Efficient simulations of high-spin black holes with a new gauge

We adopt the standard 3+1 form of the spacetime metric $g_{ab}$,

$\displaystyle ds^{2}=-\alpha^{2}\,dt^{2}+\gamma_{ij}\left(\beta^{i}\,dt+dx^{i}\right)\left(\beta^{j}\,dt+dx^{j}\right),$

(1)

where $\alpha$ is the lapse, $\beta^{i}$ the shift, and $\gamma_{ij}$ the spatial metric. In vacuum, the spacetime metric $g_{ab}$ and its time derivative $\partial_{t}g_{ab}$ must satisfy the Hamiltonian and momentum constraints,

	$\displaystyle R+K^{2}-K_{ij}K^{ij}$	$\displaystyle=0,$		(2)
	$\displaystyle D_{j}(K^{ij}-\gamma^{ij}K)$	$\displaystyle=0,$		(3)

where $R$ is the spatial Ricci scalar, $D_{i}$ the spatial covariant derivative, $K_{ij}$ the extrinsic curvature, and $K=K^{i}{}_{i}$ the trace of the extrinsic curvature.

The spatial metric and extrinsic curvature are split using a conformal decomposition as

	$\displaystyle\gamma_{ij}$	$\displaystyle=\psi^{4}\bar{\gamma}_{ij},$		(4)
	$\displaystyle K_{ij}$	$\displaystyle=A_{ij}+\frac{1}{3}\gamma_{ij}K,$		(5)

where $\psi$ is the conformal factor, $\bar{\gamma}_{ij}$ the conformal metric, and $A_{ij}$ the traceless part of $K_{ij}$. $A_{ij}$ is further decomposed as

	$\displaystyle A_{ij}$	$\displaystyle=\psi^{-2}\bar{A}_{ij},$		(6)
	$\displaystyle\bar{A}^{ij}$	$\displaystyle=\frac{\psi^{6}}{2\alpha}\left((\bar{L}\beta)^{ij}-\bar{u}^{ij}\right),$		(7)

where $\bar{u}_{ij}\equiv\partial_{t}\bar{\gamma}_{ij}$ (note that $\bar{\gamma}^{ij}\bar{u}_{ij}=0$ to uniquely fix $\bar{u}_{ij}$ [50]), and the vector gradient part $(\bar{L}\beta)^{ij}$ is defined as

$\displaystyle(\bar{L}\beta)^{ij}\equiv\bar{D}^{i}\beta^{j}+\bar{D}^{j}\beta^{i}-\frac{2}{3}\bar{\gamma}^{ij}\bar{D}_{k}\beta^{k},$

(8)

with $\bar{D}_{j}$ the covariant derivative associated with $\bar{\gamma}_{ij}$.

In the extended conformal thin-sandwich formalism, $\bar{\gamma}_{ij}$, $\bar{u}_{ij}$, $K$, and $\partial_{t}K$ are freely specifiable. The elliptic solver in SpEC [51] computes $\psi$, $\alpha$, and $\beta^{i}$ by solving

	$\displaystyle\bar{D}_{j}\left(\frac{\psi^{6}}{2\alpha}(\bar{L}\beta)^{ij}\right)-\bar{D}_{j}\left(\frac{\psi^{6}}{2\alpha}\bar{u}^{ij}\right)-\frac{2}{3}\psi^{6}\bar{D}^{i}K=0,$		(9)
	$\displaystyle\bar{D}^{2}\psi-\frac{1}{8}\psi\bar{R}-\frac{1}{12}\psi^{5}K^{2}+\frac{1}{8}\psi^{-7}\bar{A}_{ij}\bar{A}^{ij}=0,$		(10)
	$\displaystyle\bar{D}^{2}(\alpha\psi)-\alpha\psi\left(\frac{7}{8}\psi^{-8}\bar{A}_{ij}\bar{A}^{ij}+\frac{5}{12}\psi^{4}K^{2}+\frac{1}{8}\bar{R}\right)$
	$\displaystyle\hskip 33.96698pt+\psi^{5}(\partial_{t}K-\beta^{k}\partial_{k}K)=0,$		(11)

where $\bar{R}$ is the conformal Ricci scalar. Equations (4) and (5) are then used to compute $\gamma_{ij}$ and $K_{ij}$.

