On the electron energy distribution index of Swift Gamma-Ray Burst Afterglows

The electron energy distribution index, p, is a fundamental parameter of the synchrotron emission from a range of astronomical sources. Here we examine one such source of synchrotron emission, Gamma-Ray Burst afterglows observed by the Swift satellite. Within the framework of the blast wave model, we examine the constraints placed on the distribution of p by the observed X-ray spectral indices and parametrise the distribution. We find that the observed distribution of spectral indices are inconsistent with an underlying distribution of p composed of a single discrete value but consistent with a Gaussian distribution centred at p = 2.36 and having a width of 0.59. Furthermore, accepting that the underlying distribution is a Gaussian, we find the majority (>94%) of GRB afterglows in our sample have cooling break frequencies less than the X-ray frequency.


INTRODUCTION
The afterglow emission of Gamma-Ray Bursts (GRBs) is generally described by the blast wave model (Rees & Mészáros 1992;Mészáros et al. 1998) which details the temporal and spectral behaviour of the emission that is created by external shocks when a collimated ultra-relativistic jet ploughs into the circumburst medium, driving a blast wave ahead of it. Fundamental to this model is the electron energy distribution index, p; a characteristic parameter of the process by which the electrons are accelerated to relativistic speeds and by which they radiate via synchrotron emission. This acceleration mechanism, common to many astronomical jet sources (as well as particle acceleration in the solar wind and supernovae, and the acceleration of cosmic rays), is thought to be Fermi diffusive shock acceleration (Fermi 1954) due to the passage of an external shock (Blandford & Ostriker 1978;Rieger et al. 2007) after which the energy of the electrons, E, follows a power-law distribution, dN ∝ E −p dE, with a cut-off at low energies. This is consistent with recent PIC simulations (Spitkovsky 2008) and Monte Carlo models (Achterberg et al. 2001;Ellison & Double 2002;Lemoine & Pelletier 2003) but at odds with others . The blast wave model describes how synchrotron emission from relativistic electrons produces a smoothlybroken, broad band spectrum that is well characterised by a peak flux and three, time evolving, break frequencies (peak frequency, ν m ; cooling frequency, ν c ; synchrotron self-absorption frequency, ν a ) as well as the electron energy distribution index, p (Sari et al. 1998;Granot & Sari 2002). The spectrum is divided into four regimes by the three break frequencies and within each regime the spectrum is asymptotically described by F ν ∝ ν −β , where the spectral index, β, is a function of p only. By comparing the observed X-ray spectra to the predicted asymptotic values of the synchrotron spectra, we can extract information about the electron energy distribution index, p, which is dependent only on the underlying micro-physics of the acceleration process. Some (semi-)analytical calculations and simulations indicate that there is a nearly universal value of ∼ 2.2 − 2.4 (e.g., Kirk et al. 2000;Achterberg et al. 2001;Spitkovsky 2008) though other studies suggest that there is a large range of possible values for p of 1.5 − 4 (Baring 2004).
Observationally, different methods have been applied to samples of BATSE, BeppoSAX and Swift bursts which reached the conclusion that the observed range of p values is not consistent with a single central value of p (Chevalier & Li 2000;Panaitescu & Kumar 2002;Shen et al. 2006;Starling et al. 2008;Curran et al. 2009). The latter three showed that the width of the parent distribution is σ p ∼ 0.3 − 0.5. However, in the studies so far there have been some limitations: multi-wavelength studies (Chevalier & Li 2000;Panaitescu & Kumar 2002;Starling et al. 2008;Curran et al. 2009) suffer from limited samples ( 10 sources each) with sufficient temporal and spectral observations, while those studies that rely on X-ray afterglows alone are subject to a large uncertainty because the position of the cooling frequency, relative to the X-ray, is unknown. The only X-ray study of Swift afterglows so far (Shen et al. 2006) used a very limited sample of spectral indices (∼ 30), dictated by the number of GRBs observed by Swift at the time, to estimate the distribution of p. Neither, they did not take a statistical approach to the position of the cooling frequency relative to the X-ray regime, as we do here. We interpret a much larger (∼ 300) and, statistically, more significant sample of Swift observed GRB afterglows, to constrain the distribution of the values of electron energy distribution index, p. In §2 we introduce our method while in §3 we present the results of our Monte Carlo analyses and their implications in the overall context of GRB observations and particle acceleration in general. We summarise our findings in §4. All uncertainties are quoted at the 1σ confidence level.

