ESTIMATING REDSHIFTS FOR LONG GAMMA-RAY BURSTS

The redshifts have long been an important issue for long-duration gamma-ray bursts (GRBs). Since the first X-ray, optical, and radio counterparts were discovered in 1997 (Costa et al. 1997; van Paradijs et al. 1997; Frail et al. 1997), GRBs are confirmed to be in galaxies at cosmological distances. However, until now, only a small fraction of the bursts have their redshifts measured. Even Swift, over its first four years of operation, has only roughly 30% of its bursts with spectroscopic/photometric redshifts. These redshifts may have complex selection biases as a function of redshift relating to the difficulty of getting redshifts for faint and distant bursts as well as the distribution of bands used for spectra and photometry. One illustration of the importance of the selection effects is that the average of redshift value for Swift GRBs have been dropped from 〈z〉 = 2.8 (Jakobsson et al. 2006c) to 〈z〉 = 2.1 (Jakobsson 2008d) while the average redshift value of earlier satellites is much lower yet (Bagoly et al. 2008). As the real redshift distribution of GRBs cannot be changed by that much for the past few years, the only reason for it would be changing selection effects in the spectroscopic observations.

These systematic biases will greatly affect various analyses. As long-duration GRBs come from the collapse of very massive fast-rotating stars with very short main-sequence lifetimes (Woosley & Bloom 2006), their rate density will provide a measurement of the massive star formation rate of our universe. If we only deal with bursts with known spectroscopic redshifts, the derived star formation rate may be biased by the spectroscopic-redshift selection effects. The same problem comes with demographic studies of the GRB luminosity function. If only bursts with spectroscopic redshifts are used in constructing the luminosity functions, then the selection biases may distort the derived luminosity function. A third important problem is the recognition of the highest redshift bursts (say, with z > 7), against which the current spectroscopic methods may be heavily biased. Swift is expected to have 5%–10% of its bursts with z > 7 (Bromm & Loeb 2006), but the spectroscopic methods have only identified one such burst GRB090423, and that just a few month ago (Tanvir et al. 2009; Salvaterra et al. 2009).

A solution to these problems is to use the luminosity relations that work on luminosity indicators measured directly from the prompt gamma-ray emission. These luminosity relations are equations that connect a measurable property of the burst itself (the luminosity indicators) to the burst's energy (E_γ or E_γ,iso) or the peak luminosity (L). The observed fluence or peak flux can then be used (with the inverse square law) to derive a distance to the burst and (for some fiducial cosmology) the burst redshift. The advantage of this procedure is that it applies to all long class GRBs, not just to the minority selected to have a measured spectroscopic redshift. The disadvantage is that the uncertainties on the derived redshifts are much larger than those of the spectroscopic redshifts, and so the method can be used primarily for statistical or demographic purposes. The situation is similar to the case where photometric redshifts of distant galaxies in the Hubble Deep Field do not have the accuracy of the spectroscopic redshifts, yet nevertheless these photometric redshifts can be done for all the galaxies and are the cornerstone of all statistical studies of the field.

The six luminosity relations we use in this article are: τ_lag (spectral lag)–L relation (Norris et al. 2000), V (Variability)–L relation (Fenimore & Ramirez-Ruiz 2000), E_peak (peak energy of the spectrum)–L relation (Schaefer 2003), E_peak–E_γ relation (also called Ghirlanda relation; Ghirlanda et al. 2004), τ_RT (minimum rise time)–L relation (Schaefer 2002), and N_peak (number of peaks in the light curve)–L relation (Schaefer 2002). It was claimed by Butler et al. (2007) that the E_peak–E_γ,iso relation (also called Amati relation; Amati et al. 2002) differs somewhat between Swift and pre-Swift data, and hence the relation might be caused by the detection threshold effect of the instrument instead of the GRB itself. We tested this claim on the relations we are using, as shown in Section 3.

In this article, we will demonstrate a method to use these luminosity relations to estimate the GRBs redshifts. Here, we will only consider bursts with known spectroscopic redshifts, but our derivation of z_best is based on the burst properties, and the effect of the known spectroscopic redshifts are negligible. This effect is tested later in the article. By applying our method to the bursts with accurate spectroscopic redshifts measured, we can compare our estimated redshifts with the measured spectroscopic redshifts. This will be the key test of our methods and the accuracy of our derived redshifts, as well as the applicability of our methods to bursts without spectroscopic redshift. With confidence in the reliability of the method, we can then apply it to all Swift bursts including those with no spectroscopic redshifts. Also, we will be able to apply this method to future bursts, and provide the community with a rapid notification about their redshifts (Xiao & Schaefer 2009).

In this paper, we adopted five luminosity indicators with six luminosity relations (with two relations for indicator E_peak: E_peak–L relation and E_peak–E_γ,iso relation). The details of the calculations for each of these indicators are described in Schaefer (2007a). In this paper, we constructed the calculation independently, and made some significant improvements in the calculation of τ_lag, τ_RT, and N_peak. The five indicators are listed as below.

