SIMULATIONS OF GAMMA-RAY BURST JETS IN A STRATIFIED EXTERNAL MEDIUM: DYNAMICS, AFTERGLOW LIGHT CURVES, JET BREAKS, AND RADIO CALORIMETRY

The dynamics of gamma-ray burst (GRB) outflows depends on the density distribution of the ambient medium as well as on the structure of the relativistic expanding ejecta (e.g., Meszaros et al. 1998). Up to the deceleration epoch, where most of the energy is transferred to the shocked external medium, the dynamics is regulated by the local radial structure of the ejecta, while at later times (as the blastwave decelerates) it mainly depends on its global angular structure. In the absence of characteristic scales, self-similar, spherically symmetric solutions exist (Blandford & McKee 1976, hereafter Blandford–McKee) and they are widely used to interpret observational data on GRB afterglows. However, even the simplest departure from this ideal model could drastically modify the afterglow behavior. Anisotropies in the GRB outflow, for example, affect the afterglow light curve when the mean jet energy per solid angle within the visible region evolves significantly. As the jet decelerates, the relativistic beaming weakens and the visible region increases. If the outflow is collimated into a narrow jet with reasonably sharp edges, this occurs at the time when the bulk Lorentz factor Γ equals the inverse of the jet half-opening angle θ₀. A simple analytic calculation using the usual scaling laws leads then to a steepening of the afterglow flux decay rate, known as a jet break (Rhoads 1997; Sari et al. 1999; Kumar & Panaitescu 2000). It is however clear from numerical studies that such simple scalings do not provide an accurate description of the afterglow (Granot et al. 2001; Zhang & MacFadyen 2009; Meliani & Keppens 2010; van Eerten et al. 2010b; Wygoda et al. 2011; van Eerten & MacFadyen 2011; van Eerten et al. 2012). Such numerical studies have so far been limited to the case of a uniform external density while the interaction of relativistic GRB jets with a non-uniform medium remains poorly understood.

Motivated by this, here we study the dynamics of two-dimensional axially symmetric impulsive jets propagating in a spherically symmetric stratified medium of rest-mass density ρ = Ar^−k and the resulting afterglow emission. Since long duration GRBs (Gehrels et al. 2009) have massive star progenitors whose winds are expected to modify their immediate surroundings (Chevalier et al. 2004; Ramirez-Ruiz et al. 2005; van Marle et al. 2008; Mimica & Giannios 2011), we consider both steady and time varying stellar winds as possible surrounding or external media for the GRB jet evolution. The case k = 2 corresponds to a stellar wind for a massive star progenitor (Chevalier & Li 2000; Panaitescu & Kumar 2000; Ramirez-Ruiz et al. 2001; Wu et al. 2005) with a constant ratio of its pre-explosion mass loss rate $\dot{M}_w$ and wind velocity v_w, in which case ρ = Ar⁻², where $A =\dot{M}_w/(4\pi v_w)$. However, since the dependence of $\dot{M}_w$ and v_w on the time t_w before the stellar explosion that triggers the GRB is highly uncertain, it is worth considering other values of k. For example, if $\dot{M}_w\propto t_w^a$ and v_w∝t^b_w then the location of a wind element at the time of the explosion is r = t_wv_w(t_w)∝t^{1 + b}_w so that t_w∝r^{1/(1 + b)} and we have $\dot{M}_w\propto r^{a/(1+b)}$, v_w∝r^{b/(1 + b)}, and ρ∝r^{−2 + (a − b)/(1 + b)}. For a constant wind velocity (b = 0) this gives k = 2 − a, which corresponds to k = 2 for a = 0 (constant wind mass flux) and k = 1 for a = 1 (linearly increasing mass flux with time).

A brief description of our numerical methods and initial conditions for both jet and external medium models is given in Section 2. Detailed hydrodynamic simulations of GRB jets interacting with k = 1, 2 stratified media are presented in Sections 3 and 4, where Section 3 is devoted to the jet dynamics and the resulting afterglow emission is discussed in Section 4. For completeness and comparison, the interaction with a constant-density medium (k = 0) is also discussed, although the reader is referred to De Colle et al. (2012) for a review of the current state of hydrodynamical modeling with k = 0. Our conclusions are summarized in Section 5.

2.1. Code Description and Initial Conditions

To study the dynamics of a GRB jet propagating in a stratified external medium, we carry out a set of two-dimensional simulations using the special relativistic hydrodynamic (SRHD) version of the adaptive mesh refinement (AMR) code Mezcal (De Colle et al. 2012). The Mezcal code integrates the SRHD equations by using a second-order (in space and time, except in shocks where it reduces to first order in space by a minmod limiter) upwind scheme based on the relativistic Harten, Lax, and van Leer (HLL) method (Schneider et al. 1993). The equation of state (EOS), relating enthalpy to pressure and density, is taken from Ryu et al. (2006), which approximates the exact Synge (1971) EOS with an error of 0.5%. This EOS properly recovers the correct values of the adiabatic index Γ in the ultrarelativistic (Γ = 4/3) and Newtonian (Γ = 5/3) regimes. The reader is referred to De Colle et al. (2012) for a detailed description of the code and an extensive list of numerical tests.

For the initial conditions, we use a conical wedge of half-opening angle θ₀, within which the initial radial profiles of pressure, density, and Lorentz factor in the post-shock region are taken from the spherical Blandford–McKee self-similar solutions for a stratified medium:

$$\begin{equation} \rho = A_k r^{-k} \; . \end{equation} \tag{ 1 }$$

