A roadmap to gamma-ray bursts: new developments and applications to cosmology

Gamma-ray bursts are the most powerful explosions in the universe and are mainly placed at very large redshifts, up to $z\simeq 9$. In this short review, we first discuss gamma-ray burst classification and morphological properties. We then report the likely relations between gamma-ray bursts and other astronomical objects, such as black holes, supernovae, neutron stars, etc., discussing in detail gamma-ray burst progenitors. We classify long and short gamma-ray bursts, working out their timescales, and introduce the standard fireball model. Afterwards, we focus on direct applications of gamma-ray bursts to cosmology and underline under which conditions such sources would act as perfect standard candles if correlations between photometric and spectroscopic properties were not jeopardized by the \emph{circularity problem}. In this respect, we underline how the shortage of low-$z$ gamma-ray bursts prevents anchor gamma-ray bursts with primary distance indicators. Moreover, we analyze in detail the most adopted gamma-ray burst correlations, highlighting their main differences. We therefore show calibration techniques, comparing such treatments with non-calibration scenarios. For completeness, we discuss the physical properties of the correlation scatters and systematics occurring during experimental computations. Finally, we develop the most recent statistical methods, star formation rate and high-redshift gamma-ray burst excess and show the most recent constraints got from experimental analyses.

Gamma-ray bursts are the most powerful explosions in the universe and are mainly placed at very large redshifts, up to z 9. In this short review, we first discuss gamma-ray burst classification and morphological properties. We then report the likely relations between gamma-ray bursts and other astronomical objects, such as black holes, supernovae, neutron stars, etc., discussing in detail gamma-ray burst progenitors. We classify long and short gamma-ray bursts, working out their timescales, and introduce the standard fireball model. Afterwards, we focus on direct applications of gamma-ray bursts to cosmology and underline under which conditions such sources would act as perfect standard candles if correlations between photometric and spectroscopic properties were not jeopardized by the circularity problem. In this respect, we underline how the shortage of low-z gamma-ray bursts prevents anchor gamma-ray bursts with primary distance indicators. Moreover, we analyze in detail the most adopted gamma-ray burst correlations, highlighting their main differences. We therefore show calibration techniques, comparing such treatments with noncalibration scenarios. For completeness, we discuss the physical properties of the correlation scatters and systematics occurring during experimental computations. Finally, we develop the most recent statistical methods, star formation rate and high-redshift gamma-ray burst excess and show the most recent constraints got from experimental analyses. Gamma-ray bursts (GRBs) represent powerful extra-galactic transient that emit in γ-rays [1,2]. They are commonly associated with death of massive stars or with binary compact object mergers. As expected, due to their enormous luminosity, after the aforementioned processes, there would be the presence a newborn stellar mass black hole (BH) that provides particle accelerations and emits a relativistic collimated outflow, in the form of jets. At the same time, this new system furnishes non-thermal emissions at almost all wavelengths. The above picture lies on the standard model describing GRBs and requires isotropic energies in the range 10 44 -10 47 J, or 10 51 -10 54 erg, mostly larger than the brightest supernova (SN) emission, lying on 10 42 J, or 10 49 erg [3,4]. Thereby the need of singling out GRB progenitors is essential to disclose their fundamental properties as well as the physical conditions that permit relativistic jet to born and accelerate. Even though a clear landscape for GRB progenitor is still unclear, in view of their duration it is plausible to classify GRBs into long and short ones.
Clearly, in our Precision Cosmology epoch GRBs could open new windows 1 toward the universe description at intermediate redshifts 2 [5][6][7], i.e., much larger than SN ones [8]. Thus, several new observations have been developed, with always better accuracy, trying to standardize GRBs and to handle their emissions in analogy to SNe. In general, the most tricky challenge for cosmology is measuring distances and arguing luminosity in the cosmic scenario, understanding from astronomical emission at which distance the emitter is placed [9].
Unfortunately, this is not exactly the case of GRBs that are not standard candles, i.e., they do not provide the above requirement on distance and luminosity [10,11]. In fact, their highly variable γ-ray emission, mostly evident during the prompt phase, is thought to be associated with jet internal energy dissipation. However, the jet kinematics, among all its speed, collimation, energy, magnetization, etc., are all properties not well clarified, as well as energy dissipation mechanisms and/or shock acceleration efficiency. Hence, it is hard to relate luminosity to GRB distances as their microphysics is not well understood. Although the above caveats plague the overall GRB scenario, both short and long GRBs have relativistic outflows and share analogous properties 3 and many attempts have been spent to standardize GRBs for both clarifying their nature, internal structure and origin together with employing these objects for cosmological purposes [12,13].
In this review, we first introduce the concept of GRB and their main observable quantities. As stated above, according to time duration, we introduce the role of the t 90 duration to classify GRBs following the standard guidelines and underline the issues related to such a classification, e.g. ultra-long GRBs and X-ray flashes. To this end, we introduce the concepts of GRB progenitor, showing quantitatively the physical reasons that limit GRBs to be fully considered as genuine standard candles. However, we also emphasize how using luminosity correlations found in prompt and afterglow phases would be useful to characterize some sort of standardization technique. In this respect, we portray the main observable quantities coming from GRBs and deeply introduce the standard picture of GRB formation and evolution, dubbed the fireball model.
From all the above aspects, we expect GRBs to able to reconcile the cosmic expansion history at small and intermediate redshifts, connecting de facto late with early times, trying to open new windows toward the comprehension of cosmology. We therefore explain how GRBs serve as complementary probes to frame DE and cosmic expansion throughout the universe evolution, together with other standard candles, e.g. type Ia SN (SNeIa), baryonic acoustic oscillation (BAO), cosmic microwave background (CMB), Hubble differential data, etc. We show how to combine such data sets with GRBs and write the main features of experimental analysis for cosmological purposes. Great emphasis will be devoted to the circularity problem that essentially plagues cosmology with GRBs. Once introduced, we also underline strategies that do not take into account its role for fitting cosmological models with GRBs.
Hence, we provide how to challenge the standard cosmological model, namely the ΛCDM paradigm, with GRBs. To do so, we provide the main and evident features of cosmology with GRBs by showing how to perform error analyses, Bayesian treatments and how to handle systematics for several GRB correlations. We therefore develop model dependent and independent techniques of calibrations and report a few numerical outcomes related to GRBs, showing the most recent cosmological bounds, found with distinct procedures.
The review is split as follows. After this short introduction, in Sect. II, we classify GRBs and we report the most interesting properties, among all the classification, the progenitors and the main observable quantities coming from GRBs. In Sect. III, we work out the standard GRB model, namely the fireball paradigm. Here, we also discuss about particle and radiative processes, giving emphasis to the possible emissions coming from GRBs. In Sect. IV, we start introducing the concept of cosmology with GRBs. We thus highlight distance indicators and the concept of standard candles. In Sect. V, we explain in detail the experimental tools useful for getting Baysian analysis with GRBs. Finally, in Sect. VI, we provocatively report the concept of standardizing GRBs to permit those objects to be used in the same manner than other probes. Several issues have been raised in Sect. VII, although the likely most serious one, the circularity problem, is described in detail in Sect. VIII, where we also stress the opposite view in which one can also avoid calibration. Last but not least, we report the most recent developments of cosmology with GRBs in Sect. IX, while we conclude our journey in Sect. X with our final outlooks and perspectives of this work.

II. GRB CLASSIFICATION AND PROPERTIES
To achieve a recognized GRB classification, the strategy is to take into account the most relevant astronomical properties of such objects. Thus, as the most prominent GRB component is represented by the prompt γ-ray emission, it is straightforward to use it to define GRB classes based on similarity criteria.
The prompt γ-ray emission is characterized by highly-variable and multi-peaked light curves composed of either overlapping or distinct pulses with variable duration. The duration of these pulses spreads within a wide time range. Since the duration is not fixed a priori, it is natural to wonder whether one can arbitrarily define a time in which the above measures can be got. Hence, it is a consolidate convention to take the total burst duration in a time interval, whereas the dotted one is for SGRBs. The dashed horizontal lines mark mean HRs for both classes. The solid t 90 histogram represents the raw data whilst the dotted one shows the error-convolved data. Credit from Ref. [3].
dubbed t 90 , evaluated in the observer frame over which the 90%, from 5% to 95%, of the total background-subtracted counts are experimentally detectable.
In view of such a property, one gets a plausible classification, as we report below.

A. Classification: short and long GRBs
The light curve analysis of the first BATSE GRB catalogue evidenced a clear bimodal distribution of the t 90 duration, separated at roughly 2 s, and in the hardness ratio (HR), namely the ratio of the total counts of the hard 100-300 keV energy band over the softer 50-100 keV band [3,14,15].
The significance of such a classification scheme has been strengthened with the full 2704 GRBs detected by BATSE and later GRB missions, providing strong evidence for two GRB progenitor channels (see e.g. Fig. 1). However, a significant overlap in the distributions of SGRBs and LGRBs suggests that a more robust classification scheme based on physical properties is still missing.

B. Intermediate GRBs?
We ended up the previous subsection, asserting the need of a more robust classification order. This scheme is veritably challenged by the existence of an intermediate class of SGRBs with extended emission (SGRBEEs), characterized by an initial short duration and spectrally-hard γ-ray pulse, followed by a softer emission lasting up to tens of seconds [15,16]. Depending on the sensitivity and energy range of the GRB alert instrument and based on the above classification scheme, a SGRBEE could be classified as short or long. A GRB detector with low sensitivity at low-energy ranges in γ-ray could detect only the initial hard part of the burst (resulting in an SGRB), whereas a GRB detector with a higher sensitivity extending down to lower energies could detect also the softer extended emission (falling in the LGRB class).
A possible explanation to the origin of this extended emission involves a highly magnetized neutron star (NS) dipole spin-down emission (see Ref. [17] and Secs. II D 2 and III A 3).
C. Ultra-long GRBs and X-ray flashes Further, the detection of rare events characterized by extremely long-lived prompt emissions lasting 10 3 s, named ultra-long GRBs (ULGRBs), represents an additional classification threat, since it is still unclear whether ULGRBs represent a distinct class of LGRBs [18], or whether they are the high-end tail of the t 90 distribution of LGRBs [19].
Finally, it has been reported the existence of extragalactic transient X-ray sources, dubbed X-ray flashes (XRFs), with spatial distribution, spectral and temporal characteristics similar to LGRBs [20,21]. The distinguishing properties of XRFs are a) their observed prompt emission spectrum that peaks at energies which are an order of magnitude lower than those of standard LGRBs; b) their time integrated flux in the 2-30 keV X-ray band greater than that in the 30-400 keVγ-ray band.
In view of these hazy results, classifying GRBs through t 90 and HR criteria only turns out to be puzzling, since the measured t 90 varies with energy range. Thus, the definition of a novel GRB classification scheme requires multiwavelength criteria to better understand the physical properties behind the GRB emission.
In this respect, attempts to re-categorize GRBs from the popular long/short classes have been made in Ref. [22], introducing alternative classes of Type I and Type II GRBs. According to this scheme, the Type I class includes short/hard GRBs and SGRBEEs with no SN association, typically found in regions of their host galaxy with low star formation, and very likely originating in compact star mergers (see details in the next Sec. II D). On the other hand, the Type II class includes long and relatively soft GRBs with SN association, usually found in star forming regions within irregular host galaxies, and thus associated with young stellar populations and likely originating in the core-collapses of massive stars (again, see details in the next Sec. II D).
Though the above scheme seems to be promising, further research on this issue is still ongoing. Therefore, for historical reasons in Sec. II D we pursue the description of the progenitor systems keeping the bimodal classification in LGRBs and SGRBs.

D. Progenitors and open questions
Beside the above discussion, the working definition of LGRBs and SGRBs suggests the existence of two different progenitor channels. Summarizing,

I:
LGRBs could arise from the core-collapse of a massive star or collapsar [23], II: SGRBs could originate from the binary neutron star-black hole (NS-BH) or NS-NS mergers [24].
The huge observed isotropic equivalent energy release of ∼ 10 49 -10 55 erg implies that: for LGRBs, up to ∼ 10 M are converted into radiation during the prompt emission duration of ∼ 100 s, whereas for SGRBs up to ∼ 1 M are converted into radiation within ∼ 1 s [25]. The energy reservoir and the efficiency of the involved physical processes in producing the emitted energy represent a stringent requirement, specially for LGRBs. 4 .
The commonly called jets substantially alleviate this issue by reducing the GRB energy release by jet's correction factor f = 1 − cos θ. Jets can be thought, in an oversimplified picture, as outflows of relativistic matter ejected into a double-cone structure of opening angle θ. In general the jet correction is poorly constrained because requires very challenging measurements of θ and the observer viewing angle relative to the jet axis. This makes troublesome to distinguish between geometric and dynamical effects. Indeed, very soft GRBs could be bursts viewed off-axis, whereas low luminosity GRBs may be the result of large jet opening angles [15].
Measurements of θ can be obtained by the predicted signature of the achromatic jet breaks, observable in the afterglow light curve at all frequencies. This feature can be explained by the dynamics of the GRB ejecta as follows. At the beginning, at high values of the bulk Lorentz Γ factor 5 , the ejecta is narrowly beamed into the jets while its Lorentz factor is Γ −1 < θ and, regardless the hydrodynamic evolution, a GRB is observed only from a small fraction of the ejecta [15]. As the ejecta decelerates, Γ decreases below θ −1 , the beaming angle becomes larger, and a larger portion of the ejecta becomes observable. Continuous deceleration leads to the point that the entire surface of the jet is observable and the jet begins to spread sideways, producing a break in the light curve across the entire afterglow spectrum [27,28]. The sharpness of this break and the change in the afterglow decay rate depend on how long the jet remains collimated and on the jet radial density profile and energy distribution [29,30]. The time of the jet break is related to the jet opening angle, the bulk Lorentz factor and the density of the circumburst medium (CBM). The above description has two effects: • an "on-axis" observer detects the prompt emission and then, as the jet decelerates, the afterglow emission and finally a jet-break due to the faster spreading of the emitted radiation; • an "off-axis" observer cannot detect the prompt emission, but detects an orphan afterglow, namely an afterglow without a preceding GRB.
In the pre-Swift era, simultaneous breaks in the optical and near-infrared (NIR) afterglow light curves were frequently interpreted as jet breaks. Nevertheless, the improved temporal and spectral coverage of GRB afterglows, specially in X-rays by Swift, have revealed within the first few hours after the prompt emission a complex structure made of flares, plateaus and chromatic breaks [31][32][33][34]. The detected achromatic breaks are observed in a few cases. The absence of jet-break signatures in most GRB afterglows has been interpreted as due to the over-simplified assumption homogenous jets with sharp edges, whereas more complex models now include structured jets that produce several chromatic jet-breaks, or much smoother breaks, or jets that can keep their structure for longer than previously thought making difficult to detect breaks without a wide temporal coverage [30].
Besides the jet modeling issue, any GRB model has to deal with features like very luminous X-ray flares occurring up to a few 10 4 s after the GRB trigger and with shape and spectra similar to those flares observed during prompt emission and extended plateau phases that last for a few hours during the early afterglow evolution [33,34]. Both features imply an extended central-engine activity with a continuous source of energy injection lasting the above 10 4 s. In the standard picture, such long-lived energy injection requires the accretion of a significant mass onto the central BH via very large (∼ 1 M ) and low-viscosity (α < 10 −2 ) accretion disk formed at the core collapse time, or via fall-back material continuously replenishing the accretion disk [35].

