Systematics in the GRBs Hubble diagram

# Systematics in the Gamma Ray Bursts Hubble diagram

V.F. Cardone111Corresponding author : winnyenodrac@gmail.com, M. Perillo, S. Capozziello
I.N.A.F. - Osservatorio Astronomico di Roma, via Frascati 33, 00040 - Monte Porzio Catone (Roma), Italy
Dipartimento di Fisica “E.R. Caianiello”, Università di Salerno, Via Ponte Don Melilo, 84081 Fisciano (Sa), Italy
I.N.F.N. - Sez. di Napoli, Compl. Univ. Monte S. Angelo, Ed. G, Via Cinthia, 80126 Napoli, Italy
Dipartimento di Scienze Fisiche, Università degli Studi di Napoli Federico II”, Complesso Universitario
di Monte Sant’Angelo, Edificio N, via Cinthia, 80126 - Napoli, Italy
Accepted xxx, Received yyy, in original form zzz
###### Abstract

Thanks to their enormous energy release which allows to detect them up to very high redshift, Gamma Rays Bursts (GRBs) have recently attracted a lot of interest to probe the Hubble diagram (HD) deep into the matter dominated era and hence complement Type Ia Supernoave (SNeIa). However, lacking a local GRBs sample, calibrating the scaling relations proposed as an equivalent to the Phillips law to standardize GRBs is not an easy task because of the need to estimate the GRBs luminosity distance in a model independent way. We consider here three different calibration methods based on the use of a fiducial CDM model, on cosmographic parameters and on the local regression on SNeIa. We find that the calibration coefficients and the intrinsic scatter do not significantly depend on the adopted calibration procedure. We then investigate the evolution of these parameters with the redshift finding no statistically motivated improvement in the likelihood so that the no evolution assumption is actually a well founded working hypothesis. Under this assumption, we then consider possible systematics effects on the HDs introduced by the calibration method, the averaging procedure and the homogeneity of the sample arguing against any significant bias. We nevertheless stress that a larger GRBs sample with smaller uncertainties is needed to definitely conclude that the different systematics considered here have indeed a negligible impact on the HDs thus strengthening the use of GRBs as cosmological tools.

###### keywords:
gamma - rays burst : general – distance scale – cosmological parameters

## 1 Introduction

