Evolutions and Calibrations of Long GammaRay Bursts Luminosity Correlations Revisited
Abstract
Luminosity correlations of long Gammaray bursts (GRB) are extensively proposed as an effective complementarity to trace the Hubble diagram of Universe at high redshifts, which is of great importance to explore properties of dark energy. Recently, several empirical luminosity correlations have been statistically proposed from GRB observations. However, to treat GRB as the distance indicator, there are two key issues: the redshift evolution of luminosity correlations and their calibrations. In this paper, we choose the Amati relation, the correlation between the peak spectra energy and the equivalent isotropic energy of GRBs (), as an example, and find that the current GRB dataset implies that there could be a evolution of the luminosity correlation with respect to the redshift. Therefore, we propose an extended Amati relation with two extra redshiftdependent terms to correct the redshift evolution of GRB relation. Secondly, we carefully check the reliability of the calibration method using the lowredshift GRB data. Importantly, we find that the lowredshift calibration method does not take whole correlations between and coefficients into account. Neglecting these correlation information can break the degeneracies and obtain the biased constraint on which is very sensitive to values of parameters for the calibration. A small shift of parameters of “calibrated” relation could significantly change the final constraint on in the lowredshift calibration method. Finally, we simulate several GRB samples with different statistical errors and find that, in order to correctly recover the fiducial value of using the lowredshift calibration method, we need a large number of GRB samples with high precisions.
Subject headings:
cosmological parameters gammaray burst: general1. Introduction
Since the accelerating expansion of the universe was revealed by observations of type Ia supernovae (SN) [39, 37], SN plays an essential role to study the expansion history of the Universe and the nature of dark energy. However, due to the limited intrinsic luminosity and the extinction from the interstellar medium, the maximum redshift of the SN we can currently observe is about . Approximately, for most other popular probes at low redshift, such as baryon acoustic oscillations (BAO) and weak lensing, the bottleneck of redshift range in the near future observations is about . Whereas fluctuations of the cosmic microwave background (CMB) provide cosmological information at the last scattering surface (). Therefore, we do not have many effective methods to observe the evolution of our Universe at , usually referred as the “cosmological desert”.
Gammaray bursts (GRB) are the most energetic explosions after the big bang in the Universe. The equivalent isotropic energy radiated in a few seconds can reach up to erg (see e.g., [38, 34, 35, 25] for recent reviews). Thanks to the extreme brightness and the immunity to dust extinction of highenergy photons in the gammaray band, GRB are detectable up to redshift [40, 44, 7, 45]. Therefore, in the literatures, GRB are widely proposed, as a complementary data of SN, to trace the Hubble diagram of the Universe at high redshifts. However, different from the consistent luminosities of SN which allows us to qualify SN as an ideal distance indicator, the central engine mechanism of explosions of GRB has not yet clearly been understood. Therefore, drawing cosmological implications from GRB observations is quite intractable. Recently, many methods have been proposed to achieve some progresses, in order to treat GRB as the distance indicator [41, 9, 16, 29, 14, 42, 30, 27, 52, 51, 33, 28, 32, 46].
Luminosity or energy correlations of GRB are empirical connections between measurable properties of the prompt gammaray emission and the luminosity or energy of GRB. In recent years, several empirical luminosity correlations have been statistically proposed from observations, such as various popular twovariable relations: the correlation between spectrum lag and isotropic peak luminosity () [36]; the correlation between time variability and isotropic peak luminosity () [12]; the correlation between the peak energy of spectrum and isotropic equivalent energy (, the Amati relation) [3]; the correlation between peak energy and collimationcorrected energy () [17]; the correlation between peak energy and isotropic peak luminosity () [41, 53]; the correlation between minimum rise time of light curve and isotropic peak luminosity () [42]; and some wellknown multivariable relations, such as the correlation between , , and the break time of the optical afterglow light curves [29]. The general expression for luminosity correlations is
(1) 
where (or , , ), (or , ). Determinations of coefficients and of correlations could tell us whether these empirical luminosity correlations are independent on the redshift, which is very crucial for treating GRB as the distance indicator to investigate the evolution of our Universe. In addition, some new models for GRB have been proposed in resent years, such as the spectrotemporal multicomponent model [19, 20, 21, 22, 23, 24]. Interestingly, this model could result in a new timeresolved relation: the relation between the luminosity of the nonthermal component, , and its corresponding spectral peak energy in the rest frame, (i.e., relation).
More importantly, determinations of coefficients can also show the reliability of the calibration of GRB luminosity correlations directly. In the literatures, there are two methods to calibrate the luminosity correlations: “selfcalibration” method, which is based on the global fitting technique using all GRB data to constrain the cosmological parameters and the coefficients of luminosity correlations simultaneously [41, 27]; “lowredshift calibration” method, which is using the lowredshift GRB data to obtain the constraints on coefficients of luminosity correlations, and directly extrapolate them to high redshifts to constrain the cosmological parameters using the highredshift GRB data [30, 52, 51, 11, 33, 32, 50].
Apparently, these two methods are seriously dependent on the assumption that empirical luminosity correlations are universal and do not evolve with respect to redshift. In practice, for the Amati relation, Li [26] first revealed that the systematically significant variations of the intercept and the slope with the cosmological redshift. Later, Basilakos & Perivolaropoulos [5] investigated the same issue for several other luminosity relations. Due to the limited quality of GRB samples, there was no statistically significant evidence for evolution of the calibration parameters. In addition, as a cosmological probe, the evolution of GRB luminosity correlations may result in overestimate or underestimate of cosmological parameters [10]. More recently, in Lin et al. [31], they reinvestigated the Amati relation using lowredshift () and highredshift () GRB, and found that the coefficients of Amati relation in lowredshift GRB differs from those of highredshift GRB at more than confidence level. Therefore, they finally concluded that long gammaray bursts might not be standard candles. Following on this topic, Lin et al. [32] further examined the possible redshift dependence of several other luminosity correlations and found that only for the relation, the lowredshift GRB could give the similar constraints on the coefficients with those from the highredshift GRB within confidence level.
In this paper, we choose the Amati relation as an example to investigate the evolutions and calibrations of GRB luminosity correlation with the latest GRB observational data [33]. Firstly, we divide the whole GRB sample into two bins ( and ). In each bin, we use the GRB data to constrain the coefficients in the CDM framework to check whether they are independent on the redshift. Then, we propose an extended Amati relation by introducing two extra terms to characterize the evolutions with respect to the redshift. Secondly, we carefully check the reliability of the calibration method. In practice, we use the lowredshift GRB data () to calibrate the coefficients of GRB luminosity correlation and directly extrapolate them to highredshift events. Then we use this “calibrated” Amati relation to constrain the matter density parameter with the highredshift GRB data. Finally, we simulate several samples of GRB data with high precisions to investigate that in which conditions the GRB data calibrated by the lowredshift calibration method can be used for cosmological studies.
2. Evolutions of Luminosity Correlations
2.1. Amati Relation & GRB Data
Although for recent years, many empirical luminosity correlations have been statistically concluded from long GRB observations, the Amati correlation is the most widely used one among them for cosmological studies. The original version of the Amati relation is expressed as
(2) 
where to correct the redshift dilation of the spectrum [48]. The isotropic equivalent energy can be calculated from the bolometric fluences [42]:
(3) 
The uncertainty of propagates from the uncertainties of and . is calculated from the observed peak photon flux in the rest frame keV energy band by assuming the Band spectrum [4]. To calculate the luminosity distance , we assume a flat concordance CDM model with the Hubble constant: and the matter density parameter: obtained from the latest 2015 results [1]. Hence, in Eq. (3) can be written as
(4) 
relation  bins  N  

