\section{Introduction}\label{sec:1}

In the cosmological concordance model, the Universe is approximated by $\sim30\%$ of baryonic and cold dark matter and by $\sim70\%$ of an exotic form of constant \emph{dark energy} \citep[DE, see e.g.,][hereafter P20]{Planck2018}. In particular to speed up the Universe today, DE counterbalances the action of gravity through  a negative pressure. Recent observations at small redshifts favor a cosmological constant contribution, $\Lambda$, to evolving DE, albeit  recent observations provided by \citet{Planck2018} seem to show unexpected tensions among observables, not yet understood within the concordance paradigm, dubbed the $\Lambda$CDM model. In the latter, the corresponding equation of state (EoS) turns out to be \emph{exactly} $w\simeq -1$ (\citetalias{Planck2018}; \citealt{Riess2007}). Any deviations from $w=-1$ would lead to more complicated versions of DE fluids or to
extended and/or modified theories of gravity
still objects of debate \citep{Bronstein1933,Tsujikawa2013,Sahni2014,Ding2015,orlando1}.

Several methods have been proposed to investigate a possible evolution of $w(z)$ with the redshift, mainly involving supernovae (SNe) Ia \citep{Phillips1993,1998Natur.391...51P,Perlmutter1999,Riess1998,Schmidt1998} and/or standard candles and rulers. However, the shortage of data at high redshifts ($z\geq1$) inevitably brings large uncertainties over $w$ measurements \citep{Union2.1}, providing an unvoidable degeneracy problem among cosmological models. For example, the simplest $\Lambda$CDM generalization is based on a first-order approximation, $w(z) = w_0 + w_a z/(1+z)$ \citep{Chevallier2001,Linder2003} fully degenerating with several other  redshift-binned parameterizations within the sphere $z<1$ \citep[see, e.g.,][]{King2014}.

Gamma ray bursts (GRBs) have the advantage over other cosmological probes and rulers to homogeneously cover a large interval of redshift up to $z \approx 9$ \citep{Cucchiara2011,2012ApJ...749...68S} and their redshift distribution peaks at $z\sim2$--$2.5$ \citep{Coward2013}, where the farthest SN Ia has been detected \citep{Rodney2015}. In the last decades, the analysis of larger samples has led to various phenomenological correlations between GRB photometric and spectroscopic properties,  suggesting  possible cosmological applications \citep[see e.g.][and references therein]{Schaefer2007,Amati2008,CapozzielloIzzo2008,Dainotti2008,Izzo2009,AmatiDellaValle2013,Wei2014}. Recently \citet{Izzo2015} (hereafter I15) developed a method for measuring the cosmological parameters from a sample of $60$ GRBs, which minimizes the cosmological ``circularity'' problem affecting the computation of the luminosity distances in all GRB correlations\footnote{For alternative techniques see e.g. \citep{orlando2,orlando3,orlando4}.}. This technique is based on the so-called \textit{Combo relation},\footnote{The name \textit{Combo} states for \emph{combination} since it combines $E_{\rm p,i}$--$E_{\rm\gamma,iso}$ \citep{Amati2008} and $E_{\rm\gamma,iso}$--$E_{\rm X,iso}$--$E_{\rm p,i}$ relations \citep{Bernardini2012}.} which is characterized by a small data scatter and involves prompt and afterglow GRB parameters \citep{Bernardini2012,Margutti2012}. The ``candle'' is provided by the plateau X-ray afterglow luminosity $L_0$, which is related to its rest-frame duration $\tau$, to the late power-law decay index $\alpha$ of the X-ray afterglow luminosity, and to the rest-frame peak energy of the GRB prompt emission $E_{\rm p,i}$ (see Fig.~\ref{fig:MDV}). The existence of this relation has been later on confirmed by \citet{2016ApJ...825L..20D} where a similar relation based on the use of the peak luminosity has been derived.
\begin{figure}
\centering
\includegraphics[width=\hsize,clip]{GRB060418_plotOk}
\caption{The rest-frame $0.3$--$10$~keV light curve of GRB 060418A. The total fit (red curve) is composed of a steep decay (dashed purple power-law) and a plateau + late power-law decay (dot-dashed blue curve). The black dots are the data filtered by the flares (blue dots). The vertical black line indicates $\tau$; an arrow marks the luminosity $L_0$.}
\label{fig:MDV}
\end{figure}

In this work, adopting an extended GRB sample, characterized by a complete data set in gamma- and X-rays, we fix limits over DE's evolution. Our general strategy remarks the same of \citetalias{Izzo2015} and aims at getting  novel bounds over the whole matter content, i.e. $\Omega_m$. Afterwards, we fix the cosmological constant density, $\Omega_\Lambda$, testing the validity of the $\Lambda$CDM model. Our findings adopt the estimated values of $L_0$ as distance indicator and can be generalized to any DE scenarios. To this end, we constrain the evolution of $w(z)$ through a piecewise formulation over the GRB redshift intervals.  Theoretical discussions of our results are prompted, leading to an overall agreement with the $\Lambda$CDM model at small and intermediate redshifts. Outside our time, slight deviations from $w=-1$ are not excluded \emph{a priori}. Precisely, from our new GRB sample we get tighter constraints on DE's evolution, showing that at $z\leq 1.2$ DE turns out to be compatible with a cosmological constant, but beyond $z\simeq1.2$ it does not since the case $w<-1$ cannot be fully excluded. These deviations are however severely affected by uncertainties and so further space missions, e.g. {\it Space Variable Objects Monitor} \citep[SVOM,][]{SVOM}, {\it enhanced X-ray Timing and Polarimetry} \citep[eXTP,][]{eXTP}, and {\it Transient High-Energy Sky and Early Universe Surveyor} \citep[THESEUS,][]{2018AdSpR..62..191A}, will significantly enrich the GRBs database to remove the ambiguity over the form of $w(z)$ as we discuss later in the manuscript.

