© ESO 2020

1. Introduction

Modern observations increasingly shed light on the properties of the Universe’s constituents (Riess et al. 1998; Perlmutter et al. 1999; Planck Collaboration VI 2020; Verde et al. 2019). Even though data are more and more accurate, current understanding of cosmology suggests that most of the Universe’s content is in the form of dark energy (∼68%) and dark matter (∼23%). Unveiling whether dark energy is a cosmological constant¹ and understanding dark matter’s nature² are two of today’s most complex challenges (Brax et al. 2009; Bertone & Tait 2018; Luongo & Muccino 2018). For decades, disclosing Universe dynamics and clustering properties of its constituents has lead to the search for three quantities: the current Hubble rate, the present deceleration parameter value, and the amount of visible and invisible matter (Martel et al. 1998; Bousso & Polchinski 2000; Garriga et al. 2000; Scherrer 2018; Bengaly 2020). Although data from standard indicators have progressively exceeded the redshift limit of our local Universe (z ≃ 0) measurements are usually jeopardized by postulating the model a priori. Thus, model-independent treatments that permit us to fix bounds over the dark energy and dark matter contributions have become much more relevant (Sahni & Starobinsky 2000; Sahni et al. 2008, 2014; Aviles et al. 2012; Seikel et al. 2012; Yahya et al. 2014; de la Cruz-Dombriz et al. 2016; Dunsby & Luongo 2016), as they are fascinating topics providing subtle results. The purpose of these prescriptions is to frame out all dynamics, without postulating the cosmological model a priori. Among several treatments, cosmography aims to fit the observable quantities of interest by means of the Taylor series up to a certain order (Capozziello et al. 2013; Aviles et al. 2014). Cosmographic expansions intrinsically own the idea of model-independently fitting data by fixing limits over scale factor derivatives, expanding around z ≃ 0. The elegance of cosmography clashes with current surveys that span beyond redshift intervals where cosmography turns out to be predictive (Cattoën & Visser 2008). Due to this restriction over z, the need for a high-redshift cosmography that involves redshifts outside the local Universe is at an all-time high (Cattoën & Visser 2008). While a self-consistent approach to high-redshift cosmography is far out of reach, one requires high-redshift standard indicators alongside current approaches to fit data at z >  1. A possible relevant high-redshift indicator is the presence of gamma-ray bursts (GRBs), which represent the most powerful explosions in the Universe (Meszaros 2006). Their use in cosmology would be a major improvement on currently known outcomes, since they likely originate from black hole formation, emit huge amounts of energy up to 10⁵⁴ erg, and are observable from larger distances than supernovae (SNe). Understanding if GRBs can be suitable variants for standard candles, meaning, if they can be considered as standard indicators, is an open challenge for astrophysics due to the well-known circularity problem. Moreover, matching universe kinematics with GRB data would give new insights into cosmodynamics. However, relating cosmography to GRBs is a difficult task for two reasons: (1) GRBs are not yet pure standard candles, and (2) cosmography is accurate for z ≪ 1 only, whereas GRBs are objects that severely exceed the limit z ≃ 0. To fix these issues, we must (a) find a procedure to calibrate GRBs³, (b) manipulate GRB data points and make them suitable for truncated cosmographic series (see, e.g., Wang et al. 2015). On the one hand, in order to explore the first issue, one requires a self-consistent GRB theoretical model that is intrinsically capable of connecting all GRB photometric and spectral properties. Despite this caveat being far from clarified, there are consolidate GRB correlations that try to relate photometric and spectral properties. On the other hand, to enable cosmography to be predictive, two methods are currently available: auxiliary variables and rational approximations of cosmic distances. Both of them permit one to exceed the limit z ≪ 1 for fitting procedures.

Hence, the aim of this work is twofold. First, we considered the mostly used GRB correlations with the lowest data-point spread: Amati, Ghirlanda, Yonetoku, and combo correlations. Then, we overcame the circularity problem by calibrating the above correlations through a purely model-independent technique at low redshift based on Bézier’s approximants and first presented in Amati et al. (2019). Once the correlations were calibrated, our first task was to model-independently enable GRBs to be quasi-standard candles. We thus considered three approaches to fit the Universe’s dynamics by means of kinematics. To do so, we employed the first approach purely based on genuine low-redshift cosmography. The other two options resulted from a high-redshift cosmography scenario and consist of using auxiliary variables and optimal rational approximations of the Taylor expansions inferred from Padé polynomials, respectively. We therefore worked with two hierarchical scenarios, initially involving the sets 𝒜₁ ≡ {h₀, q₀, j₀}, and 𝒜₂ ≡ {h₀, q₀, j₀, s₀} later. This was done by fixing the orders of our expansions. We performed our statistical analyses using the most recent Pantheon survey of SNeIa, GRB data sets from each of the aforementioned correlations, and the BOSS measurements for Baryonic acoustic oscillations (BAO). We performed a Monte Carlo analysis based on the widely adopted Metropolis-Hastings algorithm, with corresponding contours reported up to 2 − σ confidence levels. Although preliminary, our analyses show that the use of GRBs does not fully confirm the hypothesis that dark energy is a pure cosmological constant. The discrepancies are discussed in terms of evolving dark energy or by considering the effects of spatial curvature that have been conventionally neglected in our treatments. In this respect, an in-depth discussion on systematics and a comparison with previous works have been reported, indicating the goodness of our outcomes and the consequences on cosmography.

The paper is structured as follows. In Sect. 2, we highlight the main features concerning the GRB correlations. In Sect. 3, we propose the cosmographic method, and we split its use between low- and high-redshift cosmology, respectively. We thus introduce the auxiliary variable and rational approximant methods. In Sect. 4, we describe the Bézier calibration, overcoming the circularity problem previously discussed throughout the text. Afterward, in the same section, we provide details of our experimental analysis, describing each test performed here. In Sect. 5, we summarize our results, and we propose a theoretical description of our outcomes by means of a slightly evolving scalar field that mimicks dark energy. Finally, in Sect. 6, we state our conclusions to and perspectives on our approaches.

2. GRB correlations

Gamma-ray burst correlations are of utmost importance since they relate observable quantities among them and represent strategies to enable GRBs to be used as distance indicators (for details, please refer to the following works: Amati et al. 2002, 2008; Ghirlanda et al. 2004a; Schaefer 2007; Capozziello & Izzo 2008; Dainotti et al. 2008; Bernardini et al. 2012; Amati & Della Valle 2013; Wei et al. 2014; Izzo et al. 2015; Demianski et al. 2017a,b). Calibration procedures thus become essential to healing the circularity problem, plaguing each correlation. The circularity issue consists of the fact that GRB distances are typically obtained through photometry only, and that they depend on the cosmological parameters via H(z), which is the cosmological model that describes Universe dynamics at large scales. So, using GRBs to disclose the Universe’s expansion history but calibrating it with a model would a priori represent a tricky circularity. This issue is not easy to overcome because there is no low-redshift set of GRBs to achieve a cosmology-independent calibration. The common approach to overcoming circularity is to forcibly add further data at lower redshifts (see, e.g., Demianski et al. 2017a, and references therein).

