© ESO, 2017

1. Introduction

Starting at the end of the 1990s, observations of high-redshift supernovae of type Ia (SNIa) revealed that the Universe is now expanding at an accelerated rate (see e.g. Perlmutter et al. 1998, 1999; Riess et al. 1998, 2007; Schmidt et al. 1998; Astier et al. 2006; Amanullah et al. 2010). This surprising result has been independently confirmed by statistical analyses of observations of small-scale temperature anisotropies of the cosmic microwave background (CMB) radiation Spergel et al. (2007), Planck Collaboration XVI (2014), Planck Collaboration XIII (2016). It is usually assumed that the observed accelerated expansion is caused by the so-called dark energy, a cosmic medium with unusual properties. The pressure of dark energy p_de is negative and is related to the positive energy density of dark energy ϵ_de by p_de = wϵ_de, where the proportionality coefficient, that is, the equation of state (EOS), w, is negative (w< −1/3). According to current estimates, about 70% of the matter energy in the Universe is in the form of dark energy, so that today dark energy is the dominant component in the Universe. The nature of dark energy is, however, not known. The models of dark energy proposed so far can be divided into at least three groups: a) a non-zero cosmological constant (see for instance Carroll 2001), in this case w = −1, b) a potential energy of some not yet discovered scalar field (see for instance Sahni et al. 2003; Alam et al. 2003), or c) effects connected with the inhomogeneous distribution of matter and averaging procedures (see for instance Clarkson & Maartens 2010). In the last two cases, in general, w ≠ −1 and is not constant, but depends on the redshift z. To test whether and how w changes with redshift, it is necessary to use more distant objects. It is commonly argued that since the dark energy density term becomes subdominant (with respect to the dark matter) at z ≳ 0.5Riess et al. (2004), it is not important to study its EOS at earlier epochs. However, this argument is unreliable: even within the simplest model, the dark energy contributes nearly ≃10% of the overall cosmic energy density at z ≃ 2 and strongly alters the deceleration parameter with the cosmological constant. Moreover, for several observables the sensitivity to the dark energy equation of state increases at high redshifts. In Fig. 1 we explore this aspect following a simplified approach, considering the modulus of distance μ(z) as observable: we fixed a flat Λ cold dark matter (ΛCDM) fiducial cosmological model, constructing the corresponding μ_fid(z,θ), and plot the percentage error on the distance modulus with respect to different corresponding functions evaluated in the framework of a flat CPL quintessence model Chevallier & Polarski (2001), Linder (2003). We selected Ω_m = 0.3 and h = 0.7 and varied the dark EOS parameters w₀ and w₁. It is worth noting that a higher sensitivity is reached only for z ≳ 2. Therefore, investigating the cosmic expansion of the Universe also beyond these redshifts remains a fundamental probe of dark energy. In this scenario, new possibilities opened up when gamma-ray bursts (GRBs) were discovered at higher redshifts. The current record is at z = 9.4 (Tanvir et al. 2009; Salvaterra et al. 2009; Cucchiara 2011). It is worth noting that the photometric redshift on GRB 090429B is quite high, especially on the low side; GRB 090423, for which a spectroscopic redshift is available, is much better determined. GRBs span a redshift range better suited for probing dark energy than the SNIa range, as shown in Fig. 2, where we compare the distribution in redshift of our sample of 162 long GRBs/XRFs with the Union 2.1 SNIa dataset.

Fig. 1 Dependence of the distance modulus on EOS of dark energy. Upper panel: distance modulus μ(z) for different values of the EOS parameters. The black line represents the standard flat ΛCDM model with Ωm = 0.3, h = 0.7, w0 = −1, w1 = 0. The other curves correspond to different values of the CPL parameters w0 and w1. Bottom panel: z dependence of the percentage error in the distance modulus between the fiducial ΛCDM model and different flat CPL models, as in the left panel.

Fig. 2 Redshift distribution of our sample of 162 long GRBs/XRFs and the Union 2.1 SNIa dataset.

Gamma-ray bursts are enigmatic objects, however. First of all, the mechanism that is responsible for releasing the high amounts of energy that a typical GRB emits is not yet completely known, and only some aspects of the progenitor models are well established, in particular that long GRBs are produced by core-collapse events, see for instance Meszaros (2006) and Vedrenne & Atteia (2009). Despite these difficulties, GRBs are promising objects for studying the expansion rate of the Universe at high redshifts. One of the most important aspects of the observational property of long GRBs is that they show several correlations between spectral and intensity properties (luminosity, radiated energy). Here we consider the correlation between the observed photon energy of the peak spectral flux, E_p,i, which corresponds to the peak in the νF_ν spectra, and the isotropic equivalent radiated energy E_iso (e.g., Amati et al. 2002; Amati 2006), (1)where a and b are constants, and E_p,i is the spectral peak energy in the GRB cosmological rest-frame, which can be derived from the observer frame quantity, E_p, by E_p,i = E_p(1 + z). This correlation not only provides constraints for the model of the prompt emission, but also naturally suggests that GRBs can be used as distance indicators. The isotropic equivalent energy E_iso can be calculated from the bolometric fluence S_bolo as (2)where d_L is the luminosity distance and cp denotes the set of parameters that specify the background cosmological model. It is clear that to be able to use GRBs as distance indicators, it is necessary to consistently calibrate this correlation. Unfortunately, owing to the lack of GRBs at very low redshifts, the calibration of GRBs is more difficult than that of SNIa, and several calibration procedures have been proposed so far (see for instance Cardone et al. 2008; Demianski et al. 2011, 2012; Gao et al. 2012; Postnikov et al. 2014; Dai et al. 2004; Liu & Wei 2015; Lin et al. 2015). We here apply a local regression technique to determine the correlation parameters a and b, using the recently updated SNIa sample and to construct a new calibrated GRB Hubble diagram that can be used for cosmological investigations. We then use this calibrated GRB Hubble diagram to investigate the cosmological parameters through the Markov chain Monte Carlo technique (MCMC), which simultaneously computes the full posterior probability density functions of all the parameters. The structure of the paper is as follows. In Sect. 2 we describe the methods used to fit the E_p,i–E_iso correlation and construct the calibrated GRB Hubble diagram. In Sect. 3 we present our cosmological constraints and investigate the possible redshift dependence and Malmquist-like bias effects. In Sect. 4, as an additional self-consistency check, we apply the Bayesian method for the non-calibrated E_p,i–E_iso correlation and simultaneously extract the correlation coefficients and the cosmological parameters of the model. Section 5 is devoted to the discussion of our main results and conclusions.

