THE HYDRODYNAMICS OF GAMMA-RAY BURST REMNANTS

Relativistic jets are common in the astrophysical environment. Objects known or suspected to produce them include radio galaxies and quasars (Begelman et al. 1984), microquasars (Mirabel & Rodriguez 1999), and gamma-ray bursts (GRBs; Gorosabel et al. 2006). An important difference between jets of GRBs and the better studied radio jets of quasars or microquasars is that active quasars often inject energy over extended periods of time into the jet while GRB sources are impulsive. Although quasar jets remain highly collimated throughout their lifetimes, GRB jets decelerate and expand significantly once they become non-relativistic. Expansion into a uniform medium has been well studied (Ayal & Piran 2001), but the interaction of a GRB remnant with a nonuniform medium remains poorly understood.

Much of our effort in this paper is therefore dedicated to determining how the morphology and dynamics of young GRB remnants is modified by the character of the circumburst medium. Some of the questions at the forefront of attention include the effects of the external medium and the degree to which GRB remnant dynamics and structures are modified by the presence of an accompanying supernova (SN). We address both of these issues here. Because massive stars are expected to have their close-in surroundings modified by the progenitor winds, we consider both free winds and shocked winds as possible surrounding media for the GRB remnant evolution. Detailed hydrodynamic simulations of this interaction are presented in Sections 4 and 6, while a brief description of the numerical methods and the initial models is given in Section 2. For completeness, the interaction with a constant density medium is discussed in Section 3. The role of SN explosions in shaping the evolution and morphology of GRB remnants is discussed in Section 5. The effects of a nonspherical circumburst medium are briefly addressed in Section 4. Discussion and conclusions are presented in Section 7.

2.1. The Underlying Dynamics

For simplicity, let us consider a uniform GRB jet with sharp edges and a half-opening angle θ_θ, with an initial value of θ₀ (other jet structures are discussed in Section 2.2). The typical angular size, R_⊥, of the jet at t < t_θ, where t_θ is the time at which the jet's Lorentz factor γ drops to θ⁻¹₀, is the same as for a spherical flow:

$$\begin{equation} R_\perp (t) \propto \left({E_{\rm iso} \over A}\right)^{1 \over 2(4-k)}t^{{5-k \over 2(4-k)}}\propto \left({E \over A}\right)^{{1\over 2(3-k)}}t_{\rm j}^{{-1\over 2(4-k)}}t^{{5-k \over 2(4-k)}} \end{equation} \tag{ 1 }$$

for an external density profile ρ_ext = Ar^−k (Granot et al. 2005). Here, E is the true kinetic energy content of the jet, and E_iso = f⁻¹_bE is the isotropic equivalent energy where f_b ≈ θ²₀/2 is the beaming factor. At these early times, the flow is described by the Blandford & McKee (1976) self-similar solution, which provides an accurate expression for the temporal evolution of the observed image size:

$$\begin{equation} R_\perp (t) = 4\times 10^{16}\left({E_{\rm iso,52}\over n_0}\right)^{1/8}\left({t_{\rm days} \over 1+z}\right)^{5/8}\;{\rm cm} \end{equation} \tag{ 2 }$$

for a constant density medium (k = 0) with n = 1 n₀ cm⁻³ (ρ = 1.67 × 10⁻²⁴ρ₀ g cm⁻³), and

$$\begin{equation} R_\perp (t) =2.5\times 10^{16}\left({E_{\rm iso,52} \over A_\ast }\right)^{1/4}\left({t_{\rm days} \over 1+z}\right)^{3/4}\;{\rm cm} \end{equation} \tag{ 3 }$$

for a stellar wind environment (k = 2) with A_* = A/5 × 10¹¹ g cm⁻¹.

Notice that even in the most optimistic cases, the characteristic size of a GRB image is only of order 1 μ about a day after the GRB at the Hubble distance (Granot & Loeb 2001), and so it cannot be resolved by existing telescopes. Obviously, the challenge is made easier for nearby sources, where the late size of the observable remnant could extend over a wide, possibly resolvable angle (Oren et al. 2004; Granot & Ramirez-Ruiz 2004; Granot et al. 2005). At late times t>t_NR, where t_NR is the non-relativistic transition time, the jet is expected to gradually approach the Sedov–Taylor self-similar solution, asymptotically reaching R_⊥(t) ∝ (E/A)^1/(5−k)t^2/(5−k). At t_θ < t < t_NR there is, however, a large uncertainty in the hydrodynamic evolution of the jet, and in particular in its rate of sideways expansion.

To illustrate the importance of sideways expansion on the evolution of the observed image, let us consider two extreme assumptions: (1) relativistic lateral expansion in the comoving frame (Rhoads 1999; Sari et al. 1999) for which θ_θ ≈ max(θ₀, γ⁻¹) so that at t_θ < t < t_NR we have γ ≈ θ⁻¹_θ ≈ θ⁻¹₀exp(−R/R_θ), and (2) little or no lateral expansion, θ_θ ≈ θ₀ for t < t_NR, in which case appreciable lateral expansion occurs only when the jet becomes sub-relativistic and gradually approaches spherical symmetry.