Two-dimensional simulations with k  =  0 (homogeneous medium), k = 1, and k = 2 (corresponding to a steady stellar wind medium) are then evolved to study the lateral expansion and deceleration of the jet.

To accurately study the dynamics near the jet break time, an initial shock Lorentz factor of $\Gamma _{\rm sh,0} =\sqrt{2}\times 20$ and an initial half-opening jet angle θ₀ = 0.2 rad are selected, so that Γ_{sh, 0} ≫ θ⁻¹₀. The isotropic equivalent energy is taken to be E_iso = 10⁵³ erg, corresponding to a total jet energy content of E_jet = E_iso(1 − cos θ₀) ∼ 2 × 10⁵¹ erg. The ambient medium is assumed to have a density ρ₀ = A₀ = 1.67 × 10⁻²⁴ g cm⁻³ (for the case k = 0, which corresponds to ρ₀ = n₀m_pc² with n₀ = 1 cm⁻³) and a pressure p = ηρ₀c², with η = 10⁻¹⁰. The value of η has no bearing on the outcome of the simulation as long as the Mach number remains large, i.e., $\mathcal {M} \sim \eta ^{-1/2} v_{\rm sh}/c \gg 1$, where v_sh is the shock velocity. As the simulation continues to evolve well into the Newtonian regime, this condition can be expressed as v_sh ≫ 3 (η/10⁻¹⁰)^1/2 km s⁻¹. The density profiles in the cases k = 1, 2 are fixed here by assuming the jet break radius (in the lab frame) to be the same for all k: R_j(k) = R_j(k = 0). This can be rewritten (Blandford & McKee 1976) as

$$\begin{equation} R_{\rm j} = \left(\frac{(17-4k)E_{\rm iso}}{8\pi A_k \Gamma _{\rm j}^2 c^2} \right)^{1/(3-k)} = \left(\frac{17 E_{\rm iso}}{8\pi A_0 \Gamma _{\rm j}^2 c^2} \right)^{1/3}\, , \end{equation} \tag{ 2 }$$

where $\Gamma _{\rm j}=\sqrt{2}/\theta _0$.

From Equation (2) we have A_k = A₀R^k_j(17 − 4k)/17, so that the density of the ambient medium is given by

$$\begin{equation} \rho = \frac{17 - 4 k}{17}A_0 \left(\frac{r}{R_{\rm j}}\right)^{-k} \; , \end{equation} \tag{ 3 }$$

which guarantees R_j to remain unchanged for varying k. With this constraint, the value of the density at the jet break radius ρ(r = R_j) differs, compared with the k = 0 case, by factors of 13/17 and 9/17 for k = 1 and k = 2, respectively. In the simulations presented in this paper, R_j = 9.655 × 10¹⁷ cm, corresponding to a jet break time of t_j = R_j/c ≈ 372 days.

The case k = 2 corresponds to a steady spherically symmetric wind with $\dot{M} = 1.45 \times 10^{-5}$ (v/10³ km s⁻¹) M_☉ yr⁻¹, typical for Wolf–Rayet stars (e.g., Chiosi & Maeder 1986).

The jet is expected to begin decelerating to non-relativistic speeds at

$$\begin{equation} t_{\rm NR} \approx \frac{L_{\rm Sedov}}{c} = \left(\frac{(3-k)E_{\rm iso}}{4 \pi A_k c^2} \right)^{1/(3-k)} \; \end{equation} \tag{ 4 }$$

corresponding to t_NR ≈ 970, 3800, and 11, 000 days (in the lab frame) for k = 0, 1, 2, respectively.

The simulations with k = 0, 1 employ a spherical computational domain of radial and angular size (L_r, L_θ) = (1.1 × 10¹⁹ cm, π/2) while the simulation with k = 2 uses (L_r, L_θ) = (2.2 × 10¹⁹ cm, π/2). The inner boundaries are located at (1.8, 1.2, 0.3) × 10¹⁷ cm for k = (0, 1, 2), respectively. The AMR code uses a basic grid of (100, 6) cells in the (r, θ) directions, and 15 (k = 0, 1) or 16 (k = 2) levels of refinement, corresponding to a maximum resolution of (Δr_min, Δθ_min) = (6.71 × 10¹² cm, 1.60 × 10⁻⁵ rad). To keep the resolution of the relativistic thin shell Δ∝t^{4 − k} approximately constant, the maximum number of levels of refinement N_levels is decreased with time (De Colle et al. 2012) as N_levels = max [7, N_{levels, 0} − (4 − k)log (t/t₀)/log (2)]. The simulations are halted after 150 years. We also carried out a higher resolution simulation (for the k = 2 case) using a basic grid of (1000, 16) cells in the (r, θ) directions and 14 levels of refinement. The light curves computed from this simulation are very similar to those obtained from the lower resolution run, implying that convergence has been achieved.

The Mezcal code is parallelized using the "Message Passing Interface" library, enabling the highest resolution simulation to be run in about two weeks on a local supercomputer with 160 processors and the low resolution in about a quarter of that time.

2.2. Afterglow Radiation

To compute the afterglow radiation, we use the method described in De Colle et al. (2012). As the main goal of the current calculations is to study the effect of the jet dynamics on the afterglow light curves, a simple model is employed to calculate the emanating radiation. We assume synchrotron to be the primary emitting mechanism, while ignoring self-absorption and inverse Compton scattering. The self-absorption frequency, in particular, is given by v_sa = 2.04 × 10⁹ Hz and v_sa = 8.23 × 10⁹t^−3/5_obs Hz for k = 0 and k = 2, respectively, for our choice of the parameters (e.g., Granot & Sari 2002). A multi-dimensional ray tracing code, necessary to properly handle the self-absorption process, is currently under development and will be presented elsewhere.

