When Did the Remnant of GW170817 Collapse to a Black Hole?

For more than 25 yr, binary neutron-star (BNS) mergers have been hypothesized as the progenitors of short-hard gamma-ray bursts (GRBs; Eichler et al. 1989; Narayan et al. 1992; also see Baiotti & Rezzolla 2017; Paschalidis 2017; Nathanail 2018b for reviews). The neutron-rich dynamical ejecta in such mergers was also predicted to be the ideal site for r-process nucleosynthesis (Lattimer & Schramm 1976), the radioactive decay of which would power isotropic and bright kilonova emission (e.g., Li & Paczyński 1998; also see Fernández & Metzger 2016 for reviews). Both of these predictions and the detection of gravitational waves (GWs) from a BNS merger were confirmed in the detection of GWs from GW170817 (Abbott et al. 2017b) and the almost coincident prompt gamma-ray emission from the short GRB 170817A (Abbott et al. 2017a). In less than 11 hr, the optical counterpart was localized (Coulter et al. 2017) in the nearby (≃40 Mpc) elliptical galaxy NGC 4993, and a large-scale observational campaign was launched to obtain the most detailed observations in the UV–optical–NIR energy bands of the kilonova emission (Abbott et al. 2017c).

The total mass of the merged object, determined by restricting the spins of the two coalescing neutron stars to be between that inferred from observed BNSs, is M_tot = ${2.74}_{-0.01}^{+0.04}$ M_⊙ (Abbott et al. 2017b). This is significantly larger than the highest values of masses measured for neutron stars, which are 1.97 ± 0.04 M_⊙ (Demorest et al. 2010) and 2.01 ± 0.04 M_⊙ (Antoniadis et al. 2013). Furthermore, M_tot well exceeds the maximum neutron-star mass predicted by most of the popular equations of state (EOSs). Therefore, it is expected that the merger remnant, after having ejected a (small) fraction of the total rest mass, will ultimately collapse to a black hole (BH).

Whether this occurs promptly or not is an important question. Prompt collapse to BH formation in a BNS merger can be halted for a short time due to the formation of a hypermassive neutron star (HMNS) that is supported against gravitational collapse by differential rotation, which allows for equilibrium solutions with masses that are substantially larger (Baumgarte et al. 2000; Weih et al. 2018) than those allowed by uniform rotation (Breu & Rezzolla 2016). In fact, using full general relativistic simulations of both equal- (e.g., Shibata & Taniguchi 2006; Baiotti et al. 2008) and unequal-mass binaries (e.g., Rezzolla et al. 2010; Hotokezaka et al. 2011) and for a number of EOSs, it was shown that if ${M}_{\mathrm{tot}}\lesssim 2.7\,{M}_{\odot }$, then the remnant could be a long-lived HMNS with a collapse time t_coll ≳ 10 ms.

The fact that the remnant does not collapse promptly to a BH suggests that in such BNS mergers, the arrival time of the GW and the prompt gamma-ray emission will not be coincident, and the GRB onset will be delayed by at least the collapse time, such that ${t}_{\mathrm{del}}\gt {t}_{\mathrm{coll}}$. Of course, in this line of reasoning, we are assuming that the relativistic jet is indeed launched by the BH; this is our working assumption hereafter and is supported by numerous numerical-relativity simulations that have shown that magnetically confined jet structures may form once a torus is present around the BH (Rezzolla et al. 2011; Dionysopoulou et al. 2015; Paschalidis et al. 2015; Kawamura et al. 2016; Ruiz et al. 2016).

With the exception of the case of a prompt collapse, determining the collapse time of the remnant for a binary of a given total mass is extremely challenging, as it depends on the unknown EOS, the rotational profile of the HMNS, and the numerous dissipation processes (radiative and not) that operate to distribute angular momentum and bring the star to uniform rotation. Moreover, self-consistent numerical-relativity simulations—either in pure hydrodynamics or in magnetohydrodynamics (MHD)—are not particularly accurate after the merger because of the appearance of strong shocks, which inevitably reduce the convergence order to first (Rezzolla & Zanotti 2013). As a result, most of the three-dimensional numerical-relativity simulations are limited to post-merger timescales of t ∼ 10–100 ms (Baiotti & Rezzolla 2017; Paschalidis 2017), and longer integration times become computationally unfeasible; such a timescale is also the one that has been explored in practice to assess the stability of the HMNS. If the actual collapse time from numerical simulations is rather uncertain (see Radice et al. 2018a for a systematic discussion), the upper limit for the lifetime of the remnant is on somewhat firmer ground and has been estimated to be ∼10⁴ s (Ravi & Lasky 2014). While this upper limit is quite generous, it has been deduced when considering a variety of possible scenarios for the evolution of the remnant.

Overall, numerical-relativity simulations agree in finding that after the merger, a double-core structure forms first, which then quickly relaxes to the newly formed HMNS that has a highly nonaxisymmetric bar-like structure that rotates very rapidly. This generates a substantial time-varying quadrupole moment and hence a strong emission of GWs over the next several tens of milliseconds that, alone, could cause the remnant to collapse. The law of differential rotation has very specific properties, with the angular velocity profile being characterized by a slowly rotating core and an envelope that rotates at quasi-Keplerian frequencies (Shibata & Taniguchi 2006; Kastaun & Galeazzi 2015); more importantly, this behavior has been found to be essentially universal (Hanauske et al. 2017b), i.e., only weakly dependent on the EOS.

Differential rotation in the HMNS is, however, expected to be damped on longer timescales, i.e., ≳100 ms, either by direct coupling with seed magnetic fields or via other dissipative effects. In the first case, magnetic fields can be amplified dynamically after merger because of a number of instabilities that are expected to develop, with the Kelvin–Helmholtz (Price & Rosswog 2006) and magnetorotational (MRI; Siegel et al. 2013) instabilities being those that give the fastest exponential growth. Despite some computationally very expensive attempts having been made to measure with confidence the magnetic-field strength in the HMNS (Kiuchi et al. 2015), the latter is still quite uncertain; overall, an amplification of at least 3 orders of magnitude is expected, thus bringing the initial magnetic fields from their canonical value of B ≲ 10¹² G to bulk values of B ≳ 10¹⁵ G. In the second case, viscous effects, either “direct,” i.e., in terms of pure shear and bulk viscosity of nuclear matter (Alford et al. 2018), or “effective,” i.e., in terms of the development of the MRI or neutrino transport, could significantly affect the distribution of angular momentum in the HMNS (Duez et al. 2004; Radice 2017; Shibata & Kiuchi 2017).

Here we use the properties of the kilonova emission and the delay time between the GW chirp signal and prompt gamma-ray emission onset to constrain the collapse time of the remnant. Earlier works either constrained the remnant collapse time to a broad range or provided only approximate estimates. For example, using a constant mean mass-ejection rate, Metzger et al. (2018) showed that a model involving a rapidly spinning HMNS with an ordered surface magnetic-field strength of ≈10¹⁴ G and a lifetime ${t}_{\mathrm{coll}}\sim 0.1\mbox{--}1\,{\rm{s}}$ can simultaneously explain the velocity, total mass, and electron fraction of the blue-ejecta mass. On the other hand, Granot et al. (2017) argued that a lifetime ${t}_{\mathrm{coll}}\sim 1\,{\rm{s}}$ is favored over prompt collapse, as it would naturally explain the delay time of t_del ≈ 1.74 s, given that the relativistic jet launched after collapse would then have to clear ∼1 s worth of neutrino-driven wind and dynamical ejecta.

In this work, we present a comprehensive account of the different mass-ejection channels and provide analytical mass-ejection rates derived from the results of several numerical simulations from the literature. In particular, based on the time-dependent rates, we argue that in order to produce the requisite mass of M_ej,blue ≈ 0.025 M_⊙ in the lanthanide-poor ejecta that gave rise to the rapidly fading bluer emission in the UV and optical at early times, the collapse time of the remnant cannot be much smaller than t_coll ≃ 1 s. More importantly, we arrive at a similar conclusion from an independent line of argument, where we model the dynamical evolution of the relativistic jet launched after the remnant collapses and bores its way out of the dynamical ejecta. In particular, we show that the observed delay time of t_del ≈ 1.74 s between the GW and gamma-ray emission can be satisfied if the remnant collapsed after about 1 s.

The rest of the paper is structured as follows. In Section 2 we discuss the different merger remnants obtained in BNS mergers and outline the basic properties that determine the merger outcome—prompt or delayed collapse. A detailed summary of the different mass-ejection channels, along with simple analytic expressions describing the mass-ejection rates and the total ejected masses as a function of time, are presented in Section 3. The collapse time of the remnant obtained from the condition of producing the requisite amount of blue-ejecta mass is calculated in Section 4, while in Section 5 we present the semi-analytic theory that describes the dynamical evolution of the jet-cocoon system and calculate the jet-breakout time. These results are then used, along with the mass-ejection rates, to constrain in Section 6 the collapse time of the remnant in GW170817. Finally, we discuss the implications of our results in Section 7 and provide our conclusions in Section 8.

In what follows, we review and summarize the two different lines of argument, stemming either from considerations of the maximum mass of self-gravitating equilibrium configurations (Section 2.1) or from the observations of the kilonova emission (Section 2.2), that support the idea that the remnant in GW170817 collapsed to a BH sometime after the merger (see, e.g., Margalit & Metzger 2017; Shibata et al. 2017; Rezzolla et al. 2018; Ruiz et al. 2018).

We start by recalling that there are four possible outcomes for the remnant of a BNS merger (e.g., Baiotti & Rezzolla 2017): (i) prompt collapse to a BH, (ii) an HMNS that is supported by differential rotation and thermal pressure, (iii) a supramassive neutron star (SMNS) that is supported by uniform rotation, and (iv) a stable neutron star.³ Each outcome depends on the total mass of the system and the underlying EOS. While the former can be obtained with higher certainty by measuring the masses of each of the neutron stars from the GW detection, the relatively larger uncertainty in the latter allows various possibilities for the remnant. However, here we assume that for events that are associated with short GRBs, the remnant must collapse to a BH. This condition is necessitated by the fact that a BH is most likely needed (e.g., Rezzolla et al. 2011) to power the relativistic outflow, which then produces a short GRB due to internal dissipation in the jet. This would naturally invalidate the existence of any kind of long-lived stable neutron star. Nevertheless, the possibility that an SMNS with magnetar-like surface dipole magnetic fields (i.e., B_s ≳ 10¹⁴ G) powers the GRB (Zhang & Mészáros 2001; Gao & Fan 2006; Metzger et al. 2008) before its ultimate collapse to a BH cannot be excluded. Indeed, the existence of such a protomagnetar could help explain some puzzling long-term and sustained X-ray emission associated with the afterglows of a number of short GRBs (Ciolfi & Siegel 2015; Rezzolla & Kumar 2015), for which the delay time is much larger (i.e., $\sim {10}^{2}\,{\rm{s}}$).

