# Model-Independent Determination of $H_{0}$ and $\Omega_{K,0}$ using Time-
Delay Galaxy Lenses and Gamma-Ray Bursts

Shen-Shi Du,1,5 Jun-Jie Wei,2,3 Zhi-Qiang You,4,5 Zu-Cheng Chen,4,5 Zong-Hong
Zhu,1 En-Wei Liang6
1School of Physics and Technology, Wuhan University, Wuhan, Hubei 430072,
China
2Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210023,
China
3School of Astronomy and Space Sciences, University of Science and Technology
of China, Hefei 230026, China
4Department of Astronomy, Beijing Normal University, Beijing 100875, China
5Advanced Institute of Natural Sciences, Beijing Normal University, Zhuhai
519087, China
6Guangxi Key Laboratory for Relativistic Astrophysics, School of Physical
Science and Technology, Guangxi University, Nanning 530004, China
duss.physics<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Combining the ‘time-delay distance’ ($D_{\Delta t}$) measurements from galaxy
lenses and other distance indicators provides model-independent determinations
of the Hubble constant ($H_{0}$) and spatial curvature ($\Omega_{K,0}$), only
based on the validity of the Friedmann-Lemaître-Robertson-Walker (FLRW) metric
and geometrical optics. To take the full merit of combining $D_{\Delta t}$
measurements in constraining $H_{0}$, we use gamma-ray burst (GRB) distances
to extend the redshift coverage of lensing systems much higher than that of
Type Ia Supernovae (SNe Ia) and even higher than quasars, whilst the general
cosmography with a curvature component is implemented for the GRB distance
parametrizations. Combining Lensing+GRB yields $H_{0}=71.5^{+4.4}_{-3.0}$ km
s-1Mpc-1 and $\Omega_{K,0}=-0.07^{+0.13}_{-0.06}$ (1$\sigma$). A flat-universe
prior gives slightly an improved $H_{0}=70.9^{+4.2}_{-2.9}$ km s-1Mpc-1. When
combining Lensing+GRB+SN Ia, the error bar $\Delta H_{0}$ falls by 25%,
whereas $\Omega_{K,0}$ is not improved due to the degeneracy between SN Ia
absolute magnitude, $M_{B}$, and $H_{0}$ along with the mismatch between the
SN Ia and GRB Hubble diagrams at $z\gtrsim 1.4$. Future increment of GRB
observations can help to moderately eliminate the $M_{B}-H_{0}$ degeneracy in
SN Ia distances and ameliorate the restrictions on cosmographic parameters
along with $\Omega_{K,0}$ when combining Lensing+SN Ia+GRB. We conclude that
there is no evidence of significant deviation from a (an) flat (accelerating)
universe and $H_{0}$ is currently determined at 3% precision. The measurements
show great potential to arbitrate the $H_{0}$ tension between the local
distance ladder and cosmic microwave background measurements and provide a
relevant consistency test of the FLRW metric.

###### keywords:

gravitational lensing: strong – gamma-ray burst: general – cosmological
parameters

††pubyear: 2022††pagerange: Model-Independent Determination of $H_{0}$ and
$\Omega_{K,0}$ using Time-Delay Galaxy Lenses and Gamma-Ray Bursts–A

##  1 Introduction

The precise value of Hubble constant ($H_{0}$) has emerged as one of the most
vital cracks in modern cosmology, with beyond $5\sigma$ tension between the
local (late-universe) determination from Cepheid variable-calibrated Type Ia
supernovae (SNe Ia; Riess et al. 2019; Riess et al. 2021) and the high-
redshift (early-universe) measurement by _Planck_ cosmic microwave background
(CMB; Planck Collaboration et al. 2020) observations based on the $\Lambda$
Cold Dark Matter ($\Lambda$CDM) cosmological model. This issue has spurred
intense debate of beyond-$\Lambda$CDM dark energy component or early-universe
new physics as well as unaccounted-for systematic effects in either or both of
the measurements (e.g., Riess et al. 2011; Marra et al. 2013; Verde et al.
2013; Riess et al. 2016; Bernal et al. 2016; Riess 2019; Verde et al. 2019;
Freedman et al. 2019; D’Amico et al. 2021; Di Valentino et al. 2021).
Meanwhile, various recent measurements of present-day spatial curvature
($\Omega_{K,0}$) have been shown to mutually disagree with each other, leading
to another potential cosmic crisis (Di Valentino et al., 2020). New,
independent measurements of the $H_{0}$ and $\Omega_{K,0}$ are pursued to
provide relevant clues to such discrepancies (see Di Valentino et al. 2021 for
a recent review).

Time delay cosmography has become a mature field with the combined strengths
of modern monitoring and lens modeling of quasar strong gravitational lensing,
which provides a novel measurement of $H_{0}$ using the time delays between
multiple gravitationally- and geometrically-induced images (Refsdal 1964;
Blandford & Narayan 1992; Treu & Marshall 2016; Bonvin et al. 2017; Suyu et
al. 2018; Oguri 2019; Wong et al. 2020; Birrer et al. 2020; Denzel et al.
2021). The statistical combination of six measurements on different quasar
time-delay lenses had reached a 2.4% measurement of $H_{0}=73.3^{+1.7}_{-1.8}$
km s-1Mpc-1 in the flat $\Lambda$CDM model (Wong et al., 2020); Considering
the mass-sheet transformations, which leave the lensing observables unchanged,
as the dominant source of residual uncertainty in $H_{0}$ estimation, Birrer
et al. (2020) obtained $74.5^{+5.6}_{-6.1}$ km s-1Mpc-1 inferring from the
TDCOSMO sample with seven lenses. The measurements are not only independent of
the early- and late-universe estimates, but are also proposed to control the
systematic errors – around 40 strong lensing systems can measure $H_{0}$ at a
sub-percent level (Jee et al., 2016; Shajib et al., 2019), thus raising the
prospect of adjudicating the tension between CMB and the local distance
ladder. Nevertheless, time-delay cosmography relies on the assumption of a
specific background cosmology, and the inferred value of $H_{0}$ could be
significantly different with using different cosmological models (see Table 6
in Wong et al. 2020). Thus, this measurement is strongly model-dependent.

A more intriguing application of time-delay galaxy lenses (TDGLs) is combining
the TDGLs and such reliable distance indicators from other probes to constrain
both the $H_{0}$ and $\Omega_{K,0}$ via the distance sum rule (DSR; Collett et
al. 2019), merely assuming that the light propagates along null geodesics in
the Friedmann-Lemaître-Robertson-Walker (FLRW) metric. This method only relies
on the validity of the classical FLRW metric and geometrical optics, and
infers $H_{0}$ and $\Omega_{K,0}$ with TDGLs providing the “time-delay
distance" ($D_{\Delta t}$) and with the distance indicators extending the
redshifts to the lens and source to cover that of the lensing systems well.
Additionally, the measurement of $\Omega_{K,0}$ from two or more lens-source
pairs provides a consistency test of the FLRW metric (Räsänen et al., 2015).
Such measurements has been implemented in combinations with, for instance SNe
Ia (Collett et al., 2019), multi-messenger gravitational-wave “standard
sirens" (Liao, 2019), and ultra-compact structure in radio sources (Qi et al.,
2020). We also notice that Wei & Melia (2020) proposed using quasars (QSOs) as
distance indicators up to $z\sim 5$ to match the distances of lensing systems
that can be detectable by the upcoming Rubin Observatory Legacy Survey of
Space and Time (LSST; Oguri & Marshall 2010; Ivezić et al. 2019), which is
expected to take full advantage of the lensing catalog in hand. While recent
works were indicative of the $\sim 4\sigma$ tension between the quasars and
flat $\Lambda$CDM at high redshifts (e.g., $z\geq 1.5$) and showed the
vigilance of using quasars as standardizable candles to constrain cosmological
parameters (e.g., Risaliti & Lusso 2019; Lusso et al. 2020; Yang et al. 2020;
Velten & Gomes 2020; Khadka & Ratra 2021). This prompts us to further explore
robust high-redshift distance indicators to consolidate the independent
measurements of $H_{0}$ and $\Omega_{K,0}$ in conjunction with TDGLs via the
DSR.

Gamma-ray bursts (GRBs) are the most luminous astrophysical phenomenon, which
are expected to be emerged as early as death of the first-generation stars
($z\sim 20$). The central engine and radiation mechanisms as well as the jet
composition of GRBs, which are responsible for the energetic $\gamma$-ray
emission, remain inconclusive. Some aspects of the progenitor models have been
well investigated in light of the fact that several long-duration GRBs (LGRBs)
are observed in association with core-collapse supernova explosions (see e.g.,
Zhang 2018). Despite these difficulties, the high luminosity of GRBs and the
immunity of $\gamma-$rays from dust extinctions contribute to their high
detectabilities at distances that are far beyond that SNe Ia (or even quasars)
are available, thus GRBs are promising objects serving as cosmological probes
to provide meaningful cosmological constraints over a wide redshift range (Dai
et al., 2004; Liang & Zhang, 2005; Ghirlanda et al., 2006; Amati et al., 2008;
Amati et al., 2019; Wang et al., 2015; Demianski et al., 2017, 2021; Khadka et
al., 2021). One of the most important observational properties of LGRBs is the
phenomenological relations between the observed spectral and intensity
quantities (radiated $\gamma-$ray energy or luminosity) during the prompt
phase, among which the Amati correlation (Amati et al., 2002) between the $\nu
F_{\nu}$ spectral peak energy, $E_{\rm p}$, and the bolometric isotropic-
equivalent radiated energy, $E_{\rm iso}$, has been extensively explored as
the distance indicators of LGRBs (e.g., Amati et al. 2008; Amati et al. 2019;
Demianski et al. 2017, 2021). The main obstacles that preclude the
cosmological use of GRBs come from the potential selection effects, the
“circularity problem" in calibrating the GRB luminosity correlations, and
their possible redshift evolutionary effects (see a recent review by Wang et
al. 2015). Several works have investigated the observational selection effects
and show a general conclusion that such effects, even though still debated,
might be minor (e.g., Amati 2006; Butler et al. 2009; Nava et al. 2012;
Demianski et al. 2017). The circularity has been effectively overcome with
several reliable calibration procedures (Wang et al., 2015; Amati et al.,
2019). We will elaborate on the impact of possible redshift evolutionary
effects on the Amati correlation when modelling the GRB distance indicators.