2. Standardizing GRBs and constructing the Hubble diagram

The GRBs νF_ν spectra are well modeled by a phenomenological smoothly broken power law, characterized by a low-energy index, α, a high-energy index, β, and a break energy E₀. They show a peak corresponding to a specific (and observable) value of the photon energy E_p = E₀(2 + α). For GRBs with measured spectrum and redshift it is possible to evaluate the intrinsic peak energy, E_p,i = E_p(1 + z) and the isotropic equivalent radiated energy

where N(E) is the Band function: (3)E_iso and E_p,i span several orders of magnitude (therefore GRBs cannot be considered standard candles), and show distributions approximated by Gaussians plus a tail at low energies. A strong correlation between these two quantities was initially discovered in a small sample ofBeppo SAXGRBs with known redshifts (Amati et al. 2002) and confirmed afterward by HETE-2 and Swift observations Lamb et al. (2005), Amati (2006). Several analyses of the E_p,i-E_iso plane of GRBs showed that different classes of GRBs exhibit different behaviors, and while normal long GRBs and X-ray flashes (XRF, i.e., particularly soft bursts) follow this correlation, with the exception of the two peculiar subenergetic GRBs 980425 and 031203 see, e.g., Amati et al. (2007), for a discussion on possible explanations, short GRBs do not. These facts may depend on the different emission mechanisms and/or geometry involved in different classes of GRBs and makes this correlation a useful tool to distinguish between them (Amati 2006; Antonelli et al. 2009). Although it was the first strong correlation discovered for the GRB observables, until recent years, the E_p,i–E_iso correlation was never used for cosmology because it exhibits a significant dispersion around the best-fit power law: the residuals follow a Gaussian with a value of σ_log Ep,i ≃ 0.2. This type of additional Poissonian scatter is typically quantified by performing a maximum likelihood analysis that takes the variance and the errors on dependent and independent variables into account. By measuring E_p,i in keV and E_iso in 10⁵² erg, this method gives σ_log Ep,i = 0.19 ± 0.02 (Amati et al. 2008, 2009). However, the recent increase in the efficiency of GRB discoveries combined with the fact that the E_p,i–E_iso correlation requires only two parameters that are directly inferred from observations (this minimizes the effects of systematics and increases the number of GRBs that can be used) makes this correlation an interesting tool for cosmology. Despite the very large number of bursts consistent with this correlation, its physical origin is still debated. Some authors claimed that it is strongly affected by instrumental selection effects (Nakar & Piran 2005; Band & Preece 2005; Butler et al. 2007; Shahmoradi & Nemiroff 2011). However, many other studies found that such instrumental selection biases, even if they may affect the sample, cannot be responsible for the existence of the spectral-energy correlations (Ghirlanda et al. 2008; Nava et al. 2012). Moreover, a recent time-resolved spectral analysis of GRBs that were detected by the BeppoSAX and Fermi satellite showed that E_p,i correlates with the luminosity (e.g., Ghirlanda et al. 2010; Frontera et al. 2012) and radiated energy (e.g., Basak et al. 2013) also during the temporal evolution of the bursts (and the correlation between the spectral peak energy and the evolving flux has been pointed out by, e.g., Golenetskii et al. 1983 based on Konus-WIND data). The time-resolved E_p,i-luminosity and E_p,i–E_iso correlations within individual GRBs and that their average slope is consistent with that of correlations defined by the time-averaged spectral properties of different bursts strongly supports the physical origin of the E_p,i–E_iso correlation. It is very difficult to explain them as a selection or instrumental effect (e.g., Dichiara & Amati 2013; Basak et al. 2013), and the predominant emission mechanism in GRB prompt emission produces a correlation between the spectral peak energy and intensity (which can be characterized either as luminosity, peak luminosity, or radiated energy). In addition to its existence and slope, an important property of the E_p,i–E_iso correlation is its dispersion. As shown by several works that were based on the so-called jet-breaks, in the optical afterglow light curves of some GRBs (e.g., Ghirlanda et al. 2004, 2006), 50% of the extra-Poissonian scatter of the correlation is sometimes probably a result of the distribution of jet opening angles. However, these estimates of jet opening angles are still unconfirmed and model-dependent and can be made only for a small number of GRBs. Other factors that may contribute to the dispersion of the E_p,i–E_iso and other E_p,i-intensity correlations include the jet structure, viewing angles, detectors sensitivity, and energy band (but see Amati et al. 2009), the diversity of shock microphysics parameters, and the magnetization within the emitting ejecta. At the current observational and theoretical status, it is very difficult to quantify these single contributions, which seem to act randomly in scattering the data around the best fit power-law (e.g., Ghirlanda et al. 2008; Amati et al. 2008). Importantly, it has been shown (e.g., Amati et al. 2009) that a small fraction (5–10%) of the scatter depends on the cosmological model and parameters assumed for computing E_iso, which makes this correlation a potential tool for cosmology.

