Tossing Black Hole Spin Axes

Gravitational wave (GW) astronomy began 7 yr ago with the discovery of the first binary black hole (BH+BH) merger (GW150914; Abbott et al. 2016). Since then, more than 85 such BH+BH merger events have been reported (Abbott et al. 2022a), besides mergers of double neutron stars (NS+NS) and mixed BH+NS events (Abbott et al. 2017a, 2021a). The merger events of BH/NS binaries are the most powerful energy sources known in the universe and bring new opportunities for insight into fundamental aspects of physics, including new tests of gravity in a highly relativistic regime (Abbott et al. 2021b), creation of energetic bursts of electromagnetic radiation (Abbott et al. 2017b; Ajello et al. 2018), and production of heavy chemical elements that decay and power an optical transient (Coulter et al. 2017; Drout et al. 2017; Smartt et al. 2017; Metzger 2019; Watson et al. 2019). Their origin remains controversial and has sparked a plethora of new research on massive binary star evolution in order to understand the formation process of such compact object mergers (e.g., Rodriguez et al. 2016a; Kruckow et al. 2016; Mandel & de Mink 2016; Marchant et al. 2016, 2019; Zevin et al. 2019; Bavera et al. 2020; Belczynski et al. 2020; Klencki et al. 2021).

The detected BH+BH mergers reveal individual BH component masses in the range from a few solar masses (M_⊙) to more than 100 M_⊙, thereby challenging both of the hypothesized mass gaps: at the lower and upper end of the BH mass spectrum, respectively (Bailyn et al. 1998; Renzo et al. 2020; Abbott et al. 2021d). The other fundamental parameter describing an astrophysical BH is spin. The dimensionless spin magnitude is given by χ = cJ/GM², where J is the BH spin angular momentum, M is its mass, and c and G are the speed of light and the gravitational constant, respectively.

Unfortunately, it is difficult to directly measure the magnitudes of the individual BH spins (χ₁ and χ₂), or their tilt (misalignment) angles (Θ_1,2) with respect to the orbital angular momentum vector of the binary ( L ) during the inspiral of the merger. Only the sum of their projected spins along the orbital angular momentum is well measurable, i.e., the effective inspiral spin

$$\begin{eqnarray}\begin{array}{r}\begin{array}{c}\begin{array}{ccc}{\chi }_{{\rm{e}}{\rm{f}}{\rm{f}}} & \equiv & \displaystyle \frac{\left({M}_{1}{{\boldsymbol{\chi }}}_{1}+{M}_{2}{{\boldsymbol{\chi }}}_{2}\right)}{{M}_{T}}\cdot \displaystyle \frac{{\boldsymbol{L}}}{| {\boldsymbol{L}}| }\\ & = & \displaystyle \frac{{\chi }_{1}\cos {{\rm{\Theta }}}_{1}+q\,{\chi }_{2}\cos {{\rm{\Theta }}}_{2}}{1+q},\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 1 }$$

where M_T = M₁ + M₂ is the total mass of the two BHs and q ≡ M₂/M₁ ≤ 1 is their mass ratio. Whereas misaligned BH spins are subject to precession, χ_eff is constant during the long inspiral prior to detection (Gerosa et al. 2015) and therefore an important diagnostic tool for the origin of BH+BH systems.

