GRB 120521C AT z ∼ 6 AND THE PROPERTIES OF HIGH-REDSHIFT γ-RAY BURSTS

Long duration γ-ray bursts (GRBs) are known to be associated with the violent deaths of massive stars (e.g., Woosley & Bloom 2006). In conjunction with the large luminosities of their afterglows, they can therefore serve as powerful probes of the high-redshift universe (Inoue et al. 2007), providing clues to the formation environments of the first stars, the ionization and metal enrichment history of the universe, and the properties of galaxies that are otherwise too faint to study through direct imaging and spectroscopy (Totani et al. 2006; Tanvir et al. 2012; Chornock et al. 2013). Furthermore, modeling of multi-wavelength afterglow data allows us to constrain the densities and structure of massive star environments on parsec scales, as well as the energies of the explosions and the degree of ejecta collimation.

To use GRBs as effective probes of star-formation in the re-ionization era (z ≳ 6; Fan et al. 2002, 2006), it is important to understand whether there is any evolution in the properties of their progenitors with redshift. This is best achieved by studying the afterglows of the highest-redshift events to determine their explosion energy, circumburst density and degree of collimation, and by comparing these properties with those of their lower-redshift counterparts. In the long term, such studies have the potential to uncover the contribution of Population III stars, which have been speculated to be highly energetic (E_iso ∼ 10⁵²–10⁵⁷ erg) with relatively long durations (T₉₀ ∼ 1000 s; e.g., Fryer et al. 2001; Bromm et al. 2003; Heger et al. 2003; Mészáros & Rees 2010; Suwa & Ioka 2011; Toma et al. 2011; Wang et al. 2012).

At present, there are only three GRBs with spectroscopically confirmed redshifts of z ≳ 6: GRB 050904 at z = 6.29 (Tagliaferri et al. 2005; Haislip et al. 2006; Kawai et al. 2006), GRB 080913 at z = 6.70 (Greiner et al. 2009), and GRB 090423 at z = 8.23 (Salvaterra et al. 2009; Tanvir et al. 2009). In addition, GRB 090429B has an inferred photometric redshift of z ∼ 9.4 (Cucchiara et al. 2011). To fully determine the physical properties of a GRB and its environment requires multi-wavelength observations spanning the radio through to the X-rays; only two of the z ≳ 6 events have radio detections: GRB 050904 (Frail et al. 2006; Gou et al. 2007) and GRB 090423 (Tanvir et al. 2009; Chandra et al. 2010).

Previous studies of GRB 050904 have found a high circumburst density (n ∼ 10²–10³ cm⁻³; Frail et al. 2006; Gou et al. 2007), a high isotropic-equivalent γ-ray energy (E_{γ, iso} ≈ 10⁵⁴ erg; Cusumano et al. 2006), a large isotropic-equivalent kinetic energy (E_{K, iso} ≈ few × 10⁵³ erg; Frail et al. 2006; Gou et al. 2007), and no evidence for host extinction (A_V ≲ 0.1 mag; Gou et al. 2007; Zafar et al. 2010, although see also Stratta et al. 2007, 2011). A jet break at t_jet ≈ 3 d (Tagliaferri et al. 2005) indicates a beaming-corrected γ-ray energy of 8 × 10⁵¹ erg and kinetic energy of E_K ≈ 2 × 10⁵¹ erg, the latter being one of the largest known (Gou et al. 2007). GRB 090423 has an inferred density of n ≲ 1 cm⁻³ (Chandra et al. 2010), large isotropic-equivalent γ-ray energy (E_γ ≳ 10⁵³ erg) and kinetic energy (E_{K, iso} ≳ 3 × 10⁵³ erg), and no host extinction (A_V ≲ 0.1 mag; Tanvir et al. 2009). No jet break was seen for this event, resulting in a claim of E_K ≳ 7 × 10⁵¹ erg, even larger than for GRB 050904.

Whereas individual studies of these two GRBs have been undertaken, they employed different implementations of afterglow synchrotron models and their results cannot be compared directly. Here we report multi-wavelength observations of GRB 120521C and deduce a photometric redshift of z ≈ 6, making this the third high-redshift GRB with multi-wavelength data from radio to X-rays. The availability of well-sampled light curves spanning several orders of magnitude in frequency and time allow us to perform broadband afterglow modeling, and thereby to determine the energetics of the explosion, the density profile of the circumburst environment, the microphysical parameters of the relativistic shocks, and the collimation of the ejecta. We additionally re-analyze all available afterglow data for GRBs 050904 and 090423, enabling us to compare the three high-redshift GRBs in a uniform manner. Finally, we compare the properties of the high-redshift GRBs to those of bursts at z ∼ 1 to investigate whether high-redshift GRBs exhibit evidence for an evolution in the progenitor population or favor different environments than their lower-redshift counterparts. We present our observations and analysis for GRB 120521C in Section 2 and determine a photometric redshift for this event in Section 3. We describe the theoretical model employed and our multi-wavelength modeling software in Section 4 and present our broadband afterglow model for GRB 120521C in Section 5. We apply our modeling code to re-derive the properties of GRBs 050904 and 090423 in Section 6 and compare the results to those obtained for GRB 120521C and to lower-redshift events in Section 7. We present our conclusions in Section 8. We use the standard cosmological parameters, Ω_m = 0.27, Ω_Λ = 0.73 and H₀ = 71 km s⁻¹ Mpc⁻¹. All magnitudes are in the AB system, unless stated otherwise.

GRB 120521C was discovered with the Swift Burst Alert Telescope (BAT; Barthelmy et al. 2005) on 2012 May 21 at 23:22:07 UT (Baumgartner et al. 2012). The burst duration was T₉₀ = (26.7 ± 0.4) s, with a fluence of F_γ = (1.1 ± 0.1) × 10⁻⁶ erg cm⁻² (15–150 keV; Markwardt et al. 2012). The Swift X-Ray Telescope (XRT; Burrows et al. 2005) began observing the field 69 s after the BAT trigger, leading to the detection of an X-ray afterglow at coordinates R.A.(J2000) = 14^h17^m0873, Decl.(J2000) = +42°08'410, with an uncertainty radius of 16 (90% containment).⁷ XRT continued observing the afterglow for 1.5 days in photon counting (PC) mode, with the last detection at about 0.5 days.

2.1. X-Rays

2.2. Optical and Near-IR

We obtained riz-band imaging of the XRT error circle beginning about 40 min after the BAT trigger using ACAM on the William Herschel Telescope (WHT) and MOSCA on the Nordic Optical Telescope (NOT). We analyzed the data using standard procedures within IRAF⁸ and astrometrically aligned and photometrically calibrated the images using Sloan Digital Sky Survey (SDSS) stars in the field. We found a brightening point source in the WHT z-band images within the revised XRT error circle at the position R.A.(J2000) = 14^h17^m0882, Decl.(J2000) = +42°08'416, with z = 23.5 ± 0.3 mag⁹ (at Δt ≈ 0.04 d), i ≳ 23.8 mag (3σ), and r ≳ 24.3 mag (3σ; Table 2).

Table 2. Optical and Near-infrared Observations of GRB 120521C

Δt	Telescope	Instrument	Band	Frequency	Flux densitya	Uncertaintya	Detection?
(days)				(Hz)	(mJy)	(mJy)	(1 = Yes)
0.0316	WHT	ACAM	R	4.81e+14	0.000585	0.000195	0
0.0372	WHT	ACAM	I	3.93e+14	0.00109	0.000362	0
0.0379	NOT		R	4.81e+14	0.000702	0.000234	0
0.0405	WHT	ACAM	z	3.46e+14	0.00146	0.000408	1
0.0433	NOT		I	3.93e+14	0.00135	8.00e-05	0
0.106	WHT	ACAM	z	3.46e+14	0.00444	0.000555	1
0.108	WHT	ACAM	z	3.46e+14	0.00369	0.000669	1
0.109	WHT	ACAM	z	3.46e+14	0.00476	0.000615	1
0.111	WHT	ACAM	z	3.46e+14	0.0036	0.000625	1
0.112	WHT	ACAM	z	3.46e+14	0.00402	0.000651	1
0.115	WHT	ACAM	z	3.46e+14	0.00313	0.000717	1
0.117	WHT	ACAM	z	3.46e+14	0.00398	0.000653	1
0.119	WHT	ACAM	z	3.46e+14	0.00253	0.000748	1
0.12	WHT	ACAM	z	3.46e+14	0.00408	0.000635	1
0.122	WHT	ACAM	z	3.46e+14	0.0031	0.000725	1
0.124	WHT	ACAM	z	3.46e+14	0.00301	0.000649	1
0.126	WHT	ACAM	z	3.46e+14	0.00300	0.00068	1
0.208	PAIRITEL		K	1.37e+14	0.255	0.0848	0
0.208	PAIRITEL		H	1.84e+14	0.0932	0.031	0
0.208	PAIRITEL		J	2.38e+14	0.0633	0.0211	0
0.282	UKIRT	WFCAM	K	1.37e+14	0.0125	0.00134	1
0.318	UKIRT	WFCAM	J	2.38e+14	0.0112	0.00108	1
0.321	Gemini-North	GMOS	z	3.46e+14	0.00632	0.000316	1
0.324	Gemini-North	GMOS	z	3.46e+14	0.00664	0.000332	1
0.326	Gemini-North	GMOS	z	3.46e+14	0.00659	0.000329	1
0.329	Gemini-North	GMOS	z	3.46e+14	0.00601	0.000301	1
0.332	Gemini-North	GMOS	z	3.46e+14	0.00686	0.000343	1
0.334	Gemini-North	GMOS	z	3.46e+14	0.00627	0.000313	1
0.336	Gemini-North	GMOS	z	3.46e+14	0.00623	0.000311	1
0.339	Gemini-North	GMOS	z	3.46e+14	0.00553	0.000277	1
0.341	Gemini-North	GMOS	z	3.46e+14	0.00647	0.000323	1
0.344	Gemini-North	GMOS	z	3.46e+14	0.00604	0.000302	1
0.347	Gemini-North	GMOS	z	3.46e+14	0.00593	0.000296	1
0.349	Gemini-North	GMOS	z	3.46e+14	0.00594	0.000297	1
0.352	Gemini-North	GMOS	z	3.46e+14	0.00619	0.00031	1
0.354	Gemini-North	GMOS	z	3.46e+14	0.00569	0.000284	1
0.356	UKIRT	WFCAM	H	1.84e+14	0.0126	0.00135	1
0.514	Keck	LRIS	g	6.29e+14	0.000114	3.8e-05	0
0.516	Keck	LRIS	I	3.93e+14	0.000453	0.000151	0
0.579	Gemini-North	GMOS	I	3.93e+14	0.000495	0.000165	0
0.586	Gemini-North	GMOS	z	3.46e+14	0.00433	0.000374	1
1.05	WHT	ACAM	z	3.46e+14	0.00191	0.000108	1

Note. ^aNot corrected for Galactic extinction.

Download table as:  ASCIITypeset image

Given the red color of the afterglow, r − z ≳ 0.8 mag, we considered this to be a possible high redshift source, and thus triggered a sequence of optical and infrared imaging with the Gemini-North Multi-Object Spectrograph (GMOS) on Gemini-North (iz), the Low Resolution Imaging Spectrometer (LRIS) on the W. M. Keck telescope (gI), and the Wide-Field Camera (WFCAM) on the United Kingdom Infrared Telescope (UKIRT; JHK). We reduced the data in the standard manner, using the instrument pipelines for GMOS and WFCAM. We performed aperture photometry using the Graphical Astronomy and Image Analysis tool (GAIA). We placed the aperture with reference to the GMOS z-band image with the highest signal-to-noise detection of the afterglow, and used an aperture size appropriate to the seeing FWHM. We determined the level and variance of the sky background from about 20 same-sized apertures placed on sky regions proximate to the burst location. We calibrated the optical photometry to SDSS and the JHK photometry using Two Micron All Sky Survey stars in the field.

