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

The gravitational wave (GW) signals detected by LIGO (Abbott et al., 2016b, c, 2017a, 2017b, 2017c) 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_{\ast}$) 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 $\chi_{\rm 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_{\ast}$ is imprinted on

Here $\vec{a_{1}}$ and $\vec{a_{2}}$ are the dimensionless spins of the BHs and $\hat{L}$ 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_{\rm\tau}$, 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 $\chi_{\rm 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.

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_{\odot}$ helium star in orbit around a BH of $M_{1}=15M_{\odot}$. 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\times 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_{\ast}R_{\ast}$, where $d_{\ast}=2$, $3$ or $5$ and $R_{\ast}$ 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_{\odot}$ pre-SN helium star with $R_{\ast}=0.76R_{\odot}$. We considered the innermost $3M_{\odot}$ 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: $\mathrm{E_{\rm SN}}=\zeta\ \mathrm{E_{\rm G}}$, where $E_{\rm G}=2.3\times 10^{52}{\rm 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: $\zeta=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_{\ast}$. 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 $\zeta\lesssim 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 $\zeta$ 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.

In this Letter we have explored within the classical binary scenario how the spin of LIGO sources vary when the remnant BH is formed in weak SN explosions instead of direct collapse. Our key findings are summarized below.

• •

  The final mass of the newly formed BH depends on the explosion energy. Its mass varies from $M_{2}\approx 0.8M_{\ast}$ for $\zeta=0.1$ to $M_{2}\approx 0.3M_{\ast}$ for $\zeta=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 $\zeta$. 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 $\zeta$ 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) in the context of both direct collapse solutions and current LIGO observations of binary BHs. Shown are $\chi_{\rm eff}$ as a function of the chirp mass, $\mathcal{M}$, of the resulting BH binary 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_{\ast}$) are assumed to be orbiting around a BH with $M_{1}=qM_{\ast}$ and $a_{1}$=0. The dependence of $\chi_{\rm eff}$ with $\mathcal{M}$ is obtained by varying $q$ from $0.53$ to $1$ in all cases, while the dependence of $\chi_{\rm eff}$ at a fixed $\mathcal{M}$ is obtained by changing $d$ from $2R_{\ast}$ to $5R_{\ast}$ 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, c, 2017a, 2017b, 2017c).

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 $\zeta$ while $a_{2}$ and $a_{1}$ (to a lesser extent) depend primarily on $d$ (Figure 4). For the specific $28M_{\odot}$ + $15M_{\odot}$ system studied here, we show that changes in $\zeta$ alone can produce systems like LVT151012 ($\zeta=0.1$, $d_{\ast}=5$) or GW170608 ($\zeta=0.9$, $d_{\ast}=5$). For a fixed SN energy of $\zeta=0.9$, changes in the initial separation can, on the other hand, yield systems like GW170608 ($d_{\ast}=5$) or GW151226 ($d_{\ast}=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 ($\chi_{\rm eff}$,$\mathcal{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.