For relativistic lateral expansion, θ_θ ∼ 1 at t_NR = t_NR(E), where R_NR(E) = ct_NR(E) = [(3 − k)E/4πAc²]^1/(3−k) and the jet radius will be similar to that of the Sedov–Taylor solution, R_ST(E, t) ∼ (Et²/A)^1/(5−k). In this case, one expects the flow to approach spherical symmetry only after a few dynamical times. This is probably not the case, as clearly illustrated by the results of numerical simulations (Granot et al. 2001; Cannizzo et al. 2004; Granot 2007) showing only modest lateral expansion as long as the jet is relativistic. If lateral expansion is neglected, the jet becomes sub-relativistic only at

$$\begin{equation} R_{\rm NR}(E_{\rm iso})=ct_{\rm NR}(E_{\rm iso})=\left[{(3-k)E_{\rm iso} \over 4\pi Ac^2}\right]^{1/(3-k)}. \end{equation} \tag{ 4 }$$

For expansion in a constant density medium (k = 0), Equation (4) can be rewritten as

$$\begin{equation} R_{\rm NR}(E_{\rm iso})=0.3 \left({E_{51}\over \rho _0}\right)^{1/3} \left({f_b^{-1} \over 30}\right)^{1/3}\;{\rm pc}, \end{equation} \tag{ 5 }$$

while for a wind medium (k = 2), it becomes

$$\begin{equation} R_{\rm NR}(E_{\rm iso})= 2 \left({E_{51}\over A_\ast }\right) \left({f_b^{-1} \over 30}\right)\;{\rm pc}. \end{equation} \tag{ 6 }$$

From Equation (4) it follows that R_NR(E_iso) is a factor of ∼(E_iso/E)^1/(3−k) = f^−1/(3−k)_b ∼ θ^−2/(3−k)₀ larger than R_NR(E) = ct_NR(E) and a factor of ∼f^−1/(5−k)_b ∼ θ^−2/(5−k)₀ larger than R_ST[E, t_NR(E_iso)]. Thus, Equation (4) simply states that the jet keeps its original opening angle, θ_θ ≈ θ₀, until t_NR(E_iso), and hence at this time the jet is still far from being spherical. Thus, once the jet becomes sub-relativistic, we expect it to expand sideways significantly and become roughly spherical only when it has increased its radius by a factor of ψ. This should occur roughly at a time t_sph when R_ST(E, t_sph) = ψR_NR(E_iso):

$$\begin{equation} {t_{\rm sph} \over t_{\rm NR}}(E_{\rm iso})\approx f_b^{-1/2}\psi ^{(5-k)/2}\approx \sqrt{2}\,\theta _0^{-1}\psi ^{(5-k)/2}. \end{equation} \tag{ 7 }$$

This is a factor of ∼f^−1/2_b ≈ 14(θ₀/0.1)⁻¹ larger than the expected transition time for relativistic lateral expansion in the comoving frame. In the sections that follow we present quantitative estimates of t_sph for GRB jets expanding in variety of circumburst environments.

2.2. Initial GRB Model

Common to all calculations is the initiation of the GRB explosion as two identical blobs expanding in opposite directions into the circumburst medium. Calculations were done in two dimensions in cylindrical geometry using the piecewise parabolic method adaptive mesh refinement code FLASH (ver 2.4). Both blobs and the circumburst medium are modeled by a cold, γ = 5/3, ideal gas. The initial configuration is as follows. The computational domain, as illustrated in Figure 1, is a unprolonged cylinder in which the ejecta move along the symmetry axis. In the inner region of each of the pancakes, the ejecta mass, M_j, is distributed uniformly. In all runs v_j = 0.3c, ΔR_j/R_j ≈ 0.5, θ_j ≈ 0.5, M_j ≈ 2E_j/v²_j, and R_j ≈ R_NR(E_iso).

Zoom In Zoom Out Reset image size

Without a detailed understanding of the exact shape and energy distribution of the ejecta, we have only an approximate description of how to construct the initial conditions. However, as clearly illustrated by Ayal & Piran (2001), the late time evolution of the ejecta is rather insensitive to uncertainties in the initial conditions. This stems from the fact that at late times the mass of the remnant is dominated by the circumburst gas, which washes out any variations in the initial conditions of the ejecta. We have considered various initial densities, angular widths, and shapes of the collimated ejecta and found that these are indeed unimportant in determining the late morphology of the remnant provided that the total kinetic energy content per solid angle, $\mathcal {E}(\theta)$, is concentrated within some finite angular size, θ_c, and sharply drops outside of it. This is consistent with uniform jet structure models (Rhoads 1999; Sari et al. 1999), where $\mathcal {E}$ is assumed to be uniform only within some finite half-opening angle), or with Gaussian jets (Zhang & Mészáros 2002), where $\mathcal {E}(\theta) \propto \exp (-\theta ^2/2\theta _c^2)$.

An alternative jet model that has been proposed in the literature is one with a universal structure (Rossi et al. 2002; Zhang & Mészáros 2002), where $\mathcal {E}$ is assumed to vary smoothly with θ outside of some narrow core angle, typically with equal energy per decade in θ: $\mathcal {E} \propto \theta ^{-2}$. In this case, the wings contain more energy than the core (by about an order of magnitude) and, as a result, dominate the late time hydrodynamical evolution of the remnant. This class of models is discussed in Section 3.1.