In our work, we employed a new technique that is capable of reconstructing low-redshift correlations in a model-independent way. First, we summarize below the most suitable correlations employed in the literature, such as the ones in which the spread of each GRB point is not as large as it is in other correlations. Later, we describe how to calibrate them by means of Bézier polynomials.

2.1. The Amati correlation

The most investigated correlation in GRB literature is the E_p − E_iso or Amati correlation (Amati et al. 2002, 2008; Amati & Della Valle 2013; Demianski et al. 2017a; Dainotti & Amati 2018):

$$\begin{array}{c}\hfill log(E_{\mathrm{p}}/\mathrm{keV})=a_0+a_1[log(E_{\mathrm{iso}}/\mathrm{erg})-52],\end{array}$$(1)

where a₀ and a₁ are the two calibration constants. This correlation involves prompt emission observables: the rest-frame peak energy E_p of the νF_ν energy spectrum and the isotropic radiated energy in γ-rays E_iso = 4πd_LS_bolo(1 + z)⁻¹, in which d_L is the luminosity distance to the source, S_bolo the bolometric fluence computed from the integral of the νF_ν spectrum in the rest-frame 1−10⁴ keV energy band, and the factor (1 + z)⁻¹ transforms the observed GRB duration into the rest-frame one. The Amati correlation is characterized by an extra source of variability σ_a, due to the contribution of hidden variables that we cannot observe directly (D’Agostini 2005).

2.2. The Ghirlanda correlation

Theoretical and observational arguments are in favor of the possible jetted nature of GRBs (Kumar & Zhang 2015). In this respect, the GRB prompt emission energy has to be corrected for a collimation factor f = 1 − cos θ_jet, where θ_jet is the jet opening angle inferred from an achromatic jet break observed in the GRB afterglow light curve at a characteristic time t_jet (see, e.g., Sari et al. 1999; Frail et al. 2001). However, the estimate on θ_jet depends on specific assumptions on the circumburst medium. We thus chose the case of a homogeneous medium (Frail et al. 2001). By correcting the radiated energy as E_γ = f ⋅ E_iso, we obtain the E_p − E_γ, or “Ghirlanda” correlation (Ghirlanda et al. 2004b, 2007):

$$\begin{array}{c}\hfill log(E_{\mathrm{p}}/\mathrm{keV})=b_0+b_1[log(E_{\gamma }/\mathrm{erg})-50],\end{array}$$(2)

where b₀ and b₁ are the two calibration constants. The correlation features an extra scatter: σ_b.

2.3. The Yonetoku correlation

The Yonetoku correlation involves E_p and the peak luminosity L_p being computed from the time interval of 1 s around the most intense peak of the burst light curve. The peak luminosity is computed as L_p = 4πd_LF_p, where F_p is the observed 1 s peak flux in the rest frame 30−10⁴ keV energy band (Yonetoku et al. 2004; Ito et al. 2019). The L_p − E_p, or “Yonetoku" correlation is written as

$$\begin{array}{c}\hfill log(L_{\mathrm{p}}/\mathrm{erg}{\mathrm{s}}^{-1})-52=m_0+m_{1}log(E_{\mathrm{p}}/\mathrm{keV}),\end{array}$$(3)

where m₀ and m₁ are the two calibration constants. The Yonetoku correlation features an extra scatter: σ_m.

2.4. Combo correlation

The combo correlation represents a hybrid technique since it relates the prompt emission quantity E_p, inferred from the GRB γ-ray spectrum, with the observables from the X-ray afterglow light curve, such as the rest-frame 0.3–10 keV plateau luminosity L₀, its rest-frame duration τ, and the late power-law decay index α (Izzo et al. 2015). Its expression is given by

$$\begin{array}{c}\hfill log\left(\frac{L_0}{\mathrm{erg}{\mathrm{s}}^{-1}}\right)=q_0+q_{1}log\left(\frac{E_{\mathrm{p}}}{\mathrm{keV}}\right)-log\left(\frac{\tau /\mathrm{s}}{|1+\alpha |}\right),\end{array}$$(4)

where q₀ and q₁ are the two calibration constants. The correlation is characterized by an extra scatter: σ_q.

3. Calibration of GRB correlations

Calibrating the above correlations permits us to use GRBs as distance indicators (Wang & Wang 2019; Montiel et al. 2020; Wei & Melia 2020). Clearly, we cannot consider GRBs as standard candles, although we aim to alleviate the circularity problem by adopting model-independent calibrations by means of data set that is independent enough of H(z), for which the luminosity distance d_L is not computed by assuming an a priori cosmological model.

In this paper, we consider a robust method for calibrating GRB correlations that we proposed in Amati et al. (2019). In particular, we propose a “differential age method”. Spectroscopic measurements of the age difference Δt and redshift difference Δz of couples of passively evolving galaxies imply that Δz/Δt ≡ dz/dt. Under the hypothesis that galaxies form at the same redshift, the Hubble rate can be written as H(z) = − (1+z)⁻¹Δz/Δt. So, using the updated sample of 31 Hubble data points, we approximate H(z) using model-independent Bézier parametric curves.

These curves are stable at the lower degrees of control points. They can be rotated and translated by performing the operations on the points and assuming a degree n. They are model-independent reconstructions of a given function, and they formally read:

$$\begin{array}{c}\hfill H_n(z)={\displaystyle \sum _{d=0}^n{\beta }_d{h}_n^d(z),h_n^d(z)\equiv \frac{n!{(z/z_m)}^d}{d!(n-d)!}{(1-\frac{z}{z_m})}^{n-d}.}\end{array}$$(5)

The coefficients β_d of the linear combination of Bernstein basis polynomials $h_n^d(z)$ are positive in the range of 0 ≤ z/z_m ≤ 1, where z_m is the maximum z of our Hubble data set. Aside from cases of constant (n = 0) and linear growth (n = 1) and oscillatory behaviors (n >  2), the only possible combination of Bézier polynomials leading to a nonlinear monotonic growing function up to z_m is n = 2, for which the Hubble function writes as

$$\begin{array}{c}\hfill H_2(z)={\beta }_0{(1-\frac{z}{z_m})}^2+2{\beta }_1\left(\frac{z}{z_m}\right)(1-\frac{z}{z_m})+{\beta }_2{\left(\frac{z}{z_m}\right)}^2,\end{array}$$(6)

where, in units of km s⁻¹ Mpc⁻¹, we have β₀ = 67.76 ± 3.68, β₁ = 103.3 ± 11.1, and β₂ = 208.4 ± 14.3. Our results are portrayed in Fig. 1, where H₂(z) is compared to the OHD data points.

Fig. 1. H2(z) function (solid thick blue) fitting OHD data (black points), with the 1σ (blue curves and shaded area) and 3σ (blue dashed curves) confidence regions. Reproduced from Amati et al. (2019).