The paper is structured as follows. In Section~\ref{sec:2}, we extend the sample of $60$ GRBs in \citetalias{Izzo2015} with $114$ additional GRBs characterized by a complete data set in gamma- and X-rays up to December 2018. In Section~\ref{sec:3}, we provide new constraints on the cosmological parameters $\Omega_m$ and $\Omega_\Lambda$ using the estimated values of $L_0$ as distance indicator. In particular, in Section~\ref{sec:3.4}, we constrain the evolution of $w(z)$ by considering a piecewise formulation over the GRB redshift interval. The results of our analysis are then discussed in Section~\ref{sec:7},
In Section~\ref{sec:8} we draw our conclusions.

\section{The updated GRB sample}\label{sec:2}

\subsection{Brief summary on the relation and the updated sample}\label{sec:2.1}

The Combo relation writes as \citepalias[see][for details]{Izzo2015}
\begin{equation}
\label{eq:no2}
\log \left(\frac{L_0}{{\rm erg/s}}\right) = \log \left(\frac{A}{{\rm erg/s}}\right) + \gamma \log \left(\frac{E_{\rm p,i}}{{\rm keV}}\right) - \log \left(\frac{\tau/{\rm s}}{|1+\alpha|}\right)\ ,
\end{equation}
where $\gamma$ is the slope and $A$ the normalization.
For each GRB, the rest-frame peak energy $E_{\rm p,i}$ is inferred from the $\nu F(\nu)$ GRB spectrum, while $L_0$, $\tau$, and $\alpha$ by fitting the rest-frame $0.3$--$10$~keV flare-filtered afterglow luminosity light curves with the function $L(t)=(1+t/\tau)^\alpha$ introduced in \citet{RuffiniMuccino2014} (hereafter R14, see, e.g., Fig.~\ref{fig:MDV}).\footnote{The flare-filtered luminosity light curves are iteratively fitted with the function in \citetalias{RuffiniMuccino2014}: at every iteration, data points with the largest positive residual are discarded, until a final fit with a p-value $>0.3$ is obtained.}

From our sample we exclude light curves with $\alpha>-1$, which may have a change in slope beyond the XRT time coverage and/or be polluted by a late flaring activity.
We obtain a sample of {\bf$174$} GRBs, of which $60$ from \citetalias{Izzo2015} and {\bf$114$} introduced in this work (see Table~\ref{tab:no}). The Combo relation within the flat $\Lambda$CDM model is displayed in Fig.~\ref{fig:no}. Following \citet{Dago2005}, we obtain the best-fit parameters: $\log[A/({\rm erg/s})]=49.66\pm0.20$, $\gamma=0.84\pm0.08$, and an extra-scatter $\sigma_{\rm ex}=0.37\pm0.02$. The corresponding Spearman's coefficient is $\rho_s=0.92$ and the p-value from the two-sided Student's $t$-distribution is $p=1.17\times10^{-48}$.
\begin{figure}
\centering
\includegraphics[width=\hsize,clip]{Combo_updated}
\caption{The Combo relation from the sample of $174$ GRBs (black circles) in the flat $\Lambda$CDM model with indicative values $H_0=71$ km/s/Mpc, $\Omega_m=0.27$, and $\Omega_\Lambda=0.73$. The best-fit (solid gray line) and the $1$- and $3$--$\sigma$ of extra-scatter (dotted and dashed gray lines) are also displayed.}
\label{fig:no}
\end{figure}

\subsection{The calibration}\label{sec:2.2}

The Combo relation is calibrated through a two steps method minimizing the use of SNe Ia \citepalias{Izzo2015}. This method is based on:
\begin{itemize}
\item[(1)]{the determination of $\gamma$ and $\sigma_{\rm ex}$ from various, small sub-samples of GRBs lying almost at the same redshift;}
\item[(2)]{the determination of $A$ from SNe Ia located at the same redshift of the nearest GRBs of the sample.}
\end{itemize}
The use of SNe Ia in the step (2) is limited to the lowest redshift where the effect of the cosmology on the distance modulus of the calibrating SNe Ia is negligible (see Fig.~$3$ in \citetalias{Izzo2015}).

\subsubsection{The determination of $\gamma$ and $\sigma_{ex}$}\label{sec:2.2.1}

