The Hubble Diagram: Jump from Supernovae to Gamma-Ray Bursts

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

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).

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 (Amati et al., 2002; Amati et al., 2008, 2019) is the equation

where

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 Aghanim et al. (2020) 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.

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}}$ (Scolnic et al., 2018). 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²²2https://docs.google.com/spreadsheets/d/1jbOQxUlwEYd8qWU68mbih_PBPOHp0iUs-3kgphUIvAs/edit?usp=sharing.

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.

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

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 (Byrd et al., 1987) to minimize the $\chi^{2}$ function.

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

  $$\mu^{SN}(z)=5\log\frac{cz}{H_{0}}+25+\sum\limits_{i=1}^{p}a_{i}\log^{i}(1+z);$$

  (4)

• •

  simple polylogarithmic function

  $$\mu^{SN}(z)=\sum\limits_{i=0}^{p}a_{i}\log^{i}z;$$

  (5)

• •

  shifted polylogarithmic function

  $$\mu^{SN}(z)=\sum\limits_{i=0}^{p}a_{i}\log^{i}(1+z).$$

  (6)

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.

The polynomials given by equation (5) 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 (e.g., Amati et al., 2019). We have chosen the model with $p=2$, because it minimises the Akaike Information Criterion (AIC) (Akaike, 1974). 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.

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 Shirokov et al. (2020b). 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.

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.

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)

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

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

The authors declare no conflict of interest.

The following abbreviations are used in this manuscript: