THE FINAL SPIN FROM BINARY BLACK HOLES IN QUASI-CIRCULAR ORBITS

According to the predictions of general relativity, binary systems of compact objects are the most efficient emitters of gravitational waves (GWs). Indeed, Advanced LIGO has recently detected the GW signal from a black hole (BH) binary with masses M₁ ≈ 36 M_⊙ and M₂ ≈ 29 M_⊙ (Abbott et al. 2016), at a (luminosity) distance of ∼410 Mpc. In general, Advanced LIGO and other terrestrial interferometers, such as Advanced Virgo and KAGRA, target BH binaries with a variety of masses (up to a few hundred M_⊙, if they exist; Belczynski et al. 2014, 2016). More massive BH binaries are targeted by existing pulsar-timing arrays (in the mass range 10⁸–10¹⁰ M_⊙; Manchester & IPTA 2013) and by future spaceborne interferometers such as eLISA (in the mass range 10⁴–10⁷ M_⊙; Klein et al. 2016).

One of the obvious difficulties of observing BH binaries with terrestrial interferometers is that only the final part of the inspiral and the merger/ringdown are in band. This is where the perturbative post-Newtonian (PN) techniques valid earlier in the inspiral become inaccurate, preventing the extraction of the source's physical parameters. Hence, to obtain the full gravitational waveforms, it is necessary to resort to numerical-relativity (NR) simulations. In practice, even under the reasonable assumption that BH binaries near the merger have been circularized by earlier GW emission, the space of parameters to be probed (the mass ratio q and the spin vectors ${{\boldsymbol{S}}}_{1},\;{{\boldsymbol{S}}}_{2}$, i.e., seven parameters) is too large to be handled by NR simulations alone.

To ensure a sufficient coverage of the parameter space, semianalytical techniques allowing faster waveform production are employed, e.g., the spin-effective-one-body (sEOB) model (Buonanno & Damour 1999; Damour 2001; Barausse & Buonanno 2010) or "hybrid" waveforms (Ajith et al. 2008; Khan et al. 2015), which combine results from NR simulations with PN and quasi-normal-mode calculations. These techniques are faster, but require great care when modeling the merger and the transition to the ringdown. Indeed, although the ringdown can be modeled via a linear superposition of quasi-normal modes, their frequencies depend on the remnant BH's mass and spin, which, in turn, depend on the initial binary parameters.

This relation between the binary's initial and final states is highly non-trivial because it encodes the details of the strong-field, highly relativistic merger, which is only accessible via NR calculations. Yet, a number of approaches to predict analytically or semianalytically the remnant's final-spin magnitude and direction have been proposed. These range from modeling the GW fluxes throughout the binary's evolution within the EOB model (e.g., by Damour & Nagar 2007 for nonspinning BHs) to approaches that combine information from PN theory, the extreme mass-ratio limit (EMRL), symmetry arguments, and fits to NR data, to provide "formulae" for the final spin (Buonanno et al. 2008; Kesden 2008; Rezzolla et al. 2008a, 2008b, 2008c; Tichy & Marronetti 2008; Barausse & Rezzolla 2009; Healy et al. 2014). Similar formulae have also been derived for the remnant's final mass (Kesden 2008; Tichy & Marronetti 2008; Barausse et al. 2012; Healy et al. 2014), which differs from the binary's total mass by the energy emitted in GWs. Again, a common problem in these attempts is the difficulty to cover with sufficient accuracy the seven-dimensional parameter space of quasi-circular BH binaries. Indeed, while most of these formulae formally cover the whole parameter space, they can be rather inaccurate, especially for BHs with almost extremal spins.

By combining results from NR and information from the EMRL and PN theory, we here derive a new formula for the spin magnitude and direction for the merger remnant from quasi-circular BH binaries with arbitrary masses and spins. We calibrate our formula against a catalog of 619 recently published NR simulations (Chu et al. 2009; Hannam et al. 2010; Kelly et al. 2011; Nakano et al. 2011; Pollney & Reisswig 2011; Sperhake et al. 2011; Buchman et al. 2012; Lovelace et al. 2012, 2015; Hemberger et al. 2013; Hinder et al. 2013; Kelly & Baker 2013; Pekowsky et al. 2013; Healy et al. 2014; Lousto & Zlochower 2014b; Scheel et al. 2015; Szilágyi et al. 2015; Zlochower & Lousto 2015; Husa et al. 2016; SXS collaboration 2016)⁵ and validate it by comparing its results to self-force calculations and plunge-merger-ringdown fluxes for nonspinning binaries with small mass ratios, as well as to a set of 248 NR simulations not included in the calibration data set (Karan et al. 2016).