After having approximated H(z) with Eq. (6), we impose a curvature parameter Ω_k = 0 in agreement with current indications on spatial curvature made by Planck Collaboration VI (2020), so that the luminosity distance becomes

$$\begin{array}{c}\hfill d_{\mathrm{cal}}(z)=c(1+z){\displaystyle {\int }_0^z{\displaystyle \frac{\mathrm{d}{z}^{\prime }}{H_2(z^{\prime })}},}\end{array}$$(7)

enabling the calibration of E_iso, E_γ, L_p, and L₀ in the above GRB correlations and the determination of their calibration constants. Results are shown in Fig. 2 and summarized in Table 1.

Fig. 2. Best fits (solid lines), their 1σ (dark shaded areas) and 3σ (light shaded areas) dispersions, and associated data sets (black points) of the calibrated GRB correlations considered in this work. Top-left: Amati correlation; top-right: Ghirlanda correlation; bottom-left: Yonetoku correlation; bottom-right: combo correlation.

From the calibrated correlations, we are now in a position to compute their corresponding distance moduli from the standard definition: ${\mu }_{\mathrm{GRB}}^{\mathrm{obs}}=25+5log(d_{\mathrm{cal}}/\mathrm{Mpc})$.

• For the Amati correlation, from Eq. (1), we obtain

  $$\begin{array}{cc}& {\mu }_{\mathrm{A}}^{\mathrm{obs}}=32.55+\frac{5}{2}[\frac{1}{a_1}log\left(\frac{E_{\mathrm{p}}}{\mathrm{keV}}\right)-\frac{a_0}{a_1}\hfill \\ \hfill & -log\left(\frac{4\pi S_{\mathrm{bolo}}}{{\mathrm{erg}}^{-1}{\mathrm{cm}}^{-2}}\right)+log(1+z)],\hfill \end{array}$$(8)

  where the number in the formula takes into account the normalization of E_iso and the conversion from Mpc to cm units.

• For the Ghirlanda correlation, from Eq. (2), we obtain

  $$\begin{array}{cc}& {\mu }_{\mathrm{G}}^{\mathrm{obs}}=27.55+\frac{5}{2}[\frac{1}{b_1}log\left(\frac{E_{\mathrm{p}}}{\mathrm{keV}}\right)-\frac{b_0}{b_1}\hfill \\ \hfill & -log\left(\frac{4\pi f{S}_{\mathrm{bolo}}}{{\mathrm{erg}}^{-1}\mathrm{cm}{}^{-2}}\right)+log(1+z)],\hfill \end{array}$$(9)

  where the number in the formula takes into account the normalization of E_γ and the conversion from Mpc to cm units.

• For the Yonetoku correlation, from Eq. (3), we obtain

  $$\begin{array}{cc}& {\mu }_{\mathrm{Y}}^{\mathrm{obs}}=32.55+\frac{5}{2}[m_0+m_{1}log\left(\frac{E_{\mathrm{p}}}{\mathrm{keV}}\right)\hfill \\ \hfill & -log\left(\frac{4\pi F_{\mathrm{p}}}{{\mathrm{erg}}^{-1}\mathrm{cm}{}^{-2}{\mathrm{s}}^{-1}}\right)+log(1+z)],\hfill \end{array}$$(10)

  where the number in the formula takes into account the normalization of L_p and the conversion from Mpc to cm units.

• For the combo correlation, from Eq. (4), we write

  $$\begin{array}{cc}& {\mu }_{\mathrm{C}}^{\mathrm{obs}}=-97.45+\frac{5}{2}[q_0+q_{1}log\left(\frac{E_{\mathrm{p}}}{\mathrm{keV}}\right)-log\left(\frac{\tau /\mathrm{s}}{|1+\alpha |}\right)\hfill \\ \hfill & -log\left(\frac{4\pi F_0}{{\mathrm{erg}}^{-1}\mathrm{cm}{}^{-2}{\mathrm{s}}^{-1}}\right)],\hfill \end{array}$$(11)

  where the luminosity distance is expressed in cm.

The above four GRB Hubble diagrams can be employed for the fits described in the following sections.

3.1. Connection with cosmokinematics

By redefining a(t) in terms of the redshift, z = a⁻¹ − 1 and employing the cosmological principle only, one can relate the Hubble parameter to the luminosity distance d_L(z) by

$$\begin{array}{c}\hfill H(z)={\left\{\frac{d}{\mathrm{d}z}\left[\frac{d_{\mathrm{L}}(z)}{1+z}\right]\right\}}^{-1},\end{array}$$(12)

which results in

$$\begin{array}{cc}& d_L(z)=\frac{z}{H_0}{\displaystyle \sum _{n=0}^N\frac{{\alpha }_n}{n!}{z}^n,}\hfill \end{array}$$(13)

$$\begin{array}{cc}& H(z)={\displaystyle \sum _{m=0}^M\frac{H^{(m)}}{m!}{z}^m.}\hfill \end{array}$$(14)

By construction, fixing a given expansion order N in d_L means that the corresponding one for H is N − 1. Thus, taking N = 4, more explicitly, we have

$$\begin{array}{c}\hfill d_{\mathrm{L}}^{(4)}\simeq \frac{z}{H_0}({\alpha }_0+{\alpha }_{1}z+{\alpha }_2\frac{z^2}{2}+{\alpha }_3\frac{z^3}{6}),\end{array}$$(15)

or in terms of H:

$$\begin{array}{c}\hfill H^{(3)}\simeq H_0(1+H_{1}z+H_2{\displaystyle \frac{z^2}{2}}+H_3{\displaystyle \frac{z^3}{6}}).\end{array}$$(16)

The derivatives that enter as coefficients α_i, H_i are the cosmographic sets that can be directly compared with cosmic data and observations. Thus, we formally have:

$$\begin{array}{cc}& {\alpha }_i\equiv {\alpha }_i(q_0,j_0,s_0),\hfill \end{array}$$(17)

$$\begin{array}{cc}& H_i\equiv H_i(q_0,j_0,s_0),\hfill \end{array}$$(18)

where up to the fourth order in d_L, the cosmographic set, q₀, j₀, s₀, is immediately provided by current values of the derivatives of a(t) by:

$$\begin{array}{cc}& H(t)=+\frac{1}{a}\frac{\mathit{da}}{\mathrm{d}t},\hfill \end{array}$$(19)

$$\begin{array}{cc}& q(t)=-\frac{1}{a}\frac{d^{2}a}{\mathrm{d}{t}^2}{\left[\frac{1}{a}\frac{\mathit{da}}{\mathrm{d}t}\right]}^{-2},\hfill \end{array}$$(20)

$$\begin{array}{cc}& j(t)=+\frac{1}{a}\frac{d^{3}a}{\mathrm{d}{t}^3}{\left[\frac{1}{a}\frac{\mathit{da}}{\mathrm{d}t}\right]}^{-3},\hfill \end{array}$$(21)