The microphysics processes responsible for field amplification and particle acceleration are parameterized here by assuming that the magnetic field everywhere in the shocked region holds a fraction _B = 0.1 of the local internal energy density in the flow, while the non-thermal electrons just behind the shock hold a fraction _e = 0.1 of the internal energy, and have a power-law energy distribution, N(γ_e)∝γ^−p_e, with p = 2.5. We also assume the source to be at a redshift of z = 1, corresponding to a luminosity distance of d_L = 2.05 × 10²⁸ cm. The local emissivity $P^{\prime }_{\nu ^{\prime }}$ is taken to be a broken power law,

$$\begin{equation} \frac{P^{\prime }_{\nu ^{\prime }}}{P^{\prime }_{\nu ^{\prime },{\rm max}}} = \left\lbrace \begin{array}{@{}ll}(\nu ^{\prime }/\nu ^{\prime }_m)^{1/3} & \nu ^{\prime }<\nu ^{\prime }_m<\nu ^{\prime }_c \, , \\ \ms (\nu ^{\prime }/\nu ^{\prime }_c)^{1/3} & \nu ^{\prime }<\nu ^{\prime }_c<\nu ^{\prime }_m \, , \\ \ms (\nu ^{\prime }/\nu ^{\prime }_m)^{(1-p)/2} & \nu ^{\prime }_m<\nu ^{\prime }<\nu ^{\prime }_c \, , \\ \ms (\nu ^{\prime }/\nu ^{\prime }_c)^{-1/2} & \nu ^{\prime }_c<\nu ^{\prime }<\nu ^{\prime }_m \, , \\ \ms (\nu ^{\prime }/\nu ^{\prime }_m)^{(1-p)/2}(\nu ^{\prime }/\nu ^{\prime }_c)^{-1/2} & \nu ^{\prime }>\max (\nu ^{\prime }_m,\nu ^{\prime }_c)\, , \end{array}\right. \end{equation} \tag{ 5 }$$

with the following flux normalization and break frequencies,

$$\begin{equation} P^{\prime }_{\nu ^{\prime },{\rm max}} = 0.88\,\frac{512\sqrt{2\pi }}{27}\left(\frac{p-1}{3p-1}\right) \frac{q_e^3}{m_ec^2}(\epsilon _Be^{\prime })^{1/2}n^{\prime }_e \, , \end{equation} \tag{ 6 }$$

$$\begin{equation} \nu ^{\prime }_m = \frac{3\sqrt{2\pi }}{8}\left(\frac{p-2}{p-1}\right)^2\frac{q_e}{m_e^2c^5} \epsilon _B^{1/2}\epsilon _e^2(e^{\prime })^{5/2}(n^{\prime }_e)^{-2} \, , \end{equation} \tag{ 7 }$$

$$\begin{equation} \nu ^{\prime }_c = \frac{27\sqrt{2\pi }}{128}\frac{q_e m_e c}{\sigma _T^2} (\epsilon _Be^{\prime })^{-3/2}\left(\frac{\Gamma }{t_z}\right)^2 \, . \end{equation} \tag{ 8 }$$

In our simple prescription for electron cooling, which is similar to the one used by Granot et al. (2001) and Zhang & MacFadyen (2009), the electrons are assumed to have cooled at their current local cooling rate over the dynamical time, which is in turn approximated as t'_dyn ≈ t_z/Γ, so that the expression in Equation (8) is simply derived from

$$\begin{equation} \gamma _c = \frac{3m_e c}{4\sigma _T\epsilon _Be^{\prime }t^{\prime }_{\rm dyn}} \approx \frac{3m_e c\Gamma }{4\sigma _T\epsilon _Be^{\prime }t_z}, \quad \quad \nu ^{\prime }_c = \frac{3q_e B^{\prime }\gamma _c^2}{16m_e c}\, . \end{equation} \tag{ 9 }$$

In the simulations, 1000 outputs (i.e., "snapshots" of the dynamics), spaced logarithmically in time, are saved. At each snapshot the values of the hydrodynamic variables are provided at the center of a computational cell. Once the emissivity is computed for a particular cell, the flux (in the observed frame) is assigned to a particular observed time and position on the sky (see De Colle et al. 2012 for more details).

The intrinsic limitations of the radiation method are described in detail in De Colle et al. (2012). In brief, the light curve converges at all frequencies and times except for short t_obs (corresponding to regions with shock Lorentz factor ≳ 10) where neither the very high resolution used in this paper is sufficient to fully resolve the thin ultrarelativistic post-shock region. Furthermore, the differences between our treatment of the electron cooling and the results presented by Granot & Sari (2002), together with tests of the radiation code, are described in detail in De Colle et al. (2012). The simulation with k = 0, in particular, gives afterglow light curves that are nearly identical to those computed by Zhang & MacFadyen (2009).

