Kicking gravitational wave detectors with recoiling black holes

TL;DR

This paper presents new numerical simulations of black hole mergers with high recoil velocities, introduces a novel method to analyze recoil dependence on merger phase, and assesses the detectability of these events with gravitational wave detectors.

Contribution

It provides eight new simulations of maximally spinning black hole mergers, introduces an invariant waveform analysis method, and evaluates the observability of high recoil events.

Findings

Maximum recoil velocities up to ~4700 km/s were achieved.

The waveform peak amplitude serves as an effective phase reference.

Highly recoiling black holes can be distinguished by gravitational wave detectors.

Abstract

Binary black holes emit gravitational radiation with net linear momentum leading to a retreat of the final remnant black hole that can reach up to $\sim5,000$ km/s. Full numerical relativity simulations are the only tool to accurately compute these recoils since they are largely produced when the black hole horizons are about to merge and they are strongly dependent on their spin orientations at that moment. We present eight new numerical simulations of BBH in the hangup-kick configuration family, leading to the maximum recoil. Black holes are equal mass and near maximally spinning ($|\vec{S}_{1,2}|/m_{1,2}^2=0.97$). Depending on their phase at merger, this family leads to $\sim\pm4,700$ km/s and all intermediate values of the recoil along the orbital angular momentum of the binary system. We introduce a new invariant method to evaluate the recoil dependence on the merger phase via the waveform peak amplitude used as a reference phase angle and compare it with previous definitions. We also compute the mismatch between these hangup-kick waveforms to infer their observable differentiability by gravitational wave detectors, such as advanced LIGO, finding currently reachable signal-to-noise ratios, hence allowing for the identification of highly recoiling black holes having otherwise essentially the same binary parameters.

Tables

Table 1. Table 1: Initial simple proper distance and spins of the BHs. The initial coordinate separation in all cases is D / m = 9 𝐷 𝑚 9 D/m=9 and mass ratio q = 1 𝑞 1 q=1 . ADM masses are between 0.9880 and 0.9884. Each simulation in the series can be uniquely described by the azimuthal spin angle, φ 𝜑 \varphi . Both spins have magnitude 0.97 and polar angle θ 𝜃 \theta of 50.98 ∘ superscript 50.98 50.98^{\circ} , and the angle Δ ​ φ Δ 𝜑 \Delta\varphi between α → 1 subscript → 𝛼 1 \vec{\alpha}_{1} and α → 2 subscript → 𝛼 2 \vec{\alpha}_{2} is 180 ∘ superscript 180 180^{\circ} . In terms of the spin components α → 1 = ( α 1 ​ x , α 1 ​ y , α 1 ​ z ) = ( − α 2 ​ x , − α 2 ​ y , α 2 ​ z ) subscript → 𝛼 1 subscript 𝛼 1 𝑥 subscript 𝛼 1 𝑦 subscript 𝛼 1 𝑧 subscript 𝛼 2 𝑥 subscript 𝛼 2 𝑦 subscript 𝛼 2 𝑧 \vec{\alpha}_{1}=(\alpha_{1x},\alpha_{1y},\alpha_{1z})=(-\alpha_{2x},-\alpha_{2y},\alpha_{2z}) .

$\phi$	$m{\Omega }_{22}$	$d/m$	${\alpha }_{2x}$	${\alpha }_{2y}$	${\alpha }_{2z}$
0	0.01032	12.5183	0.7536	0.0000	0.6107
30	0.01044	12.4045	0.6527	0.3768	0.6107
60	0.01053	12.2011	0.3768	0.6527	0.6107
90	0.01055	12.2128	0.0000	0.7536	0.6107
120	0.01051	12.3744	-0.3768	0.6527	0.6107
150	0.01046	12.4913	-0.6527	0.3768	0.6107
203	0.01046	12.4455	-0.6953	-0.2908	0.6107
291	0.01029	12.3250	0.2663	-0.7050	0.6107

Table 2. Table 2: Final properties of the remnant BH. The final mass M f / m subscript 𝑀 𝑓 𝑚 M_{f}/m , final spin α f subscript 𝛼 𝑓 \alpha_{f} , recoil velocity in km/s, and peak luminosity in ergs/s are given. The number of orbits before merger and time of peak luminosity are also given. Δ ​ ϕ Δ italic-ϕ \Delta\phi representing the relative phases with respect to the φ = 0 𝜑 0 \varphi=0 case.