The observational evidences accumulated in these less than fifteen years, from the anisotropy and polarization spectra of the cosmic microwave background radiation (CMBR, de Bernardis et al. (2000); Brown et al. (2009); Komatsu et al. (2011), the large scale structure traced by galaxy redshift surveys (Dodelson et al., 2002; Percival et al., 2002; Szalay et al., 2003; Hawkins et al., 2003), the matter power spectrum (Tegmark et al., 2006; Percival et al., 2007) with the imprints of the Baryonic Acoustic Oscillations (BAO, Eisenstein et al. (2005); Percival et al. (2010)) and the Hubble diagram of TypeIa Supernovae (SNeIa, Kowalski et al. (2008); Hicken et al. (2009); Kessler et al. (2009)), definitely support the cosmological picture of a spatially flat universe with a subcritical matter content and undergoing a phase of accelerated expansion. While the observational scenario is now firmly established, the search for the motivating theory is, on the contrary, still in its infancy notwithstanding the many efforts and the plethora of models proposed along these years. Ironically, the problem here is not the lack of a well established theory, but the presence of too many viable candidates, ranging from the classical cosmological constant (Carroll et al., 1992; Sahni & Starobinsky, 2000) to scalar fields (Peebles & Ratra, 2003; Copeland et al., 2006) and higher order gravity theories (Capozziello & Francaviglia, 2008; Sotiriou & Faraoni, 2010; Nojiri & Odintsov, 2008; De Felice & Tsujikawa, 2010), all of them being more or less able to fit the available data.

As often in science, adding further data is the best strategy to put order in the confusing abundance of theoretical models. In particular, pushing the observed Hubble diagram to higher redshift would allow to trace the universe background evolution up to the transition regime from the dark energy driven speed up to the decelerated matter dominated phase. Moreover, being the distance modulus log - linearly related to the luminosity distance and depending this latter on the dark energy equation of state through a double integration, one has to go to large in order to discriminate among the prediction of different models when these predict similar curves at lower redshift. Unfortunately, SNeIa are not ideally suited to this task with their present day Hubble diagram going back to and not extending further than even for excellent space based projects such as SNAP (Aldering et al., 2004).

Thanks to their enormous almost instantaneous energy release, Gamma Ray Bursts (GRBs) stand out as ideal candidates to go deeper in redshift, the farthest one yet being at (Salvaterra et al., 2009). The wide range spanned by their peak energy makes them everything but standard candles, but the existence of many observationally motivated correlations between redshift dependent quantities and rest frame properties (Amati et al., 2008; Fenimore & Ramirez - Ruiz, 2000; Norris et al., 2000; Ghirlanda et al., 2004; Liang & Zhang, 2005) offers the intriguing possibility of turning GRBs into standardizeable candles just as SNeIa. The use of these scaling relations allows to infer the GRB distance modulus with an error mainly depending on the correlation intrinsic scatter. Combining the estimates from different correlations, Schaefer (2007) first derived the GRBs Hubble diagram (hereafter, HD) for 69 objects, while Cardone et al. (2009) used a different calibration method and add a further correlation to update the GRBs HD. Many attempts on using GRBs as cosmological tools have since then been performed (see, e.g., Firmani et al. 2006, Liang et al. 2009, Qi & Lu 2009, Izzo et al. 2009 and refs. therein) showing the potential of GRBs as cosmological probes.

It is worth stressing that the possibility offered by GRBs to track the HD deep into the matter dominated epoch does not come for free. Two main problems are actually still to be fully addressed. First, lacking a local GRBs sample, all the above correlations have to be calibrated assuming a fiducial cosmological model to estimate the redshift dependent quantities. As a consequence, the so called circularity problem comes out, that is to say one wants to use GRBs scaling relations to constrain the underlying cosmology, but needs the underlying cosmology to get the scaling relations. Different strategies have been proposed to break this circularity so that our first aim here is to investigate whether they are indeed viable solutions and to what extent the residual problem impact the HD derivation.

A well behaved distance indicator should be not only visible to high and possess scaling relations with as small as possible intrinsic scatter, but its physics should be well understood from a theoretical point of view. On the contrary, there is up to now not any definitive understanding of the GRBs scaling relations so that, as a dangerous drawback, one can not anticipate whether the calibration parameters are redshift dependent or not. Lacking any theoretical model, we will therefore address here this problem in a phenomenological way adopting two different parameterizations for their evolution with allowing for a change in the zeropoint only or in both zeropoint and slope. It is worth noting that such an analysis has been severely hampered before by the low statistics in the GRBs sample, a problem that we partially overcome here thanks to the use of the 115 GRBs catalog recently assembled by Xiao & Schaefer (2010).

The plan of the paper is as follows. In Sect. 2, we briefly review the different GRBs 2D correlations known insofar and present the Bayesian motivated method we will use in the following to calibrate them using three different approaches to estimate the GRBs luminosity distances. The constraints thus obtained and their dependence on the adopted luminosity distance determination are discussed in Sect. 3, while, in Sect. 4, we address the problem of their dependence on the redshift. The issues related to the HD derivation are then discussed in Sect. 5 where we consider the impact on the HD of the luminosity distance estimate method, of the the averaging procedure and the homogeneity of the sample. A summary of the main results and a discussion of them is finally given in Sect. 6, while some complementary material is presented in Appendix.

## 2 GRBs scaling relations

It is instructive to start this analysis considering the general case of two observable quantities related by a power - law relation which, in a log - log plane, reads

 logy=alogx+b . (1)

Calibrating such a relation means determining the slope , the zeropoint and the scatter of the points around the best fit relation. In order to be useful for distance determination, one of the two quantities, say , should refer to a directly observed quantity, while the other one must depend on the redshift . Setting with a directly measurable redshift independent quantity and the luminosity distance, we get

 logy=logκ+2logdL(z)=alogx+b

so that one can then estimate the distance modulus as :

 μ=25+5logdL(z)=25+(5/2)(alogx+b−logκ) .

In order to perform such an estimate, a two steps procedure has to be implemented. First, one has to select a sample of low redshift objects with known distance and fit the scaling relation to the data thus estimated to infer the calibration parameters . Second, one has to assume that such calibration parameters do not change with the redshift so that a measurement of and the use of the above scaling relation, with the same at all , are sufficient to infer the distance modulus.

Both these steps are daunting tasks for GRBs because of the lack of low redshift objects to be used in the calibration procedure. In order to overcome this problem, different strategies can be implemented, but their impact on the final derivation of the GRBs Hubble diagram have not been investigated in detail which is our aim here.

### 2.1 2D empirical correlations

Their high luminosity offering the possibility to be detected at very large makes GRBs promising candidates to trace the Hubble diagram in the matter dominated era so that a great effort has been devoted during the late years to look for reliable and narrow empirical correlations. We limit here our attention only to two dimensional (hereafter, 2D) correlations since they can be investigated relying on a larger number of GRBs. These involve a wide range of GRBs properties related to both the energy spectrum and the light curve which are then correlated with the isotropic luminosity or the emitted collimation corrected energy . Neither nor are directly measurable quantities since they depend on the luminosity distance . Indeed, it is :

 L=4πd2L(z)Pbolo , (2)
 Eγ=4πd2L(z)SboloFbeam(1+z)−1 , (3)

where and are the bolometric peak flux and fluence, respectively, while is the beaming factor with the rest frame time of achromatic break in the afterglow light curve. Note that the bolometric quantities are related to the observed ones as (Schaefer, 2007) :

 Pbolo=P × ∫104/(1+z)1/(1+z)EΦ(E)dE∫EmaxEminEΦ(E)dE , (4)
 Sbolo=S × ∫104/(1+z)1/(1+z)EΦ(E)dE∫EmaxEminEΦ(E)dE , (5)

with and the observed peak energy and fluence and the detection thresholds of the observing instrument. We model the energy spectrum as a smoothly broken power - law, i.e. (Band et al., 1993) :

 Φ(E)=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩AEαexp[−(2+α)EEpeak]EEpeak≤α−β2+αBEβotherwise (6)

The GRBs 2D correlations are obtained setting or in Eq.(1), while different choices are available for . In particular, we will consider the following possibilities : (i.) with the peak energy (in keV) of the fluence spectrum; (ii.) with (in s) the time offset between the arrival of the low and high energy photons; (iii.) with (in s) the shortest time over which the light curve rises by half the peak flux of the burst; (iv.) with the variability which quantifies the smoothness of the light curve itself. Note that all these quantities are expressed in the GRB rest frame which motivates the redshift term, while the further scaling constant is introduced to minimize the correlation among errors.

The combination of and gives rise to the different correlation we will consider222Hereafter, we will refer to the correlation obtained setting and in Eq.(1) as the  -  correlation or scaling law., namely the  -  (Ghirlanda et al., 2004; 2006, )), the  -  (Schaefer, 2003; Yonekotu et al., 2004),  -  (Norris et al., 2000),  -  (Schaefer, 2007) and  -  (Fenimore & Ramirez - Ruiz, 2000; Riechart et al., 2001; Schaefer, 2007). These five correlations has been used in Schaefer (2007, hereafter S07) to derive a combined Hubble diagram for 69 GRBs, later updated by Cardone et al. (2009, hereafter CCD09) where also a sixth correlation between the break time and the luminosity at break time of the X - ray afterglow (Dainotti et al., 2008, 2010, 2011) has been considered.

