Jet Dynamics in Compact Object Mergers: GW170817 Likely Had a Successful Jet

The first detection of merging neutron stars by gravitational waves (GW170817; Abbott et al. 2017c) was accompanied by a wealth of electromagnetic counterparts. Notable was the first unambiguous detection of kilonova emission, which demonstrated that a significant amount of mass (≳0.04 M_⊙) had been ejected with a velocity of ∼0.2c and kinetic energy of ∼10⁵¹ erg (e.g., Abbott et al. 2017a; Coulter et al. 2017; Kasliwal et al. 2017; Shappee et al. 2017; Soares-Santos et al. 2017; Valenti et al. 2017).

Approximately 1.7 s after the merger, a low-energy short-duration gamma-ray burst (GRB) was also detected (Abbott et al. 2017b). Currently there are two leading theories concerning the origin of the gamma rays: either they are produced by shock breakout from some injected energy from a central engine (Bromberg et al. 2018; Gottlieb et al. 2018), or they are produced by less relativistic material outside the core of a traditional short-GRB jet observed off-axis (a "structured jet"; Lazzati et al. 2017, 2018; Kathirgamaraju et al. 2018). Either of these models requires relativistic material to escape through the v ∼ 0.2c ejecta, at least some fraction of which is relatively spherically distributed (Hotokezaka et al. 2013).

The shock breakout scenario is sometimes envisioned as being associated with a failed GRB jet (e.g., Gottlieb et al. 2018). In this scenario, despite the failure of the jet to escape the surrounding ejecta, because significant energy was injected into the flow in the form of a cocoon, a shock continues to propagate forward, eventually breaking out of the ejecta, accelerating to relativistic velocities and releasing gamma rays.

The thermal energy imparted to the cocoon may also be important for the quasi-thermal kilonova emission in neutron star mergers. For GW170817 in particular, there are multiple observed components to the thermal transient (Shappee et al. 2017): a blue component with faster velocities (∼0.2c), a red component with slower velocities (∼0.1c), and possibly an intermediate component (see Villar et al. 2017). It has been suggested that the blue component (particularly at early times) may be partially or fully powered by thermal energy in the cocoon (Kasliwal et al. 2017; Piro & Kollmeier 2018), rather than by the radioactive decay of neutron-rich elements, which is believed to power the majority of the thermal kilonova emission (Metzger et al. 2010b; Barnes & Kasen 2013).

The jet–ejecta interaction in neutron star mergers is somewhat analogous to a similar scenario in long-duration GRBs; there, the GRB jet needs to punch its way through a collapsing star (MacFadyen & Woosley 1999; Lazzati et al. 2012; Barnes et al. 2018). However, the scenario is not equivalent. The dynamical ejecta from a neutron star merger is a homologously expanding outflow, as opposed to the collapsar, which has a fixed size. This difference qualitatively affects the dynamics of the interaction. In this work, we perform analytical and numerical calculations of the jet–ejecta interaction, to develop new scalings for the conditions under which jets can break out of the ejecta, the resulting breakout time of the jet, and the thermalized jet energy as a function of time.

We find several important results. First, assuming a narrow engine whose energy is subdominant relative to the total kinetic energy in the dynamical ejecta, and assuming that the jet is initiated at early times compared with the breakout time, shock breakout of the cocoon only occurs if the jet is also successful. This is because the thermal energy from a failed jet is absorbed into the dynamical ejecta; because the ejecta is expanding, a low-energy shock never makes it to the edge of the outflow. Second, for a fixed total energy injection, the condition for success is independent of the duration of the engine; this is a consequence of scale invariance and of our assumption that the delay between the merger and the engine activation is negligible. Therefore, success or failure of the jet only depends on the total energy injected and the solid angle within which it is injected. Third, a consequence of this is that, for very energetic jets, the breakout time is inversely proportional to the jet power. Fourth, we find two distinct regimes, "early breakout" and "late breakout." "Late breakout" jets can break out after the engine turns off (a similar result appears to have been found for the long GRB case; see Lazzati et al. 2012). Fifth, the amount of energy that is thermalized can be estimated straightforwardly from the above considerations, leading to a robust prediction for the thermal energy in the cocoon at late times. For typical short-GRB parameters, this thermal energy is very low.

Analytical scalings for this problem are derived in Section 2. The numerical setup is briefly described in Section 3. Results of numerical calculations and comparison with the analytical scalings are presented in Section 4. Finally, these results are discussed in the context of short GRBs and GW170817 observations in Section 5.