$\phi$	$\Delta {\varphi }_{peak}$	$\Delta {\varphi }_{traj}$	$2N_{orbits}$	$M_f/m$	${\alpha }_f$	$V_{recoil}$	peak Lum.	$t_{Hpeak}$
0	0	0	14.0095	0.9251	0.8525	-4014	5.603$\times {10}^{56}$	860.5
30	33.05	29.85	13.9915	0.9217	0.8461	-4622	6.076$\times {10}^{56}$	860.1
60	65.99	77.16	13.9859	0.9200	0.8446	-3882	6.228$\times {10}^{56}$	859.2
90	86.17	120.87	13.9689	0.9215	0.8496	-1846	5.811$\times {10}^{56}$	859.8
120	106.44	143.48	13.9685	0.9244	0.8550	531	5.390$\times {10}^{56}$	851.4
150	142.53	160.00	14.0011	0.9260	0.8565	2553	5.326$\times {10}^{56}$	857.3
203	203.51	201.85	13.9950	0.9225	0.8475	4579	5.965$\times {10}^{56}$	860.4
291	264.09	320.23	13.9673	0.9245	0.8536	186	5.487$\times {10}^{56}$	861.3

Table 3. Table 3: A fit A 1 ​ cos ⁡ ( Δ ​ ϕ − ϕ 1 ) + A 3 ​ cos ⁡ ( 3 ​ [ Δ ​ ϕ − ϕ 3 ] ) subscript 𝐴 1 Δ italic-ϕ subscript italic-ϕ 1 subscript 𝐴 3 3 delimited-[] Δ italic-ϕ subscript italic-ϕ 3 A_{1}\cos(\Delta\phi-\phi_{1})+A_{3}\cos(3[\Delta\phi-\phi_{3}]) to the V r ​ e ​ c ​ o ​ i ​ l subscript 𝑉 𝑟 𝑒 𝑐 𝑜 𝑖 𝑙 V_{recoil} with 4 degrees of freedom. For the three different evaluations of Δ ​ ϕ = φ = Δ italic-ϕ 𝜑 absent \Delta\phi=\varphi= initial angle (in red circles), Δ ​ ϕ t ​ r ​ a ​ j = Δ subscript italic-ϕ 𝑡 𝑟 𝑎 𝑗 absent \Delta\phi_{traj}= trajectory angle as defined in Lousto and Zlochower ( 2013 ) (in magenta triangles)), and Δ ​ ϕ p ​ e ​ a ​ k Δ subscript italic-ϕ 𝑝 𝑒 𝑎 𝑘 \Delta\phi_{peak} from the waveform phase at the peak amplitude (in blue squares), as defined in this paper.

Parameters	Initial angle	Standard	Error	Trajectory angle	Standard	Error	Waveform phase	Standard	Error
$A_1$	4569.47	$\pm 3.825$	(0.083%)	4678.96	$\pm 408.2$	(8.724%)	4678.88	$\pm 513.0$	(10.96%)
${\varphi }_1$	0.4353	$\pm 0.0008$	(0.074%)	0.7960	$\pm 0.0731$	(9.432%)	0.2447	$\pm 0.1094$	(8.253%)
$A_3$	152.22	$\pm 3.822$	(2.511%)	10.0268	$\pm 388.4$	(3873%)	9.96288	$\pm 551.1$	(5531%)
${\varphi }_3$	0.8814	$\pm 0.00840$	(0.147%)	0.0617	$\pm 12.1$	(741%)	0.7434	$\pm 16.55$	(883%)

Table 4. Table 4: Mode pair contributions to the recoil velocity in the z 𝑧 z -direction as in Eq. III.2 for the near maximal, near zero, and near minimal recoil configurations. Only pairs with contributions > 10 absent 10 >10 km/s are included here.

${\ell }_1$	$m_1$	${\ell }_2$	$m_2$	$V\left(\phi ={30}^{\circ }\right)$	$V\left(\phi ={291}^{\circ }\right)$	$V\left(\phi ={203}^{\circ }\right)$
2	2	2	2	9122.37	6779.79	4818.65
2	-2	2	-2	-4893.59	-6865.65	-9019.23
2	-2	3	-2	-228.74	-435.80	-507.18
2	2	3	2	521.46	334.62	227.70
3	2	3	2	26.51	14.06	10.35
3	-2	3	-2	-10.09	-25.43	-25.50
4	4	4	4	85.99	47.51	20.80
4	-4	4	-4	-21.93	-32.93	-84.98

Table 5. Table 5: Convergence of key quantities for the φ = 291 ∘ 𝜑 superscript 291 \varphi=291^{\circ} system with three resolutions. Richardson extrapolation is used to determine the convergence order and infinitely extrapolated values. Recoil velocities are given in km/s and peak luminosities are ergs/s. The fifth row shows the difference between the extrapolated and N144 values, and the sixth row shows the percent difference between the two. There is an exception for the quantity in the last column, ϕ h 22 , peak subscript italic-ϕ subscript ℎ 22 peak \phi_{h_{22},\mathrm{peak}} . If we were to take the phase at a fixed time near peak for each resolution, we would observe and order of convergence between 5 and 6. However, since we take the phase at the peak for each resolution, and the time of peak is already convergent at an order of 5.5, we observe higher than normal convergence for the phase when measured this way.

	$V_{recoil}$	${\alpha }_f$	$M_f/m$	${10}^{-56}\cdot {ℒ}_{\mathrm{peak}}$	${\left|r{h}_{22}\right|}_{\mathrm{peak}}$	$t_{h_{22},\mathrm{peak}}$	${\varphi }_{h_{22},\mathrm{peak}}$