Figure 1 (top panel) shows the distribution of measured individual mean values of χ_eff of all BH+BH mergers and represents a conundrum due to the combination of two characteristics: (i) the distribution is asymmetric (Abbott et al. 2021c; Roulet et al. 2021; Zevin et al. 2021; Abbott et al. 2022b); and (ii) there is apparently a significant number of systems having χ_eff < 0 (although caution should be taken for the accuracy of individual values of χ_eff, see the Appendix). On the one hand, a perfectly symmetric distribution of χ_eff centered on χ_eff = 0 would be expected if all events originated from binaries that assemble pairs of BHs with random individual spin directions (e.g., via dynamical interactions in globular clusters, Mandel & O'Shaughnessy 2010; Rodriguez et al. 2016b). An asymmetric χ_eff distribution therefore points to, at least, a substantial contribution from BH+BH systems produced in isolated binaries. On the other hand, an isolated binary star origin traditionally predicts close-to-aligned BH component spins (Kalogera 2000; Farr et al. 2017) and according to the current school of thought, negative effective spins (requiring at least one misalignment angle component Θ_i > 90°) can only result if very large asymmetric momentum kicks are imparted onto BHs in their formation process (Wysocki et al. 2018; Callister et al. 2021a; Fragione et al. 2021), following the collapse of a massive stellar core—see Section 2.2.1. However, large BH kicks are incompatible with observations of Galactic BH binary systems (Mandel 2016; Mirabel 2017). The only other way to produce χ_eff < 0 is to prevent alignment of stellar spins prior to the formation of the second-born BH (see, however, Section 5.1).

Zoom In Zoom Out Reset image size

Finally, it should be mentioned that BH+BH mergers could also have an origin from primordial BHs, population III BH binaries, or dynamical assembly in active galactic nuclei disks (for a review, see, e.g., Mandel & Farmer 2022, and references therein). However, in those scenarios the expected χ_eff distribution is less clear. For further discussion on the interpretation of the empirical χ_eff data, see, e.g., Stevenson et al. (2017), Qin et al. (2018), Bavera et al. (2020), Olejak & Belczynski (2021), Callister et al. (2022).

Here a different mechanism is investigated in which the formation of a BH in a supernova (SN) is accompanied by the tossing of its spin axis in a random (isotropic) direction. There is unambiguous observational evidence from high-precision timing of binary radio pulsars (Breton et al. 2008; Desvignes et al. 2019) clearly demonstrating that such a process of spin-axis tossing is at work during the formation of NSs (Tauris et al. 2017). In fact, in two out of two NS+NS systems in which the young NS is observable as a radio pulsar, its misalignment angle has been measured to be very significant: 50° ± 04 and 104 ± 9°, respectively ¹ (Breton et al. 2008; Desvignes et al. 2019). Further evidence of tossing of the BH spin axis during the SN comes from recent measurements of a large spin–orbit misalignment (>40°) in the X-ray binary MAXI J1820+070 (Poutanen et al. 2022). Hence, although the process responsible for the tossing of the BH spin axis following stellar core collapse remains unsettled (Section 5.3), there are accumulating pieces of evidence that the process of spin-axis tossing is a generic process in all core-collapse SNe (although, as we shall discuss, the mechanism may be different for NSs and low-mass BHs, formed via fallback of SN ejecta, compared to massive BHs formed out of direct core collapse).

In Section 2, analytical SN kinematics is introduced for post-SN orbits without BH tossing, and sample results of BH+BH binaries are given. In Section 3, the applied method of Monte Carlo simulations with and without BH tossing is presented, along with a description of the chosen range of input parameter values and their probability density functions (PDFs). The results of these simulations, and a comparison to BH+BH observational data, are presented and discussed in Section 4. Additional discussions are given in Section 5 and our conclusions are summarized in Section 6.

2.1. Analytical Description

The dynamical effects of SNe in close binaries have been studied in detail both analytically and numerically in a number of papers (e.g., Flannery & van den Heuvel 1975; Hills 1983; Brandt & Podsiadlowski 1995; Kalogera 1996; Tauris & Takens 1998; Tauris & van den Heuvel 2023). Here we follow the prescriptions of Hills (1983) and Tauris & van den Heuvel (2023), and consider the explosion of a star (a stripped helium star of mass, M_He) in a circular binary, producing the second-born BH with a mass, M_BH,2. The companion star is in all cases assumed to be the first-born BH with a mass, M_BH,1. The change in the orbital semimajor axis (i.e., the ratio of final post-SN to initial pre-SN value) can be expressed by

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{{\rm{f}}}}{{a}_{{\rm{i}}}}=\left[\displaystyle \frac{1-{\rm{\Delta }}M/M}{1-2{\rm{\Delta }}M/M-{\left(w/{v}_{\mathrm{rel}}\right)}^{2}-2\cos \theta (w/{v}_{\mathrm{rel}})}\right],\end{eqnarray} \tag{ 2 }$$

where M = M_He + M_BH,1 is the total mass of the pre-SN system and ${v}_{\mathrm{rel}}=\sqrt{G({M}_{\mathrm{He}}+{M}_{\mathrm{BH},1})/{a}_{{\rm{i}}}}$ is the relative velocity between the two stars. The kick angle, θ is defined as the angle between the kick velocity vector, w , and the pre-SN orbital velocity vector of the exploding star, v _He, in the center-of-mass rest frame. Assuming a circular pre-SN binary is a good approximation for close binaries given the tidal interactions (Zahn 1977; Hut 1981) during Roche-lobe overflow (RLO) prior to the SN. Solving for the denominator being equal to zero in Equation (2) yields the critical angle, θ_crit, so that θ < θ_crit will result in a disruption of the orbit.

The eccentricity of the post-SN system can be evaluated directly from the post-SN orbital angular momentum, L_orb,f, and is given by

$$\begin{eqnarray}&&e=\sqrt{1+\displaystyle \frac{2\,{E}_{\mathrm{orb},{\rm{f}}}\,{L}_{\mathrm{orb},{\rm{f}}}^{2}}{{\mu }_{{\rm{f}}}\,{G}^{2}{M}_{\mathrm{BH},1}^{2}{M}_{\mathrm{BH},2}^{2}}},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{L}_{\mathrm{orb},{\rm{f}}}={a}_{{\rm{i}}}\,{\mu }_{{\rm{f}}}\,\sqrt{{\left({v}_{\mathrm{rel}}+w\cos \theta \right)}^{2}+{\left(w\sin \theta \sin \phi \right)}^{2}}.\end{eqnarray} \tag{ 4 }$$

Here μ_f and E_orb,f = −G M_BH,1 M_BH,2/2a_f are the post-SN reduced mass and orbital energy, respectively. For a definition of the kick angles, θ and ϕ, see Figure 13 in Tauris et al. (2017).