In this section, we present a quantitative discussion of how the GRB remnant morphology is modified by expansion into a uniform medium. Expansion into a constant density medium is expected in a variety of progenitor models, in particular those related to short GRBs (Fryer et al. 1999; Bloom et al. 1999; Belczynski et al. 2006; Lee & Ramirez-Ruiz 2007). The discussion largely follows that of Ayal & Piran (2001), although analytical solutions are derived here to illustrate what may not be obvious from earlier derivations.

In the absence of characteristic scales in stellar ejecta and in the ambient medium, self-similar, spherically symmetric solutions exist, and they are widely used to interpret observational data of GRB remnants (e.g., Waxman et al. 1998). However, as argued in Section 2.1, near the non-relativistic transition the remnant is far from being spherical and thus the Sedov–Taylor solution fails to provide an accurate description of the evolution of the GRB ejecta. A simple estimate for t_sph can be, however, obtained by assuming that the expansion of each of the collimated blobs is accurately described by a self-similar, spherically symmetric solution. Under these conditions, the problem reduces to depositing a finite amount of energy E at two different locations separated by a distance R_j ∼ R_NR(E_iso). This is valid as long as the external medium is uniform around the explosion site. The evolution of the shock radius for each of the blobs will then follow:

$$\begin{equation} R_{\rm ST}(t)=\xi \left({E \over \rho _0}\right)^{1/5}\;t^{2/5} \end{equation} \tag{ 8 }$$

with ξ = 1.17 for γ = 5/3. The ratio between the remnant width and height,

$$\begin{equation} \varsigma ={R_{\rm ST}(t)\over R_{\rm ST}(t) + R_{\rm j}}, \end{equation} \tag{ 9 }$$

will approach unity as the two blobs expand and merge. The GRB remnant will then become nearly spherical in shape after a time

$$\begin{equation} t_{\rm sph} \approx 243 \left({E_{51} \over \rho _0}\right)^{1/3}\left({f_b^{-1} \over 30}\right)^{5/6}\left({\varsigma _{0.9} \over 1- \varsigma _{0.9}}\right)^{5/2}\;{\rm yr}, \end{equation} \tag{ 10 }$$

when ς = 0.9ς_0.9. Obviously, the above calculation is only sketchy and should be taken as an order of magnitude figure, as it not only assumes that the mass of the collimated ejecta is negligible with respect to the swept-up gas but also neglects the presence of shocked material throughout the interaction region.

Detailed hydrodynamic simulations of the evolution of a GRB remnant in a uniform medium are presented in Figure 2, where the density contours of the expanding collimated ejecta at various times in its hydrodynamical evolution are plotted. As the collimated ejecta collides with the ambient medium, a bow shock forms. The shock propagates in the direction of motion but also perpendicular to it and, over time, wraps around the expanding ejecta. Although initially the remnant may be highly nonspherical, the ratio between its width and height will approach unity as the two blobs expand, merge into a single structure, and then finally become spherical in shape. The resultant structure will not be perfectly spherical as some density inhomogeneities around the equator resulting from the encounter remain visible. It will be, however, difficult to distinguish it (based on morphology alone) from an SN remnant after about

$$\begin{equation} t_{\rm sph} \approx 3 \times 10^3\left({E_{51}\over \rho _0}\right)^{1/3}\;{\rm yr}, \end{equation} \tag{ 11 }$$

when ς ∼ 0.9. The governing parameters of the late evolution of the remnant are the initial energy of the jet, E, and the density of the ambient medium, ρ (Ayal & Piran 2001). As illustrated in Figure 3, these two initial parameters can also determine the early evolution of the remnant for a fixed f_b.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Structured Jets

The morphology and dynamics of young GRB remnants, as we are introducing in this paper, depend on the distribution of the ambient medium and on the structure of the GRB ejecta, whose angular energy distribution we assume to be mainly concentrated along the jet's polar angle. In this section, we relax this assumption and consider the hydrodynamical evolution of jets with additional ejected material moving at large polar angles. Such complex jet structures are predicted in the context of neutrino (Levinson & Eichler 2000; Aloy et al. 2005; Rosswog & Ramirez-Ruiz 2003; Rosswog et al. 2003; Dessart et al. 2009) and hydromagnetically (Vlahakis et al. 2003; Komissarov et al. 2009; Lyubarsky 2009) driven jets and in the context of the cocoon in the collapsar model (Ramirez-Ruiz et al. 2002; Zhang et al. 2003; Morsony et al. 2007; Mizuta & Aloy 2009).

For this discussion, we shall assume that the energy content of the GRB ejecta does not sharply drop outside θ_j but instead smoothly varies with θ as a power law with equal energy per decade in θ: $\mathcal {E} \propto \theta ^{-2}$ · $\mathcal {E} \propto \theta ^{-2}$. The evolution of a structured GRB remnant since the onset of the non-relativistic phase is shown in Figure 4 for three different angle-integrated energies, E_w, of the $\mathcal {E}\propto \theta ^{-2}$ wings. The initial parameters of the GRB ejecta for θ ⩽ θ_j are chosen so that they are equal to those displayed in Figure 3. Compared to Figure 3, clear structural differences appear when the GRB outflow has grown significantly in size and the wings start to overtake the laterally expanding jet component. Three illustrative cases are depicted with E_w/E_j: 0.25 (left panel), 1.0 (middle panel), and 4.0 (right panel). In this case, the outflow ejected at large angles drives a blast wave that eventually sweeps up all the laterally expanding, jet-shocked medium. Not surprisingly, the ejection of material at large angles carrying a non-negligible fraction of the total GRB energy significantly modifies the simple estimate given by Equation (11) and limits our ability to decipher the presence of a beamed component. In such cases, the GRB dynamics will be accurately described by the classical Sedov–Taylor SN remnant evolution.

