1 Introduction↩︎

The Hubble diagram (HD) is a plot of object redshifts \(z\) with respect to object distances \(d\). The HD is a well-known and widely-used practical cosmological test [1]–[3]. Based on cosmological models, one usually describe theoretical the redshift–distance relation as a parametric function \(d(z,\mathbf{p})\), where \(\mathbf{p}\) is a parameter vector. Thus, cosmology-independent determination of distances to objects with known redshift gives one a unique opportunity to verify, compare and probe parameters of cosmological models. In case of the standard \(\Lambda\)CDM cosmological model, for the HD \(\mathbf{p}\) may be considered, for example, as \((H_0, \Omega_M)\).

Usually, the HD is built up using the so-called standard candles (SC), i.e. the objects with theoretically or empirically known absolute brightness (or magnitude). Measured visible magnitudes of SCs directly give one required distances. At the end of the 20th century, the HD was constructed for Type Ia supernovae (SNe Ia) standard candles. Thus, the accelerated expansion of Universe within the standard Friedmann–Lemaitre–Robertson–Walker (FLRW) [2] cosmological model was discovered. This led to the introduction of dark energy into the standard cosmological model (SCM) [4], [5]. However, there is still a wide discussion on cosmological models and their parameter values [6]–[9]. The \(w\)CDM model, where the \(\Lambda\)CDM model is a special case of one, is often considered as an alternative to SCM. The \(w\)CDM model is defined as the FLRW model that contains two cosmological non-interacting fluids, having equation of state of the cold matter \(p = 0\), and the quintessence (dark energy) \(p = w\rho c^2\) (with \(w < 0\)) [2].

The current limit of SN observations is about \(z\sim{}\) 2–3, while the long gamma-ray bursts (LGRBs) up to \(z\approx 10\) and are already seen [10], [11]. This makes them promising objects to prolong the HD significantly further into the Universe. The LGRB sources are related to explosions of massive core-collapse SN in distant galaxies [12], though up to now there is no satisfactory theory of the LGRB radiation origins [13], [14]. There are a number of studies that suggest using LGRBs as SCs [15]–[25]. We suppose that LGRB HD can be used for probing cosmological parameters \(\mathbf{p}\) and comparing cosmological models, as well as the SNe HD [3], [26], [27].

The Amati relation is an observed linear-like correlation between LGRB spectrum parameters, including redshift and distance to LGRB host galaxy, in logarithmic plot [18]. Despite the fact that individual LGRBs may not satisfy the Amati relation due to both physical and observational factors, we expect that for a statistically large ensemble of LGRBs, the correlation is on average correct. So the relation gives one an statistical opportunity to measure the distance to LGRBs with known redshift independently of a cosmological model. However, it depends on two unknown parameters that have to be calibrated observationally. For calibration of the correlation, the LGRBs with known distances and redshifts are needed so that the first step of our study is in determining distances \(d\) for a subsample of LGRBs in near galaxies by cosmology-independent methods.

In this study, we try to calibrate the Swift LGRBs ¹ as standard candles by using \(\Lambda\)CDM as a basis, i.e. obtaining distances from the model. However, this approach involves the circularity problem [28] as for proper calibration of cosmologically-independent distance measurement should be used. So next we calibrate LGRBs as standard candles by using the Pantheon SNe Ia [29] as a basis.

We decided to approximate the SN HD by using a smooth elementary mathematical function \(d^{SN}(z)\) that can be directly used to get the distances to sources of LGRBs with known \(z\). The LGRBs should be near enough, so they can lie inside representative SN sample, where the distance error \(\sigma_{d^{SN}(z)}\) is small enough. Corresponding cosmological background of the Amati relation and introduction in mathematical approach of its calibration was described in [3]: in this paper, we apply the proposed idea of a cosmological model-independent calibration of gamma-ray bursts using a modified formula and advanced statistical methods.

The second step is in determining the Amati relation parameters via selected near LGRBs (i.e. LGRB calibrating), finding distances to all LGRBs with known redshift via the calibrated Amati relation, and plotting the HD for them. In fact, SNe are also calibrated SCs via cepheids that are also calibrated (via parallaxes) so that our approach lies in the frameworks of extending the cosmological distance scale. Since the supernova catalogues data are tied to a fixed value of the Hubble constant \(H_0=70\) km/s/Mpc, this approach is not completely cosmologically-independent, all values are obtained with an accuracy of the scale factor \(H_0\) [30].