METHOD
Our general method is to constrain the electron energy distribution index, p, from the values of the Xray spectral indices, β X or β, observed by the Swift XRT (Burrows et al. 2005) and detailed by Evans et al. (2009, Table 7). A normalised histogram of the spectral indices of these 301 GRB spectra is plotted in figure 1. We derive p from the spectral index as opposed to the temporal index because for a given spectral index, assuming the asymptotic limit, there are only two possible values of p depending on whether the cooling frequency, ν c is less than or greater than the X-ray regime, ν X , while for a given temporal index there are multiple possible values which are model dependent (e.g., the simple blast wave model (Rees & Mészáros 1992;Mészáros et al. 1998), or modifications thereof (e.g., Granot & Kumar 2006;Genet et al. 2007;Ghisellini et al. 2007;Uhm & Beloborodov 2007); for a discussion on the choice of X-ray spectral index to derive p see Curran et al. 2009). In accordance with synchrotron emission predicted by the blast wave model, we ascribe the behaviour of the unabsorbed X-ray spectrum to be a single power law where the flux goes as: F ν (ν) ∝ ν −β and β is the spectral index. Under the standard assumptions of slow cooling and adiabatic expansion of the blast wave, the electron energy distribution index is related, in the asymptotic limit, to the spectral index by either p = 2β (ν c < ν X ) or p = 2β + 1 (ν c > ν X ) (e.g., Granot & Sari 2002); implying a difference between the slopes of the two spectral regimes of ∆β = 0.5. Throughout we use the regime probability, X, as the probability that the cooling frequency is less than the frequency of the X-ray regime (i.e., ν c < ν X ) and 1−X is the probability that ν c > ν X . We neglect the case where the cooling frequency may be passing through the X-ray regime (since there is no sign of spectral evolution in the sample) and the cases where the peak frequency, ν m , or self absorption frequency, ν a , is greater than the X-ray regime as this is not observed in late time afterglows.
To parameterise the underlying distribution of the electron energy distribution index, p, from the X-ray spectral indices observed by Swift we use a maximum likelihood Monte Carlo method. This method uses a maximum likelihood fit to return the most likely parameters of the assumed underlying model, the errors on which are estimated via a Monte Carlo error analysis. Another Monte Carlo analysis tests the probability that the observed distribution of spectral indices could be obtained from an underlying distribution of p described by the most likely parameters.
In this method we first assume a model for the underlying distribution of p which we transform into a distribution in spectral index space, via the regime probability, X. We convolve this distribution with the measured probability of the data set and calculate the loglikelihood of the parameters of the underlying model (see appendix A). To estimate the most likely parameters of the underlying model and the regime probability, X, we minimise the log-likelihood (equates to maximising the -Normalised histogram of the data (301 measurements) overlaid with high-resolution normalised histograms of the synthesized data sets (10 4 × 301 data points) from the most likely parameters of a single discrete p (grey line) and a Gaussian distribution of p (blue line) as detailed in Table 1. likelihood) using the simulated annealing method ( § 10.9 of Press et al. 1992, and references therein). Uncertainties of the fit parameters are estimated via a Monte Carlo analysis, whereby the observed data are randomly perturbed within their (asymmetric) Gaussian errors and refit multiple times (10 4 ); the standard deviation of returned most likely parameters is used as a measure of the uncertainties. We then find the probability that the observed distribution of spectral indices could have originated from an underlying distribution of p based on the most likely parameters estimated from the log-likelihood minimisation. We do so by generating 10 4 synthetic data sets drawn randomly from the underlying probability distribution of p, randomly transformed into spectral index space via the regime probability, X, and further randomly perturbed within the (asymmetric) Gaussian errors of each of our original data points. Each of these synthetic data sets is fit as above and the value of the loglikelihoods recorded; one would expect that if the data are consistent with the underlying model, at the 1σ level, the original log-likelihood value should fall within the 1σ distribution of the synthetic values. With 10 4 synthetic data sets we can measure the percentage to an accuracy of two decimal places and rule out chance agreements to the 4σ level; the ∼ 2 × 10 5 synthetic data sets required to rule out chance agreements at the 5σ level was considered too costly, computationally.
There are two underlying models, or hypotheses, regarding the data that we want to test: that the observed distribution of spectral indices, β, can be obtained from an underlying distribution of p composed of i) a single discrete value (SDp) and ii) a Gaussian distribution (GDp). Details of these models and their log-likelihoods are discussed in the Appendix, A.