Zoom In Zoom Out Reset image size

If the progenitors of GRBs are massive stars then there is an analogy to the explosions of core collapse SNe, for which there is abundant evidence that they interact with the winds from the progenitor stars. In most SN cases, the radial range that is observed is only out to a few parsecs, such that the mass loss characteristics have not changed significantly during the time that mass is supplied to the wind (Chevalier & Li 2000). The density in the wind depends on the type of progenitor. Red supergiant (RSG) stars, which are thought to be the progenitors of Type II SNe, have slow dense winds. Wolf-Rayet stars, which are believed to be the progenitors of Type Ib/c SNe and possibly of long GRBs (Woosley 1993; MacFadyen & Woosley 1999), have faster, lower-density winds. The winds from Wolf-Rayet stars are characterized by mass-loss rates $\dot{M}\approx 10^{-5}\,M_\odot \;{\rm yr}^{-1}$ and velocities v_w ≈ 10³ km s^-1 (Chiosi & Maeder 1986). In a steady, spherically symmetric wind, the stellar density is ρ_ext = Ar⁻², where $A=5\times 10^{11}(\dot{M}/10^{-5}\,M_\odot \;{\rm yr}^{-1})(v_w/10^3\;{\rm km \; s^{-1}})^{-1}\;{\rm g\;cm^{-3}}$. Note that for this choice of stellar wind parameters A_* = 1.

For this discussion, we shall first assume that the stellar wind is effectively spherical. The evolution of a GRB remnant in a stellar wind since the onset of the non-relativistic phase is summarized in Figure 5. Similar resulting structures to those described in Section 3 are clearly seen. A bow shock forms as each blob collides with the stellar wind, which eventually wraps around the ejecta before the two expanding shells collide to form a single structure. However, because in a wind medium the swept-up mass increases only linearly with radius, the GRB remnant decelerates much more slowly than in a uniform medium. Moreover, in a wind medium, resistance to sideways expansion is increased. This is because the bow shock, as it wraps around the ejecta, encounters a steadily increasing ambient pressure. As a result, the remnant will become roughly spherical only after a time,

$$\begin{equation} t_{\rm sph} \approx 5\times 10^4\;\left({E_{51}\over A_\ast }\right)\;{\rm yr}, \end{equation} \tag{ 12 }$$

when ς ∼ 0.9. Beyond this point, the evolution will evolve into a classical Sedov–Taylor SN remnant evolution.

Zoom In Zoom Out Reset image size

The estimate given by Equation (12) could be inaccurate for a number of reasons. Depending upon the wind history of the progenitor star and the properties of the surrounding interstellar medium (ISM), the density structure around R_sph ≈ 30 pc could be quite complicated. The non-steady nature of the winds in massive stars together with the relatively large ISM pressure expected in star-forming regions leaves open the possibility of interaction with denser material at much earlier times. In this case, the GRB remnant will start being decelerated by the external medium at a smaller radius than it would expanding into a free 1/r² wind. Much of our effort in Section 6 will therefore be dedicated to determining the contribution of the presupernova ejecta of Wolf-Rayet stars to the circumburst environment, and describing how this external matter can affect the observable characteristics of GRB remnants.

Large-scale density gradients in the ambient medium could result in asymmetric, nonradial distortions. In this case, t_sph could be larger than that given by Equation (12). For example, Figure 6 shows the evolution of a GRB remnant in an asymmetric wind, under the assumption that the progenitor star experiences nonspherical mass loss close to critical rotation: in other words, a scenario in which a slower and denser wind is confined to the equatorial plane. To compute the latitudinal dependence of the wind properties of a star close to critical rotation ideally requires multi-dimensional models of the star and its outflowing atmosphere, which are not available. Langer (1998), however, argued that the stellar flux and the radius might still vary only weakly from pole to the equator in very luminous stars. We therefore applied equations similar to those found by Bjorkman & Cassinelli (1993) for winds of rotating stars in the limit of large distance from the star:

$$\begin{equation} v_{\infty }(\theta) = \zeta v_{\rm esc} \left(1 - \Omega \, \sin \theta \right)^{\varphi }\;, \end{equation} \tag{ 13 }$$

where we set the parameters defined in Bjorkman & Cassinelli (1993) to ζ = 1, φ = 0.5, Ω = v_rot/v_crit = 0.8, and $v_{\rm crit}=v_{\rm esc}/\sqrt{2}=[GM_*(1-\kappa)R_*)]^{1/2}$, with M_* and R_* being the mass and radius of the star, and κ standing for the ratio L/L_Edd of stellar to Eddington luminosity. Under the above conditions, the GRB remnant expands more quickly and easily into the lower density wind at the poles, producing an increasingly asymmetric double-lobed structure.

Zoom In Zoom Out Reset image size