All scalings in this section will consider pure hydrodynamics. When the word "breakout" is used in this study, it refers to a shock rapidly accelerating at the outer layers of the ejecta. This could refer to the breakout of the head of the jet, or the breakout of a more spherical shock wave associated with the cocoon. Either way, it is defined purely in terms of hydrodynamics (the radiation itself is not considered in this study).

2.1. Success versus Failure

There are multiple components to the ejecta in binary neutron star mergers, including tidal tails (e.g., Rosswog et al. 1999), quasi-spherical ejecta produced when the neutron stars first collide and shock (e.g., Hotokezaka et al. 2013), and accretion disk winds (e.g., Metzger et al. 2008). The jet is most likely to directly interact with the quasi-spherical dynamical ejecta (the case of interaction with a wind has been studied previously by Murguia-Berthier et al. 2014, 2017). This can be modeled as a homologously expanding medium, which is qualitatively distinct from the stationary medium typically considered in the context of active galactic nuclei or long GRB jets. One notable distinction is that the initial condition of a homologous outflow has no physical length scale (so long as one purely considers hydrodynamics), whereas a long GRB progenitor has a fixed radius. The GRB engine has only one timescale, the duration of the engine (assuming that there is no significant delay between the ejection of the outflow and injection of the jet). Therefore, any physical length or timescale in the solution is proportional to the engine duration; a calculation with the same energy injected but double the engine duration would yield an identical solution, but with all physical lengths and timescales multiplied by 2.

The jet–ejecta interaction in neutron star mergers is thus less like the jet–collapsar interaction in long GRBs and more like a pulsar–ejecta interaction (Chevalier & Fransson 1992; Kasen et al. 2016), in which the only relevant dynamical timescale is the pulsar spin-down time. The main distinction between this problem and the pulsar–ejecta interaction is that the energy is now injected within a small opening angle.

Given this, the success or failure of the GRB jet cannot depend on the engine duration explicitly, but only on the total injected energy, and presumably the opening angle within which the energy is injected. We write the ability of the jet to break out of the dynamical ejecta in terms of E_crit, the minimum energy required for successful breakout:

$$\begin{eqnarray}&&{E}_{j}\gt {E}_{\mathrm{crit}}\,\,(\mathrm{successful}\ \mathrm{breakout}).\end{eqnarray} \tag{ 1 }$$

If the total energy injected is below E_crit, then all of the energy is adiabatically absorbed into the kinetic energy of the dynamical ejecta. It is reasonable to expect that E_crit must be some fraction of the kinetic energy of the dynamical outflow, because it is the only other relevant energy scale in this problem (ignoring M₀c²; see below). This assumption is in agreement with the scalings calculated by Margalit et al. (2018). It is also consistent with analytical and numerical calculations that assume a spherical injection (Chevalier & Fransson 1992; Kasen et al. 2016). Exactly what this fraction is should depend on θ_j, the opening angle of the injected jet (actually the half-opening angle as defined in this study). A likely scaling is

$$\begin{eqnarray}&&{E}_{\mathrm{crit}}\propto {\theta }_{j}^{2}{E}_{\mathrm{ej}},\end{eqnarray} \tag{ 2 }$$

based on the amount of energy subtended by the opening angle θ_j. Prior to performing numerical calculations, however, it is not obvious whether this scaling is correct, as the opening angle of the cocoon is not necessarily proportional to the opening angle of the injected jet. Moreover, there is a second energy scale in the problem, M_ejc². For highly collimated jets the speed of the head of a jet propagating through a medium is set by the momentum flux of the jet, while the speed of the cocoon spreading laterally is set by the energy content of the cocoon (Begelman & Cioffi 1989). In this case one can show that the jet propagation will depend on both E_ej and M_ejc² independently. For our problem, however, the mean ejecta velocities from neutron star mergers live in a narrow range of ∼0.1c–0.3c, set by the escape velocity from a neutron star. Thus, E_ej and M_ejc² are roughly proportional to each other. For simplicity, we will thus present our results in terms of E_ej, largely keeping the mean ejecta velocity ∼(E_ej/M_ej)^1/2 fixed. Additional tests using different mean ejecta velocities of ∼0.1c–0.4c show that our conclusions are all essentially unchanged within this range (Section 4).

2.2. Breakout Time