The LGRB-supernova connection
The possible connection between LGRB and massive progenitor stars has been speculated long before the first afterglow detection [23,36]. The first observational evidence came with the association between the broad line Type Ic (Ic-BL) SN 1998bw and the low-luminous LGRB 980425 at z = 0.0085 and with lack of an optical afterglow [37]. Later on, this association was also confirmed between the Type Ic-BL SN 2003dh, temporally and spatially coincident with the standard more luminous long GRB 030329 at z = 0.1685, with an optical afterglow light curve comparable with other cosmological GRBs [38].
The launch of Swift has increased the sample of GRB-SN pairs, both spectroscopic, at z 0.5 and most of them with isotropic-equivalent γ-ray energies E iso < 10 49 erg, and photometric, in the form of SN bumps appearing in the optical afterglows 10 − 30 days 6 after the GRB, at z 0.5 and E iso ≈ 10 51 − 10 52 erg [36]. Most of the GRB-SN pairs belong to this second kind, very likely at the hand of a selection effect: the more common low-luminosity LGRBs per unit volume are not detectable at high redshift, whereas luminous LGRBs, with higher detectability at high redshift are observed from a larger volumetric area [39].
SNe Ic associated to some long GRBs are characterized by no hydrogen (H) and no weak helium (He) lines [36]. Their occurrence close to star-forming regions offer very strong evidences that long GRBs could be associated with massive star death [36]. In this regard, the best progenitor candidates are the Wolf-Rayet stars, very massive stars with an hydrogen envelope largely depleted, endowed with a fast rotation [23,40]. Within the collapsar model, very massive stars are able to fuse material in their centers all the way to iron (Fe). At this point they cannot continue to generate energy through fusion and collapse mechanisms forming a BH. Matter from the star around the core rains down towards the center and swirls into a high-density accretion disk. In this picture, the core carries high angular momentum to form a pair of relativistic jetsout along the rotational axis where the matter density is much lower than in the accretion disk. Jets propagate through the stellar envelope at velocities approaching the speed of light, creating a relativistic shock wave at the front [15,41]. If the star is not surrounded by a thick, diffuse H envelope, the leading shock accelerates as the stellar matter density decreases. Thus, by the time it reaches the star surface, Γ ≥ 100 is attained and the energy is released in the form of γ-ray photons [15,41].
The collapsar model attempts to explain the time structure of GRBs prompt emission, through the modulation of the jets by their interaction with the surrounding medium, which could produce the variable Lorentz factor needful for internal shock occurrence [23]. As the relativistic jet propagation through the stellar envelope of a collapsing star proceeds, its collimation was shown to occur analytically and numerically [15]. Another prediction of this model is the prolonged activity of the central engine which can potentially contribute to the GRB afterglow [23,40,41]. This occurs because the jet and the disk are inefficient at ejecting all the matter in the equatorial plane of the pre-collapse star and some continues to fall back and accrete [23,40,41].
The SNe associated with LGRBs appear to belong to the bright tail of type Ic SNe and can be considered as a "subclass" of SNe Ic, alternatively addressed as hypernovae, in order to emphasize the extremely high energy involved in these explosions. Remarkably, the SNe associated with both low-and high-luminous (XRFs and normal LGRBs, respectively) share very similar spectra and their peak luminosities span only 2 orders of magnitude, whereas the associated GRBs isotropic luminosities span six orders of magnitude [36]. Another distinctive feature of the GRB-SN pairs is the high photospheric expansion velocity, up to 0.1c [36]. In this scenario, one has to fit also the class of ULGRBs. The spectroscopic detection of the SN 2011kl coincident with the ULGRB 111209A [42] favors a common core-collapse origin for LGRBs and ULGRBs. This SN exhibited a peculiar, very blue and featureless spectral shape, which was unlike other SNe Ic associated with LGRBs, but more alike the newly discovered class of superluminous SNe [43]. Other ULGRBs have either been too far or too dust-extincted to secure any detection of an underlying SN, whereas other cases evidenced the indicative flattening from a rising SN in their optical and NIR light curves at 10 − 20 days after the GRB trigger [36].
In this picture, however, exceptions to the LGRB-SN association have been found from deep optical observations in two nearby bursts, GRB 060505 and GRB 060614, for which the hypothetical accompanying SN would have been a hundred times fainter than SN 1998bw [44][45][46].
To conclude, ULGRBs and SNless LGRBs give evidences for the existence of further progenitor channels for LGRBs.

SGRBs, macronovae and gravitational waves
The Swift satellite has enabled rapid and precise localisations and an increase in the number of X-ray and optical afterglow detection of both LGRBs and SGRBs. However, SGRBs have less luminous afterglows than those of LGRBs and this fact makes difficult to obtain optical spectra and a precise burst location to plan optical follow-up to search for host galaxy associations. The lack of any associated core-collapse SNe, the typically large offsets of the GRB position with respect to galaxy center, and the frequent association with galaxies with no ongoing star formation, provide evidences in support of a compact binary merger progenitor scenario [47].
The proposed progenitors for SGRBs are NS-NS and/or NS-BH binary mergers [24,[48][49][50]. These mergers take place as binary orbits decay due to gravitational radiation emission [51]. A merger releases 5 × 10 53 erg, but most of this energy is due to low energy neutrinos and gravitational waves. So, there is enough energy available to produce a GRB, notwithstanding how a merger generates the relativistic wind required to power a burst is still object of speculations and not-well understood. It has been argued that about one over thousand of these neutrinos annihilates and produces pairs that in turn produces γ-rays via νν → e + e − → γγ, but it has been pointed out that a large fraction of the neutrinos would be swallowed by the newly-born BH [15].
A further confirmation to the binary merger scenario consists in the detection of the so-called macronova (MN). The MN emission originates in NS-BH or NS-NS mergers from a fast-moving, rapidly-cooling ejected debris of neutron-rich radioactive species that decay to form transient emission and create atomic nuclei heavier than iron through neutron capture process, named r-process [52]. The opacities of these produced heavy elements lead to a dim MN emission, requiring deep follow-up observations down to NIR bands. The first indication of a MN, in the form of a re-brightening detected approximately 9 days after the GRB trigger, has been obtained by extensive follow-up of the SGRB 130603B, one of the nearest and brightest SGRBs ever detected [53]. An MN emission accompanies also the nearest SGRB ever detected, SGRB 160821B [54,55], and the recently detected SGRB 200522A [56]. For a list of other MN emissions, see Ref. [57].
In the binary merger scenario, SGRBs are expected to be significant sources of gravitational waves (GWs). The smoking gun occurred on 17 August 2017, when the Advanced LIGO and Virgo detectors observed the event GW 170817, unambiguously detected in spatial and temporal coincident with the SGRB 170817A independently measured by the Fermi Gamma-ray Burst Monitor, and the Anti-Coincidence Shield for the Spectrometer for the International Gamma-Ray Astrophysics Laboratory [58] 7 . As a further confirmation on the nature of the progenitor system of SGRB 170817A, an intense observing campaign from radio to X-ray wavelengths over the following days and weeks after the trigger led to the spectroscopic identification of a MN emission, dubbed AT 2017gfo [59].
The observation of SGRB, GW and MN emission has improved our understanding of the physical properties related to the binary merger, such as the mass of the compact object, the ejected mass, and the details of the CBM surrounding the merger site.

E. Observable quantities from GRBs
Understanding GRB physics passes through the experimental evidence of the energy that can be collected from detectors. In particular, we can start discussing about GRB prompt emission. It is typically observed in the hard-X (above ∼ 5 keV) and γ-ray energy domain.
The operative duration of the prompt emission is due to the previously defined t 90 . Within this time interval, and also within any sub-interval with enough photons to perform a significant analysis 8 , the observed spectral energy distribution (SED) of GRBs is non-thermal, and it is best fitted by a phenomenological model composed of a smoothly joined broken power-law called Band model [60] (see Fig. 2). Its functional form is where typical power-law index values are −1.5 α 0 (with an average α −1) and −2.5 ≤ β ≤ −1.5 (with an average β −2), while the peak energy at the maximum of the of the E 2 N E (or EF E ) spectrum lies within 100 keV≤ E obs p few MeV (with an average of E obs p 200 keV). Finally, K is the normalization constant with units of photons cm −2 s −1 keV −1 . In some cases the SED is also best fitted by a power-law model 9 or by a power-law plus an exponential cutoff. However, these models are purely mathematical, i.e. not yet physically linked to GRB intrinsic properties. Hence, fitting data with them does not provide any insight about the emission physical origin, but may be useful for the classification scheme of GRBs and for comparing the fitted results with the predictions of different theoretical models.
In the recent years, with a much broader spectral coverage enabled by detectors such as Fermi, evidences for more complicated broad band spectra fitted by a combined Band+thermal model, have been found in an increasing number of bursts [61][62][63][64], where the peak of the thermal component is always observed below E obs p . However, the search of the best-fit model in describing GRB prompt emission spectra depends on the analysis method. Typically, a significant spectral analysis is performed when enough photons are collected. For weak bursts, only time-integrated spectral analyses can be done and this implies that important time-dependent features may be lost or averaged, leading to a wrong theoretical interpretation. Another issue is that the chosen spectral model is convolved with the detector response and, because of the non-linearity of the detector response matrix, this procedure cannot be inverted. Therefore, two different models can equally provide a similar minimal difference between the model and the detected counts spectrum and lead to different theoretical interpretations.
From the fit of the time-integrated prompt emission spectrum, one can get the flux F (in units of erg cm −2 s −1 ) on a detector energy bandpass E min -E max as where κ is a constant, commonly used to convert the energy, expressed in keV, to erg. 7 The GW signal, originating from the shell elliptical galaxy NGC 4993, had a duration of ∼ 100 s. By the characteristics in intensity and frequency, GW 170817 has been unambiguously associated with the inspiraling of a binary NS-SN merger of total mass 2.82 +0.09 −0.47 M , which is consistent with the masses of all known binary NS systems. 8 Typically dubbed time-integrated and time-resolved analyses, respectively. 9 However, in this case the spectral break is very likely below or above the detector bandpass. To compute the total energy emitted by a GRB in all wavelengths, a bolometric spectrum is needed. However, the GRB prompt emission triggers γ-rays detectors in a given energy bandpass, therefore a limited part of the spectrum is available, instead of a bolometric one. Moreover, GRBs are cosmological sources spread over a wide redshift range, so for GRBs observed by the same detector, the measured energy range corresponds to different energy bands in the cosmological rest frame of the sources.
To standardize all GRBs, fluxes are computed in the fixed rest-frame band 1-10 4 keV, which is a range larger than that of most of the γ-ray detectors. The "bolometric" time-integrated flux is then given by and the total isotropically-emitted energy and luminosity are, respectively where the factor (1 + z) −1 corrects the t 90 duration from the observer frame to the GRB cosmological rest-frame. In a similar way, the peak luminosity L p , computed from the observed peak flux F p within the time interval of 1 s around the most intense peak of the burst light curve and in the rest frame 30-10 4 keV energy band 10 is given by The luminosity distance d L depends upon the cosmological models adopted as backgrounds and can be related to the continuity equation recast as that relates the total energy density ρ and pressure P to the barotropic factor ω(z) ≡ P/ρ of a given cosmological model. For a two component flat background cosmology composed of standard pressure-less matter with ω = 0 and a generic DE component with ω(z) (dubbed generically ωCDM), the luminosity distance is then given by 11 where H 0 is the Hubble constant, Ω m and Ω x are the cosmological density parameters of matter and DE, respectively, and f x (z) is given by For the concordance paradigm, namely the ΛCDM model, the DE equation of state is w(z) ≡ −1 corresponding to a cosmological constant Λ. Thus, f x ≡ 1 and Ω x ≡ Ω Λ . In the following, the choice w(z) ≡ −1 is adopted, unless otherwise specified. The above isotropic energy output can be corrected for the beaming (see Sec. III), once the jet opening angle θ is known, leading to beam corrected energy It is important to stress that the prompt emission is not limited to the γ-rays and that, differently from the afterglow emission starting ∼ 100 s after the GRB trigger, current information in other energy bands is extremely difficult to observe without fast triggering. Observations at lower energies (optical and X-rays) have been enabled only for GRBs with a precursor or a very long prompt emission duration, which gave the possibility of performing fast pointing to the source during the prompt phase [66].
Regarding the GeV energy domain, a delayed (with respect to the trigger), long lived emission ( 10 2 s), and separate lightcurve [67] with a decaying luminosity as a power law in time, L GeV ∝ t −1.2 has been observed [67]. These distinctive features point towards a separate origin of the GeV with respect to the lower energy photons.
After ∼ 100 s since the trigger the prompt emission starts to decay in flux and, in many cases, this feature is caught by X-ray detectors Swift-XRT within 0.3-10 keV energy band. In general X-ray afterglow light curves show complex behaviors [15] consisting of (see Fig. 3): (1) an early steep decay, interpreted as the tail of the prompt emission at large angles, followed by a very shallow decay, called the plateau, usually accompanied by spectral parameter variations, and (2) a final decay, less steep than the first one.
The X-ray afterglow is also characterized by the presence of a flaring activity [15]. The observed behavior of these flares, the rapid rise and exponential decay together with a fluence comparable in some cases to the prompt emission, points out that the same mechanism for the prompt emission is responsible for the flaring activity [15]. Concluding, as above already stressed, these X-ray afterglow features are important to understand the nature of GRB progenitors.

Timescales and characteristic energy as observable signature of GRBs
There are other GRB observable quantities often employed in the literature, e.g., to construct GRB correlations (see Sec. VI A). They span from characteristic energies to timescales measured in several wavelengths. More specifically, a selection of them is summarized below: • t b , the time at which the late X-ray afterglow power-law decline suddenly steepens due to the slowing down of the jet until the relativistic beaming angle roughly equals the jet opening angle θ.
• τ lag , the time lag is computed as the difference of arrival times to the observer of the high energy photons (100 − 300 keV) and low energy photons (25 − 50 keV). 12 • t X , the rest-frame time, defined by a broken power-law fit of the X-ray afterglow light curve, at which a late power-law decay after the plateau phase is established.
• τ , the rest-frame time marking the end of the plateau phase, defined from a fit of the X-ray afterglow with a smooth function given in Ref. [68].
• F X and F 0 are the observed X-ray fluxes respective to t X and τ , whereas the corresponding rest-frame 0.3-10 keV luminosities L X and L 0 are computes as follows 13 • V , the variability of the GRB light curve. It is computed by taking the difference between the observed light curve and its smoothed version, squaring this difference, summing these squared differences over time intervals, and appropriately normalizing the resulting sum.

III. THEORY OF GRB PROGENITORS
GRBs require progenitor systems able to guarantee enough energy for their powerful explosions to occur and emission mechanisms that can explain the above discussed spectral features. Although essential to better understand 12 The time lag is historically computed in these energy bands which are the BATSE energy channels 3 and 1, respectively. 13 As a convention, the X-ray luminosities are computed in a rest-frame energy band with similar extrema with respect to the observed one; with this prescription their expression are simple, as portrayed in Eq. (11).
where we used the fact that X-ray data are observed by the Swift-XRT in the 0.3-10 keV energy band and the SED is in general a power-law spectrum with N the physics of GRBs, neither a clear evidence for consolidate classes of suitable progenitors, nor a definitive GRB model have been yet established, as above stressed. However, observations, in the form of GRB spectra and light curves (see Sec. II E) and correlations between observable quantities (see Secs. VI A and VII), enhanced our comprehension of these phenomena and led to a general agreement on a few aspects below listed [69]: -GRB progenitors harbor a BH 14 which acts as a central engine powering the GRB emission.
-The burst energy must be gravitational and it is released in a very short time and from a compact region.
-Substantial part of this energy is converted into kinetic energy and a relativistic jetted outflow is formed.
-The acceleration process and the role played by magnetic fields are still unclear.
-The dissipation of part of the kinetic energy produces the observed prompt emission.
-The thermal emission of the prompt emission, may be the relic of the photons emitted during the initial explosion, whose energy has not been converted into kinetic form.
-Afterward, relativistic jets interact with the CBM, gradual energy conversion occurs and the afterglow (from X-ray down to radio) is produced.
The observed spectra have a considerable amount of γ-ray photons. Photons with high energy E 1 annihilate with those at a low energy E 2 and produce e + e − pairs if √ E 1 E 2 m e c 2 (up to an angular factor), where m e is the electron mass. If GRBs were not relativistic sources, the observed light curve variability time scale of δt ≈ 10 ms would imply that their emission would originate from a very compact region not larger than R = cδt ≈ 3000 km. For typical values of the luminosity distance d L ≈ 3Gpc ≈ 10 22 cm and fluence S ≈ 10 −7 erg cm −2 (energy at the detector per unit area) of GRBs, the opacity for pair creation is enourmous and it is given by [69] τ where f e ± is the fraction of photons with energies sufficient to produce pairs and σ T is the Thomson cross-section. Such a large optical depth would imply that that the source must be optically thick leading to a thermal spectrum. On the contrary, observations indicate that GRB spectra are typically non-thermal, pointing to the opposite conclusion that their source must be optically thin. This issue is named compactness problem [69]. However, the problem is only apparent, once relativistic effects are taken into account. In fact, the causality limit of a source moving relativistically with bulk Lorentz factor Γ 1 towards the observer is R ≤ Γ 2 cδt. Consequently, the observed photons are blue-shifted and their energy at the source is lower by a factor ≈ 1/Γ, which may be insufficient for pair production. This leads to a decrease in the opacity, by a factor Γ −2(β+1) , where the β is the high-energy power-law index of photon spectrum of the burst. For Γ 100 one obtain the optically thin condition of the source. Ultra-relativistic expansion of GRBs is unprecedented in astrophysics. There are indications that relativistic jets in active galactic nuclei have Γ ∼ 2-10, but some GRBs have Γ 100. These large expansion velocities in GRB outflows find confirmations from the radio scintillation observed in their afterglows, and also from the apparent observation of self-absorption in the radio spectrum of the afterglow, where it is possible to obtain independent estimates of the dimensions of the afterglow relic [15].

A. The fireball model
The GRB standard model considers a homogeneous fireball, [69]. For a pure radiation fireball, a large fraction of the initial energy released by the newly-formed BH is converted directly into photons. Close to the BH, at a radius r 0 larger than the Schwarzschild radius, R S = 2GM/c 2 , the photon temperature is where a is the radiation constant, and the luminosity L and the radius r 0 are expressed, respectively, as L 52 = L/10 52 erg/s and r 0,7 = r 0 /10 7 cm. In the following, to understand the order of magnitude of the key physical parameters characterizing GRBs, we use the notation Q x = Q/10 x , where the quantity Q is given in cgs units. The temperature T 0 is above the threshold for pair production, hence a large number of e ± pairs are created via photon-photon interactions, leading to a fully thermalized pairs-photons plasma with the opacity in Eq. (12). 15 GRB luminosities are many orders of magnitude above the Eddington luminosity, L E = 1.25×10 38 (M/M ) erg s −1 , therefore, the radiation pressure is much larger than self gravity and the fireball expands under its own pressure up to Γ ≈ 10 2 -10 3 [48,49]. Since the final kinetic energy cannot exceed the initial explosion energy E tot , the maximum attainable Lorentz factor is defined as Γ max = E tot /M c 2 and depends upon the amount of baryons (baryon load) of rest mass M within the fireball [69].

Photon-dominated scenario
The simplest scenario considers a photon-dominated expanding shell of width δr "instantaneously" releasing its energy. From here on, prime symbols indicate quantities measured in the comoving frame of the shell, namely from an observer within it. On the other hand, r is the radial coordinate of the laboratory frame, a frame outside the shell where the observer is sitting on the central engine.
Enforcing energy and entropy conservation laws, the shell keeps accelerating up to Γ max η, which is attained at at the so-called dissipation radius r s ∼ ηr 0 ; beyond it, most of the internal energy of the shell has been converted into the kinetic one, so the flow no longer accelerates and it coasts. Thus, the fireball obeys the following scaling laws of the shell comoving temperature, Lorentz factor and comoving volume, respectively from which follows that as the shell accelerates (as Γ increases with r), its internal energy drops (as T decreases with r). Finally, the evolution the observed temperature is given by

Internal shock scenario and photospheric emission
In the case of LGRBs, the progenitor continuously emits energy at a rate L, over a longer duration t r 0 /c, and ejects mass at a rateṀ = L/(ηc 2 ). In this case, the scaling laws for the instantaneous release are still valid, provided that E is replaced by L and M byṀ , and a further equation for the mass conservation of the baryons (within the spherical symmetry assumption) is required [72] where n p (r) is the comoving number density of baryons and m p is the proton mass. For longer activity of the inner engine, fluctuations in the energy emission rate would result in the propagation of independent shells, each of them with analogous thickness r 0 and dynamics. For two consecutive shells with a difference in their Lorentz factors δΓ ∼ η or velocities δv ≈ c/(2η 2 ), collisions become possible after a typical time t col = r 0 /δv and an observer frame radius [73] r col = vt col ct col 2η 2 r 0 .
Above r col , which is a factor η larger than r s , collisions occur, dissipate the kinetic energy and convert it into the observed radiation [74,75]. The advantages of the internal shock scenario are listed below.
1. Light curve variability. The time delay between the photons produced by the collisions and a photon emitted from the central engine towards the observer, i.e., δt ob r col /(2η 2 c) ∼ r 0 /c, is similar to the central engine variability and can explain the observed variability ( 1 ms).
2. Particle acceleration. Shell collisions generate internal shock waves, which can accelerate particles to high energies via Fermi mechanism and produce γ-rays.
However, the internal shock scenario manifests some drawbacks.
1. Efficiency. From the energy and momentum conservation, the kinetic energy dissipation is highly efficient only if two shells have masses m 1 m 2 and Lorentz factors Γ 1 /Γ 2 1. The average over several collisions leads to a low global efficiency of 1-10% [81,82], which contrasts to the much higher efficiency ∼ 50% inferred from afterglow measurements [31,33]. Higher efficiency up to the ∼ 60% can be attained by considering larger contrasts Γ 1 /Γ 2 10 [82]. However, these Lorentz factor contrasts unlikely occur within the traditional collapsar or the merger scenarios.
2. Observed spectra. This model does not explain the observed spectra and needs further assumptions on how the dissipated energy produces photons (i.e., involving standard radiative processes such as synchrotron emission or Compton scattering).

Magnetized (or Poynting-flux dominated) outflows
The Poynting-flux dominated model speculates that the gravitational energy produces very strong magnetic fields, which may be crucial in the jet formation of GRBs, similarly to the Active Galactic Nuclei (AGN), where magnetic energy is converted into particle acceleration via Blandford-Znajek [83] or Blandford-Payne [84] mechanisms. The idea behind this model is that the collapse of a white dwarf (WD) induced by accretion from a massive star, the core collapse of a massive star, or NS merger does not form immediately a BH, but rather a rapidly-spinning (with a period of ∼ 1 ms) and highly-magnetised NS (with a magnetic field B 10 15 G) NS, known as magnetar [70]. The maximum amount of magnetic energy that can be stored is ∼ 2 × 10 52 erg and it can be extracted in a short timescale of ∼ 10 s and drives a jet along the polar axis of the NS powering the prompt emission [71]. The decay of rotational or magnetic energy may continue to power late time flaring or afterglow emission. The dipole radiation naturally produces a plateau phase up to the dipole spin-down time scale [15].
In this model, the magnetic field is essentially toroidal (i.e., B ⊥ β) and its polarity in the flow changes on small scale defined by the light cylinder in the central engine. The total luminosity is given by L = L k + L M , where L k = ΓṀ c 2 is the kinetic part and L M = 4πr 2 c[B 2 /(4π)] is the magnetic part [85,86]. The key parameter is the magnetization σ ≡ L M /L k = B 2 /(4πΓ 2 nm p c 2 ), which plays a similar role to the baryon load in the classical model and defines the maximum attainable Lorentz factor Γ max ≈ σ 3/2 , whereas, during the acceleration phase, one gets Γ(r) ∝ r 1/3 [85,86].
In this model the rapid variability observed in GRBs and the low efficiency in dissipating the kinetic energy via shock waves in highly magnetized plasmas are still open issues. Recent recipes suggest that central engine variability leads to the ejection of magnetized plasma shells which expand due to internal magnetic pressure gradient and collide at a distance r col . The ordered magnetic field lines of the ejecta get distorted and fast reconnection occurs. The induced relativistic turbulence may be able to overcome the low efficiency difficulty of the classical internal shock scenario [87].

Particle acceleration
To produce the non-thermal GRB spectra, part of the kinetic energy needs to be dissipated and used to increase the random motion of the outflow particles and/or accelerates some fraction of them to a non-thermal distribution. Once accelerated, these high-energy particles emit non-thermal photons. The most widely proposed particle acceleration mechanism within the internal shock scenario is the Fermi mechanism [88]. In this process, the accelerated particles cross the shock fronts, and during each crossing their energy increases at a constant rate ∆E/E ∼ 1. The accelerated particles have a power-law energy distribution N (E) ∝ E −δ with index δ ≈ 2.0-2.4 [89].
Dissipation mechanisms in magnetized outflows have been longly discussed (see Ref. [90] and references therein). Further, particles may also be accelerated via Fermi mechanism in shock waves, but it has been pointed out that in highly magnetized plasma this process may be inefficient [91].
Recently, with mounting evidences of thermal components in GRB spectra [63,99], the photospheric model has acquired growing relevance [61,100,101]. This model is not in constrast with and has to be considered complementary to the synchrotron+SSC emission, which originates from a different region of the outflow.

Synchrotron emission
The synchrotron emission has been extensively studied for interpreting first the non-thermal emission in AGNs and then the GRB afterglow emission [102]. Regarding the explanation of the GRB prompt emission spectra [74,[93][94][95]103], the synchrotron emission model has several advantages: 1) requires energetic particles and strong magnetic fields, both expected in shock waves; 2) has a broad-band spectrum with characteristic peak, associated with the observed peak energy; 3) for typical parameters, energetic electrons radiate nearly 100% of their energy.
A source at redshift z, expanding with Γ = 1 − β 2 and at an angle θ with respect to the observer, emits photons which are seen with a Doppler boost D = [Γ(1 − β cos θ)] −1 . In the comoving frame, electrons move in a magnetic field B and, thus, have random Lorentz factor γ e . Their typical energy is [102] Typical GRB peak energies ε ob ≈ 200 keV require strong magnetic fields and very energetic electrons, both feasible for Poynting flux-dominated outflows or photon-dominated outflows where strong magnetic fields may be generated via Weibel instabilities [104]. 16 On the other hand, strong magnetic fields imply the comoving cooling time of the electrons to be t cool t d ∼ R/(Γc). Thus the expected synchrotron spectrum below the peak energy would be F ν ∝ ν −1/2 (or N E ∝ E −3/2 ) [106,107], which is inconsistent with the average low energy spectral slope α = −1 (see Fig. 2) and, hence, the value α = −3/2 is called "synchrotron line of death". To overcome this problem, electrons must cool slowly, leading to a spectrum below the peak given by F ν ∝ ν 1/3 (or N E ∝ E −2/3 ), which is roughly consistent with the observations. However, the condition t cool t d leads to high values of γ e , whereas B would be very low and, in order to explain the observed flux, the electron energy would be several orders of magnitude higher than that stored in the magnetic field [108]. To overcome this, the inverse Compton contribution has to be significant, producing ∼ TeV emission. To avoid a substantial increase of the total energy budget, the emission radius should be R 10 17 cm, but cannot explain the rapid variability observed [108]. Suggested modifications (and drawbacks) to the synchrotron scenario can be found in the literature [25,98,[109][110][111][112][113].

Photospheric emission
For r ph r s , a large fraction of the kinetic energy is dissipated below the photosphere [114]. The produced non-thermal photons cannot directly escape and are advected with the flow until the transparency. Within the flow, multiple Compton scatterings occur and modify the synchrotron spectrum of the heated electrons, which rapidly cool, mainly by IC scattering. The electron distribution becomes quasi-Maxwellian, with a temperature determined by the balance between heating (external and by direct Compton scattering of energetic photons), and cooling (adiabatic and radiative) [115]. Finally, the photon field is modified by the scattering from the quasi-Maxwellian electron distribution [114].
Further, the thermal photons of the fireball contribute as seed for IC scattering, hence the non-thermal electrons, heated by energy dissipation below the photosphere, rapidly cool and reach a quasi-steady state distribution [101]. The result is a two temperature plasma, with electron temperature T e > T ph . If dissipation processes occur at intermediate optical depth τ ∼few-few tens, the resulting spectrum is: 1) similar to the Rayleigh-Jeans part of the thermal spectrum, for T < T ph ; 3) an exponential cutoff, for for T > T e .
The spectral slope obtained in the above 2) is similar to the high energy spectral slope in GRB spectra, β ∼ −2, thus, it could be concluded that E p is associated with T ph . However, recently Fermi data has shown thermal peaks at lower energies than E p , which points rather to the more natural interpretation that the thermal peak is associated with T ph and suggests that E p may be associated with T e or with the synchrotron emission. Moreover, if dissipation occurs at τ 10 2 , the resulting spectra is thermal-like. On the other hand, for τ a few a more complex spectrum forms, with the main contribution coming from synchrotron photons (emitted by the electrons) below the thermal peak and above it from multiple IC scatterings (leading to a nearly flat energy spectrum) [115]. All the above discussions is viable for dissipation processes from highly magnetized plasma as well [86,116].
However, also the above model suffers two major drawbacks, since it cannot explain 1) low energy spectral slopes less steep than the Rayleigh-Jeans part of a Planck spectrum; 2) the observed GeV emission, which may originate from some dissipation above the photosphere.

