JET COLLIMATION IN THE EJECTA OF DOUBLE NEUTRON STAR MERGERS: A NEW CANONICAL PICTURE OF SHORT GAMMA-RAY BURSTS

Recent afterglow observations of short gamma-ray bursts (SGRBs) have provided various information about their environments which can be interpreted as circumstantial evidence linking SGRBs with mergers of compact binaries such as double neutron stars (NS–NS; Paczynski 1986; Goodman 1986; Eichler et al. 1989) and black hole–neutron star (BH–NS; see Berger 2013 for a latest review). However, the compact binary merger scenario is challenged by the detection of jet breaks in the afterglow of some SGRBs and the deduced small jet opening angle of ≲ 10° (Soderberg et al. 2006; Burrows et al. 2006; Nicuesa Guelbenzu et al. 2011; Fong et al. 2012, 2013). The formation of such a collimated jet in compact binary mergers has not been clarified yet (see, e.g., Aloy et al. 2005; Bucciantini et al. 2012).

One of the most interesting features in the latest numerical relativity simulations (Hotokezaka et al. 2013a) is that NS–NS mergers in general are accompanied by a substantial amount of dynamical mass ejection. Interestingly, the excess in the near-IR band observed by the Hubble Space Telescope in Swift SGRB 130603B (Tanvir et al. 2013; Berger et al. 2013) can be explained by the kilonova/macronova model (Li & Paczyński 1998; Metzger et al. 2010; Kasen et al. 2013; Barnes & Kasen 2013; Grossman et al. 2013; Tanaka & Hotokezaka 2013) provided that a large amount of mass ≳ 2 × 10⁻²M_☉ is ejected in the NS–NS merger and it is powered by the radioactivity of r-process nuclei (Hotokezaka et al. 2013c; Tanvir et al. 2013; Piran et al. 2014). Such massive ejecta will have a large impact on the dynamics of the jet and the observed collimation could be naturally explained by their interactions.

In this Letter, we numerically investigate the jet propagation in the material ejected by double NS mergers based on a scenario indicated both by our latest numerical relativity simulations and observations of SGRB 130603B. The scenario is summarized as follows (see Figure 1).

• 1.  
  According to the latest numerical relativity simulations adopting equations of state (EOSs) that are compatible with the recent discovery of massive NSs with M ∼ 2 M_☉ (Demorest et al. 2010; Antoniadis et al. 2013), a hypermassive neutron star (HMNS) is the canonical outcome after an NS–NS merger for the typical binary mass (2.6–2.8 M_☉; Sekiguchi et al. 2011; Hotokezaka et al. 2013a; Bauswein et al. 2013).
• 2.  
  During and after the merger a large amount of mass O(0.01 M_☉) is ejected (phase (II)). The size of this ejecta is required to explain the kilonova candidate associated with SGRB 130603B. According to our numerical relativity simulations (Hotokezaka et al. 2013a), the morphology of the ejecta is quasi-spherical in the case of HMNS formation. In particular, the regions along the rotational axis are contaminated significantly by the mass ejection.
• 3.  
  Such a large amount of mass can be ejected only if the EOS of the NS matter is relatively soft (Hotokezaka et al. 2013a, 2013b; Bauswein et al. 2013). In this case, the massive NS formed after the merger is expected to collapse to a BH within several tens of milliseconds (phase (III)), forming a massive torus around it.
• 4.  
  After the formation of the BH–torus system, a jet that propagates through the expanding merger ejecta is launched (phase (IV)). An SGRB will be produced only if the jet successfully breaks out of the ejecta.

Zoom In Zoom Out Reset image size

Note that our scenario is different from that explored by previous studies (Aloy et al. 2005) based on Newtonian studies (Rosswog et al. 1999), in which the mass ejection is not isotropic but is concentrated along the orbital plane. In this case, there will be little interaction between the jet and ejecta, and no collimation by the ejecta is expected. Indeed, Aloy et al. (2005) found no strong collimation by the disk wind (see also Levinson & Eichler 2000), since their simulations were carried out in rather dilute ejecta (<10⁻³M_☉).

