A BEAMING-INDEPENDENT ESTIMATE OF THE ENERGY DISTRIBUTION OF LONG GAMMA-RAY BURSTS: INITIAL RESULTS AND FUTURE PROSPECTS

The energy budget of gamma-ray bursts (GRBs) provides fundamental insight into the nature of the explosions, the resulting ejecta properties, and the identity of the central compact remnant ("engine"). While the isotropic-equivalent γ-ray energy (E_{γ, iso}) can be easily determined from a measurement of the burst fluence and redshift, a complete accounting of the energy budget requires detailed observations of the afterglow emission. The afterglow observations provide a measure of the isotropic-equivalent blast-wave kinetic energy (E_{K, iso}), as well as the explosion geometry (quantified by a jet opening angle, θ_j). The resulting beaming corrections, f⁻¹_b ≡ 1 − cos(θ_j), can be substantial, approaching three orders of magnitude in some cases (Frail et al. 2001; Panaitescu & Kumar 2002; Berger et al. 2003a; Bloom et al. 2003). To properly determine E_{K, iso} and f_b it is essential to observe the afterglows from radio to X-rays over timescales of hours to weeks, clearly a challenging task. This is particularly a problem for the subset of "dark" GRBs for which the lack of detected optical emission, or large extinction, prevent a determination of E_{K, iso} and likely f_b (e.g., Berger et al. 2002; Piro et al. 2002).

Over the past decade, detailed afterglow observations have been obtained at a great cost of telescope time for about 20 long-duration GRBs, with the basic result that the beaming corrections are large and diverse, leading to typical true energies of E_γ ∼ E_K ∼ 10⁵¹ erg (Frail et al. 2001; Panaitescu & Kumar 2001, 2002; Berger et al. 2003a, 2003b; Bloom et al. 2003). More recently, it has been recognized that some nearby long GRBs have much lower energies, E_iso  ∼  10⁴⁹–10⁵⁰ erg, and appear to be quasi-isotropic (Kulkarni et al. 1998; Soderberg et al. 2004b; Soderberg et al. 2006). Similarly, some bursts appear to have large beaming-corrected energies of ∼10⁵² erg (Cenko et al. 2010a, 2010b). The existence of these highly energetic bursts depends at least in part on the ability to correctly infer their large beaming corrections. Indeed, the inference of jet opening angles from breaks in the afterglow light curves has become controversial in recent years due to conflicting trends in optical and X-ray light curves (Liang et al. 2008; Racusin et al. 2009). Similarly, in some cases a two-component jet has been inferred, with a narrow core dominating the γ-ray emission and a wider component dominating the afterglow emission (Berger et al. 2003b; Racusin et al. 2008). Numerical simulations suggest that off-axis viewing angles can also lead to shallow breaks that may be missed or misinterpreted (van Eerten et al. 2010).

In addition to potential difficulties with the inference of f_b, the γ-ray and kinetic energies measured from the early afterglow emission only pertain to the relativistic ejecta. The existence of a substantial component of mildly relativistic ejecta can only be determined from observations at late times when such putative material can refresh the forward shock. Clearly, the existence of substantial energy in a slow ejecta component will place crucial constraints on the activity lifetime of the central engine.

Such late-time observations also have the added advantage that they probe the blast wave when it has decelerated to non-relativistic velocities and hence roughly approaches isotropy (Frail et al. 2000; Livio & Waxman 2000). This allows us to use the well-established Sedov–Taylor self-similar solution, with negligible beaming corrections, to estimate the total kinetic energy of both the decelerated ejecta and any additional initially non-relativistic material. Since the peak of the afterglow spectrum on these timescales is located in the radio band, the lack of optical afterglow emission (e.g., due to extinction) does not have an effect on the ability to determine E_K.

This approach was first exploited by Frail et al. (2000) to model the late-time radio afterglow emission of GRB 970508 (at δt ≳ 100 days) from which the kinetic energy was inferred to be E_K ∼ 5 × 10⁵⁰ erg. Berger et al. (2004) used the same approach to model the radio afterglow emission of GRB 980703 on timescales of ≳ 40 days, and to re-model GRB 970508. They found kinetic energies of E_K ∼ 3 × 10⁵¹ erg for both bursts. Finally, Frail et al. (2005) modeled the radio emission from GRB 030329 at δt ≳ 50 days and found E_K ∼ 10⁵¹ erg. Only three bursts have been studied in this fashion so far because only those events have well-sampled radio light curves on the relevant timescales of δt ≳ 100 days.