The free variables $\bar{u}_{ij}$ and $\partial_{t}K$ are typically set to 0 to construct quasi-equilibrium initial data. SpEC sets the remaining free variables $\bar{\gamma}_{ij}$ and $K$ using a superposition of two single BH spacetimes blended together by a Gaussian weight function [47]. We define $\gamma^{\rho}_{ij}$ and $K^{\rho}$ to refer to these quantities for a boosted spinning BH, where $\rho=A,B$ labels each BH. $\bar{\gamma}_{ij}$ and $K$ are chosen to be [47]

	$\displaystyle\bar{\gamma}_{ij}$	$\displaystyle\equiv\eta_{ij}+\sum_{\rho}e^{-r^{2}_{\rho}/w^{2}_{\rho}}(\gamma^{\rho}_{ij}-\eta_{ij}),$		(12)
	$\displaystyle K$	$\displaystyle\equiv\sum_{\rho}e^{-r^{2}_{\rho}/w^{2}_{\rho}}K^{\rho},$		(13)

where $\eta_{ij}$ is the 3D flat metric, $r_{\rho}$ is the Euclidean distance from BH $\rho$, and $w_{\rho}$ controls the falloff of BH $\rho$’s contribution. We use $w_{\rho}$ equal to 3/5 of the Euclidean distance between the BBH’s L1 Lagrange point (Euclidean center-of-mass) and BH $\rho$’s center. This choice of $w_{\rho}$ is wider than a BH’s size but still relatively far from the companion BH.

The most common choice for the free data at each BH is Kerr-Schild (KS), though using harmonic-Kerr has much promise for low-spin binaries [48]. In this paper, we use a Kerr-Schild-like gauge where the horizon is spherical at each black hole to set the free data. We will discuss the KS and KS-like gauges in Sec. III.

SpEC evolves the initial data using the first-order generalized harmonic (GH) system [44]. (See Refs. [52, 53, 54] for more details on the GH systems.) The coordinates $x^{a}$ (referred to as generalized harmonic coordinates) satisfy the inhomogeneous wave equation,

$\displaystyle H^{a}=\nabla_{b}\nabla^{b}x^{a}=-\Gamma^{a},$

(14)

where $\Gamma^{a}\equiv g^{bc}\Gamma^{a}_{bc}$ is the trace of the Christoffel symbol, and $\nabla_{a}$ the $g_{ab}$-compatible covariant derivative. The gauge source function $H^{a}=H^{a}(x^{b},g_{cd})$ is any arbitrary function dependent only on $x^{b}$ and $g_{cd}$ (but not derivatives of $g_{ab}$). In these coordinates, the vacuum Einstein equations can be cast into a manifestly hyperbolic form,

	$\displaystyle g^{cd}\partial_{c}\partial_{d}g_{ab}=$	$\displaystyle-2\nabla_{(a}H_{b)}$
		$\displaystyle+2g^{cd}g^{ef}(\partial_{e}g_{ca}\partial_{f}g_{db}-\Gamma_{ace}\Gamma_{bdf}),$		(15)

where $\Gamma_{abc}=g_{ad}\Gamma^{d}_{bc}$. After expanding Eq. (15) into a first-order representation (done analogously to expanding the covariant scalar field system [55, 56]) and adding constraint damping terms (see [57, 54, 58, 14, 44] for detailed discussions on constraint damping), we arrive at the GH evolution equations implemented in SpEC,

	$\displaystyle\partial_{t}g_{ab}$	$\displaystyle=-\alpha\Pi_{ab}-\gamma_{1}\beta^{i}\Phi_{iab}+\left(1+\gamma_{1}\right)\beta^{k}\partial_{k}g_{ab},$		(16)
	$\displaystyle\partial_{t}\Pi_{ab}$	$\displaystyle=2\alpha g^{cd}\left(g^{ij}\Phi_{ica}\Phi_{jdb}-\Pi_{ca}\Pi_{db}-g^{ef}\Gamma_{ace}\Gamma_{bdf}\right)$
		$\displaystyle-2\alpha\nabla_{(a}H_{b)}-\frac{1}{2}\alpha n^{c}n^{d}\Pi_{cd}\Pi_{ab}-\alpha n^{c}\Pi_{ci}g^{ij}\Phi_{jab}$
		$\displaystyle+\alpha\gamma_{0}\left[2\delta^{c}{}_{(a}n_{b)}-(1+\gamma_{3})g_{ab}n^{c}\right]\left(H_{c}+\Gamma_{c}\right)$
		$\displaystyle-\gamma_{1}\gamma_{2}\beta^{i}\Phi_{iab}+\beta^{k}\partial_{k}\Pi_{ab}-\alpha g^{ki}\partial_{k}\Phi_{iab}$
		$\displaystyle+\gamma_{1}\gamma_{2}\beta^{k}\partial_{k}g_{ab},$		(17)
	$\displaystyle\partial_{t}\Phi_{iab}$	$\displaystyle=\frac{1}{2}\alpha n^{c}n^{d}\Phi_{icd}\Pi_{ab}+\alpha g^{jk}n^{c}\Phi_{ijc}\Phi_{kab}-\alpha\gamma_{2}\Phi_{iab}$
		$\displaystyle+\beta^{k}\partial_{k}\Phi_{iab}-\alpha\partial_{i}\Pi_{ab}+\alpha\gamma_{2}\partial_{i}g_{ab},$		(18)

