Thermal and non-thermal emission from the cocoon of a gamma-ray burst jet

TL;DR

This paper uses hydrodynamic simulations to study the thermal and non-thermal radiation emitted by the cocoon formed around a relativistic gamma-ray burst jet, predicting observable signatures in X-ray and MeV bands.

Contribution

It provides the first detailed modeling of cocoon emission spectra and lightcurves, including thermal and inverse-Compton components, linking simulations to observable GRB features.

Findings

Cocoon emission peaks in X-ray band with a duration of a few minutes.

Cocoon luminosity is comparable to GRB afterglow during plateau phase.

Inverse-Compton spectrum extends from 10 keV to MeV, offering a way to measure jet properties.

Abstract

We present hydrodynamic simulations of the hot cocoon produced when a relativistic jet passes through the gamma-ray burst (GRB) progenitor star and its environment, and we compute the lightcurve and spectrum of the radiation emitted by the cocoon. The radiation from the cocoon has a nearly thermal spectrum with a peak in the X-ray band, and it lasts for a few minutes in the observer frame; the cocoon radiation starts at roughly the same time as when $\gamma$-rays from a burst trigger detectors aboard GRB satellites. The isotropic cocoon luminosity ($\sim 10^{47}$ erg s$^{-1}$) is of the same order of magnitude as the X-ray luminosity of a typical long-GRB afterglow during the plateau phase. This radiation should be identifiable in the Swift data because of its nearly thermal spectrum which is distinct from the somewhat brighter power-law component. The detection of this thermal component would provide information regarding the size and density stratification of the GRB progenitor star. Photons from the cocoon are also inverse-Compton (IC) scattered by electrons in the relativistic jet. We present the IC lightcurve and spectrum, by post-processing the results of the numerical simulations. The IC spectrum lies in 10 keV--MeV band for typical GRB parameters. The detection of this IC component would provide an independent measurement of GRB jet Lorentz factor and it would also help to determine the jet magnetisation parameter.

Figures

Equations

$$\delta \tau$$

$$\sim$$

$${\vec{F}}^{\prime }=-\frac{c\lambda }{k_{\mathrm{es}}{\rho }^{\prime }}{\nabla }^{\prime }{e}_{\mathrm{rad}}^{\prime },$$

$$\lambda =\frac{2+\xi }{6+3\xi +{\xi }^2},$$

$$\xi =\frac{|{\nabla }^{\prime }{e}_{\mathrm{rad}}^{\prime }|}{k_{\mathrm{es}}{\rho }^{\prime }{e}_{\mathrm{rad}}^{\prime }}.$$

$${\vec{F}}^{\prime }=-\frac{c}{3k_{\mathrm{es}}{\rho }^{\prime }}{\nabla }^{\prime }{e}_{\mathrm{rad}}^{\prime },$$

$${\vec{F}}^{\prime }=-c{\vec{e}}_{\mathrm{rad}}^{\prime },$$

$$\frac{\partial e_{\mathrm{gas}+\mathrm{rad}}^{\prime }}{\partial t^{\prime }}=-{\nabla }^{\prime }\cdot F_{\mathrm{rad}}^{\prime },$$

$$\frac{\partial e_{\mathrm{gas}+\mathrm{rad}}^{\prime }}{\partial t^{\prime }}=-\frac{1}{r^{\prime }}\frac{\partial (r^{\prime }{F}_r^{\prime })}{\partial r^{\prime }}-\frac{{\partial }^{\prime }{F}_z^{\prime }}{\partial r^{\prime }}=$$

$$=-\frac{1}{r^{\prime }}\frac{\partial }{\partial r^{\prime }}\left(r^{\prime }\frac{c\lambda }{{\rho }^{\prime }{k}_{\mathrm{es}}}\frac{\partial e_{\mathrm{rad}}^{\prime }}{\partial r^{\prime }}\right)-\frac{\partial }{\partial z^{\prime }}\left(\frac{c\lambda }{{\rho }^{\prime }{k}_{\mathrm{es}}}\frac{\partial e_{\mathrm{rad}}^{\prime }}{\partial z^{\prime }}\right),$$

$$e_{\mathrm{gas}+\mathrm{rad},i,j}^{\prime ,n+1}$$

$${\int }_z^{\infty }{\kappa }_{\mathrm{es}}{\rho }^{\prime }\Gamma \mathrm{d}z=1.$$

$$t_{\mathrm{obs}}=t_k-z_j/c,$$

$$d{F}_{{\nu }_{\mathrm{obs}}}^{\mathrm{obs}}=I_{\nu }^{\mathrm{obs}}\cos {\theta }_{d}d{\Omega }_d\approx I_{\nu }^{\mathrm{obs}}\frac{d{S}_{\perp }}{d^2},$$

$$I_{{\nu }^{\prime }}^{\prime }\approx \frac{F^{\prime }}{\pi }\frac{B_{{\nu }^{\prime }}(T^{\prime })}{\int B_{{\nu }^{\prime }}(T^{\prime })d{\nu }^{\prime }}=ϵB_{{\nu }^{\prime }}(T^{\prime }),$$

$$e^{\prime }=a{T}^{\prime 4}+\frac{3}{2}\frac{k_B{\rho }^{\prime }{T}^{\prime }}{\mu m_p},$$

$$I_{{\nu }_{\mathrm{obs}}}^{\mathrm{obs}}=I_{{\nu }^{\prime }}^{\prime }{\left(\frac{{\nu }_{\mathrm{obs}}}{{\nu }^{\prime }}\right)}^3,{\nu }_{\mathrm{obs}}=\frac{{\nu }^{\prime }}{\Gamma (1-\beta \cos \alpha )}.$$

$$T_{\mathrm{obs}}=\frac{T^{\prime }}{\Gamma (1-\beta \cos \alpha )}.$$

$$F_{{\nu }_{\mathrm{obs}}}^{\mathrm{obs}}=\frac{1}{d^2}\int d{S}_{\perp }ϵB_{{\nu }_{\mathrm{obs}}}(T_{\mathrm{obs}})\delta \left(t_{\mathrm{obs}}-t+\frac{z}{c}\right).$$

$$F_{{\nu }_{\mathrm{obs}}}^{\mathrm{obs}}(t_{\mathrm{obs}})=\frac{1}{d^2\Delta t_{\mathrm{obs}}}\sum _{j,k}ϵB_{{\nu }_{\mathrm{obs}}}(T_{\mathrm{obs},\mathrm{j},\mathrm{k}})\mathrm{d}{S}_{\perp ,j,k}\mathrm{d}{t}_j,$$

$$L_{{\nu }_{\mathrm{obs}}}(t_{\mathrm{obs}})=4\pi d^2{F}_{{\nu }_{\mathrm{obs}}}^{\mathrm{obs}}.$$

$${\vec{R}}_0(t_0)=\vec{R}(r,z)-c(t-t_0)\vec{e}(\theta ,\varphi ).$$