If the kick applied to the second-born BH is directed out of the orbital plane of the pre-SN system (i.e., if ϕ ≠ {0°, ±180°}), then the spin axis of the companion star (first-born BH) will be tilted with respect to post-SN orbital angular momentum vector. This effect gives rise to the geodetic precession seen in several radio pulsar binaries (in which there has been no mass transfer since the SN explosion) where the spin axis of the main-sequence star, white dwarf or NS companion is measured to be tilted (Kaspi et al. 1996; Kramer et al. 2006; Ferdman et al. 2013; Krishnan et al. 2020). This misalignment angle can be calculated as

$$\begin{eqnarray}&&\delta ={\cos }^{-1}\left(\displaystyle \frac{{v}_{\mathrm{rel}}+w\cos \theta }{\sqrt{{\left({v}_{\mathrm{rel}}+w\cos \theta \right)}^{2}+{\left(w\sin \theta \sin \phi \right)}^{2}}}\right).\end{eqnarray} \tag{ 5 }$$

If the misalignment angle is large (δ > 90°), then the new orientation of the orbit will cause retrograde spin of one or both binary components to the sense of orbital revolution. It is trivial to see that a second critical kick angle, ${\theta }_{\mathrm{retro}}={\cos }^{-1}(-{v}_{\mathrm{rel}}/w)$, exists such that post-SN orbits will be retrograde if θ > θ_retro > 90° (independent of ϕ, see below).

2.2. Analytical Results

Figures 2 and 3 show the misalignment angle, δ as a function of the two kick angles (θ and ϕ) resulting from the production of a 12 + 8 M_⊙ BH+BH system, originating from a progenitor 12 + 10 M_⊙ BH+He star binary with a fixed pre-SN orbital period of P_orb = 3.0 days and applying a fixed kick velocity of w = 500 km s⁻¹. The locations of the two critical angles, θ_crit = 1054 and θ_retro = 1458 can be clearly seen in both figures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows a variety of post-SN parameters calculated from the same progenitor system as in Figures 2 and 3, except that here the second kick angle is kept constant, ϕ = 60°, and the kick velocity, w, is variable. The color coding in each panel refers to either: misalignment angle (top left), orbital period (top right), eccentricity (bottom left), or merger time (bottom right). The systems above the black curved line have a first-born BH with a retrograde spin because the kick tilts the orbit to δ > 90°.

Zoom In Zoom Out Reset image size

Figure 5 shows examples of effective spin, χ_eff as a function of kick velocity, w, and kick angle, θ (color-coded curves from red to blue with increasing value of θ), imparted on the second-born BH for the same resulting 12 + 8 M_⊙ BH+BH system as in Figures 2–4. The four different panels are calculated for different pre-SN orbital periods, P_orb,i, and assumptions of fixed BH component spins, χ_1,2.

Zoom In Zoom Out Reset image size

2.2.1. Negative Effective Spins from Large Kicks

An important finding from Figure 5 is that producing tight systems with χ_eff < 0 that become BH+BH mergers within a Hubble time requires a kick of w > v_rel, often corresponding to w > 400 ∼ 600 km s⁻¹. (Notice that here retrograde systems are equivalent to χ_eff < 0 because Θ₂ ≃ Θ₁ = δ, see Section 2.2.2.) It can thus be concluded that if BH spin-axis tossing is not at work then, to produce BH+BH mergers with χ_eff < 0 in isolated systems, either very large secondary BH kicks are needed or pre-SN BH and stellar spins are not aligned—two conditions that are both highly questionable (Section 5.1).

To better demonstrate the effects on post-SN orbits and misalignment angles, Θ₁ = δ of the first-born BHs due to kicks, Figures 2–5 display the results based on a relatively low-mass BH+BH system (M_BH,1 = 12.0 M_⊙, M_He = 10.0 M_⊙, and M_BH,2 = 8.0 M_⊙). Even in that case, it is seen that a minimum kick velocity of w > 400 km s⁻¹ (in a backward direction, θ → 180°) is needed to obtain a retrograde post-SN orbit (i.e., Θ₁ = δ > 90° and thus ${\chi }_{1}\cos {{\rm{\Theta }}}_{1}\lt 0$ and χ_eff < 0).

2.2.2. Spin Directions of Progenitor Systems

For isolated systems, the tilt angles of the two BH spins (Θ₁ and Θ₂, see Equation (1)) are assumed to be related to δ as follows. The tilt angle of the spin axis of the first-born BH (${{\boldsymbol{\chi }}}_{1}$), by the time the system merges, is equal to the resulting misalignment angle, δ from the second SN in the BH+BH system, i.e., Θ₁ = δ. This inference follows the simple argument (e.g., Hut 1981; Hills 1983; Bhattacharya & van den Heuvel 1991; Guillemot & Tauris 2014; Biryukov & Abolmasov 2021; Tauris & van den Heuvel 2023) that at the time of the second SN, the first-born BH is assumed to possess prograde spin with its spin vector (${{\boldsymbol{\chi }}}_{1}$) parallel to the pre-SN orbital angular momentum vector as a result of previous mass-transfer episodes from the progenitor of the collapsing stripped helium star onto the first-born BH (Section 5.1).