The breakout time can be estimated as in the pulsar problem (Chevalier & Soker 1989), though with geometric corrections due to the jet opening angle. We will proceed in a little bit less detail, to calculate the scalings but without calculating the dimensionless coefficients in front, which would depend on the detailed structure of the dynamical ejecta. Our numerical results will provide a calibration of these dimensionless coefficients.

The engine injects a jet with opening angle θ_j, power L, and duration T, with total injected energy E_j = LT. In principle, there is an additional input parameter, the delay between the ejection of the cloud and initiation of the engine. All of our numerical calculations work in the limit where this timescale is very short compared to T. In this limit, there is only one input timescale in the setup of this problem: the engine duration, T. Under these specific assumptions, every length and timescale in the solution will be proportional to T. Specifically, one should be able to write the breakout time as

$$\begin{eqnarray}&&{t}_{\mathrm{bo}}={Tf}({E}_{j}),\end{eqnarray} \tag{ 3 }$$

for some dimensionless function f of the injected energy.

One could envision a thought experiment where the duration is essentially infinite, so that the engine is only specified by a constant power. In this case, the total input energy is effectively infinite, so there is always a successful breakout. After breakout, however, one is free to shut off the engine at any time, and therefore, for a given power, the minimal engine duration T necessary for breakout is given by the relationship LT = E_crit. One can turn this formula around to estimate the breakout time: t_bo ∼ E_crit/L or

$$\begin{eqnarray}&&{t}_{\mathrm{bo}}\sim {{TE}}_{\mathrm{crit}}/{E}_{j}.\end{eqnarray} \tag{ 4 }$$

This breakout time is consistent with the value estimated by Margalit et al. (2018). As we shall show, this expression works for jets that break out before the engine has begun to turn off, so that the jet power can be considered a constant.

In this work, we also find a "late breakout" regime in which the jet breaks out after the engine is turned off, t_bo > T (something similar has been seen previously in long GRB scenarios; see Lazzati et al. 2012). In this case, essentially all of E_j has been delivered to the cocoon prior to breakout, which continues to collimate the jet. The jetted outflow has a characteristic terminal velocity v_sh that takes time to overtake the fastest-moving ejecta. Since the breakout can occur well after the engine turns off, this regime is analogous to a collimated explosion in a homologous outflow. We will calibrate the breakout time for this regime empirically in Section 4.

2.3. Thermal Energy in the Cocoon

Some fraction of the engine power will be converted into thermal energy of the cocoon as the jet propagates through the ejecta. This, in addition to thermalized energy from radioactive decay, can potentially help power the thermal electromagnetic counterpart to neutron star mergers. The engine likely outputs a primarily kinetic (or magnetic) flow that thermalizes in a shock when colliding with the ejecta. Imagine that some fraction η of the engine power is instantly thermalized in this shock. The thermal energy as a function of time is then given by the equation

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{Th}}}{{dt}}=\eta L(t)-{E}_{\mathrm{Th}}/t,\end{eqnarray} \tag{ 5 }$$

where the second term is adiabatic conversion to kinetic energy in the homologous phase (a radiation-dominated equation of state is assumed). Note that before the jet and cocoon themselves reach homology, the actual work done on the ejecta is not the homologous expression −E_Th/t, but rather p_cdV/dt (where p_c and dV/dt are the pressure and rate of change of volume of the cocoon, respectively); this in turn depends on how the jet head and cocoon propagate through the ejecta. In Equation (5) we are absorbing this physics into the factor η, which encapsulates the thermalization of the injected jet energy in the cocoon and the expansion of a cocoon-driven shock into the ejecta. We will show below that analytical considerations result in a value of η ≈ 0.8, close to the value we find in our numerical calculations.

Equation (5) can be straightforwardly solved:

$$\begin{eqnarray}&&{E}_{\mathrm{Th}}(t)=\displaystyle \frac{\eta }{t}{\int }_{0}^{t}t^{\prime} L(t^{\prime} ){dt}^{\prime} .\end{eqnarray} \tag{ 6 }$$

If the jet breaks out, however, the forward shock detaches and the engine no longer efficiently thermalizes its energy. To first approximation, the first term in Equation (5) then vanishes, and the thermal energy drops as t⁻¹:

$$\begin{eqnarray}&&{E}_{\mathrm{Th}}(t)={E}_{\mathrm{Th}}({t}_{\mathrm{bo}})\displaystyle \frac{{t}_{\mathrm{bo}}}{t}\ \ \ (t\gt {t}_{\mathrm{bo}}).\end{eqnarray} \tag{ 7 }$$