At almost every stage of this study it is required to find best-fitting parameter estimation. In case, e.g., of the Amati parameters, the best-fitting function \(y=f(x)\) should take into account errors in data points both of \(x\) and \(y\). Because of this and lack of LGRB statistics we find common linear least-square method to be unsuitable for our study. Instead we use the Theil-Sen estimator (also known as the single median method) for the linear regression purposes [31]. In this method, the estimation for the slope is defined as the median slope among all of the possible pairs of dots. We also need a method to perform the parameter estimation (the curve fitting) routines for the arbitrary function \(f\). For this we use the trust region reflective algorithm [32], which is the non-linear least-squares method implemented in the curve_fit function of the scipy.optimize Python library [33]. We can interpret both of these methods as pipelines that return the estimated parameters based on the input data. Our purpose is to also estimate the errors and covariances of the parameters. To do this, we utilise the Monte-Carlo sampling uncertainty propagation approach [34], [35]. Applied to our pipelines, this approach allows us to correctly propagate the errors directly from the input data to the estimated parameters.

In the paper we have calibrated the LGRBs from the selected subsample as SCs via supernovae, as it is described above, and have built the corresponding LGRB HD, which is the main result of this work. Nowadays, there are still not enough statistics and theory to get reliable values of cosmological parameters from the LGRB HD. However, we find our method perspective and promising one. We also compare \(\Lambda\)CDM-based and SN-based Amati parameters and find no significant difference between them.

2 Materials and Methods↩︎

2.1 The Hubble diagram as a basic cosmological test↩︎

The HD is the redshift–distance relation, which can be represented as a table of points \((z_i,d_i)\) on the plot for each \(i\)-th object. Usually, the luminosity distance \(d_L\) is used. This value depends on the visible magnitude \(m\), absolute magnitude \(M\), and distance modulus \(\mu\) as \[\label{mu}m-M=\mu=5\log (d_L / 1\,\mathrm{Mpc})+25 \, .\tag{1}\] In fact, a HD is plotting the \((z_i,\mu_i)\) relation.

In general, cosmological models give theoretical assumptions on the function \(d_L(z, \mathbf{p})\), where \(\mathbf{p}\) are model parameters. For example, in \(w\)CMD model \(\mathbf{p} = \{H_0,\,w,\,\Omega_w,\,\Omega_k\}\), and in \(\Lambda\)CMD model \(\mathbf{p} = \{H_0,\,w=-1,\,\Omega_w=\Omega_\Lambda,\,\Omega_k=0.0\}\). The HD is one of the cosmological tests that allow one to find the model parameters via fitting the observational data. For this one, the redshift \(z\) and distance \(d_L\) should be determined independently (e.g., \(z\) directly via spectrum, and \(d_L\) via various indirect methods).

2.2 Gamma-ray bursts as standard candles↩︎

Standard candles are a group of objects with the known (theoretically or empirically) typical absolute magnitude \(M\). They allow one to find luminosity distance or distance modulus directly from Equation (1 ) by using the visible magnitude \(m\), or bolometric fluxes.

In case of LGRBs, the situation is a little different. The Amati relation for LGRBs [15], [18], [19] is the equation \[\log E_\mathrm{iso} = a\log E_{\mathrm{p},i}+b \, ,\] where

\(E_{\mathrm{p},i}=E_\mathrm{p}(1+z)\) is the rest frame spectral peak energy, where \(E_\mathrm{p}\) is the observed LGRB parameter determined via fitting the LGRB spectrum by the cut-off power law (CPL) model function \(N(E)\), \[N(E)\left[\frac{\mathrm{photons}}{\mathrm{KeV}\,\mathrm{s}\,\mathrm{cm}^2}\right] = AE^\alpha \exp \big( -(2 + \alpha)E/ E_\mathrm{p} \big) \, ,\] where \(\alpha\) and \(E_\mathrm{p}\) are the CPL model parameters of a LGRB from the Swift database;

\(E_\mathrm{iso}=4\pi d_L^2 \cdot S_\mathrm{bolo} / (1+z)\) is the isotropic equivalent radiated energy in gamma-rays [36]. The distance \(d_L\) and observed integral fluence \(S_\mathrm{bolo}\) are determined as quantities transferred per a unit energy frame area and that are corrected for the instrumental (observed) spectral energy range, and source redshift. The correction is performed by the equation \[S_\mathrm{bolo}= S_\mathrm{obs} \dfrac{\int_\frac{1}{1+z}^\frac{10^4}{1+z} EN(E)\;\mathrm{d}E }{\int_{E_\mathrm{min}}^{E_\mathrm{max}} E N(E)\;\mathrm{d}E} \, ,\] where \(S_\mathrm{obs}\) is the observed fluence, and \(\{E_\mathrm{min},\,E_\mathrm{max}\}\) is the instrumental spectral energy range, which is \(\{15, 150\}\,\mathrm{keV}\) for the Swift observatory; normalised example of such spectra is shown in Fig. 1;

