# Model-independent way to determine the Hubble constant and the curvature
from phase shift of gravitational waves with DECIGO

Tonghua Liu11affiliationmark:  , Shuo Cao22affiliationmark:
3$\ast$3$\ast$affiliationmark:  , Marek Biesiada44affiliationmark:  , Yilong
Zhang22affiliationmark:  33affiliationmark:  , Jieci Wang55affiliationmark:
1\. School of Physics and Optoelectronic Engineering, Yangtze University,
Jingzhou 434023, China;
2\. School of Physics and Astronomy, Beijing Normal University, Beijing
100875, China;
3\. Institute for Frontiers in Astronomy and Astrophysics, Beijing Normal
University, Beijing 102206, China;
4\. National Centre for Nuclear Research, Pasteura 7, PL-02-093 Warsaw,
Poland;
5\. Department of Physics, and Collaborative Innovation Center for Quantum
Effects and Applications, Hunan Normal University, Changsha 410081, China
$\ast$ Corresponding author:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this Letter, we propose a model-independent method to determine the Hubble
constant and curvature simultaneously taking advantage of the possibilities of
future space-borne gravitational wave (GW) detector DECIGO in combination with
the radio quasars as standard rulers. Similarly to the redshift drift in the
electromagnetic domain, accelerating expansion of the Universe causes a
characteristic phase correction to the gravitational waveform detectable by
DECIGO. Hence, one would be able to extract the Hubble parameter $H(z)$. This
could be used to recover distance-redshift relation supported by the data not
relying on any specific cosmological model. Assuming the FLRW metric, and
using intermediate luminosity radio quasars as standard rulers one achieves an
interesting opportunity to directly assess $H_{0}$ and $\Omega_{k}$
parameters. To test this method we simulated a set of acceleration parameters
achievable by future DECIGO. Based on the existing sample of 120 intermediate-
luminosity radio-quasars calibrated as standard rulers, we simulated much
bigger samples of such standard rulers possible to obtain with VLBI. In the
case of $(N=100)$ of radio quasars, which is the size of currently available
sample, the precision of cosmological parameters determined would be
$\sigma_{H_{0}}=2.74$ $\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ and
$\sigma_{\Omega_{k}}=0.175$. In the optimistic scenario $(N=1000)$ achievable
by VLBI, the precision of $H_{0}$ would be improved to $1\%$, which is
comparable to the result of $\sigma_{H_{0}}=0.54$
$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ from _Planck_ 2018 TT, TE,
EE+lowE+lensing data, and the precision of $\Omega_{k}$ would be 0.050. Our
results demonstrate that such combined analysis, possible in the future, could
be helpful to solve the current cosmological issues concerning the Hubble
tension and cosmic curvature tension.

###### Subject headings:

cosmological parameters – cosmology: observations

##  1\. Introduction

Over the past few decades, the observations of type Ia supernova (SN Ia)
revealed that the expansion of the universe is accelerating [51, 52]. Based on
the cosmological principle (supporting homogeneous and isotropic universe) and
Einstein’s general theory of relativity (GR), it is generally accepted that
this acceleration is caused by dark energy with negative pressure. The
simplest solution is to use the cosmological constant term $\Lambda$ as a
model of dark energy driving the accelerated expansion of the universe. With
addition of cold dark matter, it forms the $\Lambda$CDM model. This model has
been supported by the vast majority of astronomical observations over the past
few decades [43, 6, 4, 3]. However, with the accumulation of precise
astrophysical observations, some problems gradually emerged. For instance,
there is a $4.4\sigma$ tension between the Hubble constant ($H_{0}$)
measurements inferred within the $\Lambda$CDM model from cosmic microwave
background (CMB) anisotropy (both temperature and polarization) data [44] and
the local measurements of the Hubble constant by the Supernova $H_{0}$ for the
Equation of State collaboration (SH0ES) [53]. This inconsistency may be caused
by unknown systematic errors in astrophysical observations or unrevealed new
physics which is significantly different from $\Lambda$CDM [19, 15]. Moreover,
some works suggested that besides the $H_{0}$ tension, spatial curvature of
the universe is still an open issue [15, 16, 22]. More specifically, combining
the Planck temperature and polarization power spectra data, Planck
Collaboration et al. [44] claimed that closed universe
($\Omega_{k}=-0.044^{+0.018}_{-0.015}$) was supported by pure CMB data.
However, the combination of Planck lensing data and low redshift baryon
acoustic oscillations (BAO), revealed that flat universe with
$\Omega_{k}=0.0007\pm 0.0019$ was preferred. It should be emphasized that both
methods depend, to some extent, on a particular cosmological model assumed,
i.e. the non-flat $\Lambda$CDM. To better understand inconsistency between the
Hubble constant and the curvature of universe inferred from local measurements
and Planck observations, it is necessary to seek a completely model-
independent method that can constrain both $H_{0}$ and $\Omega_{k}$
simultaneously [11, 68, 49, 36, 35].

Traditionally, if the Hubble parameter $H(z)$ and the luminosity distance
$D_{L}(z)$ or angular diameter distance $D_{A}(z)$ are known at the same
redshift $z$, one can directly measure $H_{0}$ and $\Omega_{k}$
simultaneously. The difficulty, however, is that possibilities of measuring
Hubble parameters in the electromagnetic window are limited. Current cosmic
chronometer sample providing such measurements via differential ages of
passively evolving galaxies consists of only of 31 $H(z)$ data. Fortunately,
gravitational wave (GW) observations provide such a chance to obtain $H(z)$
directly. Recently, the discovery of the first GW event, GW150914 [1], caused
a great excitement in astronomy and physics communities. GW signals from
inspiraling and merging compact binaries enable absolute measurements of the
luminosity distance $D_{L}(z)$ [56], hence they are often referred to as
standard sirens. Subsequent detection of the GW170817 event [2] powered by the
merger of binary neutron stars (NSs) was accompanied by electromagnetic
counterpart detected in gamma-rays, X-rays, in the optical and radio bands. If
the redshifts to standard sirens can be independently determined from their
electromagnetic (EM) counterparts, such GW signals can be used for cosmology.
More importantly, if the adiabatic inspiral phase could be monitored long
enough, the accelerating expansion of the universe would cause an additional
phase shift (more details in Sec. 2.1) in the gravitational waveform, which
would allow us to measure the cosmic acceleration directly [72]. This idea is
similar to the redshift drift [55, 38] but requires shorter cadence of
observations. Current LIGO-Virgo-Kagra (LVK) network of detectors has
negligible possibility to detect cosmic acceleration because chirp signals in
their frequency range are present from seconds to a couple of minutes. The
future generation of space-borne DECi-hertz Interferometer Gravitational-wave
Observatory (DECIGO) will be a game changer, because the inspiral signals,
ultimately to be detected by LVK ground based detectors, could be discovered
and monitored a few years before the final coalescence. Thus, a cosmological
model-independent way to provide $H(z)$ with covering a wide redshift range in
GW domain is possible [58].