In addition to the contributions to the afterglow radiation computed by post-processing the results of the hydrodynamics simulations, contributions from earlier lab frame times are included, corresponding to the blast-wave decelerating from $\Gamma _1 = \Gamma (\chi = 1) =\Gamma _{\rm sh}/\sqrt{2} = 200$ to Γ₁ = 20. Here χ(r/R_sh) = 1 + 2(4 − k)Γ²_sh(1 − r/R_sh) is a self-similar variable that quantifies the distance from the shock front (Blandford & McKee 1976). These are computed using the same conical wedge taken out of the Blandford–McKee self-similar solution that is used for initializing our simulations. The mapping of the Blandford–McKee solution is implemented by using a high-resolution grid, starting at the position of the shock front (which varies with time) and sampling the Blandford–McKee solution at intervals of fixed ΔΓ = 0.01. The values of the proper density ρ, internal energy density e_int, four-velocity u, and self-similar variable χ replace those coming from the simulations and are taken from the Blandford–McKee self-similar solution at the relevant lab frame time. In order to calculate the contributions to the observed radiation, the mapped jet radial structure is subsequently integrated over all angles (0 ⩽ θ ⩽ θ₀; 0 ⩽ ϕ ⩽ 2π). This procedure provides a reasonable description of the afterglow radiation at earlier times and it is significantly more accurate than ignoring the contributions from lab frame times preceding the start of the simulation.

Detailed hydrodynamic simulations of the evolution of a GRB jet in a stratified medium with k = 0, 1, 2 are presented in Figures 1 and 2 where the density and velocity contours of the expanding ejecta at various times are plotted. A transient phase caused by the sharp lateral discontinuity in the initial conditions is observed in all cases as the shock expands laterally and a rarefaction front moves toward the jet axis. This initial phase, during which shearing instabilities are observed to be prominent at the contact discontinuity (separating the original Blandford–McKee wedge material and the later shocked external medium), lasts for about a dynamical timescale and is followed by the establishment of an egg-like bow shock structure that persists throughout the simulations. The velocity quadrivector (Figure 2) shows strong stratification in the θ direction. The expansion velocity of the jet remains mainly radial at most angles, with a non-relativistic angular component being prominent at large angles. The substructures seen in the velocity quadrivector along the z-axis at late times (generated by the convergence of turbulent flow) carry a small fraction of the energy and have a negligible effect on the light curves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Similar resulting bow shock structures are observed for k = 0, 1, and 2. However, because the rate at which mass is swept-up is larger for smaller values of k, the bow shock lateral expansion augments with increasing k. As clearly seen in Figure 3, the ratio between the bow shock width and height as the ejecta expand changes with k. This can be understood as follows. Small values of k correspond to a larger increase in the swept-up external mass and larger decrease in the Lorenz factor. For the spherical case, in particular, M(< R)∝R^{3 − k} and Γ∝R^{−(3 − k)/2}, and the same trend should persist for the non-spherical case.

Zoom In Zoom Out Reset image size

The velocity quadrivector, u = Γv/c = Γβ, of the expanding jets are shown in Figure 4 for three different angle-integrated quantities: mass, energy, and emissivity. The mean value of u is larger when weighted over the energy or emissivity than over the shocked rest mass until t ≲ 10 × t_j. This clearly illustrates, in agreement with previous analytical and numerical results limited to the case k = 0 (Granot et al. 2001; Zhang & MacFadyen 2009), that during the relativistic phase, most of the shocked rest mass resides in relatively slow material at the edges of the jet, while most of the energy is stored in the fastest moving material near the head of the jet.

Zoom In Zoom Out Reset image size

As illustrated in Figure 4, the Blandford–McKee and Sedov–Taylor self-similar solutions fail to provide an adequate description of the jet dynamics at t_j ≲ t ≲ t_NR(E_iso) with the disagreement becoming less pronounced before t_j and after t_NR(E_iso). Between these two limiting cases, −dlog u/dlog t evolves at early and late times between the two asymptotic slope values, as seen in the bottom panel of Figure 4. The evolution of −dlog u/dlog t is, however, non-monotonic as it first increases above (3 − k)/2 and only then decreases down to (3 − k)/(5 − k). This behavior is mainly caused by the faster decrease in Γ compared to a spherical flow at t > t_j due to the lateral expansion of the jet. It also relates to the fact that the Blandford–McKee solution depends on E_iso while the corresponding Sedov–Taylor solution uses the jet's true energy, E_jet and, as a result, the ratio of u(t) for these two limiting cases is ∼θ⁻¹₀ at t = t_NR(E_jet) and ∼θ^{−2/(5 − k)}₀ at t = t_NR(E_iso) ∼ θ^{−2/(3 − k)}₀t_NR(E_jet).

Figures 5 and 6 show the resulting R_⊥(t), R_∥(t), and θ_j(t) for k = 0, 1, 2 and different recipes for estimating the transverse, parallel, and angular size scales within the jet (e.g., when averaged over mass, energy, and emissivity). For all values of k the early lateral spreading of the jet, which starts around t ∼ t_j, is observed to initially involve only a modest fraction of the total energy, with the bulk of the energy reaching angles well above θ₀ at significantly later times.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For k = 0, previous numerical simulations and analytical models assuming a small lateral expansion for t ∼ t_NR (e.g., Granot et al. 2005) have shown that spherical symmetry is approached on timescales much larger than t_NR. In particular, Figure 5 shows that the growth of R_∥ is essentially stalled at t ∼ t_NR while R_⊥ continues to grow as the flow gradually approaches spherical symmetry. This effect is less pronounced for increasing k, since R_∥ continues to increase even after t_NR(E_iso), albeit more slowly. This contributes to the faster growth in θ_j for lower k-values at late times, contrary to the opposite situation at early times (t ≲ t_j). This causes GRB jets expanding into steeper density profiles to approach spherical symmetry at progressively later times as argued by Ramirez-Ruiz & MacFadyen (2010) for k = 2.