• 1.  
  The spectral lag, τ_lag, is the delay time between the soft and hard light curves of a burst. By convention, we use the soft and hard energy bands to be those of BATSE channels 1 and 3, or Swift channels 2 and 4, covering roughly 25–50 keV and 100–350 keV. For τ_lag, by shifting the hard and soft light curves of a GRB, and calculating the cross-correlation between them, we are able to get a cross-correlation versus offset plot. The offset corresponding with the peak value of the cross-correlation is the lag time we need. The calculation is simple and easy for bright bursts, while for those faint ones, since the plot has significant scatters, the offset with the peak correlation (i.e., τ_lag) is hard to evaluate under these noisy conditions. To find the offset when the cross-correlation achieves its peak value, we need to make a reasonable fit to the peak region of the cross-correlation. As the shape and scatter of the plot varies from burst to burst, we cannot simply fit it with some specific function. If we fit it with a parabola, then any asymmetry in the cross-correlation (as is often seen) will incorrectly shift the peak in the model by an amount depending on the range of offset included in the fit. And if we fit it with a high-order polynomial function, then the high-order terms will be unstable as they are trying to follow noise and regions far from the peak. What we did was to fit the cross-correlation by polynomials with different orders (normally from 3 to 9), and choose the one which fits best in the very central region (around the peak) of the curve. Two examples are shown in Figure 1.A bootstrap procedure was previously used to calculate the uncertainties on the τ_lag (Norris 2002). In our work, we are using a simple propagation method. The uncertainties on the cross-correlation amplitude points are calculated by simply evaluating the rms scatters of these data points around the fitted curves. We are then able to generate the uncertainties for each of the fitting parameters, and the coefficient errors. Then the uncertainty for our τ_lag value, $\sigma _{\tau _{\rm lag}}$, is generated by propagation of the uncertainties on the fitting parameters and the coefficients.
• 2.  
  Variability (V) measures whether a light curve is spiky or smooth, and can be obtained by calculating the normalized variance of the original light curve around the smoothed light curve. The calculation of V is as shown in Schaefer (2007a):

  $$\begin{equation} V = \big\langle [(C-C_{\rm smooth})^2-\sigma _{C}^2]/C_{\rm smooth,max}^2\big\rangle, \end{equation} \tag{ 1 }$$

  where C is the count per time bin in the background subtracted light curve, with an uncertainty of σ_C. The time duration of the burst T_up is calculated as the summation of the time with the light curve brighter than 10% of its peak flux, and C_smooth is the count in the smoothed light curve, with a box-smoothing width to be 30% of T_up. C_smooth,max is the peak value of C_smooth.From Equation (1), the 1σ uncertainty of V can be propagated from the observational uncertainty of C in each time bin. To get rid of the cross-correlation effect between C and C_smooth, for each C_i, we calculate it with

  $$\begin{equation} x_i = C_i-C_{i,{\rm smooth}} = C_i - \frac{1}{N}(\sum _{j=i-N/2}^{i-1}C_j+\sum _{j=i+1}^{i+N/2}C_j), \end{equation} \tag{ 2 }$$

  where N is the box smooth bin (i.e., 30% of T_up).If we neglect the uncertainty of σ_C, the 1σ uncertainty of V can be propagated from Equation (1) and the uncertainties of each x_i and C_smooth,max.The variability values (V) in Table 1 are about 10 times larger than those in Schaefer (2007a). In this article, we are strictly following the definition and calculation of Equation (1). By checking the details of the calculation, we found that in Schaefer (2007a), C_max was used in the denominator instead of C_smooth,max, which caused the difference. The values of variability are very sensitive to the slight change on its calculation function, and this might be one of the reasons for the large scatter of the variability–luminosity relation.
• 3.  
  E_peak is the photon energy at which the νF_ν spectrum is the brightest. By fitting the GRB spectrum with a smoothly broken power law (Band et al. 1993) as

  $$\begin{eqnarray} \hspace{-10pt}\Phi (E){=}\left\lbrace \begin{array}{@{}l@{}l} A E^{\alpha } e^{-(2+\alpha)E/E_{\rm peak}}, & \mbox{\,if $E \le [(\alpha - \beta)/(2+\alpha)]E_{\rm peak}$} \\ B E^{\beta }, & \mbox{otherwise} \end{array}\!\!\!, \right. \nonumber \\ \end{eqnarray} \tag{ 3 }$$

  the peak energy E_peak and corresponding parameters α (the asymptotic power-law index for photon energies below the break) and β (the power-law index for photon energies above the break) can be obtained. Here, Φ is the usual differential photon spectrum (dN/dE) as a function of the photon energy (E). A GRB spectrum is also able to be fit with a power law with exponential cutoff model:

  $$\begin{equation} {\it dN}/{\it dE} = {\it A E}^{\alpha } e^{-(2+\alpha)E/E_{\rm peak}}, \end{equation} \tag{ 4 }$$

  with E_peak and power-law index α being the fitting parameters. Here, A and B are normalizing constants to indicate the brightness and constructed to ensure the continuity of the model spectrum.In this article, we are using E_peak values of the bursts from various of sources. All the values we used and the sources are listed in Table 2.
