Intermediate mass-ratio black hole binaries: Applicability of small mass-ratio perturbation theory

The dynamics of merging black hole binaries with small mass-ratios (SMR) can be expressed as a series in the symmetric mass-ratio, the post-adiabatic expansion. Using numerical relativity simulations, we compute the first three terms in this series for quasicircular non-spinning binaries. We recover the leading term predicted by SMR perturbation theory, and obtain a robust prediction of the 1-post-adiabatic term. The 2-post-adiabatic term is found to be small, indicating that the 1PA approximation may be quite accurate even for nearly comparable mass binaries. We also estimate the range of applicability for SMR and post-Newtonian series for non-spinning, quasi-circular inspirals.

The dynamics of merging black hole binaries with small mass-ratios (SMR) can be expressed as a series in the symmetric mass-ratio, the post-adiabatic expansion. Using numerical relativity simulations, we compute the first three terms in this series for quasicircular non-spinning binaries. We recover the leading term predicted by SMR perturbation theory, and obtain a robust prediction of the 1-post-adiabatic term. The 2-post-adiabatic term is found to be small, indicating that the 1PA approximation may be quite accurate even for nearly comparable mass binaries. We also estimate the range of applicability for SMR and post-Newtonian series for non-spinning, quasi-circular inspirals.
Inspiraling and merging binary black holes are the most numerous source of gravitational waves observed by the LIGO and Virgo gravitational wave detectors [1,2] and are one of the key science targets for the space-based LISA mission [3]. The mass-ratio q ≡ m 2 /m 1 ≤ 1 is one of the key parameters in the dynamics of these systems. For LIGO and Virgo observations published thus far, q has been close to unity [4][5][6] with GW190412 [7] the first system with clearly unequal masses (q ∼ 0.28).
In the future, observations of binaries with lower q are expected: With additional detections during continued operation of the detectors [8], the minimal observed q can only decrease. As third generation ground based detectors delve deeper into the low frequency regime they will become sensitive to the capture of stellar mass BHs by intermediate mass black holes (IMBHs) with massratios down to q ∼ 10 −3 [9]. The space-based gravitational wave observatory LISA will observe the mergers of massive black holes (MBHs) of millions of solar masses. While the majority of these are expected to have q 0.1, there could a significant tail of events down to q ∼ 0.01 [10,11]. LISA will also be sensitive to mergers of IMBHs with MBHs (q ∼ 10 −3 ) and extreme massratio inspirals (q ∼ 10 −5 ) as sensitive probes of black hole physics [12].
The modeling of inspiraling binaries at all mass-ratios is therefore of paramount importance for detection and analysis of GW sources. The three primary modeling approaches are post-Newtonian slow-velocity perturbation theory [13], numerical relativity (NR), i.e. direct numerical integration of the full non-linear Einstein equations [14], and small mass-ratio (SMR) perturbation theory [15]. Effective one body methods [16] provide a means to combine and resum information from all three approaches and also from newer developments like post-Minkowski expansions [17].
This article focuses on the SMR and NR approaches. The SMR approximation expands the dynamics of a coalescing binary in powers of q or the symmetric mass-ratio . At leading order, the secondary follows a geodesic in the background space-time generated by the primary. The impact of the secondary's mass on the dynamics can be included as an effective force term, the gravitational self-force. Over the last two decades much progress has been made in the developing this GSF formalism (see [18] for a recent review), but calculation of the next-to-leading order contribution to the orbital phasing has not yet been achieved. While the main motivation for SMR lies in extreme-mass-ratio inspirals, there is increasing evidence [19][20][21] that the SMR may be applicable even at comparable masses.
Numerical relativity directly solves the full non-linear Einstein equations [14]. The vast number of simulations performed to date are at comparable masses, with only very few simulations at q 0.1 (see, e.g. [22], but note [23,24] for simulations at q = 1/18 to q = 1/128). The limited coverage in q has two causes. First, the number of orbits the binary spends in the strong field region grows ∝ ν −1 . Second, because of the Courant-limit on the time-step of the numerical simulations, the number of time-steps per orbit increases ∝ q −1 . Combined, these effects cause an increase in computational cost at least quadratically in mass-ratio. The need for higher numerical resolution to resolve the ever smaller secondary (as q → 0), and to preserve phase-accuracy over the increasingly longer inspiral will increase computational cost further.
Given the expectation of binaries at all mass-ratios, an essential question is how to model intermediate mass-ratio binaries at small separation: post-Newtonian theory is not accurate close to merger owing to the high velocities; numerical relativity simulations are limited to large massratios, q 0.1; and the SMR approximation is presently only available at leading order in q, and thus may be inaccurate at intermediate mass-ratios. This letter investigates the existence of a "mass-ratio gap" where none of the modeling approaches is applicable. We analyse NR simulations at mass-ratios 0.1 ≤ q ≤ 1 computed with the SpEC code [25,26] and extract the first three terms in the arXiv:2006.12036v1 [gr-qc] 22 Jun 2020 SMR expansion of the orbital phasing. Analysing these terms, we conclude that SMR results at next-to-leading order (once they have been computed) can likely bridge the mass-ratio gap up to mass-ratios q large enough to be covered by numerical relativity. Our results, however, apply only to the simplest possible situation, black hole binaries without spin and without orbital eccentricity.
Methodology.-We use geometric units such that c = G = 1 and examine the orbital phase extracted from the gravitational radiation at future null infinity, Here h 22 is the spin-weight s = −2 spherical harmonic ( , m) = (2, 2) mode of the complex gravitational wave strain h = h + + ih × . The current work focuses on nonprecessing binaries where Eq. (1) is sufficient. The orbital phase can also be defined in terms of the Weyl scalar ψ 4 , and we use this as part of our error estimate.
Introducing the orbital frequency, we consider the orbital phase as a function of the orbital frequency, φ(Ω). In the SMR approximation, φ can be calculated by a two timescale expansion [27] leading to a power series in the mass-ratio, known as the postadiabatic (PA) expansion, Here, φ nPA are functions of M Ω, where M ≡ m 1 + m 2 is the total mass of the binary. Alternatively, one can consider φ nPA as functions of m 1 Ω and/or expand in q as the small parameter (cf. Fig. 3 below).
We use numerical relativity simulations from the SpECcode, which utilizes the quasi-local angular momentum formalism to control the black hole spins [49][50][51][52], iterative eccentricity reduction to achieve orbital eccentricities e 10 −4 [53,54], and solves the Einstein evolution equations in the generalized harmonic formulation [55][56][57][58] with constraint damping and minimally reflective outer boundary conditions [58][59][60] (see [25] for more details). Because of the use of spectral methods and a dual-frame approach [61] SpEC achieves very high accuracies even for long inspiral simulations that cover a comparatively large range in orbital frequencies. Gravitational radiation is extracted using the Regge-Wheeler-Zerilli formalism, extrapolated to future null infinity [25,62], and corrected for center of mass drifts [63].
For our investigation we utilize 55 NR simulations of nonspinning quasicircular inspirals from the public SXS catalog [22,64] with mass-ratios q ∈ [0.1, 1]. The initial orbital frequency is in the range M Ω ∼ 0.015 . . . 0.02. Simulations with smaller q tend to start at the higher frequencies, to achieve an computationally manageable overall duration of the simulations. The simulations are all available at at least two numerical resolutions.
The orbital phase φ NR (M Ω) is determined by locally fitting a low order polynomial in t to φ NR (t). The width of the fitting window is variable such that at low frequencies it encompasses several radial oscillations of any residual eccentricity in the simulations, while at larger frequencies it is small enough to avoid systematic bias due to the rapidly changing frequency. The constant of integration when integrating Eq.  NR simulations (with error bar), whereas the red line is the leading order SMR result computed by Eq. (4). The agreement between the two is quite remarkable, and is a first indication that the PA expansion of the phase in the mass-ratio is well-behaved for comparable mass-ratios.
At higher frequencies, M Ω 0.055, we find an apparently systematic deviation between NR and the SMR result. This deviation may arise from a breakdown of the PA expansion near the last stable orbit as the binary transitions from inspiral to plunge. Studies of this transition regime [65,66] lead to order ν −1/5 corrections to Eq. (3). Including such a term in our fit does indeed eliminate the systematic deviation at M Ω 0.055. However, the additional term is nearly degenerate with the 0PA and 1PA terms at low frequencies making it impossible to get robust numerical results for φ 1PA and higher. Therefore, we proceed in our analysis without such transition terms.
Given how well the numerically extracted 0PA term agrees with its SMR prediction, we henceforth set it to the SMR value when fitting for the higher order PA terms. Figure 2 shows the 1PA and 2PA term obtained from the NR simulations, together with the 0PA term already discussed in Fig. 1. The coefficients φ nPA are of comparable magnitude in the frequency range covered by our analysis, suggesting that the PA series is convergent at equal masses. Moreover, for frequencies M Ω 0.05, the 2PA coefficient is almost consistent with zero, i.e. the 0PA and 1PA terms already capture essentially all variation due to mass-ratio in the numerical data at these frequencies. In fact, "goodness-of-fit" indicators, such as the adjusted R 2 value, show only marginal improvements when adding terms to the fit beyond the 1PA coefficient.
The lower panel of Fig. 2 provides a different view on the importance of terms beyond φ 1PA : For each of the 55 NR simulations, this panel plots i.e. the contribution of all terms n ≥ 2 in Eq. (3), with overall ν scaling compensated. All R 2+ can be bounded independent of mass-ratio by an envelope function, consisting of the known 3.5PN terms of φ 2PA and a higher order polynomial in M Ω fitted by eye.
So far, we have expanded in symmetric mass-ratio ν, while scaling orbital frequencies by total mass M , cf. Eq. (3). One can also use the mass-ratio q = m 2 /m 1 as the small parameter, and/or scale orbital frequency by the large body's mass m 1 . This yields four variations, all of which agree at the leading 0PA-order. Figure 3 presents the results for the 1PA and 2PA contributions. In all four cases, the extracted 1PA and 2PA coefficients remain of similar magnitude, implying that the expansion is not dominated by higher order terms. However, the 2PA term is remarkably small only when expanding using the symmetric mass ratio ν and total mass M . The choice ν, M is indeed preferred as it is invariant under exchange of the two bodies 1 ↔ 2 [19].
Discussion.-In this Letter we performed the first comparison between numerical relativity (NR) and small mass-ratio (SMR) expanded results for a gauge invariant quantity that includes both dissipative and conservative effects, namely the accumulated orbital phase as a function of orbital frequency φ(M Ω). We have successfully extracted the post-adiabatic (PA) expansion of this quantity as a power series in the mass-ratio from non-spinning quasi-circular NR simulations.
The leading "adiabatic" (0PA) term agrees with the result from SMR calculations. In addition we obtain a robust determination of the 1PA term, serving as a concrete prediction for the ongoing SMR calculation of this term, which requires the dissipative part of the second order gravitational self-force. We also estimate the 2PA term φ 2PA from the NR data. Its amplitude is comparable to φ 0PA and φ 1PA for the frequency-range considered here, indicating that the PA-expansion remains wellbehaved. In particular, when the PA series is expanded in powers of the symmetric mass-ratio while keeping the total mass fixed, the 2PA and higher order terms are consistent with zero within the numerical accuracy for 0.015 M Ω 0.05. For higher frequencies (approaching the last stable circular orbit), we find indications of a transition-regime to plunge where the series in integer powers of ν is no longer applicable. Figure 4: Region of applicability of different approximation techniques for non-spinning quasi-circular binary black hole inspiral. The shaded regions indicate ranges within which the cumulative orbital phase-error is less than π/4, π/8 and π/16 radians, respectively.
Our analysis allows us to delineate the regions of applicability of SMR, NR and PN in a quantitative way, as shown in Fig. 4: Assuming φ 1PA will become available through GSF calculations, the envelope to the R 2+ in Fig. 2 gives a bound on the secular contributions of higher PA terms. The red shaded areas in Fig. 4 show the largest M Ω interval that can be covered such that the total accumulated phase-error due to ≥ 2PA terms is below a certain value. The region of applicability of SMR increases toward smaller mass-ratios, but is still non-negligible even at comparable masses. The post-Newtonian errors are estimated by fits against φ NR (M Ω), cf. top panel of Fig. 1. The green shaded areas indicate regions where the cumulative 3.5-PN phase-error for the entire inspiral up to the given frequency is below a certain value. Finally, the blue shaded area indicates the region covered by the NR simulations used here. These simulations have phase-accuracy better than the π/16 contour line, indicating that the usability of NR is not limited by accuracy but rather by the length of the simulations. The three modelling approaches deliver complementary information, covering different regions of the parameter space. The region of validity of each method depends on the desired accuracy, and it also depends on the use of the waveforms: For GW astronomy, only the accuracy within the portion of the inspiral within the frequency band of the relevant GW detectors is be important, and this will depend on the total mass of the binary. Moreover, the needed accuracy will depend on the signal-to-noise ratio at which it is observed.
We note that the adiabatic φ 0PA term is never accurate enough in the metric of Fig. 4, because φ 1PA contributes tens of radians in the frequency range considered, independent of the mass-ratio. This underlines the importance of calculating the 1PA term (and therefore the second order gravitational self-force) for modelling binaries of any mass-ratio. Furthermore, the application of the 1PA approximation for low frequencies is limited by a (M Ω) −1/3 divergence of the 2PA term. This motivates the development of models that incorporate both SMR and PN results, e.g. using effective-one-body theory [67].
The results in this paper come with two important caveats. First, our results only apply to non-spinning quasi-circular black hole binaries. Adding spin or eccentricity could make the convergence of the PA series significantly worse, and future studies are needed to explore the full parameter space. Even for non-spinning quasi-circular case, NR simulations at smaller mass-ratio are needed to investigate the transition to plunge, as well as longer simulations, to extend our analysis to smaller frequencies.
Second, the current analysis applies only to the inspiral, since the PA expansion is known to breakdown at the last stable orbit. Our results motivate the development of 1PA accurate models that also include plunge, merger, and ringdown, as has previously been done at 0PA order [68].