In this section we show how the E_p,i–E_iso correlation can be calibrated to standardize long GRBs and to build a GRB Hubble diagram, which we use to investigate different cosmological models at very high redshifts (see also Amati & Della Valle 2013; Wang et al. 2015, 2016; Lin et al. 2015, 2016b).

2.1. GRB data

We used a sample of 162 long GRBs/XRFs as of September 2013 taken from the updated compilation of spectral and intensity parameters of GRBs by Sawant & Amati (2016). The redshift distribution of this sample covers a broad range, 0.03 ≤ z ≤ 9.3, which means that it extends far beyond the SNIa range (z ≤ 1.7). These data are of long GRBs/XRFs that are characterized by firm measurements of redshifts and the restframe peak energy E_p,i. The main contributions come from the joint detections by Swift/BAT and Fermi/GBM or Konus-WIND, except for the small fraction of events for which Swift/BAT can directly provide E_p,i when it is in the (15−150) keV interval. For events detected by more than one of these detectors, the uncertainties on the E_p,i and E_iso take the measurements and uncertainties provided by each individual detector into account. In Table 4 we report for each GRB the redshift, the restframe spectral peak energy, E_p,i, and the isotropic-equivalent radiated energy, E_iso. As detailed in (Sawant & Amati, in prep.), the criterion behind selecting the observations from a particular mission is based on the conditions summarized below.

• 1.

  We preferred observations with exposure times of at leasttwo-thirds of the whole event duration. Hence Konus-WINDand Fermi/GBM were chosen whenever available. For Konus-WIND, the observations were reported in the official literature(see e.g. Ulanov et al. 2005),and also in GCN archives when needed. For Fermi/GBM, theobservations were taken from Gruber et al.(2012), from several other papers (e.g., Ghirlanda etal. 2004, 2005; Friedman& Bloom 2005, etc.), and from GCNs.

• 2.

  The Swift/BAT observations were chosen when no other preferred mission (Konus-WIND, Fermi/GBM) was able to provide spectral parameters and the value of E_p,obs was within the energy band of this instrument. In particular, the E_p,i values derived from BAT spectral analysis were conservatively taken from the results reported by the BAT team (Sakamoto et al. 2008a,b). Other BAT E_p,i values reported in the literature were not considered because they were not confirmed by Sakamoto et al. (2008a,b), by a refined analysis (see e.g. Cabrera et al. 2007), or because they were based on speculative methods Butler et al. (2007).

When available, the Band model Band et al. (1993) was considered as the cut-off power law, which sometimes overestimates the value of E_p,i. Finally, the error on any value was assumed to be not less than 10% to account for the instrumental capabilities and calibration uncertainties.

2.2. Cosmologically independent calibration: local regression of SNIa

The lack of nearby GRBs creates the so-called circularity problem: GRBs can be used as cosmological tools through the E_p,i–E_iso correlation, which is based on the cosmological model assumed for the computation of E_iso, however. In principle, this problem could be solved in several ways: it is possible, for instance, to simultaneously constrain the calibration parameters (a,b,σ_int) ∈ G and the set of cosmological parameters θ ∈ H by considering a likelihood function defined in the space G ⊗ H, which allows simultaneously fitting the parameters (e.g., Diaferio et al. 2011). This procedure is implemented in Sect. 5 and compared with the local regression technique. Alternatively, it has been proposed that the problem might be avoided by considering a sufficiently large number of GRBs within a small redshift bin centered on any zGhirlanda et al. (2006). However, this method, even if quite promising for the future, is currently unrealistic because of the limited number of GRBs with measured redshifts. In this section we implement a procedure for calibrating the E_p,i–E_iso relation in a way independent of the cosmological model by applying a local regression technique to estimate the distance modulus μ(z) from the recently updated SNIa sample, called Union2.1 (see also Kodama et al. 2008; Liang et al. 2008). Originally implemented by Cleveland & Devlin (1988), this technique combines the simplicity of linear regression with the flexibility of nonlinear regression to localized subsets of the data, and reconstructs a function describing their behavior in the neighborhood of any point z₀. A low-degree polynomial is fitted to a subsample containing a neighborhood of z₀, by using weighted least-squares with a weight function that quickly decreases with the distance from z₀. The local regression procedure can be schematically sketched as follows:

• 1.

  set the redshift z where μ(z) has to be reconstructed;