For the spin tilt angle of the second-born BH, we consider two cases: In the case of no tossing: Θ₂ ≃ Θ₁ = δ. This follows from the assumption that tidal torques will align the spin axis of the collapsing helium star with that of the orbital angular momentum vector prior to the SN. However, if the secondary BH is subject to tossing of its spin axis during its birth in the SN after core collapse, then Θ₂ will follow a range of different (random) values (Section 3.1.5).

3.1. Method and Modeling

After introducing the kinematic effects of the second SN on a sample system in Section 2, we now investigate the outcome on a large BH+BH population as a function of various input distributions using Monte Carlo simulations. The assumptions about the input PDFs are kept relatively simple to better capture the effect of spin-axis tossing—the main scope of this paper. The outlined method can readily be applied to input PDFs obtained from population synthesis.

Here again, the starting point is a population of pre-SN binaries composed of the first-born BH and a companion star (the progenitor of the second-born BH) on the verge of core collapse. In the following, we first discuss the input distributions (PDFs) and refer to sample figures showing the outcome of the simulations along the way. A more systematic description of the displayed results is given in Section 4 and discussed further in Section 5.

3.1.1. Initial Masses

The masses of the first-born BH and the collapsing helium star companion were either chosen at fixed values or from a range of masses recovering most of the BH+BH events in the GWTC data (Abbott et al. 2022b). In the majority of simulations demonstrated here (most panels in Figures 7 and 9–11), the mass component values were simply chosen from a flat PDF such that the resulting BH masses are in the range M_BH,1 = 8–70 M_⊙ and M_BH,2 = 3–45 M_⊙, thereby reproducing the bulk of detected BH+BH mergers (see Figure 1 for GWTC data). The distribution of masses of the secondary (collapsing) stars are thus chosen in the range between M_He = 3.75 and 56.25 M_⊙ because it is assumed here that 20% of the mass of the collapsing star is lost during its collapse to a BH. This fraction may be on the high side. However, as the results are only weakly dependent on the BH mass distribution, using a mass-loss fraction of, e.g., 5% yields very similar results. For the same reason, the exact shape of the overall BH mass PDF is not very important here and it justifies the simple use of a flat PDF for the BH masses.

In the case of investigating BH spin-axis tossing effects for fixed BH masses, default pre-SN values of 37 and 30 M_⊙ are chosen for the first-born BH and the collapsing helium star, respectively (e.g., Figures 8(A), (B), and (E)), to produce final BH+BH mergers with mass components of M_BH,1 = 37 M_⊙ and M_BH,2 = 24 M_⊙. Simulations with other fixed BH mass values, including more extreme mass ratios, are also shown in Figure 8.

Finally, simulations were also performed assuming that ome (or all) of the BH+BH systems have their mass components inverted (i.e., M_BH,1 = 24 M_⊙ and M_BH,2 = 37 M_⊙, Figure 8(F)), to reflect an extreme extrapolation of a subpopulation of BH+BH progenitor systems undergoing mass reversal during their evolution—such that the initially least massive (i.e., secondary star) of the two ZAMS stars accretes sufficient mass to end up producing the first-born BH; and the initial primary star then produces the second-born BH. For the general process of mass reversal and SN order in close binaries, see, e.g., Wellstein & Langer (1999), Tauris & Sennels (2000), and Zapartas et al. (2017).

3.1.2. Initial Orbital Separation

The pre-SN orbital separation, a₀ is chosen from a flat PDF between 4 < a₀ < 40 R_⊙. The lower limit is set to fit the size of the collapsing helium star within its Roche lobe, and the upper limit is chosen such that the majority of the surviving post-SN systems will remain in orbits tight enough to let the BH+BH system merge within a Hubble time (∼13.8 Gyr). As demonstrated in the top panels of Figure 8, the value of a₀ is not important for the distribution of χ_eff (unless BH kicks are extreme and pre-SN orbits are very wide)—only the number of systems merging within a Hubble time decreases with increasing a₀ (compare panels (A) and (B)).

3.1.3. BH Kick

The kick imparted onto the newborn BH is assumed to have a random (isotropic) direction with θ ∈ [0; 180°] and ϕ ∈ [0; 360°]. This is obtained with a PDF of $P(\theta )\,\,=0.5\sin (\theta )$ and a flat PDF for ϕ. In most cases, a kick magnitude of w = 50 km s⁻¹ is applied. This modest kick value is chosen to match constraints from BHs in Galactic X-ray binaries (Mandel 2016; Mirabel 2017) and also because the main emphasis here is to demonstrate the effect of the tossing of the BH spin axis. For comparison, however, simulations were also performed using w = 500 km s⁻¹, and in a few cases without any kick (w = 0, top row in Figure 11). It is found, perhaps somewhat surprisingly, that the distribution of χ_eff does not depend much on w (e.g., in Figure 7, compare panels (C) and (D) with panels (E) and (F)). The main effect is simply that larger SN kicks destroy more binaries. For this reason, and to avoid too many free parameters, the default value is set to w = 50 km s⁻¹. The reason for the (almost) kick invariance is that for a large part of the BH+BH progenitor systems the applied kick is still relatively small compared to the pre-SN orbital velocity (w/v_rel ≪ 1).