$$\begin{array}{cc}& s(t)=+\frac{1}{a}\frac{d^{4}a}{\mathrm{d}{t}^4}{\left[\frac{1}{a}\frac{\mathit{da}}{\mathrm{d}t}\right]}^{-4},\hfill \end{array}$$(22)

often referred to as cosmographic parameters. The above terms are the Hubble deceleration, jerk, and snap parameters, respectively (Cattoën & Visser 2007, 2008).

Cosmography is the part of cosmology that, using those coefficients, aims to reconstruct the Universe’s kinematics in a model-independent way. For these reasons, we very often refer to this approach as cosmokinetics. Equations (19)–(22) lead to errors that are produced by truncation, and mathematically the convergence decreases when additional higher order terms are involved. This fact clearly compromises the overall accuracy, and, motivated by this issue, we need to study different sets of parameters with fixed orders. An example would be to employ particular hierarchies among coefficients. We describe the hierarchical procedure in the next sections, and we highlight how cosmokinematics can be formulated at low and high redshifts by means of different reformulations of the redshift variable.

3.2. Low-redshift kinematics

The luminosity distance can be connected to the cosmographic series by directly computing the coefficients that enter the previous expansions given by Eq. (17). We supposed above that the Universe is spatially flat, meaning that it has a vanishing curvature (k = 0) and that the cosmological principle holds. This leads to the fact that we did not consider the energy momentum budget of the Universe or its fundamental compositions. We are, however, interested in writing up the coefficients in Eq. (17), taking α₀ = 1, and

$$\begin{array}{cc}& {\alpha }_1=\frac{1}{2}(1-q_0),\hfill \end{array}$$(23)

$$\begin{array}{cc}& {\alpha }_2=-\frac{1}{6}(1-q_0-3q_0^2+j_0),\hfill \end{array}$$(24)

$$\begin{array}{cc}& {\alpha }_3={\displaystyle \frac{1}{24}}(2-2q_0-15q_0^2-15q_0^3+5j_0+10q_0{j}_0+s_0),\hfill \end{array}$$(25)

and

$$\begin{array}{cc}& H_1=1+q_0,\hfill \end{array}$$(26)

$$\begin{array}{cc}& H_2={\displaystyle \frac{1}{2}}(j_0-q_0^2),\hfill \end{array}$$(27)

$$\begin{array}{cc}& H_3={\displaystyle \frac{1}{6}}[-3q_0^2-3q_0^3+j_0(3+4q_0)+s_0],\hfill \end{array}$$(28)

corresponding to the direct derivative of the functions involved, meaning d_L and H, with respect to time.

3.3. High-redshift cosmography

Although powerful, the above formalism suffers from shortcomings due to the convergence at higher redshifts. The violation of z ≃ 0 implies that series lose their meanings because they are computed at our time only. In other words, the standard cosmographic approach fails to be predictive if one employs data at higher redshift domains, which is exactly the case of GRBs. It is, in principle, inconvenient to obtain information on the evolution of cosmographic coefficients when one uses data exceeding the previous redshifts. What is trickier with the approaches to high-redshift cosmography is that one has to write up new techniques to rephrase the definition of cosmological distances.

To further analyze this severe restriction, one can perform a direct high-redshift cosmography by inventing smaller convergence radii and/or theoretical series that are different from Taylor ones, making it possible to heal the problem of convergence arising for cosmography at high redshifts (Capozziello et al. 2020).

In summary, the strategies used to address the issues of high-redshift cosmography are, so far: (1) extending the limited convergence radii of Taylor series by changing variables of expansion (Capozziello et al. 2019) and (2) changing the way in which the expansions are performed, through smooth functions, typically ratios of polynomials, which better converge at higher redshift (Aviles et al. 2014). Below, we analyze the first and second approaches in detail and develop the use of auxiliary variables and rational approximations, respectively.

In the case of auxiliary variables, one employs a tricky method in which the expansion variable is reformulated, under the condition that the new variable satisfies a peculiar construction obeying a few mathematical conditions. The idea is commonly to restrict the convergence radius up to a given value, substituting z, which goes up to infinity, with a new variable that does not. In such a way, the convergence radius is severely reduced in a new Taylor series (reformulated in terms of an independent variable): hereafter y ≡ ℱ(z). All prototypes are functions of the redshift z and obey the following requirements as a minimum:

1. ℱ(z)|_{z = 0} = 0,

2. ℱ(z)|_{z = 0} < ∞.

A more suitable ℱ(z) might address the additional requirement:

1. ℱ(z)|_{z = −1} < ∞.

Following this formula, we considered a first variable under the form y₁ = z/(1 + z), satisfying the first two conditions above. Although appealing, y₁ was categorically ruled out for its inconsistencies in fitting cosmic data (for details, see, e.g., Busti et al. 2015; Capozziello et al. 2018; Rodrigues Filho & Barboza 2018). Thus, a second, more refined one that even includes the further condition 3. can be written as

$$\begin{array}{c}\hfill y_2=arctanz.\end{array}$$(29)

The second one has shown to better perform cosmographic analyses, and, in fact, satisfies the whole set of conditions in toto, since y₂(0) = 0 and y₂(∞) = π/2.

We can compute the coefficients of Eq. (17), using y₂ by

$$\begin{array}{cc}& {\alpha }_1={\displaystyle \frac{1}{2}}(1-q_0),\hfill \end{array}$$(30)

$$\begin{array}{cc}& {\alpha }_2={\displaystyle \frac{1}{6}}(1-j_0+q_0+3q_0^2),\hfill \end{array}$$(31)

$$\begin{array}{cc}& {\alpha }_3={\displaystyle \frac{1}{24}}(10+5j_0-10q_0+10j_0{q}_0-15q_0^2-15q_0^3+s_0),\hfill \end{array}$$(32)

and

$$\begin{array}{cc}& H_1=1+q_0,\hfill \end{array}$$(33)

$$\begin{array}{cc}& H_2={\displaystyle \frac{1}{2}}(j_0-q_0^2),\hfill \end{array}$$(34)

$$\begin{array}{cc}& H_3={\displaystyle \frac{1}{6}}[2+2q_0+3q_0^2+3q_0^3-j_0(3+4q_0)-s_0].\hfill \end{array}$$(35)

Next, we discuss rational approximations. Another way to extrapolate information from improved versions of the series is to extend the convergence radius by changing the way in which expansions are effectively performed. To do so, one can consider rational approximations, meaning well-established prescriptions in extending Taylor expansions in case of a numerical pathology of slopes.

Optimizing the Taylor series with rational approximants leads to the construction of new polynomials that approach infinity at larger z better than the Taylor series. This procedure guarantees mathematical stability of the new series if data points exceed the low redshift limit and better adapts to combined surveys involved in cosmographic computations. Among all the possible choices, Padé and Chebyshev polynomials represent outstanding examples that have been proposed so far (Capozziello et al. 2018). In this work, we consider the Padé polynomials, firstly introduced in Gruber & Luongo (2014). This technique of approximations turns out to be a bookkeeping device to keep the calculations manageable for the cosmography convergence issue.