\begin{figure}
\centering
\includegraphics[width=\hsize,clip]{Combo_slope_MDV}
\caption{The Combo sub-samples (see legend). The solid gray best-fit line and the dotted (dashed) gray lines at $1$--$\sigma$ ($3$--$\sigma$) are obtained from the calibrated values of $\gamma$ and $\sigma_{\rm ex}$ as inferred in Sec.~\ref{sec:2.2.1}. The normalization $A$ is constrained to the value determined in Sec.~\ref{sec:2.2.2}.}
\label{fig:no2}
\end{figure}

\begin{table}
\centering
\caption{Best-fit parameters of the seven sub-samples at average redshift $\langle z\rangle$: the slope $\gamma_z$, the normalization $A_{\rm z}$, and the extrascatter $\sigma_{\rm ex,z}$ are shown.}
\begin{threeparttable}
\begin{tabular}{c|ccc}
\hline\hline
$\langle z\rangle$  &   $\gamma_z$       &  $\log \left(A_{\rm z}/{\rm cm}^2/{\rm s}\right)$    & $\sigma_{\rm ex,z}$ \\
\hline
\hline
$0.54\pm0.01$	&	$0.81\pm0.49$	&	$-7.38\pm1.08$	                                    &	$0.29\pm0.10$	\\
$1.18\pm0.07$	&	$0.83\pm0.31$	&	$-7.93\pm0.77$	                                    &	$0.26\pm0.08$	\\
$1.46\pm0.05$	&	$0.80\pm0.32$	&	$-8.29\pm0.79$	                                    &	$0.26\pm0.08$	\\
$1.70\pm0.06$	&	$0.94\pm0.32$	&	$-8.92\pm0.84$	                                    &	$0.19\pm0.08$	\\
$2.05\pm0.05$	&	$1.05\pm0.32$	&	$-9.44\pm0.90$							            & $0.34\pm0.08$ \\
$2.27\pm0.07$	&	$0.78\pm0.20$	&	$-8.57\pm0.54$							           & $0.16\pm0.06$	\\
$2.69\pm0.08$	&	$1.02\pm0.34$	&	$-9.38\pm0.87$							           & $0.34\pm0.10$ \\
\hline												\hline
\end{tabular}
\end{threeparttable}
\label{tab:no3}
\end{table}

To select the GRB sub-samples at the same $z$ with well constrained best-fit parameters, we impose that the corresponding relations must span at least $2$ order of magnitude in the Combo plane.
Since in each sub-sample the GRB luminosity distances $d_l$ are approximately the same, instead of using $L_0=4\pi d_l^2 F_0$, we employ the rest-frame $0.3$--$10$ keV energy flux $F_0$, rendering this procedure cosmology-independent.
We found seven sub-samples, shown in Fig.~\ref{fig:no2}.
Their best-fit slope parameters $\gamma_z$ (see Table~\ref{tab:no3}) exhibit no evident trend with $z$ within the errors, implying that the Combo relation does not result from an evolution effect. This feature allows us to assume that the above seven sub-samples do follow the same correlation with the same $\gamma$ and $\sigma_{\rm ex}$ but different normalizations. Therefore, we perform a simultaneous fit of the above sub-samples requiring $\gamma$ and $\sigma_{\rm ex}$ to be the same for all of them. The best-fit values, which minimize the sum of the seven log-likelihood functions, are $\gamma=0.90\pm0.13$ and $\sigma_{\rm ex}=0.28\pm0.03$ (see Fig.~\ref{fig:no2}).

\subsubsection{The calibration of $A$}\label{sec:2.2.2}

The Combo relation distance modulus is \citepalias[see, e.g.,][]{Izzo2015},
\begin{align}
\nonumber
\mu_{\rm GRB}=-97.45+\frac{5}{2}&\left[\log\left(\frac{A}{{\rm erg/s}}\right)+\gamma\log \left(\frac{E_{\rm p,i}}{{\rm keV}}\right) -\log\left(\frac{\tau/{\rm s}}{|1+\alpha|}\right)\right. \\
\label{eq:no3}
&\left. -\log\left(\frac{F_0}{\textnormal{erg/cm$^2$/s}}\right)-\log4\pi
\right]\,.
\end{align}
To calibrate $A$ we shall use the nearest GRBs of our sample: 1) 171205A at $z=0.0368$, 2) 180728A at $z=0.117$, and 3) 130702A at $z=0.145$ and 161219B at $z=0.1475$, almost at the same distance.
In Eq.~\ref{eq:no3} we replace $\mu_{\rm GRB}$ with the average distance modulus $\langle\mu_{\rm SNIa}\rangle$ computed from SNe Ia at the same $z$ of the above four sources \citep[][]{Union2,Union2.1}, and we use the value of $\gamma$ from Sec.~\ref{sec:2.2.1} and the values of $F_0$, $E_{\rm p,i}$, $\tau$, and $\alpha$ for each of the four GRBs from Table~\ref{tab:no}. We infer: 1) $\langle\mu_{\rm SNIa}\rangle=35.97\pm0.17$ and $\log[A/{\rm ( erg/s)}]=48.61\pm0.42$ for GRB 171205A; 2) $\langle\mu_{\rm SNIa}\rangle=38.68\pm0.09$ and $\log[A/{\rm ( erg/s)}]=48.97\pm0.25$ for GRB 180728A; 3) $\langle\mu_{\rm SNIa}\rangle=39.21\pm0.24$ and $\log[A/{\rm ( erg/s)}]=49.52\pm0.21$ for GRB 130702A and $\log[A/{\rm ( erg/s)}]=50.09\pm1.12$ for GRB 161219B. Indeed, in estimating the normalization, we select the values of $A$ coming from GRBs 130702A and 161219B because they represent a small set of bursts at the same $z$.\footnote{ Single GRB estimates are disfavored because they might be above/below the relation and, thus, over/underestimate the value of $A$. Small sample estimates in general average these effects.} Their weighted average value is $\log[A/{\rm ( erg/s)}]=49.54\pm0.21$.