3.1.4. Spin Magnitudes

The spin magnitudes of the two BH components (χ₁ and χ₂) are drawn from four different PDFs plotted in Figure 6 (see details in Table 1), based on different expectation values for the first- and second-born BH in an isolated binary system (see Section 5.1 for further discussion). Whereas the first-born BHs have been argued to be preferentially slowly spinning (χ₁ ≲ 0.1, Fuller & Ma 2019; Bavera et al. 2020), the second-born BHs are expected to spin much faster due to tidal interactions between the first-born BHs and the helium stars that collapse and produce the second-born BHs (Kushnir et al. 2016; Hotokezaka & Piran 2017; Zaldarriaga et al. 2018).

Zoom In Zoom Out Reset image size

Table 1. Normalized PDFs of BH Component Spins, χ_i (i = 1, 2) for Simulating χ_eff of BH+BH Systems That Merge

Name	PDF	Shape	Average	Median
planck1	$\displaystyle P({\chi }_{i})=\displaystyle \frac{\text{125,000}\,{\chi }_{i}^{3}}{\exp (30\,{\chi }_{i})-1}$	Strong peak near χi ∼ 0.1	〈χi 〉 = 0.128	χi,med = 0.117
power2	$\displaystyle P({\chi }_{i})=\displaystyle \frac{17.4{\left({\chi }_{i}/0.2\right)}^{0.8}}{{\left(1+{\chi }_{i}/0.2\right)}^{3.2}}$	Concentrated at small χi	〈χi 〉 = 0.288	χi,med = 0.214
planck2	$P({\chi }_{i})=\displaystyle \frac{656\,{\chi }_{i}^{3}}{\exp (8\,{\chi }_{i})-1}$	Broad peak χi ∼ 0.2 − 0.6	〈χi 〉 = 0.450	χi,med = 0.426
high	$P({\chi }_{i})=\left(\displaystyle \frac{24}{7}\right)\displaystyle \frac{1}{{\left(2-{\chi }_{i}\right)}^{4}}$	Increasing with χi	〈χi 〉 = 0.714	χi,med = 0.789

Note. See plots of the PDFs in Figure 6 and further discussion in Section 5.1.

Download table as:  ASCIITypeset image

The two distributions "planck1" and "planck2" have median χ values of χ₁ ∼ 0.1 and χ₂ ∼ 0.4, respectively, and were selected to mimic the empirical distributions inferred in Figure 17 in Abbott et al. (2022b) based on GWTC data of BH+BH mergers. These are the distributions applied in the default results shown in Figures 7(A) and (B). The effects on the χ_eff distribution, however, were primarily investigated here using the alternative PDFs "power2" and "high" for χ₁ and χ₂, respectively (Figures 7(C)–(F), and Figure 8) because these distributions likely better reflect expectations from close binary star evolution (Section 5.1). Finally, the effects on χ_eff by applying fixed or a narrow range of values of χ₁ and/or χ₂ were also investigated, see Figures 9–11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.5. BH Spin-axis Tossing

The PDF for the degree of tossing of the spin axis of the second-born BH is assumed to follow an isotropic distribution of post-SN spin axis directions: $P({{\rm{\Theta }}}_{2})=0.5\sin ({{\rm{\Theta }}}_{2})$, where Θ₂ ∈ [0; 180°]. In the case of a core collapse without the effect of BH spin-axis tossing, it is assumed that Θ₂ = Θ₁ = δ (see Section 2.2.2).

A number of histograms are presented in Figures 7–11 to display the resulting χ_eff distributions of BH+BH systems produced from Monte Carlo simulations with (green bins) and without (red bins) BH tossing, following the methodology and details outlined in Sections 2 and 3. The total number of simulated pre-SN systems is quoted in each panel of all figures. Underneath is stated the number of post-SN systems which survive and produce a BH+BH merger within a Hubble time (GW systems). A number of post-SN systems are disrupted (also stated in the legend) if the applied SN kick value is large (w = 500 km s⁻¹). The rest of the surviving systems end up in BH+BH binaries that are too wide to merge within a Hubble time. (Very few systems ${ \mathcal O }(\lt {10}^{-5})$ actually produce direct mergers if the second-born BH is kicked hard directly toward the first-born BH such that their closest spatial separation is less than the sum of their event horizon radii.)