• 2.

  sort the SNIa Union2.1 sample by increasing value of | z−z_i | and select the first , where α is a user-selected value and the total number of SNIa;

• 3.

  apply the weight function (4)where u = | z−z_i | /Δ and Δ the highest value of the | z−z_i | over the previously selected subset;

• 4.

  fit a first-order polynomial to the data previously selected and weighted, and use the zeroth-order term as the best-fit value of the modulus of distance μ(z);

• 5.

  evaluate the error σ_μ as the root mean square of the weighted residuals with respect to the best-fit value.

It is worth stressing that both the choice of the weight function and the order of the fitting polynomial are somewhat arbitrary. Similarly, the value of n to be used must not be too low to make up a statistically valuable sample, but also not too high to prevent the use of a low-order polynomial. To check our local regression routine, we simulated a large catalog with the same redshift and error distribution as the Union2.1 survey. We adopted a quintessence cosmological model with a constant EOS, w, and fixed values of the cosmological parameters (Ω_M,w,h). For each redshift value in the Union2.1 sample, we extracted the corresponding modulus of distance from a Gaussian distribution centered on the theoretically predicted value and with the standard deviation σ = 0.15. To this value, we added the value of the error, corresponding to the same relative uncertainty as the data in the Union sample. This simulated catalog was used as input to the local regression routine, and the reconstructed μ(z) values were compared to the input ones.

Fig. 3 Best-fit curves for the Ep,i–Eiso correlation relation superimposed on the data. The solid and dashed lines refer to the results obtained with the maximum likelihood (Reichart likelihood) and weighted χ2 estimator, respectively.

Defining the percentage deviation with μ_rec. and μ_w the local regression estimate and the input values, respectively, and averaging over 500 simulations, we found that the choice α = 0.02 gives (δμ/μ)_rms ≃ 0.3% with | ϵ | ≤ 1% independently on the redshift z. This result implies that the local regression method allows correctly reconstructing the underlying distance modulus regardless of redshift z< 1.4 from the Union SNIa sample. We also tested this results in different cosmological backgrounds by adopting an evolving EOS and averaging over five realizations of the mock catalog. With this efficient way of estimating μ(z) at redshift z in a model-independent way, we can now fit the E_p,i–E_iso correlation relation, using the reconstructed μ(z). We considered only GRBs with z ≤ 1.414 to cover the same redshift range as is spanned by the SNIa data. To standardize the E_p,i–E_iso relation as expressed by Eq. (1), we need to fit a data array { x_i,y_i } with uncertainties { σ_x,i,σ_y,i }, to a straight line (5)and determine the parameters (a,b). We expect a certain amount of intrinsic additional Poissonian scatter, σ_int, around the best-fit line that has to be taken into account and determined together with (a,b) by the fitting procedure. We used a likelihood, implemented by Reichart et al. (2001), that solved this problem of fitting data that are affected by extrinsic scatter in addition to the intrinsic uncertainties along both axes: (6)where the sum is over the objects in the sample. We note that this maximization was performed in the two-parameter space (a,σ_int) since b may be calculated analytically by solving the equation , we obtained (7)To quantitatively estimate the goodness of this fit, we used the median and root mean square of the best-fit residuals, defined as δ = y_obs−y_fit. To quantify the uncertainties of some fit parameter p_i, we evaluated the marginalized likelihood ℒ_i(p_i) by integrating over the other parameters. The median value for the parameter p_i was then found by solving (8)The 68% (95%) confidence range (p_i,l,p_i,h) was then found by solving e.g., D’Agostini (2005)(9)with ε = 0.68 and ε = 0.95 for the 68% and 95% confidence level. The maximum likelihood values of a and σ_int are and . In Fig. 3 we show the correlation log E_p,i−log E_iso. The solid line is the best fit obtained using the Reichart likelihood, and the dashed line is the best fit obtained by the weighted χ² method. The marginalized likelihood functions are shown in Fig. 4.

Fig. 4 Marginalized likelihood functions: the likelihood function ℒa is obtained by marginalizing over σint; and the likelihood function, ℒσint, is obtained by marginalizing over a.

As noted above, b can be evaluated analytically. We obtained b = 52.53 ± 0.02. It is worth noting that in the literature results for the inverse correlation are commonly reported, that is, the correlation E_iso–E_p,i: using the local regression technique and the Reichart likelihood, we also obtained this inverse calibration. We obtained , . In Fig. 5 we plot the best-fit curves for the E_iso–E_p,i correlation relation superimposed on the data.

Fig. 5 Best-fit curves for the Eiso–Ep,i correlation relation superimposed on the data. The solid and dashed lines refer to the results obtained with the maximum likelihood (Reichart likelihood) and weighted χ2 estimator, respectively.

Fig. 6 Calibrated GRB Hubble diagram (black filled diamond) up to very high values of redshift, as constructed by applying the local regression technique with the Reichart likelihood: we show overplotted the behavior of the theoretical distance modulus μ(z) = 25 + 5log dL(z) corresponding to the favored fitted CPL model (red line), with w0 = −0.29, w1 = −0.12, h = 0.74, Ωm = 0.24, and the standard ΛCDM model (gray line) with Ωm = 0.33, and h = 0.74. They are defined by Eqs. (12) and (21). The blue filled points at lower redshift are the SNIa data.

2.2.1. Constructing the Hubble diagram