However, the kinetic energy can still be estimated using the same methodology even from fragmentary late-time radio observations. Such an approach will naturally result in larger uncertainties for each burst, but it can be applied to a much larger sample of events. Here we present such an analysis for 24 long-duration GRBs with radio observations at ≳ 100 days, but with only 1–3 data points (at 1.4 to 8.5 GHz) per burst. Using these observations, we infer robust ranges for the kinetic energy of each burst and for the population as a whole. The plan of the paper is as follows. The radio observations are summarized in Section 2. The model for synchrotron emission from a Sedov–Taylor blast wave, and the various assumptions we employ are presented in Section 3. In Section 4 we detail the resulting kinetic energies and the range for the overall sample, and we compare these results to multi-wavelength analyses of early afterglows in Section 5. We conclude with a discussion of future prospects.

We use radio observations of 24 long GRBs at δt ≳ 100 days, since on those timescales the blast wave is expected to become non-relativistic and roughly isotropic (Livio & Waxman 2000) and the peak of the afterglow emission is at or below the centimeter band. This has been confirmed with detailed data in the case of GRBs 970508, 980703, and 030329 (Frail et al. 2000, 2005; Berger et al. 2004). We restrict the analysis to GRBs with a known redshift and with early-time detections, which, for the case of a single detection or upper limit, allow us to infer that the peak of the spectrum has transitioned below our observing frequency.

The observations are primarily from the Very Large Array (VLA¹), with the exception of GRBs 980425 and 011121 which were observed with the Australia Telescope Compact Array (ATCA). The data were obtained between 1997 and 2009 as part of a long-term GRB radio program (e.g., Frail et al. 2003a).

For the purpose of our analysis, we separate the bursts into three categories based on the quality of the data. In Group A are three bursts with late-time detections at multiple frequencies that constrain the peak of the synchrotron spectrum (the same three bursts that have been studied in detail by Frail et al. 2000; 2005; Berger et al. 2004). In Group B are 11 bursts with single-frequency detections, while Group C consists of 10 GRBs with late-time non-detections. The VLA measurements and relevant burst properties are listed in Table 1.

Our modeling of the radio data follows the methodology of Frail et al. (2000) and Berger et al. (2004) for the case of a uniform density medium.² For the typical expected parameters of long GRBs, the initially collimated blast wave approaches spherical symmetry and decelerates to non-relativistic velocity on similar timescales, t_s ≈ 150(E_{K, iso, 52}/n_e)^1/4t^1/4_{j, d} days and t_NR ≈ 40(E_{K, iso, 52}/n_e)^1/4t^1/4_{j, d} days, respectively (Livio & Waxman 2000); here, n_e is the circumburst density in units of cm⁻³ and t_j is the "jet break" time at which the jet begins to expand sideways (i.e., Γ(t_j) ∼ θ⁻¹_j, where Γ is the bulk Lorentz factor). In this paper, we assume that the blast wave has transitioned to the non-relativistic isotropic phase by the time of our observations and subsequently check for self-consistency.

The blast wave dynamics in the non-relativistic phase are described by the Sedov–Taylor self-similar solution with r(t)∝(E_STt²/n)^1/5. To calculate the synchrotron emission emerging from the shock-heated material, we make the usual assumptions: (1) the electrons are accelerated to a power-law energy distribution, N(γ)∝γ^−p for γ > γ_m, where γ_m is the minimum Lorentz factor; (2) the value of p is 2.2 as inferred from several bursts (e.g., Panaitescu & Kumar 2001, 2002; Yost et al. 2003); and (3) the energy densities in the magnetic field and electrons are constant fractions (_B and _e, respectively) of the shock energy density. Accounting for synchrotron emissivity and self-absorption, and including the appropriate redshift transformations, the flux observed at frequency ν and time t is given by (Frail et al. 2000; Berger et al. 2004)

$$\begin{equation} F_\nu = F_0(t/t_0)^{\alpha _F}[(1+z)\nu ]^{5/2}(1-e^{-\tau _\nu }) f_3(\nu /\nu _m)f^{-1}_2(\nu /\nu _m), \end{equation} \tag{ 1 }$$

where the optical depth is given by

$$\begin{equation} \tau _\nu = \tau _0(t/t_0)^{\alpha _{\tau }}[(1+z)\nu ]^{-(p+4)/2} f_2(\nu /\nu _m), \end{equation} \tag{ 2 }$$

