Off-axis Prompt X-Ray Transients from the Cocoon of Short Gamma-Ray Bursts

The discovery of gravitational waves (GWs) from mergers of binary black holes (BHs) has opened a new window of study of the universe (Abbott et al. 2016a, 2016b). Each event has been accompanied by a massive search for electromagnetic (EM) counterparts, despite the lack of a general consensus for a production mechanism of such counterparts to BH–BH mergers (Connaughton et al. 2016; Loeb 2016; Perna et al. 2016; Zhang 2016; de Mink & King 2017). The identification of an EM counterpart would possibly lead to the identification of the host galaxy and to a redshift measurement. Moreover, EM signals carry complementary information that would allow for better constraints on the source properties.

Unlike BH–BH mergers, the coalescence of two neutron stars (NSs), or an NS and a BH, is expected to be accompanied by EM radiation: these events are believed to be the progenitors of short gamma-ray bursts (SGRBs; Eichler et al. 1989; Nakar 2007; Berger 2014). However, despite decades of indirect evidence pointing to this association, only a detection of simultaneous GWs and gamma-ray radiation would constitute a smoking gun that would confirm the association and finally solve the long-standing mystery of the origin of SGRBs. Additionally, the EM counterpart would allow for a better localization of the merger sites and hence help to constrain the evolutionary scenario which led to the binary formation. By combining the parameters inferred from the GW radiation (i.e., masses) with the energetics inferred from the SGRB EM emission, further information can be gained on the mass of the ejecta (Giacomazzo et al. 2013).

The question of the simultaneous observability of GWs and EM radiation is hence of paramount importance. However, while GWs are only moderately anisotropic, the γ-ray emission is likely produced within a collimated, relativistic outflow, which reduces the probability of seeing it in association with a GW event.⁸ The longer wavelength emission, if produced by the same relativistic jet as in the standard afterglow scenario, is also expected to be collimated, at least until the jet has slowed down to trans-relativistic speeds, weeks to months after the event (Rossi et al. 2002). Observations of SGRBs so far have estimated the average typical opening angle of the jet to be $\sim 16^\circ$ or less (Fong et al. 2015; Ghirlanda et al. 2016), making the probability of seeing a jet on axis (and hence a "standard" SGRB) less than 10%.

Given the above, additional sources of EM emission from NS–NS mergers become especially relevant. In the optical/near infrared band, an important contribution can be provided by the kilo/macronova, a transient phenomenon triggered by the radioactive decay of r-process nuclei in the neutron-rich material ejected during the NS–NS merger (Li & Paczyński 1998; Metzger et al. 2010; Metzger & Berger 2012; Kasen et al. 2013, 2015; Kawaguchi et al. 2016). At all other wavelengths, the best prospects for detection come from "side" emission (SE) from the jet.⁹ Even though it is significantly weaker than the on-axis emission, SE is potentially very important when in association with a GW event, as the distance to which NS–NS events can be detected with advanced LIGO is only about 200 Mpc, even after reaching design sensitivity (it is 65–115 Mpc for advanced Virgo; currently the limit is 80–120 Mpc for advanced LIGO and 20–60 Mpc for advanced Virgo; Abbott et al. 2016c).¹⁰ Therefore, the fainter SE significantly increases the chances of detectability for off-axis events.

A study of SE from the jet becomes especially timely during the GW era. Observer lines of sight at large angles are much more likely in association with a GW event from a compact binary merger than are on-axis observations of SGRBs. The probability density function of LIGO/Virgo detections with respect to the off-axis angle peaks at about 30° (Schutz 2011). The off-axis emission from long GRB jets has been widely studied and has taken various names in the literature, such as structured jet (e.g., Rossi et al. 2002), off-axis emission (e.g., Granot et al. 2005; Salafia et al. 2016), sheath (Kathirgamaraju et al. 2017), and cocoon/wide-angle (e.g., Lazzati et al. 2017). Some theories explain the phenomenon of X-ray flashes (an X-ray-dominated sub-class of long GRBs) as the off-axis emission from collapsars (Gotthelf et al. 1996; Yamazaki et al. 2002; Fynbo et al. 2004; Lamb et al. 2005; Sakamoto et al. 2005; Guidorzi et al. 2009); alternative theories, though, exist (see, e.g., Ciolfi 2016).