After studying the dynamics of the jet in the presence of the expanding ejecta, we discuss the canonical model used to explain a particular event, SGRB 130603B. With the observationally consistent parameter set, we show that relativistic jets successfully break out of the dynamical ejecta and travel with the required collimation angle.

For constructing ejecta profile models, the results from numerical relativity are employed as reference. We first analyze results in Hotokezaka et al. (2013a), and then fit the ejecta profile along the pole using the following formulae:

$$\begin{equation} \rho (t_{\rm {i}},r) = \rho _{0}(t_{\rm {i}}) \left (\frac{r}{r_{0}} \right)^{-n}, \end{equation} \tag{ 1 }$$

$$\begin{equation} r_{\rm {max}} (t_{\rm {i}}) = v_{\rm {max}} (t_{\rm {i}}-t_{0}) + r_{\rm {max0}}, \end{equation} \tag{ 2 }$$

$$\begin{equation} v(t_{\rm {i}},r) = v_{\rm {max}} \left (\frac{r}{r_{\rm {max}}} \right). \end{equation} \tag{ 3 }$$

In the above expressions, t_i, r, ρ, and v denote the onset time of jet injection (measured from the merger time), radius, rest-mass density, and velocity of ejecta, respectively. The other variables, n, v_max, r₀, and t₀, are the fitting parameters. The power-law index of the density distribution (n) is more or less dependent on the dynamics of merger, which is in the range 3 < n < 4. We choose the middle of this range n = 3.5 in this study. v_max denotes the velocity at the dynamical ejecta front (we set v_max = 0.4c). t₀ denotes the snapshot time at which we refer to the result of the numerical relativity merger simulations. We set t₀ = 10 ms, since the morphology of the ejecta has been determined by that time and the outer ejecta continues to be in the homologous expansion phase (Rosswog et al. 2014). The location of the forward shock wave at t₀ is denoted by r_max0, which is set to be r_max0 = 1.3 × 10⁸ cm. The rest-mass density ρ₀(t_i) can be expressed as a function of ejecta mass (M_ej) as

$$\begin{equation} \rho _{0} (t_{\rm {i}}) = \frac{(n-3) M_{\rm {ej}}}{ 4 \pi r_{0}^3 } \left \lbrace \left (\frac{r_{\rm {esc}}}{r_{0}} \right)^{3-n} - \left (\frac{r_{\rm {max}}}{r_{0}} \right)^{3-n} \right \rbrace ^{-1}, \end{equation} \tag{ 4 }$$

where

$$\begin{equation} r_{\rm {esc}} = \left (\frac{2GM_{\rm {c}} r_{\rm {max}}^2 }{v_{\rm {max}}^{2}} \right)^{\frac{1}{3}}, \end{equation} \tag{ 5 }$$

M_c denotes the central remnant mass, which is chosen to be M_c ≡ 2.7M_☉, and r_esc denotes the escape radius, which is defined as $v(t_{\rm {i}},r_{\rm {esc}}) \equiv \sqrt{2GM_{\rm {c}}/r_{\rm {esc}}}$. The pressure of the ejecta is set to p = K_efρ^4/3 with $K_{\rm {ef}}=2.6 \times 10^{15}\, \rm {g}^{-1/3} \,\rm {cm}^3\, \rm {s}^{-2}$, which is cold enough not to affect the jet and ejecta dynamics.

From these formulae, we determine the ejecta profile as a function of t_i and M_ej. We first examine the case of M_ej = 10⁻²M_☉ (see Table 1), which is the approximate value of the required mass to explain the kilonova associated with SGRB 130603B (Hotokezaka et al. 2013c); we then study the dependence on M_ej (M–M3, M–M2-2, M–M1). t_i corresponds to the time of jet injection, which is supposed to be the operation timing of the central engine. For this there are no observational constraints. We set t_i = 50 ms to be the reference value, since our numerical relativity simulations predict that the lifetime of HMNS is likely to be several tens of milliseconds to explain the large mass of the ejecta M_ej ∼ 10⁻²M_☉ as well as the large mass of the torus surrounding the black hole (BH). For comparison, we study the t_i = 500 ms case for one model (M–ti500; see Table 1).