After fitting the correlation and estimating its parameters, we used them to construct the GRB Hubble diagram. We recall that the luminosity distance of a GRB with the redshift z may be computed as (10)The uncertainty of d_L(z) was then estimated through the propagation of the measurement errors on the involved quantities. In particular, recalling that our correlation relation can be written as a linear relation, as in Eq. (5), where y = E_iso is the distance dependent quantity, while x is not, the error on the distance dependent quantity y was estimated as (11)where σ_b is properly evaluated through the Eq. (27), which implicitly defines b as a function of a and σ_int, and is then added in quadrature to the uncertainties of the other terms entering Eq. (10) to obtain the total uncertainty. The distance modulus μ(z) is easily obtained from its standard definition (12)with its uncertainty obtained again by error propagation: (13)We finally estimated the distance modulus for each ith GRB in the sample at redshift z_i to build the Hubble diagram plotted in Fig. 6. We refer to this data set as the calibrated GRB Hubble diagram below since to compute the distances Hubble diagram we relied on the calibration based on the SNIa Hubble diagram. The derived distance moduli are divided into two subsets, listed in Tables 4 and 5. Table 4 lists GRBs with z ≤ 1.46, the same redshift range as for known SNIa, and Table 5 lists GRBs with z ≥ 1.47. In Fig. 7 we finally compare the GRB Hubble diagram (black points) with the SNIa Hubble diagram (blue points) and with BAO data (red points).

Fig. 7 GRB Hubble diagram (black points) is compared with the SNIa Hubble diagram (blue points) and with BAO data (red points).

3. Cosmological constraints derived from the calibrated GRB Hubble diagram

Here we illustrate the possibilities of using the GRB Hubble diagram to constrain the cosmological models and investigate the dark energy EOS. Within the standard homogeneous and isotropic cosmology, the dark energy appears in the Friedman equations: where a is the scale factor, H = ȧ/a the Hubble parameter, ρ_M the density of matter, ρ_de the density of dark energy, p_de its pressure, and the dot denotes the time derivative. The continuity equation for any component of the cosmological fluid is (16)where the energy density is ρ_i, the pressure p_i, and the EOS of the ith component is defined by . The standard non-relativistic matter has w = 0, and the cosmological constant has w = −1. The dark energy EOS and other constituents of the Universe determine the Hubble function H(z) and any derivations of it that are needed to obtain the observable quantities. When only matter and dark energy are present, the Hubble function is given by (17)where , and w(z,θ) describes the dark energy EOS, characterized by n parameters θ = (θ₁,θ₂,...,θ_n). We limit our analysis below to the CPL parametrization, (18)where w₀ and w₁ are constant parameters and they represent the w(z) present value and its overall time evolution, respectively Chevallier & Polarski (2001), Linder (2003). If we introduce the dimensionless Hubble parameter E(z,θ): (19)we can define the luminosity distance and the modulus of distance. Actually the luminosity distance id given by: (20)where , d_M(z,θ) is the transverse comoving distance and it is defined as (21)being d_C(z,θ) the comoving distance: (22)therefore we can define the modulus of distance μ_th(z,θ): (23)The standard ΛCDM model corresponds to w₀ = −1, w₁ = 0.

To constrain the parameters specifying different cosmological models, we maximized the likelihood function ℒ(θ) ∝ exp [−χ²(θ)/2], where θ indicates the set of cosmological parameters and the χ²(θ) was defined as usual by (24)Here, μ_obs and μ_th are the observed and theoretically predicted values of the distance modulus. The parameter space is efficiently sampled by using the MCMC method, thus running five parallel chains and using the Gelman-Rubin test to check the convergence Gelman & Rubin (1992). The confidence levels are estimated from the likelihood quantiles. We recall that we performed the analysis assuming a non-zero spatial curvature (not flat ΛCDM), and only in this case did we take w = −1. To alleviate the strong degeneracy of the curvature parameters and any EOS parameters, we added a Gaussian prior on Ω_k, centered on the value provided by the Planck Collaboration XVI (2014), , and with a dispersion ten times of , where . We investigated the CPL parameters assuming a flat universe. In a forthcoming paper we will present a full cosmological analysis using the high-redshift GRB Hubble diagram to test different cosmological models, where several parametrizations of the dark energy EOS will be used and also different dark energy scenarios, for instance the scalar field quintessence. In Table 1 we summarize the results of our analysis. There are indications that the dark energy EOS is evolving.

Fig. 8 Regions of confidence for the marginalized likelihood function ℒ(w0,w1). On the axes we also plot the box-and-whisker diagrams for the w0 and w1 parameters: the bottom and top of the diagrams are, as commonly done in the box-plots, the 25th and 75th percentile (the lower and upper quartiles, respectively), and the band near the middle of the box is the 50th percentile (the median). It is evident that the ΛCDM model, corresponding to w0 = −1 and w1 = 0 is not favored.

The joint probability for different sets of parameters that characterize the CPL EOS is shown in Fig. 8.

4. E_iso and E_p,i correlation