Notice that the total number of systems simulated in each panel is simply adjusted to optimize the height of the most frequent bins in the histograms. Similarly, the plotted empirical GWTC-3 data bins have been multiplied with a constant scale factor of 100 for clarity. Some panels include solid curves that represent (blue) the mean estimate of the PDF of χ_eff from the GWTC-3 data (Abbott et al. 2022b), (gray) the 5% and 95% quantiles, and (black) the blue curve normalized to the number of simulated systems. An example is shown in Figure 7(A), which may be considered the default result for BH spin-axis tossing. Note, the black curve is not a fit to the simulated distribution shown as the green histogram—it is the normalized empirical data (but the agreement is remarkable).

The main result of the investigation in this work is that the observed distribution of χ_eff cannot be reproduced from isolated binaries without tossing of the BH spin axis (as clearly demonstrated in Figures 7(B), (D), and (F) and Figures 9 and 11 (panels (B), (D), and (F)), but see also Section 2.2.1. On the other hand, simulations including tossing of the BH spin axis reproduce the observed distribution of χ_eff quite well—even for a wide range of input parameter values. The resulting χ_eff distribution is basically invariant to SN kicks, exact pre-SN orbital separations, and only weakly dependent on the input distributions of BH masses. This is demonstrated by comparing simulation results, e.g., in Figures 7(C) and (E), and Figures 8 and 11 (panel (A)), as will be discussed further below in more detail. The few observed BH+BH systems with very large values of χ_eff ∼ 0.7 are somewhat outliers and may be recovered using different PDFs for the BH spins (see Figure 7(C)), or these systems may simply have a distinct origin (Section 5.2). We now discuss the results of each of the histograms.

Figure 7 directly compares the resulting χ_eff distribution including (left column) and excluding (right column) the process of BH spin-axis tossing. Here it is clearly demonstrated that only by including BH tossing the observed χ_eff distribution can be reproduced. As an example, panel (A) shows a simulation of the expected distribution of χ_eff following a scenario with BH tossing. The plot is based on 25,000 computations of the second SN, from which 22,254 BH+BH systems are produced that will merge within a Hubble time. The agreement with the observed data (blue bins and black curve PDF) is striking—see also Figure 11(A). Figure 7(C) is similar to panel (A), except that other spin distributions (power2 and high, see Table 1) are applied to enable production of slightly larger χ_eff values. Panels (B) and (D) show similar simulations without BH tossing, which irrefutably cannot reproduce the observational data. The bottom row shows that the χ_eff distribution is basically invariant to BH kick magnitude (Section 3.1.3). Unless absurdly large, applying a kick of e.g., w = 500 km s⁻¹ only increases the number of systems disrupted due to the SN.

Figure 8 assumes fixed BH masses. In panels (A), (B), and (E), the BH masses are equal to the average BH masses of M₁ = 37 M_⊙ and M₂ = 24 M_⊙ from GWTC-3, see Figure 1) and the simulations here show that the resulting χ_eff distribution to a large degree is invariant to pre-SN orbital separation, a₀. In panel (A), a constant pre-SN orbital separation value of a₀ = 4 R_⊙ is applied, whereas a₀ = 40 R_⊙ is applied in panel (B). Despite a factor of 10 in difference in the choice of a₀, the only major difference in the outcome is that fewer systems merge within a Hubble time if a₀ = 40 R_⊙, as expected. (Notice that no post-SN systems will merge within a Hubble time if a₀ > 54 R_⊙ and w = 50 km s⁻¹.) Panels (C)–(F) show that the χ_eff distribution depends only somewhat weakly on the exact choice of BH masses.

Figure 9, and panels (B), (D), and (F) of both Figures 7 and 11, clearly demonstrate that reproducing the observed χ_eff distribution without BH spin-axis tossing is impossible, no matter the choice of input distributions of BH masses, their spins and applied kicks (but see the caveat in Section 2.2.1). To boost the kinematic effect of SNe, in terms of producing small or negative values of χ_eff, very large kicks of w = 500 km s⁻¹ were applied in all panels in Figure 9. Despite this, the simulated data which neglects BH spin-axis tossing cannot reproduce the observed data.

4.1. Spin of First-born BH

If the far majority of the observed BH+BH systems do have an origin from isolated binaries (as supported by simulations presented here), a strong mismatch between observations and simulations is obtained in cases where the first-born BH is assumed to be fairly rapidly spinning (see, e.g., Figures 10(E) and (F)). This result is in agreement with expectations from binary stellar evolution theory (Heger et al. 2005; Fuller & Ma 2019) predicting a small spin for the first-born BH (χ₁ ≲ 0.1). Indeed, it is demonstrated here (Figure 7(A) and Figure 11(A)), both applying planck1 and planck2 spin PDFs for χ₁ and χ₂, see Table 1 and Figure 6) that individual spin magnitudes of order χ₁ ∼ 0.1 and χ₂ ∼ 0.4 reproduce the empirical data quite well (see Section 5.1 for further discussion). This is in agreement with a statistical analysis of the currently measured data (see Abbott et al. 2022b, and their Figure 17). It is also seen from the empirical data that the average sum of the two spin magnitudes should not be too large (∣χ₁∣ + ∣χ₂∣ ≲ 0.8).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Obviously, one can obtain desired final χ_eff values in a narrow band by selecting χ₁ and χ₂ values in a narrow interval and applying small kicks (see, e.g., Figures 11(C) and (D)).