\begin{figure}
\centering
\includegraphics[width=\hsize,clip]{Combo_updated_hubblesolo}
\caption{The GRB Hubble diagram of the Combo relation.}
\label{fig:no5a}
\end{figure}

\section{Testing the standard $\Lambda$CDM model}\label{sec:3}

GRB distance moduli $\mu_{\rm GRB}$ and their uncertainties $\sigma\mu_{\rm GRB}$ (see Fig.~\ref{fig:no5a}) depend upon the set of observables $(E_{\rm p,i},\tau,\alpha,F_0)$ and therelation parameters. They can be compared with the theoretical distance moduli $\mu_{\rm th}=-97.45+5\log\{d_l[z,\Omega_m,\Omega_\Lambda,w(z),H_0]/{\rm cm}\}$, where the luminosity distance $d_l$ is defined as \citep[see, e.g.,][]{Goobar1995}
\begin{equation}
\label{eq:no5}
d_l = \frac{c}{H_0}\frac{(1+z)}{\sqrt{| \Omega_k |}} F_k\left(\int_0^z \frac{\sqrt{| \Omega_k |} dz^\prime}{\sqrt{E\left[z^\prime,\Omega_m,\Omega_\Lambda, w\right]}}\right)\ ,
\end{equation}
with $E\left[z,\Omega_m,\Omega_\Lambda,w\right]=\Omega_m(1+z)^3+\Omega_\Lambda(1+z)^{3(1+w)}+\Omega_k(1+z)^2$, $\Omega_k = 1-\Omega_m-\Omega_\Lambda$, and $F_k(x)=\sinh(x)$ for $\Omega_k>0$, $F_k(x)=x$ for $\Omega_k = 0$, and $F_k(x)=\sin(x)\}$ for $\Omega_k < 0$.

\begin{table}
\setlength{\tabcolsep}{0.75em}
\renewcommand{\arraystretch}{1.2}
\centering
\caption{Best-fit parameters with $1$--$\sigma$ uncertainties for the various cosmological cases discussed in this work. The last column lists the values of the $\chi^2$. $H_0$ is fixed to the values given by \citetalias{Planck2018} and \citetalias{2019ApJ...876...85R} (see text).}
\small
\begin{tabular}{lccccc}
\hline\hline
$H_0$ &  $\Omega_m$       &  $\Omega_{\Lambda}$  & $\Omega_k$    & $w$   &  $\chi^2$         \\
\hline
\hline
\multicolumn{6}{c}{Flat $\Lambda$CDM (${\rm DOF}=173$)}\\
\cline{1-6}
P20	   &   $0.32^{+0.05}_{-0.05}$   &  $0.68^{+0.05}_{-0.05}$   & $0$   &  $-1$   &  $165.54$  \\
R19   &   $0.22^{+0.04}_{-0.03}$   &  $0.78^{+0.04}_{-0.03}$   & $0$   &  $-1$   &  $171.32$  \\
\hline
\hline
\multicolumn{6}{c}{$\Lambda$CDM (${\rm DOF}=172$)} \\
\cline{1-6}
P20	   &   $0.34^{+0.08}_{-0.07}$   &  $0.91^{+0.22}_{-0.35}$   & $-0.24^{+0.24}_{-0.35}$  &  $-1$    &  $164.38$  \\
R19   &   $0.24^{+0.06}_{-0.05}$   &  $1.01^{+0.15}_{-0.25}$   & $-0.24^{+0.16}_{-0.25}$  &  $-1$   &  $169.23$ \\
\hline
\hline
\end{tabular}
\label{tab:no5}
\end{table}

\begin{figure*}
\centering
\includegraphics[width=0.9\hsize,clip]{Combo_updated_contour2}
\caption{Best-fit parameters (black circles) and confidence regions ($1$-, $2$-, and $3$-$\sigma$ from the inner/darker to the outer/lighter) in the $\Omega_m$--$\Omega_\Lambda$ plane for the selected $H_0$. The gray regions indicate the set of parameters for which no Big Bang occurs while the dashed lines those of a flat Universe case.}
\label{fig:no5}
\end{figure*}

We first constrain the Hubble constant $H_0$ with two recent estimates, in tension at $4.4$-$\sigma$: $H_0=(67.4\pm0.5)$~km/s/Mpc \citepalias{Planck2018} and $H_0=(74.03\pm1.42)$~km/s/Mpc \citep[][hereafter R19]{2019ApJ...876...85R}.
Then we compare $\mu_{\rm GRB}$ and $\mu_{\rm th}$ distributions and find the best-fit cosmological parameters by performing a $\chi^2$ statistic test.