We detected the afterglow in all filters redward of z-band and obtained non-detections with deep limits in the optical filters (gri) at the level of F_ν ≲ 0.45 μJy (3σ; Figure 1 and Table 2). On the other hand, the infrared colors were relatively blue: J − H = 0.13 ± 0.21 mag and J − K = 0.12 ± 0.21 mag. This suggested that reddening due to dust was negligible, and that the red r − z color was due to the Lyα break falling within the z-band, implying a photometric redshift of z ∼ 6. We perform a full analysis to determine a photometric redshift in Section 3.

Zoom In Zoom Out Reset image size

The Swift UV/Optical Telescope (UVOT) began observing the field 77 s after the burst. No optical counterpart was detected at the location of the X-ray afterglow (Oates & Baumgartner 2012). We performed photometry using the HEASOFT task uvotsource at the location of the NIR afterglow, and report our derived upper limits in Table 3.

Table 3. Swift UVOT Observations of GRB 120521C

Δt	Filter	Frequency	3σ Flux Upper Limita
(days)		(Hz)	(mJy)
1.5907e-02	B	6.9250e+14	2.8387e-02
1.6055e-02	UVM2	1.3450e+15	1.3113e-02
1.4277e-02	U	8.5630e+14	9.5506e-03
1.6770e-02	V	5.5500e+14	5.4590e-02
1.7337e-02	UVW1	1.1570e+15	9.5866e-03
1.6487e-02	UVW2	1.4750e+15	8.7097e-03
1.3334e-02	WHITE	8.6400e+14	3.7121e-03
1.0321e-01	B	6.9250e+14	1.3541e-02
7.4747e-02	UVM2	1.3450e+15	1.0309e-02
1.4464e-01	U	8.5630e+14	8.6864e-03
2.0687e-01	V	5.5500e+14	4.6748e-02
1.4259e-01	UVW1	1.1570e+15	3.8169e-03
2.0528e-01	UVW2	1.4750e+15	1.9738e-03
1.0784e-01	WHITE	8.6400e+14	2.9149e-03
5.7824e-01	WHITE	8.6400e+14	8.9244e-04
1.5257e+00	UVM2	1.3450e+15	1.5728e-03

Note. ^aNot corrected for Galactic extinction.

Download table as:  ASCIITypeset image

We obtained spectroscopic observations of the afterglow with Gemini-North/GMOS beginning 1.03 d post-burst for a total exposure of 3600 s, by which time the source had faded to z ≈ 23.2 mag. We used the R400 grism and a slit width of 1'', providing a wavelength coverage of 5850–10140 Å and a resolution of R ≈ 1900. The data were reduced using the GMOS pipeline. A faint trace of the afterglow was visible at the red end of the spectrum. The trace disappears around 8700 Å, which unfortunately coincides with the gap between the GMOS CCDs. Assuming this break is due to Lyα, we deduce z ≈ 6.15, consistent with the red r − z color. We plot the extracted spectrum in Figure 2, adaptively re-binned to produce approximately the same noise in each bin.

Zoom In Zoom Out Reset image size

2.3. Radio

We observed GRB 120521C with the Karl G. Jansky Very Large Array (VLA) beginning on 2012 May 22.12 UT at mean frequencies of 5.8 GHz (lower and upper sideband frequencies set at 4.9 and 6.7 GHz, respectively) and 21.8 GHz (lower and upper sideband frequencies of 19.1 and 24.4 GHz, respectively). We employed 3C286 as a flux and bandpass calibrator and interleaved observations of J1419+3821 for calculating time-dependent antenna gains. All observations utilized the VLA WIDAR correlator (Perley et al. 2011). We excised radio frequency interference from the data, resulting in final effective bandwidths of ≈1.5 GHz at 5.8 GHz and ≈1.75 GHz at 21.8 GHz. We performed all data calibration and analysis with the Astronomical Image Processing System (AIPS; Greisen 2003) using standard procedures for VLA data reduction.

In our first epoch at 21.8 GHz (0.15 d after the burst), we did not detect any significant radio emission within the refined Swift XRT error circle to a 3σ limit of 50 μJy (Table 4). However, we detected a radio source in the second epoch at 1.15 d after the burst (Figure 3). This source subsequently faded, confirming it as the radio afterglow. We also detected the afterglow at 6.7 GHz in our observations taken between 4.25 and 29.25 d after the burst; however, we did not find significant radio emission at 4.9 GHz. We treat these two side-bands separately in our analysis, but for simplicity, we show side-band averaged images in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. VLA Observations of GRB 120521C

Δt	VLA	Frequency	Integration Time	Integrated Flux	Uncertainty	Detection?
(days)	Configuration	(GHz)	(min)	density (μJy)	(μJy)	(1 = Yes)
0.15	CnB	4.9	15.28	41.7	13.9	0
		6.7	15.28	48.0	16.0	0
		21.8	15.07	50.7	16.9	0
1.15	CnB	4.9	10.12	51.0	17.0	0
		6.7	10.12	57.3	19.1	0
		21.8	15.07	112	18.5	1
4.25	B	4.9	15.27	41.1	13.7	0
		6.7	15.27	54.5	14.3	1
		21.8	14.52	66.5	18.6	1
7.25	B	4.9	15.12	39.9	13.3	0
		6.7	15.12	48.8	14.2	1
		21.8	12.97	65.8	18.3	1
12.27	B	21.8	32.95	30.6	10.2	0
14.27	B	21.8	32.68	38.4	12.8	0
13.27a	B	21.8	...	26.2	9.2	1
29.25	B	4.9	24.87	35.7	11.9	0
		6.7	24.87	29.1	9.7	1
174.66	A	4.9	46.43	28.5	9.5	0
	A	6.7	46.43	23.4	7.8	0

Note. ^aWeighted sum of data at 12.27 and 14.27 d.

Download table as:  ASCIITypeset image

We used the AIPS task JMFIT to determine the positional centroid and integrated flux of the radio afterglow in each epoch by fitting a Gaussian at the position of the source and fixing the source size to the restoring beam shape. The weighted mean position of the source, determined by combining all 21.8 GHz detections is R.A.(J2000) = 14^h17^m08803 ± 0002, Decl.(J2000) = +42°08' 4121 ± 003 (1σ). We summarize the results of the radio observations in Table 4. GRB 120521C was also observed by the Arcminute Microkelvin Imager Large Array at 15.75 GHz (AMI-LA; Staley et al. 2013), and we include the reported upper limits in our analysis.

To determine a photometric redshift, we interpolate the optical and NIR observations to a common time. To minimize this interpolation, we select a time of 8.1 hr after the burst when we obtained near-simultaneous zJHK photometry. We perform a weighted sum of the GMOS z-band observations at 7.7 hr <Δt < 8.5 hr and find F_ν = 6.22 ± 0.05 μJy at Δt ≈ 8.1 h. Since the NIR light curves are not well-sampled before 1 d, we use the z-band light curve to extrapolate the NIR fluxes. We first fit the z-band light curve with a broken power-law of the form $F_{\nu } = F_{\rm b} (({(t/t_{\rm b})^{-s\alpha _1}+(t/t_{\rm b})^{-s\alpha _2}})/{2} )^{-1/s}$, where t_b is the break time, F_b is the flux at the break time, α₁ and α₂ are the temporal decay rates before and after the break, respectively, and s is the sharpness of the break.¹⁰ We use the Python function curve_fit to estimate these model parameters and the associated covariance matrix. Our best-fit parameters are t_b = (0.34 ± 0.07) d, F_b = 6.89 μJy, α₁ = 0.83 ± 0.31, α₂ = −1.38 ± 0.43, and s = 1.7 ± 1.6 (Figure 5). Using this model to extrapolate the JHK photometry, we obtain F_ν = 11.1 ± 1.1 μJy, 12.8 ± 1.4 μJy, and 12.4 ± 1.3 μJy, at J, H, and K band, respectively, at the common time of 8.1 hr. The uncertainties are statistical only and do not include the systematic uncertainties introduced by the interpolation, which are less than 2%.

Zoom In Zoom Out Reset image size

After obtaining NIR fluxes at a common time, we build a composite model for the afterglow spectral energy distribution (SED). We use a sight-line-averaged model for the optical depth of the intergalactic medium (IGM) as described by Madau (1995), accounting for Lyα absorption by neutral hydrogen along the line of sight and photoelectric absorption by intervening systems. We also include Lyα absorption by the host galaxy, for which we assume a column of log (N_H/cm⁻²) = 21.1, the mean value for GRBs at z ∼ 1 (Fynbo et al. 2009). The free parameters in our model are the redshift of the GRB, the extinction along the line of sight within the host galaxy (A_V), and the spectral index (β) of the afterglow SED, F_ν∝ν^β. In order to not bias our results, we assume a uniform prior for the redshift and the extinction. We further use the distribution of extinction-corrected spectral slopes, β_ox from Greiner et al. (2011) as a prior on β. We use a Markov Chain Monte Carlo (MCMC) algorithm to explore the parameter space, integrating the model over the filter bandpasses and computing the likelihood of the model by comparing the resulting fluxes with the observed values. Details of our MCMC implementation are described in Section 4.2.

We find $z=5.93^{+0.11}_{-0.14}$, $\beta =-0.16^{+0.34}_{-0.25}$, and $A_{\rm V}=0.11^{+0.22}_{-0.10}$ mag, where the uncertainties correspond to 68% credible intervals about the median.¹¹ The parameters of the highest-likelihood model are z = 6.03, β = −0.34, and A_V = 0 mag, consistent with the 68% credible intervals derived from the posterior density functions (Table 5). We note that the median values differ from the highest-likelihood values. This is a standard feature of Monte Carlo analyses whenever the likelihood function is asymmetric about the highest-likelihood point. In this case, this occurs because the extinction is constrained to be positive, resulting in a truncation of parameter space. The best-fit model and a model with the median parameters are plotted in Figure 6, while the full posterior density function for the redshift is shown in Figure 7. We can rule out a redshift of z ≲ 5.6 at 99.7% confidence. The corresponding 99.7% confidence upper limit is z ≲ 6.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. Parameters from Optical/NIR SED Modeling of GRB 120521C

Parameter	Best-fit	68% Credible Regions
z	6.03	5.93$^{+0.11}_{-0.14}$
β	−0.34	$-0.16^{+0.34}_{-0.25}$
AV	0	0.11$^{+0.22}_{-0.10}$

Download table as:  ASCIITypeset image

We note that this constraint on the redshift relies on the assumed prior for β. Using broadband modeling we can locate the synchrotron break frequencies (explained in the next section) and thereby constrain β independent of the redshift. Therefore, in the subsequent multi-wavelength modeling we leave the redshift as a free parameter and fit for it along with the parameters of the explosion. For the optical and NIR frequencies, we integrate the model over the filter bandpasses to take into account absorption by the intervening IGM and the interstellar medium (ISM) of the host galaxy.

4.1. Synchrotron Model

In the standard synchrotron model of GRB afterglows, the SED consists of multiple power-law segments delineated by "break-frequencies," namely, the synchrotron cooling frequency (ν_c), the typical synchrotron frequency (ν_m), and the self-absorption frequency (ν_a). The location and evolution of these break frequencies, and the overall normalization of the spectrum depend upon the physical parameters of the explosion: the energy (E_{K, iso}), the circumburst density (n₀, or the normalized mass-loss rate in a wind environment, A_*), the power-law index of the electron energy distribution (p), the fraction of the blastwave energy transferred to relativistic electrons (_e) and to the magnetic fields (_B), and the half-angle of the collimated outflow (θ_jet). For further details of the synchrotron model, see Sari et al. (1998).