Since the rate of lateral spreading of the jet increases as Γ decreases (see, e.g., Equation (2) of Granot 2007) and Γ(R_j) = θ⁻¹₀ is the same for all k, then the jet lateral spreading is expected to increase with k for R ≲ R_j (where Γ(R) decreases with k for a given R), while the opposite should hold for R ≳ R_j (where, for a given R, Γ(R) increases with k). Such a behavior is also seen in analytic models (Granot 2007; Granot & Piran 2012).

Figure 5 also plots the temporal evolution of ΓR_⊥/R_∥ ≈ Γθ_j, which is observed to approach unity at t ≫ t_j. This should be compared with the results of semi-analytic models (e.g., Rhoads 1997; Sari et al. 1999; Kumar & Panaitescu 2000). These models predict Γθ_j ≈ 1 at t ≳ t_j, and Γ to decrease rapidly with lab frame time t, which is not observed here. In the simulations, Γ decreases rather slowly with t (as a power law). The jet angular size θ_j (see Figure 6), on the other hand, is observed to increase only logarithmically with t for all k until the flow becomes non-relativistic.

As shown in Figure 6, the weighted mean of θ_j over the emissivity (and to a slightly lesser extent over the energy) remains practically constant until t/t_j ∼ a few, while the weighted mean over the shocked rest mass is significantly larger, in accord with earlier results (Granot et al. 2001; Piran & Granot 2001; Zhang & MacFadyen 2009). This indicates that, as argued before, a large fraction of the swept-up external rest mass is concentrated at the edges of the jet, while most of the energy and emission lies near the head. Moreover, it implies that (as discussed above and seen in the temporal evolution of δ depicted at the bottom of Figure 5) the lateral expansion at early times, t ≲ t_j, is significantly faster for larger values of k, while the situation is reversed at late times.

Figure 7 plots the temporal evolution of the energy (excluding rest energy) per solid angle, = dE/dΩ, as a function of the angle θ from the jet symmetry axis, for k = 0, 1, 2. At t ≳ 50 yr the energy distribution appears nearly spherical for all k's. At earlier times, a clear k-dependence trend is observed, where the energy spreads to larger solid angles faster for a more stratified medium, but a correlation is less evident when one compares (θ) for different k-values at the same four velocity u rather than the same lab frame time t.

Zoom In Zoom Out Reset image size

Abundant confirmation is provided here that the dynamics of GRB jets are greatly modified by the radial profile of the surrounding circumburst density. Most analytic formalisms (e.g., Rhoads 1999) derive an exponential lateral spreading with lab frame time or radius at t > t_j, which ultimately erases all information about the initial jet opening angle and relies solely on the true energy content of the jet: E_jet. No exponential lateral expansion is observed in our study for k = 1, 2, consistent with previous numerical work for expansion in a constant density medium (Granot et al. 2001; Zhang & MacFadyen 2009; van Eerten & MacFadyen 2011). As illustrated in Figure 6, the evolution of the jet's angular scale containing a constant fraction of the total energy is logarithmic and is not self-similar as it retains memory of the initial jet opening angle. The deviation from the expected self-similar exponential lateral expansion behavior (Gruzinov 2007) might be at least partly due to u rapidly decreasing with the polar angle θ from the jet symmetry axis, so that the flow is no longer ultrarelativistic (u ≫ 1) as it has been previously assumed. Even with the expectation that such a self-similar solution would be only very slowly attained (Gruzinov 2007), the maximal Lorentz factor at the head of the jet in this formalism is predicted to decrease exponentially with time, which appears to be inconsistent with our numerical results.

The resolution of this apparent inconsistency between analytic models and numerical simulations can be attributed to the modest values of θ₀ used in the simulations, which result in the breakdown of the analytic models that assume Γ ≫ 1 and θ_j ≪ 1 soon after the jet starts spreading sideways Γ < θ⁻¹₀ and before it can reach a phase of exponential lateral expansion (Wygoda et al. 2011; Granot & Piran 2012). In the small region in which the analytical models are valid, 1 ≪ Γ < θ⁻¹₀, there is reasonable agreement with simulation results (Wygoda et al. 2011). A generalization of these analytic models to any values of Γ or θ_j (Granot & Piran 2012) shows reasonable agreement with the results of simulations from the early ultrarelativistic stage to the late Newtonian stage. Such generalized analytic models predict that if the jet is initially extremely narrow then there should still be an early phase of exponential lateral spreading. However, these models make the simplifying approximation of a uniform jet, while in practice u quickly drops with θ. This causes a breakdown of the u ≫ 1 assumption used to derive the self-similar solution, which is only slowly attained even under ideal conditions (Gruzinov 2007).

