Extending superposed harmonic initial data to higher spin
EExtending superposed harmonic initial data to higher spin
Sizheng Ma, ∗ Matthew Giesler,
2, 1
Mark Scheel, and Vijay Varma
3, 2, 1, † TAPIR 350-17, California Institute of Technology, 1200 E California Boulevard, Pasadena, CA 91125, USA Cornell Center for Astrophysics and Planetary Science, Cornell University, Ithaca, New York 14853, USA Department of Physics, Cornell University, Ithaca, New York 14853, USA (Dated: February 15, 2021)Numerical simulations of binary black holes are accompanied by an initial spurious burst ofgravitational radiation (called ‘junk radiation’) caused by a failure of the initial data to describe asnapshot of an inspiral that started at an infinite time in the past. A previous study showed thatthe superposed harmonic (SH) initial data gives rise to significantly smaller junk radiation. However,it is difficult to construct SH initial data for black holes with dimensionless spin χ (cid:38) .
7. We hereprovide a class of spatial coordinate transformations that extend SH to higher spin. The new spatialcoordinate system, which we refer to as superposed modified harmonic (SMH), is characterized by acontinuous parameter — Kerr-Schild and harmonic spatial coordinates are only two special cases ofthis new gauge. We compare SMH with the superposed Kerr-Schild (SKS) initial data by evolvingseveral binary black hole systems with χ = 0 . .
9. We find that the new initial data still leadsto less junk radiation and only small changes of black hole parameters (e.g. mass and spin). We alsofind that the volume-weighted constraint violations for the new initial data converge with resolutionduring the junk stage ( t (cid:46) M ), which means there are fewer high-frequency components inwaveforms at outer regions. I. INTRODUCTION
The detection of GW150914 [1] and other binary com-pact objects [2–6] has opened a new era in astrophysics.With the improvement of detector sensitivity, more andmore events are expected to be detected in the near fu-ture [7]. Therefore, an accurate modeling of coalescingbinaries is crucial for data analysis. Numerical relativity(NR) remains the only ab initio method to simulate thecoalescence of binary black hole (BBH) systems. WithNR, one can obtain the entire BBH waveform includinginspiral, merger, and ringdown. Moreover, gravitationalwave models [8–15] used to analyze detector data areultimately calibrated against NR.Numerical simulations of BBHs are based on splittingthe Einstein equation into constraint and evolution parts,where the constraint equations provide the initial datato evolve. However, the constructed initial data doesnot exactly correspond to a quasi-equilibrium state of aninspiral that started at an infinite time in the past. Forexample, the tidal distortion of a BH is not fully recovered,and the initial data do not usually include gravitationalradiation already present. As a result, once the evolutionbegins, the system relaxes into a quasi-equilibrium state,and gives rise to a pulse of spurious radiation, which isreferred to as ‘junk radiation’. Several attempts havebeen made to reduce junk radiation, by introducing PNcorrections [16–21], or by using a curved conformal metric[18, 22, 23].Recently, Varma et al. [23] carried out a systematicstudy of initial data and its effects on junk radiation andcomputational efficiency of the subsequent time evolution. ∗ [email protected] † Klarman fellow
The simulations studied in Varma et al. were performedwith an NR code: the Spectral Einstein Code (SpEC)[24], where the construction of initial data is based on theExtended Conformal Thin Sandwich (XCTS) formulation[25, 26]. Within this formalism, several free fields, includ-ing the conformal metric, must be provided. Differentchoices of the free fields generate different physical initialdata; the data still correspond to two black holes with thesame desired mass ratio and spins, but the initial tidaldistortions and strong-field dynamics differ. Varma et al. showed that the junk radiation and efficiency of the sub-sequent evolution depend on the given free fields. In par-ticular, choosing the initial data based on two superposedblack holes in time-independent harmonic coordinates[27] (heretofore called superposed harmonic (SH) data)leads to less junk radiation than superposed Kerr-Schild(SKS) initial data [22], which is typically used in SpECsimulations [28]. Varma et al. also found that SH initialdata has higher computational efficiency. However, SHinitial data works well only for BHs with dimensionlessspin χ (cid:46) .
7. For high-spin BHs, the horizons becomeso highly deformed that it is difficult to construct initialdata (cf. Fig. 10 in Ref. [23]).For both SH and SKS initial data, the conformal spa-tial metric and the trace of the extrinsic curvature aredetermined by superposing the analytic solutions for twosingle Kerr black holes. The difference is that SKS usesthe Kerr metric in Kerr-Schild coordinates, and SH usesthe Kerr metric in time-independent harmonic coordi-nates [27]. It may be surprising that making a differentcoordinate choice — the choice of coordinates for thesingle-BH analytic solution — leads to a different physi-cal BBH solution. The reason is that the superpositionof two single-BH solutions does not solve the Einsteinequations for a BBH and is used to compute only some ofthe fields; the remaining fields are computed by solving a r X i v : . [ g r- q c ] F e b constraints and by quasi-equilibrium conditions. For asingle black hole, following the complete initial data pro-cedure (including solving the constraints numerically) forboth SKS and SH would result in the same physical Kerrmetric but in different coordinates.In this paper, we extend SH to higher spins by usinga spatial coordinate map to transform the free data forthe single-BH conformal metric, while retaining harmonictime slicing for this single-BH conformal metric. Thecoordinate transformation defines a class of spatial co-ordinate systems that are characterized by a continuousparameter α . We refer to these coordinates as the modi-fied harmonic (MH) coordinate system. MH coordinatesare purely harmonic with α = 1 and correspond to spatialKS when α = 0. Similar to the cases of SKS and SH, aninitial data for a BBH system can also be constructedby superposing two single Kerr black holes in MH co-ordinates. We refer to this initial data as superposedmodified harmonic (SMH). For the BBH systems with χ > .
7, a value of α < α as close to 1 as possibleso that SMH data still shares the desirable properties ofSH initial data.This paper is organized as follows. In Sec. II, we providesome basic information about how we compute initialdata and evolve BBH systems. In Sec. III, we comparethe behavior of different single-BH coordinate systems.In particular, in Sec. III E we explicitly point out thenumerical reason that SH does not work for high-spin BHs.This immediately leads to a class of spatial coordinatetransformations, defined in Sec. III F, that can cure thenumerical issues. We then use the MH coordinate systemto construct initial data for BBHs (i.e., SMH) with χ = 0 . . II. BBH INITIAL DATA AND EVOLUTION
Following the discussions in Ref. [23], we use the XCTSformulation to construct initial data for a binary blackhole system. Within this formalism, one can freely specifythe conformal metric ¯ g ij , trace of extrinsic curvature K , and their time derivatives ∂ t ¯ g ij and ∂ t K . To obtainquasi-equilibrium initial data, we choose ∂ t ¯ g ij = 0 , ∂ t K = 0 . (1)The construction of the other free fields, ¯ g ij and K , isbased on the 3-metric g βij and the trace of extrinsic cur-vature K β of two single boosted Kerr BHs , where thesuperscript β = 1 , g ij = f ij + X β =1 e − r β /w β ( g βij − f ij ) , (2) K = X β =1 e − r β /w β K β , (3)where f ij is the flat 3-metric, and r β is the Euclideancoordinate distance from the center of each BH [29]. Notethat each metric is weighted by a Gaussian with width w β = 0 . d L β , (4)where d L β is the Euclidean distance between the Newto-nian L Lagrange point and the center of the black holelabeled by β . Here g βij and K β correspond to the Kerrsolution expressed either in the KS, harmonic, or MHcoordinate systems. BBH initial data constructed fromthe two Kerr solutions in the aforementioned coordinatesare referred to as SKS, SH, and superposed modifiedharmonic (SMH), respectively.After specifying the free fields, the initial data are com-pleted by solving a set of coupled elliptic equations thatensure satisfaction of the constraints and an additionalquasi-equilibrium condition. Additionally, these ellip-tic equations require boundary conditions. At the outerboundary (typically chosen to be 10 M from the sources),we impose asymptotic flatness [cf. Eq. (11)—(13) in Ref.[23]], and at each inner boundary we enforce an apparenthorizon condition [cf. Eq. (15)—(24) in Ref. [23]]. Aftergenerating initial data in the XCTS formalism, we alsoneed to specify the initial gauge for time evolution. Herewe use the most common choice for SpEC simulations: ∂ t N = ∂ t N i = 0 in a corotating frame, where N is thelapse function and N i is the shift vector. It was shownthat the damped harmonic gauge [30] is the most suitablefor mergers, so we do a smooth gauge transformation ona time scale of ∼ M during the early inspiral, to trans-form from the initial gauge to the better suited dampedharmonic gauge. III. MODIFIED HARMONIC COORDINATESYSTEM
In this section, we aim to investigate the reason thatmakes the harmonic coordinates problematic for high-spinBHs. We begin with a brief review of KS coordinatesin Sec. III A. Then in Sec. III B, we outline a methodthat can be used to study the numerical behavior ofKerr metric in different coordinate systems. It is thenapplied to KS spatial coordinates with harmonic slicingin Sec. III C, and to harmonic coordinates in Sec. III D.Those analyses allow us to explicitly show the numericalproblem with using harmonic coordinates for high-spinBHs, as discussed in Sec. III E. Finally in Sec. III F, weprovide a coordinate map to fix the problem.
A. Kerr in Kerr-Schild coordinates
For a stationary Kerr BH with mass M and angularmomentum χM in the z direction, the metric in KScoordinates x µ KS = ( t KS , x KS , y KS , z KS ) is given by [31] ds = g µν dx µ KS dx ν KS = ( η µν + 2 Hl µ l ν ) dx µ KS dx ν KS , (5)where η µν is the Minkowski metric, H is a scalar function,and l µ is a null covariant vector. The expressions for H and l µ are not used here but can be found in Ref. [31].With KS coordinates, the radial Boyer-Lindquist coordi-nate r can be written as [31] r = 12 ( x + y + z − a )+ (cid:20)
14 ( x + y + z − a ) + a z (cid:21) / , (6)or equivalently x + y r + a + z r = 1 . (7)Here we have used a = χM for the sake of conciseness.The outer and inner horizons of the BH are located at r ± = M ± p M − a . (8) B. Transforming from KS to a different coordinatesystem
Now we introduce a new coordinate system x µ =( t, x, y, z ), which are related to the KS coordinates x µ KS through dt KS dx KS dy KS dz KS = (cid:18) b C (cid:19) dtdxdydz , (9)where b is a 3D vector, and C is a 3 × t KS1 . Note that we here keep the formsof b and C generic, so that our present discussion can beapplied to different coordinate systems.With the Jacobian at hand, we could transform the Kerrmetric into the new coordinates, and study the numericalfeatures of each metric component, such as the problem-atic behavior of harmonic coordinates for high-spin blackholes, but this usually involves very complicated calcula-tions. However, since g KS µν can be decomposed into two Equivalently, ( x KS , y KS , z KS ) are independent of t . pieces [Eq. (5)], it is simpler to study the transformationsof η µν . In the new coordinates, we have η µν = (cid:18) − I (cid:19) → (cid:18) b C (cid:19) T (cid:18) − I (cid:19) (cid:18) b C (cid:19) = (cid:18) − − b − b T C T C − b T b (cid:19) , (10)where I is the three-dimensional identity matrix. Boththe 3-metric ( C T C − b T b ) and the shift vector − b aboveare modified by the vector b . Any numerically problematicterm in b might cause difficulty to resolve the metric inthe new coordinates. Below, we focus on the z componentof b , b z , at the inner boundary r = r + , and study itsnumerical behavior for high-spin black holes (especiallywhen a → M ) with several coordinates. C. Kerr-Schild spatial coordinates with harmonicslicing
We first apply our discussion in Sec. III B to a mixedcoordinate system: KS spatial coordinates together withharmonic temporal slicing, then we have C KSHS = I , (11a) b KSHS = 2 Mr − r − ∇ r, (11b)where the subscript ‘KSHS’ stands for Kerr-Schild spatialcoordinates with Harmonic Slicing; and r is the radialBoyer-Lindquist coordinate. Note that Eq. (11b) is theresult of [27] t KSHS = t KS − Z Mr − r − dr . (12)We refer the reader to Appendix A for the detailed ex-pression of ∇ r . The z component of b KSHS at the innerboundary r = r + is given by (as a → M ) b z KSHS = M z KS r + ( az KS ) . (13) D. Harmonic coordinates
Let us turn our attention to harmonic coordinates x µ H =( t H , x H , y H , z H ), where the spatial coordinates also becomeharmonic. For such a coordinate system, we have [27]( r − M ) = 12 ( x + y + z − a )+ (cid:20)
14 ( x + y + z − a ) + a z (cid:21) / , (14) We have checked that the same problematic terms also occur inthe Hl µ l ν piece of Eq. (5). and x + y ( r − M ) + a + z ( r − M ) = 1 , (15)where the subscript ‘H’ stands for harmonic coordinates.The harmonic slicing implies b H = 2 Mr − r − ∇ r, (16)with z component of b H at r = r + given by (as a → M ) b z H = M z H ( r + − M ) + ( az H ) . (17)Expressions for the 3 × C H ) ij = ∂x i KS /∂x j H , along with additional details, can be found in AppendixA. E. Problematic behavior of harmonic coordinates
In SpEC, the Legendre polynomials are used to numer-ically expand b z H and b z KSHS as functions of cos θ , definedby cos θ = z H p x + y + z . Here θ is the polar angle in harmonic coordinates andis not to be confused with the angular Boyer-Lindquistcoordinate. As a test, we first represent b z H [Eq. (17)] withtwenty Legendre-Gauss collocation points and a BH spinof a = 0 . M . The results of this test are shown in Fig. 1.From Fig. 1, we see that the function b z H is difficult toresolve using Legendre polynomials. This is the primaryreason that harmonic coordinates fail to accurately repre-sent high-spin BH initial data. Note that increasing theresolution to l (cid:38)
60 (for a single BH) eventually allowsus to resolve b z H , but in practice requiring such high res-olution is computationally prohibitive; furthermore, therequired resolution increases rapidly as the spin increases.Previous studies have shown success in high-spin BBHsimulations with SKS initial data up to spins of χ = 0 . b z KSHS (see Sec. III C) in which the time coordinate is harmonicbut the spatial coordinates are Kerr-Schild. Again, werepresent b z KSHS with twenty Legendre-Gauss collocationpoints and a BH spin of a = 0 . M , as shown in Fig.1. The representation is much better than the case ofharmonic coordinates. And we also confirm that withsuch mixed coordinates, BBH initial data can be indeedextended to higher spins. However, as we show later, theydo not lead to a smaller amount of junk radiation thanSKS initial data. cos θ U b z MH ( α = 0 . b z H b z KSHS b z MH ( α = 0 . Fit b z H Fit b z KSHS
Fit
FIG. 1. The function b z in KSHS [Eq. (13)], harmonic [Eq.(17)] and MH [Eq. (20)] coordinates with α = 0 .
7. Solidlines represent b z , whereas triangles represent the Legendre-Gauss collocation approximation to each function b z using 20Legendre polynomials. The spin of the BH is a = 0 . M . b z is better approximated by a fixed number ( l = 20) of Legendrepolynomials for MH than for harmonic coordinates. Looking more closely at Fig. 1, b z KSHS has fewer struc-tures than b z H , which makes b z KSHS easier to represent byLegendre polynomials. The main reason is that b z H ∼ z − in the case of high-spin BHs [since ( r + − M ) (cid:28) az H as a → M . See the denominators of Eq. (13) and (17)] . Inthe next subsection, we provide a new spatial coordinatesystem where the dependence of b z on cos θ is reduced. F. Modified Harmonic coordinates
We have seen that b z H is sensitive to cos θ for high-spin BHs. To reduce such dependence, we define a moregeneral coordinate system x + y ( r − αM ) + a + z ( r − αM ) = 1 , t MH = t H , (18)which leads to( r − αM ) = 12 ( x + y + z − a )+ (cid:20)
14 ( x + y + z − a ) + a z (cid:21) / . (19)Here we introduce a new constant parameter α . As men-tioned earlier, we refer to this new choice of spatial coordi-nates as the modified harmonic (MH) coordinate system.MH coordinates become harmonic ( spatial ) coordinateswhen α = 1 [Eq. (15)] and become KS ( spatial ) coordi-nates when α = 0 [Eq. (7)]. Meanwhile, the time slicingof MH coordinates is the same as in harmonic, regardlessof the value of α . With this new coordinate system, theradius of the outer horizon along the spin direction is(1 − α ) M + √ M − a . For a → M , this radius goes to M for KS coordinates ( α = 0) and it goes to zero forharmonic coordinates ( α = 1). Therefore, the horizonwith harmonic coordinates is highly compressed in thespin direction. However, if we let α be a number smallerthan, but still close to 1, the horizon will be less distorted.On the other hand, since α is close to 1, we can expectthat it still shares some similar properties (e.g. less junkradiation) with harmonic coordinates.As in Sec. III E, we use the function b z as an exampleto see the improvement offered by MH coordinates. Inthe MH coordinate system, we have b z MH = M z MH ( r + − αM ) + ( az MH ) . (20)Now ( r + − M ) is replaced by ( r + − αM ) . If α < r + − αM ) is not as small as ( r + − M ) , so one cananticipate that b z MH should also have fewer structuresthan b z H , just like b z KSHS . As a result, it should be easierfor Legendre polynomials to represent b z MH . To test this,we repeat our previous process and represent it with thesame set of angular Legendre-Gauss collocation points,by taking α = 0 . a = 0 . M . Fig. 1 shows thecomparison. As expected, b z MH is less sensitive to cos θ than b z H , and the representation in Legendre polynomialsshows an enormous improvement. IV. RESULTS
In this section we investigate the numerical behavior ofBBHs evolved starting with SMH initial data, comparedto evolution of SKS data. We pick four cases, as sum-marized in Table I. To make comparisons, we considerconstraint violations, computational efficiency, changesof BH parameters (mass and spin), and junk radiation.For the first three factors, we show the general featuresof SMH by focusing on Case I. For junk radiation, westudy all cases. For each simulation, we evolve with threeresolutions (labeled Lev 1,2,3 in order of increasing reso-lution). The resolution is chosen by specifying differentnumerical error tolerances to the adaptive mesh refine-ment (AMR) algorithm [33]. The orbital eccentricity isiteratively reduced to below ∼ − [34]. A. Constraint violations and computationalefficiency
Figure 2 shows the evolution of the volume-weightedgeneralized harmonic constraint energy N volume , which isgiven by [Eq. (53) of Ref. [35]] N volume = s R V F ( x ) d x R V d x , (21) -9 -8 -7 -6 -5 C o n s t r a i n t SKS Lev3SKS Lev2SKS Lev1 SMH α = 0 . Lev3SMH α = 0 . Lev2SMH α = 0 . Lev1 t/M
FIG. 2. The volume-weighted generalized harmonic constraintenergy for evolutions of Case I, with both SKS (dotted lines)and SMH (solid lines) initial data. Three resolutions areshown, labeled ’Lev1’ (red), ’Lev2’ (blue), and ’Lev3’ (black)in order of decreasing AMR tolerance. At the beginning,BHs of SMH initial data are more distorted on the grid sothe constraints are worse. However, as the gauge transitionproceeds, the constraints decay quickly. During most of thejunk stage (25 M (cid:46) t (cid:46) M ), the constraints of SMH initialdata are smaller than SKS by an order of magnitude. Theyalso converge with resolution. After the junk stage, SKS andSMH finally become comparable. -6 -5 -4 -3 -2 -1 C o n s t r a i n t SKS Lev3SKS Lev2SKS Lev1 SMH α = 0 . Lev3SMH α = 0 . Lev2SMH α = 0 . Lev1 t/M
FIG. 3. Same as Fig. 2, except that L norm is used. with F ( x ) the generalized harmonic constraint energy at x . For the first ∼ M of evolution, the constraints ofSMH are much larger than those of SKS. This is becauseBHs with SMH initial data are more distorted than SKS,and the metric is more difficult to resolve; however, themetric is much easier to resolve for SMH than SH (whichis not shown because even constructing the initial datafor SH is problematic with a spin of χ = 0 . TABLE I. A summary of parameters (mass ratio q and dimen-sionless spins χ ) for four simulations, where the spins of CaseII are chosen randomly. The orbital angular momentum ispointing along (0 , , α for MH coordinates.Simulation label q χ χ α Case I 1 (0 , , .
8) (0 , , .
8) 0.9Case II 1 (0 . , . , .
50) (0 . , . , .
46) 0.9Case III 2 (0 , , .
7) (0 , , .
8) 0.8Case IV 1 (0 , , .
9) (0 , , .
9) 0.7 N g r i d SKSSMH α = 0 . t/M C P U - h FIG. 4. Computational efficiency of evolutions of SMH ( α =0 .
9) and SKS initial data for Case I, with the highest resolution.The upper panel is the total number of grid points as a functionof time. At the beginning, the SMH initial data requires manymore grid points to meet the error tolerance. As the gaugetransition to damped harmonic gauge proceeds (on a timescale of ∼ M ), the BHs become less distorted, so AMRgradually drops points. At the same time, several concentricspherical shells around each of the BHs are dropped, whichleads to discontinuous jumps in the number of grid points. Inthe end, evolutions of SMH initial data has fewer collocationpoints than for SKS. The lower panel is the accumulated CPUhours versus time. The SMH initial data is extremely slow atthe beginning. As the collocation points and subdomains areadjusted, it speeds up. The total CPU hours for evolutions ofboth initial data sets are similar. more, at slightly later times, constraints decrease rapidly.During the junk stage ( t (cid:46) M ), the constraints forthe evolution of SMH initial data are smaller than thoseof SKS by an order of magnitude. After the junk leavesthe system, the evolution of SMH initial data is still alittle bit better than that of SKS initial data, althoughconstraints of SKS and SMH become similar at late times( t (cid:38) M for Lev 3 and t (cid:38) M for Lev 2). Lookingmore closely at the junk stage, we find the constraints forevolutions of SMH initial data converge with resolution,which was also observed for SH with low-spin BHs [23].However, Fig. 3 shows the norm of the constraint energy t/M ∆ t / M SKSSMH α = 0 . FIG. 5. The time step as a function of evolution time. Theresolution is Lev 3. Initially, the time step for evolutions ofSKS initial data is larger than for SMH. However, after severaljumps due to the shell-dropping algorithm, SMH eventuallyhas a larger time step than SKS. determined using a pointwise L norm over grid pointsrather than an integral over the volume, as given by N pointwise = vuuut N P i =1 ( F i ) N , (22)where the subscript i stands for the index of grid point,and N is the total number of grid points. For the point-wise norm the convergence is not as good as for thevolume-weighted norm. This is because the pointwisenorm gives larger weight to the interior regions near theBHs where there are more points, whereas the volumenorm gives larger weight to the exterior wave zone whichcovers more volume. The difference beween Figs. 2 and 3illustrates that the improvement of the constraints inthe case of SMH mainly comes from the outer region,where the high-frequency components in the waveformsare smaller (i.e., less junk radiation). Figure 3 also showsthat the pointwise norms ( L norm) for evolutions of bothinitial data sets become comparable much earlier than thevolume norms ( t ∼ M ). This is because the pointwisenorms are monitored by AMR, and therefore their valuesremain consistent with the numerical error tolerance inAMR during the evolution as AMR makes changes to thegrid resolution.To understand how the computational efficiency of theevolution depends on the initial data, in Fig. 4 we show thetotal number of grid points in the computational domainas a function of time. At the beginning, SMH needs manymore points than SKS. As the gauge gradually transformsto the damped harmonic gauge, the BHs become lessdistorted and AMR decides to drop grid points. Duringthe evolution, there are two factors that mainly controlthe number of grid points. One is AMR, which adjustsgrid points based on the numerical error tolerance. The L e v
1e 4 ∆ M A irr /M SKSSMH α = 0 . L e v
1e 50 100 200 300 400 500 t/M L e v
1e 5 L e v
1e 4 ∆ χ A SKSSMH α = 0 . L e v
1e 40 100 200 300 400 500 t/M L e v
1e 4
FIG. 6. The evolution of irreducible mass (left) and dimensionless spin (right) of the first BH for Case I, with three resolutions.The quantities shown are deviations from their values at t = 0. Evolutions of SMH initial data have fewer oscillations than SKS.Deviations of three parameters for both initial data sets are on the same order. other one is the domain decomposition [36]. SpEC splitsthe entire computation region into various subdomains. Inparticular, there are a series of concentric spherical shellsaround each BH. The subdomain boundaries are fixed inthe “grid frame”, the frame in which the BHs do not move,but these boundaries do move in the “inertial frame”, theframe in which the BHs orbit and approach each other [37].As the separation between the BHs decreases, the inertial-frame widths of the subdomains between them decreasesas well. During the evolution, the inertial-frame widths ofthe spherical shells are monitored. Once one of the shellsbecomes sufficiently squeezed, the algorithm drops one ofthe shells and redistributes the computational domain. InRef. [23], the authors pointed out that evolutions of SHinitial data are faster than for SKS initial data. However,that statement is not true at very early times, when SHstarts with more spherical shells and more grid points,which leads to low speed. The evolution of SH initial datathen gradually speeds up after several spherical shellsare dropped, and eventually becomes faster than thecorresponding evolution of SKS data. Our simulation hereis similar. In Fig. 4, AMR modifies N grid smoothly, whilethe discontinuous jump is caused by the shell-droppingalgorithm. For each BH, we have six spherical shellsinitially. However, four of them are dropped during thefirst ∼ M . In the end, the number of grid pointsfor evolutions of SMH is smaller than for evolutions ofSKS. This not only improves the computational efficiencyof each time step, but also increases the time step ∆ t allowed by the Courant limit (∆ t ∼ N − ). As shownin Fig. 5, the time step for SMH jumps several timesbecause of the shell-dropping algorithm. In the end, ∆ t for SMH is larger than the one for SKS. Both N grid and∆ t contributes to the high speed of evolutions of SH andSMH initial data. And we have checked that the increaseof ∆ t plays the major role in the speed increase.The bottom panel of Fig. 4 shows the accumulatedCPU hours of the simulation. At first, the evolution ofSMH is extremely slow. Once several shells are dropped,the simulation gradually speeds up. This suggests thatboth SH and SMH initial data start with more shellsthan necessary. Therefore, it might be possible to furtherimprove the computational efficiency solely by reducingthe number of shells. B. Junk radiation and changes in parameters
Since the BHs in the initial data are not in true quasi-equilibrium, the masses and spins of BHs relax once theevolution begins, resulting in slight deviations from theirinitial values. In Figure 6, we show the change of irre-ducible mass ∆ M irr ( t ) = | M irr ( t ) − M irr ( t = 0) | and thechange of spin ∆ χ ( t ) = | χ ( t ) − χ ( t = 0) | as functions oftime, for three resolutions. We can see the variations areon the same order for both SMH and SKS initial data, butSMH has smaller oscillations. With the highest resolution,the deviation of SMH is smaller by a factor of ∼ . − h lm , for Case I, II and III listed in TableI (Case IV will be discussed later). Note that the lineargrowth of h for Case II appears because only the initialpart of the waveform is shown; over the entire evolution, | r h / M |
1e 1
Case I
SMH α = 0 . SKS
Case II
SMH α = 0 . SKS
Case III
SMH α = 0 . SKS | r h / M |
1e 5 0.00.51.0 1e 2 0.000.250.500.751.00 1e 20123 | r h / M |
1e 5 024 1e 3 0.00.51.0 1e 2200 0 200 4000246 | r h / M |
1e 3 200 0 200 4000246 1e 3 200 0 200 400024 1e 3t/M
FIG. 7. Mode amplitudes of waveforms for Case I, II and III with the highest resolution. Columns correspond to three cases,and rows are for different modes. For SMH initial data, we pick α = 0 . α = 0 . h for Case II appears because only the initial part of the waveform is shown. Over the entire evolution,the mode is oscillatory. In general, the junk radiation of SMH initial data leaves the system faster. It is also smaller than thejunk radiation of SKS for most of the modes. However, there are some modes, such as h , that have the same peak as SKS. the mode is oscillatory.We can see that the junk radiation of evolutions ofSMH initial data is less than for SKS for most of themodes. In general, the junk radiation leaves the systemfaster for SMH initial data than for SKS. However, thedecrease of junk radiation for SMH is not as significantas SH for low-spin BHs [23]. Some modes of SMH initialdata, such as h , are similar to SKS. Comparing CasesII and III, we note that the junk radiation of α = 0 . α = 0 . α = 0 . α = 1). Notethat Case II has similar junk radiation as Case III whenboth cases are evolved from SKS initial data; this suggeststhat the difference in junk radiation between Cases II andIII seen in Figure 7 is probably not due to differences inparameters like the mass ratio.For Case IV, a BBH system with dimensionless spins0.9, we need to decrease α to 0.7, since for that large of | r h / M | SMH α = 0 . SKS
200 100 0 100 200 300 400 500 t/M | r h / M | FIG. 8. The h and h modes for the highest resolution ofCase IV, an equal-mass BBH system with larger spins. Thespins for both BHs are (0 , , . spin α = 0 . h and h . We can seethe junk radiation for SMH initial data is still less thanfor SKS. But the improvement is not as good as othercases. For modes other than h and h , we do not seeimprovements. The main reason appears to be that α =0 . α = 1, so that the benefit of SHinitial data is reduced. In addition, in Fig. 9 we comparethe accumulated CPU hours for evolutions of both initialdata sets. We can see the initial computational efficiencyfor SMH initial data is much lower, but it graduallycatches up after several shells are dropped. For evolutionsof only a few orbits, the expense of evolving SMH initialdata may not be worth the extra computational cost.But for evolutions of many orbits, the extra cost at thebeginning of the evolution will be comparatively small.In most of the evolutions shown here, shortly afterthe beginning of the simulation several spherical shellsaround each BH are dropped, leading to a smaller numberof grid points, a larger time step, and overall greater com-putational efficiency. However, for a general evolution,we are not always ‘lucky’ enough to gain this efficiency,since the current algorithm for dropping spherical shells aims only to avoid narrow shells rather than to speed upthe simulation. To improve the computational efficiencyfor all simulations, we could start with fewer sphericalshells at t = 0. However, the benefit of this change islimited without changing the shell-dropping algorithm.One workaround is to use smaller α , which speeds up thesimulation, but if α deviates too much from α = 1, wecannot have less junk radiation. Therefore, we suggestthat the algorithm that divides the domain in to subdo-mains should be modified to account for computationalefficiency during the evolution, or a better algorithmshould be developed to initialize subdomains. Given suchfuture algorithmic improvements, we could potentiallyrun high-spin BBH evolution with larger α , which canlead to less junk radiation. t/M C P U - h SMH α = 0 . SKS
FIG. 9. The accumulated CPU hours for evolutions of SMHand SKS initial data as functions of time. The BBH system isCase IV, and we plot results for the highest resolution. Theinitial computational efficiency of SMH initial data is muchlower than for SKS, but after a short time both evolutionsproceed at the same number of CPU hours per simulationtime.
V. CONCLUSION
In this paper, we extended SH initial data [23] to higher-spin BBHs by introducing a class of spatial coordinatesystems that represent a time-independent slicing of a sin-gle Kerr black hole and are characterized by a continuousparameter α . This coordinate representation of Kerr isused to supply free data for the initial-value problem forBBH systems; we call the resulting initial-value solutionSMH initial data. The harmonic ( α = 1) and KS ( α = 0)coordinate representations of Kerr are only two specialcases of our new representation. The coordinate shapeof the horizon becomes less spherical and more distortedfor larger α . Therefore for high-spin BHs, we pick α < α should be close to 1 so that SMH initial data stillhas the desirable properties of SH initial data as shownin Ref. [23], such as less junk radiation. We have testedthat for SMH initial data with α = 0, i.e, harmonic timeslicing with KS spatial coordinates, there is more junkradiation than for SKS initial data.We have evolved four BBH systems with dimensionlessspins 0.8 or 0.9 starting from SMH initial data with α between 0.7 and 0.9, and we compared with evolutions ofthe same system starting from SKS initial data. The firstthree cases, all with dimensionless spins 0.8, representdifferent situations: a non-precessing system with equalmasses, a precessing system with random spin directions,and a non-precessing system with unequal masses. Ingeneral, the junk radiation of SMH initial data leaves thesystem faster than that of SKS. For most gravitationalwave modes, the SMH initial data leads to less junkradiation. The exceptions, like the h mode and the h mode for Case III, have bursts with amplitudes similar toSKS. Furthermore, α = 0 . α = 0 . L -norm constraints do not have suchconvergence. Therefore, the benefit is mainly from theouter regions, where there is less junk radiation.At the beginning of the evolution for Case I, SMHrequires more collocation points than SKS to reach theerror tolerance because the horizon is distorted, hence itproceeds more slowly. At later times, SKS and SMH runat approximately the same rate, after both the compu-tational efficiency on each time slice and the size of thetime step increase for the SMH case.For Case IV, which has BHs with dimensionless spin0 .
9, we found that we needed to decrease α to 0 .
7. Wesimulated an equal-mass BBH system with equal dimen-sionless spins χ , = (0 , , .
9) and compared h and h for both SMH and SKS initial data sets. Junk radiationfor SMH is still less than for SKS, but the improvementis not as good as the case of lower spin. The comparisonof CPU hours for these two cases show that the initialcomputation efficiency for SMH initial data is much lower.But it gradually becomes the same as SKS after severalshells are dropped.We also found that the algorithm for choosing the num-ber and sizes of subdomains in SpEC could use someimprovement, particularly for the initial choice of sub-domains and the early stages of the evolution. In mostsimulations but not all, AMR eventually chooses a subdo-main distribution that increases computational efficiency.Some improvements can be gained by simply startingwith fewer spherical shells around each BH, but we find that the effects of this change are limited. Therefore, theevolution of SMH initial data for high-spin BBH will ben-efit from either an algorithm to adjust subdomain sizesbased on computational efficiency during the evolution,or a better algorithm to initialize subdomains. Thosealgorithmic improvements could allow us to run high-spinBBH evolutions with larger α , which can give rise to lessjunk radiation. ACKNOWLEDGMENTS
We want to thank Maria Okounkova, Saul Teukolskyand Harald Pfeiffer for useful discussions. M.G. is sup-ported in part by NSF Grant PHY-1912081 at Cornell.M.S. and S.M. are supported by NSF Grants No. PHY-2011961, PHY-2011968, and OAC-1931266 at Caltech.V.V. is supported by a Klarman Fellowship at Cornell.M.G. and V.V. were supported by NSF Grants No. PHY-170212 and PHY-1708213 at Caltech. S.M, M.G, M.S,and V.V. are supported by the Sherman Fairchild Foun-dation. The computations presented here were conductedon the Caltech High Performance Cluster, partially sup-ported by a grant from the Gordon and Betty MooreFoundation. This work was supported in part by NSFGrants PHY-1912081 and OAC-1931280 at Cornell.
Appendix A: Details of MH coordinates
For a Kerr BH with an arbitrary spin vector a , thetransformations between KS spatial coordinates and MHspatial coordinates are given by x KS = a + r ( r − αM ) a + ( r − αM ) x MH + αMa + ( r − αM ) ( x MH × a )+( x MH · a ) a αM ( r − αM )[ a + ( r − αM ) ] , (A1)where a = a · a , and r is the radial Boyer-Lindquistcoordinate. For α = 0, we have x KS = x MH , i.e., the iden-tity transformation. The Jacobian C ij MH = ∂x i KS /∂x j MH C ij MH = a + r ( r − αM ) a + ( r − αM ) δ ij + αMa + ( r − αM ) a k (cid:15) ijk + a i a j αM ( r − αM )[ a + ( r − αM ) ]+ M α [ a − ( r − M α ) ][ a + ( r − M α ) ] x i MH ∂ j r − M α ( r − M α )( a + ( r − M α ) ) x MH m a k (cid:15) imk ∂ j r − x m MH a m a i ∂ j r M α [ a + 3( r − M α ) ][ a + ( r − M α ) ] ( r − αM ) , (A2)where (cid:15) ijk is the Levi-Civita symbol, δ ij is the Kroneckerdelta, and the Einstein summation convention is used.For α = 1, C ij MH becomes C ij defined in Sec. III D. Bydifferentiating Eq. (19), we have ∂ i r = x MH i + ( a · x MH ) a i / ( r − αM ) r − αM ) h − x MH · x MH − a r − αM ) i . (A3) With MH coordinates, the null covariant vector l in Eq.(5) can be written as l = (cid:18) dt MH + 2 Mr − r − dr (cid:19) + ( r − αM ) x MH − a × x MH + ( a · x MH ) a / ( r − αM )( r − αM ) + a · d x MH , (A4)where the first bracket corresponds to dt KS [see Eq. (12),with t MH = t H ]. The scalar function H in Eq. (5) is givenby H = M r ( r − αM ) r ( r − αM ) + ( a · x MH ) . (A5)In addition, the lapse function N and the shift vector N i in MH coordinates are given by N − = 1 + 2 M ( r − αM ) r ( r − αM ) + ( a · x MH ) r + ( r + 2 M ) r + r − r − , (A6) N i = N r l i + N φ a j x MH k (cid:15) jki a , (A7)with N r = N M r + ρ , N φ = − N aρ Mr − r − , (A8) ρ = r + a cos θ = r + ( a · x MH ) ( r − αM ) , l i = l i , (A9)where θ is the polar Boyer-Lindquist coordinate, and l i isthe spatial component of the null covariant vector l . [1] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev.Lett. , 061102 (2016), arXiv:1602.03837 [gr-qc].[2] B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X , 031040 (2019), arXiv:1811.12907 [astro-ph.HE].[3] R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. , 101102 (2020), arXiv:2009.01075 [gr-qc].[4] https://gracedb.ligo.org/superevents/S190412m/view/ ().[5] https://gracedb.ligo.org/superevents/S190425z/view/ ().[6] https://gracedb.ligo.org/superevents/S190814bv/view/ ().[7] J. Abadie et al. (LIGO Scientific, VIRGO), Class. Quant.Grav. , 173001 (2010), arXiv:1003.2480 [astro-ph.HE].[8] V. Varma, S. E. Field, M. A. Scheel, J. Blackman,D. Gerosa, L. C. Stein, L. E. Kidder, and H. P. Pfeiffer,Phys. Rev. Research. , 033015 (2019), arXiv:1905.09300[gr-qc]. Here we do not distinguish upper and lower indices of a tensor ina Euclidean space. [9] V. Varma, S. E. Field, M. A. Scheel, J. Blackman, L. E.Kidder, and H. P. Pfeiffer, Phys. Rev. D , 064045(2019), arXiv:1812.07865 [gr-qc].[10] S. Ossokine et al. , (2020), arXiv:2004.09442 [gr-qc].[11] R. Cotesta, A. Buonanno, A. Bohé, A. Taracchini, I. Hin-der, and S. Ossokine, Phys. Rev. D , 084028 (2018),arXiv:1803.10701 [gr-qc].[12] A. Bohé et al. , Phys. Rev. D95 , 044028 (2017),arXiv:1611.03703 [gr-qc].[13] G. Pratten et al. , (2020), arXiv:2004.06503 [gr-qc].[14] C. García-Quirós, M. Colleoni, S. Husa, H. Estellés,G. Pratten, A. Ramos-Buades, M. Mateu-Lucena, andR. Jaume, (2020), arXiv:2001.10914 [gr-qc].[15] G. Pratten, S. Husa, C. Garcia-Quiros, M. Colleoni,A. Ramos-Buades, H. Estelles, and R. Jaume, (2020),arXiv:2001.11412 [gr-qc].[16] K. Alvi, Phys. Rev. D , 124013 (2000), arXiv:gr-qc/9912113.[17] N. Yunes and W. Tichy, Phys. Rev. D , 064013 (2006),arXiv:gr-qc/0601046.[18] N. K. Johnson-McDaniel, N. Yunes, W. Tichy, and B. J.Owen, Phys. Rev. D80 , 124039 (2009), arXiv:0907.0891 [gr-qc].[19] B. J. Kelly, W. Tichy, Y. Zlochower, M. Campanelli, andB. F. Whiting, Numerical relativity and data analysis.Proceedings, 3rd Annual Meeting, NRDA 2009, Potsdam,Germany, July 6-9, 2009 , Class. Quant. Grav. , 114005(2010), arXiv:0912.5311 [gr-qc].[20] G. Reifenberger and W. Tichy, Phys. Rev. D , 064003(2012), arXiv:1205.5502 [gr-qc].[21] W. Tichy, Rept. Prog. Phys. , 026901 (2017),arXiv:1610.03805 [gr-qc].[22] G. Lovelace, Numerical relativity data analysis. Proceed-ings, 2nd Meeting, NRDA 2008, Syracuse, USA, Au-gust 11-14, 2008 , Class. Quant. Grav. , 114002 (2009),arXiv:0812.3132 [gr-qc].[23] V. Varma, M. A. Scheel, and H. P. Pfeiffer, Phys. Rev. D98 , 104011 (2018), arXiv:1808.08228 [gr-qc].[24] .[25] J. W. York, Jr., Phys. Rev. Lett. , 1350 (1999), arXiv:gr-qc/9810051 [gr-qc].[26] H. P. Pfeiffer and J. W. York, Jr., Phys. Rev. D67 , 044022(2003), arXiv:gr-qc/0207095 [gr-qc].[27] G. B. Cook and M. A. Scheel, Phys. Rev. D , 4775(1997). [28] M. Boyle et al. , Class. Quant. Grav. , 195006 (2019),arXiv:1904.04831 [gr-qc].[29] G. Lovelace, R. Owen, H. P. Pfeiffer, and T. Chu, Phys.Rev. D78 , 084017 (2008), arXiv:0805.4192 [gr-qc].[30] B. Szilagyi, L. Lindblom, and M. A. Scheel, Phys. Rev.
D80 , 124010 (2009), arXiv:0909.3557 [gr-qc].[31] R. P. Kerr, Phys. Rev. Lett. , 237 (1963).[32] M. A. Scheel, M. Giesler, D. A. Hemberger, G. Lovelace,K. Kuper, M. Boyle, B. Szilágyi, and L. E. Kidder, Class.Quant. Grav. , 105009 (2015), arXiv:1412.1803 [gr-qc].[33] B. Szilágyi, Int. J. Mod. Phys. D23 , 1430014 (2014),arXiv:1405.3693 [gr-qc].[34] A. Buonanno, L. E. Kidder, A. H. Mroue, H. P. Pfeif-fer, and A. Taracchini, Phys. Rev.
D83 , 104034 (2011),arXiv:1012.1549 [gr-qc].[35] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, andO. Rinne, Class. Quant. Grav. , S447 (2006), arXiv:gr-qc/0512093 [gr-qc].[36] S. Ossokine, F. Foucart, H. P. Pfeiffer, M. Boyle, andB. Szilágyi, Class. Quant. Grav. , 245010 (2015),arXiv:1506.01689 [gr-qc].[37] M. A. Scheel, H. P. Pfeiffer, L. Lindblom, L. E. Kidder,O. Rinne, and S. A. Teukolsky, Phys. Rev. D , 104006(2006), arXiv:gr-qc/0607056.[38] D. Christodoulou and R. Ruffini, Phys. Rev. D4