low  
Amati relation:  high  
Global Fit  
Extended  low  
Amati relation:  high  
Global Fit 
To test the possible redshift dependence of the Amati luminosity correlation, we use the GRB sample compiled in Liu & Wei [33], which includes 138 wellmeasured GRB in the redshift range: . In general, we examine the evolution of the coefficients and according the following steps: Firstly, we get the equivalent isotropic energy for each GRB at redshift with the distance in Eq. (4) taken into account. Then, we divide all observed 138 GRB sample into two redshift bins: [0.0331, 1.4], [1.4, 8.1]. The reason of choosing as the threshold is that the maximal redshift of SNe Ia data samples we are using for cosmological studies is usually around . The Universe below has already been well studied by using the SNe Ia datasets [2, 43, 6]. Therefore, we treat the GRB data with and as the lowz and highz samples, respectively. In these two redshift bins, there are 59 and 79 GRB data, respectively. And then, we constrain the luminosity correlation coefficients and by maximizing the D’Agostini’s likelihood [8]:
(5) 
or, equivalently, by minimizing the :
(6) 
where is the intrinsic scatter which represents any other unknown uncertainties except for the observational statistical ones. Here, in order to test the feasibility of D’Agostini’s likelihood, we constrain the two parameters of Amati relation use all GRB data (138 GRB) in the flat CDM framework firstly. The minimum chisquared value is and . We also plot the comparison between the GRB data and the theoretical prediction of the bestfit values in Fig. 1. All these results show that the D’Agostiniâs likelihood is feasible to constrain parameters. Finally, after getting constraints on and , in Fig.2 we plot the 2 dimensional marginalized constraints in the plane of (, ) from the lowredshift (blue contours) and highredshift (red contours) GRB data, respectively.
In our analysis, we use emcee^{1}^{1}1https://pypi.python.org/pypi/emcee introduced by ForemanMackey et al. [15], a Python module that uses the Markov chain Monte Carlo (MCMC) method to get the bestfit values and their uncertainties of parameters , and by generating sample points of the probability distribution. The bestfit values with 1 errors are shown in Table 1. We also plot the results in Fig. 2 (the upper right panel), in which the blue and red contours denote the results on the parameters of Amati relation from the lowredshift and highredshift GRB subsamples, respectively. Apparently, these two constraints are not match very well. The best fit values are more than away. When considering the uncertainties, the tension between two fitting results is less than confidence level, which implies that the parameters of Amati relation might be redshiftdependent. The result is similar with some other works [31].
2.2. Redshift Evolution of Amati Relation
As the statistically obvious evolution of the luminosity correlation has been revealed, we should be careful to use the GRB data as the distance indicator for cosmological studies. Because all these investigations are based on the assumption that the luminosity correlation does not vary with respect to the redshift, i.e., both and in Eq. (2) should be constant. To solve this problem, inspired by the study for the evolution of lightcurve fitting parameters in SN observations [49, 47] and that for the evolution of the GRB intrinsic luminosity [54], we introduce two extra redshiftdependent terms to characterize the redshift evolution of the luminosity correlation. Since most redshifts of GRB are much greater than unity, rather than the linear expression in Wang & Wang [49], Wang, Li & Zhang [47] and the logarithmic parametrization in Yu et al. [54], we select two mild formulas
(7) 
to avoid extreme results at high redshifts. Therefore, the extended Amati relation can be expressed as
(8) 
where ,, and are new four coefficients in the relation which should be obtained from the GRB data.
Following the same procedures for analyzing the evolutions of and presented in the subsection 2.1, we also investigate the evolutions of luminosity correlation coefficients , , and using the GRB sample in the same lowredshift and highredshift subsamples, respectively. The bestfit values with 1 errors are shown in Table 1. Fig. 2 (the lower left six panels) shows the results regarding evolutions of luminosity correlation coefficients (, , and ). Since we have two parameters to describe the redshift dependence, compared with the result of original Amati relation, the tension between constraints of parameters from low and high GRB data are alleviated, most of which are consistent within confidence level. Of course, these two extra free parameters could bring the large uncertainties and make the constraints weaker.
3. Calibrations of Luminosity Correlations
3.1. Low Redshift Calibration
As shown in section 2.2, the tension between the constraints on parameters of the original Amati relation from the lowredshift and highredshift GRB subsamples is more than confidence level, while the constraints are consistent within confidence level for the extended Amati relation. However, this does not mean that we can safely use the GRB data for cosmological studies. We need to check the calibration of GRB luminosity correlations further. Usually, we have two calibration methods. One is to calibrate luminosity correlations with lowredshift GRB samples in the context of the flat CDM model and then extrapolate these obtained coefficients to high redshift range. The other method is the global fitting analysis using the MCMC technique in which luminosity correlation parameters and cosmological parameters are simultaneously fitted on the same weight in the context of the flat CDM model [27]. Therefore, here we use these two methods to calibrate the extended Amati relation and compare the constraints on the cosmological parameters.
In the lowredshift calibration method, we obtain the luminosity distance of the lowredshift GRB sample () by using Eq. (4). And then we calibrate the luminosity correlation by maximizing the D’Agostini’s likelihood (Eq. 2.1) for the 59 lowredshift GRBs. The constraints on the luminosity correlation coefficients are shown in Tab. 1. Then we directly extrapolate these coefficients to high redshifts. Therefore, we construct the Hubble diagram for the 79 highredshift GRBs based on the calibrated luminosity correlation and investigate cosmological implication from this Hubble diagram in the context of standard CDM model. In this framework, the luminosity distance is determined by Eq. (4), and the corresponding distance modulus is
(9) 
We use the numerical method to fit and the corresponding is
(10) 
where is the total uncertainty of distance modulus for the th GRB and it is propagated from the uncertainties of
(11) 
where
Firstly, we consider the original Amati relation. In the upper right panel of Fig. 3 (red dashed line), we show the onedimensional distribution of the matter density parameter obtained from the 79 highredshift GRB data by using the lowredshift calibration method. If we use the bestfit values and to calibrate the highredshift GRB data, we obtain the 68% C.L. constraint on the matter density parameter:
(12) 
Based on this result, we can see that this obtained constraint on is far away from the fiducial value of , which is used to constrain the parameter and from the low GRB data for calibration. The significance of this inconsistence is about confidence level.
Then, we use the extended Amati relation to do the calculations, following the same procedures. In the lower right panel of Fig. 3 (red dashed line), we also show the onedimensional distribution of obtained from the 79 highredshift GRB data by using the lowredshift calibration method. The obtained constraints on is:
(13) 
at 68% confidence level, which is consistent with previous works [50]. This result is consistent with the fiducial value of the matter density parameter very well, which also implies that the extended Amati relation might be better than the original one.
3.2. Global Fitting Method
The second method is the global fitting method, in which we calculate the matter density parameter and luminosity correlation coefficients (, , and ) simultaneously with the same weight. The constraints on these coefficients from allredshift GRB data are also listed in Table 1.
Similar with the discussions in the section 2, we start with the original Amati relation. In Table 1 and Fig. 4, we show that constraints on the parameters of Amati relation from the low GRB data in the lowredshift calibration method and from all GRB data in the global fitting method are more or less consistent within confidence level. As we know, the constraint on the matter density parameter is strongly correlated with these coefficients of GRB relation. Therefore, the constraint on using the global fitting method is quite different from that using the lowredshift calibration method:
(14) 
at confidence level, as shown in the blue line of the upper right panel of Fig. 3. This large difference implies that the obtained constraint on using the lowredshift calibration method could be biased and overestimated. On the other hand, the constraint on using the global fitting method is far away from the current bestfit values from other observational data , which might imply that the original Amati relation could not be the best relation to describe the GRB data.
Next, we move to the extended Amati relation. In Table 1 we also show the constraints on parameters of the extended Amati relation from the low GRB data in the lowredshift calibration method and from all GRB data in the global fitting method, which are completely consistent with each other at about confidence level. This also implies that the extended Amati relation could ease the tension of constraints on parameters of GRB relation obtained from different GRB data combinations, similar with the conclusion in the previous section 2. Then we compare the constraints on the matter density parameter using these two methods. In Fig. 3 we plot the twodimensional marginalized distribution with 1 and 2 contours between and , as well as the onedimensional distributions of and (blue lines). We also show the constraint on using the lowredshift calibration method (red dashed line) in the lowerright panel for comparison.
Clearly, we can see that the constraint on using the global fitting method is nearly unconstrained, while the powerful constraint on is obtained by using the lowredshift calibration method. After our careful numerical checks, we find that the reason of this huge difference is that the global fitting method includes all the correlations among parameters, such as the matter density parameter and parameters of GRB relation (as shown in the lowerleft panel of Fig. 3), during the calculations, while this lowredshift calibration method totally neglects these strong degeneracies by force and only uses the simple error propagation equation [see Eq.(11)], which only includes the standard deviations of parameters themselves, to estimate the final statistical errors for the data. This will inevitably underestimate the statistical errors of data and then overestimate the constraining power on the matter density parameter. This means that the lowredshift calibration method gives the strongly biased and overestimated constraint on the matter density parameter, which is untrustable.
In the two dimensional contour (lower left panel) of Fig. 3, we add a vertical line which denote the bestfit value obtained from the lowredshift GRB data alone in the lowredshift calibration method. This line crosses the and contours with the crossing values 0.11 0.36 and 0.05 0.76. We find that these values are quite similar with the and , lower and upper limits of the constraint on the matter density parameter by using the lowredshift calibration method (Eq. 13). This analysis implies that in the lowredshift calibration method, the strong correlations between and coefficients of the extended Amati correlation are not fully taken into account in the calculation. Neglecting this information could break the strong degeneracies between and correlations coefficients by force and then give quite stringent constraint on the matter density parameter.
Furthermore, in the lowredshift calibration method, we find that the final constraint on is very sensitive on how to calibrate the relation. Here, we assume three sets of parameters to calibrate the extended Amati relation:

(a) the bestfit values obtained from the low GRB data;

(b) the bestfit values obtained from all GRB data;

(c) , , and , which is quite similar with the first case (a), except the .
Based on the Table 1, we can easily know these three cases are consistent with each other at confidence level. Now we use the lowredshift calibration method to constrain , which is shown in the lowerright panel of Fig. 3 [(a) the red dashed line, (b) green dash dot line and (c) black dotted line, respectively]. The case (a) is the standard case in the lowredshift calibration method and obtains the stringent constraint on (see Eq. 13). But if we use case (b) to calibrate the relation, the obtained constraint on the matter density parameter is quite different from that in case (a):
(15) 
The small shift of the parameters of “calibrated” relation could significantly change the final constraint on in the lowredshift calibration method, which implies that the constraint on is very sensitive to the values of parameters for the calibration, when we do not fully take the strong correlations among parameters into account and underestimate the statistical errors. When we use the case (c), the situation becomes worse. The obtained constraint on is weaker and only has the lower limit: at 95% confidence level.
We also perform the similar check in the original Amati relation. The constraints on the matter density parameter using the cases (a) and (b) to calibrate the relation are shown in the upperright panel of Fig.3 [(a) the red dashed line and (b) green dash dot line, respectively]. We obtain the same conclusion that the constraint on is very sensitive to the values of parameters for the calibration in the lowredshift calibration method.
3.3. Simulation Results
For the time being, the number of GRB data we can use for possible cosmological studies are very small, when comparing with various type Ia supernovae data sets. There are still very large intrinsic scatter and large statistical errors on data in the current GRB data set, as we discuss above. In order to check how large biased constraint on the matter density parameter by using the lowredshift calibration method, we simulate the mock GRB data with high precisions.
To get the simulation data of GRBs, we first plot the distribution of the redshift and the observed quantity peak energy of the 138 GRBs sample [33] in Fig. 5 with the blue histogram. Then we find two distribution functions to trace the histogram of redshift and . For the redshift , we assumes the distribution function has the form
(16) 
We also assumes the distribution of peak energy is gamma distribution
(17) 
where and are parameters, and the gamma function is
(18) 
We put the distribution function of in Fig. 5 (right panel) the form of
(19) 
We also use the distribution of to generate the mock GRB data, and find that the obtained constraints on the parameters of relation and the matter density parameter are consistent with that by using the distribution of .
We use the fiducial values of the extended Amati relation: , , , , obtained from allredshift GRB data by using the global fitting calibration method in the flat CDM framework. (Different choices of fiducial values do not affect the final results.) Based on these two distribution functions, equations (16) and (19), we could obtain the simulation data for both redshift and peak energy , and then get the simulation data of bolometric fluence according to Eq. (8). For the real GRB data, we find that the uncertainties of and have the mean relative error of and , respectively [33]. Therefore, here we set the simulation data of and have the same relative errors. Furthermore, in the process of simulation, we also include the information of the intrinsic scatter and choose obtained from the calibration method using all current GRB data. Finally, the equation becomes
(20) 
where which includes the intrinsic scatter.
Firstly, we use the current precision (,) to simulate the 138 GRB data to constrain the coefficients of the extended Amati relation and obtain the similar bestfit values and error bars of coefficients with those from the current real GRB data, which means that we could safely use this mock data for further calculations.
Next, we use the lowredshift calibration method to constrain the matter density parameter from this mock GRB data. In practice, we select the lowredshift GRB samples with to constrain those four coefficients of the extended Amati relation. The results are shown in Fig. 6. The obtained bestfit values perfectly recover the fiducial ones of parameters of GRB relation. In section 3.2, we have proven that a small shift of the parameters of “calibrated” relation could significantly change the final constraint on from the real GRB data in the lowredshift calibration method. Therefore, here we want to check this using the mock GRB data. We assume the values of parameters (52.805, 2.408, 0.020, 1.939), which is a small shift from the best fit fiducial values and safely lay in the contour (red points in Fig. 6), to calibrate the extended Amati relation. Then we use the mock highredshift GRB data to constrain the matter density parameter, based on this “calibrated” relation, and obtain the 68% C.L. constraint:
(21) 
which is far away from the fiducial value . This result tells us that even a small shift of coefficients in the contour, the constraint on the matter density parameter with the current precision of GRB data could still be significantly biased using the lowredshift calibration method. The constraint on is very sensitive to the values of parameters for the calibration, when we do not fully take the strong correlations among parameters into account and underestimate the statistical errors in the lowredshift calibration method.
In Fig 7 we show the obtained biased values of matter density parameter from the mock GRB data with different precisions using the lowredshift calibration method. We can see that when GRB data have larger samples with the current precision or the relative errors of GRB data decreases, will be more and more close to the fiducial value.
4. Conclusions and discussions
Thanks to the extremely high power of the explosion, Gammaray bursts (GRBs) are proposed as a promising candidate to trace the Hubble diagram of the Universe in high redshift range. In recent years, several popular luminosity correlations have been statistically concluded from GRBs observations. People also have made great efforts to standardize them as cosmological distance indicators. In general, there are two key issues extensively discussed in the literature when we use these correlations for cosmology. The first issue is that whether luminosity correlations evolve with redshift. The other one is about methods to calibrate luminosity correlations. In this paper we take the Amati relation as an example and use the current GRB data in the redshift range [0.0331, 8.1] to investigate these two issues. Here we summarize our main conclusions in more detail:

We divide the whole GRB data into two redshift bins and constrain the coefficients and of the original Amati relation in each bin, respectively. Based on the MCMC method, we find that the constraints on both the intercept and the slope are quite different from the lowredshift and highredshift GRB data, respectively. The tension is more than level, which implies that the parameters of original Amati relation could be redshiftdependent.

We introduce two extra redshiftdependent terms to characterize the redshift evolution of the luminosity correlation. Interestingly, we find that the tension between constraints of parameters from low and high GRB data are alleviated.

Besides the evolution of Amati correlation, the calibration is also very important for using GRB data in cosmological studies. We firstly check the constraint on the matter density parameter using the lowredshift calibration method and obtain (68% C.L.) which is consistent with other works.

However, we could also use the global fitting method to calibrate the GRB data, in which we use all the GRB data to constrain coefficients of the extended Amati relation and the cosmological parameters simultaneously. In this case, we find that, due to the current poor precision of the GRB data, the constraint on is very weak, which is quite different from that by using the lowredshift calibration method.

After our careful checks, we find that the lowredshift calibration method does not take the whole correlations between and coefficients into account. Neglecting the correlation information can break the degeneracies between and coefficients. Therefore, the obtained constraint on is totally biased. A small shift of the parameters of “calibrated” relation could significantly change the final constraint on in the lowredshift calibration method, which implies that the constraint on is very sensitive to the values of parameters for the calibration.