Figure 8 shows the emerging light curves at frequencies ranging from the radio to gamma-rays (ν = 10⁹, 10¹¹, 10¹³, 10¹⁵, 10¹⁷, 10¹⁹ Hz), and corresponding spectra at different observed times t_obs, for k = 0, 1, 2, including the effects of electron cooling and the contribution from a mapped Blandford–McKee solution (with 20 ⩽ max (Γ) ⩽ 200). Figure 9 shows the light curves computed for ν = 10⁹, 10¹³, and 10¹⁷ Hz, for the two-dimensional simulation and the Blandford–McKee conical wedge, as in Figure 8, but illustrating the contributions to the light curve arising from the various evolutionary stages of the blast wave, quantified here by considering the emission from lab frame times where $\Gamma _{\rm sh}(t)/\sqrt{2}$, given by the Blandford–McKee solution, ranges between 10 and 20, 5 and 10, 2 and 5, 1 and 2, respectively. As expected, lower Lorentz factors contribute to the observed flux at later times. A slightly more subtle effect is that at the same observed time the flux at low frequencies comes from slightly later lab frame times t (corresponding to a lower Blandford–McKee Lorentz factor Γ_sh(t)). This is because there is a lower flux contribution from the sides of the jet compared with the center, as reflected by the fact that the afterglow image is more limb brightened at higher frequencies and less so at lower frequencies (Granot et al. 1999; Granot & Loeb 2001; Granot 2008), resulting in a smaller typical angular delay time (t_θ = R(t) ≈ Rθ²/2c) in the arrival of photons to the observer (which is along the jet axis in these figures). As a result, the flux at the same observed time t_obs is dominated by larger lab frame times t.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The spectra at different observer times are shown in Figure 8. For all values of k, the spectra evolves from a fast cooling (with ν_c < ν_m and F_ν∝ν^1/3, ν^−1/2, and ν^−p/2 for ν < ν_c, ν_c < ν < ν_m, and ν > ν_c, respectively) to a slow cooling regime (with ν_m < ν_c and F_ν∝ν^1/3, ν^{(1 − p)/2}, and ν^−p/2 for ν < ν_m, ν_m < ν < ν_c, and ν > ν_m, respectively). The characteristic frequency ν_m quickly drops with time with an asymptotic slope of −2.9, −2.6, −2 for k = 0, 1, 2, respectively (while one expects ν_m∝t^{−(15 − 4k)/(5 − k)}, which is relatively closed to our result), while ν_c increases at late times as ν_c∝t, that is, with a slope independent on the particular stratification of the ambient medium (for comparison, in the Sedov–Taylor regime one expects ν_c∝t^{(2k − 1)/(5 − k)}).

As shown in Section 3, a jet moving in a stratified medium (with k = 1 and k = 2) decelerates to sub-relativistic speed over larger distances with respect to a jet moving in a homogeneous medium (k = 0). The consequences of it on the light curve are particularly evident at radio frequencies (Figure 9), where the contribution from mildly and sub-relativistic material is negligible in the k = 2 case and dominant in the k = 0 up to t ∼ 10³ days.

Figures 8 and 9 show a pan-chromatic dip or flattening in the light curves at around half a day for k = 0, a third of a day for k = 1, and significantly earlier for k = 2. This feature is also seen in Figure 10, which shows the temporal index α ≡ −dlog F_ν/dlog t_obs as a function of t_obs, where the earliest value of α is larger than expected analytically for a spherical flow (or for a jet viewed along its axis, before the jet break time). Figure 9 clearly illustrates the reason for this behavior. It basically occurs at the point where the dominant contribution to the observed flux switches from the Blandford–McKee wedge with 20 ⩽ Γ_sh(t) ⩽ 200 to the simulation, which corresponds to later lab-frame times. As pointed out and calculated in De Colle et al. (2012) for the spherical case, the relaxation of the mapping of the analytic Blandford–McKee self-similar solution to the numerical solution and the finite resolution of the simulation result in a dip in the Lorentz factor that is gradually recovered as the shocked region becomes wider and thus better resolved with time. This produces a dip in the light curve that gradually goes away as the resolution of the simulation is increased (see Figures 5–7 of De Colle et al. 2012). This feature is a numerical artifact of the finite resolution of the simulation. Similar errors in the light curves were also present in previous simulations for the k = 0 case (e.g., our light curve in the case k = 0 is nearly identical to that by Zhang & MacFadyen 2009 as depicted in De Colle et al. 2012).

Zoom In Zoom Out Reset image size

A smaller contribution (although not easily quantifiable) to the pan-chromatic dip in the light curve is due to the particular initial conditions chosen in this paper. In fact, as the jet initially has sharp edges (a step function in the θ-direction), once the simulation starts there is a relaxation period occurring in the lateral direction on a dynamical timescale (as a rarefaction wave propagates from the edge of the jet toward its center). This lateral transient phase triggered by the sharp-edged jet is also imprinted in the light curves around the time of the dip or flattening, and, contrary to the limited resolution artifacts, is not expected to go away as the resolution is increased. This artifact might be less pronounced for initial conditions that are smoother in the lateral direction (e.g., a jet with an initial Gaussian angular profile).

Apart from this early-time, artificial feature, there is the expected pan-chromatic jet break that is present at all frequencies above ν_m and is observed between a day for k = 2 to several days for k = 0. These jet break features are discussed in more detail below.

4.1. Jet Breaks

Figure 10 plots the shape of the jet break, i.e., the temporal decay index of the light curve, α ≡ −dlog F_ν/dlog t_obs, as a function of observer time, t_obs, for different observed frequencies and k-values. We shall first discuss the pan-chromatic jet break features at frequencies that are above the typical synchrotron frequency at the time of the jet break, ν > ν_m(t_{obs, j}). As shown in Figure 10, the temporal decay of the light curve becomes smoother for increasing k, as derived in analytic models (Kumar & Panaitescu 2000, hereafter KP00). However, the steepening in the light curve occurs within a significantly smaller observed time period than that predicted by analytic models. Most of the increase in α occurs over a factor of ≈3–5 in time for k = 0 (compared with a decade in time predicted in KP00) and within about one decade in time for k = 2 (compared with four decades in time predicted in KP00). The relatively sharper jet break (compared with analytic expectations) in a stratified medium may permit the detection of such a jet break. We also note that there is an "overshoot" in the value of the temporal decay index α just after the jet break, which is more prominent for lower k-values (in agreement with previous results; Granot 2007). After this overshoot α gradually decreases, and there is also a noticeable curvature in the light curve as the flow becomes mildly relativistic and eventually approaches the Newtonian regime. The effects of electron cooling on the shape of the jet break appear to be rather modest in most cases.