IV. RECONCILING COSMOLOGICAL INDICATORS TO GRBS
After the first part of this review, in which we faced the main properties of GRBs, their possible theoretical background and progenitors, we are now in condition to relate our understanding directly to cosmology. In fact, for cosmological purposes it is essential to get the distance of astronomical objects and so the use of GRBs would help in computing such distances up to very large redshifts. In particular, source physical parameters mostly depend on luminosity and size and then cosmic bounds can be inferred if there exists a relation between distances and redshifts. This prerogative is intimately related to two distinct concepts, i.e. distance indicators 17 and standard candles 18 .
Below we elucidate the main properties of such objects and the most important consequences they have in observational cosmology.

A. Distance indicators
At the beginning of our review, we emphasized how distances in cosmology are relevant to compute GRB luminosity/energy. A further step consists in noticing the distances measurements are classifiable by Absolute measures, as they are computed through previously known information, e.g. trigonometric parallax.
Relative measures, as they involve empirical relations based on indirect or direct probes, e.g. Cepheids period-luminosity relation, for which the distance measures are calibrated against an absolute method to enable those measurements to be somehow anchored.
Standard cosmology shows up how to relate the redshift to metric distances in both the above cases. The machinery of dynamical distance indicators involves tightly packed all the ingredients of cosmological physics. We thus require the cosmological principle to hold in an expanding universe in the context of general relativity. Despite obvious, there is no direct analogy to classical dynamical distance indicators, as the laboratory in which measurements are got is moving as well. Precision cosmology would enrich data during incoming years, as future surveys will provide resource of data to constrain and refine our understanding about distances and cosmological parameters.
Using current data catalogs, it appears evident that GRBs can be significantly investigated once the calibration of the correlation functions are deduced from absolute confidence. Recently, techniques of non-calibration have been more often used, overcoming the problem of standardizing GRBs that are, as known, not perfect standard candles for cosmological distance tests. Later on we confront the calibration and uncalibration procedures, emphasizing how to single out the most promising treatment to handle GRBs in cosmology.

B. Standard candles
Above we stated astrophysical distances are crucial for picturing current universe. Though essential, estimating cosmic distances mainly remains a complicate prerogative. In view of the above classification, the distance estimation passes through the use of standard candles. These objects hold the fundamental property of relating the intrinsic luminosity, namely L, to some known property, enabling one to get constraints over it. Once the luminosity is known, the distance can be computed accordingly.
A standard procedure is to get measures of the energy emitted from astrophysical objects. The energy bounds are got in a precise time interval, say ∆t and by virtue of E = L · ∆t, i.e., the relation between luminosity and energy, it is possible to get distances from the energy itself, through a well-consolidate strategy, below reported.
Detectors are able to catch fractions E d of the emitted energy E, which is proportional to the ratio between the detector area A and the spherical shell 4πd 2 L in which one defines the cosmic distance d L , i.e., A general relation for d L (z) is written as where the set of free parameters to constrain is indicated by θ. Exploring a given cosmological model is equivalent to get θ. Thereby, combining the aforementioned quantities, we obtain the energy per unit detector area A and per unitary time ∆t, which defines the flux expressed by As we highlighted, the luminosity L is known for standard candles, thus one can measure F in order to get a given astrophysical object distance.

C. Classifying standard candles
We above stressed that physical laws underlying a particular astronomical object permit one to know the luminosity of standard candles. Clearly, such rules are essentially based on thermodynamic or chemical processes of a given astrophysical object. Consequently, one can classify standard candles by means of these physical laws and, according to the simplest classification scheme, we can handle at least two kinds of standard candles summarized below [117].
• Standard candles as primary distance indicators, which can be calibrated within the Milky Way galaxy.
• Standard candles as secondary distance indicators, which can be observed at larger distances than Milky Way scales. However, they require calibration, typically performed using known primary distance indicators within distant galaxies.

Primary distance indicators
The above first typology mainly includes Variable stars, i.e., among which, Cepheids, RR Lyrae and Mira. Here, the variable star type is based on the possible correlation between their period of variation, steadily measured, and their luminosity. Even though this set of stars mainly constitutes the primary indicators, further typologies are main sequence and red clump stars. Here, using the luminosity-temperature relations from the standard Hertzsprung-Russell diagram, one deduces stellar luminosity within a fairly narrow range. Last but not least, eclipsing binaries are also primary distance indicators, since their luminosity is computed by the Stefan-Boltzmann law through a direct estimate of their radius, by means of a Doppler measurement of orbital velocities combined with the light-curve data, together with the temperature, deduced from the spectrum.

Secondary distance indicators
On the other hand, the second class of standard candles are essentially based on very different indicators with respect to the first case. For instance, the prototypes of such indicators are the properties of galaxies, among all, the Tully-Fisher relation. This law matches spiral galaxy rotation speed and stellar luminosity. In particular, to argue the spiral galaxy rotaion speed one can consider, for example, the spectral line width. Another relation, widely adopted as underlying second type of indicators, is the Faber-Jackson relation. Here, it is possible to infer elliptical galaxy random stellar velocities using the total luminosity. Again, the way to get these velocities consists in the use of spectral line widths. Another quite relevant relation is the fundamental plane law, i.e., a treatment that extends the Faber-Jackson one by including surface brightness as an additional observable parameter.
Besides galaxy properties, another second typology of standard candles is represented by SNe Ia, i.e., probably the most used cosmological standard candles to accredit the late time cosmic speed up. The scenario in which they form is due to thermonuclear explosions of WDs that exceed the Chandrasekhar's limit, namely ∼ 1.4M . For such objects, we see a correlation between the time scale of the explosion and the peak luminosity. The corresponding light curves follow given shapes, in agreement with the so-called Phillips curve [118]. As stated, SNe Ia are the most fruitful standard candles. For each event, even if the luminosity is clearly different for every SN, the Phillips curve relates the B magnitude peak to the luminous decay after 15 days with an overall set of SNe distributed in the range z = 0-2.5. These redshifts span between decelerating and accelerating phases of universe's evolution, corresponding to the matter and DE dominated epochs 19 . Last but not least, these indicators are present in all galaxies, except in the arms of spiral galaxies, but their physical internal processes are still object of investigations as they are not fully-interpreted.

V. GOING AHEAD WITH STANDARD INDICATORS: THE χ 2 ANALYSIS
Using standard candles, it is possible to establish data catalogs that can be used and matched with GRB data. Hence, to experimentally fit a given model with a given set of free parameters, one requires the definition of a merit function that quantify the overall agreement between the working model with the aforementioned cosmic data. Equivalently, it is of utmost importance to get best fit parameters and corresponding estimates of error bars, together with a method to possibly measure the goodness of fit. The parameter fitting treatment commonly makes use of least-squares analyses, based on the combination among data points, say D i , a model for these data, namely y(x, θ), function of θ. Naively the simplest approach to least squares for uncorrelated data becomes where the weights w i reach the maximum variance in case w i = 1/σ 2 i , with σ i the data point errors. For correlated data, we have in which the inverse of covariance matrix, Q, has been introduced describing the degree of correlations among data.
Minimizing the χ 2 is equivalent to get suitable sets of findings that represent the best fit for our procedure. Different χ 2 values lead to probability distribution around the minimum.

A. Probability distribution
Analyzing the probability distribution, once the above treatment is worked out becomes essential. In particular, probabilities p that the observed χ 2 exceeds by chance a value χ for the correct model is clearly calculable and, in fact, Q provides a measure of the goodness of fit, as one infers it at the minimum of χ 2 . Two limiting cases, unfortunately, are possible, Q is too small or too large. The first occurrence leads to the fact that the model is either wrong or errors are underestimated and/or they do not distribute Gaussianly. The second occurrence happens when either errors are overestimated or data are correlated while rarely it could also happen the distribution is non-Gaussian.
In general, the statistical procedure suggests χ 2 roughly comparable with the data number. Consequently, using the reduced chi square, as the ratio between the chi square and the number of degrees of freedom, could be an useful trick to handle experimental workarounds.

B. The SNe Ia measurements
SNe are widely-adopted in astrophysics as standard candles. Thereby, several SN catalogs are often updated, furnishing today a large number of data points that combined with other data sets enable one to fix tighter constraints over the universe expansion history in terms of its constituents. In particular, SNe Ia are likely the most used objects that constrain DE at late times. The standard procedure makes use of the luminosity distance d L (z) and of apparent magnitude. A general relation for d L (z) has been previously written, with θ the set of free parameters of a given model. Then, we can notice that exploring a given cosmological model is equivalent to get the whole set of parameters, θ.
In particular, when one adopts a given cosmological model then an indirect requirement naturally holds: the underlying cosmological model is the most suitable one. This is clearly a limitation because this hypothesis do not always coincide with the most feasible statistical model. So, more than one scenario can lead to subtle bounds, indicating a degeneracy problem among different models. This justifies the need of analyzing different cosmological paradigms working out data set hierarchy, i.e., combining more than one data catalog. In addition, statistical criteria are also crucial to check the goodness of a given paradigm.
For SNe Ia, by virtue of Eq. (20), it is possible to relate the brightness to fluxes to get the distance modulus Neglecting error bars on z, we underline errors on µ, namely σ µ , whereas the best fit is determined by the standard maximization of the underlying likelihood function, or simply minimizing the χ 2 , provided by where the subscript min refers to the set of values that minimize the chi square function, as above requested.
Theoretical models can be therefore tested by χ 2 statistics, leading to probe DE by inferring d L in units of megaparsecs and using it by means of the apparent magnitude. Again, intertwined more than one data set with other surveys is quite essential to determine the whole set of parameters, with refined accuracy. For instance, SNe alone, as well as GRBs 20 , H 0 cannot be arguable. In fact, expanding up to the first order the luminosity distance, valid up to z 0.001, one gets that clearly vanishes at z = 0, implying H 0 cannot be constrained with SNe Ia alone. Also, a multiplicative degeneracy between H 0 and the other free parameters occurs.
Once computed the chi square statistic, the confidence regions are planes with fixed χ 2 . For example, one can get Ω M − θ i planes by marginalizing the likelihood functions over H 0 . This procedure consists in integrating the probability density p ∝ exp(−χ 2 /2) for all values of H 0 . Marginalization is a generic technique, clearly not limited to H 0 . In fact, as one desires to simultaneously constrain a few parameters and meanwhile wants to get the corresponding probability distribution regardless of the values of a given parameter, say θ , can proceed with marginalizing. Let us call θ the parameter we do not care about, the marginalized probability density, computed for example for Ω m , is given by p(Ω m ) = dθ p(Ω m , θ ).

C. BAO measurements
The BAO measurements are due to overdensity of baryonic matter due to acoustic waves. These waves propagate in the early universe [119,120] and represent standard ruler for cosmological length scale. This signature, in the largescale clustering of galaxies, constrains cosmological parameters by detection of a peak in the correlation function [121], by defining the A parameter as follows where x is the set of cosmological density parameters, E(x, z) = H(x, z)/H 0 and sinn(x) = sinh(x) for the curvature parameter Ω k > 0, sinn(x) = x for Ω k = 0, and sinn(x) = sin(x) for Ω k < 0. The A parameter has been measured from the SDSS data and reads to be A = 0.469(0.95/0.98) −0.35 ± 0.017, with z 1 = 0.35, so the χ 2 in terms of A reads The BAO corresponding angular distance measures can be defined by means of The corresponding χ 2 is given by It is clear that BAO measures are slightly model-dependent as they depend on the comoving sound horizon r s (z d ). In particular, in Eq. (28) the sound horizon depends upon the baryon drag redshift z d . This quantity requires calibration that typically is performed with CMB data, adopting a given background model, that commonly is the ΛCDM scenario. Very often the best expected values are given by z d = 1059.62 ± 0.31 and r s (z d ) = 147.41 ± 0.30 [122].