The only disadvantage of rational approximations is that the orders of Padé expansions are, in principle, not known. So, exploiting the correct expansion orders becomes of fundamental interest in developing our analyses, although the price to pay is that error bars commonly increase with respect to standard Taylor approaches. Intense investigations have been performed throughout recent years (see, e.g., Zhou et al. 2016). To do so, one can estimate the cosmographic parameters and the corresponding error bars circumscribing the analyses to the series that better provide stable results. We therefore limit our analyses below to P_2, 1 Padé approximations. Thus, provided we have Taylor expansions of f(z) under the form $f(z)=\sum _{i=0}^{\infty }{c}_i{z}^i$, with c_i = f⁽ⁱ⁾(0)/i!, it is possible to obtain the (n, m) Padé approximant of f(z) by

$$\begin{array}{c}\hfill P_{n,m}(z)={\displaystyle \left(\sum _{i=0}^n{a}_i{z}^i\right){(1+\sum _{j=1}^m{b}_j{z}^j)}^{-1},}\end{array}$$(36)

which, by construction, is rational and requires that b₀ = 1. Furthermore, it is important that f(z)−P_n, m(z) = 𝒪(z^{n + m + 1}) and the coefficients b_i come from solving the homogeneous system of linear equations $\sum _{j=1}^m{b}_j{c}_{n+k+j}=-b_0{c}_{n+k}$, valid for k = 1, ⋯, m. Once b_i are known, a_i can be obtained using the formula $a_i=\sum _{k=0}^i{b}_{i-k}{c}_k$. Finally, we get

$$\begin{array}{c}\hfill P_{2,1}(z)=\frac{z}{H_0}\left\{\frac{6(q_0-1)+[q_0(8+3q_0)-5-2j_0]z}{-2(3+z+j_{0}z)+2q_0(3+z+3q_{0}z)}\right\}.\end{array}$$(37)

We take into account low- and high-redshift cosmographic methods, and below we show how to obtain results from Monte Carlo procedures using the two aforementioned hierarchies.

4. Methods of statistical analysis

Our cosmographic fits utilize high-redshift data from GRBs. In particular, the points are in the form of Hubble diagram data, obtained from the calibrated correlations described by the numerical coefficients prompted in Table 1. To perform our fits, we considered Markov chain Monte Carlo (MCMC) simulations sampled from the widest possible parameter space by means of a modified version of the freely available Wolfram Mathematica code provided in Arjona et al. (2019). We employed the widely adopted Metropolis-Hastings algorithm. We can therefore compare our approach with several others developed in the literature, and we reduce the dependence on initial statistical distributions made in the simulations by means of the properties of the algorithm itself. The numerical procedure, in particular, is summarized in detail in what follows below for the data sets composed by supernovae, BAO, and GRBs.

4.1. Numerical procedures

We computed the best fit cosmological parameters, namely x, by minimizing the total χ² function, built up from the GRB, SNe Ia, and BAO data sets. We define it as

$$\begin{array}{c}\hfill {\chi }_{\mathrm{tot}}^2={\chi }_{\mathrm{GRB}}^2+{\chi }_{\mathrm{SN}}^2+{\chi }_{\mathrm{BAO}}^2.\end{array}$$(38)

Below, we define the χ² for each of the considered probes.

– For each of the considered GRB correlations, we define

$$\begin{array}{c}\hfill {\chi }_{\mathrm{GRB}}^2={\displaystyle \sum _{i=1}^{N_{\mathrm{GRB}}}{\left[{\displaystyle \frac{{\mu }_{\mathrm{GRB},i}^{\mathrm{obs}}-{\mu }_{\mathrm{GRB}}^{\mathrm{th}}(\mathbf{x},z_i)}{{\sigma }_{{\mu }_{\mathrm{GRB},i}}}}\right]}^2,}\end{array}$$(39)

where N_GRB is summarized in Table 1, ${\mu }_{\mathrm{GRB}}^{\mathrm{th}}$ are the theoretical GRB distance moduli for a given model.

– For SNe Ia, we computed the χ² function from the Pantheon sample, which is the largest and most recent combined sample consisting of 1048 SNe Ia ranging in 0.01 <  z <  2.3 (Scolnic et al. 2018). The SNe Ia distance modulus is parameterized as

$$\begin{array}{c}\hfill {\mu }_{\mathrm{SN}}=m_{\mathrm{B}}-(\mathcal{M}-\alpha {\mathcal{X}}_1+\beta \mathcal{C}-{\mathrm{\Delta }}_{\mathrm{M}}-{\mathrm{\Delta }}_{\mathrm{B}}),\end{array}$$(40)

where m_B and ℳ are the B-band apparent and absolute magnitudes, respectively. Furthermore, 𝒳₁ and 𝒞 are the SN light-curve stretch and the color factors, respectively; α and β are the coefficients of the luminosity-stretch and luminosity-color relationships respectively. Finally, Δ_M is a distance correction based on the host galaxy mass of SNe, and Δ_B is a distance correction based on predicted biases from simulations.

The uncertainties of each SN do not depend on ℳ, therefore, the χ² of the SN data is given by

$$\begin{array}{c}\hfill {\chi }_{\mathrm{SN}}^2={(\mathrm{\Delta }{\mu }_{\mathrm{SN}}-\mathcal{M}\mathbf{1})}^{\mathrm{T}}{\mathbf{C}}^{-1}(\mathrm{\Delta }{\mu }_{\mathrm{SN}}-\mathcal{M}\mathbf{1}),\end{array}$$(41)

where $\mathrm{\Delta }{\mu }_{\mathrm{SN}}\equiv {\mu }_{\mathrm{SN}}-{\mu }_{\mathrm{SN}}^{\mathrm{th}}(\mathbf{x},z_i)$ is the module of the vector of residuals, and C is the covariance matrix that contains statistical and systematic uncertainties on the light-curve parameters (Conley et al. 2011). The parameter ℳ can be removed from the fits by analytically marginalizing it and assuming a flat prior Goliath et al. (2001), leading to

$$\begin{array}{c}\hfill {\chi }_{\mathrm{SN},\mathcal{M}}^2=a+log\frac{e}{2\pi }-\frac{b^2}{e},\end{array}$$(42)

where $a\equiv \mathrm{\Delta }{\mathit{\mu }}_{\mathrm{SN}}^T{\mathbf{C}}^{-1}\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{SN}}$, $b\equiv \mathrm{\Delta }{\mathit{\mu }}_{\mathrm{SN}}^T{\mathbf{C}}^{-1}\mathbf{1}$, and e ≡1^TC⁻¹1. Analytical marginalizations over α and β were not performed, because they contribute to the uncertainties of each SN.