Since we here used the calibrated GRB Hubble diagram to perform the cosmological investigation described above, we discuss some arguments about the reliability of using the E_iso−E_p,i relation for cosmological tasks. For instance, the calibration method we implemented so far relies on the underlying assumption that the calibration parameters do not evolve with redshift. It is worth noting that the problem of the redshift dependence of the GRBs correlations is widely debated in the literature, with different conclusions. This problem is intimately connected to the problem of the influence of possible selection effects or biases on the observed correlations, see, for instance, Li et al. (2008), Dainotti et al. (2013a,b), Liu & Wei (2015), Wang et al. (2016), Lin et al. (2016a,b). Any answer to these fundamental questions is far from being settled until more data with known redshifts are available. In this section, however, we revisit this question from an observational point of view: we test the validity of this commonly adopted working hypothesis and search for any evidence of such a redshift dependence. We also investigate possible effects on the GRB Hubble diagram. As a first simplified approach we considered two subsamples with a comparable number of bursts, divided according to redshift: a lower redshift sample of 97 bursts with z ≤ 2, and a higher redshift sample of 67 burst with z > 2. We estimated the cosmological parameters for flat ΛCDM and CPL dark energy EOS cosmological models, considering the two subsamples separately, and compared the results. Even when the bursts belonging to these two samples experience different environment conditions, we did not find significant indication that a spurious z-evolution of the slope affects the cosmological fit: the values of the cosmological parameters derived from these two samples are statistically consistent within 1σ, as shown in Table 2. This result indirectly shows that redshift evolution dependence, even if it exists, does not undermine the reliability of the GRBs as probes of the cosmological expansion. Moreover, to investigate this redshift dependence in more detail, that is, how it affects E_iso and/or E_p,i (which we generically denote by y), we used two different approaches. First, we evaluated the Spearman rank correlation coefficient, C(z,y), taking into account that since our sample is not too large, a few points could dominate the final value of the rank correlation, which would introduce a bias. We applied a jacknife resampling method by evaluating C(z,y) for N−1 samples obtained by excluding one GRB at a time, and we adopted the mean value and the standard deviation to estimate C(z,y). We found C(z,E_p,i) = 0.299 ± 0.004, and C(z,E_iso) = 0.278 ± 0.004, which indicates a moderate evolution of the correlation. The Spearman rank correlation coefficient, however, does not include the errors of E_iso and E_p,i, so that we also implemented an alternative, and completely different, approach to determine whether and how strongly these variables are correlated with redshift. We assumed that the evolutionary functions can be parametrized by a simple power-law functions: g_iso(z) = ⁽1 + z^)k_iso and g_p(z) = ⁽1 + z^)k_p (see for instanceDainotti et al. 2013a), so that /g_iso(z) and /g_p(z) are the de-evolved quantities. In this case, the effective E_iso−E_p,i correlation can be written as a 3D correlation:(25)Calibrating this 3D relation means determining the coefficients (a, b, k_iso, and k_p) plus the intrinsic scatter σ_int. Low values for k_iso and k_p would indicate a lack of evolution, or at least negligible evolutionary effects. To fit these coefficients, we constructed a 3D Reichart likelihood, but we consider no error on the redshift of each GRB: (26)where α = ak_p. As in the 2D case, we maximized our likelihood in the space (a, k_iso, and α) since b was evaluated analytically by solving the equation , we obtain(27)We also used the MCMC method and ran five parallel chains and the Gelman-Rubin convergence test, as previously explained. We finally studied the median and 68% confidence range of k_iso and α to test whether the correlation evolves and noted that a null value for these parameters is strong evidence for a lack of any evolution. We found that , k_iso = −0.04 ± 0.1; α = −0.02 ± 0.2; σ_int = 0.35 ± 0.03, so that . We can safely conclude that the E_iso and E_p,i correlation shows, at this stage, weak indications of evolution. In Fig. 9 we plot both the de-evolved and evolved/original correlation, and the de-evolved and evolved/original Hubble diagram: we do not see any signs of evolution.

5. Fully Bayesian analysis

Fig. 9 Evolutionary effects on the Eiso – Ep,i correlation and the Hubble flow. Upper panel: de-evolved (red points) and evolved/original correlation (green poins); there is no discernable evolution. Bottom panel: de-evolved (red points) and evolved/original (green points) Hubble diagram; evolutionary effects of the correlation do not affect the GRB Hubble diagram.

In this section we simultaneously constrain the calibration parameters (a,b,σ_int) and the set of cosmological parameters by considering the same likelihood function as in Eq. (7). Our task is to determine the multidimensional probability distribution function (PDF) of the parameters { a,b,σ_int,p }, where p is the -dimensional vector of the cosmological parameters. The Amati correlation can be written in the form (28)We note that the best-fit zero point b can be analytically expressed as a function of (a,σ_int) when the cosmological parameter p_c is specified. A more computationally efficient strategy is to let b free and add it to the list of quantities to be determined, which thus sums up to . To efficiently sample the dimensional parameter space, we used the MCMC method and ran five parallel chains and used the Gelman-Rubin convergence test, as described in the previous section. Since we are mainly interested in the calibration problem and not in constraining the cosmological parameters (Ω_M,Ω_Λ,h), we considered here only the particular case of a flat ΛCDM model, as already done in the previous analysis. The analysis determines, at the same time, the cosmological parameters and the correlation coefficients, which are listed in Table 3.