N100	227.41	0.853399	0.923310	5.4062	0.475254	962.804	89.793
N120	193.35	0.853569	0.923599	5.4578	0.476050	962.595	89.800
N144	186.03	0.853642	0.923705	5.4867	0.476328	962.519	89.801
Inf. Extrap.	184.03	0.853697	0.923766	5.5235	0.476476	962.476	*
Inf. - N144	-2.00	0.000055	0.000061	0.0368	0.000148	-0.043	*
% difference	-1.09	0.0065	0.0066	0.6673	0.0311	-0.005	*
Conv. Order	8.4	4.6	5.5	3.2	5.8	5.5	*

Equations

$$V_{\mathrm{rec}}=V_1\cos (\Delta \varphi +{\varphi }_1)+V_3\cos (3\Delta \varphi +3{\varphi }_3),$$

$$V_1=$$

$${\tilde{V}}_1=\left(\frac{1+E\alpha \cos \theta }{1+F\alpha \cos \theta }\right)D\alpha \sin \theta$$

$$\frac{d{P}_z}{dt}$$

$$c_{l,m}$$

$$d_{l,m}$$

$${\Psi }_{\infty }=\frac{{\Psi }_H{\Psi }_L-{{\Psi }_M}^2}{{\Psi }_H-2{\Psi }_M+{\Psi }_L}$$

$$A=\frac{{{\Psi }_L}^2-2{\Psi }_L{\Psi }_M+{{\Psi }_M}^2}{{\Psi }_H-2{\Psi }_M+{\Psi }_L},$$

$$n=-\frac{1}{\ln \left(f\right)}\ln \left|\frac{{\Psi }_H-{\Psi }_M}{{\Psi }_L-{\Psi }_M}\right|,$$

$$\mathcal{M}=⟨a|b{⟩}_k\equiv 2{\int }_{|f|>f_{\min }}df\frac{[\tilde{a}(f){]}^{\ast }\tilde{b}(f)}{S_{h,k}(|f|)}.$$

$${\mathcal{M}}_i(\xi ,\psi )=⟨{\phi }_{ref}[{\xi }_{ref},{\psi }_{ref}]|{\phi }_i[\xi ,\psi ]⟩,$$

Full text

I Introduction

Soon after numerical relativity simulations Campanelli et al. (2007a, b) neatly revealed that astrophysical binary black holes may impart speeds of thousands of kilometers per seconds after merger on the final black hole through gravitational recoil, a search for them intensified in the astronomical community. These searches ranged from dynamical effects of their host galaxies Volonteri (2007); Loeb (2007); Holley-Bockelmann et al. (2008); sesana (2007); Blecha and Loeb (2008) leading to displacements from galaxy cores, to specific objects displaying features that could be interpreted as differential velocities of thousand of kilometers per second between narrow and broad emission lines, like CID-42 Civano et al. (2010); Blecha et al. (2012); Lanzuisi et al. (2013), J0927+2943 Bogdanović et al. (2009); Komossa et al. (2008); Vivek et al. (2009); Shields et al. (2009a); Decarli et al. (2014), J1225+1415 Jonker et al. (2010), J1050+3456 Shields et al. (2009b), and NGC1277 Shields and Bonning (2013). A systematic search was carried out and described in Eracleous et al. (2012); Lena et al. (2014); Runnoe et al. (2015, 2017). A particularly promising study of 3C186 Chiaberge et al. (2017); Lousto et al. (2017); Chiaberge et al. (2018) is currently underway. Early reviews on this field are given in Refs. Komossa (2012); Blecha et al. (2016).

Systematic numerical relativity simulations provided a method to model the recoil of the final merged black hole as a function of the precursor binary Lousto et al. (2012); Zlochower and Lousto (2015), and to determine that the maximum recoil is about 5,000 km/s for maximally spinning, equal mass, holes in the hangup kick configuration Lousto and Zlochower (2011a). Aligned spins, on the other hand, can only reach a maximum of just above 500 km/s, in an antialigned configuration with mass ratio $q\sim 2/3$ Healy et al. (2014, 2018a). While nonspinning holes only contribute about one third of this value González et al. (2007); Healy et al. (2017a). See a review of the field in Sperhake (2015). Numerical studies can also include accreting mater to determine electromagnetic counterparts of the recoil Sijacki et al. (2011); Ponce et al. (2012); Meliani et al. (2017).

