Anomalies in the gravitational recoil of eccentric black-hole mergers with unequal mass ratios
AAnomalies in the gravitational recoil of eccentric black-hole mergers with unequalmass ratios
Miren Radia, ∗ Ulrich Sperhake,
1, 2, † Emanuele Berti, ‡ and Robin Croft Department of Applied Mathematics and Theoretical Physics,Centre for Mathematical Sciences, University of Cambridge,Wilberforce Road, Cambridge CB3 0WA, United Kingdom California Institute of Technology, Pasadena, California 91125, USA Department of Physics and Astronomy, Johns Hopkins University,3400 N. Charles Street, Baltimore, Maryland, 21218, USA (Dated: January 28, 2021)The radiation of linear momentum imparts a recoil (or “kick”) to the center of mass of a mergingblack hole binary system. Recent numerical relativity calculations have shown that eccentricitycan lead to an approximate 25% increase in recoil velocities for equal-mass, spinning binaries withspins lying in the orbital plane (“superkick” configurations) [PRD (2020) 024044]. Here weinvestigate the impact of nonzero eccentricity on the kick magnitude and gravitational-wave emissionof nonspinning, unequal-mass black hole binaries. We confirm that nonzero eccentricities at mergercan lead to kicks which are larger by up to ∼
25 % relative to the quasicircular case. We also find thatthe kick velocity v has an oscillatory dependence on eccentricity, that we interpret as a consequence ofchanges in the angle between the infall direction at merger and the apoapsis (or periapsis) direction. I. INTRODUCTION
Gravitational waves (GWs) carry energy, angular mo-mentum and linear momentum away from the source withpotentially observable consequences. The radiated en-ergy corresponds to an often enormous mass deficit inthe source: for example the first ever detected black-hole(BH) binary merger GW150914 [2] radiated ∆ M ≈ M (cid:12) ,or about . of the total mass of the source. A tinyfraction of this energy is deposited into GW interferom-eters, thus enabling us to detect and characterize thesignal [3]. The angular momentum radiated in GWs re-duces the rotation rate of possible merger remnants and– at least in four spacetime dimensions – plays a criticalrole in avoiding the formation of naked singularities in theform of BHs spinning above the Kerr limit; see e.g. [4, 5].Therefore GW emission is a necessary ingredient of thetheory of general relativity, in the sense that it avoidsthe formation of spacetime singularities and preserves itspredictive power.In this paper we will focus on the radiated linear mo-mentum, which imparts a recoil (commonly referred to asa kick ) on the center of mass of the emitting system [6–8].Whereas GWs inevitably carry energy and angular mo-mentum – provided their sources do – the radiation oflinear momentum requires some degree of asymmetry,as realized in nonspherical supernova explosions or un-equal mass ratios and/or spin misalignments in binaryBH mergers. The inspiral of two equal-mass, nonspinningBHs, for example, radiates energy and angular momen-tum, whereas the emitted linear momentum is zero by ∗ [email protected] † [email protected] ‡ [email protected] symmetry. By turning these considerations around, wemay also regard the study of recoiling GW emitters asa guided search for characteristic, in some loose sense“asymmetric”, features in their orbital dynamics which, inturn, might help us to better understand astrophysicalsources through GW observations. A recoiling postmergerBH, for example, can induce a blue (or red) shift on partsof its GW signal that may be exploited in future GWobservations to directly measure BH kicks [9–11], andthe effect of kicks should be taken into account in futureringdown tests of general relativity with third-generationGW detectors to avoid systematic biases [12].For binary BH mergers, early estimates of the recoilspeeds of the remnant BH relied on a variety of approxi-mations, including post-Newtonian (PN) theory [13, 14],BH perturbation theory [15], the effective-one-body for-malism [16], the close-limit approximation [17, 18] andcombinations thereof [19]. Not long afterwards, duringthe numerical relativity (NR) gold rush, several groupsobtained more accurate results for the kick velocity fromthe merger of nonspinning BHs along quasicircular or-bits [20–22]. These calculations were followed by thediscovery that the merger of spinning BHs can lead tokick velocities ∼ , km/s when the spins lie in theorbital plane and point in opposite directions (“superkick”configurations [23–25]), and to even larger kicks of order ∼ , km/s when the spins are partially aligned withthe orbital angular momentum (“hang-up kick” configura-tions [26]). The probability of such large recoils occurringin nature depends therefore on spin alignment, and thishas been studied by several authors (see e.g. [27–31]).The possible occurrence of superkicks has importantconsequences for astrophysical BHs and their environ-ments [32–35]. It is pertinent to compare the recoil ve-locities obtained from NR simulations with the escapevelocities of various astrophysical environments [36]. Forexample, stellar-mass BH binaries are believed to form a r X i v : . [ g r- q c ] J a n dynamically in globular clusters [37]. In this case theescape velocities are generally O (10) km / s , smaller thanthe O (100) km / s kicks predicted for quasicircular, non-spinning binaries [21]. Then relativistic recoils can affectthe proportion of BH merger remnants that are retainedby globular clusters even if the BHs are nonspinning [38].At the other end of the scale, the recoil velocities of su-permassive BHs can be used to constrain theories of theirgrowth at the center of dark matter halos [39]. Kickedremnants in the accretion disk of an active galactic nucleusmay also lead to detectable electromagnetic counterpartsfor stellar-origin BH mergers [40, 41].As mentioned above a net gravitational recoil requiressome asymmetry in the system, so that the GW emis-sion is anisotropic. A natural way to accentuate theasymmetry is through the addition of orbital eccentricity.Early calculations in the close-limit approximation [18]predicted a kick proportional to e for small eccen-tricities, e (cid:46) . . More recently, numerical relativitycalculations led to the conclusion that eccentricity canlead to an approximate 25% increase in recoil velocities forsuperkick configurations with moderate eccentricities [1].The main goal of this study is to investigate the impactof nonzero eccentricity on the kick magnitude and thecorresponding GW emission of nonspinning, unequal-massBH binaries. As we shall see, the eccentricity has asubtle but significant effect on the kick magnitude, whichmanifests itself in corresponding patterns in the GWsignal, especially in subdominant multipoles.For isolated binary systems with large initial separa-tions, the emission of GWs acts to circularize the orbit bythe time the signal enters the frequency band of ground-based detectors. However, viable dynamical formationchannels of stellar-origin BH binaries could result in anonnegligible population of merging BHs that still re-tain moderate eccentricities at frequencies relevant forground-based GW detection (see e.g. [42–47]). Most ofthe events observed by the LIGO/Virgo Collaborationshow no evidence of significant eccentricities [48] but theextraordinary GW190521 event [49] is potentially consis-tent with an eccentricity as high as e ≈ . [50, 51].Orbital eccentricity is expected to be a distinguishingfeature of stellar-origin BH binaries that form dynamically,but a nonzero eccentricity is more likely at the low fre-quencies accessible by LISA, where gravitational radiationreaction has less time to circularize the binary [52–54]. Ifconfirmed, a nonzero eccentricity would hint at a possibledynamical origin for this event [50].Eccentricity is expected to play an even more promi-nent role for massive BH binaries: the dynamics of thesebinaries in stellar and gaseous environments is expected tolead to distinct, but generically nonzero, orbital eccentric-ities by the time the binaries enter the LISA sensitivitywindow (see [55] and references therein). Even largereccentricities are possible if BH binary coalescence occursthrough the interaction with a third BH [56].Our work is an exploration of the effect of large eccen-tricities near merger, and it differs in several ways from the catalog of eccentric, unequal-mass simulations pre-sented in Ref. [57]. While their study considered a largerrange of mass ratios (in our notation, / ≤ q ≤ ),they carried out fewer simulations for each value of q .The binaries in their simulations have initial eccentricitiessmaller than e = 0 . fifteen cycles before merger, andsince they start at larger orbital separations, their eccen-tricity will have further decreased by the time of merger.As we will see below, the larger initial eccentricities inour simulations allow us to highlight interesting periodici-ties in the emission of gravitational radiation and in thebehavior of the recoil velocity.The remainder of this paper is organized as follows. InSec. II we discuss our two numerical codes ( Lean and
GRChombo ), the computational framework, and thecatalog of simulations we produced for this study. InSec. III we present the main results of our simulations. InSec. IV we summarize these results and point out possibledirections for future work. In Appendix A, we detail ourtests for numerical accuracy and verify that our two codesgive comparable results. Finally, in Appendix B we discussthe tagging of cells for adaptive mesh refinement used inone of our numerical codes (
GRChombo ). Throughoutthis work we use geometrical units ( G = c = 1 ). II. COMPUTATIONAL FRAMEWORK ANDSET OF SIMULATIONSA. Numerical methods
The simulations reported in this work have been per-formed with the
GRChombo [58, 59] and
Lean [60] codes.We estimate the error budget of our simulations from bothcodes to be up to 3.5 %. Details of our convergence anal-yses are provided in Appendix A. Though different codeswere used for each sequence of configurations, we under-took comparison tests in order to ensure consistent results,and these can also be found in Appendix A. GRChombo setup
GRChombo [58] is a finite difference numerical relativ-ity code which uses the method of lines with fourth-orderRunge-Kutta timestepping. In contrast to previous stud-ies with
GRChombo we have implemented sixth-orderspatial stencils in order to improve phase accuracy [61].The Einstein equations are solved by evolving the co-variant and conformal Z4 (CCZ4) formulation [62] withthe prescription described in Sec. F of [63], namely thereplacement κ → κ /α , in order to stably evolve BHsand maintain spatial covariance. After this replacementand in the notation of Ref. [62], we use the constraintdamping parameters κ = 0 . , κ = 0 and κ = 1 in allsimulations. However, unlike Refs. [58] and [62], we usethe conformal factor defined by χ = det( γ ij ) − / , (1)where γ ij is the physical spatial metric. GRChombo is built on the
Chombo [64] library for solving par-tial differential equations with block-structured adaptivemesh refinement (AMR) which supports nontrivial meshhierarchies using Berger-Rigoutsos grid generation [65].The grid comprises a hierarchy of cell-centered Cartesianmeshes consisting of L + 1 refinement levels labeled from l = 0 , . . . , L , each with grid spacing h l = h / l . Giventhe AMR, the grid configuration changes dynamicallyduring the simulation. The regridding is controlled bythe tagging of cells for refinement in the Berger-Rigoutsosalgorithm [65], with cells being tagged if the tagging cri-terion C exceeds a specified threshold value t R . Detailsof the tagging criterion used in this work are providedin Appendix B. The Berger-Oliger scheme [65] is usedfor time stepping on the mesh hierarchy, and we take aCourant-Friedrichs-Lewy (CFL) factor of / in all simula-tions. Due to the inherent symmetry of the configurationsconsidered, we employ bitant symmetry in order to reducethe computational expense. Lean setup
The
Lean code [60] is based on the
Cactus com-putational toolkit [66] and uses the method of lineswith fourth-order Runge-Kutta timestepping and sixth-order spatial stencils for improved phase accuracy [61].The Einstein equations are implemented in the formof the Baumgarte-Shapiro-Shibata-Nakamura-Oohara-Kojima (BSSNOK) formulation [67–69] with the moving-puncture gauge [70, 71]. The
Carpet driver [72] providesAMR using the technique of “moving boxes”. We use bi-tant symmetry to exploit the symmetry of the simulationsand reduce computational expense. The computationaldomain comprises a hierarchy of L + 1 refinement levelslabeled from l = 0 , . . . l F , . . . , L , each with grid spacing h l = h / l . Before applying the symmetry, for l ≤ l F ,each level consists of a single fixed cubic grid of half-length R l = R / l , and for l F < l ≤ L , each level con-sists of two cubic components of half-length R l = 2 L − l R L centered around each BH. We adopt this notation forconsistency with that used to describe GRChombo . Thistranslates into the more conventional
Lean grid setupnotation (cf. Ref. [60]) as (cid:8) ( R , . . . , − l F R ) × (2 L − l F − R L , . . . , R L ) , h L (cid:9) . (2)A CFL factor of / is used in all simulations, and ap-parent horizons are computed with AHFinderDirect [73, 74]. In one departure from this rule, we enhance R by a factor 4/3for the simulations of sequence lq1:2 of Table I. pM pM D y x
FIG. 1. Schematic diagram of the initial BH binary setup foran arbitrary configuration in one of the sequences.
3. Initial data
For both codes, we use puncture data [75] of Bowen-York [76] type provided by the spectral solver of Ref. [77]in the form of the
Cactus thorn
TwoPunctures for
Lean , and a standalone version integrated into
GR-Chombo . In the latter case, we take advantage of theimprovements made in Ref. [78] to use spectral interpola-tion.
B. Black-hole binary configurations
We follow the construction of sequences of BH binaryconfigurations and the notation of [79]. In particular, wedenote by M and M the initial BH masses. Without lossof generality, since we are only considering unequal masses( M (cid:54) = M ), we take M > M and denote their sum by M = M + M . The reduced mass is µ = M M /M andto quantify the mass ratio, we either use q = M M , (3)or the symmetric mass ratio, η = µ/M . Finally, the totalArnowitt-Deser-Misner (ADM) mass [80] is denoted by M ADM .In order to construct a sequence for a fixed mass ratio,we first determine an initial quasicircular configuration.We specify the initial coordinate separation
D/M alongthe x axis, and the scale in the codes is fixed by choos-ing M = 0 . . Next, Eq. (65) in Ref. [81] is used tocalculate the initial tangential momentum of each BH, p = (0 , ± p, (as shown in Fig. 1). We use a Newton-Raphson method to iteratively solve for the Bowen-Yorkbare mass parameters that give the desired BH masses.The binding energy of this quasicircular configuration isthen computed using E b = M ADM − M. (4)The rest of the sequence with increasing orbital eccen-tricity is constructed by fixing the binding energy andgradually reducing the initial linear momentum parameter p . For a given configuration with fixed p , we iterativelysolve for the separation, D , and bare masses that givethe required binding energy and BH masses. The choiceto keep the binding energy constant as the momentumparameter (and thus the initial kinetic energy) is reducedmeans that the initial separation increases along the se-quence. This ensures an inspiral phase of comparableduration as the eccentricity increases. The initial orbitalangular momentum of the system is given by L = Dp [82].Even though D increases as p decreases, the initial an-gular momentum of the system monotonically decreasesas p decreases for all but the least one or two eccentricconfigurations in a sequence.We have parameterized the configurations within a se-quence by their initial tangential momentum p , but wewould like to measure the eccentricity of these configura-tions. Unfortunately, there is no gauge-invariant measureof eccentricity [83] and the ambiguity in any definitionis particularly pronounced in the late stages of inspiralfrom which our simulations start. Following Ref. [79], weuse the formalism in Ref. [84] to obtain a PN estimatefor the eccentricity. Note that this formalism has threeeccentricity parameters: e t , e r and e φ , and employs twodifferent types of coordinates: ADM-like and harmonic.The choice of which parameter and coordinate type to useis somewhat arbitrary. We mostly focus on the eccentric-ity parameter e t in harmonic coordinates as in Ref. [1].This estimate should be taken with a pinch of salt dueto the relatively small initial binary separations D in oursimulations. Furthermore, e t has an infinite gradient asa function of the initial orbital angular momentum inthe quasicircular limit (see Fig. 1 in Ref. [79]), such thatvalues of e t (cid:46) . are difficult to realize in practice, unlessthe BHs start from large initial distance. In the head-onlimit e t diverges, and a Keplerian/Newtonian interpreta-tion ceases to be valid. Despite these shortcomings, thisestimate provides us with a helpful approximation of theeccentricity and with a criterion to quantify deviationsaway from quasicircularity.The sequences considered in this work are given inTable I. Note that there are two sequences correspondingto the mass ratio q = 1 / . The sequence lq1:2 has alonger inspiral phase compared to the other sequences.For the nearly quasicircular configurations, the binarycompletes about six orbits before merger in the lq1:2 sequence, and about three orbits in all other sequences.The longer sequence of simulations was conducted in orderto identify any possible artifacts in the shorter sequencesdue to the exclusion of the earlier inspiral phase. Inaddition to the labelling of sequences in Table I, we referto individual simulations within a sequence by appending‘ -p ’ to the sequence label followed by a four digit integerwhich is given by p/M truncated appropriately: forexample, sq1:2-p0100 denotes the simulation in sequence sq1:2 with initial tangential momentum p = 0 . M . The ADM-like estimate of Ref. [84] differs only by a few percentfor e t (cid:46) . , and would not significantly alter our results. Sequence Code q E b /M r ex /M v c (km/s) sq2:3 GRChombo / − . sq1:2 Lean / − . lq1:2 Lean / − .
80 150 sq1:3
GRChombo / − . E b /M , andthe GW extraction radius r ex . For reference, we also list foreach sequence the kick velocities v c in the quasicircular limit.These values agree, within the numerical uncertainties, withthe results of Ref. [21]. C. Diagnostics
For all simulations, we have extracted values of theWeyl scalar Ψ on spheres of finite coordinate radiusgiven in Table I for each sequence. We also computed thedominant terms in the multipolar decomposition, Ψ ( t, r, θ, φ ) = ∞ (cid:88) (cid:96) =2 (cid:96) (cid:88) m = − (cid:96) ψ (cid:96),m ( t, r ) (cid:2) − Y (cid:96),m ( θ, φ ) (cid:3) , (5)where − Y (cid:96),m are the usual spin-weight − sphericalharmonics.Our main diagnostics are the energy, linear momen-tum and angular momentum radiated in GWs, which arecomputed directly from the extracted Ψ values on thespheres using standard methods. For completeness, wereproduce the formulae here.The radiated energy E rad is given by [85, 86] E rad ( t ) = lim r →∞ r π (cid:90) tt d t (cid:48) (cid:73) S r dΩ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) t (cid:48) −∞ d t (cid:48)(cid:48) Ψ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (6)The radiated linear momentum P rad is given by P rad ( t ) = lim r →∞ r π (cid:90) tt d t (cid:48) (cid:73) S r dΩ ˆ e r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) t (cid:48) −∞ d t (cid:48)(cid:48) Ψ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (7)where ˆ e r is the flat-space unit radial vector ˆ e r = (sin θ cos φ, sin θ sin φ, cos θ ) . (8)Finally, the radiated angular momentum J rad is given by J rad ( t ) = − lim r →∞ r π Re (cid:90) tt d t (cid:48) (cid:40)(cid:73) S r (cid:32)(cid:90) t (cid:48) −∞ d t (cid:48)(cid:48) ¯Ψ (cid:33) × ˆ J (cid:32)(cid:90) t (cid:48) −∞ d t (cid:48)(cid:48) (cid:90) t (cid:48)(cid:48) −∞ d t (cid:48)(cid:48)(cid:48) Ψ (cid:33) dΩ (cid:41) , (9)where the angular momentum operator ˆ J for spin weight s = − is given by ˆ J = (cid:18) Re ˆ J + , Im ˆ J + , ∂∂φ (cid:19) , (10)and ˆ J + = e i φ (cid:18) i ∂∂θ − cot θ ∂∂φ + 2i csc θ (cid:19) . (11)Additionally, we compute the radiated linear momen-tum from the multipolar amplitudes ψ (cid:96),m in Eq. (5) us-ing the formulae of Ref. [87]. From the symmetry ofour configurations, the z -component vanishes identically: P rad z = 0 . For the components in the orbital plane, wewrite P rad+ = P rad x + i P rad y . Then P rad+ ( t ) = ∞ (cid:88) ˜ (cid:96) =2 ˜ (cid:96) (cid:88) ˜ m = − ˜ (cid:96) P ˜ (cid:96), ˜ m + , (12)where P ˜ (cid:96), ˜ m + ( t ) = lim r →∞ r π (cid:90) tt d t (cid:48) (cid:40)(cid:32)(cid:90) t (cid:48) −∞ d t (cid:48)(cid:48) ψ ˜ (cid:96), ˜ m (cid:33) × (cid:32)(cid:90) t (cid:48) −∞ (cid:104) a ˜ (cid:96), ˜ m ¯ ψ ˜ (cid:96), ˜ m +1 + b ˜ (cid:96), − ˜ m ¯ ψ ˜ (cid:96) − , ˜ m +1 − b ˜ (cid:96) +1 , ˜ m +1 ¯ ψ ˜ (cid:96) +1 , ˜ m +1 (cid:105) d t (cid:48)(cid:48) (cid:19)(cid:27) , (13)and the coefficients a (cid:96),m and b (cid:96),m are given by a (cid:96),m = (cid:112) ( (cid:96) − m ) ( (cid:96) + m + 1) (cid:96) ( (cid:96) + 1) , (14) b (cid:96),m = 12 (cid:96) (cid:115) ( (cid:96) −
2) ( (cid:96) + 2) ( (cid:96) + m ) ( (cid:96) + m − (cid:96) −
1) ( (cid:96) + 1) . (15)We will find it helpful to define the partial sums, P ˜ (cid:96) + = ˜ (cid:96) (cid:88) ˜ m = − ˜ (cid:96) P ˜ (cid:96), ˜ m + , (16) P ≤ ˜ (cid:96) + = ˜ (cid:96) (cid:88) ˜ (cid:96) (cid:48) =2 P ˜ (cid:96) (cid:48) + . (17)In practice, we do not evaluate the limit in Eqs. (6), (7),(9) and (13), but rather just evaluate them at the finiteextraction radius r = r ex , as given in Table I. A discussionof the error this introduces is given in the following section.In order to exclude the spurious radiation inherentin Bowen-York initial data, we start the integration inEqs. (6), (7), (9) and (13) at t = 50 M + r ex . Therecoil velocity is computed from the radiated momentumaccording to v = − P rad M fin , (18)where M fin is the mass of the BH merger remnant. Thequantity M fin can be computed using energy balance: M fin = M ADM − ˜ E rad , (19) where ˜ E rad denotes the radiated energy including thespurious radiation. We similarly compute the spin ofthe final BH χ fin (which, by symmetry, must be in the z direction) using the radiated angular momentum: χ fin = L − J rad z M , (20)where the initial angular momentum is L = pD . For Lean simulations, we have compared M fin and χ fin with thecorresponding values derived from the apparent horizonproperties, and find agreement to within ≤ . . III. RESULTS
Using the framework summarized in the previous sec-tion, we have simulated four sequences of nonspinningBH binaries, characterized by their mass ratio (3) andby their binding energy (4). The parameters of thesesequences are listed in Table I. We have selected our massratios such that they cover the regime of maximum re-coil, realized for η = 0 . or q = 1 / . (cf. Fig. 3).Recall that sequences sq2:3 , sq1:2 and sq1:3 completeabout three orbits and sequence lq1:2 completes aboutsix orbits, respectively, in the quasicircular limit.Our main results are displayed in Fig. 2, where weplot for all sequences the total recoil speed v tot , varioustruncations of the multipolar contributions to the totalrecoil according to Eqs. (12)-(17), the total radiated GWenergy E rad and the dimensionless spin χ fin of the BHresulting from the merger.Let us first focus on the total recoil v tot , displayed ineach of the figure’s top panels as the blue solid line. Foreach mass ratio, the global maximum of the kick velocityis realized for moderate eccentricities e t ≈ . . We alsoillustrate this kick variation in Fig. 3, where the solid bluecurve shows the quasicircular kick as a function of thesymmetric mass ratio η according to Fit 3 in Table V ofRef. [88]. The velocity ranges obtained for our eccentricbinaries are overlayed as the vertical bars for each ofour sequences. The bar for each constant- η sequence isobtained by starting at the quasicircular limit on the rightof each panel in Fig. 2 and identifying the minimum andmaximum of v ( p ) , excluding the plunge regime to the leftof the global maximum.For our sequences sq2:3 , sq1:2 and lq1:2 , the magni-fication of the kick through moderate values of the orbitaleccentricity is similar to the enhancement by up to 25 %reported in Ref. [1] for the so-called superkick configu-rations [23, 24]. For sq1:3 the effect is milder, with a ∼
12 % amplification, but still well above the uncertaintyestimates of our simulations. On the other hand, as evi-denced by the oscillatory pattern of the function v ( p ) inFig. 2, appropriate nonzero values of the eccentricity canalso lead to a reduction of the maximum kick at a givenmass ratio by ∼
10 % . This overall modification of thegravitational recoil in the merger of eccentric, nonspinningBH binaries is the first main result of our study.
FIG. 2. For each sequence of simulations in Table I:
Top panel : The recoil velocity v is plotted as a function of the initialtangential momentum p/M . The individual curves represent the total kick v tot (blue, solid), the contribution to the kick from (cid:96) = 2 modes of Ψ , ψ ,m , only in Eqs. (12)-(13) v (cid:96) =2 (red, dashed), and the contributions to the kick from P ≤ ˜ (cid:96) (cid:48) + defined inEq. (17) v ˜ (cid:96) ≤ ˜ (cid:96) (cid:48) for ˜ (cid:96) (cid:48) = 2 (orange, dotted), ˜ (cid:96) (cid:48) = 3 (green, dot-dashed) and ˜ (cid:96) (cid:48) = 4 (purple, long dot-dashed). Our estimate of theeccentricity (see Sec. II B) is provided on the upper horizontal axis. Bottom panel : The final BH spin χ fin (black, solid) and theenergy radiated in GWs E rad (gold, dashed) are plotted also as functions of p/M . For both curves, the individual simulationsperformed for this analysis are shown by × symbols. Besides the global maximum, we also note a numberof local minima and maxima in the kick velocity as wevary the eccentricity in Fig. 2. For all mass ratios ( q =2 / , / , / we see about five local extrema in v ( p ) in our three short sequences, corresponding to the two upper panels and to the bottom-left panel. We notice asimilar, albeit less pronounced, oscillatory pattern in thefunctions E rad ( p ) and χ fin ( p ) for the radiated energy andfinal spin in the lower subpanels in Fig. 2. Our resultsdisplay no systematic correlation, however, between the FIG. 3. The range of recoil velocities obtained for each se-quence is plotted against the symmetric mass ratio η . Notethat for each sequence, we exclude the configurations with p < p max (i.e. the head-on limit), where p = p max is thetangential momentum that maximises the kick. The threeshort sequences are marked in gold and the long sequence ismarked in red (dashed). A fitted formula for the quasicircularkick as a function of η from Ref. [88] is also shown in blue forcomparison. extrema of the respective quantities; neither global norlocal extrema in v , E rad or χ fin coincide in magnitudeor their eccentricity values. We believe this diversityis due to the qualitatively different dependence of theradiated quantities on the GW multipoles – overlaps of different multipoles for the kick, a sum of terms ∝ ψ lm for the energy and the interaction of first and second timeintegrals for the angular momentum in Eq. (9).We added to our study the q = 1 / sequence of longerBH binary inspirals to investigate whether these anoma-lies in v = v ( p ) might merely result from ignoring in oursimulations the earlier inspiral phase. The remarkableoutcome of this test, however, is that the oscillatory be-havior in the kick as a function of eccentricity is more pronounced in the long sequence. The solid blue curve inthe bottom-right panel of Fig. 2 displays significantly morerapid oscillations in the eccentricity regime . (cid:46) e t (cid:46) . as compared to the shorter inspiral sequences. This oscil-latory behavior, and the apparent increase in the numberof oscillations as we increase the initial separation of theBHs, is the second of our results.We next attempt to gain insight into the origin ofthis behavior. For this purpose, we have computed themultipolar contributions to the total kick according toEqs. (13)-(17). The resulting velocities are displayed inFig. 2 by the additional dashed, dotted and dash-dottedcurves. Here, the curves labeled v (cid:96) =2 have been computedfrom the (cid:96) = 2 modes of Ψ ( ψ ,m only) in Eqs. (12)-(13).We computed this additional contribution (red dashedcurves in the figure) to determine whether the oscillatorybehavior is present also in the pure quadrupole signal.The answer is yes: the oscillations are clearly perceptiblein v (cid:96) =2 , even though they are a bit milder than in the total kick v tot . Considering all (cumulative) multipolarcontributions shown in Fig. 2, we notice the followingbehavior:(1) The oscillatory dependence of the kick on eccentric-ity is present at any level of truncating the multi-polar contributions in the cumulative sum (17).(2) The partial sum of the kick up to ˜ (cid:96) = 4 barelydiffers from the total kick, indicating that higher-order overlap terms do not significantly contributeto the kick.(3) The higher-order contributions ˜ (cid:96) > to the cumula-tive kick (17) systematically decrease the kick, coun-teracting the pure quadrupole contribution v (cid:96) =2 .In short, we have not identified any specific multipolesdominating the variation in the kick function v = v tot ( p ) .In our search for an explanation, we turn next to theinfall direction of the BH binary just before merger. Awell-known feature of the superkicks generated in theinspiral of BHs with opposite spins S = − S pointingin the orbital plane is the sinusoidal variation with theinitial azimuthal angle of the spin vectors: cf. Fig. 4 inRef. [89]. The initial orientation of the spins can, alterna-tively, be interpreted as a measure for the angle betweenthe in-plane spin components and the BH binary’s in-fall direction at merger [90]. The superkick is thereforecommonly determined by simulating otherwise identicalBH binary configurations for different values of this angleand fitting the resulting data with a cosine function: seee.g. Sec. III A in Ref. [1]. For the eccentric, nonspinningBH binaries considered in this work, it is the initial apsis(either a periapsis or an apoapsis) that defines a referencedirection. Unfortunately, neither the apsis nor a “binaryinfall direction” are rigorously defined quantities in thestrong-field regime of general relativity, and we considerinstead the orientation of the final kick relative to the x axis, defined by ˜ ϑ = arg( v x + i v y ) . (21)For convenience, we define ϑ = ˜ ϑ + 2 nπ, (22)where n ≥ is chosen minimally for each configurationin order to obtain ϑ as a monotonic function of the ini-tial tangential momentum p for each sequence. We willinterchangeably refer to ϑ and ˜ ϑ as the angle of the kick.Since all of our simulations start with the BHs locatedon the x axis with purely tangential initial momentum p = (0 , ± p, (Fig. 1), the x direction can be regardedas the initial direction of the apoapsis. If we furthermoreinterpret the gravitational recoil to be predominantly gen-erated by the excess beaming of the GWs in the directionof the smaller and faster BH (see Fig. 3 in Ref. [91]) dur-ing the short merger phase, the kick direction can serveas an approximate measure for the infall direction of thebinary. π π π π π π π π π π π FIG. 4. Plots involving the angle of the kick ϑ for all sequences. In the left panel, we plot the BH recoil velocity v against ϑ . Inthe right panel, we plot the location of the local extrema ϑ extrema of the left panel against the index of the extrema k countingrightwards from the global maximum on the left. We can test this prediction by computing the kick mag-nitude as a function of the angle ϑ ; if correct, we wouldexpect a periodic variation with a period close to π . Wedo not expect an exact π periodicity because the relevantperiapsis (or apoapsis) direction should be the last onebefore merger, and will shift away from the x axis duringthe inspiral due to apsidal precession – the BH analog ofMercury’s perihelion precession around the Sun. Morespecifically, we would expect deviations from a π period-icity to be more pronounced for longer inspirals, i.e. lowereccentricity and/or larger initial separations, but onlymildly dependent on the mass ratio q . Quite remarkably,all these features are borne out by the functions v = v ( ϑ ) displayed for our four sequences in the left panel of Fig. 4and the location of the extrema in this plot shown in theright panel of Fig. 4. For all sequences we observe thesame approximate π periodicity, with deviations fromthis value increasing at larger ϑ , i.e. for longer inspirals.Note also that ϑ = − π in the head-on limit, as expectedfor our initial configurations, that start with the heavierBH located on the positive x axis.While short of a rigorous proof, this result providesconsiderable evidence in favor of interpreting the oscil-latory dependence of the kick on the eccentricity as aconsequence of the corresponding variation in the infalldirection as measured relative to the last apoapsis (orperiapsis) of the eccentric binary. This interpretation alsoexplains why the longer sequence lq1:2 exhibits moreoscillations than the shorter sequences sq1:3 , sq1:2 and sq2:3 . Let us consider for this purpose two binary config-urations that only differ by a tiny amount of eccentricity δe . The longer the inspiral phase, the more time thesetwo binaries have to build up a considerable phase differ-ence and, hence, a different kick and merger GW signal.Note the potentially dramatic consequences of this be-havior for the GW emission from eccentric binaries overastrophysical timescales. For long astrophysical inspiralsretaining some eccentricity near merger, the kick and GW merger signal should exhibit critical dependence on theeccentricity. In terms of our Fig. 2, the function v = v ( e t ) would display a huge number of oscillations rather thanthe handful observed in our case, and the resulting curvewould look like a “band” rather than a single line. Withinthe band, a very small change δe t in eccentricity canproduce a finite change in the kick and merger waveform.As indicated by our analysis of the multipolar contribu-tions to the total recoil, the variations in the GW signalare of complex nature. We defer a more comprehensiveanalysis of the GW pattern to future work, but merelyillustrate with an example the type of variations thatare encountered. For this purpose, we show in Fig. 5the ( (cid:96), m ) = (2 , and (3 , multipoles of the GW sig-nal around merger for the configurations lq1:2-p0537 and lq1:2-p0567 , corresponding to a local minimum andmaximum in the kick, respectively: cf. the bottom-rightpanel of Fig. 2. In Fig. 5, the time has been shifted suchthat ∆ t = 0 corresponds to the first occurrence of a com-mon apparent horizon. The main difference perceptiblein the figure is the relative phase shift of the (3,3) moderelative to the dominant quadrupole (2,2). For the case p = 0 . M with maximal kick, the global peaks of bothmultipoles are aligned, whereas for p = 0 . M with min-imal kick, the global peak of the (2 , mode coincideswith a minimum in ( (cid:96), m ) = (3 , . We have made similarobservations for other pairs of modes such as (2 , and (2 , , and find these pairs to dominate the oscillatoryvariation in the multipolar series expansion (17). IV. CONCLUSIONS
In this paper we have studied the gravitational recoiland GW emission of sequences of nonspinning BH binarieswith mass ratios q = 2 / , / and / , and eccentricityvarying from the quasicircular to the head-on limit. Forthis purpose we have evolved configurations with - - - - -
50 0 50 100 - - - FIG. 5. The real part of the ( (cid:96), m ) = (2 , and (3 , modesof Ψ are shown as functions of time for the two binaries ofsequence lq1:2 with p/M = 0 . and with p/M = 0 . ,resulting in kick velocities of v = 128 km / s and
173 km / s ,respectively. the GRChombo and
Lean codes. Both codes yieldconvergent results for the recoil with a total error budgetof - and exhibit excellent agreement, well withinthis uncertainty estimate, for a verification configurationsimulated with both codes. In order to estimate theimpact of variations in the overall length of the inspirals,we have evolved two sequences for the case q = 1 / which complete about 3 and 6 orbits, respectively, in thequasicircular limit.The findings of our study are summarized as follows.(i) For all sequences, the total recoil reaches a globalmaximum for moderate eccentricities e ∼ . . As inthe case of the enhancement of superkicks studiedin Ref. [1], the maximum kick is enhanced by upto about 25 % relative to the value obtained forquasicircular configurations.(ii) Besides this global maximum, we observe an oscil-latory dependence of the kick v as a function ofeccentricity, with several local minima and maximain the function v = v ( e ) . Appropriate nonzero val-ues of the eccentricity can lead to a reduction of thekick by ∼
10 % relative to the quasicircular valueinstead of an increase. By splitting the kick intoseparate multipolar contributions, we notice thatthis oscillatory dependence is already present, al-beit in slightly weaker form, when we consider onlyquadrupole terms in the series expansion (12). Fur-ther contributions involving (cid:96) ≥ multipoles tendto decrease the overall kick and mildly enhance theoscillatory variation: see Fig. 2.(iii) We interpret this oscillatory variation in the kickas a consequence of changes in the angle between the infall direction at merger and the apoapsis (orperiapsis) direction. In the absence of rigorous defi-nitions for either of these directions, we approximatethis angular variation by considering the directionof the final kick and the x axis, assuming that theformer is related via relativistic GW beaming to theinfall direction and by taking into account that ourBHs start on the x axis with zero radial momentum.Displayed as a function of this angle, the kick dis-plays the expected periodic behavior with a periodclose to, but mildly deviating from π , presumablydue to periapsis precession.(iv) We have explored the dependence of this oscillatorybehavior of the recoil by simulating an additional se-quence of eccentric binaries with mass ratio q = 1 / ,but less negative binding energy, corresponding toabout 6 orbits in the quasicircular limit. We findthe oscillations in v = v ( e ) to be more pronouncedand numerous than in the shorter sequence. Weattribute this feature to the longer available timewindow during which otherwise identical binarieswith tiny differences in the initial eccentricity buildup a phase difference prior to merger. This ob-servation raises the intriguing possibility that thetotal recoil depends highly sensitively on the initialeccentricity.(v) The variations in the kick velocity are accompaniedby relative time shifts in the peak amplitudes ofsubdominant multipoles relative to the peaks of the(2,2) mode: cf. Fig. 5. For configurations with large(small) kick, the peak amplitude of subdominantmultipoles tends to be aligned (misaligned) withthe quadrupole peak.Our findings point to a variety of future investigations.While our simulations indicate an increased sensitivity ofthe GW merger signal to the initial eccentricity for largerinitial separations (i.e. longer inspirals), it is not clearhow this will be affected by the circularizing nature ofGW emission. In this context, it will also be importantto analyze in more quantitative terms the differences inthe GW signals and possible implications for parameterinference in GW observations. A thorough investigationof long eccentric inspirals on astrophysical timescales willlikely require PN methods and may benefit greatly froma multi-time scale analysis in phase space, as appliedto spin-precessing BH binaries in Refs. [92, 93] or tothe dynamics of binary systems in external gravitationalbackground potentials in Refs. [94, 95]. If there is a singleconclusion to draw from the results of this work, it isthe surprisingly rich phenomenology of the GW signals ofeccentric compact binaries – even in the absence of spins –which merits as much as it requires further investigation.0 ACKNOWLEDGMENTS
We thank Michalis Agathos, Vishal Baibhav, VitorCardoso, Thomas Helfer and Nicholas Speeney for use-ful discussions. We also thank Chris Moore, CarlosLousto and Juan Calderón Bustillo for helpful commentson this manuscript. M.R. thanks the GRChombo col-laboration [59] for their code development, and partic-ularly Katy Clough and Tiago França. M.R. is sup-ported by a Science and Technology Facilities Council(STFC) studentship. U.S. is supported by the EuropeanUnion’s H2020 ERC Consolidator Grant “Matter andstrong-field gravity: new frontiers in Einstein’s theory”Grant No. MaGRaTh–646597, and the STFC Consolida-tor Grant No. ST/P000673/1. E.B. is supported by NSFGrant No. PHY-1912550, NSF Grant No. AST-2006538,NASA ATP Grant No. 17-ATP17-0225 and NASA ATPGrant No. 19-ATP19-0051. This work has received fund-ing from the European Union’s Horizon 2020 research andinnovation programme under the Marie Skłodowska-CurieGrant No. 690904. This work was supported by the GW-verse COST Action CA16104, “Black holes, gravitationalwaves and fundamental physics”. Computational work wasperformed on the San Diego Supercomputer Center Cometand Texas Advanced Computing Center (TACC) Stam-pede2 clusters at the University of California San Diegoand the University of Texas at Austin (UT Austin) re-spectively through NSF-XSEDE Grant No. PHY-090003;the Cambridge Service for Data Driven Discovery (CSD3)system at the University of Cambridge through STFC cap-ital Grants No. ST/P002307/1 and No. ST/R002452/1,and STFC operations Grant No. ST/R00689X/1; theTACC Frontera cluster at UT Austin; and the JUWELScluster at GCS@FZJ, Germany through PRACE GrantNo. 2020225359.
Appendix A: Numerical Accuracy
As in Ref. [1], the uncertainty in our numerical resultsfor the recoil velocities has two predominant contributions:the discretization error and the finite extraction radii forthe Weyl scalar Ψ .To estimate the uncertainty arising from the latter, wehave selected a representative sample of the simulationsfrom each sequence and extrapolated the cumulative ra-diated momentum to infinity from about six extractionradii in the range r ex / ≤ r ex using a Taylor series in /r as in Ref. [96]. We report the results from the finiteextraction radii given in Table I and estimate the errorby comparing with the linear-order extrapolation. Forboth codes, we estimate the contribution from this erroris about 2% for all sequences.In order to estimate the error contribution from finitedifferencing and verify that our codes give consistentresults, we have performed simulations of sq1:2-p0100 (the binary in sequence sq1:2 with p/M = 0 . ) with bothcodes. We discuss the analyses of the convergence of each Label
L N t R b h L /M taggingR1 7 320 0.012 . M i . M i . M i GRChombo simula-tions. As explained in Sec. II A 1 and in Appendix B, the totalnumber of refinement levels is L + 1 , the number of cells alongeach dimension on the coarsest level is N , t R is the regriddingthreshold value, b is the BH tagging buffer parameter thatwe set proportional to the mass M i ( i = 1 , ) of the nearestBH for all configurations except R4, and h L denotes the gridspacing. codes separately before comparing. GRChombo convergence
For
GRChombo , we have performed the simulations of sq1:2-p0100 with resolutions h L = 3 M / , M / and M / , and we refer to the configurations correspond-ing to these resolutions as R1, R2 and R3, respectively.The full grid configurations are given in Table II and theresults of this analysis are shown in the left panel of Fig. 6.Around merger, at ( t − r ex ) /M ∼ , our results exhibitmild overconvergence in the top-left panel of Fig. 6. Theimportant results for our analysis in Fig. 2 below, how-ever, are the final kick values after the merged BH hassettled down. As can be seen from the inset, the conver-gence here is close to fifth order. From our convergenceanalysis, the difference between the result obtained fromthe R1 simulation and the more conservative fourth-orderRichardson extrapolated result leads to an estimate of thediscretization error of about . A similar error estimateis also obtained for the radiated energy, E rad . From ex-perience, we have found smaller values for the mass ratio, q < , more challenging to simulate accurately than largervalues, and we therefore feel justified in using this errorestimate (for a q = 1 / configuration) as a conservativeestimate for the error in the sq2:3 sequence simulations( q = 2 / ). We therefore used the R1 grid configurationfor this sequence with l max1 = l max2 = L = 7 (both BHs arecovered by the finest level: see Appendix B for details).For the sq1:3 simulations, we used the R4 grid con-figuration (see Table II) with l max1 = L = 7 and l max2 = L − (the larger BH is not covered by the finestlevel: see Appendix B for details). This corresponds toa resolution of h L = 3 M / . We performed a separateconvergence analysis of sq1:3-p0089 , which led to anestimated 1% discretization error.Combining both the finite extraction radius and dis-cretization errors, our estimate for the total error budgetof the GRChombo simulations is about .1 - - - - - - - - - - - - - FIG. 6. For each code, we show convergence plots for the accumulated linear momentum radiated from sq1:2-p0100 by plottingthe BH recoil velocity in the bottom panels. The Richardson extrapolated curve, v Rich4 , assuming fourth-order convergence,is also shown in the bottom panel. The grid configurations are given in Table II for
GRChombo , and in Table III for
Lean .The top panel shows the difference between the configurations along with rescalings corresponding to fourth- and fifth-orderconvergence. The inset shows a magnification of the right side of the plot: the final value of the recoil velocity is what we showin Fig. 2 below.Label
L l F R R L h L /M S1 7 4 384 1 1/20S2 7 4 384 1 1/24S3 7 4 384 1 1/32S4 7 4 384 1 1/28TABLE III. Grid configurations used for
Lean simulations.As explained in Sec. II A 2, the total number of refinementlevels is L + 1 , the number of fixed refinement levels is l F + 1 , R is the half-length of the outer grid, R L is the half-lengthof one cubic component of the innermost grid, and h L is thegrid spacing on the finest level. Lean convergence
With
Lean , we have simulated sq1:2-p0100 with res-olutions h L = M / , M / and M / . We refer tothese grid configurations as S1, S2 and S3, respectively (cf.Table III). The right panel of Fig. 6 shows convergencebetween fourth and fifth order. For simulations in sq1:2 ,we used the S2 grid configuration. From the convergenceanalysis, the difference between the result obtained fromthe S2 simulation and the fourth-order Richardson extrap-olation leads to an estimate of the discretization error ofabout . . For the lq1:2 simulations, we have undertaken a sepa-rate convergence analysis of lq1:2-p0086 using the samegrid setup as in Table III, but using higher resolutions h L /M = 1 / , / and / . We observe convergenceclose to fourth order and obtain an error estimate of from the Richardson extrapolated kick for the mediumresolution h L /M = 1 / .In summary, the Lean simulations of sequence sq1:2 are performed with resolution grid S2 of Table III and anerror budget of . , % , and those of sequence lq1:2 withgrid S4 of Table III and an error budget of .
3. Comparison between
GRChombo and
Lean
A comparison of the recoil velocity computed from
GRChombo and
Lean simulations of sq1:2-p0100 withthe grid configurations used for sq2:3 and sq1:2 runsrespectively is shown in the top panel of Fig. 7. Theeccentricity estimate for this system is e t = 0 . . Wehave chosen this configuration for two reasons. First, todetermine appropriate resolutions, we had to calibrateour codes’ accuracy at the start of our exploration, whichwe began in the regime of mild eccentricities to acquirean intuitive understanding of their behavior. Second,configurations with mild eccentricity have a longer inspiralphase than highly eccentric ones, and therefore imposea stronger requirement on phase accuracy. A mildly2 - - - FIG. 7. Comparison between
GRChombo and
Lean for theaccumulated linear momentum radiated in GWs in simulationsof sq1:2-p0100 with e t = 0 . . We compare the BH recoilvelocity (top panel) and the corresponding plus-polarized (cid:96) = m = 2 strain amplitude (bottom panel). eccentric binary is therefore ideally suited to obtain aconservative estimate of the numerical accuracy, which isrepresentative across the targeted parameter space.The final recoil velocities obtained for this configurationwith our two codes differ by about 2%, which is wellwithin the error budget of each code. We also show thequadrupole contribution h +2 , to the ‘+’ polarization straindefined by [97] h + (cid:96),m ( t, r ) − i h × (cid:96),m ( t, r )= a (cid:96),m + b (cid:96),m t + (cid:90) t d t (cid:48) (cid:90) t (cid:48) d t (cid:48)(cid:48) ψ (cid:96),m ( t (cid:48)(cid:48) , r ) , (A1)where the constants a (cid:96),m and b (cid:96),m are chosen to minimizelinear drift, in the bottom panel of the figure, to betterillustrate the agreement between the codes for these gridconfigurations.In Fig. 6, the differences between the results of dif-ferent resolutions with Lean are greater than that of
GRChombo . However, we found that
Lean entered theconvergent regime at lower resolutions than GRChombo.This is compatible with the observations of Ref. [62] thathigher resolutions were required for convergence withCCZ4 compared to BSSNOK.
Appendix B:
GRChombo tagging criterion
As explained in Sec. II A 1, the regridding is controlledby the tagging of cells for refinement in the Berger-Rigoutsos algorithm [65], with cells being tagged if thetagging criterion C exceeds the specfied threshold value t R as given in Table II. For this work, we use the taggingcriterion C = (cid:40) , if l ≥ l maxBH and r BH < ( M BH + b ) , max( C χ , C punc , C ex ) , otherwise , (B1)where l maxBH is a specifiable maximum level parameter foreach BH (so that it is not unnecessarily over resolved), r BH is the coordinate distance to the puncture, M BH isthe mass of the corresponding BH, b is a buffer parameter,and C χ , C punc , and C ex are given as follows:(i) C χ tags regions in which the gradients of the con-formal factor, χ , become steep. It is given by C χ = h l (cid:115)(cid:88) i,j ( ∂ i ∂ j χ ) , (B2)where h l is the grid spacing on refinement level l .(ii) C punc tags within spheres around each puncture inorder to ensure the horizon is suitably well-resolved.It is given by C punc = (cid:40) , if r BH < ( M BH + b )2 max( l maxBH − l − , , , otherwise . (B3)(iii) C ex ensures each sphere on which we extract theWeyl scalar, Ψ , is suitably well-resolved. It is givenby C ex = (cid:40) , if r < . r ex and l < l ex , , otherwise , (B4)where r is the coordinate distance to the centerof mass, r = r ex gives the location of the extrac-tion sphere, and l ex is a specifiable extraction levelparameter for each sphere.We also used this tagging criterion with the replacement r BH → max( x BH , y BH , z BH ) , where, e.g. x BH is the dis-tance to the puncture in the x direction. We refer to thisas “box” tagging and the original as “spherical” tagging.Naively, one might hope that C χ is sufficient to ensuresuitable refinement around the BHs, since the gradientsof χ become increasingly steep close to the punctures.However we found empirically that, without C punc , thehorizons are perturbed significantly by the refinementboundaries, leading to lower accuracy.3 [1] U. Sperhake, R. Rosca-Mead, D. Gerosa, and E. Berti,Phys. Rev. D , 024044 (2020), arXiv:1910.01598 [gr-qc].[2] B. P. Abbott et al. , Phys. Rev. Lett. , 061102 (2016),arXiv:1602.03837 [gr-qc].[3] P. R. Saulson, Gen. Rel. Grav. , 3289 (2011).[4] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.Rev. D , 041501 (2006), arXiv:gr-qc/0604012.[5] U. Sperhake, V. Cardoso, F. Pretorius, E. Berti, T. Hin-derer, and N. Yunes, Phys. Rev. Lett. , 131102 (2009),arXiv:0907.1252 [gr-qc].[6] W. B. Bonnor, M. A. Rotenberg, and L. Rosenhead, Proc.R. Soc. Lond. A Math. Phys. Sci. , 109 (1961).[7] A. Peres, Phys. Rev. , 2471 (1962).[8] J. D. Bekenstein, Astrophys. J. , 657 (1973).[9] D. Gerosa and C. J. Moore, Phys. Rev. Lett. , 011101(2016), arXiv:1606.04226 [gr-qc].[10] J. Calderón Bustillo, J. A. Clark, P. Laguna, andD. Shoemaker, Phys. Rev. Lett. , 191102 (2018),arXiv:1806.11160 [gr-qc].[11] C. O. Lousto and J. Healy, Phys. Rev. D , 104039(2019), arXiv:1908.04382 [gr-qc].[12] V. Varma, M. Isi, and S. Biscoveanu, Phys. Rev. Lett. , 101104 (2020), arXiv:2002.00296 [gr-qc].[13] M. J. Fitchett, Mon. Not. R. Astron. Soc. , 1049(1983).[14] L. Blanchet, M. S. Qusailah, and C. M. Will, Astrophys.J. , 508 (2005), arXiv:astro-ph/0507692.[15] S. A. Hughes, M. Favata, and D. E. Holz, in Conferenceon Growing Black Holes: Accretion in a CosmologicalContext (2004) arXiv:astro-ph/0408492.[16] T. Damour and A. Gopakumar, Phys. Rev. D , 124006(2006), arXiv:gr-qc/0602117.[17] C. F. Sopuerta, N. Yunes, and P. Laguna, Phys. Rev. D ,124010 (2006), [Erratum: Phys.Rev.D 75, 069903 (2007),Erratum: Phys.Rev.D 78, 049901 (2008)], arXiv:astro-ph/0608600.[18] C. F. Sopuerta, N. Yunes, and P. Laguna, Astrophys. J.Lett. , L9 (2007), arXiv:astro-ph/0611110.[19] A. Le Tiec, L. Blanchet, and C. M. Will, Class. Quant.Grav. , 012001 (2010), arXiv:0910.4594 [gr-qc].[20] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, J. R.van Meter, and M. Miller, Astrophys. J. Lett. , L93(2006), arXiv:astro-ph/0603204.[21] J. A. Gonzalez, U. Sperhake, B. Bruegmann, M. Hannam,and S. Husa, Phys. Rev. Lett. , 091101 (2007), arXiv:gr-qc/0610154.[22] F. Herrmann, I. Hinder, D. Shoemaker, and P. Laguna,Class. Quant. Grav. , S33 (2007).[23] J. A. González, M. D. Hannam, U. Sperhake, B. Brüg-mann, and S. Husa, Phys. Rev. Lett. , 231101 (2007),arXiv:gr-qc/0702052.[24] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Mer-ritt, Phys. Rev. Lett. , 231102 (2007), arXiv:gr-qc/0702133.[25] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Mer-ritt, Astrophys. J. Lett. , L5 (2007), arXiv:gr-qc/0701164 [gr-qc].[26] C. O. Lousto and Y. Zlochower, Phys. Rev. Lett. ,231102 (2011), arXiv:1108.2009 [gr-qc].[27] J. D. Schnittman and A. Buonanno, Astrophys. J. Lett. , L63 (2007), arXiv:astro-ph/0702641 [ASTRO-PH].[28] M. Dotti, M. Volonteri, A. Perego, M. Colpi,M. Ruszkowski, and F. Haardt, Mon. Not. Roy. Astron.Soc. , 682 (2010), arXiv:0910.5729 [astro-ph.HE].[29] M. Kesden, U. Sperhake, and E. Berti, Astrophys. J. ,1006 (2010), arXiv:1003.4993 [astro-ph.CO].[30] C. O. Lousto, Y. Zlochower, M. Dotti, and M. Volonteri,Phys. Rev. D85 , 084015 (2012), arXiv:1201.1923 [gr-qc].[31] E. Berti, M. Kesden, and U. Sperhake, Phys. Rev.
D85 ,124049 (2012), arXiv:1203.2920 [astro-ph.HE].[32] S. Komossa, Adv. Astron. , 364973 (2012),arXiv:1202.1977 [astro-ph].[33] M. Colpi, Space Sci. Rev. , 189 (2014),arXiv:1407.3102 [astro-ph.GA].[34] L. Blecha, D. Sijacki, L. Z. Kelley, P. Torrey, M. Vogels-berger, D. Nelson, V. Springel, G. Snyder, and L. Hern-quist, Mon. Not. Roy. Astron. Soc. , 961 (2016),arXiv:1508.01524 [astro-ph.GA].[35] L. Barack et al. , Class. Quant. Grav. , 143001 (2019),arXiv:1806.05195 [gr-qc].[36] D. Merritt, M. Milosavljevic, M. Favata, S. A. Hughes,and D. E. Holz, Astrophys. J. Lett. , L9 (2004),arXiv:astro-ph/0402057.[37] M. J. Benacquista and J. M. Downing, Living Rev. Rel. , 4 (2013), arXiv:1110.4423 [astro-ph.SR].[38] J. Morawski, M. Giersz, A. Askar, and K. Belczynski, Mon.Not. Roy. Astron. Soc. , 2168 (2018), arXiv:1802.01192[astro-ph.GA].[39] Z. Haiman, Astrophys. J. , 36 (2004), arXiv:astro-ph/0404196.[40] M. J. Graham et al. , Phys. Rev. Lett. , 251102 (2020),arXiv:2006.14122 [astro-ph.HE].[41] H.-Y. Chen, C.-J. Haster, S. Vitale, W. M. Farr, andM. Isi, arXiv:2009.14057 [astro-ph.CO] (2020).[42] J. Samsing and E. Ramirez-Ruiz, Astrophys. J. Lett. ,L14 (2017), arXiv:1703.09703 [astro-ph.HE].[43] J. Samsing, Phys. Rev. D , 103014 (2018),arXiv:1711.07452 [astro-ph.HE].[44] J. Samsing, A. Askar, and M. Giersz, Astrophys. J. ,124 (2018), arXiv:1712.06186 [astro-ph.HE].[45] C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee, K. Kre-mer, F. A. Rasio, J. Samsing, C. S. Ye, and M. Zevin,Phys. Rev. D , 123005 (2018), arXiv:1811.04926 [astro-ph.HE].[46] J. Samsing, I. Bartos, D. D’Orazio, Z. Haiman, B. Koc-sis, N. Leigh, B. Liu, M. Pessah, and H. Tagawa,arXiv:2010.09765 [astro-ph.HE] (2020).[47] H. Tagawa, B. Kocsis, Z. Haiman, I. Bartos, K. Omukai,and J. Samsing, arXiv:2010.10526 [astro-ph.HE] (2020).[48] B. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. , 149 (2019), arXiv:1907.09384 [astro-ph.HE].[49] R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. , 101102 (2020), arXiv:2009.01075 [gr-qc].[50] I. M. Romero-Shaw, P. D. Lasky, E. Thrane, and J. C.Bustillo, arXiv:2009.04771 [astro-ph.HE] (2020).[51] V. Gayathri, J. Healy, J. Lange, B. O’Brien, M. Szczep-anczyk, I. Bartos, M. Campanelli, S. Klimenko, C. Lousto,and R. O’Shaughnessy, arXiv:2009.05461 [astro-ph.HE](2020).[52] A. Nishizawa, E. Berti, A. Klein, and A. Sesana, Phys.Rev. D94 , 064020 (2016), arXiv:1605.01341 [gr-qc]. [53] K. Breivik, C. L. Rodriguez, S. L. Larson, V. Kalogera,and F. A. Rasio, Astrophys. J. Lett. , L18 (2016),arXiv:1606.09558 [astro-ph.GA].[54] A. Nishizawa, A. Sesana, E. Berti, and A. Klein, Mon.Not. Roy. Astron. Soc. , 4375 (2017), arXiv:1606.09295[astro-ph.HE].[55] C. Roedig and A. Sesana, Gravitational waves. Numericalrelativity - data analysis. Proceedings, 9th Edoardo AmaldiConference, Amaldi 9, and meeting, NRDA 2011, Cardiff,UK, July 10-15, 2011 , J. Phys. Conf. Ser. , 012035(2012), arXiv:1111.3742 [astro-ph.CO].[56] M. Bonetti, A. Sesana, F. Haardt, E. Barausse, andM. Colpi, Mon. Not. Roy. Astron. Soc. , 4044 (2019),arXiv:1812.01011 [astro-ph.GA].[57] E. Huerta et al. , Phys. Rev. D , 064003 (2019),arXiv:1901.07038 [gr-qc].[58] K. Clough, P. Figueras, H. Finkel, M. Kunesch, E. A.Lim, and S. Tunyasuvunakool, Class. Quant. Grav. ,24 (2015), arXiv:1503.03436 [gr-qc].[59] .[60] U. Sperhake, Phys. Rev. D , 104015 (2007), arXiv:gr-qc/0606079.[61] S. Husa, J. A. Gonzalez, M. Hannam, B. Bruegmann,and U. Sperhake, Class. Quant. Grav. , 105006 (2008),arXiv:0706.0740 [gr-qc].[62] D. Alic, C. Bona-Casas, C. Bona, L. Rezzolla, andC. Palenzuela, Phys. Rev. D , 064040 (2012),arXiv:1106.2254 [gr-qc].[63] D. Alic, W. Kastaun, and L. Rezzolla, Phys. Rev. D ,064049 (2013), arXiv:1307.7391 [gr-qc].[64] M. Adams et al. , Chombo Software Package for AMRApplications - Design Document , Tech. Rep. LBNL-6616E(Lawrence Berkeley National Laboratory, 2019).[65] M. Berger and I. Rigoutsos, IEEE Trans. Sys. Man &Cyber. , 1278 (1991).[66] T. Goodale, G. Allen, G. Lanfermann, J. Massó, T. Radke,E. Seidel, and J. Shalf, in Vector and Parallel Processing– VECPAR’2002, 5th International Conference, LectureNotes in Computer Science (Springer, Berlin, 2003).[67] T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor.Phys. Suppl. , 1 (1987).[68] M. Shibata and T. Nakamura, Phys. Rev. D , 5428(1995).[69] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D ,024007 (1999), arXiv:gr-qc/9810065.[70] M. Campanelli, C. Lousto, P. Marronetti, and Y. Zlo-chower, Phys. Rev. Lett. , 111101 (2006), arXiv:gr-qc/0511048.[71] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, andJ. van Meter, Phys. Rev. Lett. , 111102 (2006), arXiv:gr-qc/0511103.[72] E. Schnetter, S. H. Hawley, and I. Hawke, Class. Quant.Grav. , 1465 (2004), arXiv:gr-qc/0310042.[73] J. Thornburg, Phys. Rev. D , 4899 (1996), gr-qc/9508014.[74] J. Thornburg, Class. Quant. Grav. , 743 (2004), arXiv:gr-qc/0306056.[75] S. Brandt and B. Bruegmann, Phys. Rev. Lett. , 3606(1997), arXiv:gr-qc/9703066.[76] J. M. Bowen and J. York, James W., Phys. Rev. D ,2047 (1980).[77] M. Ansorg, B. Bruegmann, and W. Tichy, Phys. Rev. D , 064011 (2004), arXiv:gr-qc/0404056.[78] V. Paschalidis, Z. B. Etienne, R. Gold, and S. L. Shapiro,arXiv:1304.0457 [gr-qc] (2013).[79] U. Sperhake, E. Berti, V. Cardoso, J. A. Gonzalez,B. Bruegmann, and M. Ansorg, Phys. Rev. D , 064069(2008), arXiv:0710.3823 [gr-qc].[80] R. Arnowitt, S. Deser, and C. W. Misner, in Gravitationan introduction to current research , edited by L. Witten(John Wiley, New York, 1962) pp. 227–265, arXiv:gr-qc/0405109.[81] B. Bruegmann, J. A. Gonzalez, M. Hannam, S. Husa,U. Sperhake, and W. Tichy, Phys. Rev. D , 024027(2008), arXiv:gr-qc/0610128.[82] J. W. York, Jr., Initial data for collisions of black holesand other gravitational miscellany (Cambridge UniversityPress, 1989) p. 89.[83] N. Loutrel, S. Liebersbach, N. Yunes, and N. Cornish,Class. Quant. Grav. , 025004 (2019), arXiv:1810.03521[gr-qc].[84] R.-M. Memmesheimer, A. Gopakumar, and G. Schaefer,Phys. Rev. D , 104011 (2004), arXiv:gr-qc/0407049.[85] M. Campanelli and C. O. Lousto, Phys. Rev. D , 124022(1999), gr-qc/9811019.[86] C. O. Lousto and Y. Zlochower, Phys. Rev. D , 041502(2007), gr-qc/0703061, arXiv:gr-qc/0703061.[87] M. Ruiz, R. Takahashi, M. Alcubierre, and D. Nunez,Gen. Rel. Grav. , 2467 (2008), arXiv:0707.4654 [gr-qc].[88] J. Healy, C. O. Lousto, and Y. Zlochower, Phys. Rev. D , 024031 (2017), arXiv:1705.07034 [gr-qc].[89] B. Brügmann, J. A. González, M. D. Hannam, S. Husa,and U. Sperhake, Phys. Rev. D , 124047 (2008),arXiv:0707.0135 [gr-qc].[90] C. O. Lousto, H. Nakano, Y. Zlochower, and M. Cam-panelli, Phys. Rev. D , 084023 (2010), [Erratum:Phys.Rev.D 82, 129902 (2010)], arXiv:0910.3197 [gr-qc].[91] A. G. Wiseman, Phys. Rev. D , 1517 (1992).[92] M. Kesden, D. Gerosa, R. O’Shaughnessy, E. Berti,and U. Sperhake, Phys. Rev. Lett. , 081103 (2015),arXiv:1411.0674 [gr-qc].[93] D. Gerosa, M. Kesden, U. Sperhake, E. Berti, andR. O’Shaughnessy, Phys. Rev. D , 064016 (2015),arXiv:1506.03492 [gr-qc].[94] C. Hamilton and R. R. Rafikov, Mon. Not. Roy. Astron.Soc. , 5489 (2019), arXiv:1902.01344.[95] C. Hamilton and R. R. Rafikov, Mon. Not. Roy. Astron.Soc. , 5512 (2019), arXiv:1902.01345.[96] U. Sperhake, B. Bruegmann, D. Muller, and C. Sopuerta,Class. Quant. Grav. , 134004 (2011).[97] N. T. Bishop and L. Rezzolla, Living Rev. Rel.19