Recently, a larger GRBs catalog has been compiled by Xiao & Schaefer (2009, hereafter XS09) giving, for each object, the values of the different quantities (if available) needed to check the five correlations above. We use this sample here as input for our analysis referring the interested reader to XS10 for details on the catalog construction and observable quantities determination.

### 2.2 Bayesian fitting procedure

Eq.(1) is a linear relation which can be fitted to a given dataset in order to determine the two calibration parameters . Moreover, although there is still not a theoretical model explaining any of the empirical 2D correlations in terms of GRB physics, we nevertheless expect that the wide range of GRB properties makes the objects scatter around this (unknown) idealized model. As a consequence, the above linear relations will be affected by an intrinsic scatter which has to be determined together with the calibration coefficients . To this aim, in the following we will resort to a Bayesian motivated technique D’ Agostini (2005) thus maximizing the likelihood function with :

 L(a,b,σint) = 12∑ln(σ2int+σ2Yi+a2σ2Xi) (7) + 12∑(Yi−aXi−b)2σ2int+σ2Yi+a2σ2Xi

with and the sum is over the objects in the sample. Note that, actually, this maximization is performed in the two parameter space since may be estimated analytically as :

 b=[∑Yi−aXiσ2int+σ2Yi+a2σ2Xi][∑1σ2int+σ2Yi+a2σ2Xi]−1 (8)

so that we will not consider it anymore as a fit parameter.

It is worth noting that the Bayesian approach allows to find out what is the most likely set of parameters within a given theory, but does not tell us whether this model fits well or not the data. An easy way to quantitatively estimate the goodness of the fit is obtained considering the median and root mean square of the best fit residuals, defined as which we will also compute for the different 2D correlations we will consider.

The Bayesian approach used here also allows us to quantify the uncertainties on the fit parameters. To this aim, for a given parameter , we first compute the marginalized likelihood by integrating over the other parameter. The median value for the parameter is then found by solving :

 ∫pi,medpi,minLi(pi)dpi=12∫pi,maxpi,minLi(pi)dpi . (9)

The () confidence range are then found by solving :

 ∫pi,medpi,lLi(pi)dpi=1−ε2∫pi,maxpi,minLi(pi)dpi , (10)
 ∫pi,hpi,medLi(pi)dpi=1−ε2∫pi,maxpi,minLi(pi)dpi , (11)