Interestingly, the observability of these recoils with gravitational wave detectors Gerosa and Moore (2016); Calderón Bustillo et al. (2018) has been explored recently. Here we test this question in the most favorable scenario, that of the hangup-kick recoil with explicit simulations of nearly maximal spins ($\alpha=|\vec{\alpha}_{1,2}|=|\vec{S}_{1,2}|/m^{2}_{1,2}=0.97$). We compare waveforms for configurations within the hangup-kick family (See Fig. 1) leading to nearly maximally but opposed recoils and passing through essentially vanishing recoil to see the required signal-to-noise ratio to distinguish between them with the analytic advanced LIGO sensitivity curve Adv .

This paper is organized as follows, in the next section II we describe the numerical relativity techniques that we will use in the evolutions of the binary black holes. In section III we describe the results of the simulations within the hangup family with equal mass black holes and spin magnitudes $\alpha=0.97$ for eight different spin orientations. This systematic study provides a method to fit a sinusoidal dependence of the recoil velocity of the final black hole as a function of the spin orientation. A new technique to identify this relative spin orientation at merger from the waveform phase is described in subsection III.1. We also analyze in subsection III.3 the finite difference errors of our simulations by studying one member of the family with three resolutions and assess the differences with respect to its extrapolated value. We end the paper with a discussion, in section IV, of the potentially observable recoil effects on these waveforms. We evaluate the matching of our simulations with each other, taking as a reference the one with the lowest recoil velocity, to see the signal-to-noise (SNR) requirements to distinguish one from the other by advanced LIGO. We also come back to the first gravitational wave event GW150914, that we recently reanalyzed in Ref. Healy et al. (2019), to evaluate the likelihood of recoils within a different simulation family, involving one single spinning black hole.

II Numerical Techniques

The late orbital dynamics of spinning binary black holes remain a fascinating area of research since the numerical breakthroughs Pretorius (2005); Campanelli et al. (2006a); Baker et al. (2006) solved the binary black hole problem via supercomputer simulations. Among the notable spin effects (without Newtonian analogs) observed in supercomputer simulations are the hangup effect Campanelli et al. (2006b), which prompts or delays the merger of binary black holes depending on the spin-orbit coupling, $\vec{S}\cdot\vec{L}$, being positive or negative (aligned spins or antialigned spins with the orbital angular momentum $\vec{L}$); the flip-flop of individual black hole spins passing from aligned to antialigned periods with respect to the orbital angular momentum Lousto et al. (2016) the alignment instability Kesden et al. (2015) as a case of imaginary flip-flop frequencies Lousto and Healy (2016); and the total flip of the orbital angular momentum Lousto and Healy (2019) under generic retrograde orbits for intermediate mass ratio binaries $(q<1/4)$.

Perhaps one of the most notable predictions of numerical relativity are the large recoils (thousands of km/s) of the final black hole remnant Campanelli et al. (2007a), and up to 5,000 km/s Lousto and Zlochower (2011a). Those results have been obtained from simulations with spinning black holes of $\alpha=S/m^{2}=0.7,0.8,0.9$ and extrapolated to maximally spinning holes. More recently, we introduced a version of highly-spinning initial data, based on the superposition of two Kerr black holes Ruchlin et al. (2017); Healy et al. (2016), in a puncture gauge that can easily be incorporated into moving-punctures codes. In Refs. Ruchlin et al. (2017); Zlochower et al. (2017), we were able to evolve an equal-mass binary with aligned spins, and spin magnitudes of $\alpha=0.95$ and $\alpha=0.99$ respectively, using this new data and compare with the SXS results of Ref. Mroue et al. (2013), finding excellent agreement.

In order to verify the extrapolation to near maximally spinning black holes and its evolution for a precessing system (in particular here the binary has a bobbing motion), we designed a set of 8 simulations in the hangup-kick configuration with spins $\alpha=0.97$. These simulations in turn will allow us to single out the effect of recoil as a function of its merger phase and their observability with gravitational wave detectors.

In table 1 we provide the 8 configurations spanning different initial orientations of the spin projection onto the orbital plane $S_{\perp}$, with respect to the line joining each hole as described by the angle $\varphi$, and are chosen to include near maximum recoil in both z-directions ($\vec{L}$) and near zero recoil. Here $\varphi=\phi(t=0)$ and at that initial time $S_{x}=S_{\perp}\cos\varphi$ and $S_{y}=S_{\perp}\sin\varphi$ for one black hole and reversed signs for the other. The polar angle $\theta$ of the spin with respect to the orbital angular momentum $\vec{L}$ has been chosen to maximize the recoil according to the predictions in Ref. Lousto et al. (2012), i.e. reproduced here in Eqs. (2), (3); and evaluated for $\alpha=0.97$ give the value $\theta_{max}=50.98$ degrees.

We have chosen the initial separations of the binaries to guarantee around 7 orbits before merger and in order to estimate the accuracy of the finite resolution used in those simulations we perform three simulations for a member of the family (that with $\varphi=291^{\circ}$), at increasing resolutions by a factor 1.2 in order to study the convergence of the relevant quantities for our studies. Those results are reported later in subsection III.3.

