Solving the inverse cosmological calibration problem of gamma-ray bursts

We have received a new physical characteristics fitting based on actual observational data from the Swift mission's long-duration gamma-ray bursts (LGRBs). We considered such characteristics as the Amati parameters for linear correlation ($E_\text{iso}-E_{\text{p},i}$) and the $k$-correction for gravitational lensing and Malmquist bias (GLMB) effect. We used the Pantheon SN Ia catalogue and the standard $\Lambda$CDM model with a fixed Hubble constant of $H_0=70$ km/s/Mpc as the baseline for the Hubble function $\mu(z)$. In our paper, we formulated the inverse cosmological calibration problem (ICCP) in the non-parametric statistics framework. The ICCP involves fitting non-observable physical characteristics while assuming a fixed cosmological model. To solve this problem, we developed a new method that is resistant to non-Gaussian processes. This method is based on error propagation through the Monte Carlo method and the Theil-Sen method for linear regression estimate. We have demonstrated the stability and robustness of this assessment method. The parameter estimates are as follows: $a=0.92^{+0.12}_{-0.12}$, $b=50.32^{+0.33}_{-0.32}$ without considering the GLMB effect, and $a=0.63^{+0.13}_{-0.14}$, $b=50.12^{+0.33}_{-0.31}$, and $k=1.98^{+0.25}_{-0.24}$ with the effect included. The proposed method can be applied to any other calibration sample of known standard candles, a calibrated sample of LGRBs, and the Hubble function $\mu(z)$. In the future, the ICCP idea can be used as an alternative cosmological test for estimating cosmological parameters, including the GLMB effect, or even for the selection of models, providing new information about the Universe. This can be done by analysing the residual values of observational data within the Bayesian statistics paradigm.