At low frequencies, ν < ν_m(t_{obs, j}) (see Figure 10, upper panel), there is only a very modest increase in α near t_{obs, j}. On the other hand, when the break frequency ν_m sweeps past the observed frequency ν, a very sharp break is seen (i.e., increase in α). Both features are present for k = 0, and we find here that they also persist for higher k-values. Moreover, we also find that this break is sharper for smaller k-values. This is because the corresponding spectral break (at ν_m) is very sharp for our simple broken power-law spectral emissivity model and is not degraded by the contribution from multiple parts of the jet at smaller k-values (in addition ν_m decreases somewhat faster in time at t_obs > t_{obs, j} for smaller k-values). We expect that a more realistic synchrotron emissivity function would result in a significantly smoother spectral break at ν_m, which would in turn lead to a correspondingly smoother temporal break.

Figure 11 shows afterglow light curves for three different observed frequencies (ν = 10⁹, 10¹³, 10¹⁷ Hz; top to bottom panels), external density profiles (k = 0, 1, 2), and viewing angles (θ_obs/θ₀ = 0, 0.5, 1), both with and without electron cooling (left and right panels, respectively). Figure 12 shows the corresponding values of the temporal decay index α for ν = 10¹⁷ Hz. The dependence of the light curves as a function of viewing angles is in qualitative agreement with the results of van Eerten et al. (2010a) in the k = 0 case. For stratified media, Figures 11 and 12 show that the shape of the jet break is predominantly regulated by the change in viewing angle (within the initial jet aperture, 0 ⩽ θ_obs/θ₀ ⩽ 1) rather than by the external density power-law index k (in the range 0 ⩽ k ⩽ 2). For θ_obs = 0 most of the steepening occurs within a factor of ∼2–4 in time for k = 0, 1, 2 while for θ_obs/θ₀ ∼ 0.5–1 it takes ∼1–2 decades for k = 0, 1, 2. This is particularly interesting because previous analytical work have argued that the effect of varying k should be significantly larger. It can also be seen in Figure 12 that the jet induced steepening starts earlier and ends later for larger k-values and for larger viewing angles (or θ_obs/θ₀ values). Also, the overshoot in the value of α is larger for greater k-values or θ_obs/θ₀ values. The jet break time is also observed to occur later for larger viewing angles at all values of k and varies over a factor of ∼3–5 for 0 ⩽ θ_obs/θ₀ ⩽ 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The change in the jet break duration with k is due to the slower evolution of Γ with t or R ≈ ct as well as t_obs for larger k-values (Γ∝R^{(k − 3)/2}∝t^{(k − 3)/(8 − 2k)}_obs for a spherical flow). For θ_obs = 0 and ν > ν_m(t_{obs, j}), the jet break duration roughly corresponds to the time it takes the beaming cone to grow past the limb-brightened outer part of the image. If crudely neglecting lateral spreading (since most of the emission near the jet break time is from within the initial jet aperture; Piran & Granot 2001), so that the dominant effect is the "missing emission" from outside the edges of the jet (Granot 2007), and requiring that the beaming cone (of angle θ ≲ 1/Γ around the line of sight) grows by a factor of f_k, then this would correspond to a factor of ∼f^{(8 − 2k)/(3 − k)}_k in observed time. However, the resulting image is more limb brightened for smaller k-values (Granot & Loeb 2001; Granot 2008), and, as a result, one might estimate f_{k = 0} ∼ 1.3, f_{k = 1} ∼ 1.4, and f_{k = 2} ∼ 1.5, which would result in factors of ∼2, ∼3, and ∼5 in the observed time, respectively, in rough agreement with our numerical results.

As to the effect of the viewing angle for a fixed value of k, the addition to the duration of the jet break relative to θ_obs = 0 corresponds approximately to the time it takes the edge of the beaming cone (1/Γ) to grow from θ₀ to θ₀ + θ_obs. Thus, for θ_obs = θ₀ this corresponds to a factor of 2 decrease in Γ or a factor of ∼2^{(8 − 2k)/(3 − k)} increase in the observed time (i.e., factors of ∼6, ∼8, and ∼16 for k = 0, 1, and 2, respectively). This is in rough agreement with our numerical results. According to this simple estimate, the duration of the jet break for θ_obs = θ₀ and k = 2 should be a factor of ∼(2f₂)^{(8 − 2k)/(3 − k)} ∼ 3⁴ ∼ 81 in time, or almost two decades in observed time, also in agreement with the results of our calculations.

4.2. Radio Calorimetry

Figure 13 shows the radio light curves (at ν = 10⁹ Hz) for k = 0, 1, 2 from our two-dimensional numerical simulations of a double-sided jet, as well as for a spherical blast wave with the same true energy and a double-sided cone of fixed half-opening angle θ₀ calculated from a spherical blast wave with the same isotropic equivalent energy (where θ_obs = 0 in the two non-spherical cases), and the initial time of the simulation is set by the Blandford & McKee (1976) self-similar solution:

$$\begin{equation} t \cong \frac{R}{c} \cong \frac{1}{c} \left[\frac{(17-4k)E}{8 \pi A_k c^2 \Gamma _{\rm sh}^2} \right]^{1/(3-k)} \, . \end{equation} \tag{ 10 }$$