$$I_{{\nu }_0}(\theta ,\varphi )=ϵB_{{\nu }_0}\left(T=D_{\psi }{T}_{\mathrm{ph}}^{\prime }\right)$$

$$\Lambda (\alpha )=\left(\begin{array}{ccc}\cos \alpha & 0& \sin \alpha \\ 0& 1& 0\\ -\sin \alpha & 0& \cos \alpha \end{array}\right).$$

$$\{\begin{array}{ll}\cos \theta & =\sin \tilde{\theta }\cos \tilde{\varphi }\sin \alpha +\cos \tilde{\theta }\\ \tan \varphi & =\frac{\sin \tilde{\theta }\sin \tilde{\varphi }}{\sin \tilde{\theta }\cos \tilde{\varphi }\cos \alpha -\cos \tilde{\theta }\sin \alpha }\end{array}$$

$$\{\begin{array}{ll}& {\nu }_0=\Gamma (1+\beta \cos {\tilde{\theta }}^{\prime }){\nu }_0^{\prime }=D{\nu }_0^{\prime }\\ & \cos \tilde{\theta }=(\cos {\tilde{\theta }}^{\prime }+\beta )/(1+\beta \cos {\tilde{\theta }}^{\prime })\\ & \tilde{\varphi }={\tilde{\varphi }}^{\prime }\\ & I_{{\nu }_0^{\prime }}^{\prime }({\tilde{\theta }}^{\prime },{\tilde{\varphi }}^{\prime })=I_{{\nu }_0}(\tilde{\theta },\tilde{\varphi })/D^3\end{array}$$

$$\{\begin{array}{ll}& \cos {\theta }_{\mathrm{obs}}^{\prime }=(\cos \alpha -\beta )/(1-\beta \cos \alpha )\\ & {\varphi }_{\mathrm{obs}}^{\prime }=\pi \end{array}$$

$$\begin{array}{rl}\mathrm{d}{N}_{\gamma }=& \mathrm{d}{N}_{\mathrm{e}}\mathrm{d}{\nu }^{\prime }\Delta {\Omega }_{\mathrm{obs}}^{\prime }\mathrm{d}{t}^{\prime }\\ & \int \mathrm{d}{\tilde{\Omega }}^{\prime }\int \mathrm{d}{\nu }_0^{\prime }\frac{I_{{\nu }_0^{\prime }}^{\prime }({\tilde{\theta }}^{\prime },{\tilde{\varphi }}^{\prime })}{h{\nu }_0^{\prime }}\frac{{\partial }^2\sigma }{\partial {\nu }^{\prime }\partial {\Omega }_{\mathrm{obs}}^{\prime }}{\mathrm{e}}^{-\tau },\end{array}$$

$$t_{\mathrm{obs}}=t-z/c,$$

$$\nu =[\Gamma (1-\beta \cos \alpha ){]}^{-1}{\nu }^{\prime }=D_{\alpha }{\nu }^{\prime },$$

$$\begin{array}{rl}\mathrm{d}{L}_{\nu }^{\mathrm{iso}}=& \mathrm{d}{N}_{\mathrm{e}}4\pi D_{\alpha }^{2}h\nu \\ & \int \mathrm{d}{\tilde{\Omega }}^{\prime }\int \mathrm{d}{\nu }_0^{\prime }\frac{I_{{\nu }_0^{\prime }}^{\prime }({\tilde{\theta }}^{\prime },{\tilde{\varphi }}^{\prime })}{h{\nu }_0^{\prime }}\frac{{\partial }^2\sigma }{\partial {\nu }^{\prime }\partial {\Omega }_{\mathrm{obs}}^{\prime }}{\mathrm{e}}^{-\tau }\end{array}$$

Full text

††pubyear: 2017††pagerange: Thermal and non-thermal emission from the cocoon of a gamma-ray burst jet–References

1 Introduction

Long duration gamma-ray bursts (GRBs) are produced when the core of a massive star collapses to a neutron star or a black hole (for recent reviews on GRBs see, e.g., Piran, 2004; Woosley & Bloom, 2006; Fox & Mészáros, 2006; Gehrels, Ramirez-Ruiz, & Fox, 2009; Kumar & Zhang, 2015). The newly formed compact object produces a pair of relativistic jets that make their way out of the progenitor star along the polar regions. Punching their way to the stellar surface, these jets shock heat the material they encounter pushing it both sideways and along the jet’s direction. Therefore, the jet is surrounded by a hot cocoon made by this shock heated plasma, which contributes to the collimation of the jet (Ramirez-Ruiz, Celotti, & Rees, 2002). The central engine activity is known to be highly variable and long-lived, giving multiple episodes of gamma-ray emission during the prompt phase (e.g., Ramirez-Ruiz & Merloni, 2001) and sometimes sharp increases in X-ray flux (flares) at much later time from minutes to hours (e.g., Chincarini et al., 2007).

The total amount of energy deposited in the cocoon (which is equal to the work done by the jet on the medium it passes through) depends on the size of the star and the jet luminosity. The temperature of the cocoon is determined by the density of the star which controls the expansion speed of the cocoon transverse to the jet axis and hence the cocoon’s volume. The temperature, energy, and Lorentz factor of the cocoon are the main parameters that affect its luminosity. Thus, the observation of radiation from the cocoon, once it breaks out of the star, expands and later decelerates into the stellar environment, should in principle provide information about the progenitor star.

The numerous GRB simulations performed to date (e.g., Zhang, Kobayashi, & Mészáros, 2003; Mizuta et al., 2006; Morsony, Lazzati, & Begelman, 2007; Bromberg et al., 2011; Lazzati et al., 2013; López-Cámara et al., 2013; Mizuta & Ioka, 2013; Duffell & MacFadyen, 2015; Bromberg & Tchekhovskoy, 2016) have considered the hydrodynamic or magnetohydrodynamic interactions of the jet with the progenitor star, but not the radiation escaping the cocoon or the interaction of this radiation with the relativistic jet. However, it is the radiation from the cocoon that provides information on the GRB progenitor star properties, and that is a big part of the motivation for this work.

The only paper we know presenting cocoon lightcurves based on hydrodynamical simulation is by Suzuki & Shigeyama (2013). However, this paper was mainly concerned with the impact of the circumstellar medium on the early thermal X-ray emission. Recently, Nakar & Piran (2017) provided an analytic calculation of the cocoon radiation which includes the mixing of the shocked jet and shocked stellar material, but not the IC interaction between the cocoon and the jet.

In this paper, we present lightcurves and spectra of radiation escaping the cocoon by post-processing special relativistic hydrodynamic simulations which follow the evolution of the GRB jet and the associated cocoon from $\sim 10^{8}$ cm to $\sim 3\times 10^{14}$ cm. The aforementioned mixing and its effect on the emergent cocoon radiation is build into our hydrodynamical simulations. We also compute the radiation from the cocoon scattered by electrons in the relativistic jet, which produces higher energy photons. The non-thermal radio emission from the cocoon is studied in detail in a follow-up paper (De Colle et al., 2018).

The paper is organised as follows: in Section 2 we provide information about the pre-collapse Wolf-Rayet star, the numerical method employed, and the simulation set up. The jet dynamics and the hydrodynamic properties of the cocoon are discussed in Section 3. Details of the calculation of the thermodynamic parameters and the radiation are described in Section 4, where we also present lightcurves and spectra for a number of selected jet and progenitor stellar models. Section 5 describes the inverse-Compton interaction between the relativistic jet and the cocoon photons and the lightcurve of the emergent high-energy photons; the main results of the paper are discussed in Section 6. Throughout this paper, we use the convention $G_{x}=G/10^{x}$ in cgs units.

2 Methods and Initial conditions

We run two-dimensional axisymmetric simulations using the adaptive mesh refinement code Mezcal (De Colle & Raga, 2006; De Colle et al., 2012a, b, c), which solves the special relativistic, hydrodynamics equations on an adaptive grid. The simulations span six orders of magnitude in space and time, and we consider the jet crossing the star and then moving in the wind of the progenitor star.

As progenitors of GRBs, we employ two pre-supernova, Wolf-Rayet stellar models (see Figure 1). The E25 stellar model (Heger, Langer, & Woosley, 2000) is a star with an initial mass of $25M_{\odot}$ reduced to a final mass of $5.45M_{\odot}$ after losing its hydrogen and helium envelopes by massive winds. The 12TH model (Woosley & Heger, 2006) has an initial mass of $12M_{\odot}$, a final mass of $9.23M_{\odot}$ and a more extended stellar envelope. At a radius larger than the stellar radii of these models ($r>3\times 10^{10}$ cm and $r>10^{11}$ cm for the E25 and 12TH models respectively), the ambient medium density is taken to be that of the wind of the progenitor star which we assume had a mass loss rate of $\dot{M}_{w}=10^{-5}$ M⊙ yr*-1* and a velocity of $v_{w}=10^{8}$ km s*-1*.

During a period of time $t_{\rm inj}=20$ s, a conical jet with an half-opening angle of 0.2 rad is injected from a spherical, inner boundary located at $r=5\times 10^{8}$ cm. The jet luminosity is $L_{\rm jet}=2\times 10^{50}$ erg  s*-1* (corresponding to a total kinetic energy of $4\times 10^{51}$ erg). Simulations are run with three different values of the jet Lorentz factors: 10, 20, and 30. The jet internal pressure is assumed to be a small fraction ($10^{-5}$) of the rest-mass energy density, giving a negligible amount of thermal energy in the jet111Other authors (e.g. Mizuta & Ioka 2013) studied the propagation of “hot” jets, in which a large fraction of the jet energy is initialized as thermal energy. As most of the jet kinetic energy is dissipated into thermal energy in the post-shock region when the jet moves through the star, the different choices of the initial conditions should produce similar results.. The jet is switched off after 20 seconds, and subsequently the spherical, inner region (from where the jet was injected until that time) becomes part of the computational region.

A computational box with physical size $(L_{r},L_{z})=(3\times 10^{14},3\times 10^{14})$ cm (along the $r$- and $z$-axis respectively) is resolved by using a grid with 40 $\times$ 40 cells and 20 levels of refinement, corresponding to a maximum resolution of $\sim 1.4\times 10^{7}$ cm. The jet at the inner boundary is resolved by $\sim$ 7 cells in the transverse direction. The propagation of the jet is followed for 104 seconds (time measured by a stationary observer at the centre of explosion).

As the jet expands, it would be impossible to keep the same resolution during the entire duration of the simulation (as it would require employing $\approx 10^{15}$ cells!). For this reason, we adapt dynamically the grid, refining (i.e., creating four “sibling” cells from a parent cell) and derefining (i.e., replacing four sibling cells with one parent, larger cell) the cells as a function of the evolutionary stage of the simulation and the distance $R$ from the origin of the coordinate system. This is done by decreasing the maximum number of levels as $N_{\rm max}=20-3.32\times\log(t/20)$ for $t>20$ s, in such a way that the cocoon (whose size increases with time as $R\sim ct$) is resolved by a maximum of approximately $(N_{r},N_{z})=(2000\times 2000)$ cells during its full evolution.

We also run simulations with 22 levels of refinement (within a smaller computational box). Although the details of the jet propagation (i.e. the generation of instabilities at the jet/cocoon interface) depend on resolution, the calculation of the flux and light curve shown in Section 4 do not change by more than 10% when increasing the resolution from 20 to 22 levels of refinement.

Resolving the density stratification of the wind medium outside the star requires setting initially a large number of cells which remain unused until the GRB/cocoon shocks arrive at that radius. To make the simulation faster we derefine at the lowest level of refinement all the cells outside a spherical region which is expanding at the speed of light. As this “virtual” sphere expands in the environment, the mesh is refined and the values of the density and the pressure are set using the initial conditions of the simulation. In this way, we get an improvement by a factor of $\sim$ 2 in computational time.

3 Dynamics

The different stages of the dynamical evolution of the GRB/cocoon system222We describe in this Section only the dynamical evolution of the GRB jet with $\Gamma_{\rm jet}=10$ propagating through the E25 pre-supernova stellar model. The other models present a similar dynamical evolution. are shown in Figure 2. The jet is injected from a circular boundary (top, left panel of Figure 2). Due to the large stellar density the jet propagates at sub-relativistic speed inside the star, needing about five seconds to get to $r/c\sim$ 0.45 s (top, central panel of Figure 2), which corresponds to an average velocity of $v/c\sim 0.1$. A double shock structure (the “working surface”) is created at the head of the jet, with a forward shock which accelerates the stellar material, and a reverse shock which decelerates the fast-moving jet material. The hot, dense shocked plasma expands laterally (due to pressure gradients) forming a dense, hot cocoon which helps collimate the jet. In the cocoon, the interface between the shocked jet and the shocked stellar material (the “contact discontinuity”) is corrugated by the presence of Kelvin-Helmholtz instabilities.

The top, right panel of Figure 2 shows the jet breaking out from the star. As the jet arrives at the stellar envelope, the lower entropy stellar material facilitates the lateral expansion of the cocoon, which quickly encloses the star and propagates into the wind medium. The amount of material in the hot cocoon is then given by the stellar material which has crossed the forward shock when the jet was still inside the star but did not have enough time to move laterally and dissipate its thermal energy inside the star.

The bottom panels of Figure 2 show the late evolution (left to right) of the GRB jet/cocoon system. The cocoon (bottom left panel) is strongly stratified, with a much larger density close to the jet axis and decreasing by 6-7 orders of magnitude at large polar angles. The deceleration phase is clearly visible in the right, bottom panel where Rayleigh-Taylor instabilities form at the interface between the GRB cocoon and the shocked ambient medium along the equatorial plane. Along the direction of propagation of the GRB jet, the dynamical evolution is more complex. At $t=20$ s the GRB jet is switched-off, and a rarefaction wave moves at the speed of light towards the shock front (visible at $z/c=10$ s and $r/c=0$ s in the left bottom panel). As the rarefaction shock arrives to the head of the jet, the double shock structure gradually forms a thin shell which will move with constant speed up to distances $10-100$ times larger than those simulated here before decelerating.

Figure 3 shows the velocity distribution ($u=\Gamma v/c=\Gamma\beta$) at two evolutionary times. Most of the cocoon expands at sub-relativistic speeds, except the region close to the shock front, which moves with larger velocities. The Lorentz factor strongly depends on the polar angle. At $t=200$ s (top panel) the Lorentz factor drops from $\Gamma\sim 5$ (in the region $10^{\circ}\lesssim\theta\lesssim 45^{\circ}$) to sub-relativistic speeds $\Gamma\beta\lesssim 1$ at $\theta\gtrsim 75^{\circ}$. At $t=5000$ s (bottom panel of Figure 3) the cocoon has slightly decelerated (except for the region $\theta\lesssim 30^{\circ}$).

Figure 4 shows the time evolution of the shock average velocity $R/ct$, in the lab frame, as a function of the shock polar angle. As discussed above, the jet moves at non-relativistic speeds while it crosses the star. As it breaks out of the star, it accelerates during $\sim$ 10 seconds and achieves highly relativistic speeds along the jet axis (with a Lorentz factor equal to the Lorentz factor of the material injected from the inner boundary), and mildly relativistic speeds (with Lorentz factors of $\approx$ 2-5, see Figure 4) at large polar angles. The deceleration phase starts first at larger polar angles. For instance, the cocoon shock velocity at polar angles $\theta=75^{\circ}$-$90^{\circ}$ decreases from $v\sim 0.8\;c$ at $t=100$ s to $v\sim 0.7\;c$ after $t=10^{4}$ s.

Figure 5 shows the energy distribution as a function of polar angles. The energy increases during the first 20 seconds (consistently with the duration of injection of the jet from the inner boundary). Most of the energy remains collimated in a conical region of angular size $\theta\lesssim 15^{\circ}$ during the full duration of the numerical simulation. A smaller amount of energy ($10^{49}$-$10^{50}$ erg) is present in the portion of the cocoon moving at larger angles. While along the jet axis most of the energy is concentrated in a small region around the shock, at larger polar angles a significant fraction of the energy is distributed through all the cocoon’s volume. Therefore, the deceleration radius observed in the numerical simulations is smaller than one inferred by using the expression $R/c\sim Ev_{w}/(\dot{M}_{w}\Gamma^{2}c^{3})\approx 6\times 10^{3}E_{50}/\Gamma_{1}^{2}$ s.

It is very challenging, even with the excellent computational resources available these days, to run numerical simulations up to several times $10^{14}$ cm in 3D. The effect of 3D versus 2D simulations has been recently studied by, e.g., Gottlieb et al. (2018); Harrison et al. (2018). These authors showed that the presence of asymmetries in the jet (i.e., along the direction perpendicular to the jet axis) reduces the velocity of the cocoon thereby making the cocoon dimmer, with larger differences depending on the amount of asymmetry present in the jet333A 3D simulation with cylindrical symmetry would obviously recover the results of the 2D axisymmetric simulations unless the code is not, by construction, preserving symmetries (i.e. if the jet is precessing, wiggling and so on).

4 Cocoon luminosity and spectrum

4.1 Radiative diffusion

The average optical thickness of each cell in our simulation is

[TABLE]

where $N_{\rm cells}=2000$ is the maximum number of cells in the simulations along the radial and transverse directions (the cocoon terminal Lorentz factor for mass $10^{-4}M_{\odot}$ and energy of 1051 erg is $\sim 5$). However, the whole cocoon remains optically thick along the polar axis for $\sim 10^{3}$ s when the cocoon reaches $R\sim 3\times 10^{13}$ cm. The cocoon is nearly isobaric throughout the volume (except near the surface) and the pressure is dominated by radiation. Therefore, the fractional change to the thermal energy of an interior cell in a dynamical time is of order $(N_{\rm cells}\,\delta\tau)^{-1}\ll 1$. Therefore, for cells below the photosphere, diffusive radiation has a negligible effect on the dynamics, and post-processing of the hydro-simulation output to calculate cocoon radiation is expected to yield a reasonably accurate result at least for small polar angles $\lesssim 45^{\rm o}$. At large polar angles, the cocoon structure is affected by the loss of thermal energy at late evolutionary times. A full special relativistic radiation hydrodynamics calculation is necessary to properly compute the cocoon radiation seen by an observer located off-axis. In this paper, we only consider the case where the observer is on-axis at $\theta_{\rm obs}=0$.

Although radiative diffusion does not change the hydrodynamical evolution of the system, it strongly modifies the temperature structure of the cocoon surface and the resulting cocoon radiation. In our post-process procedure, radiative diffusion is coupled to the hydro-simulations as follows. Given a snapshot of the numerical simulation, the flux diffusing across the boundaries of each cell is calculated by employing the flux-limited diffusion approximation in the comoving frame of each cell

[TABLE]

where $e^{\prime}_{\rm rad}$ is the (comoving) radiation energy density (defined as $e^{\prime}_{\rm rad}=aT^{\prime,4}$, with the proper temperature defined in equation 15) of nearby cells, $\rho^{\prime}$ is the matter density, and $k_{\rm es}=0.2$ cm2 g*-1* is the electron scattering opacity for fully ionised helium (or heavier elements). The parameter $\lambda$ is the flux limiter, given by (Levermore & Pomraning, 1981):

[TABLE]

where $\xi$ is a dimensionless quantity defined as

[TABLE]

Equation (3) gives the correct scaling in the optically thin and optical thick limits, although it represents only an approximation of the intermediate case. In an optically thick medium $\xi\ll 1$ and $\lambda=1/3$, equation (3) reduces to

[TABLE]

while in the optically thin limit $\xi\gg 1$ and the flux is limited to

[TABLE]

where $\vec{e}_{\rm rad}^{\prime}=\nabla e_{\rm rad}^{\prime}/|e_{\rm rad}^{\prime}|$.

Once the flux is computed using equation (3), the effect of the radiation is included in the energy equation by solving the following equation

[TABLE]

which in cylindrical coordinates takes the following form

[TABLE]

and is discretized as

[TABLE]

where $i^{\pm}=i\pm 1/2$, $\nu_{d}=\frac{c\lambda}{\rho k_{\rm es}}$, $\Delta e^{\prime}_{i^{+}}=e^{\prime}_{{\rm rad},i+1,j}-e^{\prime}_{{\rm rad},i,j}$ and the indexes $i,j$ refer to the cell $(i,j)$ of the computational grid, with $i=1,\dots,N_{r}$, $j=1,\dots,N_{z}$ ($N_{r}$ and $N_{z}$ being the total number of cells along the $r$ and $z$ directions respectively). Finally, we notice that the stability condition restricts the time step to $\Delta t\leq(\Delta x)^{2}/4\nu_{d}$.

4.2 Calculation of luminosity and spectrum

For each simulation, we save a large number of snapshots (typically 103). After post-processing each snapshot with the radiative diffusion algorithm described above and assuming that the observer is located along the $z$-axis, we determine the position of the photosphere for each snapshot, defined as the surface corresponding to an optical depth $\tau=1$, i.e.

[TABLE]

For each cell $(j,k)$ located on the photosphere we compute the specific intensity at the location of the observer and assign it to the corresponding observing time, related to the simulation time by

[TABLE]

where $z_{j}$ is the distance from the equatorial plane and $t_{k}$ is the simulation time in the lab frame. The duration of the snapshot $k$ can be defined as $\Delta t_{k}=(t_{k+1}-t_{k-1})/2$, and the flux in the observer’s frame lasts for $\Delta t_{\rm obs}=(1-\beta_{j}\cos\alpha)\Delta t_{k}$, where $\alpha$ is the angle between the velocity vector and the $z$-axis and $\beta_{j}$ is the photosphere velocity. As our purpose is to compute the specific and total luminosity, we neglect cosmological effects and take $z=0$.

To compute the flux we add up the contributions from all cells at the photosphere (see, e.g., De Colle et al., 2012a)

[TABLE]

where d is the distance to the source, $\theta_{d}$ is the angle between the line joining the observer to the centre of the star and the line from the observer to the volume element from which radiation is being considered, and $d\Omega_{d}=dS_{\perp}/d^{2}$ is the differential solid angle of the cell as viewed by the observer. As the angular size of the source according to the observer is negligibly small, it is justified to take $\cos\theta_{d}\approx 1$. The specific intensity in the comoving frame is obtained by assuming a blackbody spectrum

[TABLE]

where $F^{\prime}$ is given by equation (3).

The comoving-frame temperature $T^{\prime}$ is computed from the values of energy (obtained at each timestep by using equation 10) and density by inverting the equation

[TABLE]

where $\mu=4/3$ is the mean molecular weight. The specific intensity in the observer’s frame is given by

[TABLE]

We employ the fact that the blackbody photon distribution function is invariant under Lorentz transformation, with the temperatures in the two frames related by the following relation

[TABLE]

We compute the total flux in the observer frame (see, e.g., De Colle et al. 2012a) by integrating over the entire photosphere

[TABLE]

Integrated and averaged over the time interval $\Delta t_{\rm obs}$, this equation reduces to

[TABLE]

where the sum extends over all $j$-th snapshots (one at each lab frame time $t_{j}$) and $k$-th cells respectively. The luminosity in the observer’s frame is

[TABLE]

4.3 Results

We present simulation results for cocoon luminosity and spectra for two different GRB progenitor star models (E25 and 12TH of Heger, Langer, & Woosley 2000; Woosley & Heger 2006). For the model E25, three different jet Lorentz factors (10, 20 and 30) are considered.

Figure 6 shows the bolometric cocoon luminosity computed with the method described in the previous section. The luminosity increases quickly to $\sim 10^{45}-10^{46}$ erg/s as the jet breaks out of the star and increases to values of $\sim 10^{47}$ erg/s after 100 s (in the observer’s frame) before dropping as $\sim t_{\rm obs}^{-3}$ for larger observing time.

The X-ray lightcurve of a sample of GRBs (Perley et al., 2014) is compared in figure 7 with the X-ray flux emitted by the cocoon. The cocoon energy is a fraction of the energy deposited by the jet into the star (i.e. the jet luminosity integrated over the time needed for the jet to cross the star), i.e. $\sim$ 1/100-1/10 of the GRB total energy (1049-1050 erg, see Figure 5). The GRB kinetic luminosity (imposed as 1050 erg/s in our simulations) is about three orders of magnitude larger than the cocoon luminosity ($\sim$ 1047 erg/s). The cocoon luminosity lasts for a much longer period of time ($\sim$ several hundred of seconds, see 6), and the total energy emitted by the cocoon is about $\sim 5\times 10^{49}$ erg. The cocoon X-ray luminosity is much smaller than GRB jet afterglow luminosity for the initial 30-50 s. Subsequently, as the GRB jet X-ray lightcurve undergoes a very steep decline for a few minutes after the end of the prompt phase, the cocoon luminosity becomes of the same order or slightly lower than the afterglows of GRBs with low isotropic energies ($E_{\rm iso}<10^{52}$ erg), and about one order of magnitude lower than GRBs with large isotropic energies ($E_{\rm iso}>10^{52}$ erg). This can be seen by comparing the cocoon luminosity with e.g. the XRT flux for GRBs in the 1051-1053 erg range (see Figure 7). Thus, the cocoon radiation should be detectable at least for some GRBs and can be identified by its distinct thermal spectrum.

Figure 6 also shows that the stellar structure is the key parameter which determines the cocoon luminosity. When a more extended stellar model is considered, the lightcurve shows a peak luminosity slightly lower and is in general dimmer at all observing times. As the jet break-out happens at later times, the luminosity increases on a longer timescale. Finally, increasing the jet Lorentz factor (while keeping fixed the total energy of the jet) does not produce a large change in luminosity (see the bottom panel of Figure 6). This result is consistent with the cocoon energy being determined by the jet ram pressure/luminosity (which is the same in all models considered here), and on the time it takes for the jet to carve out a polar cavity through the GRB progenitor star, and not on the jet Lorentz factor (the speed at which the jet moves through the star is a weak function of $\Gamma_{j}$).

Figure 8 shows the angular contribution to the spectrum at $t_{\rm obs}=30$ s for the E25, $\Gamma_{\rm jet}=10$ model. The spectrum peaks in the X-ray band were most of the energy is emitted. Actually, a plot of the X-ray luminosity (in the 0.3–10 keV band, which corresponds to Swift XRT energy coverage) would be indistinguishable from the bolometric luminosity shown in Figure 6. The spectrum is nearly thermal (a blackbody curve is also shown in the Figure for comparison), with a low-frequency spectrum $F_{\nu}\propto\nu^{2}$ and an exponential decay at large frequencies, and wider than a blackbody spectrum near the peak. Due to relativistic beaming, the flux is dominated by the emission at small angles ($\theta\lesssim 30^{\circ}$). Angles larger than $45^{\circ}$ produce a negligible flux when the jet is observed on-axis (they would dominate the cocoon’s emission for off-axis GRBs).

Figures 9 and 10 show the time evolution of the spectrum and energy peak respectively at $t_{\rm obs}=30,100$ s. The cocoon spectrum for a more extended stellar model (Figure 9, central panel) is less similar to a blackbody, and its evolution with time is more rapid. Also shown in the Figure 9 are spectra for three different values of jet Lorentz factor; the spectrum below $\lesssim$ 1 keV is nearly independent of $\Gamma_{\rm jet}$, however, at higher energies the flux is larger for cocoons formed by jets with larger $\Gamma_{\rm jet}$.

4.4 Observations of thermal X-ray radiation

Several groups searched for thermal X-ray emission after the prompt $\gamma$-ray phase (e.g. Friis & Watson, 2013; Starling et al., 2012; Sparre & Starling, 2012, and references therein). Starling et al. (2012) analysed spectra for 11 GRBs and found thermal X-ray signal for three of these (GRBs 060218, 090618 and 100316D) which are known to be associated with supernovae. In 4 additional cases the spectral fit is claimed to improve when a thermal component is included; these 4 cases are also possibly associated with supernovae. Sparre & Starling (2012) analysed Swift/XRT data for 190 GRBs and identified 6 bursts with possible blackbody components444There might have been many more bursts with a thermal afterglow component which the analysis of Sparre & Starling (2012) could not identify due to the uncertainty associated with hydrogen column density and the absorption of soft X-ray radiation. (GRBs 061021, 061110A, 081109, 090814A, 100621A and 110715A). The three new candidates identified by Sparre & Starling (2012, GRBs 061021, 061110A and 090814A) are not known to have an associated supernova. The thermal components in these 6 cases were found during the steep decline or the plateau phase of the X-ray lightcurve.

One of the models suggested for the thermal component in GRB afterglow lightcurves is that it is the radiation associated with the supernova shock breaking through the surface of the progenitor star. However, several papers (Ghisellini, Ghirlanda, & Tavecchio, 2007; Li, 2007; Chevalier & Fransson, 2008) point out a problem with this scenario, and that is that the thermal luminosity is too large to be consistent with the shock-breakout model for GRB 060218; the thermal luminosity for the other five cases are even larger than GRB 060218 (Sparre & Starling, 2012).

The observed luminosity for the thermal component of these GRB afterglows is between 1047 and 1050 erg s*-1*, and the spectral peak is at $\sim 0.5$ keV (Sparre & Starling, 2012). These values are consistent with the expectations of cocoon radiation (see Figs. 6 and 10). Measurement of the evolution of the thermal component’s peak temperature and luminosity in the future and its comparison with the cocoon simulation would make it possible to draw a firm conclusion regarding the origin of the radiation. If confirmed to be originated by the cocoon associated with the relativistic jet, then that could be used to determine the progenitor star radius and density structure.

5 Inverse-Compton scattering of cocoon radiation by the relativistic jet

As thermal photons from the cocoon diffuse out from the photosphere and run into the relativistic jet, they can get inverse-Compton (IC) scattered by electrons to much higher energies. This will produce a flash of high energy photons when the faster jet surpasses the cocoon photosphere (Kumar & Smoot, 2014). Because of the relativistic beaming along the direction of cocoon local velocity vector, only photons emitted within an angular patch of size $\sim\Gamma_{\rm c}^{-1}$ ($\Gamma_{\rm c}$ being the cocoon Lorentz factor) surrounding the jet can interact with electrons in the jet and contribute to IC luminosity. Since the angle between the jet axis and photons’ moving direction is of order $\Gamma_{\rm c}^{-1}$, IC scatterings boost cocoon photon energy to $\sim kT_{\rm c}(\Gamma_{\rm j}/\Gamma_{\rm c})^{2}\sim\Gamma_{\rm j,2.5}^{2}/\Gamma_{\rm c,1}^{2}$ MeV when electrons are cold ($\gamma_{\rm e}\sim 1$) in the jet comoving frame.

We present in this section the IC lightcurves and spectra computed by post-processing the results of the hydrodynamical simulations. To calculate the IC emission in real GRBs, we consider the following two cases. (i) In the Early-IC case, a fast jet with Lorentz factor $\Gamma_{\rm j}=300$ is launched with a delay of $t_{\rm delay}=30$ sec, owing to episodic central engine activity (e.g., Ramirez-Ruiz, Merloni, & Rees, 2001; Ramirez-Ruiz & Merloni, 2001). (ii) In the Late-IC case, a slower jet with Lorentz factor $\Gamma_{\rm j}=50$ is launched with a delay of $t_{\rm delay}=100$ sec, owing to late time central engine activity causing X-ray flares. The cocoon physical properties as well as the thermal radiation produced are only weakly dependent on the jet Lorentz factor. Thus, the approach used to compute the IC emission is also nearly independent on the prompt jet Lorentz factor, which, for numerical limitations, is smaller in our simulations than in typical GRBs.

In both cases, we assume that the central engine stays active at a low luminosity level so as to keep the jet funnel open, so the delayed jet does not interact with the cocoon hydrodynamically. If the shocked cocoon gas refills the jet funnel, then the delayed jet will be significantly decelerated if $\Gamma_{j}^{2}(M_{c}c^{2}/\Gamma_{c})>E_{j}$, where $\Gamma_{c}$ is the cocoon Lorentz factor, $M_{c}$ is the cocoon mass within the jet funnel, and $E_{j}$ is the kinetic energy of the delayed jet. Modeling the possible hydrodynamical interaction between delayed jets and the cocoon is out of the scope of the current paper. Shocks at the jet-cocoon interface may accelerate relativistic electrons, which will inverse-Compton scatter the cocoon radiation and produce additional non-thermal spectrum. These effects need to be studied in detail in a future work.

In this section, all quantities in the comoving frame of the fluid cell are denoted with a prime and unprimed quantities are measured in the lab frame (rest frame of the central engine).

5.1 Method

For each snapshot from the simulation, we first identify the position of the cocoon photosphere (equation 11). Each photospheric cell has its own diffusive flux $F_{\rm ph}^{\prime}$ (given by flux-limited diffusion in equation 3), temperature $T_{\rm ph}^{\prime}$, Lorentz factor $\Gamma$, and the angle between the fluid cell’s velocity vector and the jet axis is $\arctan(v_{r}/v_{z})$. The intensity in the comoving frame of each photospheric cell is assumed to be isotropic (equation 14) and the lab-frame intensity in an arbitrary direction is given by the corresponding Lorentz transformation.

Then, for each scattering cell located in the hypothetic jet at time $t$ at the position $\vec{R}(r,z)=(r,0,z)$, we consider light rays coming from different directions denoted by the unit vector $\vec{e}(\theta,\phi)$, where $\theta$ is the polar angle and $\phi$ is the azimuth angle. At an earlier time $t_{\rm 0}$ the position of the light ray is

[TABLE]

For each of the previous snapshots at $t_{\rm 0}<t$, we do not expect $\vec{R}_{\rm 0}(t_{\rm 0})$ to be exactly on the photosphere surface (since the snapshot times $\{t_{\rm 0}\}$ are discrete), so the closest photospheric cell is considered as where the light ray was emitted. Therefore, the intensity of the cocoon’s radiation field at any time $t$ at the position of the scattering cell $\vec{R}(r,z)$ is given by

[TABLE]

where $T_{\rm ph}^{\prime}$ is the temperature of the photospheric cell found by the ray-tracing method above, $\epsilon=F_{\rm ph}^{\prime}/[\sigma_{\rm SB}(T^{\prime})^{4}]$, $D_{\psi}=[\Gamma_{\rm 0}(1-\beta_{\rm 0}\cos\psi)]^{-1}$ is the Doppler factor, and $\psi$ is the angle between $\vec{e}(\theta,\phi)$ and the velocity vector of the photospheric cell. Then, it is straightforward to integrate light rays of different directions and frequencies to calculate the IC lightcurve and spectrum.

In the following, we describe the detailed integrating procedure presented in Lu, Kumar, & Smoot (2015). Experienced readers may jump to the results in the next subsection.

Consider a scattering cell at time $t$ and position $\vec{R}(r,z)$ with a total number of electrons denoted by $\mathrm{d}N_{\rm e}$. The cell is moving at an angle $\alpha$ from the $z$-axis with velocity $\vec{v}$ which corresponds to a Lorentz factor $\Gamma=[1-(v/c)^{2}]^{-1/2}$. We define another set of Cartesian coordinates $\tilde{x}\tilde{y}\tilde{z}$ in which $\hat{\tilde{z}}$ is parallel to the cell’s velocity vector $\vec{v}$, $\hat{\tilde{y}}$ is parallel to the original $y$-axis $\hat{y}$, and $\hat{\tilde{x}}$ is in the original $xz$ plane (which goes through the scattering cell). The relation between the two coordinate systems is described by a rotation along the $y$-axis by an angle $\alpha$ according to the right-hand rule, and the transformation matrix (in Cartesian coordinates) is

[TABLE]

In this new coordinate system, the original direction of the light ray $\vec{e}(\theta,\phi)$ becomes $\vec{\tilde{e}}({\tilde{\theta}},{\tilde{\phi}})=\Lambda\vec{e}$. This relation gives the mapping between $(\theta,\phi)$ and $({\tilde{\theta}},{\tilde{\phi}})$ as follows

[TABLE]

Then, we convert the intensity of external radiation field $I_{\nu_{\rm 0}}({\tilde{\theta}},{\tilde{\phi}})$ to the comoving frame $\tilde{x}^{\prime}\tilde{y}^{\prime}\tilde{z}^{\prime}$ of the scattering cell by using the Lorentz transformations

[TABLE]

where we have defined a new Doppler factor $D$ in the first expression. We are interested in the photons scattered towards the observer, who in the original $xyz$ frame is located on the $z$-axis and in the direction $(\theta_{\rm obs}^{\prime},\phi_{\rm obs}^{\prime})$ in the comoving $\tilde{x}^{\prime}\tilde{y}^{\prime}\tilde{z}^{\prime}$ frame given by

[TABLE]

The number of photons scattered into the solid angle $\Delta\Omega_{\rm obs}^{\prime}$ around the $(\theta_{\rm obs}^{\prime},\phi_{\rm obs}^{\prime})$ direction in a duration $\mathrm{d}t^{\prime}$ and frequency range $\mathrm{d}\nu^{\prime}$ can be obtained by integrating the radiation incoming at different frequencies $\nu_{\rm 0}^{\prime}$ and from different directions $\tilde{\Omega}^{\prime}=({\tilde{\theta}^{\prime}},{\tilde{\phi}^{\prime}})$, i.e.

[TABLE]

where $\frac{\partial^{2}\sigma}{\partial\nu^{\prime}\partial\Omega^{\prime}}$ is the differential cross section and the optical depth $\tau$ describes the attenuation by the part of the jet that lies along the incident photons’ trajectory before they enter the scattering cell. These $\mathrm{d}N_{\gamma}$ photons will arrive at the observer at time

[TABLE]

and the lab-frame frequency is given by

[TABLE]

where we have defined another Doppler factor $D_{\alpha}$. Differentiating equation (28) gives $\mathrm{d}t_{\rm obs}=(1-\beta\cos\alpha)\mathrm{d}t=\mathrm{d}t^{\prime}/D_{\alpha}$ and the Lorentz transformation of solid angles gives $\Delta\Omega_{\rm obs}^{\prime}=D_{\alpha}^{2}\Delta\Omega_{\rm obs}$ , so this scattering cell contributes a specific luminosity of

[TABLE]

where we have used $\Delta\Omega_{\rm obs}=4\pi$ for isotropic equivalent luminosity. As integrating equation (30) is computationally expensive, we made the following two simplifications. (i) When electrons are cold ($\gamma_{\rm e}^{\prime}=1$), the differential cross section is assumed to be isotropic and every electron has Thomson cross section. (ii) When electrons are hot ($\gamma_{\rm e}^{\prime}\gg 1$), we use the full angle-dependent Klein-Nishina differential cross section, but instead of the blackbody seed spectrum $I_{\nu_{\rm 0}^{\prime}}^{\prime}({\tilde{\theta}^{\prime}},{\tilde{\phi}^{\prime}})=B_{\nu_{\rm 0}^{\prime}}(D_{\psi}T_{\rm ph}^{\prime}/D)$, we use a $\delta$-function centred at $h\nu_{\rm 0}^{\prime}=2.7D_{\psi}kT_{\rm ph}^{\prime}/D$.

At last, we add up the IC luminosity from all the scattering cells in different snapshots weighted by their individual numerical timesteps.

5.2 Results

In this subsection, we present the lightcurves and spectra for the Early-IC and Late-IC cases. In the Early-IC case, the jet is launched at lab-frame time $t_{\rm delay}=30\rm\ s$ with isotropic equivalent power $L_{\rm j}^{\rm iso}=10^{52}\rm\ erg\ s^{-1}$ and Lorentz factor $\Gamma_{\rm j}=300$. We have two sub-cases corresponding to $\gamma_{e}=1$ (cold electrons) and $\gamma_{e}=10$ (hot electrons). In the Late-IC case, the jet is launched at $t_{\rm delay}=100\rm\ s$ with Lorentz factor $\Gamma_{\rm j}=50$, and we consider two luminosities $L_{\rm j}^{\rm iso}=10^{51}$ and $10^{52}\rm\ erg\ s^{-1}$. The duration of the delayed jets in all cases is $t_{\rm j}=10\rm\ s$, which is motivated by the observation that the typical ratio between the duration and peak time for X-ray flares is 0.1–0.3 (e.g. Chincarini, Moretti, et al., 2007, ApJ, 671, 1903)

The calculations in all cases are done by using an opening angle of $2/\Gamma_{j}$. The flux contributed by high-latitude ($>1/\Gamma_{j}$) regions is strongly suppressed due to relativistic beaming. We have also tested other opening angles of $1.5/\Gamma_{j}$ and $3/\Gamma_{j}$ and the differences are negligible.

We ignore the hydrodynamical interaction between the jet and cocoon and only consider IC scattering after the jet emerges from below the photosphere.

The jet emerges from the cocoon surface at radius $R_{\rm em}=10^{13}R_{\rm em,13}$ and $10^{13.5}R_{\rm em,13.5}\rm\ cm$ for the Early-IC and Late-IC cases respectively, so the peak time is

[TABLE]

The thickness of the transparent shell at the jet front that external photons can penetrate through is

[TABLE]

which needs to be compared to the size of the causally connected region given by

[TABLE]

The duration of the peak IC luminosity is given by

[TABLE]

In the Early-IC case, we have $\Delta t_{\rm peak}\simeq 0.26\rm\ s$. In the Late-IC cases with $L_{\rm j}^{\rm iso}=10^{51}$ and $10^{52}\rm\ erg\ s^{-1}$, we have $\Delta t_{\rm peak}\sim 3.7$ and $0.4\rm\ s$ respectively. We plot the lightcurves and spectra in Figure 11-14 and the main results are summarized as follows (see the captions for details):

(1) When electrons are cold, IC scattering generally taps a fraction of a few$\times 10^{-4}$ to $10^{-3}$ of the total energy of the delayed jet. This means the IC emission is likely overwhelmed by other radiation mechanisms (it is not responsible for generating the majority of the emission observed in X-ray flares.) However, when there is a modest amount of jet dissipation so that the root mean squared Lorentz factor $\sqrt{<\gamma_{\rm e}^{2}>}\gtrsim 10$ at radii between $10^{13}$ and $10^{14}\rm\ cm$, the peak IC luminosity exceeds the jet luminosity. This means that a fraction of a few to 10 percent of the jet is strongly dragged by the IC force to a lower Lorentz factor. Further jet dissipation such as internal shocks or magnetic reconnection could be triggered by this IC drag. A significant fraction of the 0.1-1 GeV emission during the prompt phase may be contributed by IC emission off cocoon photons.

(2) The spectrum is always broader than blackbody, even when electrons are in monoenergetic distribution (cold or hot). A power-law electron distribution will lead to a power-law spectrum in the high frequency part. Therefore, the IC component off cocoon emission might have been missed in previous studies looking for a thermal component.

(3) The effect of progenitor star’s mass and density profile has a much weaker effect on the IC emission than the Lorentz factor of the prompt jet. This is due to the strong dependence of IC emission on the cocoon Lorentz factor (Kumar & Smoot, 2014). It may be challenging to infer progenitor properties from the IC emission. Since the IC emission is very sensitive to the Lorentz factor of the delayed jet (peak energy $\propto\Gamma_{\rm j}^{2}$ and luminosity $\propto\Gamma_{\rm j}^{4}$), IC flux could provide an independent measurement of the $\Gamma_{\rm j}$.

Finally, we notice that internal shocks moving close to the head of the jet should produce an IC emission similar to that obtained in the case of a 30 s delay.

6 Conclusions

We have carried out numerical simulations of a long-GRB jet propagating through the progenitor star and its wind, the production of a cocoon that results from this interaction, and the spectrum and lightcurve of the emergent cocoon radiation. The dynamical evolution of the jet/cocoon was followed from $10^{8}$ cm to $\sim 3\times 10^{14}$ cm. We considered two different progenitor stars of mass 5.45 M⊙ and 9.23 M⊙ (right before the collapse) with radius of $3\times 10^{10}$ and $10^{11}$ cm respectively; they were models E25 (Heger, Langer, & Woosley, 2000) and 12TH (Woosley & Heger, 2006). The simulations were run for a luminosity $L_{\rm jet}$ = $2\times 10^{50}$ erg s*-1* and several different jet Lorentz factors ($\Gamma=10$, 20, 30). In each of these cases the jet duration was taken to be 20 s.

The cocoon emission was calculated by post-processing results of the numerical simulations. The cocoon spectrum is quasi-thermal and peaks in the X-ray band at $\sim 5$ keV ($\sim 0.5$ keV) a few seconds ($\sim 100$ s) after the cocoon emerges above the stellar surface. The bolometric luminosity of cocoon emission is $\sim 10^{47}$ erg s*-1* for about 200 s in the host galaxy rest frame ($\sim 10$ min in the observer frame for a typical redshift of 2); this luminosity is comparable to the GRB X-ray afterglow luminosity during the plateau phase which is observed in a good fraction of long-GRBs starting at about 100 s after the prompt $\gamma$-ray emission ends.

The cocoon lightcurve is nearly independent of the jet Lorentz factor when $\Gamma_{\rm jet}>20$, but depends on the stellar structure (Figure 6). When a more extended stellar model is considered, the X-ray light curve increases at later times and is dimmer, while the spectra is softer at all observing times (see Figure 9). We note that the velocity, density and energy distributions depend strongly on the jet initial conditions, in particular on the jet luminosity history (for instance, a jet with a larger luminosity can deposit larger amount of energies in the cocoon making it brighter in X-rays), on the presence of a magnetic field (which can provide extra collimation to the jet and thereby reduce the cocoon energy), on the stellar structure (as shown in this paper) and on the jet structure. In addition, the presence of large asymmetries in the jet (which need to be studied by three dimensional simulations and could be due, e.g., to precession or wiggling of the jet) would also affect the results by reducing the velocity of the cocoon then making the cocoon dimmer. Given the uncertainties in the jet characteristics, nevertheless, we showed in this paper that the thermal emission from a GRB cocoon could be detectable at least in some GRBs.

Thus, detection of a quasi-thermal component in the X-ray afterglow lightcurves of long-GRBs can be used, combined with constraints from optical observations from the associated jet-driven supernova and more detailed hydrodynamic calculations (e.g., including a more complete treament of the radiation transfer and a broad range of progenitor properties), to infer the density structure and radius of the progenitor star.

Observations of X-ray flares are usually intepreted as evidence of late energy injection. A small fraction of photons from the cocoon pass through the relativistic jet and are inverse-Compton scattered by electrons in the jet to higher energies. We assume that electrons are in monoenergetic distribution with Lorentz factor $\gamma_{\rm e}$ in the jet comoving frame. We computed the IC spectrum generated by the interaction between the cocoon photons and shocks with a delay of 30 and 100 s. For a jet of Lorentz factor 300 launched with a delay time of 30 s, the resulting IC luminosity is of $\sim 3\gamma_{\rm e}^{2}\%$ of jet luminosity and the IC spectrum peaks at $\sim 3\gamma_{\rm e}^{2}$ MeV when the jet emerges above the cocoon photosphere. Photons from the cocoon are also scattered by jets with longer delays (those responsible for X-ray flares during the GRB afterglow phase), and the resulting IC luminosity is $\sim 0.3\gamma_{\rm e}^{2}\%$ of late jet luminosity with the spectral peak at $\sim 10\gamma_{e}^{2}$ keV, if the late jet has a Lorentz factor of 50. These results on IC scatterings of cocoon photons by relativistic jets are consistent with the analytical calculations of Kumar & Smoot (2014) and would provide an indirect measurement of the GRB jet Lorentz factor and magnetization parameter if detected.

Acknowledgements

We thank Diego Lopez Cámara and Hiro Nagataki for useful discussions. F.D.C. thank the UNAM-PAPIIT grants IA103315 and IN117917. W.L. is funded by a graduate fellowship (“Named Continuing Fellowship”) at the University of Texas at Austin. ER-R acknowledges financial support from the David and Lucile Packard Foundation and UCMEXUS (CN-12-578). The simulations were performed on the Miztli supercomputer at UNAM.