Finally, the estimate given by Equation (12) would be modified if the beamed GRB is accompanied by an underlying SN, as expected in the collapsar model (MacFadyen & Woosley 1999; MacFadyen et al. 2001; Woosley & Bloom 2006). The large-angle SN outflow, responsible for exploding the star and producing the ⁵⁶Ni, would generally carry more energy and inertia than the relativistic jet itself (e.g., Soberberg et al. 2004; Mazzali et al. 2006; Kaneko et al. 2007), so that the latter always overtakes it, and sweeps up the GRB ejecta. It is to this problem that we now turn our attention.

It seems likely that GRBs originate in a very small fraction of massive stars that undergo a catastrophic energy release event toward the end of their evolution. Expressly, the association of some GRBs with type Ic SNe (e.g., Hjorth et al. 2003; Stanek et al. 2003; Pian et al. 2006) has pointed a finger at deaths of massive Wolf-Rayet stars as the cause of GRBs, or at least a subset thereof. The central engine is believed to give rise to a polar outflow with two components (MacFadyen & Woosley 1999; Ramirez-Ruiz et al. 2002; Woosley & Bloom 2006). One large-angle outflow (the SN), containing most of the energy and mass, is responsible for exploding the star and producing the ⁵⁶Ni to make the SN bright. A second outflow component (the GRB jet) occupies a narrower solid angle and probably contains less energy (which can range from comparable to much less), and most of its energy is in material with relativistic velocities (where the typical Lorentz factor of the material that carries most of the energy in this component can vary significantly among SN–GRBs).

The large-angle SN outflow, carrying more energy and inertia than the relativistic jet itself, will generally sweep up the GRB ejecta before it has been much decelerated. An order of magnitude estimate for t_sph can be obtain by assuming that the dynamics of the laterally expanding, GRB ejecta is accurately described by a self-similar, spherical solution in a 1/r² medium and that the SN outflow does not appreciably slow down. The evolution of the shock radius of the beamed remnant is then given by

$$\begin{equation} R_{\rm ST}(t)=\xi _w \left({E \over A}\right)^{1/3}\;t^{1/3} \end{equation} \tag{ 14 }$$

with ξ_w = 0.73 for γ = 5/3. The GRB ejecta, as it clears the surrounding stellar matter, will be overtaken by the large-scale outflow SN at

$$\begin{equation} t_{\rm sph} \approx 910 \left({E_{51} \over A_\ast }\right) \left({\beta _{\rm SN} \over 0.1}\right)^{-3}\;{\rm yr}, \end{equation} \tag{ 15 }$$

where β = v_SN/c ≪ 1. After this time, the merged system will quickly become spherical. Yet it is clear that this simple estimate is inadequate as a model for the real complex dynamics, which necessitates the use of hydrodynamical calculations.

The evolution of a GRB remnant accompanied by an underlying SN is shown in Figure 7 for two different SN explosion energies. The stellar wind parameters are chosen so that they are equal to those displayed in Figure 5. Compared to Figure 5, significant structural differences appear when the SN outflow has grown significantly in size and it starts to overtake the laterally expanding GRB ejecta. Two illustrative cases are depicted: E_SN: 5 × 10⁵¹ erg (left panels) and 5 × 10⁵⁰ erg (right panels). In both cases, the SN outflow carries more momentum than the beamed ejecta and drives a blast wave that eventually sweeps up all the GRB-shocked medium. This happens before the two beamed blobs collide on the equator. The merged blobs will become roughly spherical very soon after. For E_SN = 5 × 10⁵¹ erg (5 × 10⁵⁰ erg), we obtain t_sph ≈ 250 yr (3 × 10³ yr). Not surprisingly, the presence of an underlying spherical SN seriously modifies the simple estimate given by Equation (12) and limits our ability to decipher the presence of a beamed component in a GRB explosion.

Zoom In Zoom Out Reset image size

Prevailing in all these calculations is the initiation of the explosion as a pressure-driven blast wave by deposition of the explosion energy, E_SN, entirely as thermal energy. In the inner region, an ejecta mass, M_SN ≈ 5 M_☉, is distributed uniformly. This may accurately model the Sedov–Taylor stage of SN remnant evolution after the ratio of swept-up mass to the mass of the original stellar ejecta exceeds roughly 19 (Fabian et al. 1983). In most cases, deceleration of the SN outflow will begin only sometime after the GRB ejecta has been swept up. As clearly seen in Figure 8, the SN outflow at this early stage is not accurately described by the Sedov–Taylor solution. However, numerous tests show that our results are not strongly dependent upon the assumed mass of the SN ejecta or on whether kinetic energy rather than thermal energy is distributed (such that the velocity profile is linear; similar to the Sedov solution) in the inner region. None of these complications are likely to seriously modified our estimate for t_sph.

Zoom In Zoom Out Reset image size

So far we have considered either the uniform ambient medium case or the 1/r² wind case on its own. However, since the winds in massive stars are non-steady, the density structure is more complex (Wijers 2001; Chevalier et al. 2004; Ramirez-Ruiz et al. 2005; van Marle et al. 2006). The preburst stellar wind depends on the evolutionary stages prior to (and during) the Wolf-Rayet stage (Ramirez-Ruiz et al. 2001). For Galactic stars, a standard evolutionary track is to start as an O star, evolve through an RSG phase or luminous blue variable phase with considerable mass loss, and end as a Wolf-Rayet star (García-Segura et al. 1996a; García-Segura et al. 1996b). At low metallicity, the RSG phase may be absent; this may also be the case for some binary stars.