We evolve the binary black hole data sets using the LazEv Zlochower et al. (2005) implementation of the moving puncture approach Campanelli et al. (2006a) with the conformal function $W=\sqrt{\chi}=\exp(-2\phi)$ suggested by Ref. Marronetti et al. (2008). For the runs presented here, we use centered, eighth-order finite differencing in space Lousto and Zlochower (2008), a fourth-order Runge Kutta time integrator, and a 7th-order Kreiss-Oliger dissipation operator. We use a Courant factor of 0.25 in the CCZ4 formulation of the evolution equations Alic et al. (2012). Our code uses the EinsteinToolkit Löffler et al. (2012); ein / Cactus cac / Carpet Schnetter et al. (2004) infrastructure. The Carpet mesh refinement driver provides a “moving boxes” style of mesh refinement. In this approach, refined grids of fixed size are arranged about the coordinate centers of both holes. The evolution code then moves these fine grids about the computational domain by following the trajectories of the two black holes. We use AHFinderDirect Thornburg (2004) to locate apparent horizons. We measure at it the mass and the magnitude of the horizon spin using the isolated horizon (IH) algorithm detailed in Ref. Dreyer et al. (2003) and as implemented in Ref. Campanelli et al. (2007c). We measure radiated energy, linear momentum, and angular momentum, in terms of the radiative Weyl Scalar $\Psi_{4}$, using the formulas provided in Refs. Campanelli and Lousto (1999); Lousto and Zlochower (2007). We extract the radiated energy-momentum at finite radius and extrapolate to $r=\infty$ with the perturbative extrapolation described in Ref. Nakano et al. (2015). For the radiated quantities, we use all modes up to and including $\ell_{max}=6$. Quasicircular (low eccentricity) initial orbital parameters are computed using the 3rd. order post-Newtonian techniques described in Healy et al. (2017b).

III Results

We summarize the results of our BBH evolutions in Table 2 where the final black hole remnant properties and peak waveform luminosity values are reported. The modeling of remnant mass and spin for precessing binaries is given in Ref. Lousto and Zlochower (2014); Zlochower and Lousto (2015) with both quantities bearing a $\cos 2\phi$-dependence for the current family of simulations. Here, we will particularly focus on the analysis of the recoil velocities with regards to the predictions for those simulations with high spin $(\alpha=0.97)$ from the extrapolation of previous fitting formulae cfr. in equations (2) or (3).

In order to analyze the results of the present simulations, We can fit the recoil to the form

[TABLE]

where $V_{1}$, $V_{3}$, $\phi_{1}$, and $\phi_{3}$ are fitting parameters and $\Delta\phi$ is the initial phase of the spin with respect to a reference direction (in our case the y-axis).

Based on Lousto and Zlochower (2011a), we expected that the recoil would have the form

[TABLE]

where $V_{1}$ is the component of the recoil proportional to $\cos\phi$, $V_{1,1}$ arises from the “superkick” formula, and the remaining terms are proportional to linear, quadratic, and higher orders in $S_{z}/m^{2}=\alpha\cos\theta$ (the spin component in the direction of the orbital angular momentum).

A fit of the simulations reported in Lousto et al. (2012) to this ansatz (2) showed that the truncated series appears to converge very slowly with coefficients $V_{1,1}=(3677.76\pm 15.17)\ \rm km\ s^{-1}$, $V_{A}=(2481.21\pm 67.09)\ \rm km\ s^{-1}$, $V_{B}=(1792.45\pm 92.98)\ \rm km\ s^{-1}$, $V_{C}=(1506.52\pm 286.61)\ \rm km\ s^{-1}$ that have relatively large uncertainties. In what follows we will neglect the contribution of $V_{3}\sim 100$km/s; see Lousto et al. (2012) for more details.

In addition, we proposed the alternative modeling

[TABLE]

which can be thought of as a resummation of Eq. (2) with an additional term $E\alpha\cos\theta$, and fit to $D$, $E$, $F$ (where we used the prediction of Lousto and Zlochower (2011b) to model the $\tilde{V}_{1}$ for $\theta=90^{\circ}$) and found $D=(3684.73\pm 5.67)\ \rm km\ s^{-1}$, $E=0.0705\pm 0.0127$, and $F=-0.6238\pm 0.0098$. Note that $E$ is approximately $1/9$ of $F$, indicating that coefficients in this series get progressively smaller in a faster sequence than in Eq. (2).