• 4.  
  The minimum rise time in the light curve, τ_RT, was proposed for use in a luminosity relation by Schaefer (2002). The minimum rise time of a burst is taken to be the shortest time over which the light curve rises by half the peak flux of the pulse. In practice, especially for faint bursts with large Poisson noise, the rate difference between two close bins might be larger than half of the peak flux. As a result, we have to smooth the light curve before we calculate the rise time. The problem is then that if we smooth it too little, the apparent fastest rise time might be dominated by the Poisson noise, resulting in a too small rise time, and if we smooth it too much, the smoothing effect will dominate, resulting in a rise time near the smoothing time bin.As the light curves vary greatly among different bursts, there is no specific box-smoothing width that can satisfy a majority of bursts. Instead, we vary the box-smoothing width from 0 bins (that is, no binning) up to a relatively large number, say 50 bins, and for each of these smoothing widths, we generate a smoothed light curve from which a minimum rise time can be calculated. Of course, some of these light curves are over-smoothed and some are under-smoothed. In a minimum rise time versus smoothing width plot, we will have a monotonically rising curve, as shown in Figure 2. Although the shape of the curve varies among bursts, for most of the bursts there will be a region where the curve appears flat, or with a slightly increasing slope, which we call a plateau. It is easy to explain the existence of the plateau: in this region, the smoothing is enough that Poisson noise is negligible in determining the minimum rise times, while the smoothing is not so much that it determines the minimum rise time. On a plot of minimum rise time versus smoothing width, we can identify three regions: a fast rising region where the Poisson noise dominates, a nearly flat plateau region where we are seeing the real minimum rise time, following by another rising region for large box-smoothing widths where the smoothing is dominating. We can make use of this plateau region by extrapolating it back to the zero-smoothing case (where the box-smoothing width is zero), as shown in Figure 2. The intercept on the y-axis corresponds with the minimum rise time where the smoothing is effectively zero. With the extrapolation, we are side-stepping the regime where the Poisson noise dominates. As a result, the value of the intercept is just the minimum rise time we need, not affected by either the Poisson noise effect or the smoothing effect.For some of the extremely faint bursts, the Poisson noise dominant region and the smoothing effect dominant region will overlap with each other, and we are unable to find a plateau in the minimum rise time versus smooth width plot, for which the extrapolation cannot be made. Thus, our technique does not produce τ_RT values for the faintest bursts.The uncertainty of the minimum rise time is calculated by simple propagation from the fitting parameters and the uncertainties on each individual point on the minimum rise time versus smoothing width plot. The uncertainties on each of the individual points are dominated by the noise on each original data points in the light curve. The noise on the peak flux will affect our criteria, by some factor of σ_C/max(C), where C is the rate, and σ_C is the uncertainty of the rate. In addition to that, random noise on the start and stop data points for each possible rise time will affect our determination, i.e., the real rise time between two data points may be larger/smaller than half of the peak flux, however, with the random noise on the start and stop points, we take it as our rise time, which is equal to half of the peak flux. This effect is also reflected on our criteria also, by a factor of ∼2σ_C/max(C). As a result, to determine the uncertainties on each of the individual rise time on RT-smoothing width plot, we can change our criteria from 0.5 by a factor of 2.25σ_C/max(C), and record how much the resulted minimum rise time values changes. The uncertainties on the fitting parameters and the minimum rise time value can then be calculated from propagation.
• 5.  
  N_peak is defined as the number of peaks in the light curve. With C_max as the overall maximum of the background-subtracted light curve, we define a peak to be a local maximum that rises higher than C_max/4 and is also separated from all other peaks by a local minimum that is at least C_max/4 below the lower peak. In principle, N_peak is easy to count, either automatically or manually. In practice, we have the same problem of the Poisson noise and the smoothing factor effect, as what we had in the calculation of τ_RT. Faint bursts will have their unsmoothed light curves dominated by apparent peaks produced by Poisson noise, resulting in large numbers of false peaks. A random noise spike can satisfy our definition for a peak if we do not smooth the light curve, yet if we smooth it too much there will always be just one peak. Here, we adopted the same procedure of calculation as that in the calculation of τ_RT. We vary the box smooth width from 0 bins to a relatively large number, and calculate the number of peaks for each of the smoothed light curves. As in the case for τ_RT, we see a fast falling curve (where Poisson noise is contributing spurious peaks), with a plateau, where neither the Poisson noise effect nor the smoothing effect dominates. By extrapolating the plateau back to the y-axis, we get an N_peak value for an unsmoothed case of the light curve, with the effects of Poisson noise removed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Each of the luminosity indicators discussed above has one or more corresponding luminosity relations. These relations are τ_lag–L, V–L, E_peak–L, E_peak–E_γ (so called Amati relation) and E_peak–E_γ,iso (so called Ghirlanda relation), τ_RT–L, and N_peak–L relation, which are described and explained in details in Schaefer (2007a). In our redshift calculation program, we are not including Amati relation for the following two reasons. First, the Amati relation has been challenged as it returns ambiguous redshifts (Li 2007), while the E_peak–L relation (also the τ_lag–L, V–L, τ_RT–L, and N_peak–L relations) passed the same test (Schaefer & Collazzi 2007). Second, the physics of the Amati relation is nearly the same as that of Ghirlanda's E_peak–E_γ relation, except that Ghirlanda's E_peak–E_γ relation has a correction for the jet opening angle, making it much tighter. There is another luminosity relation that relates E_peak, an effective duration (called T₄₅), and the luminosity by Firmani et al. (2006). However, with more GRBs and E_peak data added in, it is realized that this luminosity relation is not making any improvement on the E_peak–L relation (Collazzi & Schaefer 2008). Some X-ray indicators and relations are also proposed as estimators of GRB redshift (Gendre et al. 2008), which is out of the scope of this paper.