We have developed Python software for broadband modeling of GRB afterglows. Our software implements the full afterglow model with smoothly connected power law segments presented in Granot & Sari (2002, henceforth GS02). The model includes synchrotron cooling and self-absorption for both ISM and wind-like environments. The full treatment of the synchrotron model including local electron cooling results in five different spectral regimes with 11 definitions of the break frequencies, corresponding to different orderings of the synchrotron frequencies. Depending on the circumburst density profile and the combination of physical parameters, the spectrum evolves from fast cooling (ν_c < ν_m) to slow cooling (ν_c > ν_m), transitioning through the various spectral regimes (Figure 2 in GS02).

Given a set of explosion parameters, we compute the location of each of the 11 break frequencies using the expressions in GS02. Owing to slightly different normalizations of the break frequencies between the five spectral regimes, a sharp transition from one spectrum to another sometimes introduces discontinuities in the light curves. This is exacerbated by the fact that the transition times between spectra are not uniquely defined (see Table 3 in GS02). To overcome this and to establish a consistent framework, we add a linear combination of all spectra through which the spectrum evolves for a given set of physical parameters, with time-dependent weights. These weights are chosen such that each spectrum dominates in its own regime of validity, while allowing for the light curves to remain smooth when break frequencies cross each other at spectral transitions. A detailed description of our weighting scheme is provided in Appendix A.

The hydrodynamics presented in GS02 assume spherical expansion. While this is a good approximation in the early phase of the afterglow evolution when the Lorentz factor of the ejecta is $\Gamma \gg \theta _{\rm jet}^{-1}$ and only a small fraction of the jet is visible to an observer on Earth, deceleration of the jet to $\Gamma \lesssim \theta _{\rm jet}^{-1}$ results in a steep decline in the observed flux density at all frequencies at later times. We account for this "jet break" by changing the evolution of the break frequencies after the break time, t_jet, using the prescription in Sari et al. (1999), smoothing over the transition with a smoothing parameter¹² (for further discussion of the jet break based on numerical simulations, see van Eerten & MacFadyen 2012 and Leventis et al. 2013).

Our software also accounts for possible contributions in the optical and NIR from the host galaxy, as well as absorption and reddening of the afterglow light by dust in the host. For the former, we add the contribution of the host to the model afterglow light curve and fit for the flux density of the host in each waveband separately.¹³ For the latter, we use the Small Magellanic Cloud (SMC) extinction curve from Pei (1992) and fit for the B-band extinction in the rest frame of the host galaxy. We use the optical B-band rather than V-band to normalize our model, since the extinction curves of Pei (1992) are normalized in B-band. We find that using a Large Magellanic Cloud extinction model does not significantly affect the derived value of A_B, and we, therefore, use the SMC model throughout for consistency. We convert A_B to A_V using A_V = 0.83A_B (Pei 1992).

Radio observations can be strongly affected by scintillation, particularly at low frequencies (below ∼15 GHz). We account for scintillation in our modeling by calculating the modulation index (the expectation value of the rms fractional change in flux density) in the direction of the source and adding the expected flux variation in quadrature to the measured uncertainty. The details of our method are described in Appendix B.

We note that several observations, particularly those in the optical/NIR, have high signal-to-noise ratios approaching ∼50, implying photometry precise to the ∼2% level. However, the relative calibration of different instruments is generally not expected to be better than about 5%. In addition, the synchrotron model is by its nature a simplification of a complex physical process, and we, therefore, cannot expect the model to accurately represent the data at the ≲ 5% level. To account for this source of systematic uncertainty, we enforce a floor of 5% on the reported uncertainties prior to fitting.

To determine the best-fit solution, we compute the likelihood function using a Gaussian error model. The likelihood function for a data set composed of both detections and non-detections is given by (e.g., Lawless 2002; Helsel 2005)

$$\begin{equation} L = \prod p(e_i)^{\delta _i} F(e_i)^{1-\delta _i} \end{equation} \tag{ 1 }$$

where e_i are the residuals (the difference between the measurement or 3σ upper limit and the predicted flux from the model), δ_i is an indicator variable (equal to 0 for an upper limit and 1 for a detection), p(e_i) is the probability density function of the residuals, and F(e_i) is the cumulative distribution function of the residuals, equal to Prob(e_i ⩽ t) for a limit t. For a Gaussian error model,

$$\begin{equation} p(e_i) = \frac{1}{\sqrt{2\pi }\sigma }e^{-{e_i^2}/{2\sigma _i^2}}, \end{equation} \tag{ 2 }$$

where σ_i are the measurement uncertainties, while

$$\begin{equation} F(e_i) = \frac{1}{2} \left[1+{\it erf}\left(\frac{e_i}{\sqrt{2}\sigma _i}\right)\right], \end{equation} \tag{ 3 }$$

where erf(x) is the error function. We determine the best-fit parameters by maximizing the likelihood function using sequential least squares programming tools available in the Python SciPy package (Jones et al. 2001).

4.2. Markov Chain Monte Carlo