5.1. BH Component Spins and Alignment

The spin magnitude of a newborn BH reflects the spin of the core of its collapsing progenitor star, which itself depends on the angular momentum transport during its previous phases of stellar evolution, including stellar winds and accretion. Efficient angular momentum transport by viscosity and magnetic torques will couple the stellar core to its envelope, thereby slowing the spin of the core as the envelope expands when it becomes a giant star. Based on such studies, it has been argued (Heger et al. 2005; Fuller & Ma 2019) that BHs from single stars, and first-born BHs in binaries, are expected to have very slow spins. This motivates the green (planck1) and blue (power2) PDFs of χ₁ shown in Figure 6 (see Table 1 for details).

The rationale for allowing the possibility χ₁ → 1 is twofold: first, the estimated spins of BHs observed in HMXBs are apparently all quite large (Miller-Jones et al. 2021) with values between 0.84 and 0.99; and, second, some of the most massive BH components (and thus in general expected to have formed first) in recorded BH+BH mergers seem to possess rather large values of χ₁ (Abbott et al. 2022b, with the usual caveat that these measurements have large uncertainties in general). This evidence suggests that either (i) isolated binary evolution may produce larger spins of the first-born BHs than predicted from current theory, (ii) the BHs in these HMXBs and BH+BH merger sources could possibly have obtained their rapid spins from chemical homogeneous evolution (Mandel & de Mink 2016; Marchant et al. 2016), or (iii) the general populations of (near-)Galactic HMXBs with rapid BH spins and BH+BH mergers are distinct (Gallegos-Garcia et al. 2022). The rapidly spinning BHs in the former population may give support to the third possibility and could therefore conveniently alleviate the problem of aligning rapidly spinning BHs (see below), following a core collapse with spin-axis tossing.

Rapid BH spins obtained from Eddington-limited accretion are not possible within the short lifetime of their massive companion stars. From similar simple arguments as in, e.g., Mandel & Fragos (2020), one can show that the ratio of orbital angular momentum of an accreted test particle (with mass m, carrying the specific orbital angular momentum from the innermost stable circular orbit (ISCO)) and the spin angular momentum of a slow-spinning accreting BH (with mass M) is given by

$$\begin{eqnarray}&&\displaystyle \frac{J}{{J}_{\mathrm{BH}}}\simeq \displaystyle \frac{1}{\chi }\ \displaystyle \frac{m}{M}\sqrt{\displaystyle \frac{{R}_{\mathrm{ISCO}}\,{c}^{2}}{{GM}}}\sim \displaystyle \frac{\sqrt{12}}{\chi }\ \displaystyle \frac{m}{M}.\end{eqnarray} \tag{ 6 }$$

Hence, for χ ≃ 0.1, we see that J/J_BH ∼ 35 m/M, and thus accreting an amount of mass equivalent to a few percent of the BH can make a significant contribution to J_BH, and therefore likely enable alignment (at least for slow-spinning BHs).

Moreover, and more importantly, tides (Hut 1981) and torques from BH–disk interactions (Bardeen–Petterson effect, Bardeen & Petterson 1975; Natarajan & Pringle 1998; King et al. 2005) have been shown to align the BH spin with the orbital angular momentum vector (but see also Stegmann & Antonini 2021). Finally, the subsequent common-envelope stage (central to most models of isolated binary star evolution) is poorly understood, and may also affect spin–orbit couplings of angular momenta from the bulk mass ejection on a short timescale.

Figures 9(E) and (F) were simulated with partial tossing, i.e., BH spin-axis tossing was implemented only for systems where M_BH,2 < 20 M_⊙ or χ₂ < 0.30 (see Section 5.3 for a discussion). The resulting χ_eff distributions do not match well with observations. This is not surprising since, in effect, these input conditions imply that a fairly large fraction of the simulations did not include tossing and therefore fails to reproduce the empirical data (see Figure 7(B)).

Figure 10 shows the effect of applying various BH spin magnitudes, in all cases including BH spin-axis tossing. Simulations with a constant large value of χ₁ = 0.8 (panels (E) and (F)) resulted in χ_eff distributions that completely fail to explain the empirical data, thereby favoring smaller values of χ₁ in general, if these systems formed from isolated binary star evolution. Panels (A) and (C) serve as sanity checks: panel (A) assumes χ₁ = 0 (i.e., the first-born BH is non-spinning at the time of the second core collapse), and panel (C) shows the effect of double BH tossing (i.e., equivalent to the situation in which no alignment occurs for the spin axis of the first-born BH during mass accretion from the companion star prior to the second SN). In both of these cases, we would expect the distribution of χ_eff to be completely symmetric with respect to χ_eff = 0. This is indeed in agreement with the presented simulated data, but in disagreement with observations, which is therefore another piece of evidence in favor of spin alignment prior to the second SN.