– The BAO is produced by the propagation of sound waves in the early Universe. These waves are observed as a peak in the large-scale structure correlation function and provide a measure of angular distance. We track these waves by means of the comoving volume variation $D_V^3(\mathbf{x},z)$ at a given z, through the BAO observable for uncorrelated data d_z(x, z), namely

$$\begin{array}{c}\hfill D_V^3(\mathbf{x},z)\equiv \frac{cz}{H(\mathbf{x},z)}{\left[\frac{d_{\mathrm{L}}(\mathbf{x},z)}{1+z}\right]}^2,d_z(\mathbf{x},z)\equiv \frac{r_{\mathrm{s}}(z_{}\text{d})}{D_V(\mathbf{x},z)}.\end{array}$$(43)

However, BAO points are slightly model dependent. The comoving sound horizon r_s(z_d) in Eq. (43) depends upon the baryon drag redshift z_d and needs to be calibrated with CMB data for a given cosmological model. For the MCMC simulations, we used Planck Collaboration XVI (2014) best fit values z_d = 1059.62 ± 0.31 and r_s(z_d) = 147.41 ± 0.30.

We used the BAO data shown in Table 2. We did not include correlated BAO from the WiggleZ data (Blake et al. 2011), because the corresponding observable depends on the cosmological parameter Ω_m. This permits largely model-independent BAO points, through which we obtain

$$\begin{array}{c}\hfill {\chi }_{\mathrm{BAO}}^2={\displaystyle \sum _{i=1}^{N_{\mathrm{BAO}}}{\left[\frac{d_{z,i}^{\mathrm{obs}}-d_z^{\mathrm{th}}(\mathbf{x},z_i)}{{\sigma }_{d_{z,i}}}\right]}^2.}\end{array}$$(44)

4.2. Technical procedures of numerical analyses

In view of our numerical analyses, we considered Amati, Ghirlanda, Yonetoku, and combo correlations. For every correlation, we investigated the cosmographic approach by adopting the strategies of expanding in Taylor series first and then with y₂ and rational approximations. As discussed, we considered up to the first order the Taylor series, whereas we employed a y₂ variable and P_2, 1 Padé rational polynomial. In doing so, we explored the ranges of parameters reported in Table 3. The intervals were built up to guarantee that the minima occurring in our computations are unique and correspond to real minima in the statistical analyses. The spread over the polynomial priors corresponds to the increased complexity in analyzing a series. The difficulty of getting results from the aforementioned definitions of z, y₂, and P_2, 1 is compatible with previous efforts made in constraining the cosmological distances. We expand up to the fourth order d_L and up to the third order d_rmz. In the latter case, however, the expansion is decreased by the term ∼z⁻¹ that enters Eq. (43). Moreover, to constrain up to the jerk parameter within hierarchy 1, we decreased the order of expansion of a unity. The same procedure is valid even for the auxiliary variable y₂. For the sake of clarity, in the case of rational approximation, we notice a favorite fitting stability when including a Taylor d_z with the Padé expansion of d_L. Thus, a mixed version of our fitting procedure seems to provide suitable results, which means that stable numerics are obtained at the end of the computations. The strategies are compatible with what has been proposed so far in the literature in the case of rational approximations. The procedures enable us to explore the widest ranges of degrees of freedom by means of the aforementioned priors. No poles in the rational approximations occur for the spanned ranges of coefficients employed here. The numerical bounds inferred are therefore reliable for checking the goodness of any cosmological paradigm. In what follows, we give details on our numerical expectations and on their consequences in view of our understanding toward the cosmological scenario.

5. Discussing numerical estimates

The strategy adopted here is to employ SNeIa, BAO, and GRB data sets in fitting the hierarchies 𝒜₁ and 𝒜₂. We did not consider using GRBs alone, although our calibrations would permit us to consider them as standardized. In fact, in support of our choice, if we considered GRBs alone, the final output over H₀ would be unconstrained in analogy to what happens for SNe. Moreover, for the sake of clarity, the controversy over the H₀ measurements is therefore unsolved in this scheme, even when adding more than one data set. This occurs since the H₀ value, thus far obtained from the OHD-based calibration procedure, turns out to be compatible with the current estimates by the Planck Collaboration VI (2020) and in agreement at the 1.49σ level with the value measured by Riess et al. (2018), respectively. Moreover, our final results indicate a non definitive concordance with H₀ with respect to the two accepted results over Hubble’s parameter today.

Besides the structural form of these correlations, Taylor outcomes are effectively stable within each hierarchy, as portrayed by the results in Table 4 and Figs. 3 and 4, respectively. In hierarchy 1, H₀ values from Amati, Ghirlanda, and Yonetoku correlations are consistent within 1 − σ with the Riess et al. (2018) prediction, meaning (73.4 ± 1.7) km s⁻¹Mpc⁻¹; the H₀ value from the combo correlation seems to be similar to Planck Collaboration VI (2020) expectations, meaning H₀ = (67.4 ± 0.5) km s⁻¹Mpc⁻¹, though it is still consistent with Riess et al. (2018) within 2 − σ. The deceleration parameter q₀ is compatible with a flat Λ cold dark matter (ΛCDM) paradigm, with a mass density parameter of Ω_m ∼ 0.3 for the combo correlation, whereas the other correlations seem to indicate smaller values. In particular, the Yonetoku correlation indicates a value quite far from the other outcomes. Quite surprisingly, the values of j₀ for every correlation turns out to be incompatible with flat ΛCDM predictions. Stronger deviations are found for the combo correlation, meaning the same correlation for which we found q₀ to be in agreement with the flat concordance model in the opposing case. In hierarchy 2, Taylor fits show great concordance with Planck Collaboration VI (2020) expectations on h₀, and clearly viable values over q₀ if compared with a flat ΛCDM model. Lower values for j₀ than those in hierarchy 1 are found instead. Even at 2σ, it is not possible to reach the limit j₀ = 1 prompted in the ΛCDM model. Within 1 − σ, the snaps s₀ are perfectly in agreement with the flat ΛCDM paradigm, which means s₀ = −0.35.

Fig. 3. Contours of our MCMC analyses for hierarchy 1 Taylor expansions. Best fits (black points), prompted with 1 − σ (darker areas) and 2 − σ (lighter areas) confidence levels, refer to SN+BAO+GRB, in particular: Amati (top left), Ghirlanda (top right), Yonetoku (bottom left), and combo (bottom right) data sets. The results indicate a slight spread when comparing q0 and j0 values, with approximatively Gaussian distributions.

Fig. 4. Contours of our MCMC analyses for hierarchy 2 Taylor expansions. Symbols, colors, and positions of the panels are the same as in Fig. 3.

The adoption of the y₂ variable enlarges h₀ significantly for hierarchies 1 and 2 (see Table 5 and Figs. 5 and 6, respectively). The Amati correlation indicates a quite surprising extreme value for h₀, meaning h₀ ∼ 0.78 for hierarchy 2. The results are, however, quite non-predictive at the level of hierarchy 1. In this case, both q₀ and j₀ are badly constrained, meaning they differ dramatically from all other fits. This signature of unpredictivity is not present as s₀ is introduced. In such a case, q₀ lies on suitable intervals, although j₀ and s₀ are very far from the ΛCDM paradigm. The main caveat with these results is that j₀ is negative, indicating no passage between matter and dark energy-dominated eras. Quite strangely, these results are therefore not exhaustive.