As an example, we follow the dynamics of a GRB remnant around a 35 M_☉ star (as calculated by Chevalier et al. 2004), which evolves (at solar metallicity) through a long-lived RSG stage with prominent consequences for the evolution of the circumstellar matter. The wind velocity in the Wolf-Rayet phase is 10³ km s⁻¹, and the mass-loss rate 10⁵ M_☉ yr⁻¹. The ISM pressure and density are assumed to be typical of the hot, low-density phase of a starburst galaxy, with P_ism ∼ 10⁷ K cm⁻³ and a density of 4 × 10⁻²⁵ g cm⁻³. When the fast Wolf-Rayet wind v_w ∼ 10³ km s⁻¹ starts blowing, it sweeps up the RSG wind material into a shell. The termination shock of the Wolf-Rayet wind is located at R_t ≈ 0.4 pc and the RSG shell at R_rsg ≈ 1.7 pc. Because the pressure in the shocked wind is nearly in equilibrium with the ISM, and the temperature ∼10⁷(v_w/10³ km s^-1)² K, the density in the bubble is ∼8 × 10²⁵(P_ism/10⁷)(v_w/10³ km s^-1)⁻² g cm^-3, independent of the mass-loss rate and the ambient density. The extent of the constant density region is ∼4R_t out to the dense RSG shell.

The resulting evolution—a beamed GRB remnant in the circumstellar medium expected around a 35  M_☉ massive stellar progenitor—is summarized in Figure 9. The presence of the sharp density gradient will only affect the dynamics of the GRB remnant when it size is comparable to, or exceeds, the scale length of the gradient. Before this time, the evolution of the remnant is similar to that depicted in Figure 5. The bow shock propagates in the direction of motion but also perpendicular to it and, over time, wraps around the expanding ejecta. In a wind bubble, the shock front will expand within the stellar wind until it reaches the sharp density discontinuity at about 1.7 pc. The encounter with the RSG shell happens before the GRB remnant has time to expand laterally, which allows for an elongation of the RSG shell in the z-axis. A less pronounced elongation is also seen perpendicular to the z-axis, which results from the collision on the equator of the two beamed blobs. Superimposed on this large-scale deformation one can also notice the familiar effects caused by the development of Rayleigh–Taylor (R-T) instability. The instability grows rather quickly and within the following 10³ yr several elongated spikes extend from the shell. At ≈2 × 10³ yr, the remnant has a mean radius of 3 pc (Figure 10) and the spikes are even more pronounced than before. Both shell and R-T spikes advance with a velocity of 300–500 km s⁻¹, which is a very small fraction of the random velocities observed in the hot cavity of the remnant. These fast motions are induced by inflection of consecutive shock waves ramming the irregular shell. The pressure in the cavity is almost uniform, whereas the density varies chaotically. The knots of ejecta formed by the R-T instability are believed to be responsible for several forms of observable radiation in some young supernova remnants, including radio synchrotron and optical emission lines. The sheared motions resulting from the instability could lead to an amplification of the magnetic field strength and enhance the bright radio remnant (Jun & Norman 1996).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The growth time of the instability may be roughly estimated in the following way. The interaction of the GRB remnant with the wind cavity leads to an increase of the shell-driving pressure by

$$\begin{equation} \Delta P={E_{\rm j} \over 2\pi R_{\rm rsg}^3}. \end{equation} \tag{ 16 }$$

The mass of the shell is equal to

$$\begin{equation} M_{\rm rsg}={4 \over 3}\pi R_{\rm rsg}^3\rho _{\rm rsg}. \end{equation} \tag{ 17 }$$

ΔP causes an acceleration g_rsg of the shell, satisfying the equation

$$\begin{equation} 4\pi R_{\rm rsg}^2 \Delta P=M_{\rm rsg} g_{\rm rsg}. \end{equation} \tag{ 18 }$$

With no other forces involved, the contact surface separating both fluids is unstable to perturbations of all wavelengths and the fluids interpenetrate. The instability grows exponentially on a characteristic time scale

$$\begin{equation} t_{\rm R-T}= \left({4\pi ^2\lambda R_{\rm rsg}^4\rho _{\rm rsg} \over 3 E_{\rm j}}\right)^{1/2}, \end{equation} \tag{ 19 }$$

where λ is the perturbation wavelength, corresponding to a growth time of

$$\begin{equation} t_{\rm R-T} \approx 10^3\left({R_{\rm rsg} \over 2\;{\rm pc}}\right)^2\left({\lambda \over 0.2\;{\rm pc}}\right)^{1/2}\rho _{{\rm rsg},-21}^{1/2}E_{51}^{1/2}\;{\rm yr}, \end{equation} \tag{ 20 }$$

where E₅₁ is the total energy of the GRB remnant and ρ_rsg=10⁻²¹ρ_rsg,−21 g cm^-3. Indeed, the observed spikes have dimensions of tenths of parsecs and grow on a time scale comparable to the estimated one. For this discussion, we have assumed that the remnant evolution is effectively adiabatic and should be modified to include the effects of radiative cooling,⁴ which are expected to become important after about few thousand years.