Using the ejecta profile obtained above for the initial conditions, we perform axisymmetric simulations of the jet propagation by employing a relativistic hydrodynamical code (Nagakura et al. 2011, 2012; Nagakura 2013). We assume that the central engine successfully operates in the vicinity of the compact remnant, and that the jet is injected with constant power from the innermost computational boundary. In these simulations, we focus only on exploring the interaction between the jet and ejecta. Therefore, the computational domain covers r_esc to 2 × 10¹⁰ cm. The canonical jet power is set to be L_j = 2 × 10⁵⁰ erg s⁻¹ for all models, which is comparable with the average jet power of SGRB 130603B (see Fong et al. (2013) for the collimation-correlated jet energy and duration of the prompt emission). We also prepare the model M–L4 for which L_j = 4 × 10⁵⁰ erg s⁻¹ to study the dependence of the jet luminosity. Throughout our simulations, we use the gamma-law EOS with γ = 4/3. The initial Lorentz factor (Γ_ini) and specific enthalpy (h_ini) are set to Γ_ini = 5 and h_ini = 20, which result in the terminal Lorentz factor Γ_term = 100. The initial jet opening angle (θ₀) is also not well constrained by observations, and hence we set θ₀ = 15° as the reference value with θ₀ = 30°, 45° to study the dependence on θ₀ (M–th30, M–th45). Note that θ₀ = 15° is larger than the opening angle of 1/Γ_ini = 1/5 ∼ 12°, so that the initial thermal expansion of the jet would not be significant (see, e.g., Mizuta & Ioka 2013). Simulations are carried out until the shock reaches the outer boundary or time is 1 s after the jet injection. Our models are summarized in Table 1.

Starting from the initial moment of jet injection at the chosen post-merger time, the jet begins to burrow through the homologously expanding ejecta with mildly relativistic velocity. In the two left panels of Figure 2, we display the density contour maps for M–ref at the time of jet breakout and the end of our simulation. At a short distance from the inner boundary, the jet structure changes from conical to cylindrical due to confinement by the dense ejecta. The small cross section of the jet head allows the shocked jet matter to escape sideways and generates a hot cocoon around the jet. Even though the density gradually decreases with the radius, the surrounding cocoon keeps the jet confined near the pole, and eventually the jet head successfully breaks out of the edge of the ejecta. The overall properties of the interaction between the ejecta and jet are very similar to those in the context of the collapsar model (Nagakura et al. 2011; Mizuta & Ioka 2013). A remarkable difference between the jet propagation in the NS–NS ejecta and the stellar mantle is that the background fluid is no longer stationary and expands with time. The jet head chases the ejecta edge from behind, and needs to catch up with it for the relativistic breakout; otherwise it would become non-relativistic ejecta and will never produce SGRBs (see below).

Zoom In Zoom Out Reset image size

For the less massive ejecta case (M–M3), the jet experiences less confinement and propagates faster than M–ref (see right panels in Figure 2). Even so, the hot cocoon is formed by the jet–ejecta interaction and works to weakly confine the jet. In order to analyze the cocoon confinement and its degree, we use the dimensionless jet luminosity parameter ($\tilde{L} \equiv \rho _{\rm {j}} h_{\rm {j}} \Gamma _{\rm {j}} / \rho _{\rm {a}}$, where ρ_a denotes the ambient density above the jet head) following the study by Bromberg et al. (2011). By employing Equations (1)–(5) and imposing the condition r_esc ≪ r_max, $\tilde{L}$ can be roughly estimated as

$$\begin{eqnarray} \tilde{L} &\sim& 10^{-3} \left (\frac{L_{\rm {j}}}{2 \times 10^{50} {\rm {erg \,s^{-1}}}} \right) \left (\frac{M_{\rm {ej}}}{10^{-2}\, M_{\odot }} \right)^{-1} \nonumber \\ &&\times \left(\frac{\theta _{0}}{15^{\circ }} \right) ^{-2} \left (\frac{t_{\rm {i}}}{50\,{\rm {ms}}} \right)^{\frac{2}{3}} \left (\frac{\epsilon _{r}}{1} \right)^{n} \left (\frac{\epsilon _{t}}{1} \right)^{3-n}, \end{eqnarray} \tag{ 6 }$$