Fig. 5. Contours of our MCMC analyses for hierarchy 1 expansions with the auxiliary variable y2. Symbols, colors, and positions of the panels are the same as in Fig. 3.

Fig. 6. Contours of our MCMC analyses for hierarchy 2 expansions with the auxiliary variable y2. Symbols, colors, and positions of the panels are the same as in Fig. 3.

Padé fits seem to improve the quality of hierarchy 1, as was expected by construction (see Table 6 and Fig. 7). In our fits, we find agreement with the Planck Collaboration VI (2020) expectations for h₀. Far from the ΛCDM approach, one notices the values of q₀ and j₀. As well as the Taylor fits, again, indications toward departures from the concordance paradigm seem to be more robust than y₂. Here, we notice that we could not go further to j₀, since introducing s₀ leads to an order ≥(3, 1). The polynomial P_3, 1 is however unconstrained at high redshift and does not perform properly as demonstrated by Capozziello et al. (2020).

Fig. 7. Contours of our MCMC analyses for hierarchy 1 Padé expansions. Symbols, colors, and positions of the panels are the same as in Fig. 3.

Before concluding, we discuss the statistical outcomes of our analysis, as summarized in Table 7. Looking at the χ² values of the cosmographic fits performed over the approximants considered in this work, we notice a possible hierarchy among the correlations themselves. In particular, by sorting the correlations by the ratio χ²/DoF for Taylor 𝒜₁ expansion (see Table 7), we immediately see that, among the considered correlations, combo has the least and closest to unity ratio, while Amati has the largest one. It is also worth mentioning that among the GRB correlations, the Padé rational approximant P_2, 1 represents a better fit with respect to the corresponding Taylor expansion 𝒜₁ for combo and Yonetoku, while for Amati and Ghirlanda it does not. This, in principle, may reflect the nature of the above couples of correlations, meaning that and Yonetoku base their standardization on luminosity measurements, while Amati and Ghirlanda base theirs on energy values, and may in principle indicate among them the best suited observable to perform cosmological fits with GRBs. However, though the best fit results are encouraging, from a statistical significance point of view, the Ghirlanda correlation has the poorest data set, meaning the number of sources fulfilling it corresponds to ∼14% of the Amati data set (see Table 1). This argument in principle may apply to the combo correlation, because, after Ghirlanda, it has the fewest data sets among the correlations considered in this work. Future updates on these two data sets will clarify this issue and strengthen the validity of the correlations. Furthermore, another argument against the Ghirlanda correlation consists of the model-dependent correction over E_iso accounting for GRB collimated emission, which is at the core of the correlation. Thus, theoretically speaking, Ghirlanda should be the least bounded correlation among all cases and for both the hierarchies. For a comparison of the results summarized in Tables 4–7 and described in recent literature see, for example, Demianski et al. (2017b).

Concluding our discussions regarding results, it is worth mentioning that our four standardized GRB empirical correlations, despite boasting high correlation coefficients, are not fully immune to systematic errors. To better explain the role of systematics in our analyses, we must first state that all likely systematic errors entering our scenario might come from two different “domains”. The first deals with the Bézier interpolation method that makes use of OHD data, while the second uses data directly from GRBs. Concerning the first case concerning OHD data, we can stress that these contributions have been extensively studied in the literature. In particular, in this case, data depend on stellar metallicity estimates, population synthesis models, progenitor biases, and the presence of an underlying young component in massive and passively evolving galaxies. Even though severe systematics are expected, several studies found that H(z) are recovered with a 3% error rate at intermediate redshifts and that systematic errors are safely kept below the statistical ones (Moresco et al. 2012, 2018; Liu et al. 2016). Indeed, neglecting the systematics, the statistical errors are incorporated in GRB distance moduli through their calibrated luminosity distances and best fit correlation parameters.

Concerning GRB data, the situation becomes more complex and requires additional and more detailed considerations. Besides the intrinsic dispersions of the GRB’s prompt emissions, other factors induce systematic errors. The first contribution comes from the detection from different γ-ray detectors characterized by various thresholds and spectroscopic sensitivity, which can spread relevant selection effects/biases in the correlation observables. The GRB detector trigger algorithms are based on photon counts and are number-biased against hard photons, which represent the minority in GRB spectra. This implies that, in general, a dim soft burst will be preferentially detected over a dim hard burst, creating apparent correlations where increasingly hard GRBs appear increasingly bright (see, e.g., Shahmoradi & Nemiroff 2011). Another contribution involves prompt emission again, in particular the fitting of GRB spectral energy distribution which appears to depend upon the data quality. A cut-off power-law function is an acceptable fit to GRB spectra if the detector does not have enough sensitivity to collect high-energy photons; however, the resulting E_p are always harder than the one inferred from the standard band function (for details, see Kaneko et al. 2006). A further source of bias arises from the intrinsic GRB luminosity, because GRB correlations require bright enough sources to have their actual redshifts determined, and, in fact, there is an evident lack of weak events at high redshifts (Amati & Della Valle 2013). For hybrid correlation, for example, combo, there is a further complication in combining X-ray observables with γ-ray ones, thus mixing properties from different rest-frame energy bands. Finally, all correlations lack the physical parameter information related to their theoretical explanation. In summary, all correlations appear to be significantly affected by selection effects in trigger and spectral analyses, and redshift determinations.

In view of the above considerations, in the past literature there have been many cautions and several studies on the possible selection effects of GRB correlations, leading, however, to contrasting results (Butler et al. 2007; Ghirlanda et al. 2008; Nava et al. 2008; Shahmoradi & Nemiroff 2011; Amati et al. 2009; Yonetoku et al. 2010). In particular, Wang (2008) noted that in four GRB correlations (including Ghirlanda and Yonetoku, which are considered in this work), though the statistical errors on the best fit correlation parameters are quite small, their corresponding χ² are very large and dominated by systematic errors. The author estimated systematic errors for all correlations by requiring that χ² ≡ DoF. Then, within ΛCDM models with different values of Ω_m, Wang (2008) noted that the derived systematic errors change by less than 3%. This result means that either the systematic uncertainties of the calibrated correlation parameters are not sensitive to cosmological parameter values, or conversely, that the inclusion (or lack thereof) of systematic errors does not significantly affect cosmological fits performed by including GRB data. In view of this result, we can quantify the influence of systematic errors on our estimates to be ≈3%. Discussions on how to go further with the concordance paradigm are the subject of the next section.

6. Going further the concordance paradigm?