with (0.95) for the () range respectively. Actually, in order to sample the parameter space, we use a Markov Chain Monte Carlo (MCMC) method running two parallel chains and using the Gelman - Rubin (1992) test to check convergence. The confidence ranges are then obtained considering the histograms of the parameters from the merged chain after burn in cut and thinning.

### 2.3 GRBs luminosity distances

A preliminary step in the analysis of all the 2D correlations hinted at above is the determination of the luminosity or the collimated energy entering as variable in the  -  scaling laws. As shown by Eqs.(2) - (3), one has to first determine the GRBs luminosity distance over a redshift range where the linear Hubble law does not hold anymore. As such, one should estimate the luminosity distance as :

 dL(z)=c(1+z)H0∫z0dz′E(z′) (12)

with the present day Hubble constant and the dimensionless Hubble parameter which depends on the adopted cosmological model thus leading to the well known circularity problem (i.e., one would like to use GRBs to probe the cosmological model, but needs the cosmological model itself to get the GRBs Hubble diagram).

Different strategies have been developed to tackle this problem. The simplest one is to assume a fiducial cosmological model and determine its parameters by fitting, e.g., the SNeIa Hubble diagram. Motivated by its agreement with a wide set of data, one usually adopt the CDM model as fiducial one thus setting :

 E2(z)=ΩM(1+z)3+ΩΛ (13)

with because of the spatial flatness assumption. In order to determine the parameters , we maximize the likelihood function with :

 χ2(ΩM,h) = NSNeIa∑i=1[μobsi−μth(zi,ΩM,h)σμi]2 (14) + (ωobsM−ωthMσωM)2+(hobs−hσh)2 .

As input data, we use the Union2 SNeIa sample (Amanullah et al., 2010) to get for objects over the redshift range and set for the matter physical density and for the Hubble constant in agreement with the results of the WMAP7 (Komatsu et al., 2011) and SHOES (Riess et al., 2009) teams, respectively. The best fit values turn out to be , while median value and confidence ranges read :

 ΩM=0.259+0.019 +0.032−0.016 −0.028  ,  h=0.723+0.025 +0.046−0.025 −0.045  .

Although the CDM model fits remarkably well the data, it is nevertheless worth stressing that a different cosmological model would give different values for and hence different values for thus impacting the estimate of the calibration parameters . Although CCD09 have shown that this effect is quite small for a large range of phenomenological dark energy models333CCD09 adopted the CPL (Chevallier & Polarski, 2001; Linder, 2003) parametrization of the dark energy equation of state, , and explored the dependence of on for different values of . While clear trends of with were obtained, the relative change is negligibly small. It is worth stressing, however, that the calibration parameters may change more if the adopted cosmological model is radically different from the CDM one (see, e.g, Diaferio et al. 2011 for an illuminating example)., it is nevertheless worth be conservative and look for model independent approaches. A first step towards this aim is to resort to cosmography (Weinberg, 1972; Visser, 2004), i.e., to expand the scale factor to the fifth order and then consider the luminosity distance as a function of the cosmographic parameters. Indeed, such a kinematic approach only relies on the validity of assumption of the Robertson - Walker metric, while no assumption on either the cosmological model or the theory of gravity is needed since the Friedmann equations are never used. We use the expansion of to the fifth order in found in Capozziello et al. (2008) and determine the Hubble constant , the deceleration , jerk , snap and lerk parameters by fitting the Union2 SNeIa data with the prior on coming from the SHOES team quoted above. Note that, in this case, the likelihood function is defined as before, but the pseudo -  is now given by Eq.(14) without the term depending on since now is no more a parameter. Using a similar MCMC algorithm (but now running three parallel chains), we find as best fit values :

 (h,q0,j0,s0,l0)=(0.741,−0.56,0.66,−0.41,3.59) ,

while median values and confidence ranges read :

 h=0.741+0.036 +0.071−0.035 −0.072  ,  q0=−0.45+0.05 +0.10−0.05 −0.14  ,
 j0=0.00+0.13 +0.79−0.12 −0.35  ,  s0=0.29+0.74 +1.57−0.20 −1.08  ,
 l0=−0.07+4.77 +12.27−4.71 −10.86  ,

in good agreement with previous results in literature (Vitagliano et al., 2010; Xu & Wang, 2010; Bouhmadi - López et al., 2010) using different datasets.