Although the calibration procedure is different (since we now fit for the cosmological parameters as well), it is nevertheless worth comparing this determination of (a,b,σ_int) with the one obtained in the previous analysis based on the SNIa sample. It is evident that although the median values change, the 95% confidence levels are in full agreement so that we cannot find any statistically significant difference. Figure 10 shows the marginalized probability distribution function of the correlation coefficients of the relation E_iso−E_p,i.

As far as the cosmological parameters are concerned, it turns out that , the range of confidence at 1σ is (0.13,0.54), and , the range of confidence at 1σ is (0.50,0.87). This result implies that , and that the range of confidence is (− 0.00730,0.00605), in a good agreement with results derived by using the calibrated GRBs Hubble diagram.

Also in this case we finally estimate the distance modulus for each ith GRB in our sample at redshift z_i, to build the fiducial Hubble diagram, by using Eqs. (1) and (10). It turns out that the fiducial and calibrated Hubble diagrams are fully statistically consistent, as shown in Fig. 11.

Fig. 10 Regions of confidence for the marginalized likelihood function ℒ(a,σ), obtained by marginalizing over b and the cosmological parameters. The gray regions indicate the 1σ (full zone) and 2σ (empty zone) regions of confidence. On the axes we also plot the box-and-whisker diagrams for the a and σint parameters: the bottom and top of the diagrams are the 25th and 75th percentile (the lower and upper quartiles, respectively), and the band near the middle of the box is the 50th percentile (the median).

6. Discussion and conclusions

The E_p,i–E_iso correlation is one of the most intriguing properties of the long GRBs, with significant implications for the use of GRBs as cosmological probes. Here we explored the Amati relation in a way independent of the cosmological model. Using the recently updated data set of 162 high-redshift GRBs, we applied a local regression technique to estimate the distance modulus using the recent Union SNIa sample (Union2.1). The derived calibration parameters are statistically fully consistent with the results of our previous work Demianski et al. (2011). Moreover, we tested the validity of the commonly adopted working hypothesis that the GRB Hubble diagram is slightly affected by redshift dependence of the E_p,i–E_iso correlation. As a first qualitative and simplified approach we considered a lower redshift sample of 97 bursts with z ≤ 2, and a higher redshift sample of 67 burst with z > 2. We estimated the cosmological parameters for flat ΛCDM and CPL dark energy EOS cosmological models, considering the two subsamples separately, and compared the results. Even when the bursts belonging to these two samples had different environment conditions, we found no significant indications that a spurious z-evolution of the slope affected the cosmological fit. Moreover, to quantify this redshift dependence, we used two different approaches. First, we evaluated the Spearman rank correlation coefficient, C(z,y) by applying a jacknife resampling method by evaluating C(z,y) for N−1 samples obtained by excluding one GRB at a time, and we adopted the mean value and the standard deviation to estimate C(z,y). C(z,E_p,i) = 0.299 ± 0.004, and C(z,E_iso) = 0.278 ± 0.004, which indicates a negligible evolution of the correlation. Moreover, we also implemented an alternative method, assuming that the redshift evolution can be parametrized by simple power-law functions: g_iso(z) = ⁽1 + z^)k_iso and g_p(z) = ⁽1 + z^)k_p and that the correlation holds for the de-evolved quantities and . In this case, we rewrote an effective 3D E_iso−E_p,i correlation, with included the evolutionary terms. Since we were interested in the implications of possible evolutionary effects of the GRB Hubble diagram and to simplify the fit, we introduced an auxiliary variable α = ak_p. A null value for k_iso and α is strong evidence for a lack of any evolution.

Fig. 11 Fiducial GRB Hubble diagram, superimposed on the calibrated Hubble diagram.

To fit the coefficient of our 3D correlation we constructed a 3D Reichart likelihood and used a MCMC method running five parallel chains and using the Gelman-Rubin convergence test. The E_iso and E_p,i correlation shows at this stage only weak indications of evolution. The derived calibration parameters were used to construct an updated calibrated GRB Hubble diagram, which we adopted as a tool to constrain the cosmological models and to investigate the dark energy EOS. In particular, we searched for any indications that w(z) ≠ −1, which reflects the possibility of a deviation from the ΛCDM cosmological model. To accomplish this task, we focused on the CPL parametrization as an example. To efficiently sample the cosmological parameter space, we again used a MCMC method. At 1σ level we found indications for a time evolution of the EOS in the considered parametrization; we conclude that the ΛCDM model is not favored even though it is not ruled out by these observations so far. Moreover, for w = −1 we also performed our analysis assuming a non-zero spatial curvature, adding a Gaussian prior on Ω_k centered on the value provided by the Planck data Planck Collaboration XVI (2014). The GRB Hubble diagram alone provides for the range of confidence at 1σ(− 0.007,0.0064). Finally, to investigate the reliability of the E_p,i–E_iso correlation in greater detail, we also used a different method to simultaneously extract the correlation coefficients and the cosmological parameters of the model from the observed quantities. To illustrate this method we assumed here, by way of an example, only the particular case of a non-flat ΛCDM model, as already done in the previous analysis. This analysis simultaneously determines the cosmological parameters and the correlation coefficients. Although the calibration procedure is different (since we now fit for the cosmological parameters as well), it is nevertheless worth comparing this determination of (a,b,σ_int) with the one obtained in the previous analysis based on the SNIa sample. It is evident that although the median values change, the 95% confidence levels are in full agreement so that we cannot find any statistically significant difference. This means that the high-redshift GRBs can be used as cosmological probes, mainly in a redshift region, , which is unexplored by SNIa and BAO samples, and that both the calibration technique based on a local regression with SNIa and a fully Bayesian approach are reliable. As a final remark we note that our results for the cosmological parameters are statistically consistent with those previously obtained by other teams, as shown in Table 6, where we display some of the most recent measurements of cosmological parameters obtained with GRBs, even if following a slightly different approach. The main peculiarity in our analysis consists of the procedure used to build up the dataset, consisting of 162 objects, as discussed in Sect. 2.1 (see also Sawant & Amati, in prep.), even more than the specific statistical analysis, which follows a Bayesian approach through the implementation of the MCMC. Because the reliability of GRBs as distance indicators has not yet been clearly proved and because it constitutes an independent topic in the field of GRB research, it is not meaningless to investigate the reliability of the E_p,i–E_iso correlation whenever new and improved datasets are available. In the near future we intend to enhance the cosmological analysis by using the high-redshift GRB Hubble diagram to test different cosmological models, where several parametrizations of the dark energy EOS are used, but a cosmographic approach will also be implemented, which will update the analysis performed in Demianski et al. (2012) to check whether with this new updated dataset the estimates of the jerk and the dark energy parameters will confirm deviations from ΛCDM cosmological model, as has been indicated in our previous analysis.