In order to investigate how large biased constraint on by using the lowredshift calibration, we simulate the mock GRB data with different precisions. We find that the mock data with current precision will give significantly biased result using the lowredshift calibration method, even we assume a small shift on the parameters of “calibrated” relation. When GRB data include larger samples with the current precision or the relative errors of GRB data decreases, the constraint on is more and more close to the fiducial value.
Acknowledgements
G.J. Wang and Z.H. Zhu are supported by the the National Science Foundation of China under grant No. 11373014. H. Yu is supported by the National Basic Research Program of China (973 Program, grant No. 2014CB845800). Z.X. Li is supported by the NSFC under grant No. 11505008. J.Q. Xia is supported by the National Youth Thousand Talents Program and the NSFC under grant No. 11422323. The research is also supported by the Strategic Priority Research Program “The Emergence of Cosmological Structures” of the Chinese Academy of Sciences, Grant No.XDB09000000, and the NSFC under grant Nos. 11633001 and 11690023.
References
 Ade et al. [2015] Ade P. A. R., et al. [Planck Collaboration], 2015, arXiv:1502.01589
 Amanullah et al. [2010] Amanullah R., et al., 2010, ApJ, 716, 712
 Amati et al. [2002] Amati L., et al., 2002, A&A, 390, 81
 Band et al. [1993] Band D., et al., 1993, ApJ, 413, 281
 Basilakos & Perivolaropoulos [2008] Basilakos S., Perivolaropoulos L., 2008, MNRAS, 391, 411
 Betoule et al. [2014] Betoule M., et al., 2014, A&A, 568, A22
 Cucchiara et al. [2011] Cucchiara A., et al., 2011, ApJ, 736, 7
 D’Agostini [2005] D’Agostini G., 2005, arXiv:physics/0511182 [physics.dataan]
 Dai, Liang & Xu [2004] Dai Z. G., Liang E. W., Xu D., 2004, Astrophys. J., 612, L101
 Dainotti et al. [2013] Dainotti M. G., Cardone V. F., Piedipalumbo E., Capozziello S., 2013, MNRAS, 436, 82
 Ding, Li & Zhu [2015] Ding X. H., Li, Z. X., Zhu Z. H., 2015, IJMPD, 24, 7
 Fenimore & RamirezRuiz [2000] Fenimore E. E., RamirezRuiz E., 2000, arXiv: astroph/0004176
 Firmani, AvilaReese, & Ghirlanda [2006] Firmani C., AvilaReese V., and Ghirlanda G., 2006, MNRAS, 372, L28
 Firmani et al. [2005] Firmani C., Ghisellini G., Ghirlanda G., AvilaReese V., 2005, MNRAS, 360, L1
 ForemanMackey et al. [2012] ForemanMackey, D., Hogg, D. W., Lang, D. and Goodman, J. 2012, arXiv:1202.3665
 Ghirlanda et al. [2004] Ghirlanda G., Ghisellini G., Lazzati D., Firmani C., 2004, ApJ, 613, L13
 Ghirlanda, Ghisellini & Lazzati [2004] Ghirlanda G., Ghisellini G., Lazzati D., 2004, ApJ, 616, 331
 Ghirlanda et al. [2006a] Ghirlanda, G., Ghisellini, G., & Firmani, C. 2006, New J. Phys, 8, 123
 Guiriec et al. [2011] Guiriec, S., Connaughton, V., Briggs, M. S., et al. 2011, ApJ, 727, L33
 Guiriec et al. [2013] Guiriec, S., Daigne, F., HascoÃ«t, R., et al. 2013, ApJ, 770, 32
 Guiriec et al. [2015a] Guiriec, S., Kouveliotou, C., Daigne, F., et al. 2015, ApJ, 807, 148
 Guiriec et al. [2015b] Guiriec, S., Mochkovitch, R., Piran, T., et al. 2015, ApJ, 814, 10
 Guiriec et al. [2016a] Guiriec, S., Gonzalez, M. M., Sacahui, J. R., et al. 2016, ApJ, 819, 79
 Guiriec et al. [2016b] Guiriec, S., Kouveliotou, C., Asano, K., et al. 2016, arXiv: 1606.07193
 Kumar & Zhang [2015] Kumar P., Zhang B., 2015, Phys. Rep. 561, 1
 Li [2007] Li L. X., 2007, MNRAS, 379, L55
 Li et al. [2008] Li H.,Xia J. Q., Liu J., Zhao G. B., Fan Z. H., and Zhang X. M., 2008, ApJ, 680, 92
 Li, Ding & Zhu [2015] Li Z., Ding X., Zhu Z. H., 2015, PRD, 91, 083010
 Liang & Zhang [2005] Liang E. W., Zhang B., 2005, ApJ, 633, 611
 Liang et al. [2008] Liang N., Xiao W. K., Liu Y., Zhang S. N., 2008, ApJ, 685, 354
 Lin et al. [2015] Lin H.N., Li X., Wang S., Chang Z., 2015, MNRAS, 453, L128
 Lin et al. [2016] Lin H. N., Li X., Chang Z., 2016, MNRAS, 459, 2501
 Liu & Wei [2015] Liu J., Wei H., 2015, Gen.Rel.Grav. 47, 141
 Mészáros [2002] Mészáros P., 2002, ARA&A, 40, 137
 Mészáros [2006] Mészáros P., 2006, Rep. Prog. Phys., 69, 2259
 Norris, Marani & Bonnell [2000] Norris J. P., Marani G. F., & Bonnell J. T., 2000, ApJ, 534, 248
 Perlmutter et al. [1999] Perlmutter S., et al., 1999, Astrophys. J., 517, 565
 Piran [1999] Piran T., 1999, Phys. Rep., 314, 575
 Riess et al. [1998] Riess A. G., et al., 1998, Astron. J., 116, 1009
 Salvaterra et al. [2009] Salvaterra R., et al., 2009, Nature, 461, 1258
 Schaefer [2002] Schaefer, B. E. 2002, ApJ, 583, L67
 Schaefer [2007] Schaefer, B. E. 2007, ApJ, 660, 16
 Suzuki et al. [2012] Suzuki N., et al., 2012, ApJ, 746, 85
 Tanvir et al. [2009] Tanvir N. R., et al., 2009, Nature, 461, 1254
 Tanvir [2013] Tanvir N. R., 2013, arXiv:1307.6156
 Wang, Dai & Liang [2015] Wang F. Y., Dai Z. G., Liang E. W., 2015, New Astronomy Review, 67, 1
 Wang, Li & Zhang [2014] Wang S., Li Y. H., & Zhang X., 2014, PRD, 89, 063524
 Wang, Qi & Dai [2011] Wang F. Y., Qi S., Dai Z. G., 2011, MNRAS, 415, 3423
 Wang & Wang [2013] Wang S., Wang Y., 2013, PRD, 88, 043511
 Wang et al. [2016] Wang J. S., Wang F. Y., Cheng K. S., Dai Z. G., 2016, A&A, 585, A68
 Wei [2010] Wei H., 2010, J. Cosmol. Astropart. Phys., 08, 020
 Wei & Zhang [2009] Wei H., Zhang S. N., 2009, Eur. Phys. J. C, 63, 139
 Yonetoku et al. [2004] Yonetoku D., Murakami T., Nakamura T., Yamazaki R., Inoue A. K., Ioka K., 2004, ApJ, 609, 935
 Yu et al. [2015] Yu, H., Wang, F. Y., Dai, Z. G., & Cheng, K. S. 2015, ApJS, 218, 13