Our new formula builds upon Barausse & Rezzolla (2009), who introduced a final-spin formula that is widely used both in the production of semianalytical waveforms (e.g., in sEOB and phenomenological waveforms) and in cosmological studies of massive BH evolution (see, e.g., Berti & Volonteri 2008; Fanidakis et al. 2011; Barausse 2012; Volonteri et al. 2013; Dubois et al. 2014; Sesana et al. 2014). Our novel prescription especially improves the accuracy of the formula by Barausse & Rezzolla (2009) for extreme mass ratios and for near-extremal spins. This is important since near-extremal spins are expected, at least in some cases, for supermassive BHs (Berti & Volonteri 2008; Fanidakis et al. 2011; Barausse 2012; Volonteri et al. 2013; Dubois et al. 2014; Sesana et al. 2014) and possibly also for stellar-mass BHs (McClintock et al. 2011). We assume G = 1 = c hereafter.

Let us first consider a BH binary with spins parallel (i.e., aligned or anti-aligned) to the orbital angular momentum ${\boldsymbol{L}}$ and denote the masses by ${M}_{\mathrm{1,2}}$ (with $q\equiv {M}_{2}/{M}_{1}\leqslant 1$) and the spin projections on the angular momentum direction by ${S}_{\mathrm{1,2}}\equiv {a}_{\mathrm{1,2}}\;{M}_{1,2}^{2}$ (${a}_{\mathrm{1,2}}$ being the dimensionless spin-parameter projections). In the EMRL q ≪ 1, the final-spin projection on the angular momentum direction must be

$$\begin{eqnarray}&&{a}_{{\rm{fin}}}={a}_{1}+\nu ({L}_{{\rm{ISCO}}}({a}_{1})-2{a}_{1}{E}_{{\rm{ISCO}}}({a}_{1}))+{ \mathcal O }({\nu }^{2}),\end{eqnarray} \tag{ 1 }$$

with $\nu \equiv q/{(1+q)}^{2}$ the symmetric mass ratio, and L_ISCO(a), E_ISCO(a), respectively, the specific (dimensionless) angular momentum and energy for a test particle at the innermost stable circular orbit (ISCO) of a Kerr BH with spin parameter a (Bardeen et al. 1972):