In this work, we propose to combine the $D_{\Delta t}$ measurements from TDGLs
and distances determined with the GRB Hubble diagram to restrict $H_{0}$ and
$\Omega_{K,0}$ via the DSR. The general cosmography with an arbitrary spatial
curvature component is implemented for the GRB distance parametrizations, with
which the bounds of $H_{0}$ and $\Omega_{K,0}$ are contributed by both the GRB
distances and $D_{\Delta t}$. To alleviate the circularity, we simultaneously
fit the GRB distance indicators and the parameters in theoretical distances in
an effective Bayesian framework, and $H_{0}$ and $\Omega_{K,0}$ can be
obtained by marginalizing over the remaining parameters. The rest of this
paper is arranged as follows. In Section 2, we depict our method to determine
the $H_{0}$ and $\Omega_{K,0}$ with the combined observations of TDGLs and the
popular high-redshift probes, including the most recently released SNe Ia,
GRBs, and QSOs, which we denote as standardizable-candle (SC) objects
hereafter. The observational data and likelihood constructions are detailed in
Section 3. We present the results and discussions in Section 4. Summary is
made in Section 5.

##  2 Methodology

###  2.1 Constraining Hubble constant and spatial curvature through distance
sum rule

In what follows, we present the method for simultaneously estimating the
$H_{0}$ and $\Omega_{K,0}$ with only the homogeneity and isotropy ansatz in
cosmology, stating that the space-time is described by the FLRW metric,

${\rm d}s^{2}=-c^{2}t^{2}+a^{2}(t)\left[\frac{{\rm d}r^{2}}{1-\kappa
r^{2}}+r^{2}d\Omega^{2}\right],$ (1)