As a further step towards a fully model independent estimate of the GRBs luminosity distances, one can use SNeIa as distance indicator based on the naive observations that a GRBs at redshift must have the same distance modulus of SNeIa having the same redshift.Interpolating therefore the SNeIa Hubble diagram gives therefore the value of for a subset of the GRBs sample with which can then be used to calibrate the 2D correlations (Kodama et al., 2008; Liang et al., 2008; Wei & Zhang, 2009). Assuming that this calibration is redshift independent, one can then build up the Hubble diagram at higher redshifts using the calibrated correlations for the remaining GRBs in the sample. CCD09 have used a similar approach based on the local regression technique (Cleveland, 1979; Cleveland & Devlin, 1988; Loader, 1999) which combines much of the simplicity of linear least squares regression with the flexibility of nonlinear regression. The basic idea relies on fitting simple models to localized subsets of the data to build up a function that describes the deterministic part of the variation in the data, point by point. Actually, one is not required to specify a global function of any form to fit a model to the data so that there is no ambiguity in the choice of the interpolating function. Indeed, at each point, a low degree polynomial is fit to a subset of the data containing only those points which are nearest to the point whose response is being estimated. The polynomial is fit using weighted least squares with a weight function which quickly decreases with the distance from the point where the model has to be recovered. We apply local regression to estimate the distance modulus from the Union2 SNeIa sample following the steps schematically sketched in CCD09 which we refer the reader to for further details and the demonstration of the reliability of the inferred luminosity distances.

## 3 Calibration parameters

While the quantities are directly observed for each GRBs, the determination of (either the luminosity or the collimated energy ) needs for the object’s luminosity distance. The three methods described above allows us to get three different values for so that it is worth investigating whether this has any significant impact on the calibration parameters for the correlations of interest. We will refer hereafter to the three samples with the quantities estimated using the luminosity distance from the fiducial CDM cosmological model, the cosmographic parameters and the local regression method as the , and samples, respectively.

As a preliminary caveat, it is worth stressing that the error on comes from two different contributions. First, there is the uncertainty obtained by propagating those on the observed quantities (i.e., the flux and the Band parameters for and the fluence , the jet angle and again the Band parameters for ) which we assume to be uncorrelated. Second, there is the uncertainty on which we estimate from the MCMC chain for the and cases, while it is output from the method for the case. It is worth noting that, given the large SNeIa sample used in the fit of CDM and cosmographic parameters and in the local regression technique, the error on the luminosity distance is actually quite small and always smaller than the one from the measurement uncertainties. It is this latter term that one should minimize in order to get better determined and/or values and hence put stronger constraints on the slope, zeropoint and intrinsic scatter.

In order to constrain the calibration parameters , we can use the Bayesian procedure described above with the three GRBs samples as input to the likelihood analysis. However, the , and samples can not contain the same number of objects since the GRBs luminosity distance (and hence the quantities) can be estimated only for with depending on the adopted method. If we rely on the fiducial CDM model, we can predict at every so that with the only caveat that one is implicitly assuming that a model fitted over the range probed by the Union2 SNeIa sample can be extrapolated to the full evolutionary history of the universe. The sample is based on the use of the cosmographic parameters fitted to the data having checked that the fifth order expansion of is reliable over the Union2 SNeIa redshift range. However, the larger is , the higher is the order one has to include in the Taylor series to get accurate results so that one must set in order to not bias the determination of the quantities because of an inaccurate distance approximation. Finally, the local regression method can only be applied to estimate over the redshift range probed by the data used for the interpolation so that one get the constraint . In order to homogenize the three samples, we therefore fit the 2D correlations using only , and GRBs with so that we can meaningfully compare the results for the three distance estimate methods.

Table 1 summarizes the constraints on the calibration parameters for the three different samples. Note that we have also included the  -  correlation (Amati et al., 2009) with since, although refers to the same quantities used in the  -  correlation, it is based on a larger number of GRBs (given that one does not need a measurement of the jet opening angle ). As a general result, we find that the fit is always quite good with reduced values close to in all the cases independently on the 2D correlation considered and the distance estimate method adopted. As such, the results from the different sample are statistically equivalent and can be safely compared.

Although the and confidence ranges for the calibration parameters for a given correlation overlap quite well for the , and samples, the best fit coefficients and the median values clearly show that the calibration based on the fiducial CDM model leads to steeper scaling laws for most of the cases. On the contrary, shallower slopes are obtained using the or samples with the  -  relation as unique exception444Actually, the  -  relation is the shallowest one and has the largest intrinsic scatter among the six 2D correlations investigated in this paper so that a different trend of the slope with the sample could also be a statistical artifact.. Although the differences in the slopes are not statistically meaningful because of the large uncertainties, it is nevertheless worth investigating whether such an effect can be ascribed to the distance determination. To this end, for a given GRB and  -  correlation, let us define the following quantities :

 ΔXY=(a1X+b1)−(a2X+b2)=(a1−a2)X+(b1−b2) ,
 ΔdLY=2log[dL1(z)/dL2(z)] .

While gives the difference in the values predicted using the correlations with coefficients and , quantifies the effect of estimating the model dependent quantity using two different luminosity distances. Should the change in the slope be the outcome of compensating the offset due to different luminosity distances assumptions, one should get . Actually, this is not the case. Considering, e.g, the  -  correlation, a weighted average gives :

 ⟨ΔXY⟩=0.19≠⟨ΔdLY⟩=0.09