The disruption of the RSG shell is, as expected, very sensitive to the total amount of energy and momentum released by the GRB explosion. To illustrate this, we study the evolution of a GRB remnant accompanied by an underlying 10⁵² erg SN in a wind bubble structure such as that illustrated in Figure 9. The large-angle SN outflow, carrying significantly more energy and inertia than the relativistic jet itself, sweeps up the GRB ejecta before it is decelerated by the RSG shell. The resultant shell, clearly apparent in the first frame of Figure 11 taken at t = 380 yr, will then be pushed outward at faster velocity than it would be in the absence of the SN ejecta. In subsequent evolutionary phases one can observe a broadening of the merged shell, accompanied by a decreasing of its density, and the development of small-amplitude density perturbations. After t = 3180 yr, the remnant has become highly irregular already when it has grown to ≈13 pc and is also rather weakly elongated in the direction of motion of the beamed ejecta (Figure 12). There are also multiple kinematic components within the remnant. Fast-moving knots and fast-moving flocculi, dense fragments of SN ejecta, expand from the explosion center with velocities of several thousand km s⁻¹. Much slower (several hundred km s⁻¹) flocculi are clearly shocked and accelerated remnants of stellar material ejected by the SN progenitor prior to the explosion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The calculations above demonstrate how the measurable properties of GRB remnants depend on the nature of the progenitor star and the medium around it. Moreover, counts of GRB remnants as a function of age may have huge selection effects, as the actual age of the GRB may be considerably less than the kinematical age estimated from the radius of the filaments divided by the expansion velocity. Although nurture makes a huge difference, hydrodynamical models of young GRB remnants are also sensitive to the structure of the ejecta, which might be rather complicated in detail as demonstrated by inspection of numerical models of SN explosions (e.g., Chevalier & Blondin 1995).

7.1. Asymmetric GRB Remnants in the Local Universe

It is obvious from the discussions in this paper that the dynamics of GRB remnants are complex, especially because of rich interactions between the ejecta and the circumburst medium. The structure of ejecta is also important, in particular for young GRB remnants. Abundant confirmation was provided of the important notion that the morphology and visibility of GRB remnants are determined largely by their circumstellar environment, from the initial density gradient created by the progenitor to the effects of large- and small-scale circumburst structures in later evolution. This distribution is expected to be nonuniform around massive GRB progenitors because of the significant mass loss and the dynamical effects of the stellar winds.

The presence of a density gradient will only affect the dynamics of a GRB when the remnant size is comparable to, or exceeds, the scale length of the gradient. Before this time, the density can be treated as approximately uniform. In the absence of characteristic scales in stellar ejecta and in the ambient medium, self-similar, spherically symmetric solutions exist, and they are widely used to interpret observational data on young GRB remnants. However, even for the simplest density distributions, we found that the resulting structure and dynamics of their remnants are very different from the standard self-similar solutions. This is mainly because at early stages the morphology of a beamed GRB remnant would be very different from that of a spherical explosion. In principle, this can be used to identify those remnants although the dynamical complexity of their surrounding circumburst medium seriously limits our ability to decipher their presence, in particular around massive star progenitors.

The values of the mean space density, ϕ_sph, of asymmetric GRB remnants expanding into a variety of circumburst environments are given in Table 1. Taking the long GRB rate expected in the local universe to be ℜ_grb = 0.3 Gpc⁻³ yr⁻¹ (Guetta & Piran 2007), an upper limit to the number density of asymmetric remnants arising from massive stellar progenitors is

$$\begin{equation} \phi _{\rm sph} \approx 5 \times 10^{-3}E_{51}A_\ast ^{-1}{f_{b}^{-1} \over 100}\;{\rm Mpc^{-3}}. \end{equation} \tag{ 21 }$$

This is a generous upper limit, since it assumes that all long GRBs occur in a free 1/r² stellar wind and that the relativistic component is energetically dominant. For comparison, an ultraconservative lower limit to the space density of long GRBs gives

$$\begin{equation} \phi _{\rm sph} \approx 10^{-5}{f_b^{-1} \over 100}\;{\rm Mpc^{-3}}. \end{equation} \tag{ 22 }$$

With the rates given by Guetta & Piran (2005), the frequency of short GRBs is about 20 times higher than that of long GRBs. The much lower external density (likely uniform) medium expected around short GRBs (e.g., Lee et al. 2005) suggests

$$\begin{equation} \phi _{\rm sph} \approx 6 \times 10^{-3}E_{51}^{1/3} \rho _{-3}^{-1/3} {f_b^{-1} \over 10}\;{\rm Mpc^{-3}}, \end{equation} \tag{ 23 }$$

where ρ₋₃ = 10⁻³ρ₀. Thus, the rate of asymmetric short GRB remnants could be much higher than Equation (22). However, the remnant's visibility would be highly biased in favor of those with massive progenitors.

Table 1. Mean Space Density of Asymmetric GRB Remnants

Remnant Type	tsph	ϕsph(ℜgrb = 1 Gpc-3 yr−1)
k = 0	3 × 103E1/351ρ−1/30 yr	3 × 103E1/351ρ−1/30f−1bℜgrb Gpc-3
k = 2	5 × 104E51A−1* yr	5 × 104E51A−1*f−1bℜgrb Gpc-3
Rrsg∼ 1 pc	102E51A−1* yr	102E51A−1*f−1bℜgrb Gpc-3
Esn ∼ Ej,51	3 × 103 yr	3 × 103f−1bℜgrb Gpc-3
Esn ⩾ Ej,51	⩽103 yr	⩽103f−1bℜgrb Gpc-3