where $g_{ab}$, $\Phi_{iab}\equiv\partial_{i}g_{ab}$, $\Pi_{ab}\equiv-n^{c}\partial_{c}g_{ab}$ are the three dynamical fields being evolved, $\gamma_{0}$, $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ the constraint damping parameters, and $n^{a}$ the future-pointing unit normal to constant-$t$ spatial hypersurfaces. See the Appendix for values of $\gamma_{0}$, $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ used in the simulations of this paper.

The simplest choice of gauge source function is the harmonic gauge, where $H^{a}=0$. The harmonic gauge dates back to Einstein’s work [59] and has been an important tool in many aspects of analytical general relativity [60, 61, 62]. Unfortunately, using a harmonic gauge condition in simulations of BBH mergers leads to explosive growth of $\gamma=\det(\gamma_{ij})$ near the apparent horizons as two BHs merge [63]. To suppress such growth, SpEC adopts the damped wave gauge or damped harmonic (DH) gauge $H^{a}=H^{a}_{\text{DH}}$ [63],

	$\displaystyle H^{a}_{\text{DH}}=\mu_{L}n^{a}\ln\left(\frac{\sqrt{\gamma}}{\alpha}\right)-\mu_{S}\frac{\beta^{i}}{\alpha}\gamma^{a}_{\ i},$		(19)
	$\displaystyle\mu_{L}=\mu_{S}=e^{-(\ln 10^{15})r^{2}/\sigma^{2}}\left[\ln\left(\frac{\sqrt{\gamma}}{\alpha}\right)\right]^{2},$		(20)

where $r$ is the Euclidean radius, and $\sigma=100M$. Note that the DH gauge reduces to the harmonic gauge to machine precision for $r\geq\sigma$.

SpEC uses the following gauge transition from the initial gauge $H_{\text{init}}^{a}$ to the DH gauge $H^{a}_{\text{DH}}$ during evolution:

$\displaystyle H^{a}=F(t)H^{a}_{\mathrm{init}}+\left[1-F(t)\right]H^{a}_{\mathrm{DH}},$

(21)

where

$\displaystyle F(t)=\begin{dcases}\exp\left[-\left(\frac{t-t_{0}}{w}\right)^{4}\right],&t\geq t_{0}\\ 1,&t<t_{0}.\end{dcases}$

(22)

Here $t_{0}$ is the start time and $w$ is the temporal width of the transition. Unless stated otherwise, we choose $t_{0}=0M$ and $w=50M$. Note that this transition function is not finely tuned. It is chosen because it decays to 0 rapidly for large $t$ and has a continuous third derivative at $t=t_{0}$.

SpEC adopts a dual-frame configuration for BBH simulations [64]. In the inertial frame, the two BHs are orbiting each other and deform as they merge. In contrast, in the grid frame, the BHs are at fixed coordinate locations and are kept approximately spherical. The two frames are related by a time-dependent analytic map determined by feedback control systems in SpEC [65].