D. Differential age and Hubble measurements
Another intriguing treatment, widely used in observational cosmology and also for calibrating GRB correlations, has been firstly proposed in Ref. [123]. The idea is to measure the Hubble rate by using galaxies, in a quite modelindependent way. In the context of GRBs, the Hubble catalog has been widely explored. For example, in Ref. [123], the core idea is to match the observational Hubble rate data (OHD) with model independent expansion of H made by Bézier polynomials. At a first glance, this differential age method (see, e.g., Refs. [124,125]) does not require any assumption over the form of H, although spatial curvature can affect the overall treatment if it varies with time, instead of being fixed 21 .
To better introduce the method, we notice it is well-known that spectroscopic measurements of the age difference ∆t and redshift difference ∆z of couples of passively evolving galaxies lead to ∆z/∆t ≡ dz/dt and so, if galaxies formed at the same time (redshift z), the Hubble rate can be approximated by Consequently, model-independent estimates may come from cosmic chronometers based on the assumption that observable Hubble rates are given by the exact formula if approximated as in Eq, (30). The χ 2 from the current 31 OHD measurements reads This procedure has the great advantage of directly considering H without passing through any cosmic distance.
E. The R parameter The CMB represents cosmic recombination epoch remnant and contains abundant early universe information. Consequently, the acoustic peak positions [120,126] can be used to characterize a given cosmological model by means of the shift parameter, defined as [127] The last scattering redshift, namely z ls , is fixed to where [128] and the χ 2 reads In analogy with BAO measures, the shift parameter is not fully model-independent.

F. Confidence levels and uncertainties
As one performs fits combining GRBs with other observable quantities, meaningful information on the best-fit parameters are achieved by computing their confidence limits or contour plots, which define the allowed parameter phase-space. These are essentially regions constructed around a set of best fit parameters got from computation. One does not mind about the number of dimensional parameter space, namely m, corresponding de facto to the number of parameters, since to make those regions compact one holds constant χ 2 boundaries, fixing the chi squared values to specific numbers. So one takes m to be the number of parameters, n the number of data and p be the confidence limit that one desires to reach. Assuming to shift by solving Q[n − m, min(χ 2 ) + ∆χ 2 ] = p, and to find the parameter region where χ 2 ≤ min(χ 2 ) + ∆χ 2 , immediately one gets the requested confidence region. Once the regions have been computed, it is needful to get uncertainties. To do so, expanding the log likelihood in Taylor series ln L = ln L(θ 0 ) (θ j − θ j0 ) + ..., we define the Hessian matrix by Since its non diagonal terms indicate correlated parameters, one can assume the errors on a given i parameter to be 1/ √ H ii . This naive representation of errors is a coarse-grained approach, dubbed conditional error, not frequently adopted in the literature. On the other hand, one can compute the Fisher information matrix, as a forecast expression for error bars with the ensamble average over observational data. In analogy to conditional errors, we write σ 2 ij ≥ (F −1 ) ij , while the marginalized errors become σ θi ≥ (F −1 ) ii . Above we underlined that the Fisher matrix is somehow related to error bars. In this respect, we mean the Fisher Information matrix enables to estimate the parameters errors before the experiment is performed. Hence, it permits to explore different experimental set ups that could optimize the experiment itself. For these reasons, the Fisher matrix is largely adopted in the literature.

G. Binning procedure
In several cases, it is useful to get constraints directly on the universe equation of state. Thus, fitting it for the latetime universe constituents is extremely important to understand the dark energy evolution. In particular, pointing out a possible variation of the equation of state of dark energy is essential to disentangle the standard model predictions from possible theoretical extensions and, in this respect, GRBs can be seen as intermediate redshift probes to disclose such an evolution.
To do that, an intriguing strategy consists in binning the dark energy equation of state, say w, in short intervals of z and then fit w in each bin, assuming it is constant in each bin. Indicating with a generic function f (z) the dark energy evolution, we have where w i is the barotropic factor within the i th redshift bin. The bin is built up by an upper boundary at z i , whereas the zeroth bin is defined as z 0 = 0. There uncorrelated sub-equations of state in every bin can be experimentally refined adding data points and, in particular, GRBs, being calibrated as we will discuss later. Several indications have shown good agreement with the standard paradigm, up to z 9, albeit relevant deviations have been found, indicating the situation is not still clear.

VI. STANDARDIZING GRBS
Being successful in standardizing GRB data is of utmost importance to characterize new data catalogs up to high redshifts. In particular, getting redshifts, or more generally spectroscopic observations, is essential for GRB-related science, as we summarized below.
(1) Computing the luminosity function for GRBs, constructing it from the prompt emission as well as afterglows.
This treatment is analogous to what we do for SNe Ia.
(2) Computing the redshift distribution of GRBs. This enables one to use GRBs as tracers for the cosmic starformation history. Consequently, spotting very high redshift GRBs will shed light on their distribution at intermediate epochs of the universe evolution.
(3) Studying the host galaxies, in particular those faint, high-redshift galaxies that are unlikely to be found and studied with other methods. Characterising the dust extinction curves of high-z galaxies.
(4) Studying GRB-selected absorption line systems and probing cosmic chemical evolution with GRBs.
(5) Studying if and how much GRBs can be used for determining the cosmological parameters of dark energy models and/or to rule out a few models. Analogously the use of GRBs can be tested in view of determining cosmographic parameters, i.e. getting model independent bounds over the cosmic evolution.

A. GRB correlations and related issues
Since the first discovery of GRBs independent groups have found different correlations, that represent a key to use GRBs for cosmological purposes. The basic idea is to intertwine different quantities of such objects among them. The observable quantities of interest are in relation with the cosmological model that lies on the background. This fact permits GRBs to be distance indicators at a first glance, but limits their use because requires to postulate the underlying cosmological model, providing a circularity in the process itself which is know as the circularity problem.
The widest majority of GRB correlations prompt the same requirement: the GRB standardization in terms of cosmological tools. Attempts for new correlations have been severely investigated, relating different observable quantities with each other. The way in which this is realized provides the theoretical interpretation behind the relation itself. In other words, evidence for a given correlation leads to interpret particular physical processes. Thus, achieving the goal of standardizing GRBs brings the certainty of getting feasible bounds on cosmological parameters. Intriguingly, a narrow set of correlations enables one to estimate GRB redshifts also. Even though this is still under speculation, in general a wide number of correlations could provide information about GRB progenitors.
More precisely, standardizing GRBs for cosmological purposes aims at reaching further hints toward progenitors of different groups of GRBs. Multi-wavelength instruments of recently-adopted satellites have significantly increased the number of GRBs that could be observed to check the validity of a given relation. So, it is even possible a few correlations may be derived from experimental evidence, instead of theoretically. Unfortunately, this could open further issues related to data processing whose outputs can be biased in the overall computations.
Going ahead, it is certainly possible to constantly observe new hints undertaking novel correlations to let free theoretical speculations that deeply search into new physics beyond the standard comprehension of GRBs.

B. Prompt emission GRB correlations
The correlations that make use of prompt emission quantities are listed below with the corresponding properties of each of them. For details on the involved quantities, see Sec. II E.

This correlation holds for
LGRBs and indicates that the more luminous bursts possess shorter time lags, i.e., L iso ∝ τ −1.25 lag [129]. It has been used as a GRB redshift indicator and to constrain cosmological parameters. However, the existence of this correlation is challenged by recent studies. For details, see Ref. [130].

This correlation holds for
LGRBs and indicates that the more luminous bursts have the more variable light curves [131]. However, the intrinsic scatter is very large and the index is still not completely settled, also in view of the fact that the variability V is different for various instruments. For details, see Ref. [130].

Amati or Ep-Eiso correlation
This correlation is of the form E p ∝ E 0.52 iso and shows that E iso is correlated with the rest-frame spectral peak energy, namely E p = E obs p (1 + z) [132]. Observations by Swift and Fermi detectors confirmed this correlation for LGRBs. An analogous E p -E iso correlation with a slope similar to that of LGRBs but a larger value of the normalization holds also for SGRBs, though with much smaller data set [133]. Moreover, the E p − E iso correlation also holds within individual GRBs using time-resolved spectra, and the slopes are consistent with the correlation from time-integrated spectra. For details, see Ref. [130].

Yonetoku or Lp-Ep correlation
This correlation reads L p ∝ E 2 p [65] and holds for both LGRBs and SGRBs [133]. Similarly to the E p -E iso correlation, the L p -E p correlation holds also within individual GRBs using time-resolved spectra. For details, see Ref. [130].

Ghirlanda or Ep-Eγ correlation
This represents a tight, less scattered, correlation between E p and E γ , valid for LGRBs [134]. One of the major drawbacks is the lack of achromatic breaks in the Swift afterglow light curves of most of GRBs. This fact limits the increase in the correlation sample. For details, see Ref. [130].

C. Prompt and afterglow emission correlations
The following correlations involve prompt and afterglow emission observables. Below we report the most common correlations. For details on the involved quantities, see Sec. II E.

Liang-Zhang or Ep-Eiso-t b correlation
It is a correlation valid for LGRBs among E iso , E p and the rest-frame break time in the optical band t b , i.e., E p ∝ E 0.52 iso t 0.64 b [135]. If we take the optical break time as the jet break time, this correlation is similar to the Ghirlanda one. However, the inclusion of additional GRBs made this correlation more scattered. For details, see Ref. [130].

Dainotti or LX-tX correlation
This correlation links the X-ray luminosity L X and rest-frame time t X , i.e., L X ∝ t −1 X , at the time when the X-ray afterglow light curve establishes a power-law decay after the plateau phase [136]. This correlation holds for LGRBs and SGRBEEs. By adding a third parameter, E iso , it has been found the new correlation of the form L X ∝ t −0.87 X E 0.88 iso [137]. However, both relations are quite scattered and seem to be a selection effect due to the flux detection limit of Swift-XRT instrument. For details, see Ref. [130].

E X iso -Eiso-Ep correlation
This is a universal correlation for both LGRBs and SGRBs which links E iso and E p to the isotropic energy of the X-ray afterglow E X iso computed in the rest-frame energy band 0.3-30 keV, i.e., E X iso ∝ E 1.00 iso E 0.60 p [138]. However, due to the fact that this correlation depends upon two cosmology-dependent quantities, it is unsuitable to constrain cosmological parameters. For details, see Ref. [139].

Combo correlation
The correlation represents the combination of the Amati correlation and the E X iso -E iso -E p correlation and, like the Amati correlation, it holds for LGRBs only [140]. It relates the prompt emission E p and the X-ray afterglow luminosity L 0 and rest-frame duration τ of the plateau phase, and the late power-law decay index α, i.e., L 0 ∝ E 0.90 The main drawback is that it is much more complicated than other correlations, as it depends upon four parameters. For details, see Ref. [141].

L-T -E correlation
This correlation connects the rest-frame end time t X and luminosity L X of the X-ray afterglow plateau phase with E iso , i.e., L X ∝ t −1.01 X E 0.84 iso [142]. This correlation is very similar to the Combo correlation but unlike it, the L-T -E one holds for a few SGRBs and requires to be corrected for the redshift evolution effects. For more details, see Ref. [142].

VII. FURTHER ISSUES RELATED TO CONSTRUCTING GRB CORRELATIONS
Despite having so many correlations reported in the literature, all of them suffer of several issues related to constructing the correlations themselves. Besides the circularity problem, we face several issues to be addressed. In the following we list the most common ones.

A. Evolution effects
GRBs are observed from a large redshift range and in principle correlation parameters may evolve with the redshift. In some cases, to estimate the cosmological parameters a correction of the kind of (1+z) −d to the energy or luminosity parameters is needed, introducing a further parameter to fit. However, though tested for subsamples of GRBs and several correlations, this issue is still ongoing [143].

B. Instrumental selection effects
A yet open issue are the instrumental selection effects that may affect the observed GRB energy or luminosity correlations. There are at least two kinds of issues due to a) the trigger threshold, i.e., the minimum photon peak flux that a burst must have in order to be detected by a given instrument, and b) the spectral analysis threshold, i.e., the minimum fluence to perform a reliable spectral analysis and determine the SED parameters.
Several analyses have been performed in the literature, using different samples and correlations of GRBs and searching for outliers and inconsistent GRBs. Several conclusions were drawn: a) some correlations may exist, though due to selection effects; b) other correlations may exist when accounting for the intrinsic scatter; c) some correlations may have statistical significance, though affected by the thresholds of GRB detectors, etc. See Ref. [144], for more details. Interestingly, using the time-resolved spectra, similar correlations were found in individual bursts, strongly supporting the fact that the correlations may be physical [64].

C. Systematic errors
Sources of systematic errors for GRBs are the sensitivity of the detectors, the differences in the estimated spectral parameters depending on detectors and/or fitting models, the lack of unknown parameters, etc. All these might dominate over the intrinsic dispersions of GRBs.
In general, due to the vast number of systematic errors, a technique to consider them in GRB fitting procedures consists in deriving those errors requiring the chi-square to be comparable to the fitted degrees of freedom ν, namely χ 2 ν, and then summing them in quadrature with the statistical errors [145].
However, studies on different GRB correlations suggest that the systematic uncertainties of correlation parameters are not sensitive to the assumptions about cosmological parameters [146].

D. Issues and interpretation of prompt emission GRB correlations
In the following, attention is given to the above described GRB prompt emission correlations, focusing on possible issues and their physical interpretation. Particular emphasis is given to the functional form of E p -E iso (or Amati), Ghirlanda and Yonetoku correlations, the most popular and most quoted correlations for prompt emission observables.
As established in Ref. [147], this correlation holds also for X-ray flares (in the rest frame energy band 0.3-10 keV) and proposed that their underlying mechanism is similar. However, this correlation is affected by evolution effects, as discussed in Ref. [130].
One of the latest proposed explanations of the L iso -τ lag correlation involves only kinematic effects [148], as the observed time-lag is τ lag ∝ D −1 and the luminosity is L iso ∝ D, where D is the Doppler boost defined in Sec. III A 5.
As discussed in Ref. [149], this correlation was constructed from a small sample of heterogeneously collected GRBs and is severely affected by sample incompleteness.
This correlation has a non-negligible scatter, thus, it is the least reliable one among all GRB prompt emission correlations [130]. The physical origin of the L iso -V correlation is still unclear. Within the internal shock scenario (see Sec. III A 2), it seems to be related to the activity of the central engine through the values of Γ and the jetopening angles [150] and, based on this interpretation, L iso (V ) is proportional to a high (low) power of Γ, hence high-luminosity pulses imply high variability prompt light curves [145].
As discussed above for the L iso -τ lag correlation, the luminosity-variability correlation as well is severely plagued by sample incompleteness [149].