We first take the flat $\Lambda$CDM model and then we leave $\Omega_k$ free to vary. In the first case we impose $\Omega_{\Lambda}=1-\Omega_m$ and we summarize our findings in Table~\ref{tab:no5}. In particular,
GRBs provide: $\Omega_m=0.32^{+0.05}_{-0.05}$ for the $H_0$ of \citetalias{Planck2018} and $\Omega_m=0.22^{+0.04}_{-0.03}$ for the $H_0$ of \citetalias{2019ApJ...876...85R}, whereas fits using the non-flat $\Lambda$CDM model provide large values of $\Omega_{\rm k}$ which are consistent within the errors with the flat $\Lambda$CDM case (see Table~\ref{tab:no5} and Fig.~\ref{fig:no5}). Our large and negative value of the curvature parameter is compatible with the recent one obtained by \citet{2020ApJ...888...99W} using only quasar data calibrated via observational Hubble data set (OHD). However, in \citet{2020ApJ...888...99W} the best-fit value $\Omega_{\rm k}=-0.918\pm0.429$ has a large uncertainty, due to the error propagations induced by anchored OHD, and it is in disagreement with the concordance model at $2.1$--$\sigma$ level and at $1.6$--$\sigma$ level with respect to our findings.

Worth noticing, the best-fit $\Omega_m$ values listed in Table~\ref{tab:no5}, though obtained by using different values of $H_0$, are consistent with each other.

\section{Piecewise reconstruction of $w(z)$ in  flat background cosmology}\label{sec:3.4}

We consider now a flat Universe with an evolving DE EoS. We adopt a redshift-binned parametrization for $w(z)$: a step function with different and constant values of $w$ in each of the redshift bins over which the entire redshift range covered by GRBs is divided \citep{Daly2004,King2014}.

The number of bins is arbitrary defined but holds a set of heuristic assumptions that we made before deciding the redshift interval for each bin. First, we have to notice that the number of bins requires the number of points that is enough to guarantee a suitable interpolation. Second, we expect the number of points increases with the redshifts since most of the GRBs of the sample are located at larger redshifts. This implies that we cannot choose tighter intervals as $z$ increases to include always the same number of data points. Rather than following this strategy, we presume to work at least with 15 points for each bins and to split the intervals in redshift differences of about $\Delta z\sim0.5$. In so doing, we assume the following four intervals for $\Delta z_i$: $\Delta z_1=0.55$, $\Delta z_2=0.63$, $\Delta z_3=0.56$ and $\Delta z_4=0.81$, whereas we assume the last interval up to $z\sim 9.5$ since in it we count 48 data points. In this respect the number of data for each interval only slightly increases, leaving unaltered the statistical significance of our fits for each interval. The bins are therefore statistically supported and reported in Table~\ref{tab:no6}.

In Eq.~\ref{eq:no5} we can set
\begin{align}
E\left[\Omega_m,f(z_{N-1}<z \leq z_N)\right]=\,\, & \Omega_m(1+z)^3+\\
&(1-\Omega_m)f(z_{N-1}<z \leq z_N)\,, \nonumber
\end{align}
with parametric  DE evolution $f(z_{N-1}<z \leq z_N)$ given by
\begin{equation}
\label{eq:fzbinned}
f(z_{N-1}<z \le z_N)=(1+z)^{3(1+\tilde{w}_N)}\prod_{i=0}^{N-1}\left[1+{\rm max}(z_i)\right]^{3(\tilde{w}_i-\tilde{w}_{i+1})},
\end{equation}
in which $N$ is the redshift bin where $z$ resides and $\tilde{w}_i$ the EoS parameter in the $i^{\mathrm{th}}$ redshift bin defined by its maximum redshift ${\rm max}(z_i)$.
This parametrization makes no assumptions on DE's nature since different parameters are introduced in every redshift bin.
However, the parameters $w_i$ are usually correlated with each other (see Eq.~\ref{eq:fzbinned}).
This means that the covariance matrix
\begin{equation}
\textbf{C}=\langle \tilde{\textbf{w}} \tilde{\textbf{w}}^{\rm T}\rangle-\langle\tilde{\textbf{w}}\rangle\langle\tilde{\textbf{w}}^{\rm T}\rangle\ ,
\end{equation}
computed from the vector $\tilde{\textbf{w}}$ with components $w_i$ and the transposed one $\tilde{\textbf{w}}^{\rm T}$, is not diagonal.
Following \citet{Huterer2005}, to uncorrelate the DE parameters $\tilde{w}_i$ we transform them through an orthogonal matrix rotation \textbf{O} into a basis that diagonalizes the inverse of the covariance matrix, i.e., ${\bf C}^{-1}={\bf O}^{\rm T}{\bf \Lambda} {\bf O}$, where ${\bf\Lambda}$ is diagonal. The new set of uncorrelated DE parameters $w_i$ is then given by applying a weight matrix defined as $\tilde{\textbf{W}}=\textbf{O}^{\rm T}{\bf \Lambda}^{1/2}\textbf{O}$ to the correlated $\tilde{w}_i$. $\tilde{\textbf{W}}$ is normalized in such a way that its rows, the weights for $\tilde{w}_i$, sum to unity.