using the and samples and

 ⟨ΔXY⟩=−0.25≠⟨ΔdLY⟩=−0.02

for the vs samples. We therefore find that the change in the slope is not induced by the different luminosity distances adopted. Actually, the use of different samples has also an impact on the intrinsic scatter determination with the sample leading to smaller best fit and median values for all the six correlations considered. Since and correlate, the change in the slope is not only due to the change in the luminosity distances, but also of the intrinsic scatter. Discriminating among these two effects and quantifying their respective contribution is not possible so that we can not draw any definitive conclusion on the impact of the calibration method on the slope of the 2D scaling relations.

## 4 Evolution with redshift

As more and more data add to the available GRBs dataset, the observational evidences for the six 2D correlations we are considering become more and more reliable. On the contrary, a clear theoretical motivation is still lacking in many cases. As a consequence, it is also not clear whether the calibration parameters evolve with the redshift or not. To investigate this issue, we consider two different possibilities for the evolution with . First, we consider the possibility that the slope is constant, but the zeropoint is evolving. In particular, we assume :

 y=B(1+z)αxA ⟶ Y=αlog(1+z)+aX+b (15)

with and . Table 2 summarizes the constraints on the parameters obtained fitting Eq.(15) to the GRBs with distances estimated from the fiducial CDM model. The constraints from the fit to the and samples are consistent within the confidence ranges so that they are reported in Appendix for completeness, but not discussed anymore here.

Comparing to the constraints in Table 1 highlights some interesting lessons. First, we note that both the best fit and median values of the slope parameter are significantly shallower than in the no evolution case. However, the confidence ranges typically overlap quite well so that, from a statistical point of view, such a result should not be overrated. Actually, adding one more parameter introduces a degeneracy between and so that one can make the best fit  -  relation shallower compensating the difference in the term with the contribute from . Indeed, is typically quite large the only exception being the  -  correlation which is also the only one with the same best fit slope for the fit with and without the evolution term. Quite surprisingly, the no evolution result (i.e., ) is consistent with the confidence range only for the  -  correlation thus arguing in favor of an evolution of the calibration parameters with . On the other hand, the reduced values are never smaller than those obtained for the fits with no evolution. Moreover, the intrinsic scatter is almost the same for both cases so that allowing for a redshift evolution does not lead to statistically preferred results. As such, we consider a most conservative option to assume that the GRBs scaling relations explored here do not evolve with .

Actually, such a conclusion is model dependent. As an alternative parametrization, we therefore allow for an evolution of the slope and not only the zeropoint of the 2D correlations. We fit the data using :

 Y=(a0+a1z)X+(b0+b1z) , (16)

i.e., we are Taylor expanding to first order the unknown dependence of the slope and zeropoint on the redshift. Note that, actually, we should limit the redshift range to very low to have a meaningful expansion, but we extrapolate this linear relation to any as a first order guess avoiding to add further unknown fitting parameters. Morevoer, since we expect are quite small, we skip to logarithmic units for these quantities in order to be more sensitive to the tiny values. Table 3 summarizes the results for the fit to the sample, while a qualitatively similar conclusions can be drawn from the those to the and samples as can be inferred from Table 6 in Appendix A.

As a general result, we find that the best fit parameters and the median values of the evolutionary coefficients are typically quite small indicating that the dependence of both the slope and the zeropoint on the redshift is quite weak, if present at all555Note that, since we are using logarithmic units, the no evolution case can not be exactly achieved corresponding to and going to . However, needless to say, obtaining, e.g., is the same as stating that there is no evolution at all. Note also that we have implicitly assumed that both and are positive. We have, however, checked that the qualitative conclusions are not affected by this assumption. . On the other hand, the slope and zeropoint at are consistent within the confidence range with the corresponding quantities in the no evolution case. Moreover, the reduced values are typically larger than in the no evolution case so that this latter assumption is statistically preferred.

## 5 GRBs Hubble diagram

Once the calibration parameters for a given  -  correlation have been obtained, it is then possible to estimate the distance modulus of a given GRB from the measured value of . Indeed, for a given , the luminosity distance is :

 d2L(z)=y/κ

with , and for , and , respectively. Using the definition of distance modulus and estimating from through the  -  correlation, we then get :

 μ(z) = 25+5logdL(z) (17) = 25+(5/2)log(y/κ) = 25+(5/2)(aX+b−logκ)