In Lazzati et al. (2017; hereafter L17) we presented calculations of the wide-angle emission of SGRBs from compact binary merger progenitors. This emission arises from a hot cocoon, which forms if the jet initially propagates through a baryon-contaminated region surrounding the merger site (Ramirez-Ruiz et al. 2002; Nakar & Piran 2017). The semi-analytical calculations in L17 were made under a number of approximations, such as the assumption of an isotropic cocoon. However, simulations by Gottlieb et al. (2017) have showed that this is not the case. Additionally, a specific, ad-hoc value of the Lorentz factor, ${\rm{\Gamma }}=10$, has been adopted.

Here we perform a full simulation of the formation of the cocoon and its interaction with the jet and unshocked ambient medium, using a nested 3D/2D approach. Simulations of the dynamical effect of the merger ejecta on a relativistic jet were pioneered by Nagakura et al. (2014) and Murguia-Berthier et al. (2014, 2017). We predict the distribution of the isotropic equivalent energy as a function of the off-axis angle. We then calculate the profile of the peak photon energy of the photospheric cocoon emission, as well as the profile of the pulse duration and time delay of the pulse from the jet launching time as a function of the off-axis angle. Our numerical methods are described in Section 2, and the results are described in Section 3. We summarize and discuss our results in Section 4.

The simulation of an SGRB jet was performed with the adaptive mesh refinement relativistic hydrodynamic (AMR–RHD) code FLASH (Fryxell et al. 2000) as modified in Morsony et al. (2007). The SGRB jet was simulated as an inflow boundary condition with similar properties as the jet in L17: ${L}_{j}={10}^{50}$ erg s⁻¹, ${\theta }_{j}=16^\circ$, ${t}_{\mathrm{eng}}=1\,{\rm{s}}$. However, it was injected with ${{\rm{\Gamma }}}_{\infty }=300$ and already mildly relativistic (${{\rm{\Gamma }}}_{0}=5$) at the inner boundary located at ${r}_{0}={10}^{8}\,\mathrm{cm}$ from the merger site. The local merger ejecta are approximated as an exponentially cutoff wind, with the density profile

$$\begin{eqnarray}&&n(r)={n}_{0}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{-2}\,{e}^{-\tfrac{r}{{r}_{0}}}\end{eqnarray} \tag{ 1 }$$

where, again, ${r}_{0}={10}^{8}\,\mathrm{cm}$, and ${n}_{0}={10}^{6}$ cm⁻³, for a total ejecta mass of 0.006 solar masses, 0.002 of which are inside the simulation domain (the remainder are located at $r\lt {r}_{0}$). Such masses and densities are comparable to the results of numerical simulations of binary NS and NS–BH mergers (Kiuchi et al. 2015; Endrizzi et al. 2016; Radice et al. 2016; Ciolfi et al. 2017).

When simulating the emission from short-duration gamma-ray bursts, it is challenging to balance the need for a large domain (the fireball becomes transparent at the photospheric radius ${r}_{\mathrm{ph}}\sim {10}^{13}$ cm), high spatial resolution (the thickness of the fireball remains of the order of one light second throughout the expansion), and full dimensionality (3D) in order to avoid the occurrence of the plug (or butterfly) instability (Gottlieb et al. 2017). Such instability is particularly harmful in short GRB simulations, as it can penetrate the fireball all the way to the back of the shell, resulting in a complete loss of the relativistic motion.