the synchrotron peak frequency, corresponding to electrons with γ = γ_m, is given by

$$\begin{equation} \nu _m = \nu _0(t/t_0)^{\alpha _m}(1+z)^{-1}, \end{equation} \tag{ 3 }$$

and the function f_l(x) is given by

$$\begin{equation} f_l(x)=\int ^x_0 F(y)y^{(p-l)/2}\,dy, \end{equation} \tag{ 4 }$$

where F(y) is an integration over Bessel functions (Rybicki & Lightman 1979). The temporal indices in the case of a uniform density medium are α_F = 11/10, α_τ = 1 − 3p/2, and α_m = −3. The normalizations are such that F₀ and τ₀ are the flux density and optical depth at a frequency of ν = 1 Hz at t = t₀, and ν₀ is the synchrotron peak frequency in the rest frame of burst at t = t₀. Furthermore, the synchrotron self-absorption frequency, ν_a, is defined by the condition τ_ν(ν_a) = 1.

We fit this synchrotron model to our radio data using F₀, τ₀, and ν₀ as free parameters. Since we have no a priori knowledge about the expected values of the synchrotron spectrum parameters we assume that they follow a flat distribution in log-space. We note that any further assumption about the distribution of these parameters will only serve to restrict the resulting energy distributions, and we therefore consider our assumed flat log-space distributions to be conservative. Indeed, using input distributions that are flat in linear space leads to energy distributions (for Group B) that are narrower by about 1–1.5 dex; a less severe effect is observed for Group A since the afterglow properties are better constrained. We also note that the use of flat linear distributions leads to lower median energies for Group B events by about a factor of three, but that this effect is much smaller than the overall allowed range of energies for each burst.

For the detected objects (Groups A and B) we retain all solutions that reproduce the measured flux density within the error bars, while for the non-detections we use 3σ as an upper bound. Moreover, for the bursts in Group A, ν₀ is constrained by the multi-frequency observations and no additional constraints are required. However, for the bursts in Groups B and C, which have only a single-frequency observation, we require that both ν_m and ν_a have values below the observing frequency since the light curves are always declining at the time of our observations.³

Using the allowed ranges of F₀, τ₀, and ν₀ we determine the set of relevant physical parameters: n_e, γ_m, and B, where B is the magnetic field strength. The radius of the blast wave, r, remains unconstrained (e.g., Frail et al. 2000; Berger et al. 2004):

$$\begin{equation} B=11.7(p+2)^{-2}F^{-2}_{0,-52}(r_{17}/d_{L,28})^4 \,\,\,{\rm G}, \end{equation} \tag{ 5 }$$

$$\begin{equation} \gamma _m=6.7(p+2)F_{0,-52}\nu ^{1/2}_{0,9}(r_{17}/d_{L,28})^{-2}, \end{equation} \tag{ 6 }$$

$$\begin{equation} n_e=3.6\times 10^{10}c_n\eta _1F^3_{0,-52}\nu ^{(1-p)/2}_{0,9} \tau _{0,32}r^{-1}_{17}(r_{17}/d_{L,28})^{-6} \,\,\,{\rm cm}^{-3}, \end{equation} \tag{ 7 }$$

$$\begin{equation} c_n=(1.67\times 10^3)^{-p}(5.4\times 10^2)^{(1-p)/2}(p+2)^2 /(p-1), \end{equation} \tag{ 8 }$$

where d_{L, 28} is the luminosity distance in units of 10²⁸ cm assuming the standard cosmological parameters (H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_M = 0.27, and Ω_Λ = 0.73), and η is the reciprocal of the thickness of the emitting shell.

The unknown radius of the blast wave can be constrained by introducing a relationship between E_ST = n_em_p(r/1.05)⁵[t_NR/(1 + z)]⁻² and the energy in the electrons and the magnetic field. We use the condition that at most half of the blast-wave energy is available for accelerating electrons and producing the magnetic field, i.e., (E_B + E_e) ≲ E_ST/2. The total energy in the accelerated electrons is E_e = [(p − 1)/(p − 2)]n_eγ_mc²V, while the energy in the magnetic field is E_B = B²V/8π, where V = 4πr²/η is the volume of the synchrotron emitting shell. The energy budget is minimized near equipartition (i.e., E_e ≈ E_B), and we use this constraint to determine the minimum required energy (e.g., Frail et al. 2000; Berger et al. 2004); this conclusion was verified with radio interferometric measurements of the size of GRB 030329 (Taylor et al. 2004; Frail et al. 2005).