where

$$\begin{equation} \epsilon _{r} \equiv r_{\rm {j}}/r_{\rm {esc}}, \end{equation} \tag{ 7 }$$

$$\begin{equation} \epsilon _{t} \equiv t/t_{\rm {i}}, \end{equation} \tag{ 8 }$$

and r_j and t denote the radius of the jet head and the time after the merger, respectively. According to Bromberg et al. (2011)⁵, the condition of cocoon confinement is $\tilde{L} \lesssim \theta _{0}^{-4/3} \sim 6 (\theta _{0}/15^{\circ })^{-4/3}$. In the vicinity of r_esc (_r ∼ 1), all models (including M–M3) satisfy the confinement condition, which indicates that the jet undergoes a collimation at least once. The cocoon pressure decreases with time because the density of ejecta has a steep radial gradient (n ∼ 3.5). Despite the weakening cocoon pressure, the opening angle of the jet becomes smaller than the initial one. In order to analyze the degree of collimation more precisely, we define the average jet opening angle as

$$\begin{equation} \theta _{\rm {ave}}(t) \equiv \frac{ \int _{r_{\rm {esc}}}^{R_{\rm {jh}}} \theta _{\rm {op}}(t,r) dr }{ R_{\rm {jh}}(t) - r_{\rm {esc}} }, \end{equation} \tag{ 9 }$$

where R_jh denotes the radius of the jet head. The jet opening angle at each radius (θ_op) is defined as the angle of the relativistic components, for which hΓ > 10. Note that if we instead employ the criterion hΓ > 100, we would obtain the incorrectly small θ_op, due to baryon pollution from numerical diffusion. Figure 3 shows the evolution of θ_ave for each model. Indeed, θ_ave is always less than θ₀, which is clear evidence of jet collimation. We also find that θ_ave after the breakout is larger than ∼θ₀/5, which is different from the results in the collapsar case (Mizuta & Ioka 2013). This may be attributed to the fact that the ejecta is not stationary with respect to the stellar mantle, and the density gradient of the ejecta is steeper than in the case of the stellar mantle.

Zoom In Zoom Out Reset image size

The initial jet opening angle is also important for the dynamics of jet propagation. In reality, it would be determined in the vicinity of the HMNS or BH from the interaction between the jet and the hot accretion disk (Aloy et al. 2005), or from pinching by magnetic fields (McKinney 2006). One of the important consequences of this study is that all models succeed in breaking out by the end of our simulation except for M–th45 (θ₀ = 45°). For the failed breakout model (M–th45), the shocked jet and ejecta cannot go sideways into the cocoon because of the large cross section of the jet and eventually expands quasi-spherically. This fact gives an interesting prediction that there may be some population of choked (failed) SGRBs or new types of events with low luminosity, which could be potential candidates for high energy neutrinos (Mészáros & Waxman 2001; Razzaque et al. 2004; Ando & Beacom 2005; Horiuchi & Ando 2008; Murase & Ioka 2013; Osorio Oliveros et al. 2013). The rate of these events is uncertain, since it depends on the jet luminosity, opening angle, ejecta mass, and the operation timing of the central engine. We also find that the delayed central engine activity tends to result in failed SGRBs or low-luminosity events since the ejecta head has already traveled farther away from the merger remnant (see r_max of M–ti500 in Table 1).

We here discuss the canonical model for SGRB 130603B based on the results of our simulations. According to de Ugarte Postigo et al. (2013) and Fong et al. (2013), SGRB 130603B has a well-collimated jet (its opening angle is ∼4°–8°) with a prompt duration ΔT₉₀ ∼ 200 ms.