Amati or Ep-Eiso correlation
Likely, the most used and investigated relation is represented by the so-called E p -E iso or Amati correlation [151] that can be here recast by Here, we have two free constants, namely a 0 and a 1 , that represent the calibration constants to determine once the relation is somehow calibrated. A possible limitation of the E p -E iso correlation is due to the extra source of variability σ a . This is thought as a direct consequence of hidden variables that contribute to the overall calibration, albeit we cannot directly observe them [152]. A possible explanation for the E p -E iso correlation considers the thermal radiation emitted when the GRB jet drills through the core of the progenitor star (see Sec. II D 1), responsible for the thermal peak in the spectrum, and the Compton scattering of this radiation by relativistic electrons outside the photosphere (see Sec. III A 5 and Ref. [153], for details).
There are claims that the Amati correlation is caused by some selection effect of observations, rather than being an intrinsic property of GRBs [149,154]. However, there is a general consensus on the fact that the correlation is real [155][156][157], though detector sensitivity affects the correlations and a weak fluence dependence may be larger than the statistical uncertainty and contribute to the dispersion of the correlation [158,159].

Ghirlanda or Ep-Eγ correlation
From theoretical and observational arguments in favor of the jetted nature of GRBs [36], the radiated GRB energy can be corrected by means of the collimation factor f = 1 − cos θ, leading to E γ = f E iso . In particular, the jet opening angle θ is evaluated at the characteristic time t b for specific assumptions on the circumburst medium, that can be assumed homogeneous [28]. The functional form adopted here for the Ghirlanda relation reads in which, as usual, b 0 and b 1 are the two free constants, fixed by means of calibration. It behooves us the extra scatter, σ b , to better constrain the relation itself. The Ghirlanda correlation shares with the E p -E iso one a similar physical interpretation (see Sec. III A 5 and Ref. [153], for details). This correlation takes into account also the jet correction in the computation of the GRB energy output (see Sec. II D).
The Ghirlanda correlation is linked to the Amati one and, thus, criticisms/analyses against/in favor of being an intrinsic property of GRBs. In Ref. [160] it was shown that as many as 33% of the BATSE bursts would not be consistent with the Ghirlanda relation, but this results depended upon the assumed distribution for the jet's correction factor f [161]. This fact limits the increase in the correlation sample. For details, see Ref. [130]. Likewise for the Amati correlation, the Ghirlanda one is statistically real but strongly affected by the thresholds of GRB detectors [159].

Yonetoku or Lp-Ep correlation
The Yonetoku or L p -E p, [65] correlation functional form here adopted reads Here, the free terms are m 0 and m 1 and require to be calibrated. Again, σ m is the extra scatter term. Detailed hydrodynamical simulations suggest that Yonetoku correlation may be due to the emission of photons from the photosphere of a relativistic jet, where the outflow becomes optically thin, whereas most of it is still optically thick (see Sec. III A 5). Quasi-thermal radiation is thus expected and the expected spectral shapes are obtained [162].
Like previous correlations, also for the Yonetoku correlation there are ongoing discussions whether it is a by-product of some selection effect or not [149,[154][155][156][157]. For this correlation, however, a weak redshift dependence has been confirmed, which may contribute to the dispersion of the correlation [158].

E. Issues and interpretation of prompt and afterglow emission GRB correlations
In analogy with the prompt emission correlations, we here focus on possible issues and physical interpretation of prompt and afterglow emission correlations. The emphasis on the functional form is given to the Combo correlation, one of the less scattered correlations for both prompt and X-ray observables without evolution effects.

Liang-Zhang or Ep-Eiso-t b correlation
The E p -E iso -t b correlation [135] has been proposed by purely considering phenomenological considerations, thus avoiding any theoretical assumption, unlike made for the Ghirlanda correlation. However, as above stated, the E p -E iso -t b somehow shares similar implications and drawbacks than the Ghirlanda one, as well as analogous physical interpretation (see Secs. II D and III A 5 and Ref. [153], for details).
As proposed in Ref. [163], this correlation, like the Ghirlanda one, appears to be affected by selection effect on E p (whereas for the Ghirlanda correlation E γ is affected as well) and suffers sample incompleteness.

Dainotti or LX-tX correlation
The Dainotti correlation, akin to the L iso -τ lag correlation, may be retrieved from kinematic effects (see Sec. III A 5 and Ref. [148], fro details), pointing out to a common origin between the two correlations.
As already discussed above, this correlation is quite scattered and seem to be a selection effect due to the flux detection limit of Swift-XRT instrument [130]. Moreover, it is also affected by evolution effects with the redshift [164].

E X iso -Eiso-Ep correlation
This correlation depends upon two cosmology-dependent quantities, although it is unsuitable to constrain cosmological parameters. Since it holds for both SGRBs and LGRBs, with different progenitor and surrounding medium properties (see Secs. II D 1 and II D 2), its physical interpretation has not been yet established. There is a speculation that it may be connected with the Γ of the outflow, which might regulate the efficiency of conversion from γ-rays to X-rays [138].
The E X iso -E iso -E p correlation utilizes the prompt emission observables E iso and E p on which the Amati correlation is based. For this reason it is straightforward to deduce that the biases and selection effects, at work for the Amati correlation, partially affect this hybrid correlation. Moreover, unlike pure afterglow correlations such as the Dainotti one, this correlation is also plagued by double truncation in the flux limit, both in the prompt and X-ray afterglow emissions, making difficult the correction for any selection effect and the use as redshift estimators and cosmological tool (see discussions in Ref. [164]).

Combo correlation
The Combo correlation is a hybrid correlation linking the prompt emission E p and the observable quantities determined from the X-ray afterglow light curve, i.e., among all the rest-frame 0.3-10 keV plateau luminosity L 0 , its rest-frame duration τ , and the late power-law decay index α [140]. For each GRB, L 0 , τ , and α can be obtained by fitting the rest-frame 0.3-10 keV flare-filtered afterglow luminosity light curves with the function L(t) = (1 + t/τ ) α . 22 The general expression is much more complicated than previous ones and reads Here, the constants k 0 and k 1 need to be determined by means of the calibration procedure. Again, the correlation is characterized by an extra scatter σ k . The Combo correlation can be explained by the external shock scenario (see Sec. III). The correlation is the result of the synchrotron emission from the electrons accelerated in a relativistic shock (see Sec. III A 5). The shock propagates through the external CBM and interacts with the magnetic field of the the turbulent plasma. Hence, the relationship among E p , L 0 and τ and the corresponding comoving quantities scale with the initial Lorentz factor of the bulk motion Γ 0 , whereas the intrinsic scatter is due to the uncertainties on the source spectral energy distribution [165].
Keeping in mind the hybrid nature and the Combo correlation, which a combination of the Amati and the E X iso -E iso -E p correlations, the same biases and selection effects at work in the E X iso -E iso -E p affect the Combo correlation as well.

5.
L-T -E correlation As already stated, this correlation shares similarities with the Combo correlation but, unlike it, needs to be corrected for the redshift evolution effects [142].
The L-T -E correlation, as well as the Combo correlation, may be explained within the magnetar scenario, which justifies the plateau phase observed in the X-ray afterglow light curves as due to the continuous energy supply from a supra-massive NS (see Sec. III A 3 and Ref. [71], for details).
Likewise the Combo correlation, which share similar features, the discussion on the same discussion on the biases and selection effects holds for the L-T -E correlation as well. In addition, in Ref. [142] it is evidenced that GRBs of the sample at high redshifts usually have relatively larger L X and E iso . This is very likely a selection effect due to the fact that very distant GRBs with small L X and E iso cannot be observed in view of the limited threshold of current detectors.

VIII. CIRCULARITY OR NOT CIRCULARITY?
Even though GRBs could be thought as indicators toward the determination of the universe's expansion history, it is remarkable to stress that small redshift GRB data are still missing. Consequently if one requires to calibrate GRBs with small redshift data, there is the strict need of other data sets, whose data points lie around z 0, that permit to perform the calibration procedure. The calibration procedure is essentially a consequence of the correlation functional forms that, by virtue of the above considerations, are commonly written as y = a x + b, with a, b real constants that depend upon the relation itself. Bearing this in mind, we focus below on different calibration strategies. In this respect, we highlight whether it is absolutely necessary or not to calibrate our correlations in order to confront cosmological models with GRB data. Below we start with such considerations, critically discussed.

A. Calibration versus non-calibration
Constraining cosmological parameters using GRBs is plagued by several conceptual and practical issues. First, all GRBs correlations by definition are built assuming an a priori background cosmology (see Secs. II E and VI A for details) and, consequently, introducing a circularity problem [166]. Second, almost the majority of GRB correlations holds for LGRBs (see Sec. VI A for details), whose observational rate falls off rapidly at low-z and in some cases such a nearby LGRBs seems to be intrinsically different from the cosmological ones [37,167]. Third, unlike SNe Ia which are calibrated with a selected sub-sample at very low redshift 23 by anchoring them to primary distance indicators as Cepheids, the shortage of low-z GRBs prevents to anchor them to primary distance indicators.
Focusing on the above circularity problem, this is essentially an epistemological issue due to the lack of very lowz GRBs and arising from the need of a background cosmology to compute the above-defined E iso , E γ and L iso entering GRB correlations [166]. For example, calibrating GRBs through the standard ΛCDM model, the estimate of cosmological parameters of any dark energy framework inevitably returns an overall agreement with the concordance model. Debates toward its use in cosmology seem to indicate that this effect could be minor, albeit it plagues cosmological constraints got from GRBs.
Possible way outs known in the literature, involve calibration techniques based on the use of SNe Ia distance moduli, cosmographic series, cosmic chronometers, etc. All these procedures represent plausible solutions, but at the same time introduce possible issues that we are going to describe below.
The above calibration procedures have to compare with the an alternative method which completely by-passes the calibration procedure [10]. This uncalibrated procedure consists in a simulataneous fit of correlation parameters together with the cosmological model parameters. This uncalibrated procedure consists in a simultaneous fit of correlation and cosmological model parameters. This procedure and the related issues are also discussed in the following.
B. Fitting procedures with calibration

SN calibration
A widely-used method to calibrate GRB correlations is through the use of SNe Ia, that span within z 2.3. In such a way, assuming this could work for any LGRBs, the GRB data points are mixed with SNe in order to build up a whole, quite large, Hubble diagram, where in the small redshift domain one has the majority of SNe, while at large z, GRBs are the most. Here, the simplest error bars on distance modulus are [168,169] where µ,i and µ,i+1 and µ i and µ i+1 are the errors and distance moduli of the SNe Ia at z i and z i+1 , respectively. For each SN catalog, we could find different interpolating functions to model the SN distribution. Thus, calibrating GRBs with SNe would seriously depend on the choice of these expressions for each catalog. Hence, GRB calibrations may turn out to be extremely sensitive to SNe Ia and the approach should be carefully handled since GRB luminosity correlations may be no longer fully independent from SN data points.

Model dependent calibration
The model dependent procedure fixes the background cosmology with a given cosmological model, where typically the dark energy evolution is assumed a priori. Since the background cosmology 24 enters the correlation functions, generically in the form y = ax + b, this means that one has to 1) assume a background cosmology, 2) fix the most suitable numerical bounds over the free coefficients of the background cosmology, and 3) calibrate the correlation.
As it appears evident, this strategy consists in determining an accredited cosmological model with particular choices of the free parameters, determined elsewhere.
This procedure is obviously strongly plagued by the circularity problem. It fixes the cosmological evolution with a given model and does not permit to constrain suitably another cosmological paradigm. In fact, if one calibrates with a generic model, say H (1) (z), any other statistical expectations on a different model, say H (2) (z), would favor the model that better matches H (1) (z). In other words, calibrating with H (1) (z) implies that the best fits are statistically argued for H (2) Another dramatic fact is that one has to constrain the free parameters of the background scenario by means of additional fits, with different data sets. This implies that an overall analysis would be plagued by error propagation between different catalogs of data and limits severely the analysis itself. To avoid other fits, one can assume exact versions of the cosmological models that should be used as backgrounds. In such a way, the corresponding error propagation reduces, albeit one does not take a real tested cosmological scenario, but rather a simplified version of it.

Model independent calibration
Calibrating correlations via model independent treatments permits to use GRBs as distance indicators, although the calibration is made by means of other standard candles. The idea of model-independent calibrations, however, enables to get the luminosity distance d L without a priori postulating the background cosmology, healing de facto the circularity problem.
A nice possibility consists of relating distances with model-independent quantities written in terms of a Taylor series expansion of the scale factor. Thus, we first notice and then we consider the following expansions where α n are the coefficients of the luminosity distance expansion and H m are the coefficients of the Hubble rate expansion. Thus, baptizing the cosmographic set, q 0 , j 0 , s 0 , as the present values of the following quantities where the scale factor a has been considered, with the requirement a ≡ (1 + z) −1 . The above quantities are named deceleration, jerk, and snap parameters, respectively, we formally have The coefficients in Eqs. (45)- (46), say α i ≡ α i (q 0 , j 0 , s 0 ) and H i ≡ H i (q 0 , j 0 , s 0 ), can be determined directly with data, without considering a cosmological model a priori. This treatment is known with the name of cosmography or cosmokinetics, i.e., the part of cosmology that reconstructs the universe's kinematics model-independently. So, at z = 0 we have and Although powerful, the above formalism suffers from shortcomings due to the convergence at higher redshifts 25 , i.e., the high GRB redshifts are very far from z = 0. 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.
Healing the convergence problem leads to a high-redshift cosmography. In this respect, several strategies have been suggested. For instance, one could 1) extend the limited convergence radii of Taylor series by changing variables of expansion, using the so-called auxiliary variables, or 2) changing the mathematical technique in which the expansions are performed, i.e. involving expansions different from Taylor ones, etc.
In the case of auxiliary variables, one employs a tricky method in which the expansion variable is reformulated as a function of the redshift itself, but with particular convergence properties. In other words, we cosmic quantities are rewritten in a more complicated function of the redshift z, namely y. Changing the redshift variable from z to y modifies accordingly the convergence radius. Formally speaking we write y ≡ F(z) [170], where we assume F(z) a generic function of the redshift. The function F(z) is properly chosen from physical prime principles. All F(z) prototypes, however, might fulfill a few mathematical conditions: 25 For the sake of completeness, this problem is not related to GRB redshifts only, but it remains an open issue of cosmography.
The first guarantees that at z = 0, our time, even y is zero. The second that the auxiliary variable does not diverge, otherwise the convergence problem would still persist. In addition, a further requirement is helpful in constructing F(z): The former condition enables y to converge at future time, as well as z. Using these hints toward the formulation of F(z), we can suggest a couple of well-consolidate examples of F(z): The second possibility is offered by rational approximations, where the expansion is thought to be a rational function, instead of a polynomial. This guarantees to optimize the Taylor series with rational approximants that better approach large z than the Taylor series, guaranteeing mathematical stability of the new series if data points exceed z = 0. Among all the possible choices, here the attention is given to the Padé polynomials, firstly introduced in Ref. [171]. This technique of approximations turns out to be a bookkeeping device to keep the calculations manageable for the cosmography convergence issue. Thus, provided we have Taylor expansions of f (z) under the and requires that b 0 = 1. Furthermore, it is important that f (z) − P n,m (z) = O(z n+m+1 ) and the coefficients b i come from solving the homogeneous system of linear equations m j=1 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 = i k=0 b i−k c k . Just for an example, we report the (2, 1) Padé polynomial as