With current technology, it is impossible to perform a simulation in 3D with the required domain and resolution. We therefore proceed as follows. We simulate the first 0.25 s of the jet evolution in 3D in a domain $-2\times {10}^{10}\leqslant x,y\leqslant 2\times {10}^{10}$ cm, ${10}^{8}\leqslant z\leqslant 4.01\times {10}^{10}$ cm. Our finest grid has a resolution of $2.4\times {10}^{6}\,\mathrm{cm}$. The density rendering of the 3D simulation at $t=0.25\,{\rm{s}}$ is shown in Figure 1 (where the jet has reached a size $\sim 1.5\times {10}^{10}$ cm). We then project the last simulation box in 2D by performing an azimuthal average around the z-axis, and use such projection as the initial condition for the subsequent evolution in 2D cylindrical coordinates (x, z). More detail about this procedure will be given in an upcoming publication (D. López-Cámara et al. 2017, in preparation). Figure 2 shows some detail of the projection procedure. With an identical color scale (shown in the left-hand side of the figure), the three panels show the azimuthally averaged density map (left panel) and the density map in the YZ and XZ planes of the full 3D simulation (top- and bottom-right panels, respectively). The azimuthally averaged properties of the outflow (such as the density shown in the left panel of Figure 2) are used as initial conditions for the subsequent 2D evolution. The 2D simulation box extends from 10⁹ to 10¹³ cm in the polar direction, and from 0 to 10¹³ cm in the equatorial plane. Joining the two simulations, we cover 5 orders of magnitude in scale, from 10⁸ to 10¹³ cm, much more than previously accomplished (Gottlieb et al. 2017; Kathirgamaraju et al. 2017). In terms of spatial resolution, the 2D simulation has a square of side $2.44\times {10}^{7}$ cm at the inner boundary, yielding the same relative resolution as the 3D simulation. At larger distances, the highest resolution element becomes a box of side $\sim 1.56\times {10}^{9}\,\mathrm{cm}$, sufficient to fully resolve the SGRB shell of thickness ∼3 × 10¹⁰ cm. Pseudocolor images of the logarithmic density and Lorentz factor of the 2D simulation at ${t}_{\mathrm{lab}}=6.7\,{\rm{s}}$ are shown in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To compute the light curve we followed the method of Lazzati et al. (2009, 2011, 2013) in Lazzati et al. 2013, considering only material moving with a speed of at least 0.4 c, corresponding to ${\rm{\Gamma }}\gt 1.1$. Despite the extended size of the simulation domain, we were not able to simulate the fireball all the way to its photosphere, and consequently we had to extrapolate some of the outflow properties to the desired radii. In particular, if we measured a comoving temperature ${\bar{T}}^{{\prime} }$ at a radius $\bar{r}$ at which the opacity to Thomson scattering is $\bar{\tau }\gg 1$, we compute the photospheric radius as ${r}_{\mathrm{ph}}=\bar{r}\sqrt{\bar{\tau }}$, and the photospheric temperature as ${T}_{\mathrm{ph}}^{{\prime} }={\bar{T}}^{{\prime} }{\bar{\tau }}^{-1/3}$. These scalings are correct as long as the outflow is in the ballistic regime and the ejecta coast at constant Lorentz factor in a self-similar fashion. To check that the extrapolation does not affect the accuracy of the results, we performed the extrapolation from the simulation results at four different times: ${t}_{\mathrm{lab}}=3.3,6.7,13.3$, and 26.6 s. The four values of the photospheric radius, energy, and comoving temperature that we obtained were consistent with each other, yielding a posteriori proof that the ballistic assumption is correct for our simulations (with some exceptions, see below). In Figures 4–6, we plot the average of the four values with symbols and show using an error bar the dispersion of the four measurements. Finally, we adopt some correction for the computed peak frequency and energy. As shown in Lazzati (2016; see also Parsotan & Lazzati 2017), the fluid temperature measured at the photosphere underestimates by a factor ∼3 the temperature of the radiation, while the energy is overestimated by a factor that ranges between a few and 10. We therefore correct our photon peak energies upwards by a factor 3 and our luminosities downward by a factor 4. More accurate results that make use of the Monte Carlo radiation transfer code MCRaT will be presented in a future publication.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows the polar profile of the kinetic energy of the ejecta (blue squares) and of the isotropic equivalent energy in radiation (red dots, bolometric). As noted above, the symbols show the average value obtained for the four different starting points of the extrapolations and the error bars indicate their dispersion (one standard deviation). We notice three structures in the polar profile that we interpret as follows. The brightest part is the core of the outflow, for $\theta \lt 15^\circ$. This is the original injected jet, the energy profile of which has been modified by the interaction with the environment from a top-hat jet to an exponentially stratified structure with e-folding angle ${\theta }_{j,\mathrm{obs}}=2^\circ$. This exponential jet is surrounded by a hot bubble that dominates at angles between 15° and 45°. We interpret it as the jet cocoon mixed with the ambient material cocoon (Nakar & Piran 2017). It is also characterized by an exponential profile with e-folding angle ${\theta }_{\mathrm{cocoon},\mathrm{obs}}=12^\circ$. Finally, there is a fairly isotropic component with a sharp cutoff at $\theta =65^\circ$, which we interpret as the ambient medium shocked by the cocoon pressure while the jet is still trapped inside the non-relativistic ejecta from the merger. Whether or not the cutoff at $\sim 65^\circ$ is physical is unclear, as the reflecting boundary condition at $\theta =90^\circ$ might affect the hydrodynamical evolution at large angles. These three components are shown in Figure 4 with dashed lines overlaid on the kinetic energy curve. In terms of energetics, of the 10⁵⁰ erg injected in the boundary condition, $5.5\times {10}^{49}$ erg remain in the jet ($\theta \leqslant 15^\circ$), $3.8\times {10}^{48}$ erg are found in the cocoon ($15\lt \theta \leqslant 45^\circ$), and 7 × 10⁴⁷ are in the shocked ambient medium. The remainder of the energy is given to slow ejecta (${\rm{\Gamma }}\lt 1.1$). Analogous components could explain the photospheric radiated energy profile, but the cocoon radiation requires a smaller e-folding angle in order to reproduce the steeper decline of the radiative energy with angle.