We first marginalize over the $\Omega_m$ parameter using the values inferred for the flat $\Lambda$CDM case (see Table~\ref{tab:no5}).
Then, we determine five redshift bins to assure a sufficient number of GRB events for each of them (see Table~\ref{tab:no6}).
The average vector $\langle\tilde{\textbf{w}}\rangle$ is calculated by letting $\tilde{\textbf{w}}$ run over the above marginalization over $\Omega_m$.
The results of the analysis are summarized in Table~\ref{tab:no6} and Fig.~\ref{fig:no6}.
It is worth noting that, for both the values of $H_0$, the $w_i$ are consistent with the value $w=-1$ for $z\lesssim1.2$.
\begin{figure*}
\centering
\includegraphics[width=0.9\hsize,clip]{Combo_updated_w_z2}
\caption{The DE EoS reconstructed evolution through the redshift-binned parametrization of $w(z)$ ($1$- and $2$--$\sigma$ from the inner/darker to the outer/lighter) for the selected $H_0$. The dashed red lines mark the value $w=-1$ in the flat $\Lambda$CDM model. The darker region shows unphysical EoS, i.e. exceeding the stiff matter regime.}
\label{fig:no6}
\end{figure*}
\begin{table}
\renewcommand{\arraystretch}{1.2}
\centering
\caption{Upper part: redshift bin $z_i$, the number of GRBs $N_i$ within, and the corresponding best-fit uncorrelated $w_i$ (with $1$--$\sigma$ attached errors) for the $H_0$ of \citetalias{Planck2018} and \citetalias{2019ApJ...876...85R}. Lower part: values of the $\chi^2$ over the DOF.
\begin{tabular}{c|c|cc}
\hline\hline
\multirow{2}{*}{$z_i$} &  \multirow{2}{*}{$N_i$}   & \multicolumn{2}{c}{$w_i$}  \\
\cline{3-4}
                      &  &   P20                                  &   R19                    \\
\hline
\hline
$0.00$--$0.55$   & $16$ &   $-1.99^{+1.10}_{-2.15}$   &   $-2.46^{+1.03}_{-1.78}$\\
$0.55$--$1.18$    & $32$&   $-2.01^{+1.86}_{-1.16}$    &   $-2.44^{+2.07}_{-1.22}$\\
$1.18$--$1.74$   & $35$ &   $-1.95^{+0.37}_{-0.31}$    &   $-2.46^{+0.36}_{-0.35}$\\
$1.74$--$2.55$    & $43$&   $-2.05^{+0.33}_{-0.31}$    &   $-2.47^{+0.23}_{-0.24}$\\
$2.55$--$9.50$    & $48$&   $-2.02^{+0.21}_{-0.20}$    &   $-2.48^{+0.21}_{-0.22}$\\
\hline\hline
$\chi^2$/DOF &  &   $168.02/169$   & $161.12/169$\\
\hline
\hline
\end{tabular}
\label{tab:no6}
\end{table}

\section{Theoretical discussions over $w$ evolution}\label{sec:7}

Our aim is to use the updated sample of $174$ GRBs alone in order to get bounds over  DE evolution. Consequently we obtain tighter constraints on the cosmological parameters, compared to previous results in \citetalias{Izzo2015}, where a smaller GRB sample has been used. This first result certifies the goodness of our Combo relation compared with other calibration functions, where the increase of data points does not forcedly lead to improve the quality of fits (Luongo \& Muccino, in publication).

Moreover, we here investigate the idea of fitting cosmological models with GRBs only, i.e. without considering surveys of standard candles. This procedure resembles the analogous method proposed by \citet{2019NatAs...3..272R}, where the authors investigated how to get constraints in cosmology by using only quasars, although a scaling between the Hubble diagram of quasars with that of SNe in the common redshift range $z=0.1$--$1.44$ has been applied. Even though GRBs are not genuine standard candles, their use in cosmology leads to promising outcomes and seems to show a particular behavior constraining $\Omega_m$. In particular, from our findings we first conclude that the higher the $H_0$ value,
\begin{itemize}
    \item[{\bf 1.}] the lower the best-fit value of $\Omega_m$, and
    \item[{\bf 2.}] the better the accuracy on cosmological parameters.
\end{itemize}
However, in all our fits the outcomes have been jeopardized by the common issues related to the use of GRBs in background cosmology \citep{orlando9}. A few of these caveats are discussed below in view of our approach.

\subsection{Issues with low-$z$ GRBs}

The shortage of very low-$z$ GRBs prevents us to obtain a precise value of $H_0$ using GRBs alone. In other words, GRBs alone are unable to fix $H_0$ and to heal somehow the net tension between recent measurements on it, see e.g. \citetalias{Planck2018} and \citetalias{2019ApJ...876...85R}. This is essential for our fits, as if one has a well-bounded $H_0$ the corresponding findings become stabler, implying that $H_0$ plays the role of \emph{an initial value} that fixes the shapes of Hubble curves used in the fit \citep{orlando5,orlando6, orlando7}.

Thus, the lacking of GRBs at small redshift  severely acts on our fits. In fact, at $z\simeq0$ we have no possibility to definitively conclude that the correct cosmological model is the $\Lambda$CDM paradigm.

This means that our approach turns out to be more predictive at intermediate redshifts, i.e. as $z>1.2$. Motivated by this fact, we fixed $H_0$ in our analyses. To do so, we considered the two most recent measurements of $H_0$ by \citetalias{Planck2018} and \citetalias{2019ApJ...876...85R} in tension between each other. It is needful to fix the kinematics of our fit for the above reasons, i.e. we do not have enough GRB sources at small redshifts.

Consequently our findings in the windows of $z<1.2$ permits to take the $\Lambda$CDM suite but leaves open the possibility that DE evolves. Other probes certify slight departures from the case $w=-1$, suggesting that our wide window found using GRBs only at $z=0$ is only due to the lack of sources in this regime.

Going further with the redshift, i.e. as $z>1.2$, leads to refined limits over the DE evolution as we discuss below. Precisely, concerning the redshift-evolving DE EoS, we find that at redshift $z\lesssim1.2$ the evolution of $w$ is almost consistent with the cosmological constant case within $1$--$\sigma$ for both the $H_0$ of \citetalias{2019ApJ...876...85R} and \citetalias{Planck2018}.

\subsection{Dark energy evolution at high $z$}

A first evident deviation from the above-discussed case arises at $z\gtrsim1.5$, where departures from the $\Lambda$CDM case lies on $\gtrsim2.4$--$\sigma$ for the $H_0$ of \citetalias{Planck2018} and on $\gtrsim4$--$\sigma$ for the $H_0$ of \citetalias{2019ApJ...876...85R}. Clearly, at $z\gtrsim1.5$ a severe observational bias exists and we could indeed observe only the brightest GRBs \citep{2020arXiv201002935B} \footnote{Further analysis on this aspect, including the systematics on the entire sample, will be discussed elsewhere.}.

It is evident that the shortage of GRBs at low redshifts has a strong influence for the estimate of the DE EoS parameters $w$. Indeed, DE has its large effects in the low-$z$ Universe ($z \lesssim 0.5$), where $\approx9\%$ of the Combo-GRBs are located. Moreover, while at low-$z$ the DE drives the evolution of the Universe, at $z\geq 1$ the matter dominates. In this light, the epoch of the transition between a DE and a matter-dominated Universes can represent a new constrain for cosmological models (Izzo et al. in preparation).

The Combo-GRB sample and analysis confirm a Universe dominated by  DE at low redshifts, but seems to show a fraction of DE to act even at higher redshifts. Although our data analysis within the flat $\Lambda$CDM model is in agreement with other recent results from different probes (e.g. \citetalias{Planck2018} and \citetalias{2019ApJ...876...85R}), our results prompt a strongly negative mean value for the evolving EoS at intermediate $z$.

In this case, DE may exhibit a phantom behavior, being a likely though troublesome DE candidate. In particular, at larger $z$ current data do not allow to exclude that $w(z)$ might vary with the redshift, possibly to values $w\lesssim-1$. This finding is somehow in line with recent results from \citet{2019NatAs...3..272R}, supporting values of the DE parameter $w<-1$, i.e. a DE density increasing with time, for $z\gtrsim1.4$. Phrasing it differently, we find that the cosmological evolution is driven by DE in lieu of a pure $\Lambda$. By construction the term that dominates is, however, strongly negative, i.e. its evolution shows $w<-1$. As a desirable consequence of this scenario, such a consideration leads to a Universe that begins at Big-Bang and terminates at Big-Rip  \citep{orlando10}.
Indeed, if phantom energy continues driving the universe outward, it could literally tear space into shreds, leading to an unvoidable big rip, as pointed out in \citet{bigrip}.

\begin{figure*}
\centering
\includegraphics[width=0.9\hsize,clip]{Combo_delta_mu2}
\caption{The rebinned GRB moduli $\mu_{\rm GRB}$ (gray circles) for the $H_0$ in \citetalias{Planck2018} and \citetalias{2019ApJ...876...85R} compared with the flat $\Lambda$CDM model $\mu_{\rm th}$ (see the horizontal black line and Tab.~\ref{tab:no5}), the non-flat $\Lambda$CDM case (see the red line and Tab.~\ref{tab:no5}), and the flat model with an evolving $w(z)$ (see the dashed blue line and Tab.~\ref{tab:no6}).}
\label{fig:no12}
\end{figure*}

\vspace{1cm}

Although appealing, this scenario is unlikely and turns out to be misleading for essentially two main reasons. First, the $\Lambda$CDM model lies inside our predictions, within the $1$--$\sigma$. Hence, the simplest framework is clearly favored than any deviations. Second, it is arguable that a phantom DE is disfavored and the fact that DE stabilizes its evolution as $z$ increases indicates the need of refined initial values, as above stated. In other words, although the Combo relation improves the quality of overall analyses with GRBs, the lacking of such sources at small redshifts influences the analysis as $z>1.2$.

Summing up the cosmological constant model is still inside the GRBs predictions got from the use of the Combo relation. Evident deviations to phantom DE and/or to slight variation from $w=-1$ may be addressed noticing that at $z\simeq0$ no enough GRB sources are available.

\subsection{Does the cosmic speed up behave like phantom?}

The overall fitting procedures that we considered are clearly dependent on the choice of the data sets and bins here used. In particular, using GRBs alone likely leads to $w<-1$, while adopting combined data sets would significantly help to refine the numerical outcomes.
However, involving more than one data set based on GRBs, although it provides consistency with \citetalias{Planck2018},  the errors are suitably large, with significant discrepancy. The case is similar to what found for quasars, since treated separately, there is no tension, but as one attempts to anchor them using another distance indicator \citep{2020ApJ...888...99W}, one gets issues on errors. Thus, the indication $w<-1$ provided by GRBs can be refined but with further error propagations induced by anchored further data surveys. A relevant approach to overcome this issue is to enlarge the whole data set adopting alternative strategies of experimental analysis \citep[see, e.g.,][]{LM2020}.

\section{Final remarks}\label{sec:8}

In this work, we considered an extended GRB sample, characterized by 174 GRBs and we fixed limits over DE's evolution using the Combo relation. We extended the approach made in \citetalias{Izzo2015} and we got refined bounds over the whole matter content, i.e. $\Omega_m$ and DE density. We first worked with the $\Lambda$CDM model and then we constrained the evolution of $w(z)$ through a piecewise formulation over seven GRB redshift interval bins. We theoretically discussed our results, showing an overall agreement with the $\Lambda$CDM model that lies inside our findings. We distinguished two regimes, $z<1.2$ and $z\geq1.2$, where GRBs turn out to be more predictive. Deviations from $w=-1$ are not excluded in both the intervals. However, we proposed that at small redshifts, deviations from $w=-1$ are jeopardized by the GRB shortage. As the redshift increases, DE's abundance is much more constrained with GRBs. Thus, beyond $z\simeq1.2$ the case $w<-1$ cannot be fully excluded from our computations.

In particular, a concise summary of the fitting of $\mu_{\rm GRB}$ with various cosmological models has been  shown in Fig.~\ref{fig:no12}, where we compared the residual GRB distance moduli $\mu_{\rm GRB}$ with respect the $\mu_{\rm th}$ from the flat $\Lambda$CDM case (see horizontal black line and Tab.~\ref{tab:no5}). In Fig.~\ref{fig:no12} we also compared these figures with the distance moduli from the non-flat $\Lambda$CDM case (see the red line and Tab.~\ref{tab:no5}), and the flat wCDM with an evolving $w(z)$ (see the dashed blue line and Tab.~\ref{tab:no6}).

From all our analyses we obtained the following main results:
\begin{itemize}
\item[1)] the Combo-GRB sample confirmed a DE-dominated Universe, responsible for the observed acceleration at low redshifts. Our data analysis within the flat $\Lambda$CDM model is in agreement with other recent results from different probes (e.g. \citetalias{Planck2018} and \citetalias{2019ApJ...876...85R});
\item[2)] as stated above, it was clear that the shortage of GRBs at low redshifts has a strong influence for the estimate of  $w$. This happens since  DE has its large effects around $z \lesssim 0.5$, where $\approx9\%$ of the Combo-GRBs are located;
\item[3)] the value of w$(z)$ agreed, within 1$-\sigma$, with the standard $\Lambda$CDM model up to $z\approx1.2$. At larger $z$, current data do not allow to exclude that $w(z)$ might vary with the redshift, possibly to values $w\lesssim-1$. This finding is somehow in line with recent results from \citet{2019NatAs...3..272R}, supporting values of the DE parameter $w<-1$, i.e. a DE density increasing with time, for $z\gtrsim1.4$.
\end{itemize}