RESULTS AND DISCUSSION
3.1. Distribution of p The results of our analysis are detailed in Table 1 which shows our most likely parameters for electron energy distribution index, p, the related Gaussian scatter, σ p , the probability that the cooling frequency is less than the -The most likely values of electron energy distribution index, p, the related Gaussian scatter, σp, the probability that the cooling frequency is less than the X-ray frequency (νc < ν X ), X, and the log-likelihood of that fit, l. Values in brackets are the average and error values from the Monte Carlo error analysis of perturbed data sets, while those in square brackets are the average and standard deviations returned from the Monte Carlo analysis of synthesized data sets. % is the percentage of synthesized data sets with a better fit than the original. X-ray frequency (ν c < ν X ), X, and the log-likelihood of that fit, l. Values in brackets are the average and error values from the Monte Carlo error analysis of perturbed data sets, while those in square brackets are the average and standard deviations returned from the Monte Carlo analysis of synthesized data sets. % is the percentage of synthesized data sets with a better fit than the original.
Though our most likely discrete value of p = 2.25 agrees well with the predicted universal value of p ∼ 2.2 − 2.3 (e.g., Kirk et al. 2000;Achterberg et al. 2001), it is comprehensively rejected by our tests; the hypothesis that the observed distribution of spectral indices, β, can be obtained from an underlying distribution of p composed of a single discrete value (SDp) is rejected at the 4σ level as all synthesized data sets had better log-likelihoods. The hypothesis that the observed distribution can be obtained from an underlying Gaussian distribution (GDp) centred at p = 2.36, having a width of 0.590 and regime probability, X = 1.000, is acceptable at the 1.5σ level (Figure 2). As a visual aid, we compare the normalised histogram of the observed data with the high-resolution normalised, average histograms of the 10 4 synthesized data sets (Figure 1). Note that the SDp model is clearly a poor fit and exhibits a secondary peak at β ∼ 0.6 due to the fact that X = 0.96 for the likelihood fit of that model to the observed spectral indices. This result confirms the results from previous small-sample GRB afterglow studies (Shen et al. 2006;Starling et al. 2008;Curran et al. 2009) as regards the non-universality of p, the central value at p ∼ 2.0 − 2.5, and the width of the distribution of σ p ∼ 0.3 − 0.5. However, our results are based on a sample of bursts an order of magnitude larger than these studies and we took a statistical approach to the position of cooling frequency relative to the X-ray, by using the regime probability, X.
Given that the value of p is not universal, it may be possible that it changes suddenly or evolves gradually with time or radius even in a single event as environmental shock parameters (e.g., magnetic field, ambient density) change or evolve (e.g., Hamilton & Petrosian 1992;Vlahos et al. 2004;Kaiser 2005). It is also possible that different components of a structured jet (Mészáros et al. 1998;Kumar & Granot 2003), multicomponent jet (Berger et al. 2003;Huang et al. 2004) or jet-cocoon (Zhang et al. 2003)   rameter, though a change or evolution of p should be observable as a change or evolution of the synchrotron spectral index, β, and no significant example of such an evolution has been observed in GRB afterglows.
3.2. Limits on regime probability, X Though previous multi-band (optical -X-ray) studies (e.g., Panaitescu & Kumar 2002;Starling et al. 2008;Curran et al. 2009) have shown that the cooling break frequencies of a number of GRB afterglows are greater than the X-ray frequency (ν c > ν X ), here, accepting that the underlying distribution of p is a Gaussian, we find that this number is consistent with zero (X = 1). From our Monte Carlo error analysis (where we refit the randomly perturbed data multiple times) we can plot the distribution of possible values of X (Figure 3; note that this is a log-log plot). The distribution of the errors is clearly not Gaussian so we should not use the standard deviation as a measure of error as we have in Table 1. We can however place a nominal 3σ lower limit on the value of X by estimating the value at which there is 99.7 percent containment; we find that this is at 0.935 < X < 0.936, compared to the nominal 1σ (68.2 percent containment) at 0.988 < X < 0.989. Any value of X less than this limit is inconsistent with the distribution of the data and would produce a clear secondary peak at β = (p − 1)/2, not observed in the data (Figure 1). To avoid this peak, an extremely wide distribution is needed, much wider than the spread of the observed data.
Hence, the upper limit on the percentage of GRBs in our sample where the cooling frequency is greater than the X-ray frequency (ν c > ν X ) is approximately 6.5 percent, or ∼ 20 out of 301 GRBs. This is a discrepancy from the ratio observed in the multi-band studies (e.g., Starling et al.; Curran et al.: 5 out of 10 and 2 out of 6, respectively) but not seriously so given the extremely low number statistics of those studies. If we create a sample of 6 random bursts from our statistical distribution of 301, there would be a non-negligible (∼ 5 percent) chance that 2, or more, bursts would have a cooling frequency greater than the X-ray frequency, consistent with the aforementioned Swift study (Curran et al.). The study of Starling et al. is based on a sample of BeppoSAX, as opposed to Swift, GRBs which may have a different limit on the regime probability, X. An investigation as to why the cooling frequency of afterglows follow this apparent limit, or its implications regarding the distribution of other parameters, is beyond the scope of this work.