C. The use of Bézier polynomials
The left term of Eq. (30) can be approximated by means of particular choices, such as using model-independent Bézier parametric curves. They are constructed to be 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 formally are defined as where we notice the linear combination of Bernstein basis polynomials h d n (z). Assuming the coefficients β d of to be positive in the range of 0 ≤ z/z m ≤ 1, where z m is the maximum z of OHD, we soon can classify those polynomials by means of the exponent n.
In particular, besides the constant case, n = 0, both linear growth, that happens for n = 1 and oscillatory regimes, say n > 2, work well. This implies that a suitable choice is n = 2. In this case we have The comparison between H 2 (z) and the OHD data points give β 0 = 67.76 ± 3.68, β 1 = 103.3 ± 11.1, and β 2 = 208.4 ± 14.3, all in units of km s −1 Mpc −1 .
After having approximated H(z) with Eq. (60), for spatially flat cosmology, Ω k = 0, the calibrating luminosity distance becomes Once the luminosity distance is written, it is possible to calibrate E iso , E γ , L p , and L 0 for the correlations that we intend to test. For E p − E iso , Ghirlanda, Yonetoku and Combo correlations we report in Table I the corresponding numerical outcomes related to the calibration process. Once calibrated, the corresponding distance moduli from Eq. (24) are computed for each correlations.

Simultaneous fits
Another relevant strategy is based on the idea to constrain the cosmological parameters together with the luminosity correlation [172,173]. In particular, the real distance modulus can be computed as where the sum is over a given number of different correlations. In particular, µ i is the best estimated distance modulus and the subscript i-th refers to the correlation, with σ µi the error bars. The uncertainty of the distance modulus for each burst is σ µ fit = ( i σ −2 µi ) −1/2 . A great advantage is that, as one computes bounds on cosmological parameters, the normalization functions and slopes of each correlations are marginalized. Consequently, we write down the χ 2 as where µ fit,i and σ µ fit,i are the fitted distance modulus and its error respectively.

Narrow calibration
Another intriguing technique consists in calibrating standard candles using GRBs in a narrow redshift range, hereafter δz. This short interval is placed near a fiducial redshift [174,175] with the great advantage that, in some cases, see e.g. Ref. [174], no low-redshift GRB sample is needful.

D. Fitting procedures without calibration
Constraints on the cosmological parameters can be obtained with an alternative method which completely by-pass the calibration procedure. It consists in taking all best data set of any GRB correlations, introducing the lowest intrinsic dispersion. The method is described below.
Any correlation between a generic energy/luminosity quantity Y and a GRB observable X has the form log Y obs = a log X + b .
The energy/luminosity quantity in general contains the information on the cosmological parameters Ω i through d L (z, Ω i ), defined by the theoretical model describing the background cosmology where F bolo may be the rest-frame bolometric fluence S bolo (1 + z) −1 for Amati-like correlations, or the bolometric observed flux F bolo for Yonetoku-like correlations. Please, notice that we only focus on these two relations just for giving an example. The same can be reformulated for other correlations, although for brevity we do not report other treatments here. The best cosmological and correlation parameters are then obtained by maximizing the log-likelihood function [152] ln where is the error in the measured value of log Y obs i , σ log Xi is the error in log X i , and σ ext is the intrinsic dispersion of the correlation.
The above treatment avoids to calibrate GRB data. This procedure has been applied to an uncalibrated E p -E iso correlation in Ref. [176] and to a uncalibrated Combo correlation in Ref. [10]. In both the cases, the correlations have been built up from samples of GRBs that have lower intrinsic dispersion. As byproducts, the resulting GRB correlations are close for different cosmological models, which can be interpreted with the fact that this procedure is model-independent. However, the application of this method seems to indicate that current GRB data are not able to put stringent constraints on cosmological parameters, though consistent with those resulting from better-established cosmological probes.
However, as hinted by the above results, this method may introduce a possible bias, namely that GRB correlation may adjust itself to the cosmology that maximize Eq. (66), rather than allowing the derivation of Ω i from a cosmology independent calibrated correlation. Moreover, it may be possible that more exotic cosmological model would lead to best-fit GRB correlations significantly different from simpler model, thus failing in providing a model-independent procedure. Therefore, this method is yet the subject of ongoing studies.