We employ the model and fitting algorithm described in Section 4 to determine the properties of GRB 120521C. The X-ray light curve displays a steep decline before ∼0.01 d, followed by a plateau phase extending to 0.25 d, neither of which can be described by the standard paradigm of the Blandford–McKee model (Blandford & McKee 1976). Such behavior is ubiquitous in the X-ray light curves of GRBs (e.g., Nousek et al. 2006; Margutti et al. 2013) and is usually attributed to the high-latitude component of the prompt emission (Kumar & Panaitescu 2000; Willingale et al. 2010) and energy injection (Nousek et al. 2006; Zhang et al. 2006; Dall'Osso et al. 2011), respectively. The models we employ only account for the emission from the afterglow blastwave shock. Therefore, we only utilize X-ray data after 0.25 d in the broadband fit.

In addition, the z-band light curve exhibits a peak at ∼8 hr. with a flux density of ≈7 μJy. If we interpret this peak as the passage of ν_m through the z-band, then ν_m should pass through 21.8 GHz at ≈200 d (evolving as t^−3/2, before a jet break) or at the very earliest around 40 days (evolving as t⁻², if we assume that a jet break occurred at 8 hours). In addition, the peak flux in the radio must be less than (in the wind model) or equal to (in the ISM model) the peak flux in optical/NIR. However, the 22 GHz radio light curve peaks before 10 d and all the radio observations are at a higher flux level than all of the optical and NIR detections. Thus, the optical/NIR and radio light curves are not compatible under the assumption that ν_m passes through z-band at 8 hr. Therefore, we do not include the z-band data before 0.25 d in our broadband fit. We return to the point of the X-ray and z-band light curves before 0.25 d in Section 5.1.

We find that an ISM model adequately explains all observations after ∼0.25 d (Figure 8). The spectrum remains in the slow cooling phase throughout, with the standard ordering of the synchrotron frequencies (ν_a < ν_m < ν_c) and with a peak flux density of F_{ν, m} ≈ 132 μJy. At Δt = 1 d, the synchrotron break frequencies are located at ν_m ≈ 5.5 × 10¹¹ Hz and ν_c ≈ 1.2 × 10¹⁶ Hz. We show the measured SED at 0.3 and 7.3 days in Figure 9, which highlights the importance of radio observations in constraining ν_m. The self-absorption frequency lies below the frequencies covered by our radio observations, ν_a ≲ 5 GHz and is therefore not fully constrained. Correspondingly, the physical parameters _e, _B, n₀, and E_{K, iso} exhibit degeneracies, with the unknown location of ν_a being the dominant source of uncertainty (Figure 10). Using the values of ν_m, ν_c, and F_{ν, max} from our best-fit model and the functional dependence of the microphysical parameters, _e, _B, n₀, and E_{K, iso} on the measured quantities ν_a ν_m, ν_c, and F_{ν, max}, we derive the following constraints: $\epsilon _{\rm e}\approx 0.15\nu _{\rm a, 9}^{{5}/{6}}$, $\epsilon _{\rm B}\approx 4.0\times 10^{-3}\nu _{\rm a, 9}^{{-5}/{2}}$, $n_{\rm 0}\approx 0.44\nu _{\rm a, 9}^{{25}/{6}}\,\rm {cm}^{-3}$, and $E_{\rm K, iso, 52}\approx 6.7\nu _{\rm a, 9}^{{-5}/{6}}$, where ν_{a, 9} is the self-absorption frequency in units of 10⁹ Hz. Imposing the restriction that _e be less than its equipartition value of 1/3, we can further restrict the self-absorption frequency to ν_a ≲ 2.7 × 10⁹ Hz. This allows us to place an upper bound on the circumburst density, n₀ ≲ 27 cm⁻³, and lower bounds on the isotropic equivalent energy, E_{K, iso, 52} ≳ 2.9 and _B ≳ 3.5 × 10⁻⁴. Similarly, imposing _B < 1/3, we can place lower bounds on the self-absorption frequency, ν_a ≳ 1.7 × 10⁸ Hz, the circumburst density, n₀ ≳ 2.8 × 10⁻⁴ cm⁻³, and _e ≳ 3.4 × 10⁻², and an upper bound on the isotropic equivalent energy, E_{K, iso, 52} ≲ 29. The parameters corresponding to the highest likelihood models are presented in Table 6, and the complete results of the Monte Carlo analysis are summarized in Table 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (45 panels)

Tarball of figure set images

Table 6. Best-fit Forward Shock Parameters

Parameter	120521Ca		090423a	050904
	ISM	Wind
z	6.04	5.70	8.23 (fixed)	6.29 (fixed)
p	2.12	2.03	2.56 (fixed)	2.07
e	$1.5\times 10^{-1}\nu _{\rm a,8}^{{5}/{6}}$ (3.4 × 10−2)	0.26	$1.3\times 10^{-1}\nu _{\rm a,8}^{{5}/{6}}$ (1.6 × 10−2)	9.1 × 10−3
b	$4.0\times 10^{-3}\nu _{\rm a,8}^{-{5}/{2}}$ (3.2 × 10−1)	2.7 × 10−3	$4.0\times 10^{-4}\nu _{\rm a,8}^{-{5}/{2}}$ (2.7 × 10−1)	2.0 × 10−2
n0	$4.4\times 10^{-1}\nu _{\rm a,8}^{{25}/{6}}$ (3.1 × 10−4)	...	$7.5\times 10^{-2}\nu _{\rm a,8}^{{25}/{6}}$ (2.4 × 10−6)	3.2 × 102
A*	...	0.81	...	...
EK, iso, 52 (erg)	$6.7\nu _{\rm a,8}^{-{5}/{6}}$ (2.9 × 101)	1.8	$7.2\times 10^{1}\nu _{\rm a,8}^{-{5}/{6}}$ (4.8 × 102)	2.4 × 102
tjet (d)	7.4	≳ 8b	16.7	1.5
θjet (deg)	2.3	≳ 10	2.5	5.4
AV (mag)	≲ 0.05	≲ 0.05	0.17	≲ 0.05
Eγ, iso (erg)	(1.9 ± 0.8) × 1053		(1.0 ± 0.3) × 1053c	(1.24 ± 0.13) × 1054d
Eγ (erg)	(1.5 ± 0.6) × 1050	≳ 2.9 × 1051	(9.5 ± 2.9) × 1049	(5.5 ± 0.6) × 1051
EK (erg)	$5.4\times 10^{49}\nu _{\rm a,8}^{-{5}/{6}}$	≳ 2.7 × 1050	$6.9\times 10^{50}\nu _{\rm a,8}^{-{5}/{6}}$	1.1 × 1052
Etot (erg)	1.8 × 1050e	≳ 3.2 × 1051	5.3 × 1051f	1.7 × 1052
$\eta _{\rm rad}=\frac{E_{\gamma }}{E_{\rm tot}}$	0.83	0.91	0.02	0.32

Notes. ^aThe best-fit values of the physical parameters, _e, _B, n₀, E_{K, iso} for GRBs 120521C and 090423 have been scaled to ν_{a, 8} = ν_a/10⁸ Hz. The values of these parameters corresponding to the highest likelihood model are given in parentheses and correspond to ν_a = 1.75 × 10⁸ Hz and ν_a = 8.6 × 10⁶ Hz for GRB 120521C and GRB 090423, respectively. ^bThe lower end of the 90% credible interval from MCMC simulations (see Table 7). The jet break time is not well constrained in the wind model for GRB 120521C. ^cvon Kienlin (2009). ^dAmati et al. (2008). ^eAssuming ν_a=1.75 × 10⁸ Hz, the best-fit value. ^fAssuming ν_a=8.6 × 10⁶ Hz, the best-fit value.

Download table as:  ASCIITypeset image

Table 7. Summary Statistics from MCMC Analyses

Parameter	120521C		090423	050904
	ISM	Wind
z	$6.01^{+0.05}_{-0.09}$	$5.71^{+0.04}_{-0.03}$	8.23 (fixed)	6.29 (fixed)
p	$2.17^{+0.09}_{-0.07}$	$2.05^{+0.04}_{-0.02}$	2.56 (fixed)	2.07 ± 0.02
e	$4.5^{+6.7}_{-2.4}\times 10^{-2}$	$0.20^{+0.09}_{-0.9}$	$2.7^{+2.0}_{-0.7}\times 10^{-2}$	$1.2^{+1.5}_{-0.5}\times 10^{-2}$
b	$0.7^{+1.5}_{-0.6}\times 10^{-2}$	$2.4^{+6.9}_{-1.7}\times 10^{-3}$	$4.8^{+9.5}_{-3.9}\times 10^{-2}$	$1.3^{+2.2}_{-1.1}\times 10^{-2}$
log n0	$-2.7^{+1.4}_{-1.0}$	...	$-4.6^{+1.1}_{-0.6}$	$2.8^{+1.1}_{-0.7}$
A*	...	$0.79^{+0.65}_{-0.44}$	...	...
EK, iso, 52 (erg)	$2.2^{+3.7}_{-1.4}\times 10^1$	$1.9^{+1.4}_{-0.9}$	$3.4^{+1.1}_{-1.4}\times 10^{2}$	$1.7^{+1.2}_{-1.0}\times 10^{2}$
tjet (d)	$6.8^{+3.8}_{-2.4}$	≳ 6a	$14.6^{+2.7}_{-2.3}$	$1.5^{+0.2}_{-0.1}$
θjet (deg)	$3.0^{+2.3}_{-1.1}$	≳ 9a	$1.5^{+0.7}_{-0.3}$	$6.2^{+3.3}_{-1.4}$
AV (mag)	<0.05	<0.05	0.15 ± 0.02	<0.05
Eγ, iso (erg)	(1.9 ± 0.8) × 1053		(1.0 ± 0.3) × 1053b	(1.24 ± 0.13) × 1054c
Eγ (erg)	$2.6^{+4.4}_{-2.0}\times 10^{50}$	≳ 2.1 × 1051	$3.2^{+2.7}_{-1.7}\times 10^{49}$	$7.4^{+4.8}_{-3.4}\times 10^{51}$
EK (erg)	$3.1^{+1.9}_{-0.9}\times 10^{50}$	≳ 5.2 × 1049a	$1.1^{+0.4}_{-0.2}\times 10^{51}$	$1.1^{+0.2}_{-0.2}\times 10^{52}$
Etot = Eγ + EK (erg)	6 × 1050	2 × 1051	1 × 1051	2 × 1052
$\eta _{\rm rad}=\frac{E_{\gamma }}{E_{\rm tot}}$	0.5	0.1d	0.03	0.4

Notes. ^aThe lower end of the 90% credible interval. The jet break time is not well constrained in the wind model for GRB 120521C. ^bvon Kienlin (2009). ^cAmati et al. (2008). ^dUsing isotropic-equivalent energies.

Download table as:  ASCIITypeset image

Our MCMC analysis allows us to constrain the redshift to $6.01^{+0.05}_{-0.09}$ (the full posterior density function is shown in Figure 6 as the blue histogram). This is consistent with the photometric redshift of $z=5.93^{+0.11}_{-0.14}$, which was based solely on the optical/NIR data and a prior on the spectral index (Section 3). At this redshift, the Swift/BAT γ-ray fluence, F_γ = (1.1 ± 0.1) × 10⁻⁶ erg cm⁻², corresponds to an isotropic energy release of E_{γ, iso} = (6.6 ± 0.6) × 10⁵² erg (104–1040 keV observer frame). Since this burst was not observed by any wide-band γ-ray satellite, we do not have information about its γ-ray spectrum outside the Swift 15–150 keV band. We, therefore, use an average K-correction based on the observed Swift/BAT fluence and computed 1–10⁴ keV rest-frame isotropic-equivalent γ-ray energies of the other z ≳ 6 GRBs: 050904, 080913, and 090423 (Sakamoto et al. 2005; Stamatikos et al. 2008; Pal'Shin et al. 2008; Palmer et al. 2009; von Kienlin 2009; Amati et al. 2008). We find that this K-correction ranges from a factor of about 1.8 (for GRBs 080913 and 090423) to 3.6 (for GRB 050904). We infer an approximate value of E_{γ, iso} = (1.9 ± 0.8) × 10⁵³ erg for GRB 120521C, where the range accounts for the uncertainty in the K-correction. Our best estimate of the kinetic energy from the Monte Carlo analysis is $E_{\rm K,iso}=(2.2^{+3.7}_{-1.4})\times 10^{53}$ erg, indicating that the radiative efficiency, $\eta _{\rm rad} = E_{\gamma , \rm iso}/(E_{\gamma ,\rm iso}+E_{\rm K, iso}) \approx 0.5$.

The 21.8 GHz radio light curve displays a plateau around 6 d at a flux level of f_{ν, m} ≈ 70 μJy (Figure 8). If we interpret this plateau as the passage of ν_m through the 21.8 GHz band, then we would expect ν_m to pass through 6.7 GHz at around 12 d with a comparable flux density and for the 21.8 GHz flux density to decline only modestly to about 50 μJy (evolving as t^{(1 − p)/2} ∼ t^−0.5). In addition, this would predict a flux density of 45 μJy at 6.7 GHz at the next epoch at Δt = 29.3 d. However, the 6.7 GHz light curve does not rise as expected, while the 21.8 GHz flux density plummets to about 26 μJy at Δt = 13.3 d. In addition, the 6.7 GHz observation at Δt = 29.3 yields a detection at barely 3 σ of 30 μJy. This behavior indicates a departure from isotropic evolution, and we find that a jet break at Δt ≈ 7 d adequately accounts for the radio observations after 10 days. The presence of a jet break means that the peak flux density of the broadband spectrum declines with time, while the break frequencies evolve faster; this explains why the 6.7 GHz flux density does not rise to the level observed at 21.8 GHz, and why the 21.8 GHz flux density rapidly declines following the plateau. Using the relation

$$\begin{equation*} \theta_{\rm jet} = 0.1\left(\frac{E_{\rm K, iso, 52}}{n_0}\right) ^{{1}/{8}} \left(\frac{t_{\rm jet}/(1+z)}{6.2\,\rm hr}\right)^{{3}/{8}} \end{equation*}$$

for the jet opening angle (Sari et al. 1999), and the distributions of E_{K, iso, 52}, n₀, z, and t_jet from our MCMC simulations (Figure 11), we find $\theta _{\rm jet} = 3.0^{+2.3}_{-1.1}$ degrees. Applying the beaming correction, E_γ = E_{γ, iso}(1 − cos θ_jet), we find $E_{\gamma } = (2.6^{+4.4}_{-2.0})\times 10^{50}$ erg. Similarly, the beaming-corrected kinetic energy is $E_{\rm K} = (3.1^{+1.9}_{-0.9})\times 10^{50}$ erg.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (10 panels)

Tarball of figure set images

The first radio detection in the 21.8 GHz band at Δt = 1.2 d (1.22 ± 0.02 mJy) is a factor of 2.7 times brighter than predicted by the model (0.45 ± 0.1 mJy, 1σ deviation from scintillation). Early-time excess radio emission in GRB afterglows has frequently been attributed to the presence of a reverse shock component (e.g., Kulkarni et al. 1999; Sari & Piran 1999; Berger et al. 2003; Soderberg & Ramirez-Ruiz 2003; Chandra et al. 2010; Laskar et al. 2013). We investigate the potential contribution of a reverse shock and derive an estimate for the Lorentz factor of the ejecta in Appendix C.

We also perform the Monte Carlo analysis detailed in Section 4.2 for a wind-like environment. The redshift distribution from the wind model is shown in Figure 7 as the green histogram. Our best-fit wind model is plotted in Figure 12. We find that the model matches the radio observations (including the first radio detection, which is missed by the ISM model), but under-predicts all X-ray data included as part of the fit. In this model, ν_a is constrained to lie between 7 and 22 GHz at Δt = 1.15 d, breaking the degeneracy encountered in the ISM model. We list the derived parameters in Tables 6 and 7. However, since the X-ray data are not fit well, we do not consider the wind model as an adequate representation of the data set.

Zoom In Zoom Out Reset image size

5.1. Potential Explanations for the z-band Peak at ≈ 8 hr

We now return to the peak in the z-band light curve at Δt ≈ 8 hr, which cannot be explained by the passage of the synchrotron peak frequency (see Section 5). One possible explanation for this peak is that the blastwave encounters a density jump, causing a long-lasting optical flare. Nakar & Granot (2007) showed that the greatest change expected in an optical light curve due to a density jump is bounded at Δα ≲ 1 (see also Gat et al. 2013), whereas the temporal behavior of the z-band flux density indicates a change of Δα ∼ 2.2. Hence, the z-band light curve is unlikely to be the result of an inhomogeneous external medium.

Another way to suppress the z-band flux before 8 hr is through absorption by neutral hydrogen in the vicinity of the progenitor. This is an attractive explanation in this case because the z-band straddles the Lyman break and the flux density in this band is therefore highly sensitive to small variations in the neutral hydrogen column along the line of sight. In particular, if the neutral hydrogen column were to decline with time due to destruction by the blastwave or by photo-ionization, it would lead to the observed behavior of the rising z-band flux density. Our first z-band detection is at ≈8 min in the rest-frame of the burst, corresponding to a distance of ∼1 AU from the progenitor, while the z-band peak occurs at ≈1.2 hr in the rest frame, corresponding to a distance of ∼8 AU. We find that an additional neutral hydrogen column of N_H ∼ 2 × 10²² cm⁻² at z = 6 would be sufficient to suppress the first z-band point to the observed flux level and the ionization of this column would therefore lead to the observed increase in flux. For a path length of ∼7 AU, this column corresponds to a density of ≈2 × 10⁸ cm⁻³ or a mass of about 10⁻⁷ M_☉ (assuming a spherical cloud). Although the requisite mass is not very large, the inferred density is four orders of magnitude higher than a typical molecular cloud in the Milky Way (Schaye 2001; McKee & Ostriker 2007). Thus, ionization of a large neutral hydrogen column along the line of sight is a feasible explanation for the rising z-band light curve only if the densities of molecular clouds at z ∼ 6 can be much greater than observed locally.

Another possible explanation for the initial rise in z-band is the injection of energy into the blastwave shock by slower-moving relativistic ejecta catching up with the decelerating blastwave. If the injection is rapid enough it could create a rising light curve at z-band, which would then be expected to break into a fading power-law if ν_m is located below z-band at the end of the injection phase. Energy injection has been frequently invoked to explain the plateau phase of GRB X-ray afterglows (e.g., Nousek et al. 2006; Zhang et al. 2006; Dall'Osso et al. 2011. The X-ray light curve of GRB 120521C indeed shows such a plateau at 0.01–0.25 d.

To test whether the X-ray and NIR light curves can result from energy injection, we use our ISM model as an anchor at Δt = t_end ≈ 8 hr, after which it is the best-fit model to the multi-wavelength data set (including the z-band and XRT observations). We then assume a period of energy injection between the start of the X-ray plateau at t_start ≈ few × 10⁻² d and t_end and use a simple power-law prescription for the energy as a function of time,

$$\begin{equation*} E_{\rm K, iso}(t) = \left\lbrace \begin{array}{@{}ll}E_{\rm K, iso,0}\left(\frac{t_{\rm start}}{t_{\rm end}}\right)^{\zeta }={{\rm const.}}, & t < t_{\rm start}\\ E_{\rm K, iso,0}\left(\frac{t}{t_{\rm end}}\right)^{\zeta }\propto t^{\zeta }, & t_{\rm start} < t < t_{\rm end} \\ E_{\rm K, iso,0} = {{\rm const.}}, & t > t_{\rm end}, \end{array}\right. \end{equation*}$$

where E_{K, iso, 0} is the total isotropic-equivalent blastwave kinetic energy after energy injection is complete. We note that the XRT light curve displays a steep decline before the plateau with α_X = −3.5 ± 0.2 at 90–345 s (Figure 8), which cannot be explained by the afterglow forward shock and is likely related to the prompt emission (see also Section 5). We therefore add an additional power-law component with a fixed slope of α_X = −3.5 to the model X-ray light curve.

We set E_{K, iso, 0} = 2.85 × 10⁵³ erg using our highest-likelihood model (values in parentheses in Table 6) and vary ζ, t_start, and t_end to obtain a good match to the X-ray and z-band light curves. We find that in general we are able to model either the X-ray plateau or the z-band rise, but not both. Our best simultaneous match to both light curves is shown in Figure 13 with the parameters, t_start ∼ 2.6 × 10³ s, t_end ∼ 1.9 × 10⁴ s, and ζ ∼ 1.25, corresponding to an increase in blastwave kinetic energy by a factor of (t_end/t_start)^ζ ∼ 12 over this period. Although the resulting light curves do not match the data perfectly, energy injection provides the most plausible explanation for the z-band peak. Finally, we note that there is some evidence for "flickering" in the form of statistically significant scatter about the overall z-band rise (Figure 5), but the observations do not sample these rapid time-scale flux variations well enough to allow us to comment on the nature or source of the variability.

Zoom In Zoom Out Reset image size

To place the physical properties of GRB 120521C derived above in the context of other high-redshift events, and to compare them in a uniform manner, we apply the above analysis to the other two GRBs at z ≳ 6 with radio to X-ray detections reported in the literature: GRB 050904 at z = 6.29 and GRB 090423 at z = 8.23.

6.1. GRB 050904

GRB 050904 was discovered with Swift/BAT on 2005 September 4 at 1:51:44 UT (Cummings et al. 2005). The burst duration was T₉₀ = 22.5 ± 10 s (Sakamoto et al. 2005), with a fluence of F_γ = (5.4 ± 0.2) × 10⁻⁶ erg cm⁻² (15–150 keV). A photometric redshift was reported by Tagliaferri et al. (2005) and Haislip et al. (2006), and spectroscopically confirmed by Kawai et al. (2006), making GRB 050904 the highest redshift GRB observed at the time.

We analyzed the XRT data for this burst in the same manner as described in Section 2.1. In our spectral modeling, we assume N_{H, MW} = 4.53 × 10²⁰ cm⁻² (Kalberla et al. 2005). The best-fit neutral hydrogen column density intrinsic to the host is $N_{\rm H, int} = 5.61^{+2.98}_{-2.44}\times 10^{22}$  cm⁻² (68% confidence intervals). In our temporally resolved spectral analysis, we find that the X-ray photon index is consistent with Γ = 2.03 ± 0.10 (68% confidence interval) for all XRT data following 490 s after the GRB trigger. We use this value of the photon index to convert the observed 0.3–10 keV light curve to a flux density at 1 keV. The X-ray data before 1.7 × 10³ s and at 3 × 10³ – 5 × 10⁴ s are dominated by multiple flares. We ignore XRT data in this time range in our analysis.

We compiled NIR observations of GRB 050904 in the Y, J, H, and K bands from the literature (Haislip et al. 2006; Gou et al. 2007), and corrected for Galactic extinction along the line of sight assuming E(B − V) = 0.061 mag (Schlafly & Finkbeiner 2011). Since z-band is located blueward of Lyman-α in the rest-frame of the GRB, flux within and blueward of this band is heavily suppressed by absorption by neutral hydrogen in the IGM and we do not include these bands in our multi-wavelength fit. This burst was observed over multiple epochs in the 8.46 GHz radio band with the VLA (Frail et al. 2006), and we use the individual observations and limits in our analysis. We list all photometry we use in our model in Table 8.

Table 8. Multi-wavelength Observations of GRB 050904

Δt	Telescope	Instrument	Band	Frequency	Flux densitya	Uncertaintya	Detection
(days)				(GHz)	(mJy)	(1σ, mJy)	(1 = Yes)
0.0202	Swift	XRT	1 keV	2.42e+17	0.00148	0.000485	1
0.128	SOAR		J	2.43e+14	0.184	0.00689	1
0.135	SOAR		J	2.43e+14	0.185	0.00695	1
0.142	SOAR		J	2.43e+14	0.146	0.00547	1
0.312	SOAR		J	2.43e+14	0.0555	0.00822	1
0.324	SOAR		Ks	1.37e+14	0.132	0.00882	1
0.408	UKIRT	WFCAM	H	1.82e+14	0.0566	0.00322	1
0.411	UKIRT	WFCAM	J	2.43e+14	0.0398	0.00226	1
0.424	UKIRT	WFCAM	K	1.37e+14	0.0755	0.00429	1
0.44	IRTF		K	1.37e+14	0.0646	0.00181	1
0.487	UKIRT	WFCAM	J	2.43e+14	0.0322	0.00214	1
0.505	VLA		X	8.46e+09	0.174	0.058	0
0.609	Swift	XRT	1 keV	2.42e+17	6.09e-05	2.37e-05	1
1.03	TNG	NICS	J	2.43e+14	0.0234	0.00322	1
1.09	VLT-UT1	ISAAC	J	2.43e+14	0.0171	0.000642	1
1.1	VLT-UT1	ISAAC	H	1.82e+14	0.0236	0.00157	1
1.12	SOAR		Y	2.91e+14	0.014	0.00379	1
1.12	VLT-UT1	ISAAC	Ks	1.37e+14	0.0324	0.00216	1
1.17	SOAR		J	2.43e+14	0.0139	0.00236	1
1.39	VLA		X	8.46e+09	0.075	0.025	0
1.91	Swift	XRT	1 keV	2.42e+17	5.37e-06	2.33e-06	1
2.09	VLT-UT1	ISAAC	J	2.43e+14	0.00797	0.00053	1
2.12	VLT-UT1	ISAAC	H	1.82e+14	0.0106	0.000706	1
2.15	VLT-UT1	ISAAC	Ks	1.37e+14	0.0144	0.000958	1
2.22	SOAR		J	2.43e+14	0.00929	0.00219	1
2.27	SOAR		Y	2.91e+14	0.00835	0.00307	1
3.1	VLT-UT1	ISAAC	J	2.43e+14	0.00345	0.000263	1
4.16	VLT-UT1	ISAAC	J	2.43e+14	0.00266	0.000204	1
5.32	VLT-UT1	ISAAC	J	2.43e+14	0.00166	0.000318	1
5.41	VLA		X	8.46e+09	0.075	0.025	0
6.22	VLA		X	8.46e+09	0.072	0.024	0
6.54	Swift	XRT	1 keV	2.42e+17	1.9e-06	8.59e-07	0
7.18	VLT-UT1	ISAAC	J	2.43e+14	0.000834	0.00167	0
20.1	VLA		X	8.46e+09	0.111	0.037	0
23.2	HST	NICMOS	F160W	1.82e+14	0.00013	2.5e-05	1
29.1	VLA		X	8.46e+09	0.09	0.03	0
33.4	VLA		X	8.46e+09	0.105	0.035	0
34.2	VLA		X	8.46e+09	0.069	0.023	0
35	VLA		X	8.46e+09	0.116	0.018	1
37.5	VLA		X	8.46e+09	0.067	0.017	1
44	VLA		X	8.46e+09	0.081	0.027	0

Notes. NIR observations are from Haislip et al. (2006), Gou et al. (2007), and Berger et al. (2007). Radio observations are from Frail et al. (2006). We report the Swift photometry included in our model (Figure 14). We do not use the Riz photometry in our model fitting (see Section 6.1) and do not list them here. ^aNot corrected for Galactic extinction.

Download table as:  ASCIITypeset image

As in previous studies of this burst (Frail et al. 2006; Gou et al. 2007), we find that an ISM model provides an adequate fit to the data. Our best-fit model is shown in Figure 14 and the corresponding physical parameters are listed in Table 6. The 8.5 GHz flux is severely suppressed by self absorption, with the self absorption frequency located around 280 GHz, above the characteristic synchrotron frequency, i.e., ν_m < ν_a. This requires a high-density circumburst environment, with n₀ ∼10³ cm⁻², while a jet break at ∼2 d is required to explain the sharp drop in the NIR light curves.

Zoom In Zoom Out Reset image size

Using MCMC analysis, we confirm the high density of the circumburst environment, $\log {(n_{\rm 0})} = 2.8^{+1.1}_{-0.7}$, with E_{K, iso}$=(1.7^{+1.2}_{-1.0})\times 10^{54}$ erg, _e$=(1.2^{+1.5}_{-0.5})\times 10^{-2}$, _B$=(1.3^{+2.2}_{-1.1})\times 10^{-2}$, and p = 2.07 ± 0.02. The values of all the parameters are consistent with those reported by Gou et al. (2007) within ∼2σ. We find a jet break time of $t_{\rm jet}=1.5^{+0.2}_{-0.1}$ d which is earlier than t_jet ∼ 3 d reported previously (Tagliaferri et al. 2005; Gou et al. 2007; Kann et al. 2007); however, our derived value of the jet opening angle, $\theta _{\rm jet}=6.2^{+3.3}_{-1.4}$ deg is consistent with the value reported by Gou et al. (2007), who also performed a full multi-wavelength analysis. We compare our derived posterior density functions for p, _e, _B, n₀, E_{K, iso}, and A_V directly with those reported by Gou et al. (2007) in Figure 15. Our distributions are similar, except that we find slightly smaller values for p. We note that we use different prescriptions for the synchrotron self-absorption frequency and evolution in the fast cooling regime. In addition, Gou et al. (2007) include the effects of inverse Compton losses, which we ignore in our model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (36 panels)

Tarball of figure set images

We find strong correlations between all four physical parameters (_e, _B, n₀, and E_{K, iso}; Figure 16). Detailed investigation using the analytical expressions for the spectra in terms of the spectral break frequencies given in GS02 reveals the cause to be multiple levels of degeneracy. For instance, the characteristic synchrotron frequency is not well constrained, since it is located below the frequencies covered by our radio observations at all times. At the same time, ν_a and the flux density at this frequency, F_{ν, a}, are not independently constrained, since this frequency lies below both the NIR and the X-rays. It is possible to change the two together in a way that leaves the NIR and X-ray light curves unchanged, without violating the radio limits. This latter degeneracy is the primary source of the observed correlations. We note that this degeneracy could have been broken with simultaneous detections in the radio and NIR.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (9 panels)

Tarball of figure set images

6.2. GRB 090423

GRB 090423 was discovered with Swift/BAT on 2009 April 23 at 7:55:19 UT (Krimm et al. 2009). The burst duration was T₉₀ = 10.3 ± 1.1 s (Palmer et al. 2009), with a fluence of F_γ = (5.9 ± 0.4) × 10⁻⁷ erg cm⁻² (15–150 keV). The afterglow was detected by Swift/XRT and ground-based near-infrared (NIR) follow-up observations, and the redshift, z = 8.26, was confirmed by NIR spectroscopy (Salvaterra et al. 2009; Tanvir et al. 2009). The burst was also observed with the Spitzer Space Telescope (Chary et al. 2009), the Combined Array for Research in Millimeter-wave Astronomy (CARMA; Chandra et al. 2010), the Plateau de Bure Interferometer (PdBI; Castro-Tirado et al. 2009; de Ugarte Postigo et al. 2012), the IRAM 30m telescope (Riechers et al. 2009), the Westerbrock Synthesis Radio Telescope (WSRT; van der Horst 2009), and the VLA (Chandra et al. 2010).

We analyzed XRT data for this burst using methods similar to GRB 050904 and GRB 120521C. We assume N_{H, MW} = 2.89 × 10²⁰ cm⁻² (Kalberla et al. 2005). The best-fit neutral hydrogen column density intrinsic to the host is $N_{\rm H, int} = (8.1^{+8.6}_{-6.5})\times 10^{22}$ cm⁻². In our temporally resolved spectral analysis, we find that the X-ray photon index is consistent with Γ = 2.03 ± 0.09 (68% confidence interval) for all XRT data following 260 s after the GRB trigger. We use this value of the photon index to convert the observed 0.3–10 keV light curve to a flux density at 1.5 keV (to facilitate comparison with Chandra et al. 2010). We compile all available photometry, together with our XRT analysis, in Table 9.

Table 9. Multi-wavelength Observations of GRB 090423

Δt	Telescope	Instrument	Band	Frequency	Flux Densitya	Uncertaintya	Detection
(days)				(GHz)	(mJy)	(1σ, mJy)	(1 = Yes)
0.0173	UKIRT	WFCAM	K	1.37e+14	0.0419	0.00209	1
0.0227	UKIRT	WFCAM	K	1.37e+14	0.0427	0.00213	1
0.0281	UKIRT	WFCAM	K	1.37e+14	0.0400	0.00200	1
0.0463	Swift	XRT	1.5 keV	3.63e+17	0.00126	0.000555	1
0.0475	Swift	XRT	1.5 keV	3.63e+17	0.000624	0.000262	1
0.0487	Swift	XRT	1.5 keV	3.63e+17	0.00121	0.000542	1
0.0498	Swift	XRT	1.5 keV	3.63e+17	0.00103	0.000452	1
0.051	Swift	XRT	1.5 keV	3.63e+17	0.000837	0.000363	1
0.0523	Swift	XRT	1.5 keV	3.63e+17	0.000886	0.000389	1
0.0536	Swift	XRT	1.5 keV	3.63e+17	0.000797	0.000351	1
0.0552	Swift	XRT	1.5 keV	3.63e+17	0.000702	0.000308	1
0.0566	Swift	XRT	1.5 keV	3.63e+17	0.000854	0.000375	1
0.0579	Swift	XRT	1.5 keV	3.63e+17	0.000982	0.000432	1
0.0593	Swift	XRT	1.5 keV	3.63e+17	0.000802	0.00035	1
0.0608	Swift	XRT	1.5 keV	3.63e+17	0.00104	0.000453	1
0.0622	Swift	XRT	1.5 keV	3.63e+17	0.000796	0.000346	1
0.0637	Swift	XRT	1.5 keV	3.63e+17	0.000765	0.000332	1
0.0644	Gemini-North	NIRI	J	2.38e+14	0.032	0.0016	1
0.0649	Swift	XRT	1.5 keV	3.63e+17	0.00137	0.000608	1
0.0661	Swift	XRT	1.5 keV	3.63e+17	0.000649	0.00028	1
0.0678	Swift	XRT	1.5 keV	3.63e+17	0.000515	0.000214	1
0.0695	Swift	XRT	1.5 keV	3.63e+17	0.000704	0.000302	1
0.0709	Swift	XRT	1.5 keV	3.63e+17	0.000883	0.000383	1
0.0727	Swift	XRT	1.5 keV	3.63e+17	0.000508	0.000217	1
0.075	Swift	XRT	1.5 keV	3.63e+17	0.000659	0.000248	1
0.0755	Gemini-North	NIRI	H	1.84e+14	0.0381	0.00191	1
0.115	Swift	XRT	1.5 keV	3.63e+17	0.000342	0.000141	1
0.118	Swift	XRT	1.5 keV	3.63e+17	0.000395	0.000165	1
0.121	Swift	XRT	1.5 keV	3.63e+17	0.000328	0.000132	1
0.124	Swift	XRT	1.5 keV	3.63e+17	0.000247	0.000103	1
0.128	Swift	XRT	1.5 keV	3.63e+17	0.00037	0.000151	1
0.131	Swift	XRT	1.5 keV	3.63e+17	0.000272	0.000108	1
0.134	Swift	XRT	1.5 keV	3.63e+17	0.000277	0.000114	1
0.14	Swift	XRT	1.5 keV	3.63e+17	0.000184	6.41e-05	1
0.185	Swift	XRT	1.5 keV	3.63e+17	0.000171	6.67e-05	1
0.19	Swift	XRT	1.5 keV	3.63e+17	0.000184	7.02e-05	1
0.195	Swift	XRT	1.5 keV	3.63e+17	0.000128	4.91e-05	1
0.202	Swift	XRT	1.5 keV	3.63e+17	0.00014	5.3e-05	1
0.209	Swift	XRT	1.5 keV	3.63e+17	0.000154	5.71e-05	1
0.363	Swift	XRT	1.5 keV	3.63e+17	6.28e-05	2.6e-05	1
0.384	PdBI			9.7e+10	0.24	0.0800	1
0.567	Swift	XRT	1.5 keV	3.63e+17	2.86e-05	1.12e-05	1
0.67	VLT-UT4	HAWKI	K	1.37e+14	0.0136	0.000681	1
0.696	ESO2.2m	GROND	K	1.37e+14	0.0115	0.00135	1
0.696	ESO2.2m	GROND	H	1.84e+14	0.0108	0.000825	1
0.696	ESO2.2m	GROND	J	2.38e+14	0.00865	0.000662	1
0.702	VLT-UT4	HAWKI	J	2.38e+14	0.00984	0.000492	1
0.781	ESO2.2m	GROND	K	1.37e+14	0.00769	0.00256	0
0.781	ESO2.2m	GROND	H	1.84e+14	0.0102	0.00068	1
0.781	ESO2.2m	GROND	J	2.38e+14	0.00789	0.000526	1
0.922	UKIRT	WFCAM	K	1.37e+14	0.0067	0.000922	1
0.934	Swift	XRT	1.5 keV	3.63e+17	1.8e-05	7.21e-06	1
1.29	PdBI			9.7e+10	0.24	0.0700	1
1.67	VLT-UT1	ISAAC	J	2.38e+14	0.00303	0.000546	1
1.87	CARMA			9.7e+10	0.54	0.18	0
2.21	VLA		X	8.46e+09	0.0927	0.0309	0
2.3	Swift	XRT	1.5 keV	3.63e+17	4.09e-06	1.72e-06	1
2.44	IRAM30m			2.5e+11	0.23	0.32	0
3.69	VLT-UT4	HAWKI	K	1.37e+14	0.00375	0.000213	1
3.72	VLT-UT4	HAWKI	J	2.38e+14	0.00211	0.000225	1
5.66	Swift	XRT	1.5 keV	3.63e+17	1.33e-06	6.04e-07	1
7.65	VLT-UT4	HAWKI	J	2.38e+14	0.0011	0.000366	0
8.29	PdBI			9.7e+10	0.24	0.0800	0
9.34	VLA		X	8.46e+09	0.0664	0.0114	1
14.3	VLA		X	8.46e+09	0.0437	0.0089	1
15.7	VLT-UT4	HAWKI	K	1.37e+14	0.00122	0.000406	0
20.7	VLA		X	8.46e+09	0.0422	0.0106	1
29.4	VLA		X	4.9e+09	0.044	0.025	0
33.1	VLA		X	8.46e+09	0.0496	0.011	1
46.3	Spitzer	IRAC	3.6 μm	8.4e+13	4.79e-05	1.3e-05	1
62	VLA		X	8.46e+09	0.0342	0.0116	0
279	Spitzer	IRAC	3.6 μm	8.4e+13	5.75e-05	1.92e-05	0

Notes. All observations except the Swift and Spitzer data are collected from Salvaterra et al. (2009); Tanvir et al. (2009); Chandra et al. (2010); Riechers et al. (2009); van der Horst (2009); and de Ugarte Postigo et al. (2012). ^aNot corrected for Galactic extinction.

Download table as:  ASCIITypeset images: 1 2

There are 134 ks of unpublished X-ray data in the Chandra archive for this GRB (PI: Garmire), taken between 16 and 42 d after the burst and distributed across five epochs. We downloaded and analyzed all available data from the Chandra archive. The GRB is marginally detected in three of the five epochs. We stacked observations taken close in time (epochs 1 and 2; epochs 3, 4, and 5) and restricted the energy range to 0.3–2 keV to increase the signal-to-noise ratio. The GRB is marginally detected in both stacks. We report the results of photometry using 15 apertures in Table 10. We convert the measured count rates into flux densities at 1.5 keV using the XRT photon index of Γ = 2.03.

Table 10. Chandra Observations of GRB 090423

Epoch	Δta	Exposure Time	Count Rateb	1.5 keV Flux
	(days)	(ks)	(10−4 s−1)	density (mJy)
1–2	16.8c	31.9	1.1 ± 0.6	(1.5 ± 0.8) × 10−7
3–5	37.8c	102.2	0.56 ± 0.26	(7.7 ± 3.6) × 10−8

Notes. ^aTime to mid-exposure. ^b0.3–2 keV, 15 (radius) aperture. ^cMean time since GRB, weighted by exposure time of individual epochs.

Download table as:  ASCIITypeset image

The 8.46 GHz radio light curve peaks at a similar flux density as does the NIR light curve, which strongly argues against a wind-like medium and suggests a constant-density environment. We also note that the millimeter observations reported in de Ugarte Postigo et al. (2012) are inconsistent with the forward-shock synchrotron model, since the flux densities at 97 GHz are much higher than in any other waveband, whereas the ISM model for GRB afterglows predicts that light curves at each frequency would reach the same peak flux density prior to the jet break. The millimeter data are shown in Figure 17 for completeness but have not been included in the analysis. This was also noted by Chandra et al. (2010), who suggested that the millimeter data and the first radio detection at 2.2 d possibly included emission from a reverse shock. We investigate this possibility further in Appendix D.

Zoom In Zoom Out Reset image size

Our best-fit model requires that the afterglow be in the slow cooling phase with the spectral ordering ν_a < ν_m < ν_c and a peak flux density of F_{ν, max} ≈ 142 μJy. At 1 day, the characteristic synchrotron frequency is ν_m ≈ 7.7 × 10¹² Hz, while the cooling break is in the X-rays, at 4.5 × 10¹⁷ Hz (1.8 keV). However, the data do not constrain ν_a. In the ISM model ν_a remains fixed before the jet break and falls as t^−0.2 after the jet break. Hence the only observational constraint on ν_a is that it is located below the radio band at all times. The model shown in Figure 17 is therefore only one of a family of models that match the data and have ν_a ≲ 8 GHz. Using the values of ν_m, ν_c and F_{ν, max} from our best-fit model and the functional dependence of the microphysical parameters, _e, _B, n₀, and E_{K, iso} on the measured quantities ν_a ν_m, ν_c, and F_{ν, max}, we derive the following constraints: $\epsilon _{\rm e}\approx 0.13\nu _{\rm a, 8}^{{5}/{6}}$, $\epsilon _{\rm B}\approx 4.0\times 10^{-4}\nu _{\rm a, 8}^{{-5}/{2}}$, $n_{\rm 0}\approx 7.5\times 10^{-2}\nu _{\rm a, 8}^{{25}/{6}}\,\rm {cm}^{-3}$, and $E_{\rm K, iso, 52}\approx 72\nu _{\rm a, 8}^{{-5}/{6}}$, where ν_{a, 8} is the self-absorption frequency in units of 10⁸ Hz. Imposing the theoretical restriction, _e < 1/3, we can further restrict the self-absorption frequency to ν_a ≲ 3.1 × 10⁸ Hz. This allows us to place an upper bound on the circumburst density, n₀ ≲ 8.3 cm⁻³, and lower bounds on the isotropic equivalent energy, E_{K, iso, 52} ≳ 28 and _B ≳ 2.4 × 10⁻⁵. Similarly, imposing _B < 1/3, we can place lower bounds on the self-absorption frequency, ν_a ≳ 6.8 × 10⁶ Hz and the circumburst density, n₀ ≳ 1.0 × 10⁻⁶ cm⁻³, and upper bounds on the isotropic equivalent energy, E_{K, iso, 52} ≲ 6.8 × 10² and _e ≳ 1.4 × 10⁻².

To further explore the degeneracies in the physical parameters of the explosion, we carried out an MCMC analysis with p fixed at our best-fit value of 2.56 (Figures 18 and 19). Our measured correlations between E_{K, iso}, n₀, _e, and _B are consistent with the expected analytic relations. We find a small amount of extinction within the host galaxy (A_V = 0.15 ± 0.02 mag), which is consistent with the low value of extinction (A_V ≲ 0.1 mag) inferred by other authors based on the X-ray and NIR observations alone (Tanvir et al. 2009; Zafar et al. 2011).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (28 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (10 panels)

Tarball of figure set images

A previous analysis of GRB 090423 claimed no jet break to ≈45 d (Chandra et al. 2010). However, our model requires a jet break at t_jet ≈ 15 d, driven by the late-time Spitzer 3.6 μm detection, as well as the radio non-detection at 62 d. In particular, ν_m passes through the radio band at 58 days while the radio light curve peaks at about 20 d, the signature of a jet break. In our model, the afterglow is optically thin at 8.46 GHz at all times. Following the jet break, the ν^1/3 part of the synchrotron spectrum transitions from t^1/2 to t^−1/3, followed by a transition to t^−p when ν_m crosses the radio band at 58 d, matching the observations. While the late-time Chandra data do not show an obvious break, the model with t_jet ≈ 15 d is consistent with the full X-ray light curve including the Chandra photometry and is required by the full model. From our MCMC analysis, we find $\theta _{\rm jet} = 1.5^{+0.7}_{-0.3}$ degrees (68% credible region). We list the best-fit parameters in Table 6 and the results of the MCMC analysis in Table 7.

Using the distribution of jet opening angles from our MCMC analysis and the isotropic-equivalent γ-ray energy, E_{γ, iso} = (1.03 ± 0.3) × 10⁵³ erg (von Kienlin 2009), we compute a beaming-corrected γ-ray energy of $E_{\gamma }=(3.2^{+2.7}_{-1.7})\times 10^{49}$ erg. The deduced value of the afterglow kinetic energy from the MCMC analysis is $E_{\rm K,iso} = (3.4^{+1.1}_{-1.4})\times 10^{54}$ erg, corresponding to a beaming-corrected energy of $E_{\rm K}=(1.1^{+0.4}_{-0.2})\times 10^{51}$ erg. Together, these results imply a low radiative efficiency, η ≡ E_γ/(E_K + E_γ) ∼ 0.03. However, we note that the value of E_K is sensitive to the upper cutoff of the prior on the E_{K, iso} and is affected by the strong correlation between E_{K, iso} and the other parameters due to the weak constraint on ν_a. In particular, lower values of the kinetic energy are allowed (with the constraint, E_{K, iso, 52} ≳ 30, corresponding to E_K ≳ 3 × 10⁵⁰ and η ∼ 0.4 for our best-fit value of θ_jet = 25). Hence, our estimate of η ∼ 0.03 should be considered a lower bound.

Having performed afterglow modeling of the three existing GRBs at z ≳ 6 with radio through X-ray data to determine the properties of the explosion and environment, we now turn to the question of how these events compare with each other, and with GRBs at lower redshifts. We compile measurements of E_γ, θ_jet, E_K, and n₀ (or A_*) for lower-redshift (z ≲ 1) events from the literature (Panaitescu & Kumar 2002; Yost et al. 2003; Friedman & Bloom 2005; Ghirlanda et al. 2007; Cenko et al. 2010, 2011). Where only a lower limit (or no information) is available for the jet opening angle, we use E_{γ, iso} as an upper bound on E_γ. This combined comparison sample includes GRBs from the pre-Swift era, as well as Swift and Fermi events.

All three z ≳ 6 GRBs presented here are well-fit by a constant density ISM model. In the case of GRBs 090423 and 120521C, the synchrotron self-absorption frequency is not directly observed and hence the best-fit model is only representative of a family of solutions. Despite this uncertainty, we are able to bound ν_a using constraints on the microphysical parameters _e, _B<1/3. We find 1.7 × 10⁸ Hz < ν_a < 2.7 × 10⁹ Hz for GRB 120521C, and 7.9 × 10⁶ Hz < ν_a < 3.2 × 10⁸ Hz for GRB 090423. The corresponding constraints on the physical parameters for these two GRBs are 2.8 × 10⁻⁴ ≲ n₀ cm³ < 27, 2.9 ≲ E_{K, iso, 52} ≲ 29, 3.4 × 10⁻² ≲ _e < 1/3, and 3.5 × 10⁻⁴ ≲ _b < 1/3 for GRB 120521C, and 1.7 × 10⁻⁶ < n₀ cm³ < 8.2, 24 ≲ E_{K, iso, 52} ≲ 5.1 × 10², 1.5 × 10⁻² ≲ _e < 1/3, and 3.3 × 10⁻⁵ < _B < 1/3 for GRB 090423. Together with the high density of n₀ ∼ 600 cm⁻³ for GRB 050904, these three high-redshift GRBs span the lowest to the highest densities inferred from GRB afterglow modeling (Figure 20).

Zoom In Zoom Out Reset image size

The light curves of all three high-redshift events display the signature of a jet break. Using the jet break time, we constrain the opening angle of the jet in each case and find θ_jet ∼ 15–6°. The median¹⁵ jet opening angle of the low-redshift sample is $\theta _{\rm jet}=7.4^{+11}_{-6.6}$ (95% confidence interval; Figure 21). Whereas this interval formally includes the measurements of θ_jet for the high-redshift sample, we note that the observed values of θ_jet for the high-redshift sample are all below the best estimate for the median of the comparison sample, suggesting that higher-redshift events may be more strongly collimated than their lower-redshift counterparts. If this difference is verified with future events, it would indicate that previous studies may have underestimated the beaming correction and therefore the rate of z ≳ 6 GRBs.

Zoom In Zoom Out Reset image size

We use the calculated values of θ_jet to compute the beaming-corrected γ-ray and kinetic energies of the high-redshift GRBs and find that both E_γ and E_K span the range of 3 × 10⁴⁹ erg to ∼10⁵² erg. We confirm previous reports that GRB 050904 is one of the most energetic GRBs ever observed (Gou et al. 2007). GRB 120521C falls in the lower half of the distribution of E_γ and E_K, whereas GRB 090423 lies at the lower end of the distribution of E_γ and near the median of the distribution of E_K. The median values of these parameters for the low-redshift sample are $E_{\gamma } = (8.1^{+11}_{-4.3})\times 10^{50}$ erg and $E_{\rm K}=(3.8^{+17}_{-2.6})\times 10^{50}$ erg. The values of E_γ and E_K for the three high-redshift GRBs span the observed distributions and present no evidence for a substantial difference from the low-redshift sample. The inferred γ-ray efficiencies (η ∼ 0.5) are also similar to the efficiencies of lower-redshift events.

From this comparison, we conclude that the existing sample of z ≳ 6 GRBs displays the same wide range of circumburst densities and beaming-corrected energies as their lower redshift counterparts (Figure 21). On the other hand, the z ≳ 6 events seem to have smaller jet opening angles than the median of the distribution at lower redshifts, suggesting that there might be some evolution in jet collimation with redshift.

We present X-ray, optical/NIR, and radio observations of GRB 120521C and use broadband modeling to deduce a redshift of $z=6.01^{+0.05}_{-0.09}$, consistent with $z\sim 5.93^{+0.11}_{-0.14}$ derived from optical/NIR SED-fitting and z ∼ 6.15 estimated from a low signal-to-noise spectrum. This is only the third GRB at z ≳ 6 for which detailed multi-wavelength observations allow us to extract the properties of the explosion. The data suggest a constant-density circumburst environment with $\log {(n_{\rm 0})}=-2.7^{+1.4}_{-1.0}$, a jet-opening angle of $\theta _{\rm jet}=3.0^{+2.3}_{-1.1}$ deg, beaming-corrected kinetic and γ-ray energies of $E_{\rm K} =(3.1^{+1.9}_{-0.9})\times 10^{50}$ erg and $E_{\gamma }=(2.6^{+4.4}_{-2.0})\times 10^{50}$ erg, and negligible extinction, A_V ≲ 0.05 mag. We also re-fit the other two GRBs at z ≳ 6 with radio detections and compare the properties of the high-redshift sample with those of their lower-redshift counterparts. We find that GRBs at z ≳ 6 exhibit a wide range of explosion energies, circumburst densities, and shock microphysical parameters. The energies and circumburst densities of these high-redshift events are comparable to those of their counterparts at z ∼ 1, and overall, they display no evidence for an evolution in the progenitor population compared to z ∼ 1 events.

We note that GRBs at z ≳ 6 may have systematically smaller jet opening angles, with a mean of θ_jet = 3.6 ± 0.7 deg, which would increase the inferred GRB rate at these redshifts by a factor of ≈4. We caution that our results are based on a small sample of three events at z ≳ 6. The primary reason for the small sample size is the historically low detection rate of GRB afterglows at radio frequencies. Like previous authors, we note that the lack of early-time radio data makes it difficult to determine the synchrotron self-absorption frequency, which in turn results in parameter degeneracies, giving rise to uncertainties in these parameters of several orders of magnitude. Rapid-response radio observations are therefore essential for studying the properties of GRBs, both at low and high redshifts. The recent refurbishment and expansion of the Very Large Array has resulted in an improvement in sensitivity by an order of magnitude, while the Atacama Large Millimeter Array promises to be an excellent facility for the study of GRBs owing to its excellent sensitivity. Detailed studies of high-redshift candidate afterglows with these facilities (e.g., the recent z = 5.913 GRB 130606A; Laskar et al. 2013) will augment this sample and help bring the study of GRBs in the reionization era into the mainstream.

We thank the anonymous referee for detailed comments on the manuscript. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Radio data for GRB 120521C were obtained under VLA project codes 12A-394 and 12A-480. Some of the data presented here were obtained at the Gemini-North Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the NSF on behalf of the Gemini-North partnership. The William Herschel Telescope is operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. We thank Ofer Yaron and S. R. Kulkarni for the Keck observations. The Berger GRB group at Harvard is supported by the National Science Foundation under grant AST-1107973. T.L. acknowledges support by NRAO. K.W. acknowledges support by STFC.

The behavior of the various spectral power law segments of a synchrotron source as outlined in GS02 is strictly valid only when the various spectral break frequencies are located far apart. However, the break frequencies evolve as a function of time and can cross, leading to transitions from one spectral shape to another. Since the normalizations of the light curves in GS02 were calculated in the asymptotic limit, spectral transitions (that occur when break frequencies approach each other and cross) lead to artificial discontinuities in the model light curves. We smooth over these glitches by adding together weighted combinations of all spectra that are accessible with the specified physical parameters (see Section 5 of GS02). For instance, in the ISM model with n₀$E_{\rm K,iso,52}^{4/7}$_B^9/7 < 18, we expect the afterglow to evolve in the order spectrum 5 → 1 → 2. Consequently, in this example, we add together a combination of spectra 5, 1, and 2 with time-varying weights such that the appropriate spectrum presents the dominant contribution in the corresponding asymptotic limit, whereas at a spectral transition (defined next), the two spectra on either side of the transition contribute equally.

For a transition from spectrum A to spectrum B, we define the transition time, t_AB as the geometric mean of the time when spectrum A ceases to be valid and the time when spectrum B first becomes valid. In the above example with the spectra evolving in the order 5 → 1 → 2, there are two transition times, denoted as t₅₁ and t₁₂, respectively.

Next, we construct weighting functions for each spectrum as follows. If a spectrum is valid in the range (− ∞, t_AB] (such as spectrum 5 in the example above), the weighting function (w_L for "left") is unity at early times, and falls as a power law near t_AB, being equal to 1/2 at t_AB:

$$\begin{equation} w_{\rm L}(t,t_{\rm AB}) = \frac{1}{1+(t/t_{\rm AB})^{\eta }}, \end{equation} \tag{ A1 }$$

where η is an ad-hoc parameter that controls the smoothness of the transition.

Similarly, if a spectrum is valid in the range [t_AB, ∞) (such as spectrum 2 in the example above), the weighting function (w_R for "right") rises as a power law at early times, is equal to 1/2 at t_AB, and asymptotes to unity as t → ∞:

$$\begin{equation} w_{\rm R} (t,t_{\rm AB})= \frac{1}{1+\left(t/t_0\right)^{-\eta }}. \end{equation} \tag{ A2 }$$

Finally, for a spectrum that is bracketed by two transition times, [t₁, t₂] (such as spectrum 1 above; note that this can be true of more than one spectrum), we define a weighting function, w_M (for "mid"):

$$\begin{equation} w_{\rm M} (t,t_1,t_2)= w_{\rm R}(t, t_1) + w_{\rm L}(t, t_2) - 1. \end{equation} \tag{ A3 }$$

The compound spectrum at any instant, F_ν(ν, t) is then computed by adding together weighted contributions from all spectra allowed under the given set of physical parameters. For instance, in the above example,

$$\begin{eqnarray} &&F_{\nu }(\nu ,t) =\nonumber\\ &&\quad{\frac{w_{\rm L}(t,t_{51}) F_{\nu }^{(5)}(\nu ,t) + w_{M}(t,t_{51},t_{12}) F_{\nu }^{(1)}(\nu ,t) + w_{\rm R}(t,t_{12}) F_{\nu }^{(2)}(\nu ,t)}{w_{\rm L}(t,t_{51}) + w_{M}(t,t_{51},t_{12}) + w_{\rm R}(t,t_{12})}}\nonumber \\ \end{eqnarray} \tag{ A4 }$$

Since these weighting functions are designed to evaluate to unity far away from a spectral transition and fall as a power law near transitions, the above expression evaluates to the correct spectral shapes in all asymptotic limits. The weighting functions for two adjoining spectra at the transition time are both equal to one half, so both neighboring spectra contribute equally at a spectral transition; this results in smooth light curves at all frequencies even across spectral transitions. Finally, we note that the index η is an arbitrary choice; we find that η = 2 (corresponding to weighting by hyperbolic tangent functions in log-space) works well and yields smooth light curves near transitions, without significantly disturbing the spectrum away from transitions.

Radio emission from a GRB afterglow traversing the Milky Way is susceptible to scintillation—scattering by inhomogeneities in the electron density distribution of the ISM along the line of sight. The phenomenon is often modeled as being produced at a scattering screen located between the source and the observer. The screen produces a speckle pattern on the detection plane, resulting in a modulation of the flux as the observer moves through the speckles. The effect of scintillation decreases above a transition frequency (ν_trans), characteristic of the general direction of the line of sight through the Galaxy (typically around 10 GHz).

The spectrum of the electron density inhomogeneities in the ISM is well-characterized by the Kolmogorov spectrum (Armstrong et al. 1995),

$$\begin{equation} \Phi _{{\rm N}_{\rm e}}(\boldsymbol {q}) = C_{\rm N}^2 q^{-11/3}, \end{equation} \tag{ B1 }$$

where $\boldsymbol {q}$ is the wave-vector and $C_{\rm N}^2$ is a normalization constant that varies from place to place with the Galaxy. The scattering measure is defined as the integral of $C_{\rm N}^2$ from the observer to the scattering screen,

$$\begin{equation} SM = \int _0^{d_{\rm scr}}\!C_{\rm N}^2(x)\,{d}(x). \end{equation} \tag{ B2 }$$

Cordes & Lazio (2002) used pulsar observations to build a model of the electron density distribution in the Galaxy. We use their model, NE2001,¹⁶ to determine the scattering measure and transition frequency along the line of sight to the GRB. We then compute the distance to the scattering screen using the formula (Cordes & Lazio 2002),

$$\begin{equation} d_{\rm scr} = 2\pi \left(\frac{n_{\rm trans}}{318.0\,{\rm GHz}}\right)^{3.4} SM^{-1.2}. \end{equation} \tag{ B3 }$$

The strength of the scattering can be quantified by a parameter, U, defined as

$$\begin{equation} U^{5/3} = \xi = 7.9\times 10^3 SM^{0.6} d_{\rm scr}^{0.5} \left(\frac{\nu }{1\,{\rm GHz}}\right)^{-1.7}, \end{equation} \tag{ B4 }$$

where ν is the observing frequency, with U ≪ 1 and U ≫ 1 corresponding to the weak and strong scattering regimes, respectively (Goodman & Narayan 2006; Walker 1998, 2001).

Having calculated U, we follow the prescription of Goodman & Narayan (2006) to compute the modulation index, m, computing the source size from the formula in Appendix A of Granot & Sari (2002). The expected scatter in the observed flux density due to scintillation is then given by

$$\begin{equation} \Delta F _{\rm scint} = m F_{\rm model}, \end{equation} \tag{ B5 }$$

where F_model is the predicted flux density from the afterglow synchrotron model. We add this uncertainty in quadrature to the flux density uncertainty in each data point prior to performing likelihood analyses.

GRB 120521C exhibits excess radio emission at 21.8 GHz at 1.15 d compared to the best-fit forward shock (FS) model (Figure 8). In Section 5 we suggested that this may be due to contribution from a reverse shock (RS). Here we discuss a self-consistent RS + FS model that accounts for this excess emission. We do not search all possible RS models exhaustively, since the excess emission is observed in only a single data point, but list a plausible model that accounts for the observations.

We begin with a general discussion of the radio light curve of reverse shocks in an ISM environment. In the standard afterglow model, the reverse shock produces a synchrotron spectrum with a characteristic synchrotron frequency (ν_{m, RS}), cooling frequency (ν_{c, RS}), self-absorption frequency (ν_{a, RS}), and overall flux normalization (F_{ν, m, RS}). At the time the reverse shock traverses the ejecta, the deceleration time (t_dec), these parameters are linked to those of the forward shock by the relations, ν_{m, RS}(t_dec) = ν_{m, FS}(t_dec)/Γ², ν_{c, RS}(t_dec) = ν_{c, FS}(t_dec), and $F_{\nu , \rm m, RS}(t_{\rm dec}) = \Gamma F_{\nu , \rm m, FS}(t_{\rm dec})$, where Γ is the Lorentz factor of the ejecta at t_dec. We use the simplest model to explain the data for GRB 120521C and assume that the ejecta are in the slow cooling regime (ν_{c, RS} > ν_{m, RS}) after t_dec, although it is possible that the opposite is true in the initial afterglow phase. At low frequencies and early times the reverse shock emission is expected to be self-absorbed (e.g., Sari & Piran 1999; Berger et al. 2003; Melandri et al. 2010) and the light curve therefore depends upon the relative ordering of ν_{m, RS}, ν_{a, RS}, and the observing frequency.

We note that the 21.8 GHz radio detection for GRB 120521C at 1.15 d, with an excess flux density of 80 μJy compared to the FS model, is preceded by a deeper non-detection at the same frequency at 0.15 d. Subtracting the FS contribution to the 21.8 GHz flux density at 0.15 d, we find an upper limit to the RS contribution at 0.15 d of ≲ 34 μJy. The light curve at 21.8 GHz is thus clearly rising between 0.15 and 1.15 d and falling thereafter, implying that it reached a peak some time between 0.15 and 1.15 d and indicating that the putative RS component is self-absorbed at this frequency at 0.15 d. Regardless of the ordering of ν_{a, RS} and ν_{m, RS}, a peak in the 21.8 GHz light curve must correspond to the passage of ν_{a, RS} through this frequency, since this is the only way to explain a late-time (t > t_dec) turn-over in an RS light curve. If we assume that ν_{m, RS} > ν_{a, RS}, then ν_{m, RS} must pass through 21.8 GHz even later than the apparent peak of the 21.8 GHz light curve at ≈1 d. Our ISM model indicates ν_{m, FS} =5.5 × 10¹¹ Hz at 1 d, implying $\Gamma (t_{\rm dec}) = \sqrt{\nu _{\rm m, FS}/\nu _{\rm m, RS}}\lesssim 5$, which is too low.

We therefore look for a self-consistent RS solution with ν_{m, RS} < ν_{a, RS} at 0.15 d. In this scenario, the light curve rises as t^5/4 prior to the passage of ν_{a, RS}, and then declines as t^{−(3p + 1)/4} ∼ t^−1.88 (using p = 2.17, the median value estimated for the FS). From the upper limit at 0.15 d we can determine the earliest time at which ν_{a, RS} can pass through 21.8 GHz. We find ν_{a, RS} =21.8 GHz at ≳ 0.66 d and F_{ν, a, RS} ≲ 0.2 μJy. This method does not allow us to precisely locate ν_{m, RS}, with the only constraint that it passes through 21.8 GHz at ≲ 0.66 d. If we additionally assume that t_dec ∼ T₉₀ ≈ 27 s, we find a solution that satisfies the relations at the deceleration time with Γ ∼ 70 and ν_{m, RS} ∼2 × 10⁸ Hz at 0.66 d. We show this combined RS + FS model in Figure 22 and note that this model obeys the NIR limits at ≈0.21 d.

Zoom In Zoom Out Reset image size

To summarize, there exists a combined RS + FS model that explains the excess flux density at 21.8 GHz at 1.15 d. Assuming that the deceleration time is of the order of T₉₀, we arrive at an initial Lorentz factor of ∼70 for this GRB, of the correct order of magnitude for GRBs (Sari & Piran 1999; Berger et al. 2003; Soderberg & Ramirez-Ruiz 2003).

The millimeter detections at a flux level of 240 μJy at ≈0.4 and 1.3 d for GRB 090423 are much brighter than expected from the forward shock alone. Based on our best-fit ISM model (Section 6.2 and Figure 17), the expected contribution of the FS to the millimeter flux density is 20 μJy and 35 μJy, respectively, corresponding to an excess flux density of 220 μJy and 205 μJy. We now consider the hypothesis that this excess is due to reverse shock emission and perform an analysis similar to that for GRB 120521C (Appendix C).

As in the case of GRB 120521C, we find that we must have ν_{a, RS} < ν_{m, RS} to avoid a low value of Γ ∼ 5. Given the millimeter data, a similar analysis to that of GRB 120521C indicates that the light curve must have peaked at ≈0.8 d with a flux density of ≈0.5 mJy (we use p = 2.56 derived from the forward shock). The data do not directly constrain ν_{m, RS}. If we assume that t_dec ∼ T₉₀ ≈ 10 s we find a solution that satisfies the relations between the RS and FS at the deceleration time (see Appendix C) with Γ ∼ 500 and ν_{m, RS} =3.5 × 10⁷ Hz at 0.80 d. In this case, the combined RS + FS model (Figure 23) over-predicts the NIR K-band observations around 0.02 d. However, we note that these observations take place at the same time as an X-ray plateau, which could result from energy injection. This would reduce the contribution of the FS to the NIR K-band light curve.

Zoom In Zoom Out Reset image size

In summary, a combined RS + FS model with Γ ∼ 500 can explain the significant excess flux density in the millimeter. The model overpredicts the NIR K-band observations at 0.01 to 0.05 d, which could potentially be explained by a lower contribution from the FS than expected, due to energy injection over this period.