We have 107 long GRBs with their spectroscopic/photometric redshifts measured, ranging from 1997 February 28 (GRB970228) to 2008 July 21 (GRB080721), observed by BATSE, Konus, HETE, and Swift. The average redshift for pre-Swift bursts is about 1.50 and that of the Swift bursts is 2.15.

The light curves of Swift GRBs are generated from the original data published on the legacy ftp site,¹ and the Swift Software ver 2.9 (HEAsoft 6.5). To generate a background-subtracted light curve of a GRB, we downloaded an event file sw00xxxxxx000.bevshsp_uf.evt.gz and a mask file sw00xxxxxx000bcbdq.hk.gz, with xxxxxx be the six digit Swift trigger number. By running a task "batbinevt," we can specify the time interval, energy bins, time bin method, output file name, and format on generating the light curve. In our work, we adopted a time interval of 0.064 s with uniform time bins, and four continuous energy bands (15–25 keV, 25–50 keV, 50–100 keV, and 100–350 keV). For the calculation of V, τ_RT, and N_peak, we have been using the light curve over the whole energy range (15–350 keV), and for the calculation of τ_lag value, we use 25–50 keV and 100–350 keV. The light curves of pre-Swift bursts are obtained from our previous work (Schaefer 2007a).

We calculated the τ_lag, V, τ_RT, and N_peak values for each of these bursts, as listed in Table 1. The first column lists the ID number of GRBs. The second column lists the satellite with the detection of the burst. Columns 3–6 show all the calculated indicator values, with the name of the indicators shown on the header row. Indicators not measured due to low signal-to-noise ratio (S/N) of the burst are represented as "..."

The values of E_peak as well as the power-law indexes α, β in Band's smoothly broken power-law model (Band et al. 1993) or E_peak and α values in a power law with exponential cutoff model are obtained from various sources, as shown in Table 2. The first column in Table 2 lists the ID number of the GRBs. The second column lists the satellite with the detection of the burst. Columns 3–6 are the values of E_peak and the power-law index α and β values, as well as the reference sources of these values. Our jet break time (t_jet) values from optical observations and their sources are also listed in Table 2, Columns 7 and 8. Bursts without measured jet break time are filled by "..." also. From the table, we see that only 33 of the bursts have their optical t_jet reported. Various jet break time in the X-ray afterglows have been reported, however, as there are usually multiple breaks for the X-ray afterglows, which are not well understood and distinguished for their causes, we are not including any of these reported t_jet values from X-ray detections. Values in square brackets for α and β are assumed values for those bursts without exact α and β values measured, which is taken to be the average value of known α and β (Schaefer et al. 1994; Krimm et al. 2009; Kaneko et al. 2006; Band et al. 1993). Some uncertainties of E_peak and t_jet are also quoted in square brackets. These uncertainties are assumed conservative assumed values, which are normally 10% of the measured E_peak or t_jet values. All of the peak flux (P) and fluence (S) values of pre-Swift bursts we use here are similar to what were used in Schaefer (2007a), and those of Swift bursts are from the data table on the Swift Web site.² All the quoted error bars in Tables 1 and 2 are converted to 1σ level.

The luminosity relations are all expressed as power laws, and we can make a linear fit on the logarithms of the redshift-corrected luminosity indicators and the logarithms of the burst luminosities. In Figure 3, we display the data and the best fits for each of these luminosity relations, for the concordance cosmological model (w = −1 and Ω_M = 0.27 in a flat universe). As there are significant uncertainties in both the horizontal and vertical axes, and some intrinsic scatters independent of both the luminosity and the redshift, we performed the ordinary least squares without any weighting (Isobe et al. 1990; Schaefer 2007a). As both the indicator and luminosity are caused by some other parameter (like the jet Lorentz factor in particular), so both indicator and luminosity are correlated. In this case, we assume that these two variables in the luminosity relations are not directly causative, and the bisector of the two ordinary least squares (Isobe et al. 1990) has been used. More details of the fitting process are referred to Schaefer (2007). As E_peak plays an important role in the calculation of P_bolo, the uncertainties on luminosity is then correlated with the uncertainties on E_peak. In our work, we ignored that effect, and simply assumed that P_bolo and E_peak have uncorrelated errors. The best-fitting functions for each of the relations are shown in Table 3.

Zoom In Zoom Out Reset image size

The V − L relation is rather scattered, as we can see from Figure 3 and Table 3. The scatter is larger than one could expected for any linear relations. As the calculation of V is not well defined, and V − L relation is our most noisy luminosity relations. In this case, we decided not to include V − L relation in our later calculation for the redshifts, although we listed all the calculated V values for the bursts in our tables. As a result, the maximum number of luminosity indicators that we can use is four, and the maximum number of luminosity relations is five.