SpEC’s domain decomposition is described in the Appendix of Ref. [66], while the adaptive mesh refinement (AMR) algorithm is described in Refs. [67, 68]. SpEC uses spherical shells around each BH. The number of spherical harmonic modes ($\ell$) used in the shells around BHs is a direct proxy for how the shape of each BH affects the computational cost of a simulation. High-spin BBH simulations use $\ell\geq 40$, which results in not only many grid points, but also close spacing between grid points. A simulation with more grid points requires more computation per time step, while a closer spacing between grid points requires a smaller time step in order to maintain stability. Both factors slow down the overall simulation. With this in mind, we seek to reduce the angular resolution needed in high-spin BBH simulations, anticipating faster simulations.

SpEC uses excision to avoid the physical singularities inside BHs. Specifically, the region within an inner boundary for each BH is excluded from the computational domain. This boundary is called an excision surface or excision boundary and lies slightly inside the apparent horizon (AH) of each BH [54, 65]. Causality prohibits any physical content in the interior region from propagating out. The excision surface has to be placed in a trapped region between the inner and outer horizons for each BH. Since the distance between the inner and outer horizons decreases as spin increases, the placement of the excision boundary becomes increasingly difficult as the spin increases. As a result, smaller time step sizes are necessary to track the apparent horizons and to keep the excision boundary inside the narrow trapped region. Thus, a gauge where horizons remain spherical for any spin should not only decrease the resolution used by AMR but also reduce the workload in tracking the apparent horizons and controlling the excision boundaries. At the outer boundary, suitable constraint-preserving boundary conditions [44] are imposed.

The Kerr metric in KS coordinates $\{t,x,y,z\}$, with spin pointing along the $z$-axis,¹¹1For spin not along the $z$-axis adding a 3D rotation suffices to determine the metric. mass $M$, and angular momentum $aM=\chi M^{2}$ is

$\displaystyle g_{ab}=\eta_{ab}+2Hl_{a}l_{b},$

(23)

where $\eta_{ab}$ is the Minkowski metric,

	$\displaystyle H$	$\displaystyle=\frac{Mr^{3}}{r^{4}+a^{2}z^{2}},$		(24)
	$\displaystyle l_{a}$	$\displaystyle=\left(1,\frac{rx+ay}{r^{2}+a^{2}},\frac{ry-ax}{r^{2}+a^{2}},\frac{z}{r}\right),$		(25)

and $r$ implicitly given by

$\displaystyle\frac{x^{2}+y^{2}}{r^{2}+a^{2}}+\frac{z^{2}}{r^{2}}=1$

(26)

is the radial coordinate in Boyer-Lindquist coordinates. Since we are only interested in astrophysical BHs, we restrict ourselves to $\chi<1$. Then for any nonzero spin, there are inner and outer horizons $r_{\pm}=(1\pm\sqrt{1-\chi^{2}})M$. In the region $r_{-}<r<r_{+}$, any object must travel radially inward, while outside this region ($r<r_{-}$ or $r>r_{+}$, excluding horizons), an object can travel both radially inward and outward. For high spins, the horizons $r_{\pm}$ become closer together and increasingly non-spherical in the KS coordinate system, requiring greater angular resolution to simulate the BHs. This suggests that we might mitigate the required resolution increase by using coordinates in which the horizons are spherical, as we now describe.

We denote the spherical KS coordinates by $\{t,\bar{x},\bar{y},\bar{z}\}$. These are related to the KS coordinates $\{t,x,y,z\}$ by

	$\displaystyle\frac{\bar{x}}{r}$	$\displaystyle=\frac{x}{\sqrt{r^{2}+a^{2}}},$		(27)
	$\displaystyle\frac{\bar{y}}{r}$	$\displaystyle=\frac{y}{\sqrt{r^{2}+a^{2}}},$		(28)
	$\displaystyle\bar{z}$	$\displaystyle=z.$		(29)

Equation (26) describes an oblate spheroid in $\{x,y,z\}$ and is equivalent to a sphere in $\{\bar{x},\bar{y},\bar{z}\}$ coordinates. That is,

$\displaystyle\bar{x}^{2}+\bar{y}^{2}+\bar{z}^{2}=r^{2}.$

(30)

The left panel of Fig. 1 shows the inner and outer horizons in KS as solid lines, and a sample excision surface as a dashed line. The right panel of Fig. 1 shows the inner and outer horizons, and excision surface but in spherical KS coordinates instead. We will abbreviate spherical KS as SphKS hereinafter.