where are the best fit coefficients for the given  -  correlation. Note that we will refer here only to the fits with no redshift evolution of the calibration parameters since the analysis in the previous section has not shown any clear evidence for a redshift dependence of . Although such a result is model dependent (since we have considered only two possible evolution parameterization) and actually limited only to the redshift range , we are nevertheless confident that the possible bias induced by the no evolution assumption is smaller than the present day measurement uncertainties. However, such an issue worths to be reconsidered when a larger and more precise sample will be available.

Eq.(17) allows to compute the central value of the distance modulus relying on the measured values of the GRBs observables, i.e., the ones entering the quantity , and the best fit coefficients of the used correlation. However, both and are known within their own uncertainties which have to be propagated to get the error on . Moreover, each correlation is affected by an intrinsic scatter which has to be taken into account into the total error budget. To this end, we adopt the procedure schematically sketched below.

1. Fix to the  - th value of the Markov Chain and estimate from Eq.(8).

2. Set the quantities needed to compute and hence randomly sampling from Gaussian distributions centred on the measured observed values and with variance equal to the measurement error.

3. Compute from Eq.(17) for all the 2000 values generated above and estimate the mean of the distribution.

4. Repeat steps and for all the points along the chain and estimate the statistical error on by symmetrizing the confidence range of the distribution thus obtained.

5. Add in quadrature the statistical error and the intrinsic scatter of the correlation to finally get the total uncertainty on the distance modulus .

Eq.(17) and the above procedure allows us to build up the Hubble diagram for all the GRBs with measured values of the quantities entering a given  -  correlation. It is then possible to both reduce the uncertainties and (partially) wash out the hidden systematic errors by averaging over the different correlations available for a given GRB. Following Schaefer (2007), we finally estimate the distance modulus for the  - th GRB in the sample at redshift as :

 μ(zi)=[∑jμj(zi)/σ2μj] × [∑j1/σ2μj]−1 (18)

with the uncertainty given by :

 σμ=[∑j1/σ2μj]−1 (19)

where the sum runs over the 2D empirical laws which can be used for the GRB considered.

As summarized in Table 1 and discussed in Sect. 3, the best fit calibration parameters depend on the method used to set the GRBs luminosity distances. Although such a dependence is not statistically meaningful being the confidence ranges for the different cases well overlapped, it is nevertheless worth investigating whether they have an impact on the derived Hubble diagram. Moreover, averaging over more than one correlation implicitly assumes that all of them are physically motivated (i.e., they are the outcome of an unknown underlying theoretical mechanism) and not the artifact of some not well understood selection effect. It is therefore also important to check what should be the effect of a possible misleading assumption. Both these issues will be discussed in the following.

### 5.1 Impact of the calibration method

Fig. 1 shows the GRBs Hubble diagrams (hereafter, HDs) obtained averaging over the six 2D correlations and using the three different calibration methods. As a guidance for the eye, the red solid line is the expected curve for the fiducial CDM model. Note that we have excluded GRBs which deviate from this line more than with the rms of for the full sample. Such a criterium is actually quite loose excluding only two objects out of 114 so that we are confident that no bias is induced by this cut.

As a general remark, we find that, notwithstanding the calibration method adopted, the GRBs HDs reasonably follow the CDM curve although with a non negligible scatter. Quite surprisingly, the scatter is significantly larger in the range because of a set of GRBs with lying systematically above the CDM prediction. One should argue for a failure of the theoretical model, but there are actually a set of points which are hard to reconcile with any reasonable dark energy model. While this could be dropped out from the sample by an iterative selection procedure, we prefer to err on the conservative side in order to introduce any bias induced by an a priori choice of a dark energy model. It is nevertheless likely that these GRBs are outliers of one or more of the 2D scaling laws used to infer their distance modulus. Looking in detail to their estimates, we find that the estimates obtained from the different correlations entering the averaging procedure are quite different, while this is not for the GRBs less deviating from the red line. Such a naive observation makes us argue in favor of a problem with the measurement of the quantities entering the correlations of interest or of a deviation from one of the best fit scaling relations. Investigating in detail this issue on a case - by - case basis would need to retrieve the original data for each GRB which is outside the aim of our work.

In order to compare the HDs from the three different calibration methods, we consider the values of with the theoretically predicted distance modulus for the fiducial CDM model. We find :

 ⟨Δμ⟩=−0.15  ,  (Δμ)med=0.07  ,  (Δμ)rms=1.15  ,
 ⟨Δμ⟩=−0.10  ,  (Δμ)med=0.16  ,  (Δμ)rms=1.15  ,
 ⟨Δμ⟩=−0.24  ,  (Δμ)med=−0.01  ,  (Δμ)rms=1.21  ,