In this study, we consider an engine model with a smooth cutoff at late times:

$$\begin{eqnarray}&&L(t)=\displaystyle \frac{{E}_{j}/T}{{(1+t/T)}^{2}},\end{eqnarray} \tag{ 8 }$$

so that Equation (6) yields

$$\begin{eqnarray}{E}_{\mathrm{Th}}(t)=\left\{\begin{array}{ll}\eta {E}_{j}\tfrac{T}{t}\left[\mathrm{ln}(1+t/T)-\tfrac{t}{T+t}\right] & t\lt {t}_{\mathrm{bo}}\\ \eta {E}_{j}\tfrac{T}{t}\left[\mathrm{ln}(1+{t}_{\mathrm{bo}}/T)-\tfrac{{t}_{\mathrm{bo}}}{T+{t}_{\mathrm{bo}}}\right] & t\gt {t}_{\mathrm{bo}}\end{array}\right..\end{eqnarray} \tag{ 9 }$$

If the shock breaks out at sufficiently early times, then one can ignore the details of how quickly the engine shuts off. However, there is an additional complication, as the engine can still deposit thermal energy (at a diminished efficiency) onto the sides of the tunnel. In this case, the first term in Equation (5) does not strictly vanish, but is instead replaced with a diminished efficiency, η₂. This would lead to the solution

$$\begin{eqnarray}{E}_{\mathrm{Th}}(t)=\left\{\begin{array}{ll}\eta {E}_{0}(t) & t\lt {t}_{\mathrm{bo}}\\ \eta {E}_{0}({t}_{\mathrm{bo}})\tfrac{{t}_{\mathrm{bo}}}{t}+{\eta }_{2}{E}_{0}(t) & t\gt {t}_{\mathrm{bo}}\end{array}\right.,\end{eqnarray} \tag{ 10 }$$

where ${E}_{0}(t)={E}_{j}\tfrac{T}{t}\left[\mathrm{ln}(1+t/T)-\tfrac{t}{T+t}\right]$. The details of this formula are of course sensitive to our choice of L(t), i.e., precisely how the power decays with time. More generally, a good estimate for the late-time thermal energy is

$$\begin{eqnarray}{E}_{\mathrm{Th}}(t)\approx \left\{\begin{array}{ll}{\eta }_{2}{E}_{j}\tfrac{T}{t} & \mathrm{early}\ \mathrm{breakout}\\ \eta {E}_{j}\tfrac{T}{t} & \mathrm{late}\ \mathrm{breakout}\end{array}\right..\end{eqnarray} \tag{ 11 }$$

In this study, we find that the transition between early and late breakout occurs around E_j ∼ 30E_crit (where E_crit is the transition between failed and successful jets) and that η = 0.8, η₂ ≈ 0.02η for our models.

Chevalier & Fransson (1992) calculate the fraction of thermalized energy for a spherically injected engine in a homologous outflow:

$$\begin{eqnarray}&&{E}_{\mathrm{Th}}=\displaystyle \frac{5-k}{11-2k}{Lt}=\displaystyle \frac{1}{2}\eta {Lt}\end{eqnarray} \tag{ 12 }$$

for a pulsar wind injected into a density profile with ρ ∝ r^−k (their Equation (2.8)). This result is weakly dependent on k and leads to a thermalization efficiency of η ≈ 0.8 for most of the propagation of the jet (using k = 2.5). Our numerical results are consistent with this value (see the analytical curves in Figure 4, which assume η = 0.8).

We numerically study the interaction of a jet with a homologous outflow using the equations of 2D axisymmetric relativistic hydrodynamics:

$$\begin{eqnarray}&&{\partial }_{\mu }(\rho {u}^{\mu })={S}_{D}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\partial }_{\mu }(\rho {{hu}}^{\mu }{u}^{\nu }+{{Pg}}^{\mu \nu })={S}^{\nu },\end{eqnarray} \tag{ 14 }$$

where ρ is proper density, ρh = ρ + + P is enthalpy density, P is pressure, is the internal energy density, u^μ is the four-velocity, and g^μν is the Minkowski metric (the equations are expressed in units such that c = 1). The source terms S_D and S^ν will be used to model the injection of mass, momentum, and energy by the central engine on small scales. The equation of state is assumed to be radiation dominated:  = 3P.

The equations are integrated using the moving-mesh hydrodynamics code, JET (Duffell & MacFadyen 2011, 2013). The JET code's moving mesh makes it effectively Lagrangian, which allows for accurate resolution of supersonic and relativistic flows. The inner and outer boundaries are also moved during the flow's evolution, so that the dynamic range covered at any given time is kept modest. The spatial resolution is given by Δθ ∼ Δr/r ∼ 0.01. This value is not exact, as zones can expand and contract, but zones are refined or de-refined so as to maintain aspect ratios of order unity. Additionally, zones are more concentrated near the symmetry axis, so that resolution there is about twice as high (Δθ ∼ Δr/r ∼ 0.005), and in practice zones near shock fronts are compressed by a factor of approximately five so that Δr/r ∼ 0.001 in the most critical regions.

Initial conditions are modeled after the numerical output of a neutron star merger simulation (Hotokezaka et al. 2013; Nagakura et al. 2014). It is an outflow described by the following:

$$\begin{eqnarray}\rho (r,\theta ,0)=\left\{\begin{array}{ll}f(\tilde{r})+{\rho }_{\mathrm{atm}}(r) & \tilde{r}\lt {R}_{0}\\ {\rho }_{\mathrm{atm}}(r) & \tilde{r}\gt {R}_{0}\end{array}\right.\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}{v}_{r}(r,\theta ,0)=\left\{\begin{array}{ll}r/{t}_{0} & \tilde{r}\lt {R}_{0}\\ 0 & \tilde{r}\gt {R}_{0}\end{array}\right.\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&f(r)={\rho }_{c}\displaystyle \frac{{(1-r/{R}_{0})}^{3/4}}{1+{(r/{R}_{1})}^{2.5}}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\rho }_{\mathrm{atm}}(r)={\rho }_{\mathrm{trans}}{e}^{-r/{R}_{0}}+{\rho }_{\mathrm{ISM}}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\tilde{r}}^{2}={r}^{2}({w}^{4/3}\,{\cos }^{2}\theta +{w}^{-2/3}\,{\sin }^{2}\theta ),\end{eqnarray} \tag{ 19 }$$

where ${\rho }_{c}=97{M}_{0}/{R}_{0}^{3}$, R₁ = 0.06R₀, ${\rho }_{\mathrm{trans}}={10}^{-8}{M}_{0}/{R}_{0}^{3}$, ${\rho }_{\mathrm{ISM}}={10}^{-20}{M}_{0}/{R}_{0}^{3}$, and R₀ = v_maxt₀, with v_max = 0.3c. The engine is initiated at a very early time t₀ = 10⁻³T, where T is the engine duration. The ejecta model above corresponds to a mass-weighted velocity of the ejecta of 0.1c and kinetic energy-weighted velocity of the ejecta of 0.18c.

In interpreting our results, the mass M₀ of the cloud is assumed to have a fiducial value of 0.07 M_⊙, which gives a kinetic energy of 10⁵¹ erg. The dynamics of the interaction are, however, scale invariant, so one could equally well assume, say, M₀ = 0.01 M_⊙ and E_ej = 1.4 × 10⁵⁰ erg. In cgs units, the external density ρ_ISM ≈ 5 × 10⁻¹⁷ g cm⁻³ is much larger than typical interstellar medium densities of ∼1 proton per cm³, but this is irrelevant for the present study, as we are only concerned with the interaction between the jet and the ejecta, so this density should be considered an artificial atmosphere, with low enough density not to affect the dynamics. The parameter w = 1.3 is the aspect ratio of the oblate cloud. The choice w = 1.3 has an impact on whether or not the jet can break out of the ejecta (Duffell et al. 2015), and therefore this affects any dimensionless coefficients reported in this study, but we expect that the overall scalings are independent of the details of the ejecta model.

The engine power is prescribed by Equation (8), L(t) = (E_j/T)/(1+t/T)², which is a constant power for t ≪ T but decays as E_jT/t² at late times. The power-law tail is nonstandard in most numerical calculations of GRB jets; typically jet power is imposed for some given duration with an instantaneous shutoff. Nature is more likely to be well modeled as a power-law decay than to have a sharp cutoff. The exponent of −2 is largely arbitrary, but it could mimic late-time power input by a central engine, e.g., magnetar spin-down or fallback accretion.

The two dimensionless variables that are varied in this study are E_j/E_ej (the ratio of injected energy to the kinetic energy of the outflow) and θ_j, the opening angle of the injection (which may be different from the opening angle of the final jet). We have also performed additional runs varying v_max in the ejecta model (Section 3). This effectively varies the energy scale E_ej while keeping the mass fixed and confirms that the scaling with E_ej at fixed M_ejc² (or vice versa) is weak over the range of mean ejecta velocities of ∼0.1c–0.4 c expected in neutron star mergers.

Figure 1 shows the four classes of solutions that we find depending on the values of the key dimensionless numbers E_j/E_ej and θ_j. These are, from left to right, successful jet breakout from the homologous ejecta on a timescale ≲T, the engine duration, failed jets, late breakout on timescales ≳T, and failed jets in which the forward shock and cocoon driven into the ejecta by the jet nonetheless break out. In what follows we delineate the regimes in which these different solutions occur and some of their key properties.

Zoom In Zoom Out Reset image size

4.1. Success versus Failure

We argued in Section 2 that the ability of a jet to break out of the ejecta should depend on the jet energy relative to a critical energy E_crit, with

$$\begin{eqnarray}&&{E}_{\mathrm{crit}}=\kappa {E}_{\mathrm{ej}}{\theta }_{j}^{2},\end{eqnarray} \tag{ 20 }$$

where κ is some dimensionless constant. Figure 2 shows success versus failure for a wide range of models with variable energy and opening angle. A successful jet is defined here as one that emerges from the outflow collimated and relativistic. A successful shock breakout is defined to occur if any shock at all makes it to the edge of the ejecta. For all except the cyan points in Figure 2, successful shock breakout necessitated a successful jet. We find the ${\theta }_{j}^{2}$ scaling to be accurate for determining success of the jet, with the coefficient κ = 0.05 in Equation (20). The value of this coefficient should depend on the details of the initial density profile and the engine power as a function of time.

Zoom In Zoom Out Reset image size

One important caveat of these numerical calculations is that they have been performed in 2D, and some of these jets fail owing to "corking" (the jet is choked because a massive cork or plug is trapped ahead of the jet, which deflects the jet away from the symmetry axis; this might not occur in 3D; see Duffell & MacFadyen 2015 for a discussion of the cork and Zhang et al. 2004 for a comparison between 2D and 3D). These are labeled as such in Figure 2. Fortunately, a scaling for E_crit(θ_j) can be discerned despite these possibly artificially choked jets.

4.2. Breakout Time

Breakout times of successful jets are shown in Figure 3. Breakout times are consistent with the following expression:

$$\begin{eqnarray}{t}_{\mathrm{bo}}=\left\{\begin{array}{ll}{\alpha }_{1}T\tfrac{{E}_{\mathrm{crit}}}{{E}_{j}} & \mathrm{early}\ \mathrm{breakout}\\ {\alpha }_{2}\tfrac{T}{\sqrt{{E}_{j}/{E}_{\mathrm{crit}}}-1} & \mathrm{late}\ \mathrm{breakout}\end{array}\right.,\end{eqnarray} \tag{ 21 }$$

where the "early breakout" scaling was estimated in Section 2.2, but the "late breakout" scaling is entirely empirical. For E_j ≳ 30E_crit, breakout times follow the "early breakout" scaling with a coefficient α₁ ≃ 3, while for E_crit < E_j ≲ 30E_crit there is more scatter, but breakout times reasonably track the "late breakout" scaling, with α₂ ≃ 9.

Zoom In Zoom Out Reset image size

4.3. Thermal Energy

The total thermal energy in the computational domain is shown as a function of time in Figure 4. In the top panel, this is compared with the analytical form given by Equation (11). For failed jets (green curves), the thermal energy follows the predicted formula precisely, with a thermalization efficiency η = 0.8. Successful jets with early breakout (pink curves) have a sudden transition at breakout, where the thermal energy drops and thermalization is no longer efficient because the jet vents through the channel carved in the ejecta. Jets with an early breakout exhibit a second increase in thermal energy due to the thermalization of the jet onto the sides of the tunnel. This comes with a reduced efficiency η₂ ≈ 0.02η, though this efficiency could be higher if the jets precess or the jet opening angle varies in time owing to variability in the accretion disk driving the jet. After the engine shuts off, i.e., for t ≳ T, the thermal energy evolution transitions to a t⁻¹ decline consistent with adiabatic conversion from thermal to kinetic energy.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 4 shows the thermal energy as a function of time in units of the ejecta kinetic energy rather than the input jet energy. This highlights that it is difficult to put a significant amount of thermal energy into the cocoon/ejecta via a jet even when the total energy of the jet approaches the initial kinetic energy of the ejecta (redder curves). This is primarily because high jet energies result in an early breakout, cutting off the thermalization and leading to most of the jet energy escaping the ejecta.

Figure 5 shows the peak thermal energy in each of our models. This peak is maximized for failed jets and late breakouts with large opening angles; for such cases, the peak thermal energy is E_peak ≈ 0.17E_j, where the factor 0.17 can be determined from the peak of Equation (9). The exact value of this coefficient at the factor of few level will depend on the precise form of L(t) and the ejecta structure. Figure 5 shows that only for very large opening angles and jet energies approaching E_ej ∼ 10⁵¹ erg can the jet deposit an appreciable thermal energy in the ejecta. More typically, E_peak ≪ E_ej.

Zoom In Zoom Out Reset image size

The energetics of short GRBs are less constrained than those of long-duration bursts, because the opening angles of short bursts are usually more uncertain. However, current estimates from afterglow observations are that typical energies are in the range of 10⁴⁹–10⁵⁰ erg, with opening angles of θ_j ∼ 0.05–0.4 (Fong et al. 2015). An independent (but model-dependent) constraint on the opening angle is provided by the relative rate of short GRBs and neutron star mergers, which favors beaming factors of ∼100 and thus θ_j ∼ 0.1.

Note that opening angles measured from afterglow observations utilize the "jet break" and may be different from the opening angle injected by the central engine, which is how θ_j is defined in the present study. Similarly, the duration of the burst might not be the same as the duration of the engine T, though we assume that a typical value of T should be in the typical observed duration range of 0.1–1 s. All of the calculations in this work assumed that the engine is initiated at early times, ≪T and ≪t_bo. This is not guaranteed to be true, and future work should investigate the consequences of relaxing this assumption. However, this assumption is physically motivated, as an accretion disk should form on millisecond timescales. Even if a hypermassive neutron star survives for ∼0.1–1 s before collapsing to a black hole, it is unlikely that much mass will remain at later times to accrete onto the black hole given the short viscous times in the disk around the hypermassive neutron star (Metzger & Fernández 2014 demonstrate this explicitly in idealized axisymmetric hydrodynamic simulations). It is possible, however, that early time outflows and jets are baryon dominated because of the strong neutrino-driven wind from the neutron star (e.g., Metzger et al. 2018). In this case, a highly relativistic outflow may emerge only at later times, effectively producing a delay in the onset of the relativistic engine even though most of the accretion energy is in fact released promptly.

For jet energies and opening angles at the upper limits suggested by observations of short GRBs (E_j ∼ 10⁵⁰ erg ∼ 0.1E_ej and θ_j ∼ 0.2–0.4), Figures 4 and 5 show that the peak thermal energy due to the jet propagating through the ejecta is ∼10⁻³E_ej−10⁻²E_ej ∼ 10⁴⁸−10⁴⁹ erg (where the latter numbers assume our fiducial E_ej = 10⁵¹ erg). Moreover, this energy is reached on a timescale comparable to the engine duration T ∼ 0.1–1 s. For comparison, the r-process produces a thermal energy of ≃1–3 MeV per nucleon on a timescale of about 1 s (Metzger et al. 2010a), with the heating then decaying as a power law at later times. For our fiducial ejecta mass of 0.07 M_⊙ this corresponds to a thermal energy from the r-process of (1–3) × 10⁵⁰ erg ≳ 0.1E_ej on a timescale of 1 s.¹ Figure 5 thus shows that r-process heating of the ejecta in neutron star mergers on second timescales dominates over jet heating in the entire parameter space of numerical models studied in this paper (usually by a large factor). In addition, the more energetic the engine, the earlier the jet breaks out, thermalizing a smaller fraction of its energy, meaning that these statements are quite robust even if one considers highly energetic engines. The r-process heating will thus set the overall luminosity of any associated thermally powered transient. Note that the dominance of r-process heating over jet heating is usually discussed in the context of day-timescale kilonova emission. Here we have shown that this is true in most cases even on second timescales. We do find, however, that for a small amount of mass along the pole the energy per unit mass from shock heating exceeds that produced by radioactive heating in the first few seconds. The angular dependence of this deposited energy may have observational consequences, but calculating these consequences explicitly would require full multidimensional radiation transport. The implications of this for kilonova light curves will be explored in future work.

For almost all of our parameter choices, we find that shock breakout of the cocoon only occurs if the jet itself is able to escape the ejecta. This is in contrast with some of the shock breakout scenarios presented for GRB 170817A (e.g., Gottlieb et al. 2018), which assume cocoon breakout after a choked/failed jet. The regime of shock breakout with a choked jet is only realized for large jet energies and large jet opening angles (see the right panel of Figure 1 and the cyan points in Figure 2), which are indeed the parameters used in Gottlieb et al. (2018). By contrast, for lower energies and opening angles, which appear to be somewhat more typical for short GRBs, the jet and cocoon both likely break out (Figure 2).

We find that the timescale for the jet to breakout of the ejecta scales inversely with the jet power for E_j ≫ E_crit (This is in part a consequence of our assumption of negligible delay between mass ejection and engine initiation). For energies just above E_crit, however, the difference between the shock velocity due to the jet and the homologous ejecta velocity is small, resulting in a late breakout on timescales that can be significantly longer than the engine duration T. In particular, reasonable engine parameters consistent with short-duration GRBs can be consistent with the 1.7 s delay between gravitational waves and gamma rays in GW170817. For example, a delay of ∼3–10 T occurs in our late breakout regime with E_j ∼ 10E_crit (Figure 3). Taking θ_j ∼ 0.1, this corresponds to E_j ∼ 10⁻²E_ej ∼ 10⁴⁹ erg, well within reasonable short-GRB parameters. The hypothesis that the delay between gamma rays and gravitational waves in GW170817 is caused by the long breakout time can be tested by future observations. For a reasonable spread in engine parameters, it is likely that future electromagnetic (EM) counterparts to neutron star mergers will be separated into two classes: bursts with a substantial delay in gamma rays, and those without. The bursts without delay would have larger ratios of jet energy to ejecta energy. Mergers producing very little dynamical ejecta and thus fainter kilonovae would have very little delay.

The calculations presented in this paper are scale-free and can be rescaled to different assumed dynamical ejecta and jet properties. We have used a kinetic energy for the quasi-spherical dynamical ejecta of E_ej = 10⁵¹ erg to anchor our discussion. This is consistent with observations of the optical–IR counterpart to GW170817 (e.g., Villar et al. 2017). Quasi-spherical dynamical ejecta with this energy can be produced in binary neutron star mergers for sufficiently small neutron star radii (e.g., Radice et al. 2016). If, however, the quasi-spherical ejecta energy is lower—which is consistent with GW170817 observations if the blue kilonova component is produced by, e.g., accretion disks or magnetar winds—a jet of a given energy and opening angle will break out more easily and will deposit less thermal energy into the cocoon as a result (Figures 3 and 5).

By using short-GRB observations to motivate engine parameters for gravitational-wave-selected events, we are assuming that at least some short GRBs are produced by binary neutron star mergers. We believe that this empirical calibration is reasonable even if short GRBs are instead due to neutron star–black hole mergers, because the mass and angular momentum content of the accretion disks are not that different. By contrast, if short GRBs are produced by an entirely different source (e.g., accretion-induced collapse of white dwarfs to neutron stars), then calibrating engine parameters of gravitational-wave-detected events to short-GRB observations may be inappropriate. It is, of course, also possible that binary neutron star mergers do produce most short GRBs, but that there is a large dispersion in short-GRB engine properties and that GW170817 just happened to be somewhat unusual relative to the bulk of the population (e.g., some short GRBs do not show any evidence for a jet break, consistent with a large opening angle of θ_j ∼ 0.4; see, e.g., GRB 050724A in Fong et al. 2015). Multimessenger observations in the coming decade will likely empirically settle these questions.

We thank Brian Metzger, Juna Kollmeier, and Tony Piro for stimulating conversations that helped to motivate much of this work. We thank Wen-fai Fong, Ore Gottlieb, Andrew MacFadyen, and Rosalba Perna for useful conversations. This work was supported, in part, by the Theoretical Astrophysics Center at UC Berkeley. This research was also funded by the Gordon and Betty Moore Foundation through grant GBMF5076. E.Q. was supported in part by a Simons Investigator award from the Simons Foundation. H.K. is supported by a DOE Computational Science Graduate Fellowship under grant no. DE-FG02-97ER25308. Numerical computations utilized the Savio computational cluster resource provided by the Berkeley Research Computing program at the University of California, Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor of Research, and Office of the CIO).

We would like to thank the anonymous referee for the very clear, helpful, and thorough review.