A. Numerical results using correlations
Observational data indicate that the cosmological expansion is currently accelerating. They also indicate that in the recent past the expansion was decelerated. The standard spatially-flat ΛCDM model [177][178][179] is the simplest model consistent with these observations [180][181][182][183]. Here, a cosmological constant Λ dominates the current energy budget and fuels the currently-accelerating cosmological expansion. In this model, above a redshift z ≈ 0.75, non-relativistic cold dark matter and baryons dominate over Λ and powered the then-decelerating cosmological expansion. While the observations are consistent with dark energy being time-and space-independent, they do not rule out slowly-evolving and weakly spatially-inhomogeneous dynamical dark energy, nor spatial flatness.
Significant constraints on cosmological parameters come from the CMB anisotropy data -that primarily probe the z ∼ 1100 part of redshift space -as well as from BAO observations -the highest of which reach to z ∼ 2.3 -and other lower-redshift SNeIa and OHD measurements. Observational data in the intermediate redshift range, between z ∼ 2.3 and ∼ 1100, are not as constraining as the lower and higher redshift data, but hold significant promise. Intermediate redshift observations include those of HII starburst galaxies that reach to z ∼ 2.4 [184][185][186][187][188][189][190], quasar angular sizes that reach to z ∼ 2.7 [188,[191][192][193][194][195], quasar X-ray and UV fluxes that reach to z ∼ 7.5 [195][196][197][198][199][200][201][202][203], as well as GRBs, that have now been detected to z = 9.4 [204]. Observed correlations between GRB photometric and spectroscopic properties that can be related to an intrinsic burst physical property would allow GRBs to be used as valuable standard candles that reach high z and probe a largely unexplored region of cosmological redshift space [see e.g. 136,140,142,145,151,[205][206][207][208], and references therein], similar to how SNeIa are used as standard candles [118] at z < 2.3. However, as stressed many times in this review, this is still a challenge for GRBs.
After it was established that GRBs were at cosmological distances, many attempts have been made to use burst correlations to constrain cosmological parameters. The first GRB Hubble diagram of a small sample of 9 bursts, obtained in Ref. [209] from the L iso -V correlation [131], led to a current non-relativistic matter energy density parameter limit of Ω m0 < 0.35 at the 1σ confidence level (for the flat ΛCDM model). Soon after, using the Ghirlanda correlation, in Ref. [210], with a sample of 12 bursts, it has been found Ω m0 = 0.35 ± 0.15 for the flat ΛCDM model, and in Ref. [211], with 14 GRBs as well as SNeIa, it has been inferred Ω m0 = 0.37 ± 0.10 and a cosmological constant energy density parameter Ω Λ = 0.87 ± 0.23, in the non-flat ΛCDM model, and Ω m0 = 0.29 ± 0.04 in the flat model. Similar constraints were obtained in Ref. [135], using the E p -E iso -t b correlation: 0.13 < Ω m0 < 0.49 and 0.50 < Ω Λ < 0.85 at 1σ confidence level in the flat ΛCDM model.
More recently, many contrasting results have been reported in the literature. In Ref. [212] cosmological parameter constraints, within the ΛCDM model and dynamical dark energy models from two different GRB data sets, were found to be different from the two data sets and also relatively broad. Similarly, in Ref. [213] it has been shown that at that time GRB data could not significantly constrain cosmological parameters. In addition, in Refs. [214,215] it has been shown that most GRB correlations have large scatter and/or their parameters differ somewhat significantly between low-and high-z GRB data sets. From the calibration of the Ghirlanda correlation, by using a SNeIa distance-redshift relation -through the (3, 2) Padé approximant-, in Ref. [214] it has been obtained Ω m0 = 0.302 ± 0.142 within the flat ΛCDM model. 26 Based on a cosmographic approach, an updated E p -E iso correlation with 162 GRBs has been used to get cosmological constraints. In Ref. [12] GRBs were calibrated with SNeIa, resulting in Ω m0 = 0.25 +0. 29 −0.12 within the flat ΛCDM model, whereas in Ref. [217] a cosmographic expansion, up to the fifth order, involving SNeIa is used to calibrate the E p -E iso correlation for GRBs, which are then used in conjunction with OHD and BAO measurements to constrain cosmographic parameters, resulting in a 1σ deviation from the ΛCDM cosmological model.
Other recent works (involving GRB data only or in conjuction with other probes) also report inconsistencies with the ΛCDM model. In Ref. [218] the E p -E iso correlation has been used, including also modeling the potential evolution of GRB observables, to conclude that calibrated GRB, SNeIa, and OHD data favor a dynamical dark energy model described by a scalar field with an exponential potential energy density. In Ref. [219], Amati, Ghirlanda, Yonetoku and Combo correlations have been calibrated in a model-independent way via OHD and jointly analyzed with SNeIa and BAO by using cosmographic methods, such as Taylor expansions, auxiliary variables and Padé approximations, to conclude that GRB do not favor the flat ΛCDM model but instead favor a mildly evolving dark energy density model. Similarly, in Ref. [220] the E p -E iso and Combo correlations have been calibrated via OHD actual and machine-learned data, and again, based on a joint analysis with SNeIa and BAO, indications against a genuine cosmological constant have been found. Analogously, in Ref. [221] different combinations of SNe Ia, quasar, and GRB data sets have been used for testing the ΛCDM model and dynamical dark energy parametrizations. It was found that GRB and quasar data sets were inconsistent with the flat ΛCDM model, in agreement with Ref. [222] for similar data. In Ref. [223] strong gravitational lensing data in conjunction with SNe Ia and GRBs have been considered and it has been found that the best-fit value of the spatial curvature parameter favored a closed universe, although a flat universe can be accommodated at the 68% confidence level.
On the other hand, some recent efforts have shown that the E p -E iso and Combo correlations calibrated using betterestablished cosmological data -such as SNe Ia or OHD measurements -provide cosmological constraints that are consistent with the flat ΛCDM model, see Table II. In Ref. [123] an updated E p -E iso correlation with 193 GRBs and a calibration based on an interpolation of the OHD data set have been considered, leading to Ω m0 = 0.397 +0.040 −0.039 in a flat ΛCDM cosmology, though the value of the mass density is higher than the one established by Ref. [182]. In Ref. [224] the E p -E iso correlation, calibrated with the latest OHD data set, has been jointly fit with CMB, BAO and SNe Ia data in a search for cosmological parameter constraints within the standard cosmological model, as well as in dynamical dark energy parametrizations, finding no evidence in favour of the alternatives to the ΛCDM model. Finally, by using the Combo correlation with 174 GRBs calibrated in a semi-model independent way, in Ref. [141] it has been found: a) for a flat ΛCDM model Ω m0 = 0.32 +0.05 −0.05 and Ω m0 = 0.22 +0.04 −0.03 for the two values of the Hubble constant H 0 of Ref. [182] and Ref. [225], respectively, and b) for a non-flat ΛCDM model Ω m0 = 0.34 +0.08 −0.07 and Ω Λ = 0.91 +0. 22 −0.35 for the H 0 of Ref. [182], and Ω m0 = 0.24 +0.06 −0.05 and Ω Λ = 1.01 +0.15 −0.25 for the H 0 of Ref. [225]. Again, by examining an uncalibrated E p -E iso correlation built up from a sample of bright F ermi-LAT GRBs [226] and another GRB sample with lower average fluence GRBs [227], in Ref. [176] cosmological parameter constraints have been obtained in a number of cosmological models, concluding that current GRB data are not able to restrictively constrain cosmological parameters, and that cosmological parameter constraints from the more-reliable GRBs are consistent with those resulting from better-established cosmological probes. In Ref. [189], a joint H(z)+BAO+quasar (QSO)+HII starburst galaxy (HIIG)+GRB fit determined Ω m0 = 0.313 ± 0.013 in the flat ΛCDM model, consistency with a cosmological constant and zero spatial curvature, though mild dark energy dynamics or a little spatial curvature are not ruled out at all.
Mixing all together, the cosmological results summarized above, obtained through GRB data, seem to be mutually inconsistent. This reflects all the efforts made so far to employ GRB as distance indicators are still affected by a certain number of issues, as we outlined previously.
First of all, we recall that GRB correlations involve a number of observable quantities affected by the so-called circularity problem [166], caused by having to compute the GRB correlations in a a priori assumed background cosmological model, being not fully model-independent [12,136,138,140,151,208,212,217]. However, even uncalibrated GRB correlations, in principle free from the circularity issue, are not able to put stringent constraints on the cosmological parameters, though consistent with those resulting from better-established cosmological probes. In addition, we recall all GRB correlations are characterized by large intrinsic dispersions, conceivably caused by unknown large systematic errors 27 [64,[143][144][145] in comparison to the case of better-established probes, such as BAO, OHD, and SNeIa, where many error sources have been better modeled. On the other hand, the influence of possible selection bias and evolution effects are currently debated [154][155][156][157][158]. One may therefore conclude that the large intrinsic dispersions of GRB correlations could be a consequence of yet undiscovered GRB intrinsic properties and/or a yet unidentified sub-class within the population of GRBs, analogously to SN populations.

B. Applications of statistical analysis with GRBs
In this section, we describe a few applications of statistical analysis using GRBs. Clearly, we focus on a particular choice and, in principle, it is possible to work out different fits and/or experimental procedures. In particular, we here propose a calibration at the very beginning, adopting the most consolidate route to handle GRBs. Above we also described the non-calibration procedure that, for brevity, we do not report here.
A widely consolidate approach is based on sampling the original catalog by means of a Monte Carlo technique, essentially built up using Markov Chain Monte Carlo simulations, that are sampled within the widest possible parameter space. Commonly, the most adopted algorithm is the Metropolis-Hastings and the standard approach to get limits uses the minimization of the total χ 2 function. As we stated above, in this review, the idea of combining more than one sample is essential in order to refine cosmological bounds on the model parameters. Hereafter we denote with x the set of parameters and include in our analysis SNe Ia and BAO data sets together with the calibrated GRB data. The former data have been obtained through calibrating the correlations. For the sake of brevity, as well as above, we only consider Amati, Ghirlanda, Yonetoku and Combo correlations. In the specific case of our three samples, i.e., GRBs, SNeIa and BAO, we combine the chi square functions by with the following recipe.
-GRB χ 2 . Here, we define where N GRB and µ th GRB are the experimental and theoretical GRB distance moduli.
-SN χ 2 . Here, by virtue of the above discussion concerning SN statistical analysis, we rewrite the chi square function in Eq. (25) by where ∆µ SN ≡ µ SN − µ th SN (x, z i ) is the module of the vector of residuals, and C the covariance matrix. In particular, we prompt the distance modulus for the most recent SN catalog, named Pantheon Sample. This represents the current largest SN sample consisting of 1048 SNe Ia lying on 0.01 < z < 2.3 [228]. The corresponding magnitudes read Here, M and m B are the B-band absolute and apparent magnitudes, respectively. The above distance moduli also depend upon other quantities required to standardize/correct the light curves of SNe Ia. The quantities X 1 and C are the light curve shape and color parameters, respectively, whereas α and β are the coefficients of the luminosity-stretch and luminosity-color relationships, respectively. ∆ M is a distance correction determined on host galaxy mass of SNe, while ∆ B a distance correction that is built up from predicted biases determined by means of simulations.
By marginalizing over M through a flat prior, it is possible to demonstrate that SN uncertainties do not depend on M and this permits one to simplify the chi square function through where -BAO χ 2 . The chi square function for BAO data is given in Eqs. (27)- (29).
Below we summarize a couple of statistical methods applied to GRB data to extract cosmological constraints.

Bézier polynomials and cosmographic series
The first method utilizes GRB data, calibrated through the above Bézier polynomials, to extract cosmological constraints by means of cosmographic model-independent series and heal de facto the circularity problem without postulating the model a priori [219]. This method has been implemented to Amati, Ghirlanda, Yonetoku and Combo GRB correlations in conjunction with SNe Ia and BAO data sets to get more stable and narrow constraints. We considered the most recent approaches to cosmography, comparing among them Taylor expansions with z and y 2 series, and Padé polynomials. Two hierarchies have been considered: hierarchy 1, up to j 0 , and hierarchy 2, up to s 0 .
Reasonable results have been found up for both hierarchies through several MCMC fits showing possible matching with the standard paradigm (see Tables III-VI). Moreover, we only partially alleviated the tension on local H 0 measurements as hierarchy 2 is considered. Taylor outcomes are quite stable within each hierarchy, as portrayed by the results in Table III, and work well with Amati, Ghirlanda and Yonetoku correlations in the sense that the corresponding numerical outcomes are consistent within 1-σ with previous findings. Again, this suggests a spatially flat ΛCDM paradigm as statistically favored model, with mass density parameter Ω m = 2(1 + q 0 )/3 ∼ 0.3 for Combo correlation, whereas the other correlations seem to indicate smaller values.
The auxiliary y 2 variable is not enough stable than Taylor expansions. It enlarges significantly h 0 , see e.g. Tab. IV and the overall results are however quite unpredictive at the level of hierarchy 1. Moreover, Padé fits seem to improve the quality of Taylor expansion hierarchy 1, as expected by construction. This is particularly evident for Combo and Yonetoku correlations, while for Amati and Ghirlanda correlations is not. It is worth noticing that to go further jerk term implies ≥ (3, 1), leading to higher orders than P 3,1 , quite unconstrained at higher redshift domains.
Quite surprisingly, our findings summarized in Tables III-VI show that the ΛCDM model is not fully confirmed using GRBs. Although this can be an indication that more refined analyses are needful, as GRBs are involved, simple indications seem to be against a genuine cosmological constant [182] and may be interpreted either with a barotropic dark energy contribution or with the need of non-zero spatial curvature [219]. Nevertheless, at this stage, our findings are in line with recent claims on tensions with the ΛCDM model [222,229,230].

Bézier polynomials and ΛCDM and ωCDM cosmological models
In a second method here summarized, the circularity problem affecting GRBs is again overcome by using Bézier polynomials to calibrate, in this case, the Amati correlation alone [123]. Unlike the previous method, now GRB data are utilized in conjuction with the SNe Ia JLA data set alone [231] and employed to explicitly constrain two different cosmological scenarios: the concordance ΛCDM model and the ωCDM model, with the dark energy equation of state parameter free to vary [123].
In the Monte Carlo integration, through the Metropolis-Hastings algorithm, H 0 has been fixed to the best-fit value obtained from the model-independent analysis over OHD data, i.e. H 0 = 67.76 km s −1 Mpc −1 . The results for Ω m and w (see Table VII) agree with previous findings making use of GRBs. The statistical performance of the models under study has been evaluated through the Akaike information criterion (AIC) criterion [232] AIC ≡ 2p − 2 ln L max , where p is the number of free parameters in the model and L max is the maximum probability function calculated at   the best-fit point, and the deviance information criterion (DIC) criterion [233] DIC ≡ 2p ef f − 2 ln L max , where p ef f = −2 ln L + 2 ln L max is the number of parameters that a data set can effectively constrain. 28 The best model is the one that minimizes the AIC and DIC values. Unlike the AIC criterion, the DIC statistics does not penalize for the total number of free parameters of the model, but only for those which are constrained by the data [234]. Differently from the previous approach, we found that the ΛCDM model is preferred with respect to the minimal ωCDM extension (see Table VII) and then conclude that no modifications of the standard paradigm are expected as intermediate redshifts are involved. However, future efforts dedicated to the use of our new technique to fix refined constraints over dynamical dark energy models are encouraged in order to fix the apparent dichotomy in the results of the two described methods.
C. The role of spatial curvature An update sample of GRBs has been developed in 2020, in which the Combo relation extracts bounds on the spatial curvature with no other probes [141], differently from previous attempts that commonly assume a spatially flat background, with the inclusion of SNe Ia and BAO.
The way in which the Combo relation is calibrated is without OHD data set, but rather invoking two step methods. In this picture, we assume [140] I. the terms k 1 and σ k are got from small GRB sub-samples with almost the same redshift; II. k 0 is determined from the use of SNe Ia limited to the lowest redshift of the GRBs of the Combo data set, in which the calibration of SNe Ia is negligible [140].
GRB sub-samples with the same z are chosen among those ones providing well constrained best-fit parameters. 28 Here, the brackets indicate the average over the posterior distribution.

169.23
Considering that in each sub-sample the GRB luminosity distances d L are quite the same, we employ the rest-frame 0.3-10 keV energy flux F 0 . This improves the dependence on the model and enables one to render our procedure cosmology-independent. In the specific example, here reported, seven sub-samples have been suggested and the best fits reported in Tab. VIII, showing no evident trends with z within the errors. This technique is used for Combo relation since previous results suggest its advantages in the fitting strategies above described 29 . One can perform a simultaneous fit of the above sub-samples with the same k 1 and σ k implying k 1 = 0.90 ± 0.13 and σ k = 0.28 ± 0.03.
The calibration of k 0 is performed by means of the nearest couple of GRBs of the employed sample with the same redshift. In particular, it is possible to take GRB 130702A at z = 0.145 and GRB 161219B at z = 0.1475. Then, µ obs C can be replaced via its average distance modulus µ SNIa = 39.21 ± 0.24, determined by SNe Ia with the same z of the above two GRBs, considering the bound over k 1 and the values of F 0 , E p , τ , and α for the two GRBs adopted throughout the computation. The computed value is k 0 = 49.54 ± 0.21.
At this stage, comparing between GRB distance moduli µ obs C , with uncertainties σµ obs C , with theoretical expectations, it is possible to get constraints over background cosmologies. In particular, for a non-flat ΛCDM model, i.e., the simplest scenario to work with, we write and we compute numerical bounds once H 0 is marginalized 30 as in Tab. IX. In particular, slightly larger estimations on matter density are got from GRBs, i.e., Ω m = 0.32 +0.05 −0.05 with H 0 of Ref. [122] and the opposite, i.e., Ω m = 0.22 +0.04 −0.03 , for the H 0 of Ref. [225]. Analogous results, i.e., compatible with the flat case, are computed using the non-flat ΛCDM model. 29 By construction, these sub-samples exhibit the same correlation with the same q 1 and σq, but involving different normalizations. 30 Constraints come from H 0 = (67.4 ± 0.5) km s −1 Mpc −1 [122] and H 0 = (74.03 ± 1.42) km s −1 Mpc −1 [225].

X. FURTHER APPLICATION OF GRBS AS PROBES OF THE HIGH-REDSHIFT UNIVERSE
Above we faced some the standard statistical methods associated with GRBs and we propose a few applications with the corresponding experimental bounds. Now we shortly review other phenomenological applications of GRBs that can shed light on the high-redshift Universe.

A. Star formation rate from GRBs
LGRBs are likely associated with core-collapse SNe [38,235]. By virtue of this fact, GRB-SN associations may provide a complementary technique that measures high-redshift star formation rate [236][237][238][239][240]. Again a problem related to calibration occurs, i.e. how to calibrate GRB event rate with star formation rate. Moreover, the luminosity function is not known a priori and its reconstruction depends upon the particular model selected for the analysis [241][242][243][244]. The typical functional structures for the luminosity functions are i) a broken power law [245], ii) a single power law with an exponential cut-off at low luminosities [242]. Reconstructed SFRs from GRBs are typically larger than those from other observations [246]. The reason behind this apparent inconsistency may reside in the fact that usually SFR at high-z is got from the observations of the brightest galaxies, whereas GRBs, in view of their high luminosities, may help in detecting faint galaxies at high-z otherwise unobserved [246,247].

B. High-redshift GRB rate excess
Although appealing, the above developments do not show why GRBs do not follow the star formation history being enhanced by hidden high-redshift mechanisms [247][248][249][250][251][252]. In particular, the star formation rate at high-redshift, namely z > 6, appears too large if confronted to star formation rate got from high-redshift galaxy surveys [253].
A natural origin of the high-redshift GRB rate excess can be found in the metallicity evolution, as LGRBs seem to prefer low-metallicity environment, as supported by recent studies that favor such a requirement 31 . Typical mass bounds on stars suggest > 30M , being responsible for BH remnants.

C. Gravitationally-lensed GRBs
Gravitationally lensed GRBs (GLGRBs) have been proposed in Ref. [49], where it was speculated that such a phenomenon would result in multiple light curves detected at different times, as due to the different light paths of the produced multiple GRB images.
Quests for GLGRBs were mostly based on strong lensing 32 and similarities among GRB light curves with identical spectra and close locations in the sky, as primary search citeria [256][257][258][259][260]. However, such searches led to null results, possibly due to Poisson noise uncertainties, affecting GRB light curves specially at low signal-to-noise ratios, which may have introduced large differences between the lensed GRB images [259]. On the other hand, some GLGRB could exhibit time delays shorter than (or comparable to) the burst duration, hence leading to unresolved (or locally separated) peaks separated by the time delay [261][262][263].
Several searches have been performed in the literature, resulting in a few or null candidates, based on different techniques and lens models, such as globular cluster with mass of ≈ Considering models where the lens is a supermassive BH [264,265], GLGRB candidates can provide an opportunity to estimate the number density of massive compact objects at cosmological distances by calculating the rate of GRB lensing. Assuming such supermassive BH lenses with mass ≈ 10 6 M as dark matter compact objects [266], the density parameter of BHs is Ω BH = 0.007 ± 0.004, or equivalently Ω BH /Ω M = 0.027 ± 0.016 [265]. In this respect, finding more GLGRBs candidates from supermassive BH may enhance our understanding of the matter content of the Universe.

XI. FINAL OUTLOOKS
In this review, we outlined some of the most recent developments toward GRB physics, properties and their applications to cosmology. In particular, the review is structured into two main parts.
In the first part, we discussed the basic demands of GRBs, including their main observable quantities, their classification scheme into LGRBs and SGRBs and the corresponding issues, emphasizing new possible ways of classification, still object of speculation. Afterwards, we gave emphasis on GRB progenitors and on their fundamental microphysics, in view of the experimental evidences characterizing prompt and afterglow emissions, etc.
LGRB connections with SNe have been explored as well along with SGRB matching with KNe and GWs. Great emphasis has been devoted to portray the standard GRB formation, working with the well consolidated fireball model. Particle acceleration and radiative processes, predicted in such a picture, have been largely reported, with particular concern on observable signatures and on the standard model frontiers. Even in this part, we illustrated that GRBs cannot be contemplated as genuine standard candles, since there is no consensus toward their internal processes that we depicted throughout the manuscript. Accordingly, their luminosity cannot be easily put in relation with their redshifts as, for instance, one does for SNe.
For these aspects and for the overall limitations above described, we developed in the second part considerable cosmological applications of GRB physics. We tried to standardize GRBs by means of the most recent techniques and accentuated GRBs are essential to reconcile small with intermediate redshift domains, opening new scenarios toward our universe comprehension. In this respect, we featured how GRBs could be used as complementary and outstanding probes to trace dark energy's evolution in support of other indicators, e.g. SNeIa, BAO, CMB, Hubble differential data, etc. Thereby we have shown a few statistical treatments related to Bayesian analysis in cosmology, able to combine GRBs with other catalogs of data, reporting the most recent cosmological constraints on dark energy models. To do so, we expounded the bristly circularity problem, burdening GRBs in cosmological set ups. In particular, we also changed perspective, showing how to avoid calibration, i.e., how not to employ the circularity. We confronted the two methods and checked which departures could be expected from the standard cosmological model through the use of GRBs in both the cases. Details on error propagation and GRB systematics have been discussed for several cosmic GRB correlations. Model dependent and independent techniques of calibrations have been likewise portrayed.
Perspectives about GRB developments will be based on clarifying the overall issues raised in this review. In particular, it is of utmost importance to shed light on how to standardize GRBs, in view of a likely self-consistent evolutionary paradigm, so far missing. With this recipe, we expect in the incoming years to improve GRB use in cosmology and get rid of circularity and greatly reduce the systematics and all the other issues that affect GRB data and challenge their use in cosmology. In particular, some models akin to those characterizing other cosmic indicators will spell out how to describe in toto GRB physics and evolution.