Recall from Sec. II.3 that as the BHs inspiral, the excision regions track the BHs. The excision boundary must be inside $r_{+}$ but outside $r_{-}$ so that all information leaves the computational domain and no boundary condition must be applied. This becomes more difficult as the spin increases, partly because the space between $r_{+}$ and $r_{-}$ decreases. We attempt to reduce the work of the control system by expanding the region between the horizons by performing a radial transformation. Note that the idea of expanding the region between horizons in the initial data is not new. For example, Ref. [69] applies a fisheye radial transformation to the quasi-isotropic coordinates to expand the horizon size. We here apply a different radial transformation to the SphKS coordinates. We refer to this gauge as wide Kerr-Schild (WKS). We continue using the notation in Sec. III.1 and denote the coordinates of WKS as {$t,\tilde{x},\tilde{y},\tilde{z}$}. We introduce a new variable $\tilde{r}$ that is related to $r$ and choose the coordinate transformation between WKS and SphKS as

	$\displaystyle\frac{\tilde{x}}{\tilde{r}}$	$\displaystyle=\frac{\bar{x}}{r},$		(31)
	$\displaystyle\frac{\tilde{y}}{\tilde{r}}$	$\displaystyle=\frac{\bar{y}}{r},$		(32)
	$\displaystyle\frac{\tilde{z}}{\tilde{r}}$	$\displaystyle=\frac{\bar{z}}{r}.$		(33)

With this convention, Eqs. (26) and (30) are equivalent to

$\displaystyle\tilde{x}^{2}+\tilde{y}^{2}+\tilde{z}^{2}=\tilde{r}^{2},$

(34)

i.e. a sphere of radius $\tilde{r}$ in WKS.

Starting with SphKS, we want a radial transformation $r\rightarrow\tilde{r}$ that keeps $r_{+}$ fixed but shrinks $r_{-}$ radially inward by some factor $b$, i.e.,

	$\displaystyle\tilde{r}(r_{+})$	$\displaystyle=r_{+},$		(35)
	$\displaystyle\tilde{r}(r_{-})$	$\displaystyle=br_{-}.$		(36)

We can achieve these relations with a quadratic,

$\displaystyle\tilde{r}(r)=\frac{1-b}{r_{+}-r_{-}}r^{2}+\frac{br_{+}-r_{-}}{r_{+}-r_{-}}r,$

(37)

with $r_{-}/r_{+}\leq b\leq 1$ to ensure monotonicity, e.g., the lower bound of $b$ is $\sim 0.75$ for a spin-0.99 BH. This $b$ is a squeezing parameter that controls how far $r_{-}$ is pushed inwards after the transformation.

Unfortunately, this quadratic has the pathology that the metric would no longer be asymptotically flat, so we use the quadratic only near the BH, smoothly transitioning to $\tilde{r}=r$ far away from the BH. Specifically,

$\displaystyle r(\tilde{r})=\frac{\sqrt{B^{2}+4A\tilde{r}}-B-2A\tilde{r}}{2A}\frac{1}{1+e^{(\tilde{r}-C)/\lambda}}+\tilde{r},$

(38)

where

	$\displaystyle A$	$\displaystyle=\frac{1-b}{r_{+}-r_{-}},$		(39)
	$\displaystyle B$	$\displaystyle=\frac{br_{+}-r_{-}}{r_{+}-r_{-}}=1-Ar_{+}.$		(40)

We write the relation as a function $r(\tilde{r})$ instead of $\tilde{r}(r)$ for easier implementation in SpEC. $C$ and $\lambda$ are parameters of the sigmoid function, controlling the center and width of transition. Theoretically, we could choose $C$ and $\lambda$ to satisfy $e^{(r_{+}-C)/\lambda}\sim 10^{-15}$ so that Eq. (37) holds exactly inside $r_{+}$ within numerical precision, but such combinations always result in large Jacobians. A quantity with large derivatives (and second derivatives) needs sufficiently high resolution to be resolved, which increases computational cost. In practice, because the starting point of WKS is broadening the region between $r_{+}$ and $r_{-}$ nonlinearly, we simply choose $C=3M$ and $\lambda=8M$. This combination of $C$ and $\lambda$ maintains stable simulations without increasing computational cost too much. Figure 2 shows $\tilde{r}/r$ versus $r$ for these $C$ and $\lambda$.