$$\begin{eqnarray}&&{E}_{{\rm{ISCO}}}(a)=\sqrt{1-\displaystyle \frac{2}{3{r}_{{\rm{ISCO}}}(a)}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{L}_{{\rm{ISCO}}}(a)=\displaystyle \frac{2}{3\sqrt{3}}[1+2\sqrt{3{r}_{{\rm{ISCO}}}(a)-2}],\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{r}_{{\rm{ISCO}}}(a)=3+{Z}_{2}-\displaystyle \frac{a}{| a| }\sqrt{(3-{Z}_{1})(3+{Z}_{1}+2{Z}_{2})},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{Z}_{1}=1+{(1-{a}^{2})}^{1/3}[{(1+a)}^{1/3}+{(1-a)}^{1/3}],\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{Z}_{2}=\sqrt{3{a}^{2}+{Z}_{1}^{2}}.\end{eqnarray} \tag{ 6 }$$

The final-spin expression of Rezzolla et al. (2008a) and Barausse & Rezzolla (2009) reproduces Equation (1) only in the special case a₁ = 0, when ${a}_{{\rm{fin}}}=2\sqrt{3}\nu +{ \mathcal O }({\nu }^{2})$. Indeed, one of the drawbacks of those early expressions is that they may yield spins a_fin > 1 for small mass ratios ν ≪ 1, in clear disagreement with Equation (1), which predicts a_fin ≤ 1, the equality holding for a₁ = 1.

To enforce the EMRL exactly, we consider the following ansatz for the final-spin projection:

$$\begin{eqnarray}{a}_{{\rm{fin}}} & = & {a}_{{\rm{tot}}}+\nu [{L}_{{\rm{ISCO}}}({a}_{{\rm{eff}}})-2{a}_{{\rm{tot}}}({E}_{{\rm{ISCO}}}({a}_{{\rm{eff}}})-1)]\\ & & +\displaystyle \sum _{i=0}^{{n}_{M}}\;\displaystyle \sum _{j=0}^{{n}_{J}}\;{k}_{{ij}}{\nu }^{2+i}\;{a}_{{\rm{eff}}}^{j},\end{eqnarray} \tag{ 7 }$$

where k_ij are free coefficients to be determined from the NR data, ${a}_{{\rm{tot}}}\equiv {({S}_{1}+{S}_{2})/({M}_{1}+{M}_{2})}^{2}=({a}_{1}+{a}_{2}{q}^{2})/{(1+q)}^{2}$ is the "total" spin parameter used in Barausse et al. (2012), while ${a}_{{\rm{eff}}}\equiv {S}_{{\rm{eff}}}/{({M}_{1}+{M}_{2})}^{2}$ is an "effective" spin parameter. In more detail, we assume ${S}_{{\rm{eff}}}=(1+\xi {M}_{2}/{M}_{1}){S}_{1}+(1+\xi {M}_{1}/{M}_{2}){S}_{2}$, which yields ${a}_{{\rm{eff}}}={a}_{{\rm{tot}}}+\xi \nu ({a}_{1}+{a}_{2})$. This choice is inspired by Damour (2001), who finds that the leading-order conservative spin–orbit dynamics depends on the spin only through S_eff with ξ = 3/4, while the leading-order conservative spin–spin dynamics depends on S_eff with ξ = 1 (see also Racine 2008; Gerosa et al. 2015b; Kesden et al. 2015). In the following, we will keep ξ as a free parameter and determine it from the NR simulations.⁶

Note that Equation (7) matches Equation (1) for ν ≪ 1 since ${a}_{{\rm{tot}}}={a}_{1}(1-2\nu )+{ \mathcal O }({\nu }^{2})$. Moreover, by singling out a_tot as the first term in Equation (7), we have isolated the "direct" contribution of the progenitor spins to the remnant's spin. However, this does not mean that all leading-order effects of the smaller BH's spin a₂ are already included. For instance, the specific energy and angular momentum at the ISCO receive corrections of ${ \mathcal O }({a}_{2}\;\nu )$ (see, e.g., Barausse & Buonanno 2010), which propagate into a term of ${ \mathcal O }({a}_{2}\;{\nu }^{2})$ in the final spin; cf. Equation (1). This effect, together with other ones, is captured by the coefficient k₀₁.

The coefficients k_0j of the ν² terms in Equation (7) also encode the information about the self-force dynamics (both dissipative and conservative) and the leading-order (in mass ratio) plunge-merger-ringdown emission. More specifically, the conservative self-force produces shifts $\nu \;{\rm{\Delta }}{E}_{{\rm{ISCO}}}$ and $\nu \;{\rm{\Delta }}{L}_{{\rm{ISCO}}}$ in the ISCO specific energy and angular momentum away from the geodesic values of Equations (2) and (3). For a nonspinning binary (a₁ = a₂ = 0) with ν ≪ 1,

$$\begin{eqnarray}&&{\rm{\Delta }}{L}_{{\rm{ISCO}}}\approx -0.802.\end{eqnarray} \tag{ 8 }$$

This follows from evaluating Equation (3c) of Le Tiec et al. (2012) at the ISCO frequency, which should include conservative self-force effects as in Equation (5) of the same reference. The plunge-merger-ringdown angular momentum flux is instead given by (Bernuzzi & Nagar 2010)

$$\begin{eqnarray}&&{\rm{\Delta }}{J}_{{\rm{MR}}}\approx 3.46\;{\nu }^{2}.\end{eqnarray} \tag{ 9 }$$

Therefore, for a nonspinning binary one expects

$$\begin{eqnarray}{a}_{{\rm{fin}}} & \approx & \;\nu {L}_{{\rm{ISCO}}}(0)-2{\nu }^{2}[{E}_{{\rm{ISCO}}}(0)-1]{L}_{{\rm{ISCO}}}(0)\\ & & +{\rm{\Delta }}{L}_{{\rm{ISCO}}}{\nu }^{2}-{\rm{\Delta }}{J}_{{\rm{MR}}}+{ \mathcal O }({\nu }^{3})\\ & \approx & \;2\sqrt{3}\nu -3.87\;{\nu }^{2}+{ \mathcal O }({\nu }^{3}),\end{eqnarray} \tag{ 10 }$$

hence k₀₀ ≈ −3.87. (Note that at ${ \mathcal O }(\nu )$ and after setting a₁ = a₂ = 0, this equation reduces to Equation (1).) However, since the transition from inspiral to plunge does not happen exactly at the ISCO when accounting for deviations from adiabaticity, but takes place smoothly around the ISCO (Ori & Thorne 2000), and since the the plunge-merger-ringdown fluxes are intrinsically approximate (as it is difficult to define unambiguously the plunge-merger-ringdown as separate from the late inspiral), we keep k₀₀ as a free parameter. As it happens, at least for n_M = 1, n_J = 2, the fitted value is k₀₀ ≈ −3.82, which is reasonably close to the one predicted by the considerations above.⁷

In principle, we could fit all the coefficients k_ij (as well as ξ) to the NR results. However, since simulations for equal-mass nonspinning BH binaries have determined the final remnant's spin with accuracy far better than for other configurations, we impose that Equation (7) with q = 1 and a₁ = a₂ = 0 yields exactly the final spin ${a}_{{\rm{fin}}}=0.68646\pm 0.00004$ measured by the NR simulations of Scheel et al. (2009). This gives the relation

$$\begin{eqnarray}&&\displaystyle \frac{\sqrt{3}}{2}+\displaystyle \sum _{i=0}^{{n}_{M}}\;\displaystyle \frac{{k}_{i0}}{{4}^{2+i}}=0.68646\pm 0.00004.\end{eqnarray} \tag{ 11 }$$

With this constraint, we fit Equation (7) to the 246 simulations for parallel-spin binaries in our calibration data set.

However, before performing the fit, it is useful to quantify the average error of the final spins calculated from NR simulations. This is possible because our calibration data set contains simulations by different groups with the same initial data. More precisely, 71 parallel-spin simulations have one or more "twins," i.e., binaries with exactly the same initial properties, so that the mean of the absolute differences between twin NR simulations can be measured to be $\delta {a}_{{\rm{fin}}}\approx 0.002$. This estimate allows not only performing a fit, but also computing its reduced chi-squared ${\chi }_{{\rm{red}}}^{2}$, thus gauging whether we are overfitting the data, which would correspond to ${\chi }_{{\rm{red}}}^{2}\lt 1$.

Since Equation (7) can be expanded to an arbitrary order via its last term, we have performed fits of the parallel-spin calibration data set with n_M = 1, n_J = 2 (6 coefficients), ${n}_{M}=3,{n}_{J}=3$ (16 coefficients), and n_M = 3, n_J = 4 (20 coefficients). The fitted coefficients are given in Table 1. Table 2 reports the mean (μ) and root-mean-square (rms, σ) of the residuals from the NR data, as well as ${\chi }_{{\rm{red}}}^{2}$, for the three aforementioned sets of coefficients, and for the formulae of Barausse & Rezzolla (2009), Husa et al. (2016), and Healy et al. (2014) (which use 4, 11, and 19 coefficients, respectively). Table 2 shows that our new formula converges when increasing the number of coefficients, although the optimal choice to avoid overfitting appears to be n_M = 1, n_J = 2. The convergence of our formula is also displayed in the left panel of Figure 1, which shows the probability distribution functions (PDFs) of the residuals, obtained as Gaussian fits. The right panel shows instead the (fitted) PDFs for our formula (for ${n}_{M}=3,{n}_{J}=4$) and for the formulae of Barausse & Rezzolla (2009), Husa et al. (2016), and Healy et al. (2014); the inset shows the actual residual distribution for our formula. Note that already with 16 coefficients our new formula has a slightly smaller rms than Healy et al. (2014), with the important advantage that it can also be used for generic binaries (see below), unlike the formulae of Husa et al. (2016) and Healy et al. (2014).

Zoom In Zoom Out Reset image size

Table 2.  The Mean and rms (μ and σ) of the Residuals ${a}_{\mathrm{fin}}^{\mathrm{num}}\mbox{--}{a}_{\mathrm{fin}}^{\mathrm{fit}}$ from the Numerical Data, as Well as ${\chi }_{{\rm{red}}}^{2}$, for Our Formula and those of Barausse & Rezzolla (2009), Husa et al. (2016), and Healy et al. (2014); also Displayed Is the Number of Coefficients in the Various Cases

Model	Coeffs.	μ	σ	${\chi }_{\mathrm{red}}^{2}$
nM = 1, nJ = 2	6	−0.000215	0.00198	0.985
nM = 3, nJ = 3	16	−0.000066	0.00168	0.712
${n}_{M}=3,{n}_{J}=4$	20	−0.000029	0.00166	0.694
Barausse & Rezzolla (2009)	4	−0.002310	0.00564	9.313
Husa et al. (2016)	11	−0.000240	0.00453	5.150
Healy et al. (2014)	19	0.000014	0.00170	0.718

Download table as:  ASCIITypeset image

To generalize Equation (7) to generic spins, we write the remnant's spin as the total spin ${\boldsymbol{S}}={{\boldsymbol{S}}}_{1}+{{\boldsymbol{S}}}_{2}$ plus an angular momentum contribution (i.e., the angular momentum at the binary's "effective" ISCO), i.e., ${{\boldsymbol{S}}}_{{\rm{fin}}}={\boldsymbol{S}}+{\rm{\Delta }}{\boldsymbol{L}}$. Since the final mass is ${M}_{{\rm{fin}}}=({M}_{1}+{M}_{2})(1-{E}_{{\rm{rad}}})$ (with ${E}_{{\rm{rad}}}\lesssim 0.1$ the mass radiated in GWs; Barausse et al. 2012), the final-spin parameter is

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{{\rm{fin}}}={{\boldsymbol{a}}}_{{\rm{tot}}}+{\boldsymbol{\ell }}\nu ,\quad {{\boldsymbol{a}}}_{{\rm{tot}}}=\displaystyle \frac{1}{{(1+q)}^{2}}({{\boldsymbol{a}}}_{1}+{{\boldsymbol{a}}}_{2}{q}^{2}),\end{eqnarray} \tag{ 12 }$$