We can fit to the recoil dependence on the initial angle $\phi$ between the spin and the y-axis. Alternatively, one can seek a reference frame, closer to merger, when most of the net recoil appears to be generated. In Refs. Lousto and Zlochower (2013); Zlochower and Lousto (2015) we have described in a totally coordinate based frame (punctures trajectories) the way to extract an instantaneous orbital plane and spin projections as displayed in Figure 3 of reference Lousto and Zlochower (2013) or Figure 1 of Zlochower and Lousto (2015). We implement here a new measure of this angle about merger with respect to the $\varphi=0$ case as a reference. We introduce the notion of using the peak amplitude phase of the waveform $\phi_{peak}$, as a measure for the reference phase of the recoil modeling and provide more detail in subsection III.1.

These fits are displayed in Fig. 2 giving rise to an estimate of the maximum recoil for these configurations in the form of the fitted amplitude of the sinusoidal dependence on $\Delta\phi$ as given by Eq. (1) to extract the leading $\cos\Delta\phi$-dependence and have a control of the nonleading $\cos 3\Delta\phi$ term. The three different evaluations of $\Delta\phi=\varphi=$ initial angle (in red circles), $\Delta\phi_{traj}=$ trajectory angle as defined in Lousto and Zlochower (2013) (in magenta triangles), and $\Delta\phi_{peak}$ from the waveform phase at the peak amplitude (in blue squares), as defined in subsection III.1 below.

Table 3 displays the measured fitting parameters and its statistical asymptotic standard errors with 4 degrees of freedom. Evaluating Eqs. (2) and (3) with the parameters for the series studied here ($\alpha=0.97$ and $\theta=50.98^{\circ}$), we find $V_{1}=4675.97\pm 64.71$ and $\tilde{V}_{1}=4678.90\pm 57.52$ km/s respectively. Comparing to the three fits given in table 3, we see excellent agreement when using $\Delta\phi_{traj}$ ($4678.96\pm 40.82$ km/s) and $\Delta\phi_{peak}$ ($4678.90\pm 57.52$ km/s).

III.1 Reference frame at peak waveform amplitude

The peak amplitude $h_{peak}^{22}$ and peak waveform frequency $\Omega_{peak}^{22}$ modeling in aligned binaries simulations was introduced in Ref. Healy and Lousto (2018). Here we use its definition to determine a reference time and hence phase of the waveform at which we can assign a recoil dependence of the form (1) and as represented in Fig. 3. We compare this gauge invariant method to determine the differential (near merger) phase dependence to the coordinate based method of Lousto and Zlochower (2013); Zlochower and Lousto (2015) that was used in the original hangup-kick work Lousto and Zlochower (2011a) and to determine the numerical coefficients in Eqs. (2) and (3). Note that the two methods defined using a (near merger) measure as reference lead to very similar results. The statistical errors of those methods appear much larger than those measured from the initial angle $\phi$ given the difficulties in measuring directions in the strong dynamical regime of the merger phase.

The notion that the phase of the waveform at peak luminosity as a reference in the strong field regime, near the merger of the two black holes, is a very interesting one, since it is amenable to be generalized in the fully precessing case. In that case one has to determine the orbital plane orientation from the direction of the maximum power of gravitational waves at the peak luminosity. Also measure the phase of the waveform along that privileged direction. Appropriate families of simulations to extract modeling parameters should then be designed. This will be the subject of a future research by the authors.

III.2 Recoil Generation

These systems provide an illustrative example of how the recoil is cumulated during late inspiral, merger, and ringdown. Due to the symmetry of these systems, the recoil of the remnant BH is solely in the $z$-direction, which is aligned with the gravitational wave extraction frame. The recoil can be calculated from individual modes of $\Psi_{4}=\sum_{l=2}^{\infty}\sum_{m=-l}^{l}A^{l,m}\left({}_{-2}Y^{l,m}(\theta,\phi)\right)$ by Eqs. (3.15), (3.18), and (3.19) in Ruiz et al. (2008):

[TABLE]

Table 4 shows the contributions to the recoil from the mode pairs of Eq. III.2 that contribute more than 10 km/s for the three simulations that appear in Fig. 4. These three simulations are the ones with the near maximal, near zero, and near minimal recoil velocities (top to bottom). To good approximation, when the amplitude of the (2,2) mode is larger than the amplitude of the (2,-2) mode, the recoil velocity will increase. This is easiest to see near merger, as in the top panel of Fig. 4, but is true throughout. In this panel, the red (2,2) dominates from late inspiral through ringdown, resulting in a near maximal recoil for these configurations. In the bottom panel, the opposite is true, the blue (2,-2) dominates over the same range, and the recoil is approximately the same, but in the opposite direction (note the $y$-axis on the right is reversed). The middle panel is interesting, in that it exhibits a late-time continuation of the orbital wobbling leading to an in-phase cancellation or anti-kick, where at first we obtain a large recoil (around 1,000 km/s) followed by another large recoil which cancels the original, resulting in a final recoil close to 0. This anti-kick can be explained again by which mode is dominating near merger. At first, the blue (2,-2) is dominating in the late inspiral, but as we approach the peak, red (2,2) dominates, producing the large positive recoil. However, during ringdown, blue (2,-2) dominates again producing the large negative recoil cancellation. Table 4 shows that the contributions of the (2,2) and (2,-2) mode with themselves produce the largest contributions to the recoil, but will always carry an opposite sign (because of the $c_{lm}$ coefficient.) For the near maximal and near minimal configurations, these two modes account for approximately $90\%$ of the kick, leaving the remaining approximately 400 km/s to the other mode pairs. Interestingly in the near zero configuration, the (2,2) and (2,-2) mode pairs only contribute 85 km/s after the cancellation, leaving the bulk of the recoil (an additional 100 km/s) to the higher mode pairs. If the same analysis were applied to an aligned system, where the spins are aligned with the orbital angular momentum, we would still obtain very large recoil contributions from the (2,2) and (2,-2) mode pairs. However, due to the symmetry, these would cancel completely (and all other mode pairs), to give a net-zero recoil in the $z$-direction.