\(a\) and \(b\) are the Amati relation parameters mentioned above that can be calibrated empirically as in this study.

If observed parameters \(z\), \(S_\mathrm{bolo}\), \(E_{\mathrm{p},i}\), and model-independent measurement of \(d_L(z)\) (at the small scales \(z\lesssim 1.5\)) are known, then best fit of the Amati relation parameters \(a\) and \(b\) can be found. Using \(a\), \(b\), \(z\), \(S_\mathrm{bolo}\), and \(E_{\mathrm{p},i}\), one can find the LGRB distances \(d^\mathrm{LGRB}_L(z)\) via the Amati relation and plot the LGRB HD.

So, to find \(a\) and \(b\), one needs to get luminosity distances \(d_L\) for near LGRBs. Although taking \(d_L(z)\) from \(\Lambda\)CMD model or from any other cosmological model is a dead-loop (a circularity problem), we have first calibrated \(a\) and \(b\) using \(\Lambda\)CMD model using parameters from [6] as a comparison basis.

Then we have chosen the way of the SN calibration. Thus, we attempted to build a function \(d_L^\mathrm{SN}(z)\) by using SNe Ia as SCs. Distance \(d_L\) was calculated as function \(d_L^\mathrm{SN}(z)\) for near LGRBs that have redshifts \(z\) low enough in order to perform an interpolation of the SN HD.

This lets one to directly get the Amati relation parameters and to use them for cosmological model-independent determination of distances \(d_L\) to all LGRBs with known observational data for \(z\), \(S_\mathrm{bolo}\), and \(E_{\mathrm{p},i}\) and to plot the HD for them.

2.3 Catalogues of SNe Ia and LGRBs↩︎

We have used the Pantheon database (containing 1,048 SNe). In this catalogue, there are two kinds of redshift, which are the redshift in terms of the cosmic microwave background radiation \(z_\mathrm{cmb}\), and the redshift in terms of the heliocentric system \(z_\mathrm{hel}\) [29]. Between these two methods, we have chosen the \(z_\mathrm{cmb}\) one. All the catalogue data is calculated in assumption of \(H_0=70\) km/s/Mpc, so we have set this value for all our study.

For our purposes, we use the sample of 174 GRBs from the Swift catalogue with the measured parameters of \((\alpha,\,E_{\mathrm{p}},\,S_{\mathrm{obs}},\,z)\), where the spectral parameter \(\alpha\), the peak energy \(E_{\mathrm{p}}\), and the observed flux \(S_{\mathrm{obs}}\) are described above, and \(z\) is the redshift. Some of the uncertainties of these parameters were missing (several \(\alpha\) errors and approximately a third of \(E_{\mathrm{p}}\) errors were not presented). To deal with this issue, we assume the points with unknown errors to have the relative uncertainties equal to the median relative uncertainties of the corresponding parameter of the sample. The redshift \(z\) values are also presented in the lack of uncertainties, as is expected. A significant proportion of the LGRB parameters has asymmetric and relatively huge errors, which leads to the inexpediency of using the common uncertainty propagation rule. Instead, we use the approach of Monte-Carlo sampling, which is thoroughly described in the Section 2.4. Our LGRB sample is represented in Google Table².

2.4 Monte-Carlo uncertainty propagation↩︎

The standard approach to the error propagation problem is known as the linear uncertainty propagation (LUP) theory. One of the best implemented error propagation software is the uncertainties package of the Python language. The errors in this approach are interpreted as the standard deviations, while the values have the sense of the mean. So, the values of variable with its error defines by the normal distribution, which denotes the possible whereabouts of the variable via distribution parameters. However, the LUP theory requires that the errors be relatively small compared to the values of variable, so that the functions in calculations are nearly linear compared to these small shifts. Also, the LUP approach is not capable of handling asymmetric errors. Thus, we cannot use it.

For that reason, we decided to use the Monte-Carlo sampling. In this approach, the values and their errors also have the sense of defining the distributions of possible locations, and we use these distributions to draw samples of the size of 10,000. To handle the asymmetric uncertainties, we interpret the values as the medians, the lower bounds are interpreted as the 0.16 quantiles, and the upper bounds as the 0.84 quantiles. With this interpretation of errors, the case of symmetric uncertainties reduces to the normal distribution with the known mean and standard deviation. In the case of asymmetric errors, we use the split-normal distribution, which results from joining the two halves of normal distributions with different standard deviations at their mode. Some of the \(E_{\mathrm{p}}\) values from our LGRB data set had remarkably huge lower uncertainties, so that trying to Monte-Carlo sample them would lead to negative \(E_{\mathrm{p}}\) values, which has absolutely no physical sense. Because of that, in the case of \(E_{\mathrm{p}}\) variables, we draw the split-normal distributions in the space of \(\log E_{\mathrm{p}}\). So, the peak energy values are drawn using the log-split-normal distributions. Taking the logarithms does not move the quantiles, so the errors in this case are not changed.