where we have reabsorbed the radiated energy E_rad in ${\boldsymbol{\ell }}\equiv {\rm{\Delta }}{\boldsymbol{L}}/[{M}_{1}{M}_{2}{(1-{E}_{{\rm{rad}}})}^{2}]$ $+{\boldsymbol{S}}[2{E}_{{\rm{rad}}}+3{E}_{{\rm{rad}}}^{2}+{ \mathcal O }{({E}_{{\rm{rad}}})}^{3}]$ $/({M}_{1}{M}_{2})$ (note that ${\boldsymbol{\ell }}$ remains finite in the test-particle limit because $| {\rm{\Delta }}{\boldsymbol{L}}| ={ \mathcal O }(\nu )={E}_{{\rm{rad}}}$ as $\nu \to 0$). By evaluating Equation (12) for parallel spins and comparing it to (7), we obtain

$$\begin{eqnarray}| {\boldsymbol{\ell }}| & = & | {L}_{{\rm{ISCO}}}({a}_{{\rm{eff}}})-2{a}_{{\rm{tot}}}({E}_{{\rm{ISCO}}}({a}_{{\rm{eff}}})-1)\\ & & +\left.\displaystyle \sum _{i=0}^{{n}_{M}}\;\displaystyle \sum _{j=0}^{{n}_{S}}\;{k}_{{ij}}{\nu }^{1+i}\;{a}_{{\rm{eff}}}^{j}\right|,\end{eqnarray} \tag{ 13 }$$