Finally, using the inferred radius for each possible solution, we require for self-consistency that β  ≲  1, where β  =  2r(1 + z)/5ct. The resulting β distributions are shown in Figure 1. These results indicate that most of the detected bursts obey the self-consistency requirement, although we reject GRBs 000926, 020819, and 021004 from the analysis since ≳ 50% of their allowed solutions lead to relativistic velocities at the time of the observations. This does not rule out that the Sedov–Taylor solution is applicable, but simply indicates that additional observations are required to narrow down the range of allowed solutions. We furthermore find that the bulk of the upper limits do not rule out relativistic expansion, and we therefore do not use these limits in the energy distribution analysis below.

Zoom In Zoom Out Reset image size

The resulting solutions for each burst can be cast in terms of a two-dimensional parameter space in n_e versus E_ST. Thus, there is a degeneracy between the two parameters, in the sense that larger densities lead to lower energies. Clearly, the bursts in Group A, for which the peak of the synchrotron spectrum is well-determined, lead to the best constraints in this two-dimensional phase space. Indeed, as shown for GRB 980703 (Figure 2), the allowed range of energies for all solutions that reproduce the observed flux density is about 10⁴⁸–10⁵³ erg, while the solutions that also satisfy the requirements that (E_B + E_e) ≲ E_ST/2 span a much narrower range of about 10⁵¹–10⁵² erg, with a roughly log-normal distribution centered on log(E_ST) ≈ 51.6.

Zoom In Zoom Out Reset image size

The bursts with only single-frequency observations (Group B) cover a much larger area in the E_ST–n_e phase-space since only an upper bound can be placed on ν_m and ν_a. To further constrain the energy we place an additional conservative constraint on the density of n_e < 100 cm⁻³, motivated by the results of detailed broadband modeling that show n_e ∼ 0.1–10 cm⁻³ (e.g., Panaitescu & Kumar 2001, 2002; Yost et al. 2003). Given the anti-correlation between density and energy, our conservative limits lead to a wider range of allowed energies than if we choose a limit of n_e < 10 cm⁻³. An example of this additional constraint for a Group B burst (GRB 010921) is shown in Figure 2. We do not place a lower bound on the density, since for the phase-space of allowed solutions this would not lead to a significant change in the energy distribution. We stress that beyond placing an upper bound on n_e no constraints have been placed on the distributions of either n_e or β since both are inferred, and not input, parameters in our model.

As noted above, we do not consider the energies for the bursts in Group C since the upper limits generally allow a wide range of solutions that are not consistent with the Sedov–Taylor formulation. This indicates that the limits are generally not deep enough to provide a meaningful constraint on the energy. Future deep radio observations may provide much better constraints (see below).

The resulting energy probability distributions for the bursts in Groups A and B are shown in Figure 3. The median energy and 90% confidence range (i.e., 5%–95% of the distribution) for each burst are listed in Table 2. We include in these ranges the small subset of solutions that lead to β values in slight excess of 1 since these are at most mildly relativistic and furthermore do not significantly change the distributions (Figure 3). We find that varying the electron power-law index over the range p = 2.1–2.5 (e.g., Curran et al. 2008) leads to a change in the median energy of only 0.1–0.2 dex (compared to our fiducial value of p = 2.2), with larger values of p leading to lower median energies. Similarly, varying the magnetic energy fraction away from equipartition to _B = 0.1 and 0.01 leads to an increase in the median energy of about 0.25 and 0.5 dex, respectively. Both of these effects are much smaller than the overall spread in energy for each burst, but they do produce minor systematic trends.

Zoom In Zoom Out Reset image size

The combined distribution for the subset of 11 bursts whose solutions are generally self-consistent is shown in Figure 4. The median and 90% confidence ranges are 7 × 10⁵¹ erg and 1.1 × 10⁵⁰–3.3 × 10⁵³ erg, respectively.

Zoom In Zoom Out Reset image size

For completeness we also include the distribution of resulting densities, shown in Figure 5. The median and 90% confidence ranges are 13 cm⁻³ and 0.6–90 cm⁻³, respectively.

Zoom In Zoom Out Reset image size