Inspired by above mentioned ideas, we will determine the curvature and the
Hubble constant simultaneously using the Hubble expansion rate from simulated
future space-borne detector DECIGO in combination with the “angular size-
distance” relation of ultra-compact structure in radio quasars (QSOs) from the
very-long-baseline interferometry (VLBI) observations. The advantage of our
work consists in assuming only the homogeneity and isotropy of the universe,
thus, it is a geometric and cosmological model-independent way for the
measurement of Hubble constant and spatial curvature simultaneously. The paper
is organized as follows: In Sec. 2, we present the methodology and simulated
data. The results and discussion are given in Sec. 3. Finally, we summarize
our findings in Sec. 4.

##  2\. Methodology and simulation of GW and QSO data.

Assuming that the universe is homogeneous and isotropic on large scales, the
space-time of our universe can be described by Friedmann-Lemaître-Robertson-
Walker (FLRW) metric [66]:

$ds^{2}=-c^{2}dt^{2}+a^{2}(t)\left[\frac{dr^{2}}{1-kr^{2}}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})\right],$ (1)

where $t$ is the cosmic time and $r$, $\theta$, $\phi$ are the comoving
spatial coordinates. The scale factor $a(t)$ is the only gravitational degree
of freedom and its evolution is determined by matter and energy content in the
universe. The dimensionless curvature $k=-1,0,+1$ corresponds to open, flat
and closed universes, respectively. With such a metric, the angular diameter
distance $D_{A}(z)$ can be expressed as [66]

${D_{A}}(z)=\left\\{\begin{array}[]{lll}\frac{D_{H}}{(1+z)\sqrt{|\Omega_{\rm
k}|}}S_{k}\left(\sqrt{|\Omega_{\rm
k}|}\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime})}\right),\\\
\frac{D_{H}}{(1+z)}S_{k}\left(\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime})}\right),\\\
\frac{D_{H}}{(1+z)\sqrt{|\Omega_{\rm k}|}}S_{k}\left(\sqrt{|\Omega_{\rm
k}|}\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime})}\right),\\\
\end{array}\right.$ (2)

where $D_{H}=c/H_{0}$ is known as the Hubble distance, the dimensionless
Hubble parameter $E(z)$ is defined as $H(z)/H_{0}$, where $H(z)$ is the
universe expansion rate at redshift $z$ and $H_{0}=H(z=0)$ is the Hubble
constant. The curvature parameter $\Omega_{k}$ is related to $k$ as
$\Omega_{k}=-c^{2}k/(a_{0}H_{0})^{2}$, where $c$ is the speed of light. For
convenience, we denoted $S_{k}(x)=\sin(x),x,\sinh(x)$ for
$\Omega_{k}<0,=0,>0$. If the angular diameter distance $D_{A}(z)$ and Hubble
parameter $H(z)$ are known, then we can constrain both $H_{0}$ and
$\Omega_{k}$. Standard procedure is to assume the cosmological model like
$\Lambda$CDM or wCDM (i.e. the cold dark matter model with dark energy not
being exactly cosmological constant) and use $H(z)$ formula suggested by the
Friedman equation to calculate theoretical value of angular luminosity
distance $D_{A}^{th}(z;\mathbf{p})$ at the distance corresponding to the
standard ruler, whose observed distance $D_{A}^{obs}(z)$ can be measured. Then
the chi-square based likelihood is maximized (or chi-square minimized) to
obtain the parameters $\mathbf{p}$, in such case $\mathbf{p}$ comprise also
cosmological model parameters like $\Omega_{m},\Omega_{\Lambda},etc.$. In this
work, we use the opportunity to measure the acceleration parameter $X(z)$
(more details see Sec.2.1, Eq. (5) and thereafter) from an additional phase
shift in the gravitational waveform based on the expected performance of the
future DECIGO space-based detector. This would enable to obtain directly the
Hubble parameter $H(z)$ at the redshift of NS-NS binary monitored with DECIGO.
Having the large sample of acceleration parameters $X(z_{i}),\;i=1,...,N$ we
use Artificial Neural Network (ANN) [65] to reconstruct the $X(z)$ function
supported only by the data, without assuming any particular cosmological
model. This reconstructed relation (with the corresponding uncertainty band)
can be used to calculate the distance Eq. (2) to every standard ruler used to
measure the distance. In this case the set of cosmological parameters
$\mathbf{p}$ comprise only $H_{0}$ and $\Omega_{k}$. In the role of standard
rulers, we will use the ultra-compact radio quasars. Detailed description of
GW data and radio quasars is given below.

###  2.1. The acceleration parameter from gravitational waveform detected by
DECIGO

In an expanding Universe, infinitesimal time intervals $dt_{o}$ measured in
the observers frame are related to the analogous time intervals $dt_{e}$ in
source (emitter) frame by the following (well known) relation:
$\frac{dt_{o}}{dt_{e}}=\frac{a(t_{o})}{a(t_{e})}=\frac{1}{a(t_{e})}=1+z.$ This
relation underlies the notion of the cosmological redshift. The next order
derivative reads:
$\frac{d^{2}t_{o}}{dt_{e}^{2}}=(1+z)^{2}\left[{\dot{a}}(t_{o})-{\dot{a}}(t_{e})\right]=(1+z)\left[(1+z)H_{0}-H(z)\right]$
and underlies the phenomenon of the redshift drift. Its order of magnitude is
the inverse of the age of the Universe, hence this correction is usually being
neglected in most cosmological considerations. In the context of GW signals
from inspiralling binaries, the relevant argument of the wave strain $h(t)$ is
the time to coalescence $\Delta t=t_{c}-t$ with $t_{c}$ denoting coalescence
time, measured in the observers frame. One can relate it to the time to
coalescence $\Delta t_{e}$ measured in the source frame by: $\Delta t=\Delta
T+X(z)\Delta T^{2}$, where $\Delta T=\Delta t_{e}(1+z)$ is the usual first
order correction due to cosmological time dilation and the second order
(Taylor expansion) factor has been denoted $X(z)$ and is called the
acceleration parameter (more details in the text following Eq.(5)).
Calculating Fourier transformed ${\tilde{h}}(f)$ and using stationary phase
approximation [72] it is straightforward to see that accelerating expansion
introduces the phase shift $\Psi_{acc}=-2\pi fX(z)\Delta T^{2}$, which is
additional with respect to the first order treatment neglecting the
acceleration. 111One should not be mislead that $\Psi_{acc}$ is simply
proportional to $f$ due to frequency dependence of $\Delta T$. Rigorous
dependence on $f$ is revealed in Eq. (4). See Fig. 1 of Kawamura et al. [27]
for an artistic illustration of this effect, i.e., GW signals coming from a
neutron star binary at a given redshift with and without the acceleration of
the expansion of the Universe. It was demonstrated [58, 41] that in a
gravitational wave signal emitted by a neutron star binary at the redshift of
$z=1$, accelerated expansion of the Universe would cause a phase delay of
about l sec during the 10-year observing time.