Butler et al. (2007) estimated the Amati's relation (E_peak–E_γ,iso) with Swift and pre-Swift data. They report an inconsistency in the relations from the two data sets, which they then attributed to differences in threshold between Swift and earlier detectors. Their claimed difference has not been reproduced by other groups (including Cabrera et al. 2007; Schaefer 2007b; Krimm et al. 2009), while their claimed threshold effects have been found to not significantly affect the observed Amati relation (Schaefer 2007b; Nava et al. 2008; Ghirlanda et al. 2008; but see Shahmoradi & Nemiroff 2009). Indeed, their analysis is based on Bayesian priors which systematically push high-E_peak values below ∼400 keV, as demonstrated by detailed comparisons with Konus, Suzaku, and RHESSI measures. Nevertheless, in this paper, we can perform yet another test to see whether the claimed threshold differences between Swift and pre-Swift bursts cause any significant change in the luminosity relations. Butler et al. (2007) includes Swift BAT bursts between GRBs 041220 and 070509, 77 of which have spectroscopic redshift measured. While in our analysis on Swift GRBs, we are including all bursts with spectroscopic redshift between 050126 and 080721. With both samples being based on largely overlapping samples selected in a nearly identical manner, we conclude that the flux limits of the two samples are essentially identical.

To this end, we have separately fitted the pre-Swift and Swift data. This test was only done for four luminosity relations (τ_lag–L, V–L, E_peak–L, and τ_RT–L), with the E_peak–E_γ relation having too few bursts, and the N_peak–L relation not being usable for the comparison as a limit. The best-fit luminosity relations are given in Table 4. The data and best-fit models are displayed in Figure 4. At first glance, we see that the difference between the two best-fit lines are small compared to the scatter in the data, and a detailed analysis is described below.

Zoom In Zoom Out Reset image size

We made a F-test for the fitting results for these four relations. First, we made a bisector linear fit on the combined data with both pre-Swift and Swift bursts, and recorded the χ² value of the fit as χ²_JOINT. Then we separate the data to two sample sets, pre-Swift and Swift, and made the same bisector linear fit separately on each set of the data. The sum of the two χ² values for the separately fitted lines are recorded as χ²_SEPARATE. Then the F value can be calculated as

$$\begin{equation} F = \frac{\chi _{\rm JOINT}^2/(N_{\rm pre}+N_{sw}-2)}{\chi _{\rm SEPARATE}^2/(N_{\rm pre}+N_{sw}-4)}, \end{equation} \tag{ 5 }$$

where N_pre is the number of pre-Swift bursts, and N_sw is the number of Swift bursts. N_pre + N_sw − 2 is the degree of freedom of the fitting on the combined data, and N_pre + N_sw − 4 is the degrees of freedom of the fitting on separated data.

These F values for each of the luminosity relations are listed in Table 4. If the pre-Swift and Swift relations differ much from each other, the separate fitting have been significantly improved over the fitting on the mixed data, the F value would be much larger than unity. Otherwise, if there is no significant difference between pre-Swift and Swift relations, the F value would be around unity. From Table 4, we see that F is rather close to unity. All of the results show that the separately fitted result is not significantly improved over the fitted results on all data mixed together, which means that luminosity relations for pre-Swift luminosity relations and Swift luminosity relations do not significantly differ from each other. And by looking at the plots in Figure 4, we see that (1) the envelope of squares and diamonds are indistinguishable, and (2) the pre-Swift and Swift best-fit lines are close to each other compared to the scatter in the data. And hence, we have no significant evidence that these four luminosity relations differ for Swift bursts.

The 1σ range on the normalization difference between Swift and pre-Swift bursts for all these four luminosity relations are also listed in Table 4. We can make an analysis with the normalization difference on Amati's relation claimed by Butler et al. (2007), which is corresponding with a 0.39 difference in log space. By making the comparison between the Butler's factor (0.39) and our normalization difference, we can exclude Butler's factor at a 2.5σ level for τ_lag–L relation, a 3.7σ level for τ_RT–L relation, a 1.3σ level for E_peak–L relation (which is not significant), and we cannot exclude the Butler's factor for V–L relation.

Our aim in this paper is to test our method of redshift calculation. We applied it to the bursts with known spectroscopic redshifts (z_spec), and if it works well, we will be able to apply it to all the long GRBs in our future work. Although we are dealing with the bursts with known z_spec, these z_spec are only involved in the fitting of luminosity relations, and this effect is negligible in our calculation. Only after we have derived our redshift z_best from the luminosity relations will we compare them with the known z_spec to test our accuracy.

Below we will describe how our method applies on one GRB step by step.

• 1.  
  First, we measure each of the luminosity indicators of the burst. The definition and method of calculations have been discussed in Sections 2 and 3. The results of the indicators for all bursts in the sample are listed in Tables 1 and 2.
• 2.  
  Next, we derive the luminosity values for each relations from Table 3. A complexity is that the luminosity relations depend on the redshift of the burst (so as to correct the luminosity indicators back to the burst rest frame), so we have to perform this calculation for an array of trial redshifts (we take it to be from redshifts of 0 to 20 at intervals of 0.005), and then we will obtain a list of luminosities (or isotropic energies for Ghirlanda's relation) values depending on the list of trial redshifts for each of the indicators. We notate each of these calculated luminosities (as a function of redshift z_trial) for the ith relation as L_i(z_trial) (or E_γ,i for the Ghirlanda's relation).