where $c$ is the speed of light, $a(t)=(1+z)^{-1}$ is the scale factor, and
the dimensionless curvature $\kappa$=1, 0, and $-1$ corresponds to a spatially
closed, flat, and open universe, respectively. The cosmic expansion rate is
defined as $H(z)\equiv{\dot{a}}(t)/a(t)$ with its present value is $H_{0}$
(hereafter, the subscript “0" represents the present value), which acts as a
normalization of the Hubble diagram. Let $D_{\rm A}(z_{s},z_{l})$ be the
angular diameter distance of a source emerging at redshift $z_{s}$ as lensed
at $z_{l}$, where subscripts “$s$" and “$l$" correspond to source and lens,
respectively. Assuming that the geometrical optics holds, the dimensionless
co-moving distance, $d(z_{s},z_{l})\equiv(1+z_{s})H_{0}D_{\rm
A}(z_{s},z_{l})/c$, is given by

$\begin{split}d(z_{s},z_{l})=\frac{1}{\sqrt{|\Omega_{K,0}|}}{\rm
sinx}\left({\sqrt{|\Omega_{K,0}|}}\int^{z_{s}}_{z_{l}}\frac{d\tilde{z}}{E(\tilde{z})}\right),\end{split}$
(2)

where $E(z)\equiv H(z)/H_{0}$, the curvature density parameter
$\Omega_{K,0}\equiv-\kappa c^{2}/a_{0}^{2}H_{0}^{2}$, and sinx is sinh (sin)
in the case of $\Omega_{K,0}>0$ ($\Omega_{K,0}<0$) and, for $\Omega_{K,0}=0$
this equation simplifies as the integral. With the denotation of $d(z)\equiv
d(0,z)$, Equation (2) links $d_{ls}\equiv d(z_{s},z_{l})$, $d_{s}\equiv
d(z_{s})$, and $d_{l}\equiv d(z_{l})$ in the FLRW metric via the distance sum
rule (Räsänen et al., 2015),

$d_{ls}=\epsilon_{1}d_{s}\sqrt{1+\Omega_{K,0}d^{2}_{l}}-\epsilon_{2}d_{l}\sqrt{1+\Omega_{K,0}d^{2}_{s}},$
(3)

where $\epsilon_{i}=\pm 1$ with the positive case is required for the one-to-
one correspondence between $z$ and $t$, with $d^{\prime}(0)>0$. We assume the
case of $\epsilon_{i}=1$, then Equation (3) can be rewritten as

$\frac{d_{ls}}{d_{s}d_{l}}=\sqrt{1/d_{l}^{2}+\Omega_{K,0}}-\sqrt{1/d_{s}^{2}+\Omega_{K,0}}.$
(4)

Strong lensing systems benefit from their very simple underlying physics and
purely gravitational origin compared to the other cosmological probes. The
measurement of time delay $\Delta t_{i,j}$ between lensed images $i$ and $j$
at the coordinates $\bm{\theta}_{i}$ and $\bm{\theta}_{j}$ and a concrete lens
model provide the measurement of time-delay distance $D_{\Delta t}$ through
(Refsdal, 1964; Schneider et al., 1992)

$\Delta t_{i,j}=\frac{D_{\Delta t}}{c}\Delta\psi_{i,j},$ (5)

where
$\Delta\psi_{i,j}=[(\bm{\theta}_{i}-\bm{\beta}_{s})^{2}/2-\psi(\bm{\theta}_{i})-(\bm{\theta}_{j}-\bm{\beta}_{s})^{2}/2+\psi(\bm{\theta}_{j})]$
represents the “Fermat potential" difference between lensed images with the
two-dimensional lensing potential $\psi$ and source position $\bm{\beta_{s}}$.
The time-delay distance is the combination of the three angular diameter
distances:

$D_{\Delta
t}=(1+z_{l})\frac{D_{A,l}D_{A,s}}{D_{A,ls}}=\frac{c}{H_{0}}\frac{d_{l}d_{s}}{d_{ls}}.$
(6)

Therefore, $H_{0}$ and $\Omega_{K,0}$ can be directly determined by comparing
Equation (4) with the time-delay distance ratio $d_{ls}/{d_{s}d_{l}}$
extracted from lensing time delay data through Equations (5) and (6).

To reduce the model dependence and take the full advantage of the whole
$D_{\Delta t}$ measurements, we make use of GRBs as distance indicators to
derive $d_{l}$ and $d_{s}$ out to very high redshifts in Equation (4) via
${{d(z)=H_{0}d_{\rm L}(z)[c(1+z)]^{-1}}}$, which holds in any space-time for
any gravity theory. Due to the lack of low-$z$ GRBs, we also combine the well-
calibrated Pantheon SNe Ia (see Section 3.2) as references to optimize the
constraints on $d_{l}$ and $d_{s}$. In principle, one needs to select a pair
of SC events that match the redshifts of the lens and source in a lensing
system. However, it is difficult for this to be fulfilled with all the
discrete data, and there are always differences in the lensing redshifts from
the nearest SC objects. This can be resolved by smoothing the evolution of
luminosity distances from all discrete events through a model-independent way,
for instance, making use of a polynomial (Räsänen et al., 2015) or non-
parametric approach like the Gaussian process (Shafieloo et al., 2012). It is
worth stressing that, the polynomial $d(z)=\sum_{i=1}^{n}c_{i}z^{i}$
conditional on $d(0)=0$ and $d^{\prime}(0)=1$ adopted in previous works (e.g.,
Collett et al. 2019; Wei & Melia 2020) as the parametrization of SC distances
might suffer from the divergence and convergence problems when fitting to the
high-redshift data sets, and the finite truncated polynomials could result in
some systematic errors. General cosmography provides an ideal tool to explore
the dynamical evolution of cosmic expansion according to the cosmological
principle and without specifying a cosmological model, intrinsically owning
the idea of fitting data using the Taylor series111Analogous approaches have
been proposed to parameterize cosmic distances, such as Padé series (Gruber &
Luongo, 2014) and Chebyshev polynomials (Capozziello et al., 2018). However,
it is very difficult to determine the approximation order due to the dearth of
data points and the qualities of cosmological measurements. expanding around
$z\simeq 0$ with the scale factor derivatives up to a certain order (Visser,
2004; Capozziello et al., 2013, 2019; Aviles et al., 2014).

To avoid the convergence and divergence problem mentioned above, we adopt the
Taylor series expansions, with the scale factor derivative up to the fourth
order, as a function of the y-redshift (Cattoën & Visser, 2008), $y=z/(1+z)$,
for the SC distance parametrizations (see Capozziello et al. 2013 for more
details). In the Appendix, we define the leading-order cosmographic parameters
to be constrained, including the deceleration $q_{0}$, jerk $j_{0}$, snap
$s_{0}$ parameters, etc., and different distances determined by Taylor series
are presented (see Capozziello et al. 2013 for more detailed descriptions
about the evolutionary status of the universe corresponding to these
parameters). Note that the SC distances encode both the $H_{0}$ and
$\Omega_{K,0}$ in the cosmography formulae. Thus our method can adequately
explore how constraining the $H_{0}$ and $\Omega_{K,0}$ are in the DSR
approach.

###  2.2 Bayesian framework

We infer the free parameters within the Bayesian framework (the methods
described below are referred to Trotta 2017). Given the observational data,
${\bm{D}}$, and prior knowledge in modelling the SC distance indicators (see
Section 3.2) and cosmographic models (i.e., the order of expansions), the
posterior probability distribution of the free parameters, ${\bm{\theta}}$, is
given by

$p({\bm{\theta}|{\bm{D}}})=\frac{{\mathcal{L}}({\bm{\theta}};{\bm{D}})\pi({\bm{\theta}})}{\int{{\mathcal{L}}({\bm{\theta}};{\bm{D}})\pi({\bm{\theta}})}d{\bm{\theta}}}.$
(7)

The term ${\mathcal{L}}({\bm{\theta}};{\bm{D}})=p({\bm{D}}|{\bm{\theta}})$ is
the _likelihood_ of observational data conditional on the knowledge of such
free parameters of hypothetic models, which permits to combine different data
sets with the summation of log-likelihoods for the parameter inference. For
the _prior_ $\pi({\bm{\theta}})$, we use uniform distributions by taking large
and physically acceptable ranges for each free parameter unless specified
otherwise, as shown in Table 1. The right-hand-side integral in Equation (7)
is the _evidence_.

The marginalized posterior probability distribution of each parameter is
obtained with implementing the dynamical nested sampling code dynesty
(Speagle, 2020)222https://github.com/joshspeagle/dynesty. Dynamical nested
sampling has desirable merits over the Markov chain Monte Carlo algorithms
with sophisticated treatment of multimodal solutions, more flexible access to
independent samples, and in particular, focusing on effectively and
simultaneously estimating the posterior and evidence while sampling from
complicated, multimodal distributions. Besides, dynamic nested sampling allows
samples to be targeted adaptively to better sample specific areas of the
posterior during the fit. In this work, for each Bayesian parameter, we report
the median of the marginalized posterior probability distribution along with
the $1\sigma$ error bar taken as half of the 16th$-$84th quantile range.

##  3 Data and likelihood constructions

###  3.1 Strong Lensing Time-Delay Distances

We compile the recent six strongly lensed quasars with time-delay measurements
from H0LiCOW collaboration (Wong et al., 2020), including B1608+656 (Suyu et
al., 2010; Jee et al., 2019), RXJ1131-1231 (Suyu et al., 2014; Chen et al.,
2019), HE0435-1223 (Wong et al., 2017; Chen et al., 2019), WFI2033-4723 (Rusu
et al., 2020), J1206+4332 (Birrer et al., 2019), and PG1115+080 (Chen et al.,
2019). All lenses except B1608+656 are blindly analyzed with respect to the
cosmological parameters. The posterior distributions of time-delay distance of
these galaxy lenses are publicly released at the H0LiCOW
website333http://www.h0licow.org. For B1608+656, the time-delay distance
likelihood function was approximated with a skewed log-normal distribution:

${\mathcal{L}}_{D_{\Delta\tau}}=\frac{1}{\sqrt{2\pi}(x-\lambda_{D})\sigma_{D}}\exp\left[{-\frac{(\ln(x-\lambda_{D})-\mu_{D})^{2}}{2\sigma^{2}_{D}}}\right],$
(8)

with the parameters $\mu_{\rm D}=7.0531$, $\sigma_{\rm D}=0.22824$, and
$\lambda_{\rm D}=4000.0$, where $x=D_{\Delta t}$ (in units of Mpc). The time
delay posterior distributions of the remaining five lensed quasars were
released with the samples of Monte-Carlo Markov Chains. We also include
DES-J0408-5354 from STRIDES collaboration (Shajib et al., 2020), which sets
the current most precise measurement of $H_{0}$ from a single time-delay lens.
For this event, we use the time-delay distance posterior distribution derived
by Shajib et al. (2020). Besides, SDSSJ0946+1006 (Gavazzi et al., 2008) is a
double-source-plane strong lensing with the presence of two sources
($z_{s1}=0.609$, $z_{s2}=2.3$) that lensed by the same foreground galaxy
($z_{l}=0.222$), offering an accurate constraint on the cosmological scaling
factor
$\xi=(d_{ls1}d_{s2}/d_{s1}d_{ls2})=(d_{ls1}/d_{l}d_{ls1})\cdot(d_{l}d_{l,s2}/d_{l,s2})$,
which is sensitive to $\Omega_{K,0}$ and with no dependence on $H_{0}$. In
SDSSJ0946+1006, the constrain result is $\xi^{-1}=1.404\pm 0.016$ (Collett &
Auger, 2014). Thus, we can well approximate the likelihood function
${\mathcal{L}}_{\xi^{-1}}$ using a Gaussian function with the mean value and
dispersion $\mu_{\xi^{-1}}\pm\sigma_{\xi^{-1}}=1.404\pm 0.016$, i.e.,

${\mathcal{L}}_{\xi^{-1}}=\frac{1}{\sqrt{2\pi}\sigma_{\xi^{-1}}}\exp\left\\{-\frac{[\xi^{-1}-\frac{T(z_{l})-T(z_{s2})}{T(z_{l})-T(z_{s1})}]}{2\sigma^{2}_{\xi^{-1}}}\right\\},$
(9)

with $T(z)=\sqrt{1/d^{2}(z)+\Omega_{K,0}}$.

Table 1: Priors used in the nested sampling and each free parameter’s descriptions. $\mathcal{U}$ represents a uniform probability distribution. Parameter | Prior | Description of Parameter
---|---|---
$M_{\rm B}$ [mag] |  $\mathcal{U}$(-25,-15) | B-band absolute peak magnitude of SNe Ia.
$A$ |  $\mathcal{U}$(0,100) | Intercept of the Amati correlation of LGRBs in logarithm.
$B$ |  $\mathcal{U}$(0,5) | Slope of Amati correlation of LGRBs in logarithm.
$\delta_{\rm grb}$ |  $\mathcal{U}$(0,2) | Intrinsic dispersion of Amati correlation of LGRBs in logarithm.
$\beta^{\prime}$ |  $\mathcal{U}$(-5,5) |  Parameter related to the intercept ($\beta$) of quasar’s $L_{\rm X}-L_{\rm UV}$ relation in logarithm.
$\gamma$ |  $\mathcal{U}$(0,2) |  Slope of quasar’s $L_{\rm X}-L_{\rm UV}$ relation in logarithm.
$\delta_{\rm qso}$ |  $\mathcal{U}$(0, 2) |  Intrinsic scatter of quasar’s $L_{\rm X}-L_{\rm UV}$ relation in logarithm.
$H_{0}$ [km$\cdot$s${}^{-1}\cdot$Mpc-1] |  $\mathcal{U}$(40, 90) | Present cosmic expansion rate.
$\Omega_{K,0}$ |  $\mathcal{U}$(-1, 1) | Present energy density of spatial curvature component.
$q_{0}$ |  $\mathcal{U}$(-2,0) | Cosmographic deceleration parameter.
$j_{0}$ |  $\mathcal{U}$(-10,10) | Cosmographic jerk parameter.
$s_{0}$ |  $\mathcal{U}$(-150,150) | Cosmographic snap parameter.
Table 2: The bounds obtained with different independent data combinations at
1$\sigma$ credibility.

Data combined | $H_{0}^{*}$ | $\Omega_{K,0}$ | $q_{0}$ | $j_{0}$ | $s_{0}$ | $M_{B}$ |  $A$/$\gamma$ |  $B$/$\beta^{\prime}$ |  $\delta_{\rm grb}$/$\delta_{\rm qso}$
---|---|---|---|---|---|---|---|---|---
Lensing+SN Ia | $75.8^{+3.7}_{-2.5}$ | $0.12^{+0.16}_{-0.14}$ | $-0.64^{+0.24}_{-0.25}$ | $2.6^{+3.9}_{-4.1}$ | $1^{+48}_{-27}$ | $-19.17^{+0.09}_{-0.08}$ |  |  |
Lensing+SN Ia+SH0ES | $74.2^{+0.5}_{-1.4}$ | $0.04^{+0.10}_{-0.10}$ | $-0.56^{+0.25}_{-0.22}$ | $1.0^{+3.8}_{-3.6}$ | $-6^{+39}_{-25}$ | $-19.24^{+0.04}_{-0.02}$ |  |  |
Lensing+SN Ia+$H_{0}^{\star}$ | $76.6^{+2.5}_{-2.7}$ | $0.11^{+0.16}_{-0.11}$ | $-0.69^{+0.29}_{-0.19}$ | $2.7^{+3.7}_{-3.9}$ | $6^{+42}_{-31}$ | $-19.14^{+0.05}_{-0.09}$ |  |  |
Lensing+SN Ia | $73.8^{+2.3}_{-1.8}$ | flat | $-0.61^{+0.23}_{-0.29}$ | $2.1^{+3.6}_{-4.7}$ | $-22^{+61}_{-12}$ | $-19.24^{+0.06}_{-0.05}$ |  |  |
Lensing+QSO | $77.6^{+2.8}_{-3.3}$ | $-0.21^{+0.41}_{-0.49}$ | $-0.41^{+0.28}_{-0.38}$ | $6.3^{+3.6}_{-2.7}$ | $34^{+41}_{-35}$ |  | $0.58^{+0.01}_{-0.01}$ | $-1.63^{+0.04}_{-0.02}$ | $0.22^{+0.01}_{-0.01}$
Lensing+QSO+$H_{0}^{\star}$ | $77.0^{+2.6}_{-2.6}$ | $-0.23^{+0.38}_{-0.50}$ | $-0.37^{+0.26}_{-0.35}$ | $7.6^{+1.8}_{-4.5}$ | $34^{+34}_{-40}$ |  | $0.58^{+0.01}_{-0.01}$ | $-1.63^{+0.04}_{-0.02}$ | $0.22^{+0.00}_{-0.00}$
Lensing+QSO | $77.9^{+2.9}_{-2.5}$ | flat | $-0.45^{+0.34}_{-0.21}$ | $7.2^{+2.3}_{-3.7}$ | $18^{+48}_{-20}$ |  | $0.58^{+0.01}_{-0.01}$ | $-1.62^{+0.03}_{-0.03}$ | $0.22^{+0.01}_{-0.01}$
Lensing+GRB | $71.5^{+4.4}_{-3.0}$ | $-0.07^{+0.13}_{-0.06}$ | $-0.39^{+0.33}_{-0.67}$ | $-1.4^{+5.2}_{-4.9}$ | $96^{+44}_{-97}$ |  | $49.19^{+0.17}_{-0.18}$ | $1.47^{+0.06}_{-0.08}$ | $0.37^{+0.02}_{-0.02}$
Lensing+GRB+$H_{0}^{\star}$ | $72.9^{+3.1}_{-3.0}$ | $-0.02^{+0.17}_{-0.07}$ | $-0.69^{+0.43}_{-0.56}$ | $0.4^{+4.6}_{-5.1}$ | $56^{+59}_{-76}$ |  | $49.20^{+0.17}_{-0.17}$ | $1.46^{+0.07}_{-0.07}$ | $0.37^{+0.02}_{-0.02}$
Lensing+GRB | $70.9^{+4.2}_{-2.9}$ | flat | $-0.43^{+0.33}_{-0.53}$ | $-0.4^{+3.8}_{-5.2}$ | $98^{+46}_{-78}$ |  | $49.21^{+0.15}_{-0.20}$ | $1.47^{+0.06}_{-0.08}$ | $0.37^{+0.03}_{-0.02}$
Lensing+GRB† | $70.4^{+3.5}_{-3.3}$ | $-0.07^{+0.03}_{-0.02}$ | $-0.43^{+0.32}_{-0.62}$ | $-3.3^{+3.4}_{-4.6}$ | $112^{+37}_{-62}$ |  | $49.31^{+0.10}_{-0.13}$ | $1.47^{+0.04}_{-0.04}$ | $0.27^{+0.01}_{-0.02}$
Lensing+GRB†+$H_{0}^{\star}$ | $71.4^{+3.2}_{-2.8}$ | $-0.07^{+0.03}_{-0.02}$ | $-0.75^{+0.51}_{-0.44}$ | $-3.8^{+4.6}_{-3.9}$ | $112^{+37}_{-72}$ |  | $49.31^{+0.11}_{-0.13}$ | $1.47^{+0.03}_{-0.05}$ | $0.27^{+0.01}_{-0.02}$
Lensing+GRB† | $69.9^{+3.0}_{-2.5}$ | flat | $-0.28^{+0.28}_{-0.32}$ | $-2.3^{+3.1}_{-3.2}$ | $133^{+16}_{-50}$ |  | $49.23^{+0.12}_{-0.09}$ | $1.47^{+0.04}_{-0.04}$ | $0.27^{+0.02}_{-0.02}$
Lensing+SN Ia+QSO | $78.0^{+2.9}_{-3.0}$ | $0.14^{+0.16}_{-0.15}$ | $-0.94^{+0.15}_{-0.11}$ | $8.9^{+1.0}_{-2.7}$ | $96^{+12}_{-39}$ | $-19.14^{+0.09}_{-0.07}$ | $0.60^{+0.01}_{-0.01}$ | $-1.56^{+0.01}_{-0.02}$ | $0.22^{+0.00}_{-0.00}$
Lensing+SN Ia+QSO+SH0ES | $74.1^{+0.8}_{-1.0}$ | $0.01^{+0.10}_{-0.10}$ | $-0.90^{+0.20}_{-0.12}$ | $7.8^{+1.8}_{-3.6}$ | $61^{+24}_{-34}$ | $-19.24^{+0.03}_{-0.03}$ | $0.60^{+0.00}_{-0.01}$ | $-1.54^{+0.01}_{-0.01}$ | $0.22^{+0.00}_{-0.00}$
Lensing+SN Ia+QSO+$H_{0}^{\star}$ | $77.0^{+2.7}_{-2.4}$ | $0.11^{+0.15}_{-0.14}$ | $-0.91^{+0.13}_{-0.14}$ | $8.5^{+1.4}_{-2.5}$ | $86^{+22}_{-31}$ | $-19.36^{+0.01}_{-0.02}$ | $0.60^{+0.00}_{-0.01}$ | $-1.52^{+0.01}_{-0.01}$ | $0.23^{+0.00}_{-0.00}$
Lensing+SN Ia+QSO | $75.4^{+2.3}_{-1.8}$ | flat | $-0.98^{+0.18}_{-0.11}$ | $9.1^{+0.9}_{-3.1}$ | $99^{+10}_{-47}$ | $-19.20^{+0.05}_{-0.06}$ | $0.59^{+0.01}_{-0.01}$ | $-1.59^{+0.01}_{-0.01}$ | $0.22^{+0.00}_{-0.00}$
Lensing+SN Ia+GRB | $75.6^{+3.3}_{-2.9}$ | $0.10^{+0.18}_{-0.12}$ | $-0.50^{+0.23}_{-0.24}$ | $0.7^{+3.7}_{-3.5}$ | $-10^{+40}_{-22}$ | $-19.18^{+0.09}_{-0.08}$ | $49.15^{+0.15}_{-0.20}$ | $1.38^{+0.07}_{-0.07}$ | $0.38^{+0.02}_{-0.02}$
Lensing+SN Ia+GRB+SH0ES | $74.0^{+0.7}_{-1.1}$ | $0.04^{+0.11}_{-0.08}$ | $-0.43^{+0.19}_{-0.23}$ | $-0.5^{+3.3}_{-2.9}$ | $-9^{+26}_{-24}$ | $-19.23^{+0.02}_{-0.03}$ | $49.15^{+0.16}_{-0.18}$ | $1.39^{+0.06}_{-0.07}$ | $0.38^{+0.02}_{-0.02}$
Lensing+SN Ia+GRB+$H_{0}^{\star}$ | $75.3^{+3.0}_{-2.5}$ | $0.10^{+0.16}_{-0.12}$ | $-0.52^{+0.28}_{-0.17}$ | $-1.2^{+4.8}_{-2.0}$ | $-17^{+38}_{-17}$ | $-19.19^{+0.08}_{-0.07}$ | $49.12^{+0.17}_{-0.17}$ | $1.40^{+0.05}_{-0.08}$ | $0.38^{+0.02}_{-0.02}$
Lensing+SN Ia+GRB | $73.3^{+2.2}_{-1.8}$ | flat | $-0.46^{+0.24}_{-0.26}$ | $-0.8^{+3.8}_{-3.7}$ | $-23^{+38}_{-18}$ | $-19.23^{+0.04}_{-0.06}$ | $49.17^{+0.15}_{-0.20}$ | $1.38^{+0.08}_{-0.06}$ | $0.38^{+0.02}_{-0.03}$
Lensing+SN Ia+GRB† | $74.6^{+2.5}_{-2.9}$ | $0.07^{+0.14}_{-0.12}$ | $-0.23^{+0.16}_{-0.17}$ | $-4.4^{+2.7}_{-2.0}$ | $-31^{+16}_{-15}$ | $-19.19^{+0.05}_{-0.10}$ | $48.98^{+0.11}_{-0.11}$ | $1.45^{+0.05}_{-0.03}$ | $0.29^{+0.02}_{-0.01}$
Lensing+SN Ia+GRB†+SH0ES | $73.5^{+1.1}_{-0.8}$ | $0.04^{+0.1}_{-0.08}$ | $-0.19^{+0.12}_{-0.18}$ | $-4.6^{+2.6}_{-1.7}$ | $-30^{+14}_{-16}$ | $-19.23^{+0.02}_{-0.04}$ | $49.00^{+0.10}_{-0.12}$ | $1.47^{+0.04}_{-0.04}$ | $0.29^{+0.02}_{-0.01}$
Lensing+SN Ia+GRB†+$H_{0}^{\star}$ | $74.8^{+2.3}_{-2.8}$ | $0.09^{+0.12}_{-0.14}$ | $-0.17^{+0.09}_{-0.22}$ | $-4.1^{+2.3}_{-2.1}$ | $-32^{+16}_{-14}$ | $-19.20^{+0.07}_{-0.08}$ | $48.98^{+0.11}_{-0.11}$ | $1.46^{+0.04}_{-0.04}$ | $0.29^{+0.02}_{-0.02}$
Lensing+SN Ia+GRB† | $73.2^{+1.6}_{-1.8}$ | flat | $-0.19^{+0.13}_{-0.20}$ | $-4.4^{+2.0}_{-2.4}$ | $-35^{+14}_{-14}$ | $-19.24^{+0.04}_{-0.06}$ | $49.00^{+0.11}_{-0.11}$ | $1.46^{+0.04}_{-0.04}$ | $0.29^{+0.02}_{-0.02}$

  * $*$

∗ In units of km$\cdot$s${}^{-1}\cdot$Mpc-1.

  * ⋆

⋆With using $H_{0}=74.5^{+5.6}_{-6.1}$ km$\cdot$s${}^{-1}\cdot$Mpc-1 as a
prior.

  * †

† With 500 mock GRBs combined.

###  3.2 Distance indicators used to derive $d_{l}$ and $d_{s}$

Figure 1: Upper panel: The derived best-fitting $E_{\rm p,z}-E_{\rm iso}$
relation assuming a flat $\Lambda$CDM cosmology with taking $H_{0}=70$
km$\cdot$s${}^{-1}\cdot$Mpc-1 and $\Omega_{\rm m,0}=0.3$ (red dashed line),
1$\delta_{\rm grb}$ (dark gray) and 3$\delta_{\rm grb}$ (light gray)
confidence regions. The black symbols are the observed data and the cyan
points denote the simulated (see section 4.2). Lower panel: The Hubble
diagrams of the standardizable candles implemented in this work, including SNe
Ia (red), quasars(blue), and GRBs(orange). Note that we adopt the quasars with
$z\leq 5$ for our analysis in this work (see Section 3.2). The term
$h_{0}=H_{0}/(100~{}\rm km\cdot s^{-1}\cdot Mpc^{-1})$.

Gamma-Ray Bursts. We compile 193 redshift-known LGRBs from Amati et al.
(2019), and implement the Amati correlation as GRB distance indicators. For
this GRB sample, $\sim 70\%$ of them lie in the $0.5<z<2.5$ range with an
average of $\bar{z}\sim 1.9$. Similar to Amati et al. (2019), we write the
Amati correlation as a power-law function of the cosmological rest-frame peak
energy, $E_{\rm p,z}$, derived from the observer-frame one, $E_{\rm p,obs}$,
via $E_{\rm p,z}=E_{\rm p,obs}(1+z)$, and the prompt $\gamma$-ray energy of
LGRBs, $E_{\rm iso}=4\pi d^{2}_{\rm L}(z,\bm{\Omega})S_{\rm bol}(1+z)^{-1}$,
where $d_{\rm L}$ is the luminosity distance (in units of cm) in terms of an
approximate parametrization or a specific cosmology with free parameters
$\bm{\Omega}$, $S_{\rm bol}$ is the observable bolometric integral flux (in
units of erg$\cdot$cm-1) within the energy band of observing instruments, and
the factor $(1+z)$ is the correction of cosmological time dilation effect.
Taking the logarithm, the $E_{\rm p,z}-E_{\rm iso}$ correlation can be
modelled as

$\log(E_{\rm iso}/{\rm erg})=A+B\;\log(E_{\rm p,z}/\rm keV),$ (10)

with the constants $A$ and $B$ being the coefficient values of Amati
correlation. Equation (10) links the theoretical distance $d_{\rm
L}(z,\bm{\Omega})$ to the observable sets {$z,E_{\rm p,z},S_{\rm bol}$} of
LGRB population. Then the likelihood function is given such that

$\begin{split}{\mathcal{L}}_{\rm
grb}=\prod_{i=1}^{193}\frac{1}{\sqrt{2\pi(\sigma^{2}_{\log(S_{\rm
bol}),i}+B^{2}\sigma^{2}_{\log(E_{\rm p,z}),i}+\delta^{2}_{\rm grb})}}\\\
\times\exp\left[-\frac{\Delta_{{\rm grb},i}^{2}}{2(\sigma^{2}_{\log(S_{\rm
bol}),i}+B^{2}\sigma^{2}_{\log(E_{\rm p,z}),i}+\delta^{2}_{\rm
grb})}\right],\end{split}$ (11)

and

$\Delta_{\rm grb}=\log(S_{\rm bol})-A-B\;\log(E_{\rm p,z})+\log(4\pi
d^{2}_{\rm L}/(1+z)),$ (12)

Here, $\sigma_{\log E_{\rm p,z}}$ ($\sigma_{\log S_{\rm bol}}$) represents the
error propagation (in logarithm) of the observed $E_{\rm p,z}$ ($S_{\rm
bol}$), $\delta_{\rm grb}$ is a free parameter accounting for the intrinsic
dispersion around the best-fitting relation of Equation (10). Throughout we
neglect the uncertainties of spectroscopically measured redshift and
K-correction for taking into account the cosmological redshifting of the GRB
spectra and the limited bandwidth of observing telescopes (Bloom et al.,
2001).

In the upper panel of Figure 1, we present the best-fitting $E_{\rm
p,z}-E_{\rm iso}$ relation (red dashed line) along with 1$\delta_{\rm grb}$
and 3$\delta_{\rm grb}$ confidence regions by assuming the flat $\Lambda$CDM
as a fiducial cosmology ($H_{0}=70$ km$\cdot$s${}^{-1}\cdot$Mpc-1,
$\Omega_{\rm m,0}=0.3$, the same as below), where $A=49.06\pm 0.32$,
$B=1.53\pm 0.12$, and $\delta_{\rm grb}=0.35\pm 0.04$ dex. The derived GRB
Hubble diagram is plotted in the lower panel (see orange points). One can see
that the GRBs distance moduli are consistent with the SN Ia ones at low
redshifts 444The absolute magnitude of all SN Ia objects is fixed to -19.35
mag, which corresponds to $H_{0}=70$ km$\cdot$s${}^{-1}\cdot$Mpc-1, see
Scolnic et al. (2018). , albeit with an intrinsic dispersion of around 0.35
_dex_ in logarithmic space. Besides, Demianski et al. (2021) parametrized the
redshift evolution of Amati parameters in the form of a power-law function and
showed that there is no significant trend systematic with redshift for the
observables, coefficients, as well as $\delta_{\rm grb}$, by cross-comparing
the calibrations in different redshift bins and observing instruments using a
local regression technique. Therefore, we can safely use this sample of LGRBs
as distance indicators.

Pantheon Supernovae Ia. The _Pantheon_ compilation released the current most
complete SN Ia sample, which consists of 1048 objects covering the redshift
range of $0.01\leq z\leq 2.3$ (Scolnic et al., 2018). As shown in the lower
panel of Figure 1, most of the SN Ia events ($\sim 80\%$) are located in the
redshift range $0.2<z<0.6$, and only 6 of them are at $z>1.4$, much
insufficient to be used to cover the redshift range of current strong lensing
systems. The observed SN Ia distance moduli (red points in this panel) are
calculated from the Spectral Adaptive Light curve Template 2 (SALT2) light-
curve fit parameters with using a modified version of the Tripp formula
(Tripp, 1998):

$\hat{\mu}_{\rm obs}=m_{\rm B}-M_{\rm B}+\alpha\cdot X-\eta\cdot C+\Delta_{\rm
M}+\Delta_{\rm B},$ (13)

where $m_{\rm B}$ is the B-band apparent magnitude, $X$ is the stretch factor,
$C$ is the color parameter, $\alpha$ and $\eta$ are the nuisance coefficients
of the stretch- and color-luminosity relations, respectively, and $M_{\rm B}$
is another nuisance parameter that represents the B-band absolute peak
magnitude of a fiducial SNe Ia with $X$ = 0 and $C$ = 0. Furthermore,
$\Delta_{\rm M}$ and $\Delta_{\rm B}$ are the distance corrections based on
the host-galaxy mass and from various biases predicted from simulations,
respectively. For the Pantheon SNe Ia, the values of $\alpha$ and $\eta$ are
calibrated by using the Bayesian Estimation Applied to Multiple Species with
Bias Corrections (BBC) method (Kessler & Scolnic, 2017), whistle correcting
those expected biases (See Marriner et al. 2011 for more details.). With the
BBC technique, Scolnic et al. (2018) reported the apparent magnitudes with
bias corrections ($m_{\rm corr}$) of all SNe Ia. Then the observed distance
modulus is given by directly subtracting $M_{\rm B}$ from $m_{\rm corr}$,
i.e., $\hat{\mu}_{\rm obs}=m_{\rm corr}-M_{\rm B}$. Using the vector of SN
distance residuals $\Delta$m$=\hat{\bm{\mu}}_{\rm obs}-{\hat{\bm{\mu}}}_{\rm
model}$, the likelihood function of the model fit is given by

$-2\ln({\mathcal{L}}_{\rm sn})=\Delta{\bf m}^{T}\cdot{\bf
C}^{-1}\cdot\Delta{\bf m}.$ (14)

Here, ${\bf C}$ is the uncertainty matrix (including both statistical and
systematic uncertainties), $\hat{\bm{\mu}}_{\rm obs}$ and
${\hat{\bm{\mu}}}_{\rm model}$ are respectively the observed and model vector
of SN distance moduli with $\hat{\mu}_{\rm model}=5\log[d_{\rm
L}({z,\bm{\Omega}})/(10~{}{\rm pc})]$.

Actually, in our method SN Ia data can be anchored near their high-redshift
end in combination with $D_{\Delta t}$ measurements from TDGLs (Taubenberger
et al., 2019; Wong et al., 2020), thus we do not take a fiducial value of
$H_{0}$ in SN Ia distances in our fitting.

Quasar Sample. Quasars are of high luminosity emerging at galaxy scales and
have served as a new cosmological probe of the high-redshift universe
(Risaliti & Lusso, 2019; Lusso et al., 2020). Recently, Lusso et al. (2020)
reported a new catalogue of 2421 optically selected quasars with the
measurements of spectroscopic redshifts and X-ray observations from either
Chandra or XMM-Newton. This high-quality sample enables set up the Hubble
diagram up to redshift $z\sim 7.5$ by using the non-linear phenomenological
relation between the X-ray and ultraviolet luminosities of quasars (Avni &
Tananbaum, 1982; Lusso & Risaliti, 2017), written logarithmically as

$\log L_{\rm X}=\beta+\gamma\log L_{\rm UV},$ (15)

where $L_{\rm X}$ and $L_{\rm UV}$ correspond to the monochromatic
luminosities at 2 keV and 2500 ${\rm\AA}$, respectively. Then the distance
indicator of quasars can be re-expressed as a function of the UV and X-ray
fluxes (both in units of erg cm-2 s-1 Hz-1), $F_{\rm UV}$ and $F_{\rm X}$,
luminosity distance, $d_{\rm L}$ at redshift $z$, i.e.,

$\log F_{\rm X}=\beta+\gamma\log F_{\rm UV}+(\gamma-1)\log 4\pi d^{2}_{\rm
L}(z).$ (16)

We can obtain the quasar distances by directly fitting the relation between
$F_{\rm X}$ and $F_{\rm UV}$ as a function of the slope $\gamma$ and intercept
$\beta$. The redshift dependence of $\gamma$ and $\beta$ has been evaluated in
a cosmology-independent way by using sub-samples in narrow redshift bins. With
this approach, Lusso et al. (2020) verified that both $\gamma$ and $\beta$ do
not show significant trend with redshift at $z\leq 5$. We thus use 2404
quasars with redshift $z\leq 5$ in the following analyses. It should be noted
that the intercept $\beta$ is very difficult to be determined due to the lack
of confident physical interpretations for the $L_{\rm X}-L_{\rm UV}$ relation
as well as its intrinsic dispersion, $\delta_{\rm qso}$. Therefore, we set
$\gamma$, $\beta$, and $\delta_{\rm qso}$ as free parameters to be
marginalized over in our fitting, which can also partly alleviate the problem
of the mismatch in the shape of Hubble diagrams of quasars and other probes
(e.g., the SNe Ia). The quasar likelihood function is given by

$\begin{split}{\mathcal{L}}_{\rm
qso}=\prod_{i=1}^{2404}\frac{1}{\sqrt{2\pi(\sigma^{2}_{\log(F_{\rm
X}),i}+\gamma^{2}\sigma^{2}_{\log(F_{\rm UV}),i}+\delta^{2}_{\rm qso})}}\\\
\times\exp\left[-\frac{\Delta_{{\rm qso},i}^{2}}{2(\sigma^{2}_{\log(F_{\rm
X}),i}+\gamma^{2}\sigma^{2}_{\log(F_{\rm UV}),i}+\delta^{2}_{\rm
qso})}\right],\end{split}$ (17)

and

$\Delta_{\rm qso}=\log(F_{\rm X,30})-\beta^{\prime}-\gamma\;\log(F_{\rm
UV,27.5})-2(\gamma-1)\log(d_{\rm L,28.5}(z,{\bm{\Omega}})),$ (18)

where $\beta^{\prime}=\beta+(\gamma-1)\log(4\pi)$, $F_{\rm X}$ and $F_{\rm
UV}$ are normalised (in logarithm) to 30 and 27.5, respectively.

The lower panel of Figure 1 also displays the quasar Hubble diagram (blue
symbols), which is derived by fitting the data with fifth-order logarithmic
polynomials (Lusso et al., 2020). We can see that quasars exhibit a relatively
large dispersion and appear to deviate from the flat $\Lambda$CDM beyond
$z\sim 2$ (see Lusso et al. 2020 for more details). We thus focus on the
cosmographic fits with the combined $D_{\Delta t}$ of TDGLs, SN Ia, and GRB
data. Nevertheless, the combined results with quasars can be employed for
comparison and validation in our model-independent measurements.

##  4 Results and Discussions

###  4.1 Current observational constraints

Figure 2: Marginalized one- and two-dimensional posterior probability
distributions (and contours) at 1$\sigma$ and 2$\sigma$ credibilities produced
by dynamical nested sampling with the combination of TDGLs, SNe Ia, and GRBs
in general cosmography with a curvature component (upper panel) and with
taking a flat universe prior (lower panel).

Figure 3: Same as in Figure 2, but with the combination of TDGLs, SNe Ia, and
quasars.

We first focus on the cosmographic fits by combining TDGLs, SNe Ia, and GRBs.
For all the free parameters, we use uninformative priors, in particular, for
the cosmographic parameters, we adopt the priors used by Li et al. (2020). In
Table 1, we give the priors and detailed descriptions of each free parameter
to be constrained. Table 2 tabulates the best-fitting results and 1$\sigma$
credible intervals produced from a random sampling of the Bayesian formulae
using the dynamic nested sampling code. We display the 1$\sigma$ and 2$\sigma$
marginalized posterior probability distributions (contours) in the
_cornerplots_ of Figure 2. Combining the lensing and low-redshift SNe Ia
yields $H_{0}=75.8^{+3.7}_{-2.5}$ km s-1Mpc-1 and
$\Omega_{K,0}=-0.12^{+0.16}_{-0.14}$, with smaller statistical uncertainties
compared to previous results (Collett et al., 2019), which can be attributed
to the fact that the constraints of $H_{0}$ and $\Omega_{K,0}$ are contributed
by both the $D_{\Delta t}$ and SC distances. This also suggests that the
cosmographic parametrizations with a spatial curvature component can, in some
degrees, consolidate the direct measurements of $H_{0}$ and $\Omega_{K,0}$ via
the DSR. Taking a flat-universe prior, we obtain $H_{0}=73.8^{+2.3}_{-1.8}$ km
s-1Mpc-1, significantly improved by 38%. From the upper panel of Figure 2, we
notice that the degeneracy between $M_{B}$ and $H_{0}$ in SN Ia distances
still exists due to the limited constraints on $H_{0}$ by current $D_{\Delta
t}$ data, which degrades the determinations of $H_{0}$ and then $\Omega_{K,0}$
when combining SNe Ia. To reduce such degeneracy effects, we utilize the prior
of $H_{0}$ with the most recently released measurement by SH0ES team
($H_{0}=73.30\pm 1.04$ km s-1Mpc-1; Riess et al. 2021), and find that the
$\Delta\Omega_{K,0}$ falls by $\sim 40\%$. From the Lensing+GRB combination,
we obtain $H_{0}=71.5^{+4.4}_{-3.0}$ km s-1Mpc-1 and
$\Omega_{K,0}=-0.07^{+0.13}_{-0.06}$, providing a 6% precision measurement of
$H_{0}$ and favoring a flat universe with $\Delta\Omega_{K,0}\approx 0.1$.

The uncertainties of combining Lensing+GRB mainly come from the poor
constraints of cosmographic parameters due to the large dispersion of current
GRB data. In the case of a flat universe, the error bar of $H_{0}$ is slightly
narrowed (by $\sim 5\%$), and the cosmographic parameters are also moderately
improved. In addition, the derived values of deceleration, jerk, and snap
parameters are in agreement with a present accelerating universe in
$\Lambda$CDM cosmology, which approximately corresponds to $q_{0}\approx-0.55$
($\Omega_{\rm m,0}=0.3$), $j_{0}=1$, and $s_{0}\approx 0$ (Capozziello et al.,
2013), although $s_{0}$ is poorly constrained.

SN Ia data possess much better qualities at relatively low redshifts relative
to GRBs (and quasars), and thus can be implemented as references to calibrate
the GRB (quasar) distances in the same redshift range. Adding the SN Ia data
into the Lensing+GRBs, we find that the uncertainty of $H_{0}$ is narrowed by
25% relative to the value determined by Lensing+GRBs and the $H_{0}-M_{B}$
degeneracy can not be improved. The constraint of $\Omega_{K,0}$ is not
improved, which could be the result of the effect of $H_{0}-M_{B}$ degeneracy
in SN Ia distances and the mismatch between SN Ia and GRB Hubble diagram at
$z\gtrsim 1.4$555This can be seen from the constrained values of the slope $B$
and snap parameter between Lensing+GRB and Lensing+SN Ia+GRB, see Table 2..
Such degeneracy effects can be seen from the result by taking the SH0ES prior
on $H_{0}$, which reduces $\Delta\Omega_{K,0}$ by 39%. In our prospects below,
we will show the effects of the uncertainties regarding to GRB distances on
the measurements of $H_{0}$ and $\Omega_{K,0}$ combining a larger sample of
realistic GRBs produced based on a Monte-Carlo simulation.

Figure 3 presents the results by combing Lensing, SNe Ia, and quasars; see
also Table 2. In the case of Lensing+QSO, we obtain a relatively higher value
of $H_{0}$ with $\Delta H_{0}\approx 3$ km s-1Mpc-1; the spatial curvature
density is loosely constrained with $\Delta\Omega_{K,0}\sim 0.5$. Compared to
Lensing+GRB, the looser constraint on $\Omega_{K,0}$ comes from that the
quasar distances ($z\lesssim 5$) have a larger dispersion and larger error
bars at high redshifts (e.g., $z\gtrsim 2.5$; see the lower panel in Figure
1), which results in larger uncertainties in the measurements of $d_{l}$ and
$d_{s}$ and then the $\Omega_{K,0}$. Also, GRB distances have a wider
available redshift coverage than the current quasars, which can provide more
meaningful constraints. Taking a flat-universe prior gives only a slightly
improved value of $H_{0}$. When adding SN Ia data into Lensing+QSO, we obtain
only a slight improvement on the limit of $H_{0}$, whereas the
$\Delta\Omega_{K,0}$ can be reduced by a factor of 3. While the cosmographic
parameters are more constraining relative to Lensing+SN Ia/QSO, their values
deviate from that of $\Lambda$CDM cosmology at beyond 1$\sigma$ credible
intervals. Meanwhile, compared to the results combining GRB distances above,
the outcomes with QSOs also suffer from the $M_{B}-H_{0}$ degeneracy in
fitting SN Ia distances and the mismatch between SN Ia and QSO Hubble
diagrams. Therefore, one should be careful when combining quasar data as
complementary distance indicators. The flat-universe assumption provides a 21%
reduction for the error bar of $H_{0}$. When imposing the SH0ES prior on
$H_{0}$, we obtained a precision similar to the SN Ia and GRB yields.

When adding GRBs to Lensing+SN Ia+QSO, we find no systematic deviation between
QSO and GRB distances. There is no significant improvement in our results
($H_{0}=77.4^{+2.9}_{-3.0}$ km s-1Mpc-1 and
$\Omega_{K,0}=0.13^{+0.17}_{-0.13}$) compared to Lensing+SNIa+QSO and/or
Lensing+SNIa+GRB. The result could be due to the mismatch between the SN Ia
and GRB/QSO distances, which is a sign of some systematic errors in QSO/GRB
distance calibration as using the current $D_{\Delta t}$ measurements for the
anchoring along with using SN Ia distances as a reference. Such uncertainties
can be further analyzed by using more high-quality time-delay lensing systems,
as which “govern” the precision and accuracy in restricting the $H_{0}$. Also,
high-quality and large-statistics SNIa/GRB data sets are needed to find out
the underlying systematical effects.

In Figure 4, we also display the measurements of $H_{0}$ with different data
combinations. The value with combining Lensing+GRB is between the CMB and
SH0ES results with a precision of $\sim$6% (see Table 2). The Lensing+SN Ia
combination results in a value favoured by the SH0ES measurement (pink shaded
region) with a $\sim$3% precision in a flat universe. However, the
combinations of quasars give a relatively high expansion rate. When combining
with SNe Ia, the values coincide with the SH0ES measurement as well as the
lensing-only result in a flat $\Lambda$CDM cosmology ($73.5\pm 1.6$
km$\cdot$s${}^{-1}\cdot$Mpc-1 from this work, see the green symbol) at
1$\sigma$ credible regions.

In our method, the constraints of $H_{0}$ and $\Omega_{K,0}$ rely on the
$d_{l}$ and $d_{s}$ provided by anchoring the pairs of SC objects on the
$D_{\Delta t}$ along with their wide redshift coverage. The measurement of
$D_{\Delta t}$ depends on the assumption in lens modelling, which dominates
the source of systematic uncertainty in measuring $H_{0}$ in time-delay
cosmography. Birrer et al. (2020) considered the mass-sheet transformations
and obtained $H_{0}=74.5^{+5.6}_{-6.1}$ km s-1Mpc-1. We also use this value as
a prior to see the impact of the systematic effects arising from lens
modelling on the limit of $\Omega_{K,0}$. The results are presented in Table
2. One can see that such effects do not significantly influence the bound of
$\Omega_{K,0}$.

Figure 4: The present values of Hubble parameter determined with different
data combinations. The thin and shaded lines denotes the 1$\sigma$ results
with taking a uninformative prior and flat-universe prior on $\Omega_{K,0}$,
respectively. The cyan shaded region denotes the CMB (early-universe)
measurement in flat $\Lambda$CDM cosmology ($H_{0}=67.4\pm 0.5$
km$\cdot$s${}^{-1}\cdot$Mpc-1, Planck Collaboration et al. 2020) and the pink
shaded region shows the (late-universe) determination by local distance ladder
($H_{0}=73.30\pm 1.04$ km$\cdot$s${}^{-1}\cdot$Mpc-1, Riess et al. 2021).

###  4.2 Future prospects

As shown above, the determinations of $H_{0}$ and $\Omega_{K,0}$ via the DSR
suffer from the lack of low-redshift GRB observations and the large intrinsic
dispersion in the GRB Hubble diagram. Such uncertainties, as well as the
effect of $M_{B}-H_{0}$ degeneracy in fitting SN Ia distances, can be improved
with larger statistics of GRB data, although the measurement of $H_{0}$ is
mainly determined by the lensing $D_{\Delta t}$ measurements. Here we perform
a Monte-Carlo simulation to investigate the constraints on the cosmographic
and spatial curvature parameters by combining the current lensing, Pantheon SN
Ia sample, and future GRB observations of the Amati quantities. Our simulation
procedure for producing the $z,E_{\rm p,z}$, and $S_{\rm bol}$ of each GRB
event is as follows. First, we fit the gamma distribution function to the
current spectroscopic redshift data set, and randomly sample the redshifts in
the range $0.001<z<2.5$, which well covers the current lensing systems.
Second, we fit the gaussian function to the $\log(E_{\rm p,z}/\rm keV)$ data,
which gives a mean 2.57 and width 0.45, and $\log(E_{\rm p,z}/\rm keV)$ of
each mock GRB is randomly drawn with subject to such a density function.
Third, we produce the mock $E_{\rm iso}$ samples using the best-fitting Amati
correlation in Equation (10) with considering an intrinsic dispersion
$\delta_{\rm grb}=0.3$ dex, which are determined in a flat $\Lambda$CDM
cosmology (see Section 3.2). The jet collimated corrections of the GRB prompt
$\gamma$-ray energy are neglected, despite the arguments that the jet opening
angles could be narrow as below 10 degrees (e.g., Goldstein et al. 2016).
Finally, we calculate the mock $S_{\rm bol}$ data from the simulated $E_{\rm
iso}$ by evaluating the luminosity distances assuming a flat $\Lambda$CDM
model and select the GRB events with a moderate detector threshold of $1\times
10^{-7}$ erg cm-2 (e.g., at the energy band 20-2000 keV). We specify a common
statistical uncertainty of 30% for the simulated samples of $\\{E_{\rm
p,z},S_{\rm bol}\\}$. By running the Monte-Carlo simulation, we produce 500
realistic GRB events for our following analysis666The expected number of GRBs
observed at the design sensitivities of current (Swift, Gehrels et al. 2004;
and Fermi, Meegan et al. 2009) and future (SVOM, Götz et al. 2009; LOFT,
Feroci et al. 2012; and THUSEUS, Amati et al. 2018) missions dedicated to
observing GRBs in the coming decade’s operations can reach $\sim 600$. , which
are shown in the upper panel of Figure 1 (cyan pints).

Figure 5: Same as in Figure 2, but with the combination of TDGLs, SNe Ia, and
500 mock GRBs.

Figure 5 presents the results by combining the 500 mock GRB events. The
1$\sigma$ credible intervals are also summarized in Table 2. We find that the
restrictions of cosmographic parameters are more stringent compared to that
with the real data and the $\Delta\Omega_{K,0}$ can be reduced to $\sim 0.02$
by using Lensing+mock GRB. The case of using a flat-universe prior further
reduces the error bar of $H_{0}$ by 14% and results in a 4% precision. When
adding the mock GRBs into Lensing+SN Ia, the $M_{B}-H_{0}$ degeneracy can be
moderately eliminated (the contours in the $M_{B}-H_{0}$ plane are narrowed
and the correlation between $M_{B}$ and $H_{0}$ tends to be weaker) and both
the uncertainties of $H_{0}$ and $\Omega_{K,0}$ are narrowed compared to
Lensing+SN Ia and/or Lensing+SN Ia+GRB. Using the prior of $H_{0}$ informed by
SH0ES, we can further eliminate the $M_{B}-H_{0}$ degeneracy in SN Ia
distances and the result of $\Omega_{K,0}$ is more constraining. In the
assumption of a flat universe, the precision of $H_{0}$ can reach $\sim 2$%.

##  5 Summary

We have provided new independent constraints on the $H_{0}$ and $\Omega_{K,0}$
combining the $D_{\Delta t}$ measurements from TDGLs and GRB luminosity
distances, only based on the validity of the FLRW metric and geometrical
optics. This method benefits from the wide redshift coverage of the absolute
distances of lensing systems by using the high-redshift LGRBs, taking the full
advantage of $D_{\Delta t}$ measurements in determining $H_{0}$. We use the
Amati correlation with 193 redshift-known LGRBs to establish the GRB Hubble
diagram up to $z\sim 8.1$, much higher than that of SNe Ia or quasars. The
$H_{0}$ and $\Omega_{K,0}$ are determined by fitting the $D_{\Delta t}$ and
GRB distance indicators through the DSR, in which the GRB distances are
parametrized using the general cosmographic parameters together with a
curvature component. Thus the limitations of $H_{0}$ and $\Omega_{K,0}$ are
contributed by both the combined $D_{\Delta t}$ and SC distances.

We present our main findings below. With the combination of Lensing+GRB, we
obtain $H_{0}=71.5^{+4.4}_{-3.0}$ km s-1Mpc-1 and
$\Omega_{K,0}=-0.07^{+0.13}_{-0.06}$. In a flat universe, the $H_{0}$ is
measured to be $H_{0}=70.9^{+4.2}_{-2.9}$ km s-1Mpc-1, with a 6% precision.
The SN Ia distances can be anchored near their high redshift end in the
combination of Lensing+SN Ia, in which the constraints of $H_{0}$ and
$\Omega_{K,0}$ suffer from the $M_{B}-H_{0}$ degeneracy in SN Ia distances due
to the limited precision of $H_{0}$ provided by $D_{\Delta t}$ measurements of
only six TDGLs. In the combination of Lensing+SN Ia+GRB, current GRB data can
neither eliminate the $M_{B}-H_{0}$ degeneracy when fitting to SN Ia distances
nor effectively improve the constraint on $H_{0}$; the mismatch between the SN
Ia and GRB Hubble diagrams due to some systematic errors in the GRB distance
calibration contributes to a part of uncertainty in restricting
$\Omega_{K,0}$. We performed a Monte-carlo simulation to produce a larger
sample of realistic GRBs, and showed that current Lensing+mock GRB could yield
inspiring improvements in measuring $H_{0}$ (with 4% precision) and
$\Omega_{K,0}$ (with an uncertainty of $\sim 0.02$). Future increment of GRB
observations could help eliminate the $M_{B}-H_{0}$ degeneracy in SN Ia data
and reconcile the SN Ia and GRB Hubble diagrams in the combination of
Lensing+SN Ia+GRB, providing more constraining results of $H_{0}$ and
$\Omega_{K,0}$. In addition, we only include one double-source-plane strong
lensing system in our sample, which provides an accurate measurement of the
cosmological scaling factor that is sensitive to the value of $\Omega_{K,0}$.
The inclusion of more such events will further improve the restrictions on
$\Omega_{K,0}$. In contrast, the resultant values of $\Omega_{K,0}$ by using
quasars are not very constraining compared to SNe Ia and/or GRBs in
conjunction with TDGLs, although they have relatively large statistics. In the
combinations with current lensing, SNe Ia, and GRB data sets, we find no
evidence of significant deviation from a zero spatial curvature or an
accelerating universe, and the values of $H_{0}$ coincide with the result of
local distance ladder with the best accuracy $\Delta H_{0}\approx 2.0$ km
s-1Mpc-1 (3% precision) in a flat universe.

To date, various solutions for the $H_{0}$ tension between the CMB and local
measurements (e.g., any new pre-combination physics or the modification of the
post-combination expansion history of the Universe) are unlikely to be
successful (see e.g., Knox & Millea 2020 and reference therein). In addition,
various independent measurements might suffer from some unknown systematic
effects (see a recent review by Di Valentino et al. 2021). More precise,
model-independent measurements of $H_{0}$ and $\Omega_{K,0}$ via the DSR with
the increases in the number of time-delay lensing systems (in which more
sophisticated lens modelling are expected) and in the statistics of reliable
standardizable candles are admired for verifying and providing some clues to
the solution of the $H_{0}$ tension. Such measurements also give clues to the
nature of dark energy and provide a consistency test of the classical FLRW
metric.

## Acknowledgements

We thank the referee for very helpful comments and suggestions. We are
grateful to L. Amati for providing the GRB data and Elisabeta Lusso for
sharing their latest released quasar sample. Z.-H. Zhu is supported by the
National Natural Science Foundation of China under Grants Nos. 12021003,
11920101003 and 11633001, and the Strategic Priority Research Program of the
Chinese Academy of Sciences, Grant No. XDB23000000. J.-J. Wei is supported by
the NSFC under grant No. U1831122, the Key Research Program of Frontier
Sciences (grant No. ZDBS-LY-7014) of Chinese Academy of Sciences, and the
Natural Science Foundation of Jiangsu Province (grant No. BK20221562). E.-W.
Liang is supported by the National Natural Science Foundation of China (grant
Nos. 12133003). Z.-C. Chen is supported by the National Natural Science
Foundation of China (Grant No. 12247176) and the China Postdoctoral Science
Foundation Fellowship No. 2022M710429.

## Data Availability

The data of strong lensing systems we use are available at
http://www.h0licow.org. The GRB and quasar sample are shared on reasonable
requests to Amati et al. (2019) and Lusso et al. (2020), respectively. The
PYTHON/JUPYTER notebook to reproduce the results of this article will be
shared on reasonable request to the corresponding author.

## References

  * Amati (2006) Amati L., 2006, MNRAS, 372, 233
  * Amati et al. (2002) Amati L., et al., 2002, A&A, 390, 81
  * Amati et al. (2008) Amati L., Guidorzi C., Frontera F., Della Valle M., Finelli F., Landi R., Montanari E., 2008, MNRAS, 391, 577
  * Amati et al. (2018) Amati L., et al., 2018, Advances in Space Research, 62, 191
  * Amati et al. (2019) Amati L., D’Agostino R., Luongo O., Muccino M., Tantalo M., 2019, MNRAS, 486, L46
  * Aviles et al. (2014) Aviles A., Bravetti A., Capozziello S., Luongo O., 2014, Phys. Rev. D, 90, 043531
  * Avni & Tananbaum (1982) Avni Y., Tananbaum H., 1982, ApJ, 262, L17
  * Bernal et al. (2016) Bernal J. L., Verde L., Riess A. G., 2016, J. Cosmology Astropart. Phys., 2016, 019
  * Birrer et al. (2019) Birrer S., et al., 2019, MNRAS, 484, 4726
  * Birrer et al. (2020) Birrer S., et al., 2020, A&A, 643, A165
  * Blandford & Narayan (1992) Blandford R. D., Narayan R., 1992, ARA&A, 30, 311
  * Bloom et al. (2001) Bloom J. S., Frail D. A., Sari R., 2001, AJ, 121, 2879
  * Bonvin et al. (2017) Bonvin V., et al., 2017, MNRAS, 465, 4914
  * Butler et al. (2009) Butler N. R., Kocevski D., Bloom J. S., 2009, ApJ, 694, 76
  * Capozziello et al. (2013) Capozziello S., De Laurentis M., Luongo O., Ruggeri A., 2013, Galaxies, 1, 216
  * Capozziello et al. (2018) Capozziello S., D’Agostino R., Luongo O., 2018, MNRAS, 476, 3924
  * Capozziello et al. (2019) Capozziello S., D’Agostino R., Luongo O., 2019, International Journal of Modern Physics D, 28, 1930016
  * Cattoën & Visser (2008) Cattoën C., Visser M., 2008, Phys. Rev. D, 78, 063501
  * Chen et al. (2019) Chen G. C. F., et al., 2019, MNRAS, 490, 1743
  * Collett & Auger (2014) Collett T. E., Auger M. W., 2014, MNRAS, 443, 969
  * Collett et al. (2019) Collett T., Montanari F., Räsänen S., 2019, Phys. Rev. Lett., 123, 231101
  * D’Amico et al. (2021) D’Amico G., Senatore L., Zhang P., Zheng H., 2021, J. Cosmology Astropart. Phys., 2021, 072
  * Dai et al. (2004) Dai Z. G., Liang E. W., Xu D., 2004, ApJ, 612, L101
  * Demianski et al. (2017) Demianski M., Piedipalumbo E., Sawant D., Amati L., 2017, A&A, 598, A112
  * Demianski et al. (2021) Demianski M., Piedipalumbo E., Sawant D., Amati L., 2021, MNRAS,
  * Denzel et al. (2021) Denzel P., Coles J. P., Saha P., Williams L. L. R., 2021, MNRAS, 501, 784
  * Di Valentino et al. (2020) Di Valentino E., Melchiorri A., Silk J., 2020, Nature Astronomy, 4, 196
  * Di Valentino et al. (2021) Di Valentino E., et al., 2021, Classical and Quantum Gravity, 38, 153001
  * Feroci et al. (2012) Feroci M., et al., 2012, Experimental Astronomy, 34, 415
  * Freedman et al. (2019) Freedman W. L., et al., 2019, ApJ, 882, 34
  * Gavazzi et al. (2008) Gavazzi R., Treu T., Koopmans L. V. E., Bolton A. S., Moustakas L. A., Burles S., Marshall P. J., 2008, ApJ, 677, 1046
  * Gehrels et al. (2004) Gehrels N., et al., 2004, ApJ, 611, 1005
  * Ghirlanda et al. (2006) Ghirlanda G., Ghisellini G., Firmani C., 2006, New Journal of Physics, 8, 123
  * Goldstein et al. (2016) Goldstein A., Connaughton V., Briggs M. S., Burns E., 2016, ApJ, 818, 18
  * Götz et al. (2009) Götz D., et al., 2009, in Meegan C., Kouveliotou C., Gehrels N., eds, American Institute of Physics Conference Series Vol. 1133, Gamma-ray Burst: Sixth Huntsville Symposium. pp 25–30 (arXiv:0906.4195), doi:10.1063/1.3155898
  * Gruber & Luongo (2014) Gruber C., Luongo O., 2014, Phys. Rev. D, 89, 103506
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, ApJ, 873, 111
  * Jee et al. (2016) Jee I., Komatsu E., Suyu S. H., Huterer D., 2016, J. Cosmology Astropart. Phys., 2016, 031
  * Jee et al. (2019) Jee I., Suyu S. H., Komatsu E., Fassnacht C. D., Hilbert S., Koopmans L. V. E., 2019, Science, 365, 1134
  * Kessler & Scolnic (2017) Kessler R., Scolnic D., 2017, ApJ, 836, 56
  * Khadka & Ratra (2021) Khadka N., Ratra B., 2021, MNRAS,
  * Khadka et al. (2021) Khadka N., Luongo O., Muccino M., Ratra B., 2021, arXiv e-prints, p. arXiv:2105.12692
  * Knox & Millea (2020) Knox L., Millea M., 2020, Phys. Rev. D, 101, 043533
  * Li et al. (2020) Li E.-K., Du M., Xu L., 2020, MNRAS, 491, 4960
  * Liang & Zhang (2005) Liang E., Zhang B., 2005, ApJ, 633, 611
  * Liao (2019) Liao K., 2019, Phys. Rev. D, 99, 083514
  * Lusso & Risaliti (2017) Lusso E., Risaliti G., 2017, A&A, 602, A79
  * Lusso et al. (2020) Lusso E., et al., 2020, A&A, 642, A150
  * Marra et al. (2013) Marra V., Amendola L., Sawicki I., Valkenburg W., 2013, Phys. Rev. Lett., 110, 241305
  * Marriner et al. (2011) Marriner J., et al., 2011, ApJ, 740, 72
  * Meegan et al. (2009) Meegan C., et al., 2009, ApJ, 702, 791
  * Nava et al. (2012) Nava L., et al., 2012, MNRAS, 421, 1256
  * Oguri (2019) Oguri M., 2019, Reports on Progress in Physics, 82, 126901
  * Oguri & Marshall (2010) Oguri M., Marshall P. J., 2010, MNRAS, 405, 2579
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Qi et al. (2020) Qi J.-Z., Zhao J.-W., Cao S., Biesiada M., Liu Y., 2020, arXiv e-prints, p. arXiv:2011.00713
  * Räsänen et al. (2015) Räsänen S., Bolejko K., Finoguenov A., 2015, Phys. Rev. Lett., 115, 101301
  * Refsdal (1964) Refsdal S., 1964, MNRAS, 128, 307
  * Riess (2019) Riess A. G., 2019, Nature Reviews Physics, 2, 10
  * Riess et al. (2011) Riess A. G., et al., 2011, ApJ, 730, 119
  * Riess et al. (2016) Riess A. G., et al., 2016, ApJ, 826, 56
  * Riess et al. (2019) Riess A. G., Casertano S., Yuan W., Macri L. M., Scolnic D., 2019, ApJ, 876, 85
  * Riess et al. (2021) Riess A. G., Casertano S., Yuan W., Bowers J. B., Macri L., Zinn J. C., Scolnic D., 2021, ApJ, 908, L6
  * Risaliti & Lusso (2019) Risaliti G., Lusso E., 2019, Nature Astronomy, 3, 272
  * Rusu et al. (2020) Rusu C. E., et al., 2020, MNRAS, 498, 1440
  * Schneider et al. (1992) Schneider P., Ehlers J., Falco E. E., 1992, Gravitational Lenses, doi:10.1007/978-3-662-03758-4.
  * Scolnic et al. (2018) Scolnic D. M., et al., 2018, ApJ, 859, 101
  * Shafieloo et al. (2012) Shafieloo A., Sahni V., Starobinsky A. A., 2012, Phys. Rev. D, 86, 103527
  * Shajib et al. (2019) Shajib A. J., et al., 2019, MNRAS, 483, 5649
  * Shajib et al. (2020) Shajib A. J., et al., 2020, MNRAS, 494, 6072
  * Speagle (2020) Speagle J. S., 2020, MNRAS, 493, 3132
  * Suyu et al. (2010) Suyu S. H., Marshall P. J., Auger M. W., Hilbert S., Blandford R. D., Koopmans L. V. E., Fassnacht C. D., Treu T., 2010, ApJ, 711, 201
  * Suyu et al. (2014) Suyu S. H., et al., 2014, ApJ, 788, L35
  * Suyu et al. (2018) Suyu S. H., Chang T.-C., Courbin F., Okumura T., 2018, Space Sci. Rev., 214, 91
  * Taubenberger et al. (2019) Taubenberger S., et al., 2019, A&A, 628, L7
  * Treu & Marshall (2016) Treu T., Marshall P. J., 2016, A&ARv, 24, 11
  * Tripp (1998) Tripp R., 1998, A&A, 331, 815
  * Trotta (2017) Trotta R., 2017, arXiv e-prints, p. arXiv:1701.01467
  * Velten & Gomes (2020) Velten H., Gomes S., 2020, Phys. Rev. D, 101, 043502
  * Verde et al. (2013) Verde L., Protopapas P., Jimenez R., 2013, Physics of the Dark Universe, 2, 166
  * Verde et al. (2019) Verde L., Treu T., Riess A. G., 2019, Nature Astronomy, 3, 891
  * Visser (2004) Visser M., 2004, Classical and Quantum Gravity, 21, 2603
  * Wang et al. (2015) Wang F. Y., Dai Z. G., Liang E. W., 2015, New Astron. Rev., 67, 1
  * Wei & Melia (2020) Wei J.-J., Melia F., 2020, ApJ, 897, 127
  * Wong et al. (2017) Wong K. C., et al., 2017, MNRAS, 465, 4895
  * Wong et al. (2020) Wong K. C., et al., 2020, MNRAS, 498, 1420
  * Yang et al. (2020) Yang T., Banerjee A., Ó Colgáin E., 2020, Phys. Rev. D, 102, 123532
  * Zhang (2018) Zhang B., 2018, The Physics of Gamma-Ray Bursts, doi:10.1017/9781139226530.

##  Appendix A Appendix

In this appendix, we present the formulas for the parametrizations of the
cosmic distances used in the general cosmography.

The leading-order cosmographic parameters are defined as the scale factor
deviations in terms of cosmic time:

$H(t)\equiv\frac{1}{a}\frac{da}{dt},q(t)\equiv-\frac{1}{aH^{2}}\frac{d^{2}a}{dt^{2}},j(t)\equiv\frac{1}{aH^{3}}\frac{d^{3}a}{dt^{3}},s(t)\equiv\frac{1}{aH^{4}}\frac{d^{4}a}{dt^{4}},...,$
(19)

which are respectively referred to the Hubble, deceleration, jerk, and snap
parameters, etc. With these definitions, we investigate the cosmic expansion
with expanding the scale factor into a Taylor series around the present time
$t_{0}$ ($z=0$),

$\begin{split}\frac{a(t)}{a(t_{0})}=1+H_{0}(t-t_{0})-\frac{1}{2}q_{0}H_{0}^{2}(t-t_{0})^{2}+\frac{1}{6}j_{0}H_{0}^{3}(t-t_{0})^{3}\\\
+\frac{1}{24}s_{0}H_{0}^{4}(t-t_{0})^{4}+...].\end{split}$ (20)