As expected, the light curves computed from a spherical blast wave with the same isotropic energy and from the two-dimensional simulation match reasonably well at early times. For the double-sided jet, it can be seen that the bump in the light curve near the non-relativistic transition time, caused because the counterjet (whose contribution is indicated by a dashed line) becomes visible, is much more prominent for low values of k and becomes significantly more modest for larger k-values. This effect is caused by the more gradual deceleration of the jet for larger k-values (as the same mass of external medium is swept-up over a larger range in radii), which causes the counterjet to become visible more gradually, resulting in a wider, lower peak flux bump. In particular, for k = 2 it amounts to a fairly modest and rather slow flattening of the light curve, which might be hard to discern observationally. This might, however, not help explain the lack of a clear flattening or rebrightening in the late radio afterglow of GRB 030329 (e.g., Pihlström et al. 2007), since in that case detailed afterglow modeling favors a uniform external density (k = 0; van der Horst et al. 2008).

Zoom In Zoom Out Reset image size

Comparison of the radio flux at late times from a double-sided jet and from a spherical blast wave with the same true energy near the non-relativistic transition time shows that they are broadly similar but may differ by up to a factor ≲ 3. For k = 0 and k = 1 the spherical analog slightly overpredicts the flux before the contribution from the counterjet becomes important, and underpredicts the flux once the emission from the counterjet becomes dominant, while for k = 2 the spherical analog consistently underpredicts the flux by up to a factor of ≲ 3. This may result in a small but not negligible error in the estimation of the true energy in the double-side jet assuming a spherical sub-relativistic flow, as is commonly done in radio calorimetry studies (Kaneko et al. 2007; Berger et al. 2003; Frail et al. 2005; Gorosabel et al. 2006; Kulkarni et al. 1998), both overestimating or underestimating the real true energy depending on the stratification of the ambient medium and the observer time.

We have studied the dynamics of GRB jets during the afterglow stage as they propagate into different external density profiles, ρ_ext = Ar^−k for k = 0, 1, 2, using detailed hydrodynamic simulations. Our main results, which relate both to the dynamics and the resulting afterglow emission, can be summarized as follows.

For the same initial half-opening angle θ₀ and external density at the jet break radius (which is defined by Γ₁(R_j) = θ⁻¹₀), the lateral spreading is initially (at R < R_j) larger for higher k-values. This arises because at the same radius (or lab frame time) the typical Lorentz factor is lower. At late times (R > R_j) the situation is reversed, and the effective jet opening angle at a fixed lab frame time is similar for different k-values. Since for higher k-values a larger range of radii is required in order to sweep-up the same amount of mass, the whole evolution extends over a much wider range of radii and times. As a result, the jet break in the afterglow light curve is smoother and more gradual, the non-relativistic transition occurs later, and the flow approaches spherical symmetry more slowly and over longer timescales. The effective jet opening angle is observed to increase only logarithmically with lab frame time (or radius) once the jet comes into lateral causal contact (i.e., when Γ drops below θ⁻¹₀).

As long as the jet is relativistic, most of the energy and emission are concentrated near the head of the jet while the slower material at the edges carries relatively little energy (even though it carries a substantial fraction of the swept-up rest mass). This holds true for all k-values. Once the jet becomes sub-relativistic, at t > t_NR(E_iso), it quickly spreads laterally and swiftly starts to approach spherical symmetry. The energy weighted mean value of u(t) is observed to be of order unity at t/t_j ∼ 2 rather than at t ∼ t_NR(E_iso), as one might naively expect. We find that there is little k-dependence on the temporal evolution of θ_j, so that irrespective of the external medium radial profile, all of the expanding jets approach spherical symmetry at similar times (∼1–1.5 decades after t_j). A similar conclusion can be reached from the calculated evolution of R_∥/R_⊥ with t/t_j.

We find that contrary to the expectations of analytic models, the shape of the jet break is affected more by the viewing angle (within the initial jet aperture, 0 ⩽ θ_obs/θ₀ ⩽ 1) than by the steepness of the external density profile (for 0 ⩽ k ⩽ 2). Larger viewing angles result in a later jet break time and a smoother jet break, extending over a wide range in time, and with a larger overshoot (initial increase in the temporal decay index α beyond its asymptotic value), which is observed to be more prominent for lower k-values. Larger k-values result in more gradual jet breaks, but the sharpness of the jet break is affected even slightly more by the viewing angle as argued above. The counterjet becomes visible around t_NR, and for k = 0 this results in a clear bump in the light curve. However, for larger k-values the jet deceleration is more gradual and as a result a wider and lower bump is produced, which becomes hard to detect for k = 2, where it reduces to a mild flattening in the light curve. This may explain the lack of a clear counterjet signature in some late time radio afterglow light curves of long duration GRBs, although the dynamical complexity of their surrounding circumburst medium seriously limits the validity of a non-evolving power-law density profile (e.g., Ramirez-Ruiz et al. 2005).

Finally, we showed that the use of a spherical blast wave for estimating the total energy of the jet, as is commonly done in radio calorimetry studies, results in an error in the estimation of the true energy content of the jet that depends on the stratification of the ambient medium (being on average larger for k = 2). In particular, in the case k = 2, the spherical blast wave analogy consistently overestimates the true energy, while for the cases k = 0 and k = 1 it produces and underestimate or an overestimate depending on whether the estimation of the jet energy is done before or after the non-relativistic transition time.

We are grateful to A. MacFadyen, W. Lee, and W. Zhang for discussions. This research was supported by the David and Lucille Packard Foundation (E.R.R. and F.D.C.), the NSF (E.R.R.) (AST-0847563), the E.R.C. advanced research grant "GRBs" and a DGAPA postdoctoral grant from UNAM (D.L.C.). We acknowledge the support by S. Dong for administrating the Pleiades supercomputer, maintained and operated by the University of California at Santa Cruz, where the numerical calculations in this paper were performed.