Figure 5 shows the profile of the peak energy of the detected spectrum as a function of the off-axis angle. We also see in this case a three-component structure with a very high-energy core (the on-axis SGRB jet) surrounded by an X-ray dominated cocoon, eventually transitioning to a predominantly soft X-ray regime at large angles, with properties independent of the off-axis angle. We notice that the error bars in the transition regions are particularly large, indicating that it is where the different components interact, and therefore where the assumption of ballistic evolution fails and our results should be taken with caution. This is particularly obvious at $\theta \sim 40^\circ$, at the cocoon boundary with the shocked ambient medium, and at $\theta \sim 70^\circ$, where the sharp cutoff in energy is observed.

Figure 6 shows the pulse duration as a function of the off-axis angle. This is computed as the angular timescale ${t}_{\mathrm{ang}}=R/c{{\rm{\Gamma }}}^{2}$. For the cocoon and shocked ambient medium material, the angular timescale is the dominant timescale, while for the jet-dominated case the width of the fireball dominates, and what is shown may be a lower limit on the duration of the transient. We notice that the angular timescale coincides with the time delay after the merger where the transient is observed, so that a prediction of this model is that the delay between the engine formation (possibly indicated by the detection of a GW signal) and the detection of the transient should be equal to the duration of the transient itself (O. S. Salafia et al. 2017, in preparation). We find that this should be a few seconds in the cocoon-dominated regime, and much longer (minutes) for the shocked ambient medium.