Acknowledgments

M.D. is grateful to the INFN for financial support through the Fondi FAI GrIV. E.P. acknowledges the support of INFN Sez. di Napoli (Iniziative Specifica QGSKY and TEONGRAV). L.A. acknowledges support by the Italian Ministry for Education, University and Research through PRIN MIUR 2009 project on “Gamma ray bursts: from progenitors to the physics of the prompt emission process”. (Prot. 2009 ERC3HT.)

References

All Tables

All Figures

Fig. 1 Dependence of the distance modulus on EOS of dark energy. Upper panel: distance modulus μ(z) for different values of the EOS parameters. The black line represents the standard flat ΛCDM model with Ωm = 0.3, h = 0.7, w0 = −1, w1 = 0. The other curves correspond to different values of the CPL parameters w0 and w1. Bottom panel: z dependence of the percentage error in the distance modulus between the fiducial ΛCDM model and different flat CPL models, as in the left panel.

In the text

	Fig. 2 Redshift distribution of our sample of 162 long GRBs/XRFs and the Union 2.1 SNIa dataset.
In the text

	Fig. 3 Best-fit curves for the Ep,i–Eiso correlation relation superimposed on the data. The solid and dashed lines refer to the results obtained with the maximum likelihood (Reichart likelihood) and weighted χ2 estimator, respectively.
In the text

	Fig. 4 Marginalized likelihood functions: the likelihood function ℒa is obtained by marginalizing over σint; and the likelihood function, ℒσint, is obtained by marginalizing over a.
In the text

	Fig. 5 Best-fit curves for the Eiso–Ep,i correlation relation superimposed on the data. The solid and dashed lines refer to the results obtained with the maximum likelihood (Reichart likelihood) and weighted χ2 estimator, respectively.
In the text

Fig. 6 Calibrated GRB Hubble diagram (black filled diamond) up to very high values of redshift, as constructed by applying the local regression technique with the Reichart likelihood: we show overplotted the behavior of the theoretical distance modulus μ(z) = 25 + 5log dL(z) corresponding to the favored fitted CPL model (red line), with w0 = −0.29, w1 = −0.12, h = 0.74, Ωm = 0.24, and the standard ΛCDM model (gray line) with Ωm = 0.33, and h = 0.74. They are defined by Eqs. (12) and (21). The blue filled points at lower redshift are the SNIa data.

In the text

	Fig. 7 GRB Hubble diagram (black points) is compared with the SNIa Hubble diagram (blue points) and with BAO data (red points).
In the text

Fig. 8 Regions of confidence for the marginalized likelihood function ℒ(w0,w1). On the axes we also plot the box-and-whisker diagrams for the w0 and w1 parameters: the bottom and top of the diagrams are, as commonly done in the box-plots, the 25th and 75th percentile (the lower and upper quartiles, respectively), and the band near the middle of the box is the 50th percentile (the median). It is evident that the ΛCDM model, corresponding to w0 = −1 and w1 = 0 is not favored.

In the text

	Fig. 9 Evolutionary effects on the Eiso – Ep,i correlation and the Hubble flow. Upper panel: de-evolved (red points) and evolved/original correlation (green poins); there is no discernable evolution. Bottom panel: de-evolved (red points) and evolved/original (green points) Hubble diagram; evolutionary effects of the correlation do not affect the GRB Hubble diagram.
In the text

Fig. 10 Regions of confidence for the marginalized likelihood function ℒ(a,σ), obtained by marginalizing over b and the cosmological parameters. The gray regions indicate the 1σ (full zone) and 2σ (empty zone) regions of confidence. On the axes we also plot the box-and-whisker diagrams for the a and σint parameters: the bottom and top of the diagrams are the 25th and 75th percentile (the lower and upper quartiles, respectively), and the band near the middle of the box is the 50th percentile (the median).

In the text

	Fig. 11 Fiducial GRB Hubble diagram, superimposed on the calibrated Hubble diagram.
In the text