INTRODUCTION
Long gamma-ray bursts (LGRB) are large-scale explosive events that occur in distant galaxies and are observed in the highest energy part of the electromagnetic spectrum.The LGRB sources are related to explosions of massive core-collapse SN in distant galaxies (Cano et al. 2017), though up to now there is no satisfactory theory of the LGRB radiation origins (Willingale & Mészáros 2017).They have a typical duration of a few seconds, although they can also be as long as milliseconds up to an hour.The idea of using LGRBs as standard candles (SCs) in cosmology is actively discussed in the last two decades (Amati et al. 2002(Amati et al. , 2008a(Amati et al. , 2019;;Demianski et al. 2017a,b;Lusso et al. 2019;Yonetoku et al. 2004;Wang et al. 2011;Wei & Wu 2017).Since LGRBs are extremely distant observable objects, using them as standard candles makes it possible to significantly extend the scale of cosmological distances.This makes it possible to use LGRBs to determine the parameters of cosmological models.One of the main cosmological tests is the distance-redshift diagram, also known as the Hubble diagram (HD).Today, the HD remains a relevant cosmological test, with which one can test various hypotheses or estimate cosmological parameters.Based on cosmological models, one usually describes theoretical the distance-redshift relation as a ★ E-mail: arhath.sis@yandex.ru† E-mail: roustique.g@gmail.com‡ E-mail: n.lovyagin@spbu.ruparametric function  (, p), where 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 (Baryshev & Teerikorpi 2012).
The HD is a well-known and widely-used practical cosmological test (Sandage 1997;Baryshev & Teerikorpi 2012;Shirokov et al. 2020c;Shirokov & Baryshev 2020).The first use of it can be considered as the discovery of the expansion of the Universe in 1929, described in the classic work of Edwin Hubble (Hubble 1929).The linearity of this Hubble law on the scales 1-300 Mpc was confirmed by Sandage at the Hale Telescope (current results are presented in Sandage, Reindl & Tammann 2010).At the end of the 20th century, the HD constructed for type Ia supernovae made it possible to discover the accelerated expansion of the Universe within the framework of the Friedmann-Lemaitre-Robertson-Walker metric.This led to the introduction of dark energy into the standard cosmological model (SCM) (Riess et al. 1998;Perlmutter et al. 1999).However, there is still a wide discussion on cosmological models and their parameter values (Aghanim et al. 2020;Riess et al. 2018;Riess 2020;Yershov et al. 2020).
Supernovae of the Ia type are classical observable standard candles, however, with current observable facilities, the limit redshifts for them of the order of  ∼ 2 ÷ 3.While LGRBs were watching up to  ∼ 10 (Amati 2018).Therefore, the inclusion of LGRBs in the scale of cosmological distances allows one to construct a HD that includes a deeper region of the Universe.Also, one can say that LGRBs are between supernovae and the cosmic microwave background radiation (CMBR), which is associated with probing cosmological parameters of very far ( ∼ 1 000), or very early Universe.Thus, the LGRB HD could be used as cosmological test that is a kind of "connecting link" between the SN and CMBR cosmology (Shirokov et al. 2020a).
We suppose that LGRB HD can be used as probing cosmological parameters p and comparing cosmological models, as well as the SNe HD.However, because of large variance and uncertainties of LGRB parameters, the LGRBs can be use as SCs in a mean sense only.The mean sense implies that we must use averaging characteristics (e.g., quantiles, or function fittings) for determining cosmological parameters.In order for LGRBs to be considered as SC, a correlation between the energy of the gamma fluence and the distance to LGRBs is required.One such correlation is the linear dependence of the isotropic equivalent radiated energy  iso on the peak energy of the spectrum  p, (Amati et al. 2008b), also called the Amati relation.The relation give one an opportunity to cosmology-independent measuring the distance by LGRBs with known redshift.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  for subsample of LGRBs in near galaxies by cosmology-independent methods.The Amati relation is a rather rough correlation, since, for example, it does not take into account the angle of the jet cone, under which we see this burst.A more accurate correlation that takes this angle into account is called the Ghirlanda relation (Ghirlanda et al. 2004(Ghirlanda et al. , 2007)).The main problem in using LGRBs for cosmological tests is getting rid of the dependence of the distances measured from the observed parameters on the choice of the cosmological model (Kodama et al. 2008).Due to this dependence, the measured distances to the LGRBs require additional calibration.
In our article Lovyagin et al. (2022), we have showed a method for calibrating LGRBs with  < 1.4 for type Ia supernovae from the Pantheon catalogue (Scolnic et al. 2018).The method is consisted in the fact that the distance modulus  -redshift  dependence for supernovae was approximated, and then used to calculate the LGRB distance moduli in the log  iso -log  p, plane.Next, using the estimated Amati parameters, we have plotted a HD for our LGRB sample from the Swift catalogue (Gehrels et al. 2004).In this paper, the LGRB calibration by supernovae is carried out in a different way.We stick to the Bayesian statistics as a way to getting new information about Universe.
Solving inverse cosmological calibration problem (ICCP) in the frameworks of fixed cosmological model, we can get fitting physical characteristics based on observational data.In case of LGRB these characteristics are the Amati parameters of linear correlation  iso −  p, (Amati et al. 2008a) and the -correction for gravitational lensing and Malmquist bias (GLMB) effect (Shirokov et al. 2020c).We try to estimate these characteristics through fitting the Swift LGRB data by using the SCM as a basis.
At the first stage, the cosmological parameters of the ΛCDM model are estimated from the HD for supernovae.The cosmological model with these parameters is then considered as fixed.Next, we estimate the parameters ,  of the linear Amati relation and the parameter .Thus, we find the optimal values of the parameters , , , i.e. such that the LGRB distance moduli, calculated by using these parameters, are corresponded to the ΛCDM model with parameters, estimated from type Ia supernovae from the Pantheon catalogue.

STANDARD COSMOLOGICAL MODEL
At the moment, the ΛCDM model with parameters  0 = 70, Ω Λ = 0.7 is accepted as the standard cosmological model (SCM), which is a special case of the Friedmann model, a non-stationary model of the Universe that satisfies the Einstein equations, which in a simplified form can be written as follows: where   is the Einstein tensor, Λ is the cosmological constant,   is the metric tensor,   is the energy-momentum tensor, and  = 8/ 4 is a constant.In the ΛCDM model, it is assumed that the Universe, in addition to baryonic matter, is also filled with cold dark matter (CDM) inaccessible to direct observation, as well as dark energy, which is related to the cosmological constant Λ from the Einstein equations 1.The content of matter (visible and dark combined) is described by such a parameter as the density divided by the critical density: For dark energy, the relative density is defined as follows: In addition to these two parameters, the curvature contribution is also introduced1 Ω k : Space-time curvature occurs if the average density of the universe differs from the critical one: Thus, there are 4 quantities to parameterise the ΛCDM model: •  0 -Hubble constant, a linear coefficient in Hubble's law that relates the distance to an extragalactic object and the speed of its receding (or approaching).Standard units: km/s Mpc ; • Ω m -density of matter (visible together with dark matter) reduced to critical density.Dimensionless value; • Ω Λ -density of dark energy reduced to critical density.Dimensionless value; • Ω k -the contribution of curvature to the total average density, Ω k = Ω m + Ω Λ − 1 (in most papers it is defined with the opposite sign).Dimensionless value.
Of these parameters, three (relative densities) are related by Equation (4), so this model actually has three degrees of freedom, not four.But usually, for the sake of clarity, the above four quantities are written out in their entirety.In this case, in the SCM framework Ω k = 0, and thus, there are two degrees of freedom.These parameters are included in the equations (5,6,7,8) relating redshift to distance.Accordingly, in order to test any cosmological model (in particular, to estimate the parameters of the ΛCDM model), we need to build a cosmological test, e.g., a HD (the distance -redshift dependence) obtained from observations.Expression for metric distance: where sn We will also use quantities such as energy distance and luminosity distance One of the main quantities, considered in this work, is the luminosity distance modulus The HD for some observed objects we will call the set of points (  ,   ).Having constructed a HD based on some catalogue of observed objects, we can estimate the cosmological parameters ( 0 , Ω m , Ω Λ , Ω k ) using approximation procedures.

OBSERVATIONAL DATA
As observational data, we must choose such objects for which the redshift and the distance modulus are measured independently of each other.Examples of such objects are type Ia supernovae, which are standard candles in cosmology.However, it is not supernovae that are of greatest interest to us, but LGRBs: in this paper, we consider the possibility of using the latter as standard candles.For further analysis, we need the initial LGRB observational data, from which we will search for their distance moduli.

Catalogue of supernovae used
We have used the Ia Pantheon supernova catalogue.It contains distance moduli and their errors, as well as redshifts for 1 048 supernovae2 .The redshift is determined in two different ways: in the heliocentric system  hel and taking into account the motion of Galaxy relative to the cosmic microwave background radiation  cmb .In our work, it was the latter values that were used.

Catalogue of gamma-ray bursts used
We have used data from 174 LGRBs observed by the Swift space observatory3 .For each LGRB we used the following four values.
•  is redshift.It is determined by identifying a LGRB with its host galaxy.The (spectroscopic) redshift error is absent in the catalogue.It is considered as negligibly small, so we do not use it in our calculations.
•  obs is the LGRB observable fluence in 10 −7 erg/cm 2 .The catalogue for these values also presents the errors corresponding to the confidence interval with the 90% error level (approximately 1.6).It is important for us that the errors must correspond to the 1 level, so these errors were divided by 1.6.
•  p is the LGRB peak energy in keV.This is the energy value at which the  2  () spectrum peak is observed.The peak energy errors are asymmetric, i.e. the upper and lower values are different.Moreover, these errors are only determined for approximately 1/3 of LGRBs, and they are missed for the remaining 2/3 of LGRBs.
•  is the dimensionless parameter to describe shape of the LGRB spectrum according to the CPL (cut-off power law) model, given by Equation ( 9).The errors for this parameter are also asymmetric, and they are given for almost all gamma-ray bursts.
As mentioned above, the errors of some quantities are not presented for all available LGRBs in our sample.In order to take into account points without errors along with points for which errors are presented, we used the following solution: all unknown errors were replaced with median relative errors for a given quantity in our sample.The original and enhanced catalogues can be found in Appendices B and C, respectively.

Photometric quantities of LGRBs
These formulas (with the exception of the CPL spectrum model) are described in more detail in Shirokov et al. (2020c).The CPL spectrum model is defined as follows: where  is a normalisation parameter.This model is used for converting from the observed fluence to the standard "bolometric" gamma fluence (from 1 to 10 4 keV) in the integration limits, the energies are given in keV.The integration limits in the denominator correspond to the spectral range of the main telescope (BAT) of the Swift space observatory (from 15 to 150 keV).
We also want to introduce a phenomenological model for taking into account the Malmquist effect and gravitational lensing in accordance with Shirokov et al. (2020c).It is carried out in the form of a parameterisation using a following power law: The parameter  characterises the change in the bolometric fluence due to the Malmquist effect and gravitational lensing.Knowing the bolometric fluence, we can find the isotropic equivalent radiated energy This value has the following meaning: it is the energy that would be emitted by a LGRB if it emitted energy isotropically (with a spherical indicatrix).
Another quantity we need is the peak energy in the source frame.It is equal to (13)

The Amati relation and the extended Amati relation
Lorenzo Amati was one of the first who discovered the correlation between the LGRB spectrum peak energy in the source frame  p, , given by Equation ( 13), and the LGRB isotropic equivalent radiated energy  iso , given by Equation ( 12) (Amati, L. et al. 2002).This dependency is called the Amati relation and can be written in this form: The coefficients  and  are called Amati parameters.Equation ( 14) can be rewritten as follows: Thus, for each LGRB, the ratio  iso / p,  is the same with the same accuracy as the Amati relation ( 14).Based on the presence of such a constant quantity for LGRBs, it follows that they could be employed as standard candles, similar to type Ia supernovae.
It follows from this that if the Amati parameters (, ) are known, then the LGRB distance moduli  can be calculated regardless of the cosmological model using the Amati relation.Indeed, substituting the expression for the isotropic equivalent radiated energy  iso with the corrected bolometric fluence, given by Equation ( 11), into the Amati relation, one can express from there the energy distance   (Equation ( 7)) and convert it into the distance modulus  (Equation ( 8)).This will result in the following expression: We will call this expression the extended Amati relation.
In accordance with equations ( 5), (7 for   ), and (8), the distance modulus is given by where  represents the corresponding cosmological model derived function.So, the distance scale depend on the value of the Hubble constant  0 .The Pantheon SNIa catalogue uses  0 = 70 km/s/Mpc to determine the absolute magnitudes of supernovae .In the subsequent part of the article, we verify our methods by varying the parameters of the cosmological model, including  0 (the correctness of the method implies that the obtained values of  0 should be equal to the initial one within the margin of error).Therefore, the distance scale and all the values of the Amati parameters obtained in this paper are tied  0 = 70 km/s/Mpc.Changing this value will, in particular, result in a change in the values of  and .

Direct problem
We will divide the problem-solving process into separate stages, providing a detailed description for each of them.
In order to solve the direct problem, one need to determine the parameters of the ΛCDM model using the provided LGRB sample by following stages.
(i) Obtaining initial data.At the initial stage, we have the Swift LGRB catalogue, described in detail in Section 3.2.It is a table of the following form as in Table 1.Each LGRB in this table is described by a set of parameters (, ,  p ,  obs ).The parameters  and  p have asymmetric errors.The symbols Δ l and Δ u mean the "bottom" and "top" errors, respectively.The symbol Δ 90% means the error corresponding to the 90% confidence interval for a normal distribution.To bring this error to a standard deviation (confidence interval 68.3%), one need to divide it by a factor 1.6 ≈ Φ −1 (1 + 0.9)/2 , where Φ is the integral function of the standard normal distribution and Φ −1 is its inverse function (quantile).As mentioned in Section 3.2, the LGRB catalogue is incomplete, because some errors are missing.In total, there are 174 LGRBs in our sample.
(ii) Extending the catalogue.We decided to fill the missing errors in the original table by calculating the median for the relative errors of the corresponding parameter.After these manipulations, we got the following sample in Table 2.
(iii) Calculation of the quantities  iso and  p, .Using some given value for the parameter  (For example,  = 0, 0.5, 1), we can use the formulas from Section 4.1 to calculate the isotropic equivalent radiated energy  iso , given by Equation ( 12), and the peak energy in the source frame  p, , given by Equation ( 13), for each LGRB.
(iv) Computing the Amati parameters.Using some of the linear regression methods, we can obtain an estimate for the parameters (, ), which are the usual parameters of a linear regression in the Amati relation ( 14).
(v) Calculation of gamma-ray burst distance modules.Given the parameter  and the computed Amati parameters (, ), we can find the LGRB distance modulus  A GRB using the extended Amati relation ( 16).
(vi) Constructing a HD.With the  A GRB values we calculated above and with the -source catalogue, we can construct a redshiftdistance modulus diagram, also called a HD.
(vii) Estimating the parameters of the ΛCDM model.Using approximation algorithms, one can approximate the theoretical dependence () (Equation ( 8)) depending on the parameters of the cosmological ΛCDM model by using a some initial sample, for instance, the SN Ia supernovae catalogue to the HD constructed for our LGRBs.
But this problem has an issue when passing from point 2 to point 3.The fact is that Equation ( 12) for  iso contains the distance   , which, in turn, depends on the choice cosmological model, what can be seen from Equations (5,6,7).Because of this, it turns out that the Amati parameters (, ), and, as a consequence, the distance moduli  A GRB of LGRBs are also dependent on the cosmological model!As a result, the parameters of the ΛCDM model ( 0 , Ω m , Ω Λ , Ω k ) that we get at the output also depend on the chosen dependency () and, as a consequence, on some initially chosen cosmological parameters ( 0 , Ω m , Ω Λ , Ω k ) 0 .This is known as a circularity problem (Kodama et al. 2008).

Addressing the circularity problem
There are several approaches to resolve the circularity problem.The first approach based on using known standard candles such as supernovae, globular clusters or others to calibrate near GRBs as a standard candles.In particular, in Lovyagin et al. (2022) we used approximation of the dependence () by the HD via type Ia supernovae from the Pantheon catalogue.The second approach is to simultaneously constrain the calibration parameters and the cosmological parameters by considering a chosen likelihood function, analogically Amati & Della Valle (2013); Demianski et al. (2021).
We also suggest here another approach based on method of iterations.Thus, the parameters ( 0 , Ω m , Ω Λ , Ω k )  obtained at the output of the algorithm can be used in the repeated solution when passing between point 2 and point 3, and as a result, new estimates for these parameters can be obtained ( 0 , Ω m , Ω Λ , Ω k ) +1 , and so on until convergence is reached.
While we assume that this approach remains robust and accurate within the framework of Bayesian statistics, we refrain from employing it in this article, leaving it as a potential future research.Instead, the article addresses the inverse cosmological calibration problem as described in Section 5.3, i.e., obtaining the Amati parameter values within the framework of the ΛCDM model with the gravitational lensing and Malmquist bias correction (Shirokov et al. 2020c), but not the parameters of the cosmological model itself.When addressing the direct problem, our focus is limited to constructing the HD (item vi in 5.1) for illustrative purposes.The scheme of the direct problem solution approach for is illustrated in Figure 1: the completed part of the work is highlighted in black, while the possible addressing the circularity problem within this context is shown in red.

Inverse cosmological calibration problem
The core of the idea of the inverse cosmological calibration problem (ICCP) is that we use the least-squares algorithms to minimise the residuals from the HD ( model −  A GRB ), while the points   are calculated from the log  p, -log  iso plot (the Amati plane) that is varied by (, , ).In this study, we fix parameters ( 0 , Ω m , Ω Λ , Ω k ) of the ΛCDM model, so optimal values for the parameters , , and  are determined using this method.
Our main objective is to solve the ICCP, which involves determining the parameters of the Amati parameters (, ) and correction parameter  using the provided LGRB sample by following stages.
(i) Using obtained initial data and extended catalogue as described in 5.1.
(ii) Obtaining ΛCDM parameters.This can be done, for example, using the Ia supernovae of the Pantheon catalogue as an input calibration sample (since the ICCP can be considered as another way to calibrate LGRBs from supernovae, it can address the circularity problem).
(iii) Determining parameters using least squares to minimise residuals of the LGRB HD.In current research we fix cosmological model parameter determined at the second stage.In the first case, we also fix the -parameter to be equal to 0 and determine the values of  and .In the second case we determine all three values ,  and .
A detailed scheme for solving the inverse problem is shown in Figure 2.
In this approach, the parameter  is considered as a correction parameter to the ΛCDM model.It is used as a correction for gravitational lensing and the Malmquist bias (GLMB-correction from (Shirokov et al. 2020c)).Determining it by solving the direct problem without addressing the circularity problem is not possible, whereas the inverse problem is free from this flaw.Assuming the accurate determination of the fundamental parameters of the ΛCDM model, fixing  = 0 should yield the Amati parameter values  and  within the uncertainty limits that match those obtained from solving the direct problem without addressing the circularity problem.

Approximation and regression algorithms used
As evident from Section 5 various stages required the usage of both linear regression (such as in calculating the Amati parameters) and non-linear regression (for example, in estimating the parameters of the ΛCDM model using the HD).
As an algorithm for linear regression, we have chosen the Theil-Sen method (Gilbert 1987).The essence of this method is that for a given set of points {(  ,   )}  the estimated slope  of the linear function  =  +  is determined as the median value of all possible slopes (  −   )/(  −   ) for each pair of points.The coefficient  of the linear relationship can then be estimated as the median of the values   −   .The main advantage of this method over the classical method of least squares is its robustness, i.e. resistance to outliers.
In the context of non-linear regression, we employed the least_squares procedure, which is implemented in the scipy.optimizepackage of the Python programming language.This routine employs a technique known as the "trust region reflective" method to identify the minimum of the variance function as follows: where  represents the vector of estimated parameters,   denotes the model's deviation from the data, and  denotes the loss function.
Choice for the identity mapping as the loss function  simplifies the problem into the classical finding of the minimum sum of squared deviations.However, this approach can introduce significant bias in parameter estimates due to data outliers.To mitigate the impact of outliers, an alternative loss function must be selected.For instance: where  serves as a conditional soft boundary, separating outliers.It serves as an illustration within this study, however, this approach suffers from a dead-loop (circularity problem) and cannot be used as a cosmological model parameter test.A potential solution using iterative methods is highlighted in red.In this study, it is stated as a task for a possible future research; instead, a solution to the inverse problem, illustrated in Figure 2, has been employed.
For the value of this constant, we have chosen doubled the standard deviation of   .

Monte Carlo method for calculating errors of indirect measurements used
We chose the Monte Carlo method as the error propagation method (Anderson 1976;Albert 2020;Gorokhov et al. 2023).To proceed to the description of the essence of this method, first one need to decide on the interpretation of the error.
The standard approach for addressing the error propagation issue is referred to as the linear uncertainty propagation (LUP) theory.In this approach, errors are treated as standard deviations, and the values are

Enhanced catalogue
The log E iso -log E p,i plane

GRB Hubble Diagram
Figure 2. Scheme for solving the inverse problem.The influence of data and parameters on the calculated values is illustrated here.From the Pantheon type Ia supernova HD (top right), we determine the ( 0 , Ω m , Ω Λ , Ω k ) parameters of the ΛCDM model, which we then assume to be fixed.Parameters (, ,  ) are required.To estimate them, we minimise the deviations in the HD for LGRBs (right).The fixed parameters of the ΛCDM model define the dependence of the distance on the redshift  (), which is shown by the blue curve in the HDs.To find  iso , one need to use the formulas from Sections 2 and 4.1.These formulas depend on the choice of  (Equation ( 11)) and on the dependency of  () (Equations (12, 7)).Therefore, the value of  iso also depends on  and fixed parameters of the ΛCDM model ( 0 , Ω m , Ω Λ , Ω k ).The Amati parameters ,  define the position of the straight line in the plane log  iso -log  p, , which we called the log  iso -log  p, plane.In this case, all three sought-for parameters , ,  affect the LGRB distance moduli calculated using the extended Amati relation 16.
regarded as the means.Thus, variable values along with their errors are determined using the normal distribution, which represents the potential range of the variable through distribution parameters.When a measured value  is provided with an error of , it signifies that this value is random and follows a normal distribution with a mean of  and a standard deviation of .
The standard formula for calculating the error of indirect measurement is derived from simple considerations.Let us have  normally distributed random variables with known expected values {  }  =1 , variances { 2  }  =1 and covariances {   } ≠  .Let  also be a linear function of these quantities: The value of  will also have a normal distribution.Its variance can be obtained as follows: or, if we treat the standard deviation as an error, If the function  is non-linear, it is represented as an expansion in a Taylor series up to the first order: This is the standard formula for calculating the error of indirect measurements.Unfortunately, this approach has limited applicability due to the following disadvantages: • Measured values and errors are treated as expected values and standard deviations of random variables that follow a normal distribution.Therefore, this method is not applicable to data that has asymmetric errors.
• It is necessary that the errors be small enough and the functions smooth enough that the inaccuracy due to representing the function in a linearised form is negligible.
In our case, the observational data exhibit asymmetric errors with large values.Also, the data transformation function  is far from linear.In addition, it cannot be guaranteed that the indirect (i.e., transformed) values will obey the normal distribution.This can be demonstrated with a simple example.
Let two quantities with errors be given:  ±  and  ±  .We want to find their ratio / and its error  / .As is known, the value, which is the ratio of two normally distributed random variables, obeys the Cauchy distribution (among astrophysicists it is better known as the Lorentz profile).The latter is a standard example of a distribution for which neither the expected value, nor the variance, nor the higherorder moments are defined.This example shows that for our case, not only the standard error propagation formula ( 23) is inapplicable, but also the interpretation of values with an error as random variables following the normal distribution with specified mean and variance is also not valid.Therefore, it is necessary to choose another interpretation of the errors and error propagation method, which would be a generalisation to arbitrary distributions.
The concepts of mean and standard deviation can be naturally generalised using quantiles.The 0.5 level quantile (median) is defined for any continuous distribution, and in the case of a normal distribution, it coincides with the expected value (mean).The quantiles 0.16 and 0.84 are the bounds of the confidence interval with the confidence level 0.68, which in the case of a normal distribution reduces to the interval with bounds ( − ,  + ).Thus, for arbitrary distributions, we will use the median as a "generalised mean" instead of the mean, and the differences of the 0.84 level quantile with the median ( 0.84 −  0.5 ) and the median with the 0.16 level quantile ( 0.5 −  0.16 ) as the top and bottom errors, respectively.
We can now move on to describing how to calculate indirect measurement (i.e., transformed values) errors using the Monte Carlo method.The concept is as follows: given that each input value is represented by a measured value and its error value, which define a form of continuous distribution of a random variable (typically a normal distribution in the simplest scenario), we can simulate this distribution using a random number generator.This simulation involves generating a sample of size  that follows the appropriate distribution.That is, each of the input values is represented as a sample of size .The problem can be represented as  different realisations of our experiment, in each of which all indirect measurements are calculated.Thus, for each indirect value, we will also have a sample of size , which will be a model of the distribution of this value as a random variable.Based on this sample, it is possible to determine the quantiles of the 0.16, 0.5 and 0.84 levels in order to estimate the value of this value and its errors.
If the initial value has symmetrical errors, then the normal distribution can be used as a model one for it.For values that have different upper and lower errors, we used a distribution consisting of two halves of different Gaussian distributions.This is also the so-called split-normal distribution.An example of such a distribution, as well as an example of using the Monte Carlo method to estimate the errors of indirect measurements, is shown in Figure 3. Another important point is that the value  p (the peak energy of the LGRB spectrum) cannot be negative.Some of the  p values from our LGRB data set had remarkably huge lower uncertainties, so that trying to create a Monte-Carlo sample could lead to negative  p values, which has absolutely no physical sense.Because of that, in the case of  p variables, we draw the split-normal distributions in the space of log  p .To prevent "tails" of distributions from falling into the negative area, we created Monte Carlo model samples not for the  p value itself, but for its logarithm.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.
This method has a number of advantages.Firstly, it is easy to implement and interpret the results.Secondly, this method is universal in the sense that it is applicable to any kind of function.Thirdly, in this method, any correlations between random variables are automatically taken into account in subsequent calculations, without the need to calculate derivatives.The only significant disadvantage is the increase in calculation time by  times (in our case  = 10 000), but we live in an era where we can afford it.

SOFTWARE IMPLEMENTATION
To solve the tasks, we wrote Python code.The corresponding repository is publicly available at https://github.com/Roustique/sngrb.The repository includes our software implementation of the following steps: • Estimating the parameters ( 0 , Ω m , Ω Λ , Ω k ) of the ΛCDM model from the HD (  ,   )  for a sample of observational objects.
• Solving the direct problem (see Section 5, Figure 1) of constructing a HD for LGRBs.
• Solving the inverse problem (see Section 5, Figure 2) of finding the Amati parameters , , and the correction parameter  based on a fixed cosmological model.
The software implementation was done utilising the following libraries.
• NumPy -for supporting arrays, including multidimensional ones, and functions that operate on them.lg(E iso /1keV) Figure 3.An example of the application of the Monte Carlo method for one of the LGRBs (GRB110422A).Upper block: original 4 parameters ,  p ,  obs , and .For the first three, Monte Carlo distributions of size  = 10 000 are presented.The vertical solid black line marks the median, and the dotted line marks the errors.The redshift  has no error and is therefore shown in the figure as a single vertical line.Bottom block: distributions of the logarithms of the parameters  p, and  iso calculated using the Monte Carlo method.The solid lines also mark the positions of the medians, and the dotted lines mark the errors.
• SciPy -we used the optimize.least_squaresprocedure for non-linear regression and stats.mstats.theilslopesfor linear regression using the Theil-Sen method (details in Section 6.1).We also used the stats module for implementing cumulative and differential distribution functions, as well as quantiles.
• Pandas -for working with catalogues.

Obtaining the ΛCDM model parameters from supernovae
Using the procedure least_squares with a loss function that is the identical mapping (written in detail in Section 6.1), we have estimated the parameters of the ΛCDM model from supernovae from the Pantheon catalogue.We have used 4 different models: in one, all 4 cosmological parameters ( 0 , Ω m , Ω Λ , Ω k ) varied, and in the others, the parameters  0 = 70 and/or Ω k = 0 were fixed.The standard model is shown in Figure 4 with a black line and its parameters, along with the parameters of other models, are presented in Table 3.The four cases considered allow us to conclude that our approximation algorithms work correctly, since the values of the cosmological parameters remain close to ( 0 = 70, Ω m = 0.3, Ω Λ = 0.7).
Since the Pantheon SN Ia catalogue data is known to be bound to  0 = 70 and the ΛCDM model is bound to Ω k = 0 we decided to choose a model with fixed values of this parameters as the standard.Other models were used to check the robustness of our methods.In this case we have obtained Ω m = 0.296 and Ω Λ = 0.704 as best-fitting values4 .

Direct problem
To solve the direct problem, we used the previously obtained ΛCDM parameters Ω m = 0.296, Ω Λ = 0.704, and  = 0.In this context, the Amati parameters  and  were determined through linear approximation employing the Theil-Sen method.More precisely, by using our Monte-Carlo sample with a size of  = 10, 000 (as detailed in Section 6.2), we acquired 10, 000 sets of  and  values using the Theil-Sen method.From these sets, we calculated the medians and quantiles for the Amati parameters:  = 0.92 +0.10 −0.11 and  = 50.48+0.28  −0.27 .With these values, we constructed the joint SN Ia and LGRB HD.The results are shown in Figure 5.The algorithm provides optimal values for the Amati parameters and extends the HD at the LGRB scales.

Inverse cosmological calibration problem
The results of our algorithm for solving the ICCP case (as detailed in Section 5.3) are shown in Figure 6.We have regarded two models: one with a fixed  = 0 and the other with a free .The input cosmological model is based on the ΛCDM parameters obtained in Section 8.1.In the ICCP case with  = 0, the results coincide within the margin of errors with the results obtained for the direct problem with the same value of  = 0. Thus, considering the HD to minimize residuals has minimal impact on the Amati parameter values.However, when we introduce variation in the parameter , its influence becomes evident, with an optimal value of  = 2 and a reduction in the correlation in the log  iso -log  p, plane.Currently, the Monte-Carlo samples has a more regular configuration, and the joint SN Ia and GRB HD appears visually more regular as well.

CONCLUSIONS
The sample of 174 LGRBs from the Swift mission was calibrated using the non-parametric statistical methods (partially described in Lovyagin et al. (2022) andGorokhov et al. (2023)), the ΛCDM model as potentially changeable basis, and the extended Amati relation with cosmological correction for gravitational lensing and Malmquist bias (GLMB) from Shirokov et al. (2020c).The calibration of the Amati parameters and the GLMB correction (, , ) was carried out in such a way that residuals in the  - diagram (the HD) were minimised, while its points   were calculated from the log  p, -log  iso plot (the Amati plane) that was varied by (, , ).This way we called as inverse cosmological calibration problem (ICCP).16)) from the ΛCDM model  () (Equation ( 8)) on the HD).Middle: the log  iso -log  p, plane.Bottom: the joint SN Ia and LGRB HD.
Thus, we suggest the ICCP as a new method resolving the circularity problem of the LGRB calibration as standard candles.In this research we used the ΛCDM model parameters found from the type Ia supernovae HD from the Pantheon catalogue, except for the Hubble constant which was fixed at  0 = 70 km/s/Mpc.Two cases were tested, one with a fixed GLMB parameter  = 0 (the case of no GLMB correction) and the second with variation of all three parameters , ,  that is of interest to us.Additionally, we have compared the case of ICCP solution with  = 0 to the solution of the direct cosmological problem (we called an estimation of  and  using linear regression in the log  iso -log  p, plane with a fixed cosmological model as a direct problem).Although the illustrative solution of the direct problem suffers from the circularity problem, it can be seen that the results are very similar, the parameters coincide up to an error.Thus, we can conclude that our methods and determination of the parameters of the ΛCDM model were accurate (including the fact that parameters of ΛCDM model fit the HD for LGRB under the conditions of this problem within the current limits of accuracy).
Adding a third degree of freedom in the form of a correction parameter  affecting the bolometric fluence, given by Equation (11), and hence the extended Amati relation ( 16), causes a significant upward bias on the factor ( + 1) 2 in estimates of the LGRB distance moduli   in the HD.While the Amati parameters  = 0.92 +0.12 −0.12 ,  = 50.32+0.33  −0.32 for  = 0 are changing in the direction of decreasing the correlation with  = 0.63 +0.13  −0.14 ,  = 50.12+0.33 −0.31 for  = 1.98 +0.25 −0.24 .All estimates of the Amati parameters are presented in Table 4.It is natural to assume that the Amati correlation of the bolometric fluence ( iso ) with its hardness ( p, ) depends on the mass of the collapsing star.Therefore, for physical reasons, the correlation should be related to the evolution of stars, and also serve as a test of alternative ideas, such as the existence of quark stars (Sokolov 2015(Sokolov , 2016(Sokolov , 2019) ) and test of gravitational theories (Baryshev 2020;Shirokov et al. 2020b).
The accuracy of the observational values of LGRBs, as well as the number of their observations at the moment, leave much to be desired.However, we conclude that LGRBs can be used as standard candles, allowing one to expand the HD up to  ≈ 10, and their position on the HD can be correctly consistent with the selected cosmological model.Apparently, the standard ΛCDM model is indeed consistent with the indirect estimates of the LGRB HD, however, this requires further research using our method and varied models as input basis, as well as increasing the sample size of LGRBs.As an application of our methods for research we suggest addressing the circularity problem analogically Amati & Della Valle (2013); Demianski et al. (2021) via iterations and solving the inverse problem to estimate parameters of cosmological models basing on the residuals values of observational data in the Bayesian statistics paradigm.Thus, the ICCP idea can be used as an alternative cosmological model test in nearest future. LGRB

Figure 1 .
Figure1.Visualisation of the direct problem and its circularity.The scheme represents the process of deriving HD from the SN-derived ΛCDM model parameters (black).It serves as an illustration within this study, however, this approach suffers from a dead-loop (circularity problem) and cannot be used as a cosmological model parameter test.A potential solution using iterative methods is highlighted in red.In this study, it is stated as a task for a possible future research; instead, a solution to the inverse problem, illustrated in Figure2, has been employed.

Table 3 .
The estimated parameters of the ΛCDM model based on supernovae from the Pantheon catalogue.The model with fixed  0 = 70 and Ω k = 0 is chosen as the standard.

Figure 4 .
Figure 4. Curves corresponding to ΛCDM models with parameters found from supernova data.The blue line is the model with all free parameters.The orange line is the model with the fixed Hubble constant  0 = 70.The green line is the model with a fixed curvature contribution Ω k = 0.The black line is the model with fixed parameters  0 = 70, Ω k = 0.

Figure 5 .
Figure 5.The results of the direct problem algorithm with fixed  = 0. Top: the Monte Carlo distributions for the parameters  and , obtained as a result of linear approximation by the Theil-Sen method in the log  iso -log  p, plane.The estimates of these parameters with their errors are shown.Middle: the log  iso -log  p, plane.Bottom: the joint SN Ia and LGRB HD.

Figure 6 .
Figure 6.The results of the ICCP algorithm.Left: the case of fixed  = 0. Right: all three parameters , ,  are varied.Top: The Monte Carlo distributions for the varied parameters, obtained as a result of solving the ICCP (minimising the deviations of  A GRB (Equation (16)) from the ΛCDM model  () (Equation (8)) on the HD).Middle: the log  iso -log  p, plane.Bottom: the joint SN Ia and LGRB HD.

Table 1 .
Form of the original catalogue of LGRBs.The full catalogue is presented in Appendix B.

Table 2 .
Form of the extended catalogue of LGRBs. obs errors are reduced to 1.The full catalogue is presented in Appendix C.