which can be generalized to precessing spins by following Barausse et al. (2012; see also Rezzolla et al. 2008a; Barausse & Rezzolla 2009) and replacing

$$\begin{eqnarray}&&{a}_{{\rm{tot}}}\to {a}_{{\rm{tot}}}(\beta ,\gamma ,q)\equiv \displaystyle \frac{| {{\boldsymbol{a}}}_{1}| \mathrm{cos}\beta +| {{\boldsymbol{a}}}_{2}| \mathrm{cos}\gamma \;{q}^{2}}{{(1+q)}^{2}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{a}_{{\rm{eff}}}\to {a}_{{\rm{eff}}}(\beta ,\gamma ,q)\equiv {a}_{{\rm{tot}}}(\beta ,\gamma )+\xi \nu ({a}_{1}\mathrm{cos}\beta +{a}_{2}\mathrm{cos}\gamma ),\end{eqnarray} \tag{ 15 }$$

with β (γ) being the angle between ${{\boldsymbol{a}}}_{1}$ (${{\boldsymbol{a}}}_{2}$) and the orbital angular momentum. Clearly, with this choice, Equation (12) matches Equation (7) for parallel spins (i.e., for β = 0, π and γ = 0, π).

Moreover, for equal masses (q = 1), the leading-order PN spin effects in the conservative sector (i.e., the leading-order spin–orbit coupling) enter the dynamics only through the combination $\hat{{\boldsymbol{L}}}\cdot {\boldsymbol{S}}/{M}^{2}={a}_{{\rm{tot}}}(\beta ,\gamma ,1)\propto {a}_{{\rm{eff}}}(\beta ,\gamma ,1)$ (see, e.g., Damour 2001; Barausse & Buonanno 2010), where a "hat" denotes a unit-norm vector. Therefore, at this approximation order, the binding energy and angular momentum at the effective ISCO depend on the spins only through a_tot(β, γ, 1) (or equivalently a_eff(β, γ, 1)), as reflected in Equations (13)–(15). Similarly, in the EMRL, the leading contributions to $| {\boldsymbol{\ell }}|$ come from the ISCO energy and angular momentum of a test particle in Kerr. By construction, $| {\boldsymbol{\ell }}|$ has the correct EMRL for parallel spins, but the EMRL is also recovered approximately for generic-spin configurations, at least at leading order in the primary-BH spin. Indeed, this happens because the ISCO angular momentum and energy for a test particle in a non-equatorial orbit in a Kerr spacetime are L_ISCO(a_tot(β, γ, 0)) and E_ISCO(a_tot(β, γ, 0)), at leading order in the spin (see discussion in Barausse et al. 2012).

Putting things together, the final-spin magnitude reads

$$\begin{eqnarray}| {{\boldsymbol{a}}}_{{\rm{fin}}}| & = & \tfrac{1}{{(1+q)}^{2}}[| {{\boldsymbol{a}}}_{1}{| }^{2}+| {{\boldsymbol{a}}}_{2}{| }^{2}{q}^{4}+2| {{\boldsymbol{a}}}_{1}| | {{\boldsymbol{a}}}_{2}| {q}^{2}\mathrm{cos}\alpha \\ & & +{2(| {{\boldsymbol{a}}}_{1}| \mathrm{cos}\beta +| {{\boldsymbol{a}}}_{2}| {q}^{2}\mathrm{cos}\gamma )| {\boldsymbol{\ell }}| q+| {\boldsymbol{\ell }}{| }^{2}{q}^{2}]}^{1/2},\end{eqnarray} \tag{ 16 }$$

where α is the angle between the two spins. In principle, the angles α, β, and γ depend on the binary separation. However, β and γ enter in our formulae only through a_tot(β, γ, q) and a_eff(β, γ, q). These combinations remain constant during the adiabatic inspiral (Apostolatos et al. 1994), if only the leading PN order in the spins (i.e., the leading-order spin–orbit coupling) is included, and either (i) the masses are equal, or (ii) only one BH is spinning, or (iii) the mass ratio is extreme (i.e., ν ≈ 0). Similarly, under the same assumptions, we can safely assume that α remains constant during the adiabatic inspiral (Apostolatos et al. 1994), i.e., the angle between the two spins is preserved by the leading-order spin–orbit coupling for equal masses, while it does enter the final-spin prediction if only one BH is spinning, or when ν ≈ 0 (indeed, the effect of the smaller BH's spin vanishes at leading order in ν, because $| {{\boldsymbol{S}}}_{2}| ={ \mathcal O }{(\nu )}^{2}$). Outside these special cases, α, β, and γ are not exactly constant. For instance, in general, α oscillates and the oscillations may even become "flip-flop" unstable between separations ${r}_{\mathrm{ud}\pm }={(\sqrt{{a}_{1}}\pm \sqrt{q{a}_{2}})}^{4}{({M}_{1}+{M}_{2})/(1-q)}^{2}$ for certain unequal-mass configurations where the primary-BH spin is aligned with the orbital angular momentum and the spin of the secondary is anti-aligned with it (Gerosa et al. 2015a; Lousto & Healy 2016). These configurations, however, are unlikely if the spins are isotropically distributed or if the spins are almost aligned with the angular momentum of a circumbinary disk due to the Bardeen–Petterson effect (Bardeen & Petterson 1975).

Therefore, we follow Barausse et al. (2012), Barausse & Rezzolla (2009), and Rezzolla et al. (2008a) and define α, β, and γ at the initial binary separation r_in

$$\begin{eqnarray}\mathrm{cos}\alpha & \equiv & {\hat{{\boldsymbol{a}}}}_{1}\cdot {\hat{{\boldsymbol{a}}}}_{2}{| }_{{r}_{{\rm{in}}}},\\ \mathrm{cos}\beta & \equiv & \hat{{\boldsymbol{L}}}\cdot {\hat{{\boldsymbol{a}}}}_{1}{| }_{{r}_{{\rm{in}}}},\\ \mathrm{cos}\gamma & \equiv & \hat{{\boldsymbol{L}}}\cdot {\hat{{\boldsymbol{a}}}}_{2}{| }_{{r}_{{\rm{in}}}}.\end{eqnarray} \tag{ 17 }$$

Indeed, Barausse & Rezzolla (2009) and Kesden et al. (2010) have verified that the final-spin predictions are robust against the initial separation r_in, i.e., in most cases the definitions (17) are justified.

A comparison between our new formula with n_J = 2, n_M = 1 and the generic-spin simulations in our calibration data set yields the residuals displayed in Figure 2. Also shown is the corresponding PDF with mean μ ≈ −0.005 and rms σ ≈ 0.007. Note that in this case we cannot reliably estimate ${\chi }_{{\rm{red}}}^{2}$, as none of the generic-spin configurations have "twins" in our calibration data set, and the NR error is expected to be larger than in parallel-spin binaries because of precession. Also shown by Figure 2 is an unattractive feature of our formula, namely, that the distribution of residuals is biased toward negative values (i.e., our formula systematically overpredicts the final spin for generic binaries). Although this bias is small, and because it follows from assuming that α, β, and γ are constant, we can amend it by replacing the angles α, β, and γ by "effective" angles α*, β*, and γ* defined as

$$\begin{eqnarray}&&{{\rm{\Theta }}}^{* }\equiv 2\mathrm{arctan}\left[(1+{\epsilon }_{{\rm{\Theta }}})\mathrm{tan}\displaystyle \frac{{\rm{\Theta }}}{2}\right]\approx {\rm{\Theta }}+{\epsilon }_{{\rm{\Theta }}}\mathrm{sin}{\rm{\Theta }},\end{eqnarray} \tag{ 18 }$$

where Θ = α, β, γ, _Θ are free coefficients to be fixed by the data, and we impose _β = _γ to make our formula symmetric under exchange of the two BHs. Clearly, for parallel spins α* = α, β* = β, and γ* = γ. A comparison with the NR data gives _α ≈ 0 and _β = _γ ≈ 0.024, where we have used the second equality of Equation (18) (the first equality gives similar results). The corresponding residual distribution has a smaller bias and is shown in the inset of Figure 2, together with a PDF with μ ≈ −0.001, σ ≈ 0.006.

Zoom In Zoom Out Reset image size

As a further "blind" test of our formula, we consider data from the recently published catalog of Karan et al. (2016) that is not already included in our calibration data set (i.e., 83 parallel-spin and 165 precessing-spin simulations). Already when using only n_M = 1 and n_J = 2, the comparison yields mean and rms residuals $\mu \approx -5\times {10}^{-5}$ and σ ≈ 1.4 × 10⁻⁴ for parallel spins, and $\mu \approx -0.004$ (μ ≈ −0.0005) and σ ≈ 3.3 × 10⁻⁴ (σ ≈ 3.5 × 10⁻⁴) for precessing spins with unadjusted (adjusted) angles β and γ.

Finally, for the final-spin direction, we follow Barausse & Rezzolla (2009) and Apostolatos et al. (1994) and note that at leading PN order in the spins (i.e., including the leading-order spin–orbit coupling alone), the GW-driven evolution in the adiabatic inspiral approximately preserves the direction of the total angular momentum ${\boldsymbol{J}}\equiv {\boldsymbol{L}}+{\boldsymbol{S}}$. Barausse & Rezzolla (2009), and later Lousto & Zlochower (2014a), verified that $\hat{{\boldsymbol{J}}}$ is approximately preserved (to within a few degrees) also in the plunge, merger, and ringdown. The only exception to this "simple precession" picture are binaries with spins almost anti-aligned with the orbital angular momentum at large separations (Apostolatos et al. 1994; Kesden et al. 2010). Indeed, when the GW emission sheds enough angular momentum that ${\boldsymbol{L}}\approx -{\boldsymbol{S}}$, these binaries undergo "transitional precession" (Apostolatos et al. 1994), whereby the direction of $\hat{{\boldsymbol{J}}}$ changes significantly on short timescales. Note that among the configurations that give rise to "simple precession" are also the "flip-flop" binaries of Lousto & Healy (2016) and Gerosa et al. (2015a). Since transitional-precession configurations comprise a small portion of the parameter space (Kesden et al. 2010), we follow Barausse & Rezzolla (2009) and assume that the final-spin direction is simply given by $\hat{{\boldsymbol{J}}}({r}_{{\rm{in}}})$, i.e., the final-spin angle θ_fin relative to the initial angular momentum is simply

$$\begin{eqnarray}&&\mathrm{cos}{\theta }_{{\rm{fin}}}=\hat{{\boldsymbol{J}}}({r}_{{\rm{in}}})\cdot \hat{{\boldsymbol{L}}}({r}_{{\rm{in}}}).\end{eqnarray} \tag{ 19 }$$

Indeed, the 157 simulations (Lousto & Zlochower 2014b; Zlochower & Lousto 2015) in our data set that report the final-spin direction confirm that the final spin is almost aligned with $\hat{{\boldsymbol{J}}}({r}_{{\rm{in}}})$, to within ∼18° in the worst case and to within 4° (6°) in 64% (78%) of the cases. The distribution of the angle between the final spin and $\hat{{\boldsymbol{J}}}({r}_{{\rm{in}}})$ is shown in Figure 3; it is unclear whether the small counts for θ_fin ≳ 10° are due to imprecisions in the formula or in the numerical simulations.

Zoom In Zoom Out Reset image size

Finally, we note that unlike other formulae for the final-spin direction (Buonanno et al. 2008; Rezzolla et al. 2008a; Tichy & Marronetti 2008), Equation (19) is valid also when ${r}_{{\rm{in}}}\gg {M}_{1}+{M}_{2}$. (This is also the case for our formula for the final-spin magnitude). This is particularly important to predict the final spin in massive BH mergers. Indeed, cosmological simulations (both numerical and semianalytical ones) cannot follow the evolution of massive BH binaries below the separation r_GW at which the GW dynamics starts driving the orbital evolution. For a binary with ${M}_{1}+{M}_{2}\sim {10}^{8}\;{M}_{\odot }$ in a gas-poor environment, ${r}_{{\rm{GW}}}\sim {10}^{-2}\;{\rm{pc}}\sim 2\times {10}^{3}({M}_{1}+{M}_{2})$, a separation at which other prescriptions for the final-spin direction become significantly inaccurate (see discussion in Barausse & Rezzolla 2009; Barausse 2010).⁸

By combining information from the test-particle limit, perturbative/self-force calculations, the PN dynamics, and an extensive set of NR simulations collected from the literature, we have constructed a novel formula for the final spin from the merger of quasi-circular BH binaries with arbitrary mass ratios and spins. When applied to parallel-spin configurations, our novel formula performs better than other expressions in the literature, and we have also tested its validity for precessing-spin binaries, which other formulae are not able to model accurately. Also, unlike models such as that of Healy et al. (2014), our formula is purely algebraic. Finally, we have used our collected NR data set to confirm that the final-spin direction is almost parallel to the initial total angular momentum direction, as first suggested by Barausse & Rezzolla (2009).

We thank Nathan Johnson-McDaniel for useful comments and Davide Gerosa for clarifications on flip-flop binaries. We acknowledge support from the European Union's Seventh Framework Programme (FP7/PEOPLE-2011-CIG) through the Marie Curie Career Integration Grant GALFORMBHS PCIG11-GA-2012-321608, from the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904, and from the ERC Synergy Grant "BlackHoleCam—Imaging the Event Horizon of Black Holes" (grant 610058).