Download table as:  ASCIITypeset image

7.2. Observational Prospects

The most difficult task at present is to relate hydrodynamical modeling to observations. A few of the observables, such as expansion rates and thicknesses of the flow structures, can be relatively easily determined from the models. However, modeling radio and X-ray emission is in general difficult, as we are still lacking an understanding of how electrons are accelerated in shocks. Very similar difficulties are encountered in modeling nonthermal X-rays. Thermal X-ray spectra are in principle easier to model, but in practice the difficulties are formidable. The reason for these difficulties is our poor understanding of a number of topics, such as the amount of electron heating in collisionless shocks, the detailed structure and composition of ejecta, their clumping, the presence of the inhomogeneous circumstellar medium, and the presence of dust.

The difficult task of interpreting observations with the help of hydrodynamical models is perhaps best illustrated by Cas A. This is a remnant of a massive star explosion and a classic prototype shell SN remnant. It has been detected throughout the whole electromagnetic spectrum (Ryle & Smith 1948; Ashworth 1980; Aharonian et al. 2001; Fesen 2001; Hwang et al. 2004). Observations in various wavelength bands probe very different components of the remnant: synchrotron radio emission gives us information about relativistic electrons, thermal X-ray emission is produced by the bulk of the shocked hot gas, much cooler gas in radiative shocks emits at optical wavelengths, and observations in infrared reveal still cooler gas and dust. However, we still do not understand what is the relationships between all these features and the remnant's hydrodynamics (e.g., Hwang et al. 2004).

To estimate the emissivity of GRB remnant without undertaking the complicated effort of calculating X-ray spectra, we have summed up the internal energy density of the gas, _int, for each zone in the simulations and plotted L_X ∼ f_Xint as a function of time in Figure 13, where f_X ∼ 0.01 is the fraction of internal energy that is radiated away in the 0.3–10 keV range. This is only a very approximate procedure and should be taken as an order of magnitude estimate at present. Superimposed on these plots, on a common scale, are all GRB afterglows with X-ray luminosity measurements covering several tens of days. There are, unfortunately, only a few curves, because such measurements can only be made on GRBs that are relatively nearby, but their light curves should be illustrative. We then compared these X-ray light curves with those of SNe and historical Galactic SN remnants. The resulting plot is striking in several ways. Despite the huge disparity in initial appearance, there are indications of a common convergence of all classes of phenomena to a common resting place: L_X ∼ 10³⁹–10⁴⁰ erg s⁻¹ about a few years after the explosive event. Moreover, it clearly illustrates that the transition from a GRB to an SN remnant appears to be rather smooth. Clearly, detailed studies relating hydrodynamical modeling to observations are needed to study the transition from a GRB into a stellar remnant.

Zoom In Zoom Out Reset image size

Figure 13 shows the schematic light curves of two GRB remnants expanding into different external density environments. For a 1/r² density profile, as in the black curve of Figure 13, the luminosity of the remnant is initially dominated by the emission of the individual beamed components as they interact with the stellar wind. The luminosity continues with a quasi-steady decay rate until the individual beamed components collide to form a single structure. The resultant light curve will then be characterized by a modest increase in luminosity. For a GRB expanding inside a pre-existing wind bubble, as in the gray curve of Figure 13, the resultant light curves can evolve much faster into luminous remnants such as Cas A (Hwang et al. 2004) or W49B (Miceli et al. 2006) due to their strong interaction with the dense circumburst medium. Interestingly, in both of these young remnants a jet feature is clearly observed (Lopez et al. 2009a, 2009b; Rest et al. 2010)

This plot summarized many of the issues outlined in this paper in which we argued that the morphology and visibility of GRB remnants are governed by the pre-existing structure of their birthplace environments. Since GRB remnants result from the impact of their ejecta with circumstellar gas, their visibility is highly biased in favor of those with massive progenitors. Many young GRBs from massive progenitors would be bright because their ejected mass is interacting with nearby gas expelled by the progenitor itself. This circumstellar gas is likely to have mass comparable to that of the accompanied SN debris and will not extend much further than a few parsecs. After less than a century, the blast wave from the GRB will pass through this relatively dense circumstellar gas. Inferring the presence of a beamed component in these hypernova scenarios would be challenging. The number density of asymmetric GRB remnants in the local universe could be far larger if they expand in a tenuous ISM as expected, for example, in the merger of two neutron stars (although there are reasons to suspect that the ejecta may not be too narrowly beamed) and may be easier to constrain directly (acknowledging the obvious trade-offs in sensitivity and angular resolution, particularly for radio and X-ray observations).

We have benefited from many discussions with C. Fryer, J. Granot, L. Lopez, T. Piran, and S. Woosley. We are especially grateful to W. Zhang and R. Parsons for countless insightful conversations. We also thank the referee for useful suggestions. The software used in this work was in part developed by the DOE-supported ASCI/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. Computations were performed on the IAS Scheide computer cluster. This work is supported by NSF: PHY-0503584 (ER-R) and DOE SciDAC: DE-FC02-01ER41176 (ER-R and AM).