CONCLUSION
We use the X-ray spectral indices of gamma-ray burst afterglows observed by the XRT aboard Swift to parameterise the underlying distribution of the electron energy distribution index, p, within the framework of the blast wave model. The electron energy distribution index is a fundamental parameter of the synchrotron emission from a range of astronomical sources and in this case of the synchrotron emission of GRB afterglows. We use a maximum likelihood Monte Carlo analysis to test two hypotheses, namely that the observed distribution of spectral indices, β, can be obtained from an underlying distribution of p composed of i) a single discrete value and ii) a Gaussian distribution. We find that the observed distribution of spectral indices are inconsistent with the first hypothesis but consistent with the second, a Gaussian distribution centred at p = 2.36 and having a width of 0.59. Furthermore, if we accept that the underlying distribution is a Gaussian, the majority ( 94 percent) of GRB afterglows in our sample have a cooling break frequency less than the X-ray frequency.
We thank the referee for their comments. PAC, PAE, MJP, MdP acknowledge support from STFC. AJvdH was supported by an appointment to the NASA Postdoctoral Program at the MSFC, administered by Oak Ridge Associated Universities through a contract with NASA.
where L(p, σ p , X) is the likelihood function. This is a numerically calculable function that is minimised to find the most likely parameters, errors on which can be estimated via a Monte Carlo error analysis. If the distribution of the electron energy distribution index, p, can be described by a single discrete value, P(p|p) = 1 (p =p) 0 (p =p) , then the probability of β is P(β|p) = X (β =p/2) 1 − X (β =p/2 − 0.5) 0 (β =p/2 & β =p/2 − 0.5) (A6) and the above convolved probabilities and likelihoods hold with σ p = 0.