From the above considerations, we concluded by remarking the important role of GRBs as probes for the low- and high-$z$ Universe: the evidence of tight correlations between their observed properties has really become crucial to understand the evolution history of the Universe at very large redshifts, where no astrophysical probe has gone before.

Future perspectives are based on missions like the incoming SVOM and the proposed eXTP and THESEUS that will provide a robust dataset of GRBs at all redshifts ($\approx100$ at $z\gtrsim5$ over a three year mission) allowing a more detailed study of the early Universe up the re-ionization era.
The application of this forthcoming dataset to further test the Combo relation can also shed light on the physics of the X-ray afterglow and provide a physical explanation behind this GRB relation \citep[see, e.g.,][]{2017MNRAS.468..570M,2018arXiv180408652S}.

\begin{acknowledgements}
The authors are thankful to the anonymous referee for constructive suggestions that helped to significantly improve the manuscript.
This work made use of data supplied by the UK \emph{Swift} Science Data Center at the University of Leicester.
The work was partially supported by the Ministry of Education and Science of the Republic of Kazakhstan,  Grant: IRN AP08052311.
MM is supported by INFN as part of the MoonLIGHT-2 experiment in the framework of the research activities of CSN2.
LA acknowledges financial contribution from the agreement ASI-INAF n.~2017-14-H.O.
We warmly thank Michael Oliver for the support in completing this paper.
\end{acknowledgements}