The Monte-Carlo approach to propagating the errors is simple yet powerful. It allows one not only to calculate the uncertainties of calculated values, but also to track and take into account the correlations between variables for free.

2.5 Best-fitting methods↩︎

To obtain the best-fit parameters of the approximation function \(d_L^{SN}(z)\) for the SN HD, one needs to minimise the functional value \[\label{chi} \chi^2=\sum_{i=1}^n \frac{1}{\sigma_i^2} \bigl(y_i-f(x_i,\mathbf{p})^2\bigr) \, ,\tag{2}\] where

• \((x_i,y_i)\) are observed table values;

• \(\sigma_i\) is error of \(i\)-th value;

• \(f\) is the model function and \(\mathbf{p} = \{p_1,p_2,\dots,p_m\}\) are the parameters.

This function utilizes the trust region reflective algorithm [32] to minimize the \(\chi^2\) function.

2.6 Interpolation function of the SN HD↩︎

In fact, we can use any smooth function to use it as a \(d_L^{SN}(z)\) function. All that we need in order to minimize the error between the real value and a mathematically predicted one for the luminous distance at any low \(z\), where we have enough supernovae to do this, is the approximation accuracy that is provided by our approach in the redshifts \(z_{min} \ll z\lesssim z_{max}\), where \(z_{min}\) and \(z_{max}\) are redshifts of the nearest and the farthest supernovae (\(z=0.01012\) and \(z\approx 2.26\), respectively).

Since in logarithmic scales (along both of the distance and redshift) the relation \(d_L(z)\) is already known, it has almost linear behaviour at low \(z\), we have decided to use a polylogarithmic function of a degree \(k\). We have tried to use three following functions:

• theoretically-inspired function \[\label{tj}\mu^{SN}(z) = 5\log\frac{cz}{H_0}+25+\sum\limits_{i=1}^p a_i\log^i(1+z);\tag{3}\]

• simple polylogarithmic function \[\label{dlsn}\mu^{SN}(z) = \sum\limits_{i=0}^p a_i\log^iz;\tag{4}\]

• shifted polylogarithmic function \[\mu^{SN}(z) = \sum\limits_{i=0}^p a_i\log^i(1+z).\]

The first function may be considered as addition of small correction to the linear Hubble law at low \(z\). The results of using other function is shown in Fig. 2.

3 Results↩︎

3.1 Approximation of the SN HD↩︎

The polynomials given by equation (4 ) for degrees \(p=1,\,2,\,3,\,4,\,5\) have been used. The result of best-fitting is shown in figure 3.

Usually, polynomials of degree 1 or degree 2 are used for HD approximation [19]. We have chosen the model with \(p=2\), because it minimises the Akaike Information Criterion (AIC) [37]. We use the obtained estimates of the parameters and their covariance matrix to generate them from a multivariate normal distribution of the size of 10,000. The corner plot of this sample is shown in the figure 4. This sample will come useful later for the Monte-Carlo uncertainty propagation. These polynomial coefficients have no physical sense, they just show the best-fit of a cosmology received by model-independent redshift–distance relation.

3.2 Amati relation parameters probing and Gamma-Ray Bursts Hubble Diagram↩︎

By analysing figure 3, we took LGRBs with \(z<1.4\) to calibrate the Amati coefficients \(a\), and \(b\). All SNe samples in this range are representative enough and their approximations via polylogarythmic function have low values of its formal errors. In total, 75 of 174 LGRBs in this range (with known \(z\) and calculable \(E_{\text{p},i}\) and \(S_\text{bolo}\)) are available. To estimate the Amati parameters \(a\) and \(b\), we use the Theil-Sen estimation, which is a reliable and robust method for linear regression. In this method, the slope (parameter \(a\)) is estimated as the median of all of the slopes between all pairs of points. The intersect (parameter \(b\)) is then estimated as the median of the values \(y_i - ax_i\). Thus, we have the pipeline that takes the LGRBs data table and returns the Amati coefficients \(a\) and \(b\). The Monte-Carlo error propagation approach can be also applied to this pipeline, so that the output would consist not only of the estimated parameters of \(a\) and \(b\), but their Monte-Carlo samples of the size of 10,000. We can further use the medians of this samples as the estimated values for the parameters, and the quantiles of 0.16 and 0.84 as the upper and lower 1\(\sigma\)-borders.