• 3.  
  With the values of the peak flux P, fluence S, E_peak, and the power-law indexes in the broken power-law model α and β (or α from the power law with exponential cutoff model), the bolometric peak flux P_bolo and fluence S_bolo can be calculated. The range for "bolometric" is set to be 1 keV–10,000 keV in the GRB rest frame, and the equations for the detailed calculation are referred to Schaefer (2007a). As a result, for each burst, we calculated P_bolo(z_trial) and S_bolo(z_trial) for each trial redshift value from 0 to 20.For those bursts with t_jet values, the jet opening angle θ_jet (in units of degrees) can be calculated as

  $$\begin{equation} \theta _{\rm jet} = 0.161\left[\frac{t_{\rm jet}}{1+z_{\rm trial}}\right]^{3/8}\left[\frac{n\eta _{\gamma }}{E_{\gamma,{\rm iso},52}}\right]^{1/8}, \end{equation} \tag{ 6 }$$

  where t_jet is the jet break time in the unit of days, n is the density of the circumburst medium in particles per cubic centimeter, η_γ is the radiative efficiency, and E_γ,iso,52 is the isotropic energy in units of 10⁵² erg (Rhoads 1997; Sari et al. 1999). We simply adopt η_γ = 0.2 and n = 3 cm⁻³ in Equation (6). The beaming factor, F_beam, is then calculated as

  $$\begin{equation} F_{\rm beam} = 1- \cos \theta _{\rm jet}. \end{equation} \tag{ 7 }$$

  From above we see that, since both θ_jet and E_γ,iso are redshift sensitive, our calculated F_beam value also varies between different z_trial values.
• 4.  
  From all the parameters above, for each of the indicators, a list of the luminosity distances can be calculated as

  $$\begin{equation} d_{L,i}(z_{\rm trial}) = \sqrt{\frac{L_{i}(z_{\rm trial})}{4\pi P_{\rm bolo}(z_{\rm trial})}}. \end{equation} \tag{ 8 }$$

  For Ghirlanda's relation, the list of luminosity distances is calculated as

  $$\begin{equation} d_{L,i}(z_{\rm trial}) = \sqrt{\frac{E_{\gamma,i}(z_{\rm trial})[1+z_{\rm trial}]}{4\pi F_{\rm beam} S_{\rm bolo}(z_{\rm trial})}}. \end{equation} \tag{ 9 }$$

  Then from each list of luminosity distance above, the distance modulus can be obtained:

  $$\begin{equation} \mu _{i}(z_{\rm trial}) = 5 \log [d_{L,i}(z_{\rm trial})] -5, \end{equation} \tag{ 10 }$$

  with d_L,i(z_trial) expressed in units of parsecs. The uncertainties are propagated strictly following the calculation.From all above, we are able to get up to four lists of measured distance moduli μ_i(z_trial) with their 1σ uncertainties: $\mu _{\tau _{\rm lag}} \pm \sigma _{\mu _{\tau _{\rm lag}}}$, $\mu _{E_{\rm peak}}\pm \sigma _{\mu _{E_{\rm peak}}}$, $\mu _{E_{\rm peak}-E_{\gamma }}\pm \sigma _{\mu _{E_{\rm peak}-E_{\gamma }}}$, and $\mu _{\tau _{RT}}\, \pm \, \sigma _{\mu _{\tau _{RT}}}$. As we have asymmetric uncertainties for E_peak, we carry the uncertainties on both directions in the calculation, and generated both plus and minus uncertainties for $\mu _{E_{\rm peak}}$ and $\mu _{E_{\rm peak}-E_{\gamma }}$. $\mu _{N_{\rm peak}}$ is also calculated, as a lower limit on the distance modulus. All of these distance moduli are a function of the assumed z_trial for 0 < z_trial < 20.
• 5.  
  Given each trial redshift, we can calculate its distance modulus directly from the cosmological model, μ_cos(z_trial). Here we adopt the concordance model, with equation of state for dark energy p = wρc², w = −1, Ω_M = 0.27, Ω_Λ = 1 − Ω_M = 0.73, and H₀ = 70 (km s⁻¹) Mpc⁻¹. In this case, the luminosity distance can be expressed as

  $$\begin{equation} d_L(z_{\rm trial}) = cH_0^{-1} (1+z_{\rm trial}) \int _{0}^{z_{\rm trial}} dz^{\prime } [(1+z^{\prime })^3 \Omega _M + \Omega _{\Lambda }]^{-1/2}. \end{equation} \tag{ 11 }$$

  From the equation above and our list of trial redshifts, a list of luminosity distances d_L(z_trial) will be calculated, and also a list of distance modulus μ_cos(z_trial) which equals 5log [d_L(z_trial)] − 5. The μ_cos(z_trial) values are only depending on the trial redshift (running from 0 to 20) and the cosmological model we choose.
