Black Hole Formation in Fallback Supernova and the Spins of LIGO Sources

Here we investigate within the context of field binary progenitors how the the spin of LIGO sources vary when the helium star-descendent black hole (BH) is formed in a failed supernova (SN) explosion rather than by direct collapse. To this end, we make use of 3d hydrodynamical simulations of fallback supernova in close binary systems with properties designed to emulate LIGO sources. By systematically varying the explosion energy and the binary properties, we are able to explore the effects that the companion has on redistributing the angular momentum of the system. We find that, unlike the mass, the spin of the newly formed BH varies only slightly with the currently theoretically unconstrained energy of the SN and is primarily determined by the initial binary separation. In contrast, variations in the initial binary separation yield sizable changes on the resultant effective spin of the system. This implies that the formation pathways of LIGO sources leading to a particular effective spin might be far less restrictive than the standard direct collapse scenario suggests.


INTRODUCTION
The gravitational wave (GW) signals detected by LIGO (Abbott et al. 2016b(Abbott et al. ,c, 2017a have uncovered a population of black holes (BHs) that is significantly more massive than the population known to reside in accreting binaries (Remillard & McClintock 2006). While there is significant debate in the community about how black hole binaries are assembled (Zwart & McMillan 2000;Kalogera et al. 2007;Sadowski et al. 2008;Postnov & Yungelson 2014;Abbott et al. 2016a;Belczynski et al. 2016;Rodriguez et al. 2016a;de Mink & Mandel 2016;Gerosa & Berti 2017;Wysocki et al. 2017), the classical scenario (Tutukov & Yungelson 1993;Voss & Tauris 2003) remains one of the leading candidates. In this channel, a wide massive binary undergoes a series of mass transfer episodes leading to a tight binary comprised of a massive helium star (M * ) and a BH (M 1 ), prior to the formation of the second BH (M 2 ). LIGO observations of the mass-weighted angular momentum perpendicular to the orbital plane χ eff , have been argued to provide constraints on this formation channel (Rodriguez et al. 2016b;Farr et al. 2017;Stevenson et al. 2017). This is because vital information on the mass transfer history of the binary and the spin of M * is imprinted on Here a 1 and a 2 are the dimensionless spins of the BHs andL is the direction of the angular momentum in the orbital plane. The angular momentum of the secondary BH is intimately linked to that of the progenitor helium star, which in turn is determined by its mass-loss history and the torque exerted by the primary black hole (Qin et al. 2018;Kushnir et al. 2016;Zaldarriaga et al. 2018). This torque can effectively drive synchronization of the stellar spin and the orbit in binaries tighter than d τ , the maximum separation allowed for synchronization within the life of the helium star (Zaldarriaga et al. 2018). The final angular momentum of the star thus provides a reasonable estimate of a 2 when the mass and angular momentum losses from the final supernova explosion are ignored. All previous works have assumed that such effects are small based on the simple expectation that LIGO BHs are formed by direct collapse.
Motivated by the fact that the nature of BH-forming supernova (SN) explosions is not to well understood (Fryer et al. 2012;Ugliano et al. 2012;Pejcha & Thompson 2015), in this Letter we explore the effects on χ eff when the second BH M 2 instead is formed by a fallback SN explosion (Moriya et al. 2010;Fryer et al. 2012;Dexter & Kasen 2013;Lovegrove & Woosley 2013;Perna et al. 2014;Batta et al. 2017;Fernández et al. 2018). For this purpose we make use of 3d hydrodynamical simulations of fallback SN in close binary systems with properties aimed at reproducing LIGO GW signals. The structure of this Letter is as follows. In Section 2 we describe the numerical formalism used to initiate the fallback SN explosion and compute the subsequent evolution of the binary. In Section 3 we describe the dynamics of the fallback material and its effect on the final spin of both BHs. Lastly, in Section 4 we present our key findings and relate them to the current population of LIGO sources.
The velocity profile of the ejecta at d * = 1 (soon after the shock emerges from the stellar envelope) is plotted in units of the escape velocity for three different initial explosion energies.The shaded region shows the initial BH mass. The subsequent growth of the BH depends on the amount of fallback material, which in turn, depends strongly on the energy injected. The dashed lines show the current mass of the BH for different explosion energies.

METHODS AND INITIAL SETUP
Here we follow the setup described in Batta et al. (2017) to study the evolution of the progenitor binary system after the birth of M 2 . We make use of a modified version of the three-dimensional smoothed particle hydrodynamics (SPH) code GADGET2 (Springel 2005), with our initial setup consisting of a tidally locked 28M helium star in orbit around a BH of M 1 = 15M . The reader is referred to Batta et al. (2017) for further details on initial particle distribution and numerical accretion. We settled for a resolution of 5 × 10 5 particles, which showed convergence for the accretion rates and properly captured the dynamics of the ejecta and the binary system.
To study the interaction of the fallback material with the newly formed BH binary we explore three sets of simulations. Each set starts with the initial binary in a circular orbit with a separation d = d * R * , where d * = 2, 3 or 5 and R * is the pre-SN star's radius. Then for each orbital separation we run simulations with at least four different SN explosion energies. In all cases we assume a 1 = 0 based on the results of Qin et al. (2018), and assume that synchronization of the stellar spin and the orbit has taken place, as is expected for the initial separations used in this analysis. Given the large uncertainties in BH natal kick estimations (Mandel 2016;Repetto & Nelemans 2015), we assumed the simplest scenario where no natal kick is applied to the recently formed BH. This combined with the synchronization of the stellar spin and the orbit, translates into BH's spins aligned with the orbital angular momentum.
The initial profile of the star was obtained from the 35OC KEPLER model calculated by Woosley & Heger (2006) of a 28M pre-SN helium star with R * = 0.76R . We considered the innermost 3M of the pre-SN star to be the newly formed BH with a 2 (t = 0) = 0, which we subsequently treat as a sink particle. After the removal of the inner core, we use a parameterized energy injection routine to mimic the supernova engine and derive the density and velocity profile of the expanding envelope. Specifically, we use a spherically symmetric kinetic piston at the inner boundary to adjust the energy injected into the envelope. In all calculations the energy is parametrized as follows: E SN = ζ E G , where E G = 2.3 × 10 52 erg is the binding energy of the pre-SN star.
The distribution that describe the ejecta is determined solely by the structure of the pre-SN star (Woosley & Weaver 1995;Matzner & McKee 1999) and is established by the hydrodynamics of the interaction. Initially, the shock propagates through the stellar material, pressurizing it and setting it into motion. Once the shock wave approaches the surface of the star, a rarefaction stage begins in which stellar material is accelerated by the entropy deposited by the shock. This stage terminates once the pressure ceases to be dynamically important and the material expands freely. Figure 1 shows the radial velocity profile of the envelope when the shock surfaces the stellar envelope for three different explosion energies: ζ = 0.1, 0.5 and 0.9. The gray area shows the initial BH mass while the dashed lines show M 2 (t) at the time the shock reaches r = R * . Despite the complicated hydrodynamical interaction, the density and pressure approach steep power laws in velocity as the material approaches homologous expansion (Figure 1). For a given progenitor structure, the ensuing ejecta will take similar distributions.
When ζ 1, energy injection fails to unbind the star such that a sizable fraction of its mass eventually fallbacks onto the newly formed BH. If this takes place in a binary system (Batta et al. 2017), a non-negligible fraction of the bound material can expand to a radius comparable or larger than the binary's separation, thus immersing the BH companion in gas. The interaction of fallback material with the binary transfers orbital angular momentum to the gas, which upon accretion onto the orbiting BHs is finally transferred into spin angular momentum. Differences in ζ result in diverse accretion histories, which ultimately regulate the BHs' final masses and spins. It is to this topic that we now draw our attention.   Flow dynamics are similar in all three simulations shown in Figure 2. First, the envelope expands to rapidly engulf the companion BH. A bow shock is created as a result of this initial interaction. It is, however, only when the slower moving material reaches the companion that the resulting torque can supply the envelope gas with sizable angular momentum. This envelope material will remain bound to the system and will form a disk around M 2 if restricted to the region within which orbiting gas is gravitationally bound to the newly formed BH. A disk, albeit lighter, also forms around M 1 , whose final mass depends sensitively on the initial separation.
The total mass bound to M 2 is the same in all simulations, yet the fraction of angular momentum accreted increases with decreasing separation. As a result, a 2 is higher for progressively more compact binaries despite the final mass of M 2 reaching similar values. This can be seen in the left panel of Figure 3, in which we show the evolution of a 2 as a function of the accreted mass in units of M * . Initially, a 2 increases as envelope material is accreted directly onto the BH. The innate angular momentum in this initial phase is determined by tidal synchronization, which increases as the binary separation decreases (Kushnir et al. 2016(Kushnir et al. , 2017Zaldarriaga et al. 2018). A transition in the evolution of a 2 is observed in the left panel in Figure 3 when material that is effectively torqued by the binary is able to form a disk and is subsequently accreted onto M 2 . This material has a higher specific angular momentum than the one initially set by tidal synchronization and, when accreted, is able to spin up the newly formed BH at a faster rate. The resultant change in slope observed in the left panel in Figure 3 due to the accretion of disk material is observed to occur earlier for smaller separations, which results in higher total spins values than those given by direct collapse of the same fallback material.
At a fixed ζ, the spin of the newly formed BH depends sensitively on d * . For d * 5, the resultant torque on the fallback material can be considerable and, as a result, a 2 can be appreciable larger than the one expected from tidal synchronization. In this case, the final mass of M 2 remains unchanged while the final spin can vary drastically. The final mass of M 2 is, on the other hand, controlled by ζ. The middle panel in Figure 3 shows the evolution of a 2 for a fixed separation d * = 3 and changing ζ. Initially, the spin evolution follows the trend expected from direct collapse. This is because the torque is unable to modify the original angular momentum of the promptly collapsing stellar material. A transition to disk accretion is seen in all cases, with the shift always occurring late in the mass accretion history of M 2 . The resultant spin is similar in all cases due to the self similarity of the mass distribution of the expanding ejecta (Matzner & McKee 1999), which results in a comparable mass ratio of directly falling stellar material to disk material for different values of ζ. For example this fraction varies from 5.3 % for ζ = 0.5 to 7.6% for ζ = 0.9 (see middle panel of Figure 3 for d * = 3). This mass ratio is mainly responsible for determining the final spin of the BH and varies only slightly with ζ.
We have discussed, in the context of the classical scenario, the effects that the binary separation and the energy of the SN have on the resulting spin of the newly formed BH. The right panel of Figure 3 provides a clear summary of our findings as it shows the final spin of M 2 as function of d and ζ. The final mass of the newly formed BH is also shown by the size of the symbols. The masses for M 2 range from 22.5M for ζ = 0.1 (M 2 /M * ≈ 0.8) to 8.2M for ζ = 0.9 (M 2 /M * ≈ 0.3). Together with the results from our simulations we also plot the expected spin obtained from the direct collapse of the pre-SN stellar profile. This formalism makes use of the KEPLER model and assumes solid body rotation determined by tidal synchronization. Then, by assuming the spherical collapse of the star, we obtain the BH's spin a 2 for different fractions f = 0.3, 0.5, 0.8, 1 of the collapsed stellar mass M * . If the entire star was to collapse directly onto a BH, this will give a final spin a 2 solely dependent on d, as predicted by the dashed line in Figure 3 labeled M * .
When ζ is small and a significant fraction of the material is promptly accreted by the BH, the simple direct collapse formalism provides an accurate description of the final spin of the newly formed BH. This can be seen by comparing the dashed line in Figure 3 labeled 0.8M * with the simulation results obtained for ζ = 0.1, which give BHs with M 2 ≈ 0.8M * and final spins that closely resemble the direct collapse ones. By contrast, when f M * M * , the final spin is significantly higher than the one predicted by direct collapse of the same enclosed material. This is because in such cases the fallback material is effectively torqued by the BH companion, which results in disk formation and consequentially higher final spin values. Binary BH formation in the classical scenario depends critically on the currently poorly constrained energy of the resulting SN, which for fallback-mediated remnant growth results in faster spinning BHs than what would have been attainable for a single star progenitor. Figure 4 shows the dependence of a 1 and M 1 on ζ and d. In contrast to M 2 , the final mass of the companion BH is only weakly altered by changes in ζ. The reason is that a comparatively small mass can be effectively restricted to the region within which the expanding envelope material is gravitationally bound to M 1 . This bound material forms a disk whose final mass depends on both d and ζ. The resultant changes in a 1 , under the assumption of a 1 (t = 0) = 0, are observed to be more pronounced when the initial binary separation changes. Although, as expected, no sizable changes take place at large separations given that only a tiny fraction of the companion's envelope can be under the gravitational influence of M 1 . The final value of a 1 shows a modest variation with SN energy with a small preference for ζ ≈ 0.5 at small separations. This indicates that although the ejecta distributions are similar for changing values of ζ, the fraction of bound material to M 1 is largest for this particular explosion energy, although its exact value is likely to change for different pre-SN progenitors. collapse. Our key findings are summarized below.

Spin Evolution of the Orbiting Black Hole
• The final mass of the newly formed BH depends on the explosion energy. Its mass varies from M 2 ≈ 0.8M * for ζ = 0.1 to M 2 ≈ 0.3M * for ζ = 0.9 ( Figure 1).
• At a fixed SN energy, the final spin increases significantly with decreasing d as a larger fraction of the fallback material is torqued by the companion. This results in similar mass BHs but with widely different spins (see left panel in Figure 3).
• Due to the self similarity of the mass distribution of the expanding ejecta, the final spin of the BH varies only slightly with ζ. This results in BHs with a wide range in masses but similar spins (see middle panel in Figure 3).
• In the presence of a companion, the final spin of a BH formed by a fallback SN explosion can be significantly higher than the one predicted by direct collapse of the same stellar material (see right panel in Figure 3).
• The spin of the BH companion, on the other hand, depends on both ζ and d. This is because its accretion history is determined by the amount of fallback material that it is able to seize (Figure 4).
In Figure 5 we present a comparison of our results (upper panel)  system. The shaded quadrilateral regions (upper and lower panel) show systems produced by the direct collapse of pre-SN helium stars of varying masses, whose structures have been taken from the KEPLER models of Woosley & Heger (2006). The final spin of M 2 is calculated using the radial stellar profile and assuming rigid body rotation of the tidally synchronized SN progenitor. The pre-SN helium stars (M * ) are assumed to be orbiting around a BH with M 1 = qM * and a 1 =0. The dependence of χ eff with M is obtained by varying q from 0.53 to 1 in all cases, while the dependence of χ eff at a fixed M is obtained by changing d from 2R * to 5R * at constant q. To facilitate comparisons, we plot as shaded ellipses the 90% credibility intervals of the GW signals measured so far (Abbott et al. 2016b(Abbott et al. ,c, 2017a. Some points should be emphasized. The current LIGO observations are inconsistent with the direct collapse of pre-SN helium stars in close binaries (Kushnir et al. 2016;Belczynski et al. 2017;Zaldarriaga et al. 2018;Hotokezaka & Piran 2017). When the assumption of direct collapse is relaxed, the mass of M 2 can be altered by small changes on the explosion energy ζ while a 2 and a 1 (to a lesser extent) depend primarily on d (Figure 4). For the specific 28M + 15M system studied here, we show that changes in ζ alone can produce systems like LVT151012 (ζ = 0.1, d * = 5) or GW170608 (ζ = 0.9, d * = 5). For a fixed SN energy of ζ = 0.9, changes in the initial separation can, on the other hand, yield systems like GW170608 (d * = 5) or GW151226 (d * = 2).
Irrespective of the exact progenitor system, the processes discussed here implies that the formation pathways of LIGO binary BHs are more complicated than the standard scenario suggests. But the effects are especially interesting for weak SN explosions taking place in close binary systems. Future LIGO observations can offer clues to the nature of the SN explosion leading to the formation of BHs, which is currently not well understood (Perna et al. 2014;Sukhbold et al. 2016;Raithel et al. 2018). For instance, GW170608 could be indicative of weak SN explosion of a more massive pre-SN progenitor system while GW151226 might arise due to direct collapse of a lighter, yet more compact progenitor system.
The properties of LIGO sources in the (χ eff ,M) plane is diverse. One appealing aspect of the classical scenario is that the great variety of binary and explosion parameters can probably help explain this diversity. Given the need for a large helium core mass in progenitors, BH formation may be favored not only by slow rotation but also by low metallicity (Izzard et al. 2004). Larger mass helium cores might have less energetic explosions but this is currently highly uncertain. Many massive stars may produce supernovae by forming neutron stars in spherically symmetric explosions, but some may fail during neutrino energy deposition, forming BHs in the centre of the star (Fryer et al. 2012;Ugliano et al. 2012;Pejcha & Thompson 2015) and possibly a wide range of weak SN explosions (Moriya et al. 2010;Fryer et al. 2012;Dexter & Kasen 2013;Lovegrove & Woosley 2013;Batta et al. 2017;Fernández et al. 2018). One expects various outcomes ranging from very massive BHs with low spins (GW150914), to lighter and faster spinning BHs (GW151226). The number density of binary BHs of different masses and spins would provide a natural test to distinguish between different stellar explosion avenues.