2.1. Evidence for Delayed Collapse: Maximum Mass

Determining the threshold mass to prompt collapse is possible via numerical simulations, and this has been exploited to set lower-bound constraints on the radii of neutron-star models (Bauswein et al. 2013, 2017; Köppel et al. 2019). In its lowest-order approximation, the threshold mass, that is the minimum mass above which a self-gravitating system (e.g., the merger remnant) will undergo prompt collapse, can be expressed in terms of the maximum mass for a nonrotating configuration, M_TOV, as M_th/M_TOV ≈ 1.415, where this expression has an uncertainty of ΔM_th = 0.05 M_⊙ (Köppel et al. 2019). Using this expression and taking a conservative value of M_TOV ≃ 2 M_⊙, we readily derive that M_th ≃ 2.82 M_⊙, and hence above the value measured for M_tot by Abbott et al. (2017b).⁴ Stated differently, even taking a conservative point of view, it is unlikely that the remnant of GW170817 collapsed to a BH promptly.

At the same time, it is also reasonable to expect that the remnant did collapse to a BH at some point after the merger, simply based on the estimates about the maximum mass that can be supported by a self-gravitating fluid configuration. We recall that Breu & Rezzolla (2016) showed that a quasi-universal relation (i.e., a relation that is essentially independent of the EOS) exists between the maximum mass that a neutron star can support via uniform rotation, M_max, and the corresponding mass for a nonrotating configuration. Exploring a very large number of EOSs, the maximum mass was found to be

$$\begin{eqnarray}&&{M}_{\max }=\left(1.20\pm 0.02\right){M}_{\mathrm{TOV}}\simeq 2.40\,{M}_{\odot }\lt {M}_{\mathrm{tot}},\end{eqnarray} \tag{ 1 }$$

again considering that M_TOV ≃ 2 M_⊙. Similar considerations can also be made for configurations that are differentially rotating, once a precise form for the law of differential rotation has been chosen. In particular, when considering the most common law of differential rotation, i.e., the j-constant law, Weih et al. (2018) found that a quasi-universal relation between the maximum mass of a differentially rotating neutron star and the maximum nonrotating mass, i.e.,

$$\begin{eqnarray}&&{M}_{\max ,\mathrm{dr}}\simeq \left(1.54\pm 0.05\right){M}_{\mathrm{TOV}}\simeq 3.08\,{M}_{\odot }\gt {M}_{\mathrm{tot}},\end{eqnarray} \tag{ 2 }$$

once again after taking M_TOV ≃ 2 M_⊙. In summary, these two results combined suggest that the remnant of GW170817 is likely to have eventually collapsed to a BH, but also that this must have happened when the remnant had lost a certain amount of differential rotation.

The timescale over which differential rotation is lost depends on the intervening process. In the presence of magnetic fields, differential rotation in an HMNS is damped by magnetic braking on the Alfvén timescale (Shapiro 2000)

$$\begin{eqnarray}\begin{array}{rcl}{t}_{{\rm{A}}} & = & \displaystyle \frac{{R}_{\mathrm{NS}}}{{v}_{{\rm{A}}}}\\ & \approx & \ 120{\left(\displaystyle \frac{{B}_{\mathrm{int}}}{{10}^{15}{\rm{G}}}\right)}^{-1}{\left(\displaystyle \frac{{R}_{\mathrm{NS}}}{12\mathrm{km}}\right)}^{-1/2}{\left(\displaystyle \frac{M}{2.7{M}_{\odot }}\right)}^{1/2}\,\mathrm{ms},\end{array}\end{eqnarray} \tag{ 3 }$$

where ${v}_{{\rm{A}}}:= {B}_{\mathrm{int}}/\sqrt{4\pi \rho }$ is the Alfvén speed in the neutron-star interior, ρ is its mean rest-mass density, and ${B}_{\mathrm{int}}$ is the internal magnetic-field strength. However, differential rotation in the inner regions of the HMNS can be removed on a much shorter timescale and of the order of a few tens of dynamical timescales, i.e., 10–15 ms (Kastaun & Galeazzi 2015; Hanauske et al. 2017b). Once differential rotation is lost, and assuming the remnant has not yet collapsed to a BH, the HMNS will effectively become an SMNS, and angular momentum will be lost over a much longer timescale, since the GW is essentially quenched. This spin-down is likely to be due to magnetic-dipole radiation, such that the spin-down timescale will be

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{sd}} & = & {f}^{-1}\,\displaystyle \frac{{{Ic}}^{3}}{2{{\rm{\Omega }}}_{0}^{2}{R}_{{}_{\mathrm{NS}}}^{6}{B}_{s}^{2}}\\ & \approx & \ \displaystyle \frac{1.1\,\times \,{10}^{4}}{f}{\left(\displaystyle \frac{{B}_{s}}{{10}^{14}{\rm{G}}}\right)}^{-2}{\left(\displaystyle \frac{{P}_{0}}{1\mathrm{ms}}\right)}^{2}{\left(\displaystyle \frac{{R}_{\mathrm{NS}}}{12\mathrm{km}}\right)}^{-6}\,{\rm{s}},\end{array}\end{eqnarray} \tag{ 4 }$$

where I ≃ 10⁴⁵ g cm² is the moment of inertia,⁵ B_s is the initial surface dipole magnetic-field strength, Ω₀ $:=$ 2π/P₀ is the initial angular spin frequency, and P₀ is the initial spin period. Setting θ_B as the inclination angle between the magnetic-dipole axis and the rotational axis, the parameter f depends on the details of the model and is given by f = χ sin² θ_B in vacuum, with χ being a constant (χ = 2/3) in Newtonian gravity, while a function of the stellar compactness in general relativity, e.g., χ ∼ 2–3 (Rezzolla et al. 2001, 2003); on the other hand, f = 1 + sin² θ_B in force-free models (Spitkovsky 2006).

Since most of the rotational energy of the star is released at the spin-down time, this timescale has been linked to the plateau phases observed in the X-ray light curves of many short-hard GRBs (e.g., Zhang & Mészáros 2001; Metzger et al. 2011; Gompertz et al. 2013; Rowlinson et al. 2013). In this scenario, the rotational spin-down of the SMNS would power the late-time energy injection into the relativistic outflow, followed by the collapse of the SMNS to a BH that would manifest as a sudden drop in the X-ray flux (however, see, e.g., Ciolfi & Siegel 2015; Rezzolla & Kumar 2015 for an alternative explanation).