Finally, we consider out to what distance the prompt EM transient from a merger observed at a given off-axis angle would be observable with the GBM on board Fermi.¹¹ Figure 4 shows the amount of energy that would be seen by the GBM instrument on board Fermi. We characterize the GBM sensitivity curve as flat within the energy range [10–150] keV and zero outside. Since the predicted photospheric temperatures are at the edge of the sensitivity band of the instrument, the detected radiation depends sensitively on the details of the spectrum. We evaluate the detected energy in two somewhat opposite and extreme cases: a purely thermal, single-temperature spectrum, and a Comptonized spectrum with photon index −2.5.¹² The latter is obtained by substituting the part of the thermal spectrum above the peak with the prescribed power-law. Comptonization of the photospheric spectrum is expected in case of trans-photospheric dissipation (Giannios 2006; Pe'er et al. 2006; Lazzati & Begelman 2010; Ryde et al. 2011; Lundman et al. 2013; Chhotray & Lazzati 2015). To determine the maximum detectable distance, we adopt a transient detection threshold of 0.7 counts per square cm per second (Meegan et al. 2009). The results are shown in Figure 7. If the transient is thermal, only fairly on-axis bursts would be detectable in the ∼200 Mpc¹³ sphere where LIGO/Virgo are expected to detect NS–NS mergers (red, thin line in Figure 7). If, instead, some sub-photospheric dissipation is present, the cocoon radiation would become detectable at wider off-axis angles, at least for bursts within ∼50 Mpc. In Figure 7, the thick blue line is computed adding a high-frequency power-law with a photon index −2.5 on top of the thermal emission. The extra energy in the power-law photons is, for all cases, significantly smaller than the energy of the blackbody spectrum. It is possible for Comptonization to add significant energy to the spectrum. We do not consider this case because it would require either strong shocks (not seen in our simulation) or some form of magnetic dissipation (which we cannot account for in our AMR–RHD simulation). An example light curve for the Comptonized pulse observed by Fermi at 30° off-axis is shown in Figure 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Lazzati et al. (2017) we discussed the possibility of detecting a short X-ray transient associated with an off-axis short GRB. The transient would be due to the expansion of the high-pressure cocoon that forms around a relativistic jet as it works its way out of a region of high ambient density. Analogous processes have been discussed for long GRBs (Ramirez-Ruiz et al. 2002; Nakar & Piran 2017) and may be responsible for the detection of X-ray flashes (Gotthelf et al. 1996; Yamazaki et al. 2002; Fynbo et al. 2004; Lamb et al. 2005; Sakamoto et al. 2005; Guidorzi et al. 2009). In L17, we made some simplifying assumptions, such as the assumption of isotropy for the cocoon material and the assumption of a Lorentz factor ${{\rm{\Gamma }}}_{\infty ,\mathrm{cocoon}}=10$ for its asymptotic expansion. Such assumptions (especially isotropy) have been cast in doubt by Gottlieb et al. (2017), who studied the afterglow emission from such a component (see also L17). In this paper, we have presented the results of a numerical simulation of an SGRB jet with properties similar to the fiducial case presented in L17. We find that the cocoon material produces an X-ray flash detectable by the GBM on board Fermi under favorable conditions (either a relatively small off-axis angle $\theta \lesssim 30^\circ$ or a distance of tens of Mpc). We also find that such transients are a few seconds long (compared to a prediction of a fraction of a second in L17) and peak at a few keV, making the detection with the GBM somewhat challenging. Qualitatively analogous results have been presented in Kathirgamaraju et al. (2017). However, their simulations are bidimensional and for a magnetized jet, so that our results are not directly comparable in detail. It should be noted that neither the three components classification in Figure 4 nor the classification as sheet discussed in Kathirgamaraju et al. (2017) can capture entirely the complexity of the interaction of a relativistic jet with the ambient medium through which it propagates. On the one hand, all components are mixed by reciprocal interaction and difficult to disentangle. On the other hand, the presence of magnetic fields, radiation drag (Chhotray et al. 2017), and mixing have additional roles in modifying the jet structure and its surrounding. Further observational work and theoretical studies are going to be fundamental for the understanding of such a complex and fascinating phenomenon. The fairly unbiased detection of mergers from GW detectors with respect of the system inclination will help us to map the angular energy distribution of the relativistic and non-relativistic outflows triggered by the merger (L17).

Contrary to the semi-analytic results in L17, our numerical results predict X-ray transients from the cocoon that are detectable only within a range of off-axis angles. Still, the possibility of detecting at least some off-axis events from their prompt cocoon emission stands, especially if sub-photospheric dissipation can add a non-thermal tail to the spectra. It should also be taken into account that different jet/ambient properties would give rise to different signals. For example, a less dense environment extended over a larger volume would produce transients with the same energy but less baryon contamination, and therefore higher temperatures (L17). More luminous engines active for a shorter time would instead produce more energetic cocoons. In extreme cases, the entire jet might be trapped in the ejecta powering a strong cocoon emission, akin to trapped jets in massive stars (Lazzati et al. 2012) and to binary merger engines that produce outflow with wider collimation (Nagakura et al. 2014). All such transients would be easier to detect with current instrumentation. Finally, we find that the radiative efficiency of the cocoon is fairly small, of the order of a few percent, and therefore the cocoon ejecta energy will be available for producing a detectable afterglow (see, e.g., L17, Gottlieb et al. 2017).

This research was in part supported by NASA ATP grant NNX17AK42G (D.L.) and NSF awards AST-1616157 (R.P.) and AST-1333514 (B.J.M.). The software used in this work was developed in part by the DOE NNSA ASC- and DOE Office of Science ASCR-supported Flash Center for Computational Science at the University of Chicago. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.