Here we focus on the two main properties of the jet: its breakout radius and opening angle. The breakout radius r_b is defined as the radius where the jet head reaches the edge of the ejecta. Broadly speaking, the spatial length of the jet (Δl_j) in SGRB 130603B is ΔT₉₀ × c ∼ 6 × 10⁹ cm. We regard r_b ≲ Δl_j to be a preferred condition for the generation of SGRBs. In this case, the central engine must be active longer than the jet breakout time t_b, so that the late parts of the jet could reach the emission region without dissipating much energy to the cocoon. The duration of the central engine can be estimated as Δt_ce ∼ t_b + (Δl_j − r_b)/c, which is ∼300 ms for M–ref (see Bromberg et al. (2012) for a comparison with long gamma-ary bursts, GRBs). By this criterion, M–th30, M–ti500 and M–M1 are discarded as candidates for SGRB 130603B.

The second property we focus on is the jet opening angle and its evolution. As shown in the previous section, the jet undergoes confinement by the ejecta and breaks out with a smaller opening angle than the initial one. The model M–M3 does not satisfy observational constraints for SGRB 130603B, because the opening angle that it reaches is too large (see Figure 3). Therefore, M–M3 may not be a good model for SGRB 130603B. Note that since θ_op includes the jet component inside the ejecta, it is not exactly equal to the observed opening angle. We check the average opening angle of the jet outside the ejecta, and it is not very different from θ_op.

According to these criteria, M–ref, M–L4, and M–M2-2 are favored candidates for SGRB 130603B. Note that if the intrinsic jet luminosity is much larger than L_j ∼ 10⁵⁰ erg s⁻¹, there is a possibility of GRBs production even for ∼0.1 M_☉ ejecta mass. Note also that if the initial jet opening angle is sufficiently small, it may not require the cocoon confinement to explain the observed small jet opening angle. However, that would become demanding with the central engine and we have yet to know the generation mechanism of such a well-collimated jet.

In this Letter, we investigate the jet propagation in the dynamical ejecta after the NS–NS merger. Similar to the collapsar model, the interaction between the jet and the merger ejecta generates the hot cocoon and the jet undergoes collimation at least by the deepest and densest layers of the ejecta, which is qualitatively consistent with the criterion $\tilde{L} \lesssim \theta _{0}^{-4/3}$. Importantly, models, except for those with quite large initial opening angles (θ₀ = 45), succeed in the breakout with an opening angle smaller than the initial one. We also, for the first time, show the possibility of the existence of some population of choked SGRBs or new types of low-luminosity events.

Using only the duration of the prompt emission, the jet opening angle, and ejecta mass, we present the canonical model for SGRB 130603B. Under the assumption of spherically symmetric ejecta, M–M2-2 model satisfies all observational constraints. In reality, however, the ejecta profile is not exactly spherically symmetric, and the mass contained in its equatorial region tends to be larger. According to this, the ejecta mass in the realistic system would be several factors larger than our spherical models. Therefore, M–ref and M–L4 could also be candidates for SGRB 130603B (Hotokezaka et al. 2013c; Tanvir et al. 2013; Piran et al. 2014).

The results of this study and of Hotokezaka et al. (2013c) suggest that the EOS of NS may be soft among several models of the EOS with maximum mass >2 M_☉ if the central engine of this SGRB is an NS–NS merger. The required condition for the central engine is that the jet should be collimated ≲ 15° before reaching the ejecta, and its lifetime should be ∼300 ms with L_j ≳ 2 × 10⁵⁰ erg s⁻¹ as the average jet power, and the time lag between merger and jet launching should not be much longer than several tens of milliseconds.

As discussed in this Letter, the cocoon confinement changes the conventional picture of jet propagation for the production of SGRBs, and reinforces the scenario of an NS–NS binary merger for SGRBs. In BH–NS mergers, the morphology of the dynamical ejecta is non-spherical, i.e., the ejecta is concentrated on the equatorial plane (see Kyutoku et al. 2013), so the jet never undergoes strong collimation unless neutrinos or magnetic-driven winds from the accretion disk provide enough baryons in the polar region.

We thank Yudai Suwa, Kenta Kiuchi, Kazumi Kahiyama, and Takashi Nakamura for fruitful discussions. This work was supported by Grant-in-Aid for the Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan (22244030, 23740160, 24000004, 24103006, 24244028, 24244036, 24740165, 25103512) and HPCI Strategic Program of Japanese MEXT. The work of K.H. is supported by a JSPS fellowship grant No. 24-1772.