for the , and samples, respectively. One could naively be surprised that the sample (whose calibration is based on the same fiducial CDM model used as reference here) does not give systematically smaller deviations. Actually, the calibration is made using only GRBs with , but the full sample is dominated by higher objects so that the mean and median of are not biased by the choice of the reference model. As a consequence, we find the rms is roughly the same for the three samples, although the mean values suggest that the calibration based on the linear regression leads to values that are preferentially larger than expected. From a different point of view, higher values at high argue in favor of a model where the transition from deceleration to acceleration takes place to a larger , which can be achieved (for a dark energy model with constant equation of state) pushing to values smaller than -1 (i.e., going to the phantom regime) or decreasing the matter content. However, such a naive argument should be reconsidered taking care of the uncertainties on the single GRBs distance modulus estimates.

Actually, from the point of view of the impact of calibration methods, we are not interested in minimizing the mean or median values of , but rather on comparing the distributions among them. Should these distributions be equal, then one could conclude that the HDs based on different calibration methods are consistent with each other. The different mean and median values seem to suggest that this is not the case. However, these distributions are quite asymmetric and wide enough to allow for a considerable overlap among them. Moreover, for each single GRB, the values of are consistent with each other within the uncertainties. We can therefore conclude that the HDs obtained using different calibration methods are consistent with each other within the uncertainties.

We, nevertheless, warn the reader that such a result is partially due to the large error bars. Should these be reduced, the consistency among the three different HDs could also gone lost. To this end, however, one should first reduce the intrinsic scatter of the correlations used in the average procedure since it is this term which typically dominates the error budget. At the present time, one should give off the correlations with the larger values (e.g, the  -  one), but this would lower the number of objects in the sample and increase the error because of averaging over a lower number of estimates. Larger GRBs samples with measured redshift and observational parameters are then needed to address the issue of the impact of calibration methods on a safer basis.

### 5.2 Impact of the averaging procedure

As yet stated above, averaging the values from different correlations helps reducing the total uncertainties and partially washes out the systematics connected to each single scaling relations. However, this does not come with no price to pay. Indeed, such a strategy implicitly assumes that all the scaling relations are physical ones and not the artifact of a not well controlled selection effect. This danger is actually hard to avoid since no clearly defined and well accepted theoretical model is presently available to describe the physics of GRBs and predict scaling relations in agreement with the empirically motivated ones.

As a first check, we compare the values obtained estimating using each single correlation. For the sample, we get the results summarized in Table 4 where we also give the number of GRBs used. The conclusions are unaltered for the and samples so that we do not discuss them.

While the median values of are roughly comparable, both and are definitely larger for the  -  and  -  correlations. It is noteworthy that these are also the only two correlations where all the GRBs can be used to infer the HD so that one could ascribe the larger to the inclusion of a higher number of outliers. On the other hand, this same effect can be a statistical artifact due to the use of a smaller GRBs sample. It is actually hard to understand whether a selection effect is at work here. Indeed, the inclusion of a GRB in a sample used for a given correlation depends mainly on observational requirements (e.g., the GRBs afterglow should last enough to measure the rise time) so that one could argue that this has nothing to do with its departure from a theoretical quantity such as the fiducial distance modulus. On the other hand, should selection effects induce a fake correlation, one can imagine that the deviations from the fiducial curve are randomly distributed leading to small .

Pending the questio of which relation is physical, we can quantify the impact of an incorrect assumption by evaluating again the distance moduli excluding the  -  and  -  correlations. Since these are the ones with the greatest rms and the largest intrinsic scatter (see Table 1), we expect them to have the greatest impact on the departure from the fiducial CDM theoretically predicted HD. We find :

 ⟨Δμ⟩=−0.11  ,  (Δμ)med=0.04  ,  (Δμ)rms=1.05  ,

for 106 GRBs, being the number smaller than before because some GRBs with determined from the  -  and/or  -  correlations only now drop out of the final sample. Compared with the case with all correlations included, the mean, median and rms are indeed smaller. However, since we average now on a lower number of correlations, the uncertainty on for each GRB is now larger so that the values with and without these two correlations are consistent within the  confidence ranges.

We must therefore conclude that, given the available GRBs dataset, including or not a given correlation in the averaging procedure should be a compromise between the need to avoid an uncertain systematic bias and the desire to ameliorate both statistics and precision.

### 5.3 Satellite dependence

The XS10 GRBs sample is made out by collecting the data available in the literature so that the final catalog is not homogenous at all. In particular, the satellite used to get both the spectral and afterglows data change from one case to another. In order to investigate whether this could have any impact on the HD, we consider again the deviations from the fiducial CDM model using only the 80 GRBs detected with the Swift satellite (excluding outliers). We get :

 ⟨Δμ⟩=−0.29  ,  (Δμ)med=0.03  ,  (Δμ)rms=1.33  ,
 ⟨Δμ⟩=−0.24  ,  (Δμ)med=0.11  ,  (Δμ)rms=1.32  ,
 ⟨Δμ⟩=−0.40