When studying evolutions using SphKS initial data we noticed that the apparent horizons become spheroids around $t=50M$ when the gauge has mostly transitioned to damped harmonic. With the change to a spheroidal horizon we observe the expected increase in the required angular resolution. Since the damped harmonic (DH) condition is generally only necessary during merger, we perform simulations where we delay the transition from $H^{a}_{\mathrm{init}}$ to $H^{a}_{\mathrm{DH}}$ in an attempt to extend how long the horizons remain (nearly) spherical. We change both $t_{0}$ and $w$ in Eq. (22) and report our results in Sec. IV.5.

We evolve a spin-0.99 BH to time $4000M$ in both KS and SphKS coordinates. We keep the domain decomposition the same by using the transformation Eqs. (27$-$29) from the SphKS domain, such that any sphere is mapped to a spheroid matching the given spin. AMR is disabled, ensuring the domain decomposition and resolution are unchanged throughout the simulations. We fix the radial resolution in both coordinates but vary the angular resolutions, represented by $\ell$. The number of angular grid points is then $2(\ell+1)^{2}$. We investigate four single BH simulations: three in KS with $\ell=22,26,30$ and one in SphKS with $\ell=22$.

In the top panel of Fig. 3, we show the $L_{2}$-norm of the metric errors $g_{ab}-g_{ab}^{\mathrm{analytic}}$ for all four simulations. The blue lines represent the evolution using KS coordinates, with $\ell=22,26,30$ corresponding to the solid, dashed and dash-dotted styles. The metric error decreases as the angular resolution increases, especially at early time (before $t<500M$). After $t\sim 1000M$, the metric error approaches a $10^{-7}$ error floor. This error floor appears because the numerical error gets reflected at the outer boundary and amplified when traveling inward [70].²²2The floor is decreased when we move the outer boundary farther out. The orange line is for the evolution using SphKS coordinates with $\ell=22$.³³3Simulations in SphKS coordinates with higher $\ell$’s yield nearly the same metric errors as the $\ell=22$ curve, so they are not shown. It reaches the same error floor at late time but is at least 10 times smaller than the metric error of the $\ell=22$ KS simulation (the blue solid curve).

In the bottom panel of Fig. 3, we show the GH constraint energy $\mathcal{E}_{c}$ for these four simulations. The constraint energy (or constraint violation) used in this paper is defined as

$\displaystyle\mathcal{E}_{c}=\displaystyle\sqrt{\frac{\int\mathcal{C}_{\mathrm{GH}}^{2}\sqrt{\gamma}d^{3}x}{\int\sqrt{\gamma}d^{3}x}},$

(41)

where

	$\displaystyle\mathcal{C}_{\mathrm{GH}}^{2}=\delta^{ab}$	$\displaystyle[\gamma^{ij}(\mathcal{C}_{ia}\mathcal{C}_{jb}+\delta^{cd}\mathcal{C}_{iac}\mathcal{C}_{jbd}+\gamma^{kl}\delta^{cd}\mathcal{C}_{ikac}\mathcal{C}_{jlbd})$
		$\displaystyle+\mathcal{F}_{a}\mathcal{F}_{b}+\mathcal{C}_{a}\mathcal{C}_{b}].$		(42)

Here, $\delta^{ab}$ is the 4D Kronecker delta and {$\mathcal{C}_{a}$, $\mathcal{F}_{a}$, $\mathcal{C}_{ia}$, $\mathcal{C}_{iab}$, $\mathcal{C}_{ijab}$} are the five constraints used in Ref. [44]. Note that our definition of $\mathcal{E}_{c}$ is different from the one in Ref. [44]. Among the simulations using KS coordinates, the constraint violation is smaller as $\ell$ increases. In contrast to the metric errors, the constraint violations do not hit an error floor because of the use of a constraint-preserving boundary condition [70]. We see that the SphKS evolution has constraint violation over a factor of $10^{3}$ smaller than the evolution in KS coordinates using the same angular resolution ($\ell=22$). Even compared to the KS $\ell=30$ case, the SphKS $\ell=22$ simulation still has a smaller constraint violation by a factor of 10.