• 6.  
  For each of the trial redshifts from 0 to 20, we have distance moduli lists of $\mu _{\tau _{\rm lag}}$, $\mu _{E_{\rm peak}}$, $\mu _{E_{\rm peak}-E_{\gamma }}$, $\mu _{\tau _{RT}}$ along with their 1σ uncertainties as well as μ_cos. We can then compare these μ_i(z_trial) with μ_cos(z_trial) in a χ² sense. Thus, $\chi ^2_{\tau _{\rm lag}} = [(\mu _{\tau _{\rm lag}}-\mu _{\rm cos})/\sigma _{\mu _{\tau _{\rm lag}}}]^2$ and so on for the other relations. We get a χ²_i versus z_trial plot, as shown in the left panel of Figure 5. Then sum over all the χ²_i, we get a χ²_total ($\chi _{\rm total}^2 = \chi _{\tau _{\rm lag}} ^2+ \chi _{E_{\rm peak}}^2 + \chi _{E_{\rm peak}-E_{\gamma }}^2 + \chi _{\tau _{RT}}^2 + \chi ^2_{N_{\rm peak}}$) versus z_trial plot, as in the right panel of Figure 5. Our best redshift (z_best) corresponds with the minimum χ²_total, where the luminosity relations and the cosmological model agree with each other best. We are also able to find the uncertainties of our z_best. By searching through the χ²_total–z_trial plot, we can find the redshifts with which the χ²_total = χ²_total,min + 1, which corresponds with the edges of the 1σ range of our redshift. Similarly, the redshifts with χ²_total = χ²_total,min + 4 gives us the 2σ range and that with χ²_total = χ²_total,min + 9 gives the 3σ range of our redshift.

Zoom In Zoom Out Reset image size

The resulting values of z_best and 1σ range of z for each burst are shown in Table 5. We have also collected the spectroscopic redshifts for these bursts (see values and references in Table 5). The last column of Table 5 lists the effective luminosity relations we used for each burst in our calculation. The number of luminosity indicators used in the redshift calculation for each of the burst is listed in Table 5. A comparison can be made between our calculated redshifts and their spectroscopic (or photometric) redshifts. The comparison plot is shown is Figure 6, and below are some of the conclusions after we made the analysis.

Zoom In Zoom Out Reset image size

• 1.  
  Of the total 115 bursts, eight have only a lower z limit that can be calculated. For the remaining 107 bursts, we took the spectroscopic redshifts z_spec as the model value, and our calculated redshift z_best as a measured value, with the 1σ uncertainty of $\sigma _{z_{\rm Plus}}$ and $\sigma _{z_{\rm Minus}}$. Then the χ² can be calculated as

  $$\begin{equation} \chi ^2 = \sum _{j=1}^{N_{\rm burst}}{(z_{j, {\rm best}} - z_{j, {\rm spec}})^2/\sigma _{j}^2}, \end{equation} \tag{ 12 }$$

  with σ_j equals to $\sigma _{j,z_{\rm Plus}}$ or $\sigma _{j,z_{\rm Minus}}$, depending on whether our z_i,best is smaller or larger than the z_i,spec. Each j represents an index number identifying the burst. The number of degrees of freedom in this comparison equals the number of bursts (N_burst), so the reduced χ² is χ²/N_burst. The reduced χ² value is 1.28, which is somewhat larger than unity. While after excluding one 3σ outliers GRB010222, whose contribution to χ² is as high as 44, our reduced χ² is equal to 0.86. This is certainly not a significant deviation from unity. So we conclude that the scatter in Figure 6 is consistent with our quoted error bars being correct.
• 2.  
  Of the 107 bursts with their z_best calculated, 73 have their z_spec falling into the 1σ range of our z. The ratio of the numbers is about 70%, which is slightly larger than the ideal case 68.5%. And for the eight bursts with only lower z limit calculated, six of them have their z_spec larger than our 1σ lower limit of z. This is another way of testing our quoted error bars, and as in the previous item, we find no significant deviation from the expected results. As such, to a close degree, we see that our derived error bars are accurate.
• 3.  
  We can test to see whether our z_best is biased high or low. For this, we calculated the average value of log₁₀(z_best/z_spec). If our result is unbiased, the average should be zero to within the error bars. We find the average is 0.01. This demonstrates that our z_best is not biased to within the 1% level.
• 4.  
  To test the accuracy of our z_best comparing with the z_spec, we calculated the rms scatter of log₁₀(z_best/z_spec). The result comes out to be 0.26, while the rms scatter of log₁₀(z_spec) is 0.30. In Schaefer (2007a), it was pointed out that the accuracy of the redshift estimation is 26% (corresponding to a log₁₀(z_best/z_spec) rms of 0.11), which is better than what we are claiming here. The reason for the larger rms scatter is, although we are dealing with the GRBs with known z_spec, we are not making any use of the z_spec in our whole calculation, and our luminosity L, luminosity distance d_L, distance modulus μ_ind, as well as the bolometric flux P_bolo and fluence S_bolo are all varying with our trial redshift. This brought us one extra degree of freedom in the calculation, which caused larger uncertainties in our result. When the known z_spec value is used (as for the Hubble diagram work in Schaefer 2007a), the scatter becomes substantially small as compared to our work in this paper. Another reason for the large rms is because it is mostly dominated by the very noise bursts (those bursts with low S/N, inaccurate measurement of luminosity indicators, and large error bars for z_best results), which are not able to be used for the notification of the redshifts. If we make a selection on bursts with relatively accurate redshift estimation, say bursts with $\sigma _{z_{\rm minus}/z_{\rm best}}<0.5$ and $\sigma _{z_{\rm plus}/z_{\rm best}}<1$, we get the rms of log(z_best/z_spec) of 0.19, which is much smaller, and the rms scatter of log₁₀(z_spec) for this subsample is 0.28.