The spin of the second-born BH in an isolated binary system is also determined by angular momentum transport, stellar evolution, and mass loss of its progenitor star. Particularly important for the second BH are tidal interactions in the post-HMXB (post-common-envelope or post-RLO) phase where a naked helium star is orbiting the first-born BH in a tight orbit. Using analytical arguments and semi-analytical calculations, one can derive approximate estimates of the spin-up of the secondary naked core prior to its collapse and thereby constrain the spin of the second-born BH (van den Heuvel & Yoon 2007; Kushnir et al. 2016; Hotokezaka & Piran 2017; Zaldarriaga et al. 2018). To mimic efficient spin-up via tides, PDFs high and planck2 (see Table 1) were applied for χ₂ in the work presented here.

5.2. BH+BH Outliers

5.3. BH Spin-axis Tossing Mechanism

In this work, strong evidence has been presented for the first time for BH spin-axis tossing to operate during the formation of BHs in SNe. This finding is based on a comparison of simulated distributions of χ_eff of BH+BH systems produced in isolated binaries with empirical data of BH+BH mergers from GWTC-3 of the LIGO-Virgo-KAGRA network (up to and including O3a and O3b). It is noteworthy that the theoretical simulations of the χ_eff distribution are only weakly dependent on a large range of input parameters, which further strengthens the evidence for BH spin-axis tossing. Future comparisons to empirical data should be followed up by further detailed and advanced statistical analysis, including selection effects, to rigidly test and confirm (or refute) the hypothesis presented here. The conclusions of this paper can be summarized as follows:

• 1.  
  Empirical data (GWTC-3) strongly disfavor the main origin of BH+BH mergers from isolated binary star evolution without BH spin-axis tossing. Including (isotropic) BH tossing, however, provides a simple solution and an excellent fit to the bulk data for a large range of input distributions of relevant physical parameters (e.g., Figure 7(A)).
• 2.  
  An implication of this hypothesis is that it jeopardizes previous ideas on the use of BH spin data as a diagnostic for double BH formation models.
• 3.  
  The simulations presented here are broadly in agreement with empirically derived (Abbott et al. 2022b) spin magnitudes with average values of 〈χ₁〉 ∼ 0.1 and 〈χ₂〉 ∼ 0.4 for the first- and second-born BHs, respectively (see PDFs planck1 and planck2 in Figure 6). However, we find that there is some indication that PDFs with a small tail extending up to χ₁ → 1, and increasing probability for larger χ₂ values as expected from binary star interactions (PDFs power2 and high, respectively), may fit the data even better.
• 4.  
  More statistical investigations (including selection effects) are needed to firmly verify and quantify the effect of BH spin-axis tossing, as well as for deriving constraints on the spin magnitudes of the first- and second-born BHs more robustly.
• 5.  
  It remains a riddle to explain the physical mechanism of the tossing itself—a process that may be generic to all core-collapse SNe (NSs and BHs).
• 6.  
  Producing BH+BH mergers with χ_eff < 0 from isolated binaries without BH tossing requires a BH SN kick of order w ∼ 500 km s⁻¹, in strong disagreement with observations of BHs in X-ray binaries. This finding is in agreement with previous work (e.g., Wysocki et al. 2018).
• 7.  
  The crucial question of the efficiency of alignment of the spin axes prior to the second SN is still an unresolved issue.

T.M.T. cordially thanks Thomas Janka for numerous discussions and explanations, Tom Callister for providing the PDFs of χ_eff data from GWTC-3, Aleksandra Olejak, Simon Stevenson, Norbert Langer, Norbert Wex, and Alejandro Vigna-Gómez for discussions, and Malte Klockmann Pein Malmberg for discussions during his BSc project on SN kinematics. Finally, acknowledgment is given to the two anonymous reviewers and the scientific editor for ApJ.

The observational data used in this work for comparison with theory were taken from the recent GWTC-3 catalog (Abbott et al. 2022a, 2022b). ² The merger events, which are double NSs (NS+NS: GW170817 and GW190425), putative mixed systems (BH+NS: GW190426_152155, GW190917, GW191219, and GW200115) or BH+BH binaries with a secondary component mass that was most likely less than 3.0 M_⊙ (GW190814 and GW200210), were removed from the sample. This leaves a total of 85 BH+BH merger events. Since the uncertainty of the measurements of χ_eff covers the range from ∣Δχ_eff∣ < 0.1 (GW191204) to ∣Δχ_eff∣ > 0.4 (GW200322), and typically ∣Δχ_eff∣ = 0.1 − 0.25, the individual values of χ_eff should be taken with caution, whereas the shape of the χ_eff distribution for the general population is reliable. Selection effects (Gerosa et al. 2018) from masses or spins, or internal relations (Callister et al. 2021b), are not considered here and are not expected to influence the main conclusions of this work, but see e.g., Mandel & Smith (2021) for discussions on priors used in the LIGO-Virgo-KAGRA analysis and the correlation between spin and mass ratio.