Our current understanding of the universe lies on the assumption of the ΛCDM paradigm, which aims to describe large-scale dynamics by means of a vacuum energy cosmological constant Λ. Its experimental estimate is, however, dramatically different from quantum field predictions. Recently, the need to revise this approach has culminated in severe discrepancies between different measures of the current Hubble rate. There is an apparent (9 ± 2.5)% tension between H₀ measured by the Planck Collaboration VI (2020) and Riess et al. (2018), respectively. Other significant evidence in favor of not yet understood new physics comes from the potential driving inflationary times, which is written in excellent approximation in terms of a Starobinsky R + αR² correction of Einstein’s gravity (Planck Collaboration VI 2020; Calzá et al. 2020). The question is thus: are there new physics beyond the concordance paradigm? In our approach, we consider GRBs of which the redshift domain is far from being primordial, as for inflationary epochs, but it lies on intermediate Hubble evolution. Consequently, comparing our results and confronting them with the ΛCDM paradigm could indicate possible hints toward deviations from the standard paradigm. Looking at our results summarized in Sect. 5, it seems that the standard paradigm is poorly constrained and fails to be predictive when GRBs are involved. In particular, while the auxiliary variables are far from being suitable benchmarks of our conclusions, both Taylor and Padé approximations show that q₀ and s₀ are roughly comparable to the ΛCDM model, whereas j₀ is not. The issue over H₀ is unsolved since many measures seem to agree with the Planck Collaboration VI (2020) expectations, but many others with those of Riess et al. (2018). We conclude that, besides the H₀ tension, the concordance paradigm is only partially recovered if one assumes the above results are not jeopardized by high systematics. In this case, we discussed above that it might not contribute significantly to the overall expectations. However, even if our indications seem to be in favor of a model practically comparable to the ΛCDM approach, but with varying j, more accurate studies are needed to check whether systematics really affect the results or not. Even though this represents a suitable possibility, another aspect to be considered could be to take a spatial curvature different from zero. This scenario arises because s can also vary according to time differently to the way it does so according to the ΛCDM framework, and our bounds do not, a priori, exclude the occurance of slight departures. This is why the speculation of non-negligible spatial curvature would need to be intensively investigated to exclude or confirm whether GRBs with different calibrated correlations are in favor of a spatially curved universe (Wei & Melia 2020).

7. Final outlooks and perspectives

In this paper, we revised the Universe’s dynamics by using high-redshift data from GRBs. To do so, we constrained cosmographic parameters by means of model-independent techniques through four different GRB calibrations. To overcome the issue of circularity, we proposed a recent new method, based on Beziér polynomials. These have been employed to set correlations to heal the circularity problem, using OHD data catalogs. The main advantage of our method is that it frames our analysis in a model-independent way, healing, de facto, the circularity problem without postulating the model a priori. With the model-independently calibrated Amati, Ghirlanda, Yonetoku, and combo correlations, several MCMC analyses were computed to obtain 1 − σ and 2 − σ confidence levels and to test the standard cosmological model. To do so, we considered the most recent approaches to cosmography, comparing, for example, Taylor expansions with z and y₂ series, and Padé polynomials. We derived limits over the Universe’s expansion history, fixing constraints over h₀, q₀, j₀, and s₀ within two hierarchies, 𝒜₁ up to j₀ and 𝒜₂ up to s₀. We used GRB data from the aforementioned four calibrated GRB correlations, together with SNeIa and BAO catalogs. Reasonable results have been found for 𝒜₁ and 𝒜₂ hierarchies through several MCMC fits, demonstrating possible matching with the standard paradigm. Moreover, we only partially alleviated the tension on local H₀ measurements where the 𝒜₂ hierarchy is considered. Quite surprisingly, our findings showed that the ΛCDM model is not fully confirmed using GRBs. Although this can be an indication that more refined analyses are required where GRBs are involved, simple indications seem to be against a genuine cosmological constant and are interpreted throughout the text, either with a barotropic dark energy contribution, or with the need for nonzero spatial curvature. Nevertheless, at this stage, our findings are in line with recent claims regarding tensions with the ΛCDM model (Yang et al. 2019; Lusso et al. 2019; Risaliti & Lusso 2019). In this respect, we also developed comments on previous works concerning the role of systematic errors. In particular, we studied how systematics can affect our results in view of the cosmographic bounds. We found that no significant deviations are expected in agreement with previous works discussed on the same topic.

Future perspectives will shed light on the role of spatial curvature. We will resolve the experimental fits in which spatial curvature is not fixed to zero at the very beginning. Moreover, we intend to adapt our approach, based on Bézier polynomials, directly to the case of dark energy paradigms. We endeavour to investigate this possibility in order to verify, using GRB data, whether or not our obtained deviations are somehow reproducible in the context of dark energy scenarios.

Acknowledgments

O. L. and M. M. acknowledge the INFN support, and in particular the Frascati National Laboratories, for Iniziative Specifiche MOONLIGHT2. O. L. acknowledges the support provided by the Ministry of Education and Science of the Republic of Kazakhstan, Program IRN: BR05236494 for support. O. L. and M. M. are thankful to Hirotaka Ito for sharing the updated GRB data of the Yonetoku correlation. The work is dedicated to the lovely memory of our colleague, Diego Leonardo Caceres Uribe, excellent scientist and good friend.

References

All Tables

All Figures

	Fig. 1. H2(z) function (solid thick blue) fitting OHD data (black points), with the 1σ (blue curves and shaded area) and 3σ (blue dashed curves) confidence regions. Reproduced from Amati et al. (2019).
In the text

	Fig. 2. Best fits (solid lines), their 1σ (dark shaded areas) and 3σ (light shaded areas) dispersions, and associated data sets (black points) of the calibrated GRB correlations considered in this work. Top-left: Amati correlation; top-right: Ghirlanda correlation; bottom-left: Yonetoku correlation; bottom-right: combo correlation.
In the text

Fig. 3. Contours of our MCMC analyses for hierarchy 1 Taylor expansions. Best fits (black points), prompted with 1 − σ (darker areas) and 2 − σ (lighter areas) confidence levels, refer to SN+BAO+GRB, in particular: Amati (top left), Ghirlanda (top right), Yonetoku (bottom left), and combo (bottom right) data sets. The results indicate a slight spread when comparing q0 and j0 values, with approximatively Gaussian distributions.

In the text

	Fig. 4. Contours of our MCMC analyses for hierarchy 2 Taylor expansions. Symbols, colors, and positions of the panels are the same as in Fig. 3.
In the text

	Fig. 5. Contours of our MCMC analyses for hierarchy 1 expansions with the auxiliary variable y2. Symbols, colors, and positions of the panels are the same as in Fig. 3.
In the text

	Fig. 6. Contours of our MCMC analyses for hierarchy 2 expansions with the auxiliary variable y2. Symbols, colors, and positions of the panels are the same as in Fig. 3.
In the text

	Fig. 7. Contours of our MCMC analyses for hierarchy 1 Padé expansions. Symbols, colors, and positions of the panels are the same as in Fig. 3.
In the text