• 5.  
  We made the same calculation with pre-Swift luminosity relations on calculating Swift redshifts and Swift luminosity relations on calculating pre-Swift redshifts, for which we call it z_best,II. Our calculation on z_best,II is totally independent of the spectroscopic redshifts, as for each burst, the luminosity relations used in the calculation are calibrated independent of the z_spec for any individual burst. The comparison between z_best and z_best,II will also show the difference between Swift and pre-Swift luminosity relations. The comparison plot between z_spec and z_best,II is shown in Figure 7. From the comparison between Figures 6 and 7 we see that the scatter and distribution of z_best and z_best,II do not differ significantly from each other. Actually, from our calculation, the average value of log₁₀(z_best,II/z_spec) is −0.02, and the rms scatter of log₁₀(z_best,II/z_spec) is 0.27, both of which are equal to those of our z_best value within error bars. This result demonstrate that those two sets of z_best values do not differ from each other, which means that these two sets of luminosity relations do not have significant difference. It also tells us that the effect of redshift involved in our calculation (in fitting luminosity relations) is negligible.
• 6.  
  We need to verify whether our result is effective in selecting high redshift bursts. If our predicted redshift is z, the possibility of a real redshift to be higher than z and lower than z are both 50%, which cannot be used as a test. However, by considering our uncertainties of z_best, if our predicted redshift is 2z, we can make a test by counting how many of the GRBs are with real redshifts larger than z. As there are not many bursts with high redshifts, a test is done on a relatively lower redshift region, where most of the GRBs are involved. We picked up all our GRBs with predicted z > 4, the total number is 12, and 0 out of 12 have their spectroscopic redshift z < 2. This result tells us that our method is actually effective in demographic studies and in picking up high z bursts, if we take into the consideration of the error bars before we make the prediction.

Zoom In Zoom Out Reset image size

From all the analysis above, we can conclude that our z_best is not biased on average, and our 1σ error bars are accurate. We can claim that our method works well on the bursts with known spectroscopic redshifts, and can be applied to all long GRBs (even without their z_spec measured).

In this paper, we developed a method to calculate the redshifts for long GRBs, using their light curves and spectra. We applied our method to bursts with known spectroscopic redshifts, detected by BATSE, HETE, Konus, and Swift. By comparing our calculated redshifts with their spectroscopic redshifts, we are able to examine the accuracy of our method.

We compared each of the luminosity relations for pre-Swift and Swift bursts by making a F-test. With the F values close to unity, we have significant evidence against any claim that the relations are caused by the detection threshold effects or any other artificial effects of the instruments.

We compared our results with the spectroscopic redshifts. We find that our z_best are not biased (with the average value of log₁₀(z_best/z_spec) equal to 0.01), and our reported 1σ error bars are good (with χ²_red = 0.86 and 70% of the z_spec fall into the 1σ region of z_best). Our accuracy on the redshifts are not as accurate as those from spectroscopy, yet nevertheless with a reasonable accuracy for demographical and statistical studies, with the rms of log₁₀(z_best/z_spec) is 0.26. The rms value is about twice as what was found in Schaefer (2007a). One of the reason is, in Schaefer (2007a), the accuracy is calculated assuming a known z_spec, and in this paper, as we are treating the unknown-redshift case, so extra degrees of freedom has been brought in the calculation, which caused the accuracy to get worse by about a factor of 2. Another reason is that the large rms is dominated by those faint and noisy burst, and for a subsample with $\sigma _{z_{\rm minus}/z_{\rm best}}<0.5$ and $\sigma _{z_{\rm plus}/z_{\rm best}}<1$, we get the rms of [log(z_best/z_spec)] of 0.19, which is much smaller. As our z_best are from the light curves, the spectra and the concordance cosmological model, it is independent of the spectroscopic redshift. As a result, our method can be applied to all long GRBs.

For Swift bursts, as we are measuring over a relatively narrow energy band (15 keV–350 keV), the uncertainties in the calculation of peak energy E_peak are large, and it becomes one of the main restrictions on our accuracy. With the launch of Fermi, we are hoping to get bursts with more accurately measured E_peak values and light curves covering a broader energy band. We are expecting a substantial improvement in the accuracy of redshifts for Fermi bursts.

For the next step, we will apply our method to all Swift long GRBs, aiming to get a nearly complete Swift GRBs redshift catalog. Such a catalog will inevitably be incomplete due to bursts with incomplete light curves and bursts too faint for their properties to be usefully measured. Our resulting redshifts will have an accuracy worse than those obtainable with optical spectroscopy, yet our accuracy will be good enough for various important statistical studies. As such, our catalog will be used for the demographic studies, without the detection threshold effect of the spectroscopic-redshift measurements. We are also trying to deal with all the future bursts, and to provide a rapid notification of the redshift on the GCN circular to the community. We are hoping to find some possible high redshift (e.g., z > 7) and possible low redshift (e.g., z < 0.3) GRBs, which will be important in the observation of the high redshift universe and the GRB–SNe connection study.

This work is supported by NASA under grant NESSF 07Astro07F-0029.