The key results of our analysis are that the median energy for the 11 bursts with self-consistent solutions is E_K ≈ 7 × 10⁵¹ erg, while the 90% confidence range is 1.1 × 10⁵⁰–3.3 × 10⁵³ erg. The median value is about a factor of three times higher than previous calorimetric measurements for GRBs 970508, 980703, and 030329, for which energies of 3 × 10⁵¹, 3 × 10⁵¹, and 10⁵¹ erg, respectively, were determined (Berger et al. 2004; Frail et al. 2005).

Similarly, the inferred energies are somewhat larger than the distributions of beaming-corrected γ-ray and kinetic energies inferred from broadband multi-epoch studies (Figure 6). From various such analyses, the median γ-ray energy is 〈E_γ〉 ≈ 8 × 10⁵⁰ erg (Frail et al. 2001; Bloom et al. 2003; Friedman & Bloom 2005), while the median kinetic energy is 〈E_K〉 ≈ 5 × 10⁵⁰ erg (e.g., Panaitescu & Kumar 2001, 2002; Yost et al. 2003); see Figure 6. In both cases the 90% range spans about 2.5 orders of magnitude, somewhat narrower than our inferred 90% confidence range for E_ST. The extension to larger energies found in our analysis mainly reflects the lack of spectral peak determinations for the bursts with single-frequency observations (see Figure 3). These large energies can be generally eliminated with a measurement of the synchrotron peak in the GHz frequency range (e.g., Group A bursts; Figure 6).

Zoom In Zoom Out Reset image size

In the context of our results, we note that recent numerical work by Zhang & MacFadyen (2009) led these authors to conclude that the timescale to reach isotropy is ∼10² yr rather than ∼1 yr as indicated by the analytic formulation of Livio & Waxman (2000), which we follow here. As a result, they note that using the Sedov–Taylor formulation may lead to an erroneous estimate of the kinetic energy. However, inspection of the resulting potential discrepancies reveals that this effect is at most a factor of two as long as self-consistency between the inferred energy and density and the transition to the Sedov–Taylor phase is ensured (see their Figure 10). The discrepancies become larger if the wrong timescale is assumed for the transition to non-relativistic expansion, but this quantity is not a free parameter. Indeed, our distributions of β values point to self-consistency for most bursts, and allow us to reject objects that are potentially still relativistic. Since the potential systematic uncertainty of about a factor of two is significantly smaller than the overall spread in allowed energy for each burst, we do not consider this to be an obstacle to our analysis, or to future work on the energy scale using late-time radio measurements.

As clearly demonstrated in Figure 3, the most constrained energy determinations require a measurement of the synchrotron spectral peak (Group A); the absence of such a constraint requires additional assumptions about the circumburst density and results in a much wider energy range. Indeed, this is the key reason for the wider range of allowed high energy solutions (≳ 10⁵² erg) compared to the results for E_γ and E_K (Figure 6). Observations of GRBs 970508, 980703, and 030329 demonstrate that the spectral peak is typically located at ∼few GHz on a timescale of ∼150 days. Thus, observations in the 1–10 GHz range on a timescale of ∼few hundred days should allow us to determine the peak flux and frequency. This will in turn provide an energy estimate with a similar level of precision to the results of early-time broadband modeling.

This is a fortuitous conclusion since with the full frequency coverage of the Expanded VLA (EVLA) it will soon be possible to cover this entire range in a few hours of observations to a sensitivity that is about an order of magnitude better than the VLA. As we demonstrated here, such a modest investment of observing time (2–3 hr per burst) can yield a robust estimate of the GRB energy distribution, regardless of the ability to measure jet opening angles. Pursuing these observations for all bursts with a measured redshift will require only ∼50–100 hr of EVLA time per year. Indeed, with such observations we should be able to constrain the energy distribution to a comparable level as existing studies within a single year given that about 30 GRBs with known redshifts just from 2009 are now available for EVLA observations (a similar sample is available from 2008 bursts). In the longer term, the large number of objects will allow us to test the energy distribution as a function of redshift, at least over the range z ∼ 1–3 where the bulk of the detected bursts occur (Berger et al. 2005; Jakobsson et al. 2006). Similarly, this approach will be particularly useful for bursts that lack detailed optical or X-ray light curves due to observational constraints or dust extinction, and for bursts with controversial estimates of the jet opening angles.

We thank Dale Frail and Eli Waxman for helpful discussions and comments on the manuscript, and the anonymous referee for helping improve the quality of the manuscript.