Because the horizons are spherical in the SphKS gauge, spacetime quantities are constructed and evolved directly in spherical domains. In the KS gauge, quantities are evolved in spheroidal domains, so a spatial map converting spheroids to spheres is necessary in the spectral calculation of derivatives. The Jacobian and Hessian of this spatial map and its inverse can introduce errors. Thus, in Fig. 3, we see the single BH simulation in SphKS provides a more accurate result than in KS.

The KS $\ell=22$ simulation takes 978 CPU hours to reach $t=4000M$, while the SphKS $\ell=22$ simulation takes 948 CPU hours on the same number of cores. However, a better comparison is to the $\ell=26$ KS case, which takes 1438 CPU hours and still has considerably larger errors. It may not be surprising that we are able to reduce the numerical error of single BH simulations by using coordinates better adapted to the geometry of the BH, but it is reassuring to have confirmation.

We evolve three pairs of non-eccentric, non-precessing, equal-mass, equal-spin BBH systems, corresponding to spin 0.9, 0.95, and 0.99, all along the $z$-axis. Each pair consists of a run with superposed KS initial data and another run with superposed SphKS initial data. They both merge after nearly the same number of orbits and at nearly the same simulation time. The initial orbital frequency $\Omega_{0}$ and the initial rate of change of separation $\dot{a}_{0}$ are tuned separately for each run, subject to a fixed initial separation $D_{0}$. We perform this tuning by eccentricity reduction [71] to achieve a negligible eccentricity $(e<0.0007)$. The specific values of these parameters, including the number of orbits, are provided in Table 1. Furthermore, we simulate each BBH run at three resolutions, Lev-1, Lev-2, and Lev-3. For Lev-$i$, the target truncation error of the AMR algorithm is $\sim 2\times 4^{-i}\times 10^{-4}$.

We focus on the spin-0.9 and spin-0.99 simulations in this paper. Comparisons of CPU times and constraint energy [Eq. (41)] between the two gauges for the spin-0.9 and spin-0.99 simulations are shown in Fig. 4. Comparison of waveforms for the spin-0.9 and spin-0.99 Lev-3 simulations are shown in Fig. 5. Additionally, the CPU times⁴⁴4Because the number of CPUs used may vary during a simulation, CPU time is a better measure of efficiency than wallclock time and we do not include wallclock time in this paper. at the end of ringdown are recorded in Table 2. We will explain and analyze these figures and tables in greater detail.

We do not show figures for the spin-0.95 case because the spin-0.95 and spin-0.99 simulations share the same qualitative behavior. However, we still provide the CPU times of the spin-0.95 case in Table 2 to show the trend that using SphKS accelerates BBH simulations more significantly as the spin increases.

We omit a discussion of the spin-0.95 BBH data and jump directly to the spin-0.99 BBH case because the efficiency and waveform comparisons are very similar.

We evolve a 25-orbit, non-precessing, non-eccentric, equal-mass, spin-0.9 BBH system using the wide Kerr-Schild (WKS) gauge (Sec. III.2) with $b=0.95$ and $b=0.9$ (recall that $b=1$ corresponds to the SphKS gauge). The parameters of initial setup are listed in Table 3. Their values are close to the previous spin-0.9 SphKS BBH run, so we compare these three simulations together.

We find that the WKS simulations share most of the properties of the SphKS simulations. Namely, there is no consistent improvement in either CPU efficiency or constraint energy by switching from SphKS to WKS. Both the strain modes Re$(rh_{22})$ and Re$(rh_{44})$ from different gauges overlap well. However, the amount of junk radiation significantly increases when $b\neq 1$. We find that the junk is approximately doubled when $b=0.9$ compared to $b=1$. Therefore, we do not recommend the use of WKS for evolutions of high-spin BBHs.

In this section, we simulate BBH systems in SphKS where we delay the transition from the initial spherical gauge to damped harmonic (DH) gauge with the hope that this further improves efficiency. We choose $t_{0}=4000M$ and $w\in\{50M,100M,400M\}$ [see Eq. (22) for the definitions of $t_{0}$ and $w$]. All of these cases share the same initial parameters with the non-delayed-transition SphKS spin-0.9 run in Table 1. All runs have negligible eccentricity ($e<0.0007$), and we only focus on the highest resolution, Lev-3.