Waveform templates used in matched-filtering technique do not account for
acceleration. Hence, with $\sim 10^{7}$ orbital cycles monitored during 10
years the accumulated phase shift due to acceleration would be detectable. As
a second-generation space-based GW detector, DECIGO [26] is designed to
improve the detection sensitivity of GW at lower frequencies, with its most
sensitive frequencies between 0.1 and 10 Hz. It can be achieved by the
proposed four clusters of spacecrafts, each of which having triangular 1000 km
Fabry-Perot Michelson interferometer. This design improves the determination
of GW sources’ position in the sky, as well as the detection sensitivity of GW
at lower frequencies. The expected sensitivity of DECIGO is $\sim 4\times
10^{-24}$ Hz-1/2 for one cluster, which enables the early detection of
inspiraling sources, and more importantly, the direct measurement of the
acceleration of the expansion of the universe. The angular resolution is $\sim
1$ arcsec2, which enables the unique identification of the host galaxy of the
binary, and improves the precision of distance measurements [58]. Therefore,
DECIGO will provide us a new perspective on studying acceleration of cosmic
expansion and on gravitational wave cosmology [zhang2022·].

Besides, benefiting from the precise time prediction of coalescence time and
high resolution of sky location of GW sources, redshift determination from EM
follow-up observations would be more reliable and frequent. DECIGO is expected
to observe $10^{6}$ GW events coming from the neutron star binaries up to the
redshift $z\sim 5$ [27, 41]. For our purpose, we follow the work of Holz &
Hughes [23], where they assumed that the GW event randomly occurs in any one
of the galaxies, and considering the large number of spectroscopic and
photometric redshift measurements in future galaxy surveys (such as
JDEM/WFIRST [20, 31], BigBOSS [57]/DESI [13, 14] or 4MOST [39]), they
concluded that 10,000 GW events from DECIGO with redshift determinations could
be obtained by EM follow-up observations. Considering the GW signal from NS-NS
binary system whose components have masses $m_{1}$, $m_{2}$, one can define
the redshifted chirp mass $\mathcal{M}_{z}\equiv M(1+z)\eta^{3/5}$, where
$M=m_{1}+m_{2}$ is the total mass of the binary system,
$\eta=m_{1}m_{2}/M^{2}$ is the symmetric mass ratio, and $z$ is the source
redshift. In the frequency domain, gravitational waveform including the
effects of the cosmic acceleration can be written as

$\tilde{h}(f)_{\mathrm{acc}}=\tilde{h}(f)e^{i\Psi_{\mathrm{acc}}(f)},$ (3)

where

$\Psi_{\mathrm{acc}}(f)=-\Psi_{N}(f)\frac{25}{768}X\left(z\right)\mathcal{M}_{z}x^{-4},$
(4)

is the phase correction due to the cosmic acceleration, where
$x\equiv(\pi\mathcal{M}_{z}f)^{2/3}$ and
$\Psi_{N}(f)\equiv\frac{3}{128}(\pi\mathcal{M}_{z}f)^{-5/3}$. Terms with
$x^{-n}$ represent the $n$ \- th Post Newtonian (PN) order, hence
$\Psi_{\mathrm{acc}}(f)$ is of “4 - PN” order and becomes more important at
lower frequencies. The acceleration parameter $X(z)$ is defined as [58, 60]

$X(z)\equiv\frac{H_{0}}{2}\left(1-\frac{H(z)}{(1+z)H_{0}}\right).$ (5)

Note such a parameter can be rewritten as $X(z)=[\dot{a}(0)-\dot{a}(z)]/2$.
Thus it directly reflects whether the universe expansion has accelerated
($X(z)>0$) or decelerated ($X(z)<0$). Angular diameter distance, given by Eq.
(2) can be expressed in terms of the acceleration parameter in the following
form:

${D_{A}}(z)=\left\\{\begin{array}[]{lll}\frac{D_{H}}{(1+z)\sqrt{|\Omega_{\rm
k}|}}f\left(\sqrt{|\Omega_{\rm
k}|}\frac{H_{0}}{1+z}\int_{0}^{z}\frac{dz^{\prime}}{H_{0}-2X(z^{\prime})}\right),\\\
\frac{c}{(1+z)^{2}}f\left(\int_{0}^{z}\frac{dz^{\prime}}{H_{0}-2X(z^{\prime})}\right),\\\
\frac{D_{H}}{(1+z)\sqrt{|\Omega_{\rm k}|}}f\left(\sqrt{|\Omega_{\rm
k}|}\frac{H_{0}}{1+z}\int_{0}^{z}\frac{dz^{\prime}}{H_{0}-2X(z^{\prime})}\right).\\\
\end{array}\right.$ (6)

The gravitational waveform without cosmic acceleration is given by

$\tilde{h}(f)=\frac{\sqrt{3}}{2}\mathcal{A}f^{-7/6}e^{i\Psi(f)}\left[\frac{5}{4}A_{\mathrm{pol},\alpha}(t(f))\right]e^{-i\left(\varphi_{\mathrm{pol},\alpha}+\varphi_{D}\right)},$
(7)

where the amplitude

$\mathcal{A}=\frac{1}{\sqrt{30}\pi^{2/3}}\frac{{\mathcal{M}_{z}}^{5/6}}{D_{L}},$
(8)

includes the luminosity distance $D_{L}$, which indicates that it can be
measured directly. For the phase $\Psi(f)$, we use the restricted-2PN [28]
waveform including spin-orbit coupling at the 1.5PN order, which is given by

$\displaystyle\Psi(f)$ $\displaystyle=$ $\displaystyle 2\pi
ft_{c}-\phi_{c}-\frac{\pi}{4}+\frac{3}{128}(\pi\mathcal{M}f)^{-5/3}$ (9)
$\displaystyle\times$
$\displaystyle\left[1-\frac{5}{84}\mathcal{S}^{2}\bar{\omega}x^{-1}-\frac{128}{3}\beta_{g}\eta^{2/5}x\right.-\frac{2355}{1462}I_{e}x^{-19/6}$
$\displaystyle+$
$\displaystyle\left(\frac{3715}{756}+\frac{55}{9}\eta\right)x-4(4\pi-\beta)x^{3/2}$
$\displaystyle+$
$\displaystyle\left(\frac{15293365}{508032}+\frac{27145}{504}\eta\right.\left.\left.+\frac{3085}{72}\eta^{2}-10\sigma\right)x^{2}\right],$

where the direction of orbital angular momentum and the direction of the
source from the Sun are also taken into account [70]. The polarisation
amplitude $A_{\mathrm{pol},\alpha}(t)$ and the polarisation phase
$\varphi_{\mathrm{pol},\alpha}(t)$ represent the waveform measured by each
detector ($\alpha=\mathrm{I},\mathrm{II}$). The Doppler phase $\varphi_{D}(t)$
denotes the difference between the phase of the wave front at the detector and
the phase of the wave front at the solar system barycenter. All these
quantities are explicitly given in Yagi & Tanaka [70].

The waveform $\tilde{h}(f)_{\mathrm{acc}}$ in Eq. (3) depends on 11
parameters:

$\theta^{i}=\left(\ln\mathcal{M}_{z},\ln\eta,\beta,t_{c},\phi_{c},\bar{\theta}_{\mathrm{S}},\bar{\phi}_{\mathrm{S}},\bar{\theta}_{\mathrm{L}},\bar{\phi}_{\mathrm{L}},D_{L},X\right),$
(10)

where $t_{c}$ and $\phi_{c}$ represent the time of and the phase at the final
coalescence, $\beta$ is the spin-orbit coupling parameter.
$(\bar{\theta}_{\mathrm{S}},\bar{\phi}_{\mathrm{S}})$ indicates the direction
of the source in the solar system barycenter frame, and
$(\bar{\theta}_{\mathrm{L}},\bar{\phi}_{\mathrm{L}})$ specifies the direction
of the orbital angular momentum. For simplicity, we set
$m_{1}=m_{2}=1.4M_{\odot}$ and take $t_{c}=\phi_{c}=\beta=0$. For each
fiducial redshift $z$, we randomly generate $10^{4}$ sets of
$(\bar{\theta}_{\mathrm{S}},\bar{\phi}_{\mathrm{S}},\bar{\theta}_{\mathrm{L}},\bar{\phi}_{\mathrm{L}})$,
and for each set, we calculate the uncertainties of the acceleration parameter
$X$ estimated, using the Fisher matrix formalism:

$\Gamma_{ij}=4Re\int_{f_{\mathrm{min}}}^{f_{\mathrm{max}}}\frac{\partial_{i}\widetilde{h}^{*}(f)\partial_{j}\widetilde{h}(f)}{S_{h}(f)}df,$
(11)

where partial derivatives are with respect to the parameters and $S_{h}(f)$ is
the DECIGO noise power spectrum, which analytical expression is [25, 27, 71]:

$\displaystyle S_{h}(f)$ $\displaystyle=$ $\displaystyle 6.53\times
10^{-49}\big{[}1+(\frac{f}{7.36Hz})^{2}\big{]}+4.45\times 10^{-51}$ (12)
$\displaystyle\times$
$\displaystyle(\frac{f}{1Hz})^{-4}\times\frac{1}{1+(\frac{f}{7.36Hz})^{2}}$
$\displaystyle+$ $\displaystyle 4.94\times
10^{-52}\times(\frac{f}{1Hz})^{-4}\rm{Hz^{-1}}.$

The lower cut-off of frequency in Eq. (10)
$f_{\mathrm{min}}=(256/5)^{-3/8}\pi^{-1}\mathcal{M}_{z}^{-5/8}\Delta
T_{obs}^{-3/8}$ is related to the observation duration $\Delta T_{obs}$. The
upper cut-off frequency of the detector $f_{\mathrm{max}}$ is set equal to 100
Hz. In addition, since there will be eight uncorrelated interferometric
signals from DECIGO, the Fisher matrix above should be multiplied by 8 [26,
71]. One can use Fisher analysis to estimate the measurement accuracy of the
acceleration parameter $X$ by marginalizing other parameters in the Fisher
matrix

$\sigma_{X}=8^{-1/2}\left[\left(\Gamma^{-1}_{jj}\right)^{1/2}\right],$ (13)

where the $\sigma_{X}$ is the root-mean-square uncertainty for the $j$-th
parameter (the acceleration parameter $X$), and the $\Gamma^{-1}_{jj}$ is the
covariance matrix element for parameter $X$ calculated from the Eq. (10).

Figure 1.— Simulated acceleration parameters $X(z)$ based on the future DECIGO
and the corresponding uncertainties $\sigma_{X}$ are shown in grey. Inset
figure shows the effect of binning. Red bars represent weighted mean of $X(z)$
in redshift bins of width $\Delta z=0.1$ with the corresponding uncertainties.
Black solid line with corresponding blue uncertainty band display the ANN
reconstructed $X(z)$ relation. The advantage of binning is noticeable.

Finally, the redshifts of GW sources are sampled according to the merger rate
of double compact objects that reflects the star formation history [17], which
is called “rest frame rates” in the cosmological
scenario222http:www.syntheticuniverse.org. In our simulation, we assume flat
$\Lambda$CDM as the proposed fiducial model with the cosmological parameters
derived by _Planck 2018_ measurements [44]. From these simulated results ANN
technique [65] was used to reconstruct the $X(z)$ function and the
corresponding uncertainty band in the redshift range $z\in[0,5]$. Machine
learning methods have become very popular in astrophysics and cosmology.
Motivated by successful applications of ANN [30, 64, 50], we chose this
method, which is a general method that could reconstruct a function from any
kind of data without assuming any particular parameterization of the function.
Moreover, this non-parametric technique does not assume Gaussianity of
measured random variables, hence it is a completely data driven approach.
Based on the future generation of space-borne GW detector DECIGO, we simulated
a dataset of $10^{4}$ measurements of the acceleration parameter $X(z)$.
However, the acceleration parameters have large uncertainties, because there
are too many unknown parameters in the waveform $\tilde{h}(f)_{\mathrm{acc}}$,
which leads to large uncertainties in the estimation from the Fisher matrix.
Therefore, we divided this sample into 50 bins and adopted ANN to reconstruct
the acceleration parameter as a function of redshift. We trained the ANN
network on the simulated $X(z)$ data, and then predicted the $X(z)$ at other
redshifts (where there were no GW data available). In general, the artificial
neural network consists of an input layer, one or more hidden layers and an
output layer. In this work, the input of the neural network is the redshift
$z$, while the output is the corresponding cosmic acceleration parameter
$X(z)$ and its respective uncertainty $\sigma_{X}$ at that redshift. More
details of data reconstruction using ANN can be found in [64, 30, 50]. Our
strategy of non-parametric reconstruction of $X(z)$ function implies that one
does not need to match the QSOs and NS-NS sources by redshift, but rather for
each QSO the corresponding value of acceleration parameter will be
reconstructed from the total sample of GW signals. The values of the
acceleration parameters $X(z)$ and their corresponding uncertainties are shown
in Fig. 1. Let us remark that individual acceleration parameters have
considerable uncertainties, hence we binned the data in the redshift bins of
width of width $\Delta z=0.1$ and in each redshift bin we calculated the
weighted mean,
$X(z)=\frac{\Sigma_{i=1}^{50}(X_{i}(z)/\sigma^{2}_{X_{i}(z)})}{\Sigma_{i=1}^{50}1/\sigma^{2}_{X_{i}(z)}}$,
$\sigma^{2}_{{X}(z)}=\frac{1}{\Sigma_{i=1}^{50}1/\sigma^{2}_{X_{i}(z)}}$. The
ANN network was subsequently trained on these $n=50$ bins. The final
uncertainties level is similar as in the recent work [59], which used
convolutional neural network (CNN) to derive the cosmic acceleration
parameters.

###  2.2. Angular diameter distances from radio quasars

Quasars, as the brightest sources in the universe, have still a great
potential to be used as useful cosmological probes, despite of high observed
dispersion and extreme variability. Observations of these objects have long
been sought to study cosmology beyond the limits of supernovae. With the
advances in observational technology over recent decades, quasars can be
observed at multiple wave-bands, and great efforts have been made to use them
as standard candles or standard rulers in cosmology, using the Baldwin effect
[5], the Broad Line Region radius-luminosity relation [62], the properties of
highly accreting quasars [63], etc. Recently, Risaliti & Lusso [54] attempted
to use quasars as standard candles through the non-linear relation between
their intrinsic UV and X-ray emission. However, the large dispersion still
remains the greatest problem in calibrating the QSOs as standard candles. In
this work, we focus on the “angular size-distance” relation of ultra-compact
structure in radio quasars both observed so far and those expected from the
future VLBI observations.

Figure 2.— Simulated sample of $1000$ angular size measurements of ultra-
compact structure in radio quasars from the future VLBI observations.

With the signal received at multiple radio telescopes across the Earth’s
surface, together with the registered correlated intensities considering the
different arrival times at various facilities, the characteristic angular size
of a distant radio quasar is defined as [21]

$\theta={2\sqrt{-\ln\Gamma\ln 2}\over\pi B},$ (14)

where $B$ represents the interferometer baseline measured in wavelengths, and
$\Gamma=S_{c}/S_{t}$ is the visibility modulus which is defined by the ratio
between the total and correlated flux densities. In a more detailed study
involving the compact structure of radio quasars, Gurvits [21] showed that the
dispersion in linear size is greatly mitigated by retaining only the sources
with a flat spectrum ($-0.38<\alpha<0.18$). Subsequently, the application of
compact radio sources to cosmology has made a breakthrough [24, 61, 32]. With
gradually refined selection technique and the elimination of possible
systematic errors, Cao et al. [7] compiled a sample of 120 intermediate-
luminosity quasars ($10^{27}W/Hz<L<10^{28}W/Hz$) selected from 613
milliarcsecond ultra-compact radio sources covering the redshift range
$0.0035<z<3.787$ in 2.29 GHz VLBI survey. More importantly, they demonstrated
that intermediate-luminosity compact radio quasars can be used to obtain
constraints on cosmological parameters. The details regarding the angular
sizes $\theta(z)$ and redshifts $z$ of 120 compact radio sources are listed in
Table 1 of Cao et al. [7], which extends the Hubble diagram to the redshift
range $0.46<z<2.76$. The angular size of the compact structure in radio QSOs
can be written as

$\theta(z)=\frac{l_{m}}{D_{A}(z;\mathbf{p})},$ (15)

where $D_{A}(z)$ is the angular diameter distance at redshift $z$, and
cosmological model parameters $\mathbf{p}$ contain the Hubble constant $H_{0}$
and other parameters entering the the dimensionless expansion rate expansion
rate $E(z)$. The parameter $l_{m}$ denotes the linear size scaling factor
describing the apparent distribution of radio brightness within the core,
which needs to be calibrated with external indicators such as SN Ia or cosmic
chronometers. In this analysis, we adopt the calibrated value of this
cosmological standard ruler at 2.29 GHz as $l_{m}=11.03\pm 0.25$ pc, through a
cosmological model-independent technique based on cosmic chronometer
measurements [8]. This sample of VLBI ultra-compact structure in radio
quasars, which was extended to multi-frequency analysis [9], has been widely
applied in various cosmological applications at high-redshifts [29, 74, 46,
47, 69, 33, 34, 37].

Figure 3.— The 2-D plots and 1-D marginalized distributions with 1-$\sigma$
and 2-$\sigma$ contours of cosmological parameters ($H_{0}$ and $\Omega_{k}$)
from the simulated QSO and GW data. Dashed lines represent the fiducial
parameters that we used in our simulation, i.e.
$H_{0}=67.4\;\rm{km~{}s^{-1}Mpc^{-1}}$ and $\Omega_{k}=0$.

Considering that currently no other sample of calibrated radio quasars is
available, we further simulate the possible future data on compact radio
quasars from VLBI surveys based on better uv-coverage. Our simulation is based
on the flat $\Lambda$CDM with $H_{0}=67.4\;\rm{km~{}s^{-1}Mpc^{-1}}$ and
$\Omega_{m}=0.308$ adopted from the _Planck_ results [44], i.e. the same
$\Lambda$CDM model as used for GW simulations. Taking the linear size of
intermediate-luminosity quasars as $l_{m}=11.03\pm 0.25$ pc and following the
redshift distribution of QSOs from [42], we simulate 1000 “angular size-
redshift” data representative of compact radio sources in the redshift range
$0.50<z<5.00$. The fractional statistical uncertainty of the angular size
$\theta$ was taken at a level of $3\%$, which is resonable for both current
and future VLBI surveys based on better uv-coverage [45]. Considering the
intrinsic variance in the size of the AGN cores, derived from high angular
resolution observations, an additional 10% systematical uncertainty in the
observed angular sizes was also assumed in our analysis. In addition, we
assume that the measured angular sizes follow a Gaussian distribution
$\theta_{\mathrm{means}}=\mathcal{N(\theta_{\mathrm{fid}},\sigma_{\theta})}$,
where $\theta_{\mathrm{fid}}$ is obtained from Eq. (15) under assumed fiducial
cosmological model. Simulated data are shown in Fig. 2. It should be remarked
that sources with $z<0.5$ are ignored, because the era of quasar formation has
ended, and the nature of quasars has changed dramatically. This problem is
further expanded and elaborated in [21], which indicates that only the high
redshift portion of the radio quasars can be used as a standard rulers. More
details of the specific procedure of QSO simulation can be found in Cao et al.
[10], Qi et al. [47, 48].

##  3\. Results and discussion

In order to obtain the best-fit values of $H_{0}$ and $\Omega_{k}$ with
corresponding uncertainties, we used the Markov Chain Monte Carlo (MCMC)
method [18] to minimize the $\chi^{2}$ objective function:

$\chi^{2}=\sum\limits_{i=1}^{1000}\frac{\left[D_{A,i}^{th}(z;\mathbf{p})-{D_{A,i}^{obs}(z)}\right]^{2}}{\sigma_{D_{A,i}^{th}}^{2}+\sigma_{D_{A,i}^{obs}}^{2}}\;,$
(16)

where the p represents the Hubble constant ${H}_{0}$ and cosmic curvature
$\Omega_{k}$, the observed angular diameter distance is inferred from angular
size-distance relation of ultra-compact structure in QSOs
$D_{A,i}^{obs}(z)=l_{m}/\theta_{i}(z)$. Theoretical counterpart,
$D_{A,i}^{th}(z;\mathbf{p})$ is derived from Eq. (6) where the ANN
reconstructed $X(z)$ function was used to calculate respective integrals
applying a simple trapezoidal rule. The uncertainties
$\sigma_{D_{A}^{th}}^{2}$ and $\sigma_{D_{A}^{obs}}^{2}$ were calculated
through the standard uncertainty propagation formula.

Combining the simulated QSO data from future VLBI observations together with
the acceleration parameters from DECIGO registered GW events, we obtained the
best-fit values of
$H_{0}=67.19^{+0.77}_{-0.74}$$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ and
$\Omega_{k}=0.022^{+0.051}_{-0.046}$ at 68.3% confidence level. The posterior
joint distribution function along with marginal distributions is shown in Fig.
3. Our results are consistent with the fiducial values of $H_{0}=67.4$
$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ and $\Omega_{k}=0$ assumed, which proves
the robustness of our simulation. Compared with the previous model-independent
works based on the cosmic chronometer measurements $H(z)$ and distances from
popular cosmological probes, our results have a higher precision than the
result of combined Pantheon SN Ia plus $H(z)$ giving $\Omega_{k}=0.63\pm 0.34$
[64], and the result of UV+X-ray quasars plus $H(z)$ yielding
$\Omega_{k}=-0.92\pm 0.43$ [67]. More importantly, such combination of radio
and GW windows may achieve constraints with much higher precision on $H_{0}$,
namely with $1\%$ uncertainty. Moreover, it creates an opportunity to estimate
the Hubble constant and cosmic curvature using the objects at much higher
redshifts ($z\sim 5$) than currently available. Therefore the future
measurements of this kind would shed light on the Hubble tension, which is
recognized as an issue of local vs. CMB probes.

The above mentioned result was obtained assuming the optimistic scenario of
having $N=1000$ intermediate luminosity radio quasars calibrated as standard
rulers. Therefore, we investigated the uncertainties of best fitted parameters
as a function of QSO sample size $N$. Let us stress that the bottleneck is the
QSO sample, while DECIGO is expected to provide enough data to extract
acceleration parameters allowing for ANN reconstruction of $X(z)$ relation.
The results are displayed in Fig. 4-5, in which two cases are considered
separately. The first is the case of $H_{0}$ and $\Omega_{k}$ treated as free
parameters. The second one is the precision of one parameter with a fixed
prior on another one. Regarding the first case, one can see that even in the
most conservative case ($N=50$), the result of $\sigma_{H_{0}}=4.36$
$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ i.e. $6\%$ precision is encouraging for
using the proposed cosmological model-independent method. When the number of
standard rulers increases to $N=300$, the precision of $\sigma_{H_{0}}=1.36$
$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ is comparable to the result of
$\sigma_{H_{0}}=1.3$ $\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ given by the SH0ES
collaboration [53]. If the sample size increases to $N=1000$, the precision on
$H_{0}$ will be improved to $\sigma_{H_{0}}=0.77$, which is comparable with
the result of $\sigma_{H_{0}}=0.54$ $\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ from
_Planck_ 2018 TT, TE, EE+lowE+lensing data [44].

Figure 4.— The precision on $H_{0}$ as a function of the number of simulated
QSOs. The blue polyline represents the case of free $\Omega_{k}$ parameter.
The red polyline represents the case of fixed $\Omega_{k}=0$ parameter. For
the illustration of the uncertainty floor observed, standard statistical
$1/\sqrt{N}$ curves have been overplotted.

The joint fit of the $\Omega_{k}$ parameter, in the optimistic scenario
($N=1000$) has the uncertainty of $\sigma_{\Omega_{k}}=0.050$. This is not
competitive to the recent results from _Planck_ lensing data, where a flat
universe was precisely constrained with a $1\sigma$ uncertainty of
$\sigma_{\Omega_{k}}=0.002$ [44]. However, some recent works [15, 16, 22]
focused on the Hubble tension and other tensions in $\Lambda$CDM cosmology
have raised the issue of the value of a spatial curvature of the Universe.
More specifically, combining the _Planck_ temperature and polarization power
spectra data, these works showed that a closed universe
($\Omega_{k}=-0.044^{+0.018}_{-0.015}$) was supported by CMB data, which
provided an explanation of anomalous high lensing amplitude [16] reported in
_Planck_ Legacy 2018. The robust conclusion of these works is to emphasize the
importance of alternative and cosmological model- independent methods to
constrain both $H_{0}$ and $\Omega_{k}$. Our method is one of this kind and
the results are encouraging. The precision would be increased if instead of
two-parameter fit, one would fit a single parameter with the fixed prior on
another one, which will be discussed below in more details.

Figure 5.— The precision on $\Omega_{k}$ as a function of the the number of
the simulated QSOs. The blue polyline represents the free $H_{0}$ parameter
case. The red polyline represents the fixed $H_{0}=67.4$
$\rm{km~{}s^{-1}Mpc^{-1}}$ parameter case. For the illustration of the
uncertainty floor observed, standard statistical $1/\sqrt{N}$ curves have been
overplotted.

Assuming the flat Universe, as shown by the red solid polyline in Fig. 4, the
constraint on $H_{0}$ will be improved significantly. In this case, the
constraint from only $N=100$ QSOs would have uncertainty of
$\sigma_{H_{0}}=1.08~{}\rm km\ s^{-1}\ Mpc^{-1}$, which is slightly better
than the result of $\sigma_{H_{0}}=1.3~{}\rm km\ s^{-1}\ Mpc^{-1}$ reported by
the SH0ES collaboration. Only $N=500$ QSOs are required to produce the results
comparable with those obtained by _Planck_ 2018 TT, TE, EE+lowE+lensing data.
In the most optimistic scenario ($N=1000$), we can get a result of
$\sigma_{H_{0}}=0.35$ $\rm km\ s^{-1}\ Mpc^{-1}$, which is at $0.5\%$ level.
All these results demonstrate that the combination of GW+QSO data proposed by
us could become an important cosmological probe, precise enough to address the
Hubble tension issue. With the fixed prior of $H_{0}=67.4~{}\rm km\ s^{-1}\
Mpc^{-1}$, the constraints on $\Omega_{k}$ are also considerably improved.
This can be understood in terms of degeneracy between $\Omega_{k}$ and
$H_{0}$, which can also be seen from the joint constraint confidence regions
of $H_{0}$ and $\Omega_{k}$ shown in Fig.3. In the most conservative case
($N=50$) the uncertainty of $\Omega_{k}$ is $\sigma_{\Omega_{k}}=0.100$, while
in the most optimistic scenario ($N=1000$), $\sigma_{\Omega_{k}}=0.023$. This
is not as good as the result from the _Planck_ 2018 data [44], yet it refers
to an alternative measurement using extragalactic objects in more local
Universe than CMB. An interesting point is that as the number of QSOs
increases, the constraint on $\Omega_{k}$ does not improve significantly after
reaching $N=500$. From the statistical point of view, when the number of
observations increases by $N$, the precision should increase by $\sqrt{N}$.
The saturation observed suggests the existence of a feature in the analysis
that prevents from going lower. We have checked how the hypothetical
assumption of zero uncertainty for the angular diameter distance influence the
uncertainty of cosmological parameters inferred. It turned out that it has
negligible effect $\sigma_{H_{0}}$ and $\sigma_{\Omega_{k}}$ were only
slightly different from the above mentioned values (at the second significant
digit). Better understanding of the origin of such an uncertainty floor would
be desired and is left for future studies. Moreover, in future applications of
our method with the real data, the influence of systematic effects cannot be
ignored and should be studied carefully.

Our method is built on two pillars: GWs as a source of $H(z)$ data and QSOs as
source of $D_{A}(z)$ data. In order to gain some understanding of relative
importance of uncertainties related to GW signals and QSO, we did the
following test. Namely, based on $N=1000$ QSO data paired with GW
reconstructed $X(z)$ values, we left the uncertainty of QSO data unchanged but
doubled the uncertainty of GW inference. The results we obtained in this case
were: $\sigma_{H_{0}}=1.18$ $\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ and
$\sigma_{\Omega_{k}}=0.070$. Similarly, we kept the uncertainty of GW
unchanged but doubled the uncertainty from QSOs. We obtained the precision of
Hubble constant to be $\sigma_{H_{0}}=1.04$
$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ and the precision of curvature parameter
to be $\sigma_{\Omega_{k}}=0.079$. Fig. 6 displays the $1\sigma$ and $2\sigma$
confidence level contours for these two cases for estimated parameters. In
conclusion, we found that the contribution of QSO and GW data uncertainties to
the precision of parameter constraints were comparable.

Figure 6.— 1D and 2D marginalized probability distributions of $H_{0}$ and
$\Omega_{k}$ as in Fig.3 but corresponding to doubling the uncertainties from
GW or QSO data respectively, while keeping the other unchanged. Dashed lines
represent the fiducial parameters that we used in our simulation, i.e.
$H_{0}=67.4\;\rm{km~{}s^{-1}Mpc^{-1}}$ and $\Omega_{k}=0$.

Finally, we tested the influence of the GW data on the fits. For this purpose
we simulated $N=1000$ GW sources, reconstructed corresponding $X(z)$ function
and analysed jointly with $N=1000$ QSOs. The result was $H_{0}=65.69\pm 1.71$
$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ and $\Omega_{k}=0.03\pm 0.142$ at 68.3%
confidence level. One can see that with an order of magnitude smaller sample
of GW signals the uncertainties raised $2-3$ times relative to the main case
studied. The precision would be at the present level attainable by local
cosmological probes. Hence, the GW signal sample size and quality are crucial.

##  4\. Conclusion

In this Letter, we proposed a cosmological model-independent method to
determine the Hubble constant and curvature parameter simultaneously based on
the future generation of space-borne GW detector DECIGO in combination with
the data regarding intermediate-luminosity radio-quasars calibrated as
standard rulers obtained from the VLBI observations. Since the accelerating
expansion of the Universe would cause an additional phase shift in the
gravitational waveform, we can extract the Hubble parameter $H(z)$ useful for
further cosmological studies. On the other hand, the future VLBI surveys will
observe a large amount of compact structures in radio quasars, which would
provide the angular diameter distances $D_{A}(z)$ based on a simple “angular
size-distance” relation. Such combination of observations performed in the
radio and GW windows creates a valuable opportunity to directly test $H_{0}$
and $\Omega_{k}$ at much higher precision.

To fully demonstrate the potential of such combination of GW+QSO data-set, we
simulated respective data representative of the expected yield of the DECIGO
and VLBI. Simulation of the future VLBI data on QSOs was based on the
properties of currently available sample of $N=120$ intermediate-luminosity
radio-quasars calibrated as standard rulers [7]. The same fiducial cosmology
was assumed in both GW and QSO simulations. In the most optimistic case of
$N=1000$ QSO calibrated as standard rulers, joint fits for $H_{0}$ and
$\Omega_{k}$ were consistent with the input simulation values and their
$68.3\%$ confidence level uncertainties were $\sigma_{H_{0}}=0.75$
$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ and $\sigma_{\Omega_{k}}=0.048$,
respectively. It means that achievable precision is at the level of $1\%$.
Besides the joint fit on two parameters simultaneously, we considered fits of
one parameter with fixed priors on the second one. In such cases precision
improved as could be expected. In particular the precision of $H_{0}$ with a
fixed prior on $\Omega_{k}$ can reach $0.5\%$ in the most optimistic case of
$N=1000$ QSOs. We also discussed the performance of the method on QSO samples
of different size (from $N=50$ to $N=1000$). In particular regarding
$\Omega_{k}$ fit the uncertainty behavior did not reflect simple $1/\sqrt{N}$
behaviour, but saturated at $N=500.$ This phenomenon remains to be fully
understood, but urges a careful consideration of systematic effects in future
applications of the proposed method to the real data.

As a final remark, there are many potential ways to improve our results. In
the future, multi-frequency VLBI observations will yield more high-quality
quasar observations based on better uv coverage [45]. Radio quasars having a
compact structure, will be measured with higher angular resolution, and lower
statistical and systematic uncertainty. On the other hand, we also look
forward to a large amount of future data, not only from the radio quasars, but
also from other astronomical probes, such as SNe Ia and BAO, allowing us to
further improve the precision of $H_{0}$ and $\Omega_{k}$ measurements. The
results we obtained are very encouraging regarding the method proposed, which
is alternative and complementary to CMB analysis. In light of current tensions
in cosmology when local probes suggest different values of important
parameters like the Hubble constant, from those obtained by CMB analysis,
every cosmological probe alternative to already existing is valuable.

The authors are grateful the referee for careful reading the manuscript and
numerous insightful comments. Liu T.-H. was supported by National Natural
Science Foundation of China under Grant No. 12203009; Chutian Scholars Program
in Hubei Province (X2023007). Cao S. was supported by Beijing Natural Science
Foundation No. 1242021; National Natural Science Foundation of China under
Grants Nos. 12021003, 11633001, and 11920101003; the Strategic Priority
Research Program of the Chinese Academy of Sciences, Grant No. XDB23000000;
and the Interdiscipline Research Funds of Beijing Normal University. Biesiada
M. was supported by Hubei Province Foreign Expert Project (2023DJC040). Wang
J. C. was supported by the National Natural Science Foundation of China under
Grant No. 12122504 and 12035005.

## References

  * Abbott et al. [2016] Abbott, B. P., et al. [LIGO Scientific and Virgo Collaborations], 2016, PRL, 116, 061102
  * Abbott et al. [2017] Abbott, B. P., et al. [LIGO Scientific Collaboration, the Virgo Collaboration], 2017, PRL, 119, 161101
  * Abbott et al. [2018] Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, PRD, 98, 043526
  * Anderson et al. [2014] Anderson, L., Aubourg, É., Bailey, S., et al. 2014, MNRAS, 441, 24. doi:10.1093/mnras/stu523
  * Baldwin [1977] Baldwin, J. A. 1977, ApJ, 214, 679
  * Caldwell et al. [1998] Caldwell, R. R., Dave, R., & Steinhardt, P. J. 1998, PRL, 80, 1582
  * Cao et al. [2017] Cao, S., Zheng, X., Biesiada, M., et al. 2017, A&A, 606, A15
  * Cao et al. [2017] Cao, S., Biesiada, M., Jackson, J., et al. 2017, JCAP, 2017, 012
  * Cao et al. [2018] Cao, S., Biesiada, M., Qi, J., et al. 2018, EPJC, 78, 749
  * Cao et al. [2019] Cao, S., Qi, J., Biesiada, M., et al. 2019, Physics of the Dark Universe, 24, 100274
  * Collett et al. [2019] Collett, T., Montanari, F., & Räsänen, S. 2019, PRL, 123, 231101
  * Cutler & Holz [2009] Cutler, C. & Holz, D. E., 2009, PRD, 80, 104009
  * DESI Collaboration [2016] DESI Collaboration, Aghamousa, A., Aguilar, J., Ahlen, S. et al. “The DESI Experiment Part I: Science,Targeting, and Survey Design” (2016) arXiv:1611.00036
  * DESI Collaboration [2023] DESI Collaboration, Adame, A. G., Aguilar, J., Ahlen, S., et al. (2023) “The Early Data Release of the Dark Energy Spectroscopic Instrument” arXiv:2306.06308
  * Di Valentino et al. [2021] Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021, Astroparticle Physics, 131, 102605
  * Di Valentino et al. [2020] Di Valentino, E., Melchiorri, A., & Silk, J. 2020, Nature Astronomy, 4, 196
  * Dominik et al. [2013] Dominik, M., Belczynski, K., Fryer, C., et al. 2013, ApJ, 779, 72
  * Foreman-Mackey et al. [2013] Foreman-Mackey, D., Hogg, D. W., Lang, D., et al. 2013, PASP, 125, 306
  * Freedman [2017] Freedman, W. L. 2017, Nature Astronomy, 1, 0169
  * Gehrels [2010] Gehrels, N., 2010, arXiv:1008.4936
  * Gurvits [1994] Gurvits, L. I. 1994, ApJ, 425, 442
  * Handley [2021] Handley, W. 2021, PRD, 103, L041301
  * Holz & Hughes [2005] Holz, D. E., & Hughes, S. A. 2005, ApJ, 629, 15
  * Jackson & Dodgson [1997] Jackson, J. C. & Dodgson, M. 1997, MNRAS, 285, 806
  * Kawamura et al. [2006] Kawamura, S., Nakamura, T., Ando., M. et al. 2006, CQG, 23, 125
  * Kawamura et al. [2011] Kawamura, S., et al. 2011, CGQ, 28, 094011
  * Kawamura et al. [2019] Kawamura, S., Nakamura, T., Ando., M. et al. 2019, IJMPD, 28, 1845001
  * Kidder et al. [1993] Kidder., L. E., Will, C. M., & Wiseman, A. G. 1993, PRD, 47, 3281
  * Li et al. [2017] Li, X., Cao, S., Zheng, X., et al. 2017, arXiv:1708.08867
  * Liu et al. [2021] Liu, T., Cao, S., Zhang, S. et al. 2021, Eur.Phys.J. C, 81, 903
  * Levi et al. [2011] Levi, M. E., Kim, A. G., Lampton, M. L., Sholl, M. J. 2011, arXiv:1105.0959
  * Lima & Alcaniz [2002] Lima, J. A. S. & Alcaniz, J. S. 2002, ApJ, 566, 15
  * Liu et al. [2020b] Liu, T., Cao, S., Zhang, J., et al. 2020b, MNRAS, 496, 708
  * Liu et al. [2021a] Liu, T., Cao, S., Zhang, S., et al. 2021a, EPJC, 81, 903
  * Liu et al. [2023] Liu, T., Liu, Z., Wang, J., et al. 2023, JCAP, 2023, 059. doi:10.1088/1475-7516/2023/07/059
  * Liu et al. [2022] Liu, T., Cao, S., Biesiada, M., et al. 2022, ApJ, 939, 37. doi:10.3847/1538-4357/ac93f3
  * Liu et al. [2023] Liu, T., Cao, S., Ma, S., et al. 2023, PLB, 838, 137687.doi:10.1016/j.physletb.2023.137687
  * Loeb [1998] Loeb, A. 1998, ApJL, 499, L111. doi:10.1086/311375
  * de Jong et al. [2019] de Jong, R.s., et al., 2019, The Messenger, 175, 3
  * Nishizawa et al. [2011] Nishizawa, A., Taruya, A., & Saito, S., 2011, PRD, 83, 084045
  * Nishizawa [2012] Nishizawa, A., et al. 2012, PRD, 85, 044047
  * Palanque-Delabrouille et al. [2016] Palanque-Delabrouille, N., Magneville, C., Yèche, C., et al. 2016, A&A, 589, C2
  * Planck Collaboration et al. [2016] Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
  * Planck Collaboration et al. [2020] Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Pushkarev & Kovalev [2015] Pushkarev, A. B. & Kovalev, Y. Y. 2015, MNRAS, 452, 4274
  * Qi et al. [2017] Qi, J.-Z., Cao, S., Biesiada, M., et al. 2017, EPJC, 77, 502
  * Qi et al. [2019a] Qi, J.-Z., Cao, S., Biesiada, M., et al. 2019a, PRD, 100, 023530
  * Qi et al. [2019b] Qi, J.-Z., Cao, S., Zhang, S., et al. 2019, MNRAS, 483, 1104
  * Qi et al. [2021] Qi, J.-Z., Zhao, J.-W., Cao, S., et al. 2021, MNRAS, 503, 2179
  * Qi et al. [2023] Qi, J.-Z., Meng, P., Zhang, J.-F.& Zhang, X., 2023, PRD, 108, 063522
  * Riess et al. [1998] Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009. doi:10.1086/300499
  * Riess et al. [2007] Riess, A. G., Strolger, L.-G., Casertano, S., et al. 2007, ApJ, 659, 98
  * Riess et al. [2019] Riess, A. G., Casertano, S., Yuan, W., et al. 2019, ApJ, 876, 85
  * Risaliti & Lusso [2019] Risaliti, G. & Lusso, E. 2019, Nature Astronomy, 3, 272
  * Sandage [1962] Sandage, A. 1962, ApJ, 136, 319. doi:10.1086/147385
  * Schutz [1986] Schutz, B. F., 1986, Nature, 323, 310
  * Schlegel et al. [2010] Schlegel, D. J., Bebek, C., Heetderks, H. et al. 2010, arXiv:0904.0468
  * Seto et al. [2001] Seto, N., Kawamura, S., & Nakamura, T. 2001, PRL, 87, 221103
  * Sun et al. [2024] Sun, M.-F., Li, J., Cao, S. & Liu, X. 2024, A&A, 682, A177
  * Takahashi & Nakamura [2005] Takahashi, R., Nakamura, T. 2005, Prog. Theor. Phys. 113, 63
  * Vishwakarma [2001] Vishwakarma, R. G. 2001, CQG, 18, 1159
  * Watson et al. [2011] Watson, D., Denney, K. D., Vestergaard, M., et al. 2011, ApJL, 740, L49
  * Wang et al. [2013] Wang, J.-M., Du, P., Valls-Gabaud, D., et al. 2013, PRL, 110, 081301
  * Wang et al. [2020a] Wang, G. J., Ma, X. J., Xia, J. Q., 2020a, MNRAS, 501, 5714
  * Wang et al. [2020b] Wang, G.-J., Ma, X.-J., Li, S.-Y., et al. 2020b, ApJS, 246, 13. doi:10.3847/1538-4365/ab620b
  * Weinberg [2008] Weinberg, S. 2008, Cosmology, by Steven Weinberg. ISBN 978-0-19-852682-7. Published by Oxford University Press, Oxford, UK, 2008
  * Wei et al. [2020] Wei, J. J., Melia, F., 2020, ApJ, 888, 99
  * Wei & Melia [2020] Wei, J.-J. & Melia, F. 2020, ApJ, 897, 127
  * Xu et al. [2018] Xu T., Cao, S., Qi, J., et al. 2018, JCAP, 06, 042
  * Yagi & Tanaka [2010] Yagi, K., & Tanaka, T. 2010, PRD, 81, 064008
  * Yagi & Seto [2011] Yagi, K. & Seto, N. 2011, PRD, 83, 044011
  * Yagi et al. [2012] Yagi, K., Nishizawa, A. & Yoo, C.-H., 2012, JCAP, 04, 031
  * Zhang et al. [2022] Zhang, Y., Cao, S., Liu, X., et al. 2022, ApJ, 931, 119
  * Zheng et al. [2017] Zheng, X., Biesiada, M., Cao, S., et al. 2017, JCAP, 2017, 030