The results of the estimation process are presented in the figure 5. "Repeated Theil-Sen estimation" means that the figure shows the average regression of 10,000 ones. In Table 1, we compare the Amati parameters, calculated through \(d_L^{\Lambda CDM}\) and \(d_L^{SN}\), where \(d_L^{\Lambda CDM}\) calculated by equation (B5) from [3]. With the estimated Amati parameters \(a\) and \(b\), it is now possible to find the distance modulus for the whole sample of LGRBs. The final SN+LGRB HD is shown in Fig. 6.

4 Discussion and conclusion↩︎

We have calibrated the near Swift LGRB sample up to \(z<1.4\) as standard candles through the Pantheon SN Ia catalogue using the Amati relation. The calibrated HD up to \(z\sim 7-8\) is shown in figure 6. This diagram is constructed naturally from the standard cosmological model at fixed value of \(H_0=70\) km/s/Mpc, and it does not take into account any systematic corrections.

The Amati relation for the near LGRBs calibrated by SN Ia and for all LGRBs calibrated by the \(\Lambda\)CDM model (Table 1) match with \(1\sigma\) level. It can be concluded that, in terms of statistic significance, there are no observed deviations from the standard cosmological model for the far LGRBs with \(z>1.4\).

However, visual analysis of the HD in the figure 6 shows a trend towards the fogging of distant LGRBs relatively to the standard model (shifting to the linear Hubble law). Possible explanations include LGRB evolution [38], [39]; gravitational lensing [40] and gravitational mesolensing that leads into increasing the brightness of objects in all spectral ranges[41]–[43], observational selection, and others. It should also be noted that the cosmic time dilation effect is taken into account in this study. This implies that our results are valid within models that include this effect. These issues require additional research, including the accumulation of the LGRB sample and the development of statistical methods.

The GRB observations in multimessenger astronomy epoch open new possibilities for testing the fundamental physics lying in the basis of the standard cosmological model: classical general relativity, cosmological principle of matter homogeneity, and the Lemaıtre space expansion nature of cosmological redshift.

The corrected HD cosmological tests can probe strong-field regime of gravitation theory, spatial distribution of galaxies, Hubble Law and time dilation of physical processes at high redshifts. Constructing the high-redshift GRB Hubble diagram and comparison of time dilation in GRB pulses, GRB afterglow and core-collapse SN light curves to test the expanding space paradigm [44]. Perspectives for performing cosmological tests in multimessenger astronomical observations of GRBs were considered and several new tests were proposed in [26], [27]. The LGRB HD is the test that could be used, in particular, to gravitational lensing and Malmquist biases testing [3], classical general relativity, cosmological principle of matter homogeneity, and the Lemaitre space expansion nature of cosmological redshift testing. The GRB HD can be combined with the gravitational wave standard sirens at intermediate redshifts [45]–[47]. We find the method of LGRBs HD construction proposed in this paper perspective for performing such tests in future.

Author contributions↩︎

Conceptualization, Nikita Lovyagin and Rustam Gainutdinov; Formal analysis, Rustam Gainutdinov; Funding acquisition, Stanislav Shirokov; Investigation, Stanislav Shirokov; Methodology, Nikita Lovyagin, Rustam Gainutdinov, Stanislav Shirokov and Vladimir Gorokhov; Project administration, Stanislav Shirokov; Resources, Rustam Gainutdinov; Software, Rustam Gainutdinov; Supervision, Stanislav Shirokov and Vladimir Gorokhov; Validation, Nikita Lovyagin and Stanislav Shirokov; Visualization, Rustam Gainutdinov; Writing – original draft, Nikita Lovyagin; Writing – review & editing, Vladimir Gorokhov.

Funding↩︎

Part of the observational data was exposured on the unique scientific facility the Big Telescope Alt-azimuthal SAO RAS and the data processing was supported under the Ministry of Science and Higher Education of the Russian Federation grant 075-15-2022-262 (13.MNPMU.21.0003)

Data availability↩︎

The codes developed in Python underlying this article are available in the repository on https://github.com/Roustique/sngrb

Acknowledgments↩︎

We thank the anonymous reviewer for important suggestions that helped us to improve the presentation of our results.

Conflict of interest↩︎

The authors declare no conflict of interest.

Abbreviations↩︎

The following abbreviations are used in this manuscript:

References↩︎