Using the FLRW geometry and the definition of redshift $z=1/a(t)-1$ with
simply assuming $a(t_{0})=1$, we then obtain the dimensionless co-moving
distance as a function of redshift with the fourth-order cosmographic
parameters,

$\displaystyle
d(z)=z-\frac{1}{2}(1+q_{0})z^{2}+\frac{1}{6}(2+4q_{0}+3q_{0}^{2}-j0+\Omega_{K,0})z^{3}$
(21)
$\displaystyle-\frac{1}{24}[6+18q_{0}-27q_{0}^{2}-15q_{0}^{3}-(9+10q_{0})j_{0}$
$\displaystyle-s_{0}+(6+6q_{0})\Omega_{K,0}]z^{4}+{\mathcal{O}}(z^{5}).$

To mitigate the convergence problem that the cosmographic expansion encounters
at higher redshift (e.g., $z>1$), it is useful to recast the $d(z)$ in the
form of the _y-redshift_ , $y=z/(1+z)$. In this manner, the distance in
$z\in(0,\infty)$ is mapped into $y\in(0,1)$, whereby the cosmographic
expansions would be well behaved in the redshift range from the local universe
to the “Big Bang". Thus, with the change of variable $z:\to y$, the
dimensionless co-moving distance becomes

$\displaystyle
d(y)=y+\frac{1}{2}(1-q_{0})y^{2}+\frac{1}{6}(2-2q_{0}+3q_{0}^{3}-j_{0}+\Omega_{K,0})y^{3}$
(22)
$\displaystyle+\frac{1}{24}[6-6q_{0}+9q_{0}^{2}-15q_{0}^{3}-(3-10q_{0})j_{0}$
$\displaystyle+s_{0}+(6-6q_{0})\Omega_{K,0}]y^{4}+{\mathcal{O}}(y^{5}).$