In the case of GW170817, if the remnant was an SMNS, then the large rotational energy released in the form of an isotropic MHD wind would have produced a spin-down luminosity of ${L}_{\mathrm{sd}}\gg {10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ for a surface dipole magnetic-field strength of B_s ≤ 10¹⁵ G. However, the observed bolometric luminosity at t ≲ 10 days post-merger was ${L}_{\mathrm{bol}}\,\lt {10}^{42}\,\mathrm{erg}\ {{\rm{s}}}^{-1}$ and cannot be explained even if the surface dipole field was as high as B_s ≤ 10¹⁶ G (e.g., Margutti et al. 2018). In addition, such a high rate of energy injection by the spinning down SMNS would have powered an exceptionally bright afterglow emission that was not observed. Both of these reasons ruled out the possibility of an SMNS remnant (e.g., Granot et al. 2017; Margalit & Metzger 2017; Pooley et al. 2018).

2.2. Evidence for Delayed Collapse: Kilonova Emission

We next consider a second and distinct line of argument, based this time on the observation of the prompt emission and properties of the outflow that powered the kilonova emission, which again suggest that the collapse to a BH for the remnant in GW170817 was delayed.

We recall that the detection of GW170817 was also rich in electromagnetic (EM) counterparts. The UV/optical/IR counterparts that followed the GW were interpreted as the result of heating due to the decay of freshly synthesized heavy elements. In addition, the onset of the EM emission, that is, of the prompt gamma-ray emission, as recorded by Fermi/GBM and INTEGRAL/SPI-ACS, was delayed by t_del = 1.74 ± 0.05 s with respect to the GW chirp signal that marks the merger of the two neutron stars (Abbott et al. 2017a). The spectrum of the rapidly fading early UV emission was well described by a blackbody, and the average velocity of the photosphere was determined to be β $:=$ v/c = 0.3–0.2 over the course of ∼0.6–1 day after the merger (Arcavi et al. 2017; Drout et al. 2017; Evans et al. 2017; Kasen et al. 2017; Kasliwal et al. 2017; Kilpatrick et al. 2017; Pian et al. 2017; Shappee et al. 2017; Smartt et al. 2017). Such high velocities and the trend toward redder optical and NIR emission that was seen after ∼2 days agree well with predictions from kilonova modeling (e.g., Metzger et al. 2010).

However, such observations also require two separate components in the outflow, one to explain the early UV or “blue” kilonova and another for the “red” kilonova that showed a gradual decline in flux over the initial ∼2–3 weeks. The difference between the two components is that the outflow that powered the blue kilonova was comprised of lanthanide-poor material and therefore had a higher electron fraction, with Y_e,blue ∼ 0.2–0.3, and lower opacity, κ_blue ≲ 1 cm² g⁻¹, to bound–bound electronic transitions. On the other hand, the outflow that powered the red kilonova was relatively lanthanide-rich, with a much lower Y_e,red ≲ 0.2 and much higher opacity with ${\kappa }_{\mathrm{red}}\sim 10\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$.

The high-opacity lanthanide-rich ejecta is expected to emerge from the dynamical ejection of material, predominantly due to tidal stripping (e.g., Rezzolla et al. 2010), as the two neutron stars merge. The ejected matter in this case lies mostly in the equatorial plane and is confined to lower latitudes. On the other hand, the lower-opacity lanthanide-poor ejecta can originate from neutrino-irradiated winds from the accretion disk (e.g., Fernández et al. 2015a; Just et al. 2015) and/or from an HMNS (e.g., Perego et al. 2014). This is illustrated in the cartoons in Figure 1, where the left panel reports a scheme of the matter distribution before the collapse of the remnant. In this case, mass is ejected dynamically due to shock heating and neutrino- and MHD-driven winds from the HMNS and the disk. This results in a two-component ejecta: a lanthanide-poor with high electron fraction (0.2 ≲ Y_e ≲ 0.3) ejecta that powers the early rapidly fading bluer emission, and a lanthanide-rich with low Y_e ≲ 0.2 ejecta that powered the longer-lasting redder emission. The increase in Y_e at higher latitudes is caused by two effects: (i) ejection of only moderately neutron-rich dynamical ejecta in the polar regions during merger and (ii) neutrino irradiation of the ejecta due to neutrinos from the HMNS. After collapse to a BH (right panel of Figure 1), an accretion torus is formed, and the BH is expected to launch a relativistic jet. Mass ejection at this point primarily occurs from the accretion torus due to MHD and viscously driven winds. As the jet traverses the blue ejecta, it slows down initially and becomes collimated by the cocoon formed by the shock-heated jet material and surrounding ejecta. As the jet clears the homologously expanding ejecta, it breaks out and produces prompt gamma-ray emission as a result of internal dissipation.

Zoom In Zoom Out Reset image size

Had the remnant of GW170817 suffered a prompt collapse, the newly formed BH would launch a relativistic jet almost immediately, apart from a negligible delay on the order of the dynamical time in assembling an accretion torus. In this scenario, the motion of the jet would be essentially unimpeded, since most of the ejected mass would be in the equatorial plane and at high colatitudes with θ > 30° (measured from the jet symmetry axis; for an example, see Figure 10 of Bovard et al. 2017, bottom row, third panel from the left). Furthermore, the jet would cross it in negligible time if any small amount of circum-merger ejecta were to be present in its path. Therefore, in the case of prompt collapse, the only time delay would be the radial time delay t_R $:=$ R_γ/2Γ²c for emission coming from material along our line of sight due to the less-than-speed-of-light motion of the relativistic flow in reaching the emission radius R_γ. However, as discussed in the previous section, and as we will show in Section 5 with more detailed calculations, t_R < t_del, thus making the prompt-collapse scenario problematic for GRB 170817A.

A prompt collapse also creates two additional tensions with the observations. First, the disk-wind component alone is not able to eject enough material in order to satisfy the total ejected mass needed for the blue component. Second, it cannot yield the high electron fractions (Y_e ∼ 0.3) at high latitudes that are needed to explain the blue kilonova. Both of these issues can be settled if the central remnant is an HMNS. First, hotter surface temperatures lead to a significantly higher level of neutrino irradiation, which ultimately raises Y_e to larger values that can be as high as Y_e ≲ 0.45. Second, it was argued by Metzger et al. (2018) that an enhancement in the mass-ejection rate and the inferred velocity of β_max,blue ∼ 0.2–0.3 can be obtained via a strongly magnetized, neutrino-irradiated wind of an HMNS that possesses magnetar-strength surface dipole fields. In this case, mass ejection is aided by magnetocentrifugal acceleration along open magnetic-field lines, whereas mass ejection due to the thermal pressure by neutrino heating alone would yield an insufficient ${M}_{\mathrm{ej},\mathrm{blue}}$ and a smaller asymptotic ejecta velocity of ${\beta }_{\max ,\mathrm{blue}}\sim 0.1$. On the other hand, if the collapse was delayed due to the formation of an HMNS, which survived for ≲1 s, then including the breakout time of the relativistic jet from the quasi-spherically expanding ejecta into the total delay time would naturally explain the delayed onset of the prompt gamma-ray emission (e.g., Granot et al. 2017).

In the following three sections, we will discuss how to set constraints on the actual time of collapse of GW170817's remnant by modeling either the ejected mass or EM delay.

As mentioned above, the UV/optical/IR observations following GW170817 are normally explained as the result of heating due to the decay of freshly synthesized heavy elements contained in the mass ejected during and after the merger of the binary. The ejected mass was found to be ${M}_{\mathrm{ej}}\approx 0.02\mbox{--}0.06\,{M}_{\odot }$, which was confirmed by several different groups (Chornock et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Nicholl et al. 2017; Perego et al. 2017; Tanaka et al. 2017; Tanvir et al. 2017; Villar et al. 2017; Waxman et al. 2017; Kawaguchi et al. 2018; Metzger et al. 2018). For convenience, we follow Kasen et al. (2017), whose kilonova analysis revealed a two-component ejecta—a lanthanide-poor ejecta with mass ${M}_{\mathrm{ej},\mathrm{blue}}\approx 0.025\,{M}_{\odot }$ and a lanthanide-rich ejecta with mass M_ej,red ≈ 0.04 M_⊙—and use these values to constrain the lifetime of the remnant and indicate the collapse time.

We recall that BNS mergers are expected to eject mass in a variety of different processes. In the last few years, thanks to numerical simulations and semi-analytic modeling, a robust picture has been drawn of the different components of ejected mass, which can be divided into three main ejection mechanisms:

• 1.  
  matter dynamically ejected (Rosswog et al. 1999; Shibata et al. 2005; Baiotti et al. 2008; Kiuchi et al. 2010; Hotokezaka et al. 2011, 2013; Rezzolla et al. 2010, 2011; Palenzuela 2013; Palenzuela et al. 2015; Sekiguchi et al. 2015; Radice et al. 2016; Lehner et al. 2016; Sekiguchi et al. 2016; Foucart et al. 2016; Bovard et al. 2017; Dietrich et al. 2017b, 2017a; Dietrich & Ujevic 2017; Radice et al. 2018b; Papenfort et al. 2018),
• 2.  
  matter ejected via neutrino-driven winds (Metzger et al. 2008, 2009; Dessart et al. 2009; Lee et al. 2010; Fernández & Metzger 2013; Perego et al. 2014; Just et al. 2015; Martin et al. 2015), and
• 3.  
  matter ejected via magnetically driven winds (Siegel et al. 2014; Fernández et al. 2015a, 2018; Ciolfi et al. 2017; Fujibayashi et al. 2017, 2018; Siegel & Metzger 2017, 2018).

We note that some works also discuss “viscosity-driven” ejected matter (Metzger et al. 2008; Goriely et al. 2011; Fernández & Metzger 2013; Fujibayashi et al. 2018; Fahlman & Fernández 2018). These viscous effects, however, are ultimately originated by magnetic turbulence driven by magnetic instabilities in either the HMNS disk or the torus; because of their magnetic origin, we will hereafter associate them with category 3.

All of these channels of ejection of matter have been explored in a large number of works and by several groups employing rather different numerical techniques and approximations (e.g., Fernández et al. 2015a; Just et al. 2015). In what follows, we discuss the mass-ejection rates for each channel separately and provide simple analytic prescriptions that provide simple and hopefully robust representations of the complex results of the numerical simulations.

3.1. Dynamical Mass Ejection

Dynamical mass ejection is a robust feature of BNS merger simulations and has been reported by a number of groups under a variety of physical conditions. The exact amount of mass ejected, however, depends both on the properties of the binary (mass, mass ratio, EOS, initial spin, magnetization, etc.) and on the criterion used to perform the measurement, which could even lead to a 100% difference (see, e.g., Bovard & Rezzolla 2017 for a discussion). Furthermore, it is presently difficult to predict with even decent precision the amount of matter to be ejected without a fully dynamical simulation, as the data do not seem to be well captured by simple analytic expressions. To counter this problem, we have exploited the data reported by Radice et al. (2018b), who presented measurements for more than 30 simulations. For the sake of the arguments made here, we simply need an average of the amount of matter ejected dynamically, which we deduce from the values for the binaries in the mass range $2.7{M}_{\odot }\,\le M\le 2.8{M}_{\odot }$, reported in Table 2 of Radice et al. (2018b). In this way, we find a mean value with a considerable variance, i.e., ${M}_{{\rm{dyn}}}\approx (1.5\pm 1.1)\times {10}^{-3}{M}_{\odot }$; although the uncertainty associated with this estimate is obviously very large, it also represents the most reasonable one at present and reflects many of the uncertainties discussed in Bovard & Rezzolla (2017) and Bovard et al. (2017).

The data reported by Radice et al. (2018b) also allow us to reconstruct a mass-ejection rate, which we approximate as given by exponential growth followed by a power-law decay fitted to the data from the numerical simulations

$$\begin{eqnarray}\displaystyle \frac{{\dot{M}}_{\mathrm{dyn}}}{{M}_{\odot }\,{{\rm{s}}}^{-1}}=\left\{\begin{array}{cc}0.054\,\exp \left(t/{t}_{\mathrm{dyn}}\right), & t\lt {t}_{\mathrm{dyn}},\\ 1.4\times {10}^{-31}{\left(t/1{\rm{s}}\right)}^{-15}, & t\geqslant {t}_{\mathrm{dyn}},\end{array}\right.\end{eqnarray} \tag{ 5 }$$

where t_dyn ≈ 10 ms is the typical timescale over which the dynamical ejection has saturated after the merger. In this case, we assume a constant fractional uncertainty of 70%. Here, and for all of the following expressions for the mass-ejection rate, we set t = 0 to be the time of the merger. These dynamical-ejection contributions are shown in Figure 2, where the top panels report the evolution of the various channels in which mass is ejected (see also below), and the bottom panels show the corresponding rates. Shown instead on the left and right are the portions of the evolution on very short timescales after the merger (left panels) and timescales comparable with the collapse time (right panels).

Zoom In Zoom Out Reset image size

3.2. Neutrino-driven Winds

After merger, neutrino-driven winds are expected to be generated from the remnant and the surrounding disk. Determining the mass-ejection rate from the various simulations performed to model this process is made difficult by the fact that only rarely is the mass ejected via neutrinos reported. Fortunately, however, a close relation exists between the neutrino-driven mass-ejection rate ${\dot{M}}_{\nu }$ and the neutrino luminosity L_ν in its dominant components (i.e., electron neutrino and electron antineutrino), and it is given by (Qian & Woosley 1996; Rosswog & Ramirez-Ruiz 2002; Dessart et al. 2009)

$$\begin{eqnarray}&&{\dot{M}}_{\nu }\approx 5\,\times \,{10}^{-4}{\left(\displaystyle \frac{{L}_{\nu }}{{10}^{52}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{5/3}{M}_{\odot }\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 6 }$$

Since most publications report the measured neutrino luminosity, Equation (6) allows us to obtain ${\dot{M}}_{\nu }$ from L_ν. At the same time, since we want to distinguish the mass ejection before and after collapse, we will separately consider the emission coming from the HMNS and the surrounding disk (see left panel of the cartoon in Figure 1).

In particular, we first consider the neutrino-driven winds coming from the evolution of the HMNS, and to this scope, we use the simulations in Fujibayashi et al. (2017), where an HMNS with a disk is evolved under the assumption of axisymmetry but in a general-relativistic framework and for a timescale of 400 ms. In particular, we use the neutrino luminosity reported in Figure 3 of Fujibayashi et al. (2017) and fit it with two power-law segments that change at 200 ms after merger, as measured by the simulations. As a result, we deduce that the mass-ejection rate associated with the HMNS is

$$\begin{eqnarray}\displaystyle \frac{{\dot{M}}_{\nu ,\mathrm{HMNS}}}{{M}_{\odot }{{\rm{s}}}^{-1}}\approx \left\{\begin{array}{cc}1.13\times {10}^{-4}\,{\left(t/1{\rm{s}}\right)}^{-0.98}\ \ & {t}_{\mathrm{dyn}}\lt t\lt 0.2\,{\rm{s}},\\ 2.02\times {10}^{-4}\,{\left(t/1{\rm{s}}\right)}^{-0.62}\ \ & t\geqslant 0.2\,{\rm{s}}.\end{array}\right.\end{eqnarray} \tag{ 7 }$$

Three remarks should be made on Equation (7). First, since the data reported by Fujibayashi et al. (2017) naturally do not distinguish the mass-ejection rates from the HMNS and the disk, the values reported in Equation (7) effectively also account for the contribution from the disk, ${\dot{M}}_{\nu ,\mathrm{DISK}}$. As we will see below, however, the latter is at least 1 order of magnitude smaller and with a faster decay; hence, it can be neglected in Equation (7). Second, the rates in Equation (7) match the result of the simulations up to t = 0.4 s and are extrapolated for longer times; this is fine, given the steep decay after $t=0.2\,{\rm{s}}$ and the fact that this is not a dominant contribution to the ejected mass. Third, Equation (7) also comes with a fractional uncertainty, which we estimate to be of the order of 20% on the basis of the discussion of the convergence of the results in Section 3.5 of Fujibayashi et al. (2017). Finally, the contribution to the mass ejection from the HMNS disk will obviously cease once the remnant collapses to a BH. We model this sudden shutoff in the mass-ejection rate via exponential decay on a dynamical timescale ${t}_{\mathrm{dyn}}\simeq 10\,\mathrm{ms}$.

The contribution to the mass ejection coming from the disk around the HMNS is best approximated by exploiting the results of numerical simulations that include neutrino-driven winds (and subsequent nucleosynthesis) coming from a system composed of a BH and an accretion torus (Fernández et al. 2015b; Just et al. 2015, 2016; Fujibayashi et al. 2017; Fernández et al. 2018; Siegel & Metzger 2017, 2018). In particular, we have made use of the results from Siegel & Metzger (2017), which are in good agreement with those of Fernández et al. (2018), and fitted the data in the upper panel of Figure 13 of Siegel & Metzger (2017) with two power laws joined at t = 0.12 s, which then yields

$$\begin{eqnarray}\displaystyle \frac{{\dot{M}}_{\nu ,\mathrm{DISK}}}{{M}_{\odot }{{\rm{s}}}^{-1}}\approx \left\{\begin{array}{cc}1.43\,\times \,{10}^{-7}{\left(t/1{\rm{s}}\right)}^{-2.0}\, & {t}_{\mathrm{dyn}}\lt t\lt 0.12\,{\rm{s}},\\ 1.22\,\times \,{10}^{-11}{\left(t/1{\rm{s}}\right)}^{-6.4}\, & t\geqslant 0.12\,{\rm{s}},\end{array}\right.\end{eqnarray} \tag{ 8 }$$

with ${\dot{M}}_{\nu ,\mathrm{DISK}}$ ≈ ${\dot{M}}_{\nu ,\mathrm{tor}}$ and ${\dot{M}}_{\nu ,\mathrm{tor}}$ the mass-ejection rate from a BH-torus system. Note that ${\dot{M}}_{\nu ,\mathrm{DISK}}$/${\dot{M}}_{\nu ,\mathrm{HMNS}}$ ≪ 1, so that the assumption made in Equation (7) that ${\dot{M}}_{\nu ,\mathrm{HMNS}}+{\dot{M}}_{\nu ,\mathrm{DISK}}\,\approx {\dot{M}}_{\nu ,\mathrm{HMNS}}$ is indeed well justified.

To model the mass ejection at times beyond those simulated by Siegel & Metzger (2017), i.e., 400 ms, we simply extrapolate in time the second segment of the power law in Equation (8). Also in this case, we attempt to estimate an error in the mass-ejection rate of Equation (8) by characterizing the numerical noise reported in Figure 13 of Siegel & Metzger (2017), deducing a lower limit on the fractional error of Equation (8) of ∼20%. All of these different neutrino-driven contributions are shown in Figure 2.

3.3. MHD–Viscosity-driven Winds

Magnetically driven winds from the HMNS are expected to be of great importance for the overall ejected mass, especially for the blue-kilonova component (Metzger et al. 2018). The magnetic energy is only a small contribution of the total energy budget during the inspiral, and magnetic-induced effects are in fact expected to be negligible sometime before the merger (Giacomazzo et al. 2009). After the merger, however, a good fraction of the kinetic energy of the system can be converted into magnetic energy via magnetic shearing and instabilities. High-resolution simulations of the remnant have shown that instabilities such as the Kelvin–Helmholtz or the MRI can develop (Siegel et al. 2013), leading to magnetic energies as large as 10⁵⁰ erg (Obergaulinger et al. 2010; Giacomazzo et al. 2015; Kiuchi et al. 2015, 2018). Simulations have also shown that copious winds are a natural outcome of a highly shearing magnetized HMNS (Siegel et al. 2014).

Despite being potentially the most important, mass ejection from magnetically driven winds has been explored much less than that from other channels. It is therefore difficult to quantify the mass-ejection rate from the few works available in the literature (Siegel et al. 2014; Dionysopoulou et al. 2015; Ciolfi et al. 2017; Kiuchi et al. 2018). Here we make use of the results of Ciolfi et al. (2017), who investigated a long-lived magnetized HMNS and showed that the dynamically ejected mass does not depend on the magnetic-field strength (see the ejected mass for the magnetized and unmagnetized binaries in Figure 23 of Ciolfi et al. 2017). On the other hand, when the dynamical ejection of matter has saturated, the magnetized HMNS does eject mass at a larger rate than in the unmagnetized case, and we have found it to be well described by a simple power law of the type ${\dot{M}}_{{\rm{B}},\mathrm{HMNS}}\propto {t}^{-0.8}$. In similar general-relativistic MHD simulations that will be presented in a separate work (E. M. Most et al. 2019, in preparation), the mass-ejection rate has been found to have a similar but different power-law dependence, i.e., ${\dot{M}}_{{\rm{B}},\mathrm{HMNS}}$ ∝ t⁻¹. Hereafter, we choose an average value from the two simulations and set the mass-ejection rate due to MHD-driven winds from the HMNS (and its surrounding disk) to be

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{M}}_{{\rm{B}},\mathrm{HMNS}}}{{M}_{\odot }{{\rm{s}}}^{-1}}\approx 5.18\,\times \,{10}^{-3}{\left(t/1{\rm{s}}\right)}^{-0.9}\qquad {t}_{\mathrm{dyn}}\leqslant t\leqslant {t}_{\mathrm{coll}}.\end{eqnarray} \tag{ 9 }$$

As is evident from Figure 2, this mass-ejection channel dominates the mass loss at early times; thus, we associate a more conservative fractional error of 30% with this estimate. At the same time, because of the very small number of works exploring the effect of these winds from the HMNS, it is prudent to also consider a different scenario in which a 30% uncertainty is introduced for the exponent in the power law; we will discuss the two scenarios separately below.

The viscously/MHD-driven winds, first from the HMNS disk and later from the BH torus, are expected to be the dominant channel for mass ejection, especially after the HMNS collapses (Lippuner et al. 2017; Fahlman & Fernández 2018; Fujibayashi et al. 2018). We model this mass ejection using the long-term simulations of Fujibayashi et al. (2018) and deduce that the rate is impressively constant over a considerably long timescale and can be approximated as (see Figure 8 of Fujibayashi et al. 2018)

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{M}}_{{\rm{B}},\mathrm{DISK}}}{{M}_{\odot }\,{{\rm{s}}}^{-1}}\simeq 0.012\qquad \quad {t}_{\mathrm{dyn}}\leqslant t\leqslant {t}_{\mathrm{coll}},\end{eqnarray} \tag{ 10 }$$

Note that as long as t < 0.4 s, ${\dot{M}}_{{\rm{B}},\mathrm{DISK}}\ll {\dot{M}}_{{\rm{B}},\mathrm{HMNS}}$, so that its contribution in Equation (9) can be safely neglected (see the right panels of Figure 2 for a graphical impression). In addition, for 0.4 s < t < t_coll, the contribution of ${\dot{M}}_{{\rm{B}},\mathrm{DISK}}$ to the total ejected mass is very small (≲9 × 10⁻³ M_⊙). Furthermore, since the rates show a rather large variance, here too we associate a conservative 30% fractional error with the estimated rate (Equation (10)).

As for the case of neutrino-driven winds, the contribution of the torus to the mass-accretion rate due to magnetic winds will be slightly different from that of the disk around the HMNS, but certainly comparable. Hence, we assume that also in the case of magnetic-related rates, ${\dot{M}}_{{\rm{B}},\mathrm{DISK}}$ ≈ ${\dot{M}}_{{\rm{B}},\mathrm{tor}}$. The assumption that the two rates ${\dot{M}}_{{\rm{B}},\mathrm{DISK}}$ and ${\dot{M}}_{{\rm{B}},\mathrm{tor}}$ are comparable is indeed confirmed by the recent simulations of Fernández et al. (2018), who performed three-dimensional general-relativistic MHD simulations of a BH-torus system for very long times, i.e., timescales of the order of $\sim 9\,{\rm{s}}$. At the same time, the new simulations also reveal that the rate (Equation (10)) is not always constant and decreases rapidly ≈1 s after the start of the simulations (thus, approximately 1 s after the collapse of the HMNS to a BH) following a power law ∝t^−2.3 (Fernández et al. 2018). Hence, we model the magnetically driven mass-ejection rate from the disk/torus as

$$\begin{eqnarray}\displaystyle \frac{{\dot{M}}_{{\rm{B}},\mathrm{DISK}}}{{M}_{\odot }\,{{\rm{s}}}^{-1}}\approx \left\{\begin{array}{cc}0.012 & {t}_{\mathrm{coll}}\lt t\lt \bar{t},\\ 0.012{\left(t/\bar{t}\right)}^{-2.3} & t\geqslant \bar{t},\end{array}\right.\end{eqnarray} \tag{ 11 }$$

where $\bar{t}:= {t}_{\mathrm{coll}}+1\,{\rm{s}}$. Also, these magnetic contributions are shown in the various panels of Figure 2.

The lifetime of the remnant from GW170817 can now be conveniently constrained by comparing the amount of ejected matter discussed in the previous sections with the amount of matter that is deduced to have been ejected on the basis of the blue-kilonova emission. We recall that the latter is expected to be the result of lanthanide-poor material and thus with a lower opacity (Arcavi et al. 2017; Drout et al. 2017; Kasen et al. 2017; Kilpatrick et al. 2017; Villar et al. 2017). Producing such material requires high temperatures and densities, such that neutrino-irradiation processes can produce material with high electron fractions (i.e., neutron poor), which, in turn, are responsible for light r-process element production with intrinsically low opacities. Such large temperatures and densities can be found in the HMNS but not in the disk surrounding it or in the torus around the BH once it is formed. On the contrary, the disk outflows with low electron fraction (i.e., neutron-rich) material are expected to produce heavy r-process elements resulting in higher opacity (Kasen et al. 2013; Tanaka & Hotokezaka 2013; Metzger & Fernández 2014; Perego et al. 2014; Lippuner et al. 2017; Tanaka et al. 2018; Fujibayashi et al. 2018).

Exploiting the mass-ejection rates reported in the previous sections, we can estimate the collapse of the HMNS by requiring that the accumulated ejected mass from the HMNS becomes equal to the blue component from the kilonova modeling, which has been estimated to be ${M}_{\mathrm{ej},\mathrm{blue}}\approx 0.025\,{M}_{\odot }$ (Kasen et al. 2017). In particular, we assume that all of the mass lost by the HMNS contributes in its entirety to the mass in the blue component, M_ej,blue. On the other hand, the mass ejected dynamically, as well as that coming from the disk and torus, will contribute to both the blue and the red components. More specifically, we assume that one-third of the dynamically ejected mass contributes to the blue component, while two-thirds of it is channeled into the red component (Bovard et al. 2017; Perego et al. 2017; Shibata & Kiuchi 2017; Radice et al. 2018b). The same percentages are also applied for the matter ejected from the HMNS disk, and these values are in agreement with the results of Shibata & Kiuchi (2017) and Fujibayashi et al. (2018). On the other hand, once the HMNS collapses to a BH, the contributions change, as only a smaller fraction of the matter ejected is able to contribute to the blue component. Following Fernández et al. (2018), we distribute 1/12 of the torus-ejected matter to the blue component and the remaining 11/12 to the red component. Ideally, these splitting factors ought to be a function of time and have an associated uncertainty; however, for simplicity, we consider them as constant here and factor in the uncertainties in the ejection rates.

In summary, the mass-ejection rates (Equations (5)–(11)) and their contributions to the total blue- and red-ejecta masses

$$\begin{eqnarray}&&{M}_{\mathrm{ej},\mathrm{blue}/\mathrm{red}}(t)=\int {\dot{M}}_{\mathrm{ej},\mathrm{blue}/\mathrm{red}}(t^{\prime} ){dt}^{\prime} ,\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{\mathrm{ej},\mathrm{blue}}(t) & = & {\eta }_{\mathrm{dyn},\mathrm{blue}}\,{\dot{M}}_{\mathrm{dyn}}(t)\\ & & +\ {\eta }_{\mathrm{HMNS},\mathrm{blue}}[{\dot{M}}_{{\rm{B}},\mathrm{HMNS}}(t)+{\dot{M}}_{\nu ,\mathrm{HMNS}}(t)]\\ & & +\ {\eta }_{\mathrm{DISK},\mathrm{blue}}[{\dot{M}}_{{\rm{B}},\mathrm{DISK}}(t)+{\dot{M}}_{\nu ,\mathrm{DISK}}(t)],\end{array}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{\mathrm{ej},\mathrm{red}}(t) & = & {\eta }_{\mathrm{dyn},\mathrm{red}}\,{\dot{M}}_{\mathrm{dyn}}(t)\\ & & +\ {\eta }_{\mathrm{DISK},\mathrm{red}}[{\dot{M}}_{{\rm{B}},\mathrm{DISK}}(t)+{\dot{M}}_{\nu ,\mathrm{DISK}}(t)],\end{array}\end{eqnarray} \tag{ 14 }$$

where, for t ≤ t_coll,

$$\begin{eqnarray}\begin{array}{rcl}{\eta }_{\mathrm{HMNS},\mathrm{blue}} & = & 1,{\eta }_{\mathrm{dyn},\mathrm{blue}}=1/10,{\eta }_{\mathrm{DISK},\mathrm{blue}}=1/3,\\ {\eta }_{\mathrm{dyn},\mathrm{red}} & = & 9/10,{\eta }_{\mathrm{DISK},\mathrm{red}}=2/3,\end{array}\end{eqnarray} \tag{ 15 }$$

whereas, for t > t_coll,

$$\begin{eqnarray}\begin{array}{rcl}{\eta }_{\mathrm{HMNS},\mathrm{blue}} & = & 0,{\eta }_{\mathrm{dyn},\mathrm{blue}}=0,{\eta }_{\mathrm{DISK},\mathrm{blue}}=1/12,\\ {\eta }_{\mathrm{dyn},\mathrm{red}} & = & 0,{\eta }_{\mathrm{DISK},\mathrm{red}}=11/12.\end{array}\end{eqnarray} \tag{ 16 }$$

In this way, we can constrain the lifetime of the remnant by determining when the ejected mass in the blue component reaches the estimated observational value, i.e., ${M}_{\mathrm{ej},\mathrm{blue}}({t}_{\mathrm{coll}})=0.025\,{M}_{\odot }$. Doing so, we deduce that the collapse time as constrained from the ejecta and the condition that ${t}_{\mathrm{coll}}\leqslant {t}_{\mathrm{del}}$ must have been at

$$\begin{eqnarray}&&{t}_{\mathrm{coll}}={1.14}_{-0.50}^{+0.60}\,{\rm{s}}.\end{eqnarray} \tag{ 17 }$$

This is illustrated in Figure 3, where we show the ejected mass in the blue (blue solid line) and red (red solid line) components, respectively. Note that we set the most likely time of collapse from the crossing of the ejected blue component with the corresponding observational estimate (blue dashed line). Note also that the ejected red component is always less than the blue component, and that the latter is almost three times larger; only at ≈2.5 s have the red and blue components been ejected in equal amounts, and at $\approx 200\,{\rm{s}}$, the red component has reached the nominal value of ${M}_{\mathrm{ej},\mathrm{red}}\approx 0.04\,{M}_{\odot }$ deduced from the observations. The very slow increase in the red component is due to the fact that after 9 s, the torus has been mostly accreted, thus providing only a very small contribution to the secularly ejected matter (see Equation (11) and Fernández et al. 2018).

Zoom In Zoom Out Reset image size

The uncertainty in t_coll reported in Equation (17) was obtained from the corresponding uncertainty in the ejected blue-component mass. More specifically, the lower limit on t_coll was obtained when the upper limit on the blue-component mass crossed the estimate inferred from kilonova observations. The upper limit, on the other hand, is limited by the delay time of t_del = 1.74 s, since t_coll < t_del. When a constant fractional error δ_i in the ejected mass M_i(t) from a given mass-ejection channel i was considered, the final uncertainty in the blue-component mass was obtained by adding in quadrature the uncertainties for the different components that contributed to the blue-ejecta mass, namely, Equation (13), such that ${\sigma }_{{M}_{\mathrm{ej},\mathrm{blue}}}(t)\,=\sqrt{{\sum }_{i}{\sigma }_{{M}_{i}}^{2}(t)}$, where ${\sigma }_{{M}_{i}}(t):= {\delta }_{i}{\eta }_{i}{M}_{i}(t)$ and η_i is their fractional contribution. Similarly, when considering an uncertainty in the power-law exponent of the mass-ejection rate (Equation (9)), with $\dot{M}$ ∝ t^α, the uncertainty in the ejected mass from this channel was calculated as ${\sigma }_{M}(t)=| (\partial M(t)/\partial \alpha ){\sigma }_{\alpha }|$, where σ_α is the uncertainty in the exponent.

We note that a more conservative lower limit on the collapse time is obtained, with ${t}_{\mathrm{coll}}={1.14}_{-0.84}^{+0.60}\,{\rm{s}}$, when we consider a 30% fractional uncertainty in the power-law exponent in the mass-ejection rate (Equation (9)) from MHD-driven winds. In this case, due to a smaller exponent in the power law, the mass-ejection rate increases at times t < 1 s, as compared to the case with a constant fractional error in the ejected mass. This causes the upper limit on the blue-component mass to cross the blue-component mass limit (as shown in Figure 3) at an earlier time. Hence, there is a smaller lower limit on the collapse time; see also Figure 7.

If, in the previous sections, we have set constraints on the collapse time of the remnant in GW170817 from considerations of the ejected mass, in this section, we set a similar constraint following a complementary route that is instead based on considerations about the time needed by the jet launched by the BH to break out. Essential in this line of argument is the velocity of the ejected matter, which, for simplicity, we assume to be a quasi-isotropic outflow with its front at radius r = R_ej(t) = β_maxct expanding into the interstellar medium at a characteristic dimensionless velocity β_max $:=$ v_max/c. The mean rest-mass density of the ejecta at any given radius r < R_ej(t) can be parameterized as a power law (e.g., Nagakura et al. 2014; Matsumoto & Kimura 2018), such that

$$\begin{eqnarray}&&{\rho }_{\mathrm{ej}}(r\lt {R}_{\mathrm{ej}},t)=\displaystyle \frac{(3-k)}{4\pi }\displaystyle \frac{{M}_{\mathrm{ej},\mathrm{blue}}(t)}{{R}_{\mathrm{ej}}^{3}(t)}{\left[\displaystyle \frac{r}{{R}_{\mathrm{ej}}(t)}\right]}^{-k},k\lt 3,\end{eqnarray} \tag{ 18 }$$

where M_ej,blue(t) is the mass of the blue ejecta at time t. Our fiducial model assumes k = 2 for the density profile of the ejecta along the rotational axis. This profile has been shown to emerge in the long-term simulations of Fujibayashi et al. (2018). Some modifications to this profile will occur if the remnant is endowed with strong global magnetic fields, which would alter the mass-ejection rate due to the MHD winds in the polar regions (e.g., Metzger et al. 2007; Siegel et al. 2014).

In our phenomenological modeling, at t = t_coll, ejection of matter from the remnant essentially stops, and the central BH launches a relativistic jet. Naturally, there may be a small delay (of the order of the dynamical time) between the termination of mass ejection and the launching of the jet as mass starts to accrete onto the BH; because of its smallness, we neglect such a delay and assume that the jet is launched at $t={t}_{\mathrm{coll}}$.

In order for the relativistic jet to reach the radius r = R_γ, where it undergoes some internal dissipation and produces the prompt emission, it must first burrow its way through the intervening ejecta. We assume that the jet is launched at relativistic speeds with Lorentz factor Γ_j and power per unit solid angle

$$\begin{eqnarray}&&\displaystyle \frac{{{dL}}_{j}}{d{\rm{\Omega }}}={r}^{2}{{\rm{\Gamma }}}_{j}^{2}{h}_{j}{\rho }_{j}{c}^{3}\simeq \displaystyle \frac{{L}_{j}}{{\rm{\Delta }}{{\rm{\Omega }}}_{j}},\end{eqnarray} \tag{ 19 }$$

where h_j and ρ_j are the jet’s specific enthalpy and rest-mass density, respectively. The last equality assumes that the jet power is approximately uniform and confined to a solid angle ΔΩ_j.

As the jet expands, its interaction with the surrounding ejecta slows it down, such that the jet’s head moves only mildly relativistically with a Lorentz factor Γ_h $:=$ (1 − ${\beta }_{h}^{2}$)^−1/2 < Γ_j, where β_h $:=$ v_h/c is its dimensionless velocity. A double-shock structure develops at the jet’s head (see, e.g., Figure 1 in Bromberg et al. 2011), where a forward shock propagates into the cold circum-merger ejecta and a reverse shock moves into the cold jet, which slows it down. A contact discontinuity separates the newly shocked jet material from the shocked circum-merger ejecta. The velocity of the jet head can be calculated by balancing the momentum flux density in the frame of the jet’s head (e.g., Begelman & Cioffi 1989; Matzner 2003; Bromberg et al. 2011),

$$\begin{eqnarray}&&{\rho }_{j}{h}_{j}{\left[{{\rm{\Gamma }}}_{j}{{\rm{\Gamma }}}_{h}({\beta }_{j}-{\beta }_{h})\right]}^{2}={\rho }_{\mathrm{ej}}{h}_{\mathrm{ej}}{\left[{{\rm{\Gamma }}}_{\mathrm{ej}}{{\rm{\Gamma }}}_{h}({\beta }_{h}-{\beta }_{\mathrm{ej}})\right]}^{2},\end{eqnarray} \tag{ 20 }$$

where Γ_ej ≃ 1 is the Lorentz factor of the ejecta and h_ej is its specific enthalpy. Here the unshocked jet and circum-merger ejecta are considered cold, and therefore they contribute negligible pressure. From this relation, it is straightforward to obtain the head velocity (e.g., Matzner 2003; Bromberg et al. 2011; Murguia-Berthier et al. 2014, 2017),

$$\begin{eqnarray}&&{\beta }_{h}=\displaystyle \frac{{\beta }_{j}+{\tilde{L}}^{-1/2}{\beta }_{\mathrm{ej}}}{1+{\tilde{L}}^{-1/2}},\end{eqnarray} \tag{ 21 }$$

where the dimensionless jet luminosity, $\tilde{L}$, which determines the jet dynamics, is given by the ratio of the jet energy density to that of the surrounding ejecta,

$$\begin{eqnarray}&&\tilde{L}:= \displaystyle \frac{{{\rm{\Gamma }}}_{j}^{2}{h}_{j}{\rho }_{j}}{{\rho }_{\mathrm{ej}}}\simeq \displaystyle \frac{{L}_{j}}{{{\rm{\Sigma }}}_{j}{\rho }_{\mathrm{ej}}{c}^{3}}.\end{eqnarray} \tag{ 22 }$$

Here Σ_j $:=$ ${z}_{h}^{2}$ΔΩ_j describes the area of the working surface at the head of the jet as it bores its way out of the ejecta, where ${z}_{h}(t):= \int {\beta }_{h}{cdt}$ is the position of the jet’s head. For a conical (i.e., uncollimated) flow, ΔΩ_j = $\pi {\theta }_{0}^{2}$, where θ₀ ≪ 1 is the initial jet half-opening angle. In this case, the jet expands ballistically, with Σ_j ∝ ${z}_{h}^{2}$, and the feedback of the hot cocoon on the jet itself is ignored (see discussion below). In principle, the jet power can vary over time due to the intermittent activity of the central BH. For simplicity, however, we assume that the jet power remains constant in time, so that L_j effectively represents the time-averaged jet power.

Because of the high-pressure shocked ejecta downstream of the forward shock, the shocked jet material and ejecta are able to move sideways away from the jet’s head until the angular size of the shocked causal region, ${{\rm{\Gamma }}}_{h}^{-1}$, becomes smaller than the instantaneous half-opening angle of the jet, i.e., ${{\rm{\Gamma }}}_{h}^{-1}\lt {\theta }_{j}(t)$. Before this condition is met, a jet-cocoon system is produced, where the pressure in the hot cocoon collimates the relativistic jet and reduces its semi-aperture from its initial value, θ₀ (see Figure 1 of Harrison et al. 2018). Here we follow the treatment of Bromberg et al. (2011), which applies for a static ambient medium into which the jet expands, to write the equations that govern the cross section of the working surface at the jet’s head. Of course, this implies that the local velocity of the ejecta, which is modeled here as expanding homologously with dimensionless velocity

$$\begin{eqnarray}&&{\beta }_{\mathrm{ej}}(r\lt {R}_{\mathrm{ej}},t)={\beta }_{\max }\left(\displaystyle \frac{r}{{R}_{\mathrm{ej}}(t)}\right),\end{eqnarray} \tag{ 23 }$$

must be significantly smaller than that of the jet’s head, i.e., ${\beta }_{\mathrm{ej}}\ll {\beta }_{h}$. The assumption of homologous expansion, as compared to a constant velocity wind, is useful here, since at a smaller radius, the ejecta velocity will also be smaller when the jet is launched. Given that the jet head will initially be subrelativistic, a slower ejecta velocity helps to preserve the condition of having β_ej < β_h throughout most of the jet propagation inside the expanding ejecta. Under this assumption, the cocoon pressure, assumed here to be dominated by radiation, depends on the energy injected by the jet, such that

$$\begin{eqnarray}&&{p}_{c}=\displaystyle \frac{1}{3}\displaystyle \frac{{E}_{c}}{{V}_{c}}=\displaystyle \frac{1}{3}\displaystyle \frac{{L}_{j}\int (1-{\beta }_{h}){dt}}{\pi {r}_{c}^{2}(t){z}_{h}(t)},\end{eqnarray} \tag{ 24 }$$

where ${r}_{c}(t):= \int {\beta }_{c}{cdt}$ is the cocoon radius expanding with dimensionless velocity

$$\begin{eqnarray}&&{\beta }_{c}:= {\left[\displaystyle \frac{{p}_{c}}{{\rho }_{\mathrm{ej}}({z}_{h}){c}^{2}}\right]}^{1/2}.\end{eqnarray} \tag{ 25 }$$

Our model assumes that the pressure in the cocoon is uniform and its geometry is that of a cylinder. Therefore, under this simplification, we study the dynamics of the jet’s head along the jet axis only, z_h(t). The jet is injected at a radius z_inj, after which it starts to be influenced by the cocoon as the pressure in the cocoon builds up. The location z_cj > z_inj, where the cocoon collimates the jet and after which the jet maintains a constant cylindrical radius, is derived from the condition of pressure equilibrium, p_c = p_j, which yields

$$\begin{eqnarray}&&{z}_{{cj}}\simeq {\left[\displaystyle \frac{{L}_{j}}{\pi {{cp}}_{c}}\right]}^{1/2}.\end{eqnarray} \tag{ 26 }$$

However, the cocoon starts collimating the jet at an even smaller distance of z ≈ z_cj/2, at which point the area of the working surface can be expressed as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{j}=\pi {r}_{j}^{2}\simeq \displaystyle \frac{1}{4}\pi {\hat{z}}^{2}{\theta }_{0}^{2}\simeq \displaystyle \frac{1}{4}\displaystyle \frac{{L}_{j}{\theta }_{0}^{2}}{{{cp}}_{c}},\end{eqnarray} \tag{ 27 }$$

where r_j is the cylindrical radius of the jet’s head.

The set of equations presented above, Equations (21)–(27), which were introduced mostly by Bromberg et al. (2011), represent a closed system that we have solved numerically to follow the evolution of the jet-cocoon system. Under the simplifying assumption that the velocities of the different components in the system remain constant, a linearized set of equations can be obtained that can be integrated analytically. Then, for a static medium, it can be further shown that the jet interacts strongly with the surrounding medium and remains collimated for $\tilde{L}\lesssim {\theta }_{0}^{-4/3}$ (Bromberg et al. 2011). In an expanding medium, this result does not strictly apply, and deviations from it are expected.

As discussed above, these equations are based on several simplifying assumptions, e.g., of uniform pressure in the cocoon and its cylindrical geometry, or the fact that only the first collimation shock, which is much stronger than the subsequent ones (Mizuno et al. 2015), is accounted for. These assumptions, however, were verified to provide a reasonable approximation in a comparison with two- and three-dimensional hydrodynamic simulations (Harrison et al. 2018).

In Figure 4, we show the dynamics of the jet-cocoon system expanding inside the quasi-isotropic ejecta from the HMNS with β_max = 0.3, where the relativistic jet is injected at $t={t}_{\mathrm{coll}}=1\,{\rm{s}}$ with a bulk Lorentz factor Γ_j = 10, corresponding to an initial half-opening angle of ${\theta }_{0}={{\rm{\Gamma }}}_{j}^{-1}\simeq 5\buildrel{\circ}\over{.} 7$. Here and in the rest of the paper, we assume a reference value for the single-sided jet power of L_j = E_j/t₉₀ = 10⁵⁰ erg s⁻¹, where E_j is the single-sided true jet energy and t₉₀ is the duration of the prompt emission within which 90% of the fluence is accumulated. In the case of GRB 170817A, the prompt gamma-ray emission lasted for t₉₀ = 2.0 ± 0.5 s (Abbott et al. 2017a; Goldstein et al. 2017), where the main hard spike had a duration of ∼0.5 s and was followed by a tail of softer emission. Many numerical simulations (e.g., Granot et al. 2018; Lazzati et al. 2018; Mooley et al. 2018; Xie et al. 2018) of this event find a true jet energy of E_j ∼ 10⁵⁰ erg that motivated the assumed value of the jet power in this work.

Zoom In Zoom Out Reset image size

Immediately after injection, the velocity of the jet’s head slows down significantly to β_h ≪ 1 as it encounters the ejected mass. The pressure in the surrounding cocoon starts to build up and the jet begins to be collimated, as seen from the reduction in the cylindrical radius r_j of the jet (see Figure 4 for $t-{t}_{\mathrm{coll}}\lesssim {10}^{-4}\,{\rm{s}}$). This also aids in accelerating the jet due to the shrinkage in the area of the working surface at the jet’s head. This phase, however, only lasts for a very short time, and, as the cocoon expands, shown by the increase in r_c, the drop in its pressure is insufficient to maintain collimation of the jet. At this point, the size of the jet’s head starts to expand, and the jet essentially propagates uncollimated. Therefore, the jet only encounters a single recollimation shock, and that is only due to the inherent simplicity of the underlying model. In principle, the jet may experience multiple recollimation shocks before breakout, which would cause r_j to oscillate. Nevertheless, the half-opening angle of the jet is now smaller than its initial value, θ_j < θ₀, and the jet maintains it until it clears the ejecta and breaks out when z_h(t) = R_ej(t) at a breakout time t_br ≈ 0.65 s in the case considered in Figure 4. After the first recollimation shock, we find that the jet head maintains an approximately constant velocity, with β_h ≃ 0.7 and constant $\tilde{L}\simeq 6\lt {\theta }_{0}^{-4/3}$, until just before breakout, after which point the jet head becomes relativistic with β_h ≃ 1.

Also shown with dashed lines in Figure 4 is the corresponding dynamics in the case in which the jet is forced to remain uncollimated; i.e., it maintains the same opening angle as the initial one. In this case, it is clear that the jet takes longer to break out (by a factor of ≃2 in time) and might even be choked if the engine shuts off before that time. Since the solid angle is fixed in this case, the radius of the working surface grows with the vertical distance traveled by the jet, i.e., r_j ∝ z_h at all times.

Figure 5 shows the jet-breakout time t_br as a function of the collapse time and how t_br depends on the different ejecta velocities β_max and jet initial half-opening angles. The mass-ejection rates are those discussed in Section 3, so that a longer t_coll results in a larger ejected mass in the path of the jet and thus to a delayed jet breakout. The sharp drop in t_br at early times simply shows that the ejecta has not had enough time to expand, and the jet breaks out almost immediately. Note that the rate at which the ejecta is expanding has a significant effect on the breakout time (both scales are logarithmic), with longer t_br in the faster-expanding ejecta. A similar effect can also be produced if the jet has a larger initial half-opening angle θ₀, which increases the area of the working surface at the jet’s head and causes the jet to slow down. The small jump at t_coll = 10⁻² s is due to the enhanced rate of mass ejection by the newly formed HMNS at t = t_dyn ∼ 10 ms post-merger (see left panels of Figure 2).

Zoom In Zoom Out Reset image size

After breakout, the jet very rapidly accelerates, and its energetic core expands with ultrarelativistic speeds, Γ_j ≫ 1, where it is surrounded by less energetic and relatively slower-moving material. This phase of the jet dynamics is outside the scope of this work and is not calculated by the one-dimensional model used above. In the next section, and for simplicity, the jet is treated as an expanding conical flow, so it is convenient to use the radial coordinate r in place of the cylindrical coordinate z.

The delay of 1.74 s in the onset of the prompt emission measured in GRB 170817A can be explained by a sum of three different timescales (e.g., Granot et al. 2017). First, and most obviously, the launching of the jet is delayed by the HMNS that will survive for a timescale ${t}_{\mathrm{coll}}$ (see Section 4), after which it will collapse to a BH. Second, the jet has to clear the ejecta before it can radiate the prompt emission photons, otherwise the Thomson optical depth of the medium will be too large and will suppress any high-energy emission. The jet clears the ejecta over the breakout time t_br that, in turn, depends on t_coll, as well as the total ejected mass, the density profile, and the jet-launching properties. Third, since GWs propagate at the speed of light (Abbott et al. 2017a) and the jet propagates at slightly lower velocities, there is always a radial time delay t_R between the two components to reach the prompt gamma-ray emission radius R_γ, i.e.,

$$\begin{eqnarray}&&{t}_{R}=(1+{z}_{c})\displaystyle \frac{{R}_{\gamma }}{2{{\rm{\Gamma }}}^{2}c}\simeq \displaystyle \frac{{R}_{\gamma }}{2{{\rm{\Gamma }}}^{2}c},\end{eqnarray} \tag{ 28 }$$

where, in the second equality, we have accounted for the fact that the considered cosmological redshift z_c of GW170817 is essentially zero. As a result of these considerations, the total delay time can be written as

$$\begin{eqnarray}&&{t}_{\mathrm{del}}=1.74\pm 0.05\,{\rm{s}}={t}_{\mathrm{coll}}+{t}_{\mathrm{br}}({t}_{\mathrm{coll}})+{t}_{R},\end{eqnarray} \tag{ 29 }$$

where we have explicitly expressed that the jet-breakout time depends on the collapse time of the remnant. The breakout time considered above is that for jet material that breaks out of the ejecta along the jet symmetry axis. In our model, we assume that after breakout, the relativistic jet core is surrounded by a sheath of lower-energy and slower-moving material (e.g., Kathirgamaraju et al. 2018) that undergoes dissipation. Alternatively, some works (e.g., Kasliwal et al. 2017; Gottlieb et al. 2018; Nakar et al. 2018) have argued that the prompt gamma-ray emission was produced by the breakout of the shock driven by the expanding cocoon from the subrelativistic circum-merger ejecta. In this case, the breakout time at an angle θ_obs away from the jet symmetry axis will be longer, and consequently, the collapse time will be shorter.

At r = R_γ, the jet experiences internal dissipation and produces the prompt gamma-ray emission. The underlying energy dissipation mechanism is still uncertain (see, e.g., Kumar & Zhang 2015 for a review). If the jet is kinetic energy–dominated, energy can be dissipated in internal collisions of mass shells ejected by the BH with variable Lorentz factors, where faster shells coming from the engine would catch up with slower shells at the dissipation radius. Alternatively, if the jet is Poynting flux–dominated, magnetic energy is the main energy reservoir, which can be tapped via magnetic reconnection or MHD instabilities.

The condition for causality dictates that R_γ ≲ 2Γ²cδt_obs, where $\delta {t}_{\mathrm{obs}}$ is the observed flux-variability time (e.g., Sari & Piran 1997). The same condition is also obtained in the internal-collision scenario, in which case Γ is the Lorentz factor of the slower-moving shell and δt_obs represents the engine-variability time. The observed variability time of the GRB light curve can be associated with either the rise time of the pulse, the decay time, or their combination, depending on their contribution to the total pulse width, which can be expressed as Δt_pulse $:=$ δt_r + δt_θ. To be more specific, consider a thin, spherical, relativistic shell that emits over a range of radii (R_γ − δR) ≤ r ≤ R_γ. The emission from this shell will be temporally spread out over the timescale $\delta {t}_{r}=\delta R/2{{\rm{\Gamma }}}^{2}c$ and will represent the rise time of the pulse. On the other hand, the decay time of the pulse occurs due to light travel-time effects across a curved emitting surface, also referred to as the “angular time,” such that $\delta {t}_{\theta }={R}_{\gamma }(1-\cos \tilde{\theta })/c$, where the term in parentheses represents the extra distance traveled by the photon that was emitted at an angle $\tilde{\theta }$ away from the line of sight. For a relativistic shell, the emission is beamed within a cone of angle $\tilde{\theta }\simeq {{\rm{\Gamma }}}^{-1}\ll 1$, which yields $\delta {t}_{\theta }\simeq {R}_{\gamma }/2{{\rm{\Gamma }}}^{2}c$. This finally yields (e.g., Genet & Granot 2009)

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{\mathrm{pulse}}=\displaystyle \frac{{R}_{\gamma }}{2{{\rm{\Gamma }}}^{2}c}\left(1+\displaystyle \frac{\delta R}{{R}_{\gamma }}\right).\end{eqnarray} \tag{ 30 }$$

The minimum variability time for the prompt gamma-ray emission in GRB 170817A, which would correspond to the rise time to the pulse peak, i.e., δt_r, measured by Fermi/GBM in the energy range of 50–300 keV, is δt_min = 0.125 ± 0.064 s, and the pulse decay time is δt_θ ≈ 0.5 s (Abbott et al. 2017a; Goldstein et al. 2017). On the other hand, INTEGRAL/SPI-ACS observed the main hard spike spread out over two time bins of $50\,\mathrm{ms}$, with a total duration of ${\rm{\Delta }}{t}_{\mathrm{pulse}}\lesssim 100\,\mathrm{ms}$ in the ∼75–2000 keV energy range (Savchenko et al. 2017). These two measurements suggest that the pulse duration is energy-dependent and that the decay time of the pulse $\delta {t}_{\theta }\lesssim 100\,\mathrm{ms}$ at high energies. Unless $\delta R/{R}_{\gamma }\approx 1$, the pulse rise time cannot be used to robustly determine the emission radius. Therefore, we associate the variability time with the angular time, δt_obs ≈ δt_θ, and use it to constrain the emission radius given an estimate of Γ.

Almost all distant GRBs are observed within the semi-aperture of the bright core, such that θ_obs ≲ θ_c, where θ_obs is the viewing angle and θ_c is the half-opening angle of the bright core (both angles are measured from the jet symmetry axis). In case the structure of the jet can be modeled as having sharp edges, then θ_c = θ_j, where θ_j is the half-opening angle of the jet. In these GRBs, compactness arguments (e.g., Lithwick & Sari 2001; Gill & Granot 2018b) dictate that Γ > 100, otherwise opacity due to γγ-annihilation (i.e., γγ → e⁻ + e⁺) would prevent the radiation of gamma-ray photons. Modeling of the radio/optical/X-ray afterglow emission from GRB 170817A, using both semi-analytical models (e.g., Gill & Granot 2018a; Lamb & Kobayashi 2018; Troja et al. 2018) and numerical simulations (e.g., Granot et al. 2018; Lazzati et al. 2018; Margutti et al. 2018; Kathirgamaraju et al. 2019) of structured jets, has revealed that the line of sight is significantly off-axis, with 20° ≲ θ_obs ≲ 28° and θ_c ∼ 5°.

The prompt gamma-ray emission of GRB 170817A at a distance of $40\,\mathrm{Mpc}$, D₄₀, with isotropic equivalent energy of ${E}_{\gamma ,\mathrm{iso}}=(5.36\pm 0.38)\times {10}^{46}{D}_{40}^{2}$ erg, was ∼3–4 orders of magnitude lower than what is typically observed from distant short-hard GRBs. This fact alone already hinted at an off-axis jet. Early analysis of the burst energetics showed that the prompt gamma-ray emission was dominated by low-energy material along the line of sight (e.g., Granot et al. 2017) and with only modest Γ. By applying the compactness arguments, it was further realized that the emission region must be moving toward the observer with Γ ≳ 5 (Matsumoto et al. 2019). This lower limit on Γ can now be used to obtain a constraint on the dissipation radius (also see Beloborodov et al. 2018),

$$\begin{eqnarray}&&{R}_{\gamma }\gtrsim 7.5\times {10}^{11}{\left(\displaystyle \frac{{\rm{\Gamma }}}{5}\right)}^{2}\left(\displaystyle \frac{\delta {t}_{\mathrm{obs}}}{0.5\,{\rm{s}}}\right)\,\mathrm{cm}.\end{eqnarray} \tag{ 31 }$$

If the decay time of the pulse is indeed as short as ∼100 ms, then the dissipation radius can be smaller by a factor of ∼5.

The angular time also yields the radial delay time t_R = δt_θ between the relativistic ejecta and the GWs to arrive at the emission radius R_γ. From Equation (29), we find that a longer t_R would yield a shorter collapse time. Therefore, to obtain a conservative estimate of t_coll, we use ${t}_{R}\approx 0.5\,{\rm{s}}$. Using this, we next numerically solve Equation (29) for a given jet luminosity and mass-ejection rate to obtain the collapse time of the HMNS as a function of t_del. This is shown in Figure 6, where the different lines refer to different values of β_max and the jet power, using as reference values k = 2 and θ₀ ≃ 57. The thin vertical line marks the t_del = 1.74 s that was observed for GRB 170817A, and it shows a collapse time of t_coll ≃ 0.75 s obtained for β_max = 0.3 and L_j = 10⁵⁰ erg s⁻¹. In this way, we deduce that the collapse time as constrained from the EM delay must have been at

$$\begin{eqnarray}&&{t}_{\mathrm{coll}}=0.82\pm 0.15\,{\rm{s}},\end{eqnarray} \tag{ 32 }$$

where the range in the variance is dictated by considering ${\beta }_{\max }=0.2$ or 0.3, respectively. This result is robust in regard to the jet power and changes only weakly as the jet power is varied by an order of magnitude. Furthermore, varying the density profile by using different values of the power-law index k in Equation (18) yields only small differences.

Zoom In Zoom Out Reset image size

Figure 7 offers a synthetic and yet comprehensive summary of the results of this paper by comparing the estimates of the collapse time of the remnant of GW170817. More specifically, the horizontal orange line reports the estimated collapse time, t_coll, as deduced from the ejected matter producing the blue-kilonova emission (see Equation (17)). Similarly, the vertical green line reports the estimate of t_coll needed for the jet to burrow its way through the intervening ejecta (see Equation (32)). In both cases, the corresponding uncertainty regions are shown with shaded regions of the corresponding color; the light orange shaded area refers to the more conservative lower limit on the collapse time obtained when considering an uncertainty in the power-law exponent in the mass-ejection rate (see Section 4). Also shown is the delay time between the GW and EM emission, t_del (black solid line), and the corresponding uncertainty (gray shaded area). Combining this information, we deduce that the best estimate for the collapse time of the remnant of GW170817 is

$$\begin{eqnarray}&&{t}_{\mathrm{coll}}={0.98}_{-0.26}^{+0.31}\,{\rm{s}}.\end{eqnarray} \tag{ 33 }$$

Zoom In Zoom Out Reset image size

A more conservative lower limit on the collapse time is obtained, with t_coll = ${0.98}_{-0.43}^{+0.31}$ s, when we consider a 30% uncertainty in the power-law exponent in the mass-ejection rate (Equation (9)) from MHD-driven winds. The level of uncertainty considered here in the power-law exponent is more conservative than that inferred from simulations, which is ≈11% and for which we find a less conservative lower limit on the collapse time as compared to that obtained from a constant fractional error on the ejected mass. Due to the small number of simulations that have explored the MHD-driven mass-ejection channel, the magnitude of the uncertainty in the power-law exponent is poorly known at best. Therefore, we adopt as a reference the estimate in Equation (33) for the collapse time in this work.

Such a survival time for the remnant is considerably longer than the one explored in three-dimensional general-relativistic simulations of BNS mergers, which so far have been limited to simulations with ${t}_{\mathrm{coll}}\sim 10\mbox{--}100$ ms. Since on timescales of ${t}_{\mathrm{coll}}\lesssim 100\,\mathrm{ms}$, the most efficient process to remove energy and angular momentum is the emission of GWs, the implication stemming from Equation (33) is that other, less efficient mechanisms need to be invoked to extract angular momentum on such long timescales. At the same time, magnetic spin-down provides a simple and compelling explanation. Using Equation (4) for a remnant mass of M ≈ 2.7 M_⊙, it is simple to deduce that a spin-down timescale of the order of the collapse time, i.e., ${t}_{\mathrm{sd}}\sim {t}_{\mathrm{coll}}\sim 1\,{\rm{s}}$, can be accomplished with magnetar-strength surface magnetic fields, B_s ≃ 10¹⁶ G. In addition, such surface magnetic fields would naturally provide an enhanced rate of mass loss in the polar region due to magnetocentrifugal emission (Metzger et al. 2018) and a high rate of expansion of the ejecta, with $0.2\lesssim {\beta }_{\max }\lesssim 0.3$. Interestingly, the collapse of the remnant when such large magnetic fields are present would also lead to an explosive outflow (Salafia et al. 2017; Nathanail 2018a; see also Nathanail et al. 2018 for a discussion of the observable signatures).

At the same time, magnetic fields in the post-merger object cannot reach the largest values allowed by energy equipartition, i.e., ${B}_{\mathrm{int}}\simeq {10}^{17}-{10}^{18}$ G. Such ultrastrong magnetic fields would, in fact, prolong the breakout time, as the jet would have to plow through more material and spend more time catching up with the expanding ejecta front. Furthermore, the duration distribution of short-hard GRBs is suggestive of a characteristic jet-breakout time of ${t}_{\mathrm{br}}\simeq 0.4\,{\rm{s}}$ (Moharana & Piran 2017); assuming that this timescale is indeed common to all short GRBs, this would limit the collapse time to t_coll ∼ 0.7 s for a typical jet power of L_j = 10⁵⁰ erg s⁻¹. The total time delay in this case would be of the order of t_del ≳ 1 s, which can be confirmed by future multimessenger detections of BNS mergers. In summary, strong (but not ultrastrong) magnetic fields represent a natural explanation for the high mass-ejection rates and long delay times between the GW and the gamma-ray emission.

Although the production of such magnetic fields still represents a serious challenge for direct numerical simulations, the existence of a long-lived HMNS can, in principle, also be confirmed by the detection of GWs after the merger. The GW signal, in fact, has clear spectroscopic features, with precise oscillation frequencies in the range ∼2–4 kHz that can be related to the EOS and the properties of the progenitors (e.g., Bauswein & Janka 2012; Takami et al. 2014; Rezzolla & Takami 2016). As mentioned earlier, in the case of GW170817, the search for GWs from the remnant in the signals from both the Advanced LIGO and Virgo detectors came out empty (Abbott et al. 2017d; but see van Putten & Della Valle 2019 for a different claim). Future detections of BNS mergers by advanced or third-generation detectors such as the Einstein Telescope (Punturo et al. 2010; Sathyaprakash et al. 2012) will provide better estimates of the collapse time of the remnant and help clarify what is the dominant process of removal of angular momentum on such long timescales.

The GW event GW170817 and its EM counterpart, GRB 170817A, have not only marked the birth of multimessenger GW astronomy but also confirmed that mergers of binary systems of neutron stars are responsible for the phenomenology of short GRBs. Unfortunately, because of its high-frequency properties during the final stages of the inspiral, the GW signal from GW170817 became invisible to the LIGO/Virgo detectors before the actual merger took place. As a result, we are unable to extract direct GW information about the object that was produced by the merger and was most likely a metastable HMNS.

Making use of the observational features of GRB 170817A, namely, that the main hard pulse of gamma-ray emission had a duration of $\sim 0.5\,{\rm{s}}$ and its onset was delayed by t_del ≈ 1.74 s with respect to the GW chirp signal, and assuming that the prompt gamma-ray emission is ultimately the result of a relativistic jet powered by a rotating BH, we have constrained the lifetime of the remnant of GW170817. In particular, we have used the properties of the kilonova emission and the delay time between the GW chirp signal and prompt gamma-ray emission onset to constrain the collapse time of the HMNS. We have argued that in order to produce the requisite mass of ${M}_{\mathrm{ej},\mathrm{blue}}\approx 0.025\,{M}_{\odot }$ in the lanthanide-poor ejecta that gave rise to the rapidly fading bluer emission in the UV and optical at early times, the collapse time of the HMNS cannot be different from about 1 s. More importantly, we have reached a very similar conclusion from an independent line of argument, where we model the dynamical evolution of the relativistic jet launched after the HMNS collapses and bores its way out of the dynamical ejecta.

To derive the estimates of the collapse time from the ejected matter, we have combined into a single systematic description the various processes leading to the ejection of matter from the remnant, which include the dynamical ejecta, the neutrino, and the MHD-driven winds from both the HMNS and the accretion torus. Estimates for the mass-ejection rates from all of these channels have been extracted from numerical-relativity simulations carried out on much shorter timescales that have been extrapolated to the timescales of interest after casting them into suitable analytic forms. In this way, we have deduced that the collapse time is constrained to be ${t}_{\mathrm{coll}}={1.14}_{-0.50}^{+0.60}\,{\rm{s}}$. Similarly, to derive the collapse time starting from the delay between the GW chirp signal and the onset of the prompt gamma-ray emission, we have computed the jet-breakout time for a fiducial set of parameters for the homologously expanding ejecta, e.g., an expansion velocity 0.2 ≤ β_max ≤ 0.3 and a jet luminosity $0.1\leqslant {L}_{j,50}\leqslant 10$. These estimates are less uncertain than those related to the ejection of mass and have allowed us to deduce that the collapse time is constrained to be t_coll = 0.82 ± 0.15 s. Combining these two results and the corresponding uncertainties, we conclude that the remnant of GW170817 must have collapsed to a rotating BH at about t_coll = ${0.98}_{-0.26}^{+0.31}$ s after merger. Interestingly, this estimate represents the first attempt to use GW and astronomical observations to constrain the properties of the gravitational collapse to a BH.

Looking at the future, the collection of new multimessenger GW detections in the coming months will provide additional important information on the delay between the GW and the EM signals and hence better constrain the collapse time of the HMNS for a given mass of the binary system. In turn, this will bring clarity to many open questions, not only on the merger of BNSs but also on the mechanism leading to the production of a short GRB. First, they will help to narrow down the search for the correct EOS, since the collapse time depends on it, with some EOSs yielding only a short-lived HMNS (i.e., ${t}_{\mathrm{coll}}\lesssim 10\,\mathrm{ms}$), while others produce a longer-lived remnant. Second, it will provide a better handle on the dominant dissipative processes; e.g., energy and angular momentum losses due to GW radiation dominate for t_coll ≲ 100 ms, whereas magnetic braking due to strong internal fields dominates thereafter. Third, since the large majority of the ejected matter leading to the blue-kilonova signal is expelled before the remnant collapses to a BH, a clear correlation between t_coll and M_ej,blue is expected and can be used to constrain the properties of the remnant. Finally, if a long-lived HMNS is produced, the collapse time will yield a constraint on the internal magnetic field and provide further insight into the physics of the dissipative effects that drive the HMNS to uniform rotation. In future work, we plan to explore several of these issues both through self-consistent numerical simulations and via semi-analytical modeling.

It is a pleasure to thank Jonathan Granot, Elias Most, and David Radice for useful discussions and input. Support comes from the ERC synergy grant “BlackHoleCam: Imaging the Event Horizon of Black Holes” (grant No. 610058), “PHAROS,” COST Action CA16214; the LOEWE-Program in HIC for FAIR; and the European Union’s Horizon 2020 Research and Innovation Programme (grant 671698; call FETHPC-1-2014, project ExaHyPE).