III.3 Convergence of the numerical simulations

Numerous convergence studies of our past simulations have been performed. In Appendix A of Ref. Healy et al. (2014), in Appendix B of Ref. Healy and Lousto (2017), and for nonspinning binaries are reported in Ref. Healy et al. (2017a). For very highly spinning black holes ($s/m^{2}=0.99$) convergence of evolutions was studied in Ref. Zlochower et al. (2017) and for ($s/m^{2}=0.95$) in Ref. Healy et al. (2018a) for unequal mass binaries. For our current simulations we studied in detail one member of the hangup kick family, that with the lowest recoil, at an initial spin orientation angle $\phi=291^{\circ}$. With three resolutions, lowering our standard resolution for the whole family by factors of 1.2. Resolutions are denoted by “NXXX”, where XXX is a number related to the resolution in the wavezone. For example, “N144”, the standard resolution for these simulations, has a wavezone resolution of M/1.44, and “N100”, has a resolution of M/1.00. We then assume that a quantity $\Psi$ behaves with resolution $h$ in the range of low $h_{L}$ to high $h_{H}$ as $\Psi(h)=\Psi(\text{exact})+A\,h^{n}$, where we compute $\Psi(h)$ at the three resolutions $h_{L},h_{M},h_{H}$. We evaluate the extrapolation to infinite resolution $\Psi(\text{exact})\to\Psi_{\infty}$ as

[TABLE]

where we also determine the constant $A$ and the convergence rate $n$. We have also assumed that the low, medium, and high resolutions are related by a common factor $f$ as $h_{M}=h_{L}/f$ and $h_{H}=h_{L}/f^{2}$, as is the case presented here with $f=1.2$.

We found roughly the expected 4th-6th order convergence as displayed in Table 5 for the values of the recoil velocity and peak luminosity as well as the final black hole mass and spin (as measured on the apparent horizon.) The results of an extrapolation to infinite resolution and the differences with respect to the standard resolution (labeled as N144) are displayed in Table 5 to provide a measure of the expected errors for the whole family of simulations. Generically, for other simulations, we monitor the accuracy by measuring the conservation of the individual horizon masses and spins during evolution as well as the level of satisfaction of the Hamiltonian and momentum constraints. All eight N144 configurations show comparable behavior in these quantities.

IV Discussion

We compute the waveforms $a$ and $b$ matching as the inner product in frequency $f$-domain

[TABLE]

where the kth detector’s noise power spectrum is $S_{h,k}(f)$ and we adopt a low-frequency cutoff $f_{\rm min}$. By construction, we maximize over both a time and phase shift between waveforms. For our analysis of GW150914, we adopt the same noise power spectrum employed in previous work Abbott et al. (2016a); Lange et al. (2018), the advanced LIGO design sensitivity noise curve. We use a reference total mass of $M_{total}=74M_{\odot}$ and $f_{\rm min}=30Hz$. This choice of $M_{total}$ starts our waveform frequencies just below 30Hz after an initial windowing function is applied. The minimal SNR needed to distinguish between the two waveforms, given the mismatch is $\mathrm{SNR}^{2}\geq\frac{1}{1-\mathcal{M}}$.

To determine if waveforms from within this family of configurations can be distinguished between different members of the family, we perform a series of matches between configurations. That is, we choose a simulation and reconstruct the gravitational wave at a given polar and azimuthal angle and use this as our reference waveform. For each of the other configurations in the series, we can then calculate the match against our reference waveform and produce a “world map” of matches. We calculate the match

[TABLE]

where $i$ runs over each configuration, and where $\xi$ and $\psi$ are the angles used to reconstruct the second waveform at a given point in the skymap: $0\leq\xi\leq\pi$, and $-\pi<\psi\leq\pi$. In Fig. 5, we chose $\varphi_{ref}=291^{\circ}$ reconstructed at $\xi_{ref}=0^{\circ}=\psi_{ref}$ and calculate the SNR from the minimum, maximum, and mean matches over the world map. We show that the last few cycles of the gravitational waveform from black holes in the the hangup-kick configuration leading to a large recoil of the final remnant of the BBH merger is potentially measurable by LIGO with reasonable SNR, i.e. around approximately 30. For comparison, the matching between different resolutions of the reference case, $\varphi=291^{\circ}$, gives us SNR of the order of 96 and 25 for N120 and N100 respectively. Extrapolation to infinite resolution of the simulations $\text{N}_{\infty}$ leads to a SNR of over 100 in order to differentiate the N144 from the $\text{N}_{\infty}$ result.

Given the spin misalignments of comparable masses BBH observed in the current detections Abbott et al. (2018), these kind of configurations seems not so unlikely to occur in nature. While the search for detecting very highly spinning black holes with gravitational wave observations continues, it is important to search for them with the appropriated highly spinning templates and our simulations can contribute to fill in this gap near maximally spinning holes and properly cover this region of BBH parameter space. Parameter estimation techniques directly using numerical relativity waveforms from catalogs have been applied successfully for GW150914Healy et al. (2019) and GW170104Healy et al. (2018b) and will be the subject of further applications for O2 LIGO-Virgo observations.

Phenomenological modeling of waveforms, such as the PhenomP Schmidt et al. (2015) mimic precession from rotating aligned cases which leads to misevaluations of the recoil. See however new attempts to take recoil into account in other waveform models Chamberlain et al. (2019); Gerosa et al. (2018). An improved analysis of GW150914 using a two spins effective one body model is provided in Abbott et al. (2016b).

In Ref. Healy et al. (2019) we have been able to use a different family of simulations of binary black holes with one single spinning hole with amplitude $\alpha=0.8$ at all different orientations covering the two dimensional space of initial $(\theta,\varphi)$. Those lead to a “world heat map” as shown in the figure 8 of Healy et al. (2019) for the likelihood $\ln{\cal L}$ to represent the signal GW150914. Bit-equivalent data to the frames used for this study is available through GWOSC (Gravitational Wave Open Science Center) Vallisneri et al. (2015), and the likelihood, $\ln{\cal L}$, is calculated using the RIFT framework Pankow et al. (2015); Lange et al. (2018). In addition to this 3-parameter space estimation, we can consider the subfamily with the mass ratio $q$ and inclination angle $\theta$ leading to the highest likelihood $\ln{\cal L}$ and use this one remaining $\phi$-parametrized subfamily to parametrize the $\phi$-dependence of the recoil. The resulting “orbits” from the interpolation of the data are displayed in Fig. 6, showing the top three $\ln{\cal L}$ families and the preference for recoils of about $-1,500$ km/s.

Ultimately, determining accurately the recoil of the final hole from a binary system is paramount to determine (given a mass ratio) the spin orientations at merger. Being able to determine the “phase” of the spin relative to the linear momentum of the holes at the merger (as determined by the maximum amplitude of radiation) allows to predict the recoil of the remnant black hole. Such determination has been performed for GW150914 Healy et al. (2019) leading to estimated recoils of around 1,500 km/s as displayed in Fig. 6. The differences this induces on the merger and ringdown phases can be estimated as well, as a consistency check and a test of the theory of gravitation.

For the source of GW150914 we were also able to estimate the inclination of the orbit from purely numerical waveforms, as displayed in Figure 9 of Ref. Healy et al. (2019). The ability to find a single maximum, not bimodal, orientation of the binary, is somewhat related to the measure of precession and this in turn is related to the spin misalignment with the orbital angular momentum that may induce large recoil velocities, those depending on the merger phase that we model in this paper for the maximum recoil configurations.

The application of this techniques that we tested in the case of the first gravitational wave signal GW150914, can be used in other detections of BBH mergers, as GW170104 and others in O2 Abbott et al. (2018) and forthcoming observations and will be the subject of a future paper by the authors.

Acknowledgements.

The authors thank R. O’Shaughnessy and Y. Zlochower for discussions on this work and H. Pfeiffer for comments on the original manuscript. The authors gratefully acknowledge the National Science Foundation (NSF) for financial support from Grants No. PHY-1912632, No. PHY-1607520, No. PHY-1707946, No. ACI-1550436, No. AST-1516150, No. ACI-1516125, No. PHY-1726215. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) [allocation TG-PHY060027N], which is supported by NSF grant No. ACI-1548562. Computational resources were also provided by the NewHorizons, BlueSky Clusters, and Green Prairies at the Rochester Institute of Technology, which were supported by NSF grants No. PHY-0722703, No. DMS-0820923, No. AST-1028087, No. PHY-1229173, and No. PHY-1726215. Computational resources were also provided by the Blue Waters sustained-petascale computing NSF projects OAC-1811228, OAC-0832606, OAC-1238993, OAC-1516247 and OAC-1515969, OAC-0725070. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.