With the relation of $d(z)=H_{0}d_{\rm L}(z)/c(1+z)$, the luminosity distance
in the FLRW universe is given by

$\displaystyle d_{\rm
L}(z)=\frac{c}{H_{0}}\\{z+\frac{1}{2}(1-q_{0})z^{2}-\frac{1}{6}(1-q_{0}-3q_{0}^{3}+j0-\Omega_{K,0})z^{3}$
(23)
$\displaystyle+\frac{1}{24}[2-2q_{0}-15q_{0}^{2}-15q_{0}^{3}+(5+10q_{0})j_{0}$
$\displaystyle+s_{0}-(2+6q_{0})\Omega_{K,0}]z^{4}+{\mathcal{O}}(z^{5})\\}.$

In terms of the _y-redshift_ , $d_{\rm L}(z)$ is converted to

$\begin{split}d_{\rm
L}(y)=\frac{c}{H_{0}}\\{y+\frac{1}{2}(3-q_{0})y^{2}+\frac{1}{6}(11-5q_{0}+3q_{0}^{2}-j0+\Omega_{K,0})y^{3}\\\
+\frac{1}{24}[50-26q_{0}+21q_{0}^{2}-15q_{0}^{3}-(7-10q_{0})j_{0}\\\
+s_{0}+(10-6q_{0})\Omega_{K,0}]y^{4}+{\mathcal{O}}(y^{5})\\}.\end{split}$ (24)

The distance modulus $\mu(z)=5\log[d_{\rm L}(z)/{\rm Mpc}]+25$ can be obtained
from the “logarithmic Hubble relation": $\ln(d_{\rm L}(y)/\,\rm{Mpc})/y=\ln
10/5[\mu(y)-25]-\ln y$, which gives

$\begin{split}\mu(y)=25+\frac{5}{\ln\,10}\\{\ln(c/H_{0})+\ln\,y+\frac{1}{2}(3-q_{0})y+\frac{1}{24}[17-2q_{0}\\\
+9q_{0}^{2}-4(j_{0}-\Omega_{K,0})]y^{2}+\frac{1}{24}[11-q_{0}+2q_{0}^{2}-10q_{0}^{3}-(1-8q_{0})j_{0}\\\
+s_{0}+(16-4q_{0})\Omega_{K,0}]y^{3}+{\mathcal{O}}(y^{4})\\}.\end{split}$ (25)