# Exploring bulk viscous unified scenarios with Gravitational Waves Standard
Sirens

Weiqiang Yang<EMAIL_ADDRESS>Department of Physics, Liaoning Normal
University, Dalian, 116029, P. R. China  Supriya Pan
<EMAIL_ADDRESS>Department of Mathematics, Presidency
University, 86/1 College Street, Kolkata 700073, India  Institute of Systems
Science, Durban University of Technology, PO Box 1334, Durban 4000, Republic
of South Africa  Eleonora Di Valentino<EMAIL_ADDRESS>School
of Mathematics and Statistics, University of Sheffield, Hounsfield Road,
Sheffield S3 7RH, United Kingdom  Celia Escamilla-Rivera
<EMAIL_ADDRESS>Instituto de Ciencias Nucleares, Universidad
Nacional Autónoma de México, Circuito Exterior C.U., A.P. 70-543, México D.F.
04510, México  Andronikos Paliathanasis<EMAIL_ADDRESS>Institute of
Systems Science, Durban University of Technology, PO Box 1334, Durban 4000,
Republic of South Africa  Departamento de Matemáticas, Universidad Católica
del Norte, Avda. Angamos 0610, 1280 Casilla, Antofagasta, Chile  Mathematical
Physics and Computational Statistics Research Laboratory, Department of
Environment, Ionian University, Zakinthos 29100, Greece

###### Abstract

We consider the unified bulk viscous scenarios and constrain them using the
Cosmic Microwave Background observations from Planck 2018 and the Pantheon
sample from Type Ia Supernovae. Then we generate the luminosity distance
measurements from ${\cal O}(10^{3})$ mock Gravitational Wave Standard Sirens
(GWSS) events for the proposed Einstein Telescope. We then combine these mock
luminosity distance measurements from the GWSS with the current cosmological
probes in order to forecast how the mock GWSS data could be effective in
constraining these bulk viscous scenarios. Our results show that a non-zero
time dependent bulk viscosity in the universe sector is strongly preferred by
the current cosmological probes and will possibly be confirmed at many
standard deviations by the future GWSS measurements. We further mention that
the addition of GWSS data can significantly reduce the uncertainties of the
key cosmological parameters obtained from the usual cosmological probes
employed in this work.

##  I Introduction

Understanding the nature of dark matter and dark energy has been a challenge
for cosmologists. The standard cosmological model, namely, the so-called
$\Lambda$-Cold Dark Matter ($\Lambda$CDM) model representing a mixture of two
non-interacting fluids $-$ a positive cosmological constant ($\Lambda>0$) and
a cold dark matter component, has undoubtedly proved its unprecedented success
by explaining a large span of astronomical data. However, this simplest
cosmological scenario has some limitations. For example, the cosmological
constant problem Weinberg (1989) and the coincidence problem Zlatev  _et al._
(1999) have already questioned the existing assumptions in the $\Lambda$CDM
model, e.g. constant energy density of the vacuum and the non-interacting
nature between $\Lambda$ and CDM. These limitations motivated the cosmologists
to find alternative cosmological scenarios beyond $\Lambda$CDM by relaxing the
above assumptions, and as a consequence, several new cosmological models were
introduced, see Copeland  _et al._ (2006); Nojiri and Odintsov (2006);
Sotiriou and Faraoni (2010); De Felice and Tsujikawa (2010); Capozziello and
De Laurentis (2011); Clifton  _et al._ (2012); Bamba  _et al._ (2012); Cai
_et al._ (2016); Nojiri  _et al._ (2017); Bahamonde  _et al._ (2021) for a
review of various dark energy and modified gravity models. Additionally, the
appearance of cosmological tensions at many standard deviations between Planck
Aghanim  _et al._ (2020a) (assuming $\Lambda$CDM in the background) and other
cosmological probes, such as distance ladders Wong  _et al._ (2020); Riess
_et al._ (2021); Shajib  _et al._ (2020); Birrer  _et al._ (2020); Freedman
(2021); Anand  _et al._ (2022); Blakeslee  _et al._ (2021); Pesce  _et al._
(2020); Shah  _et al._ (2021); de Jaeger  _et al._ (2022); Riess  _et al._
(2022) or weak lensing Heymans  _et al._ (2012); Hildebrandt  _et al._ (2017);
Asgari  _et al._ (2021); Abbott  _et al._ (2018, 2022) and galaxy cluster data
Ade  _et al._ (2016); Bocquet  _et al._ (2019); Hasselfield  _et al._ (2013)
has further weakened the confidence in the $\Lambda$CDM cosmological model
Riess (2019); Verde  _et al._ (2019); Di Valentino (2021); Di Valentino  _et
al._ (2021a, b). Thus, the list of cosmological models aiming to address the
cosmological tensions is increasing in time, see the review articles Knox and
Millea (2020); Jedamzik  _et al._ (2021); Di Valentino  _et al._ (2021c);
Perivolaropoulos and Skara (2022); Schöneberg  _et al._ (2022); Abdalla  _et
al._ (2022); Kamionkowski and Riess (2022) and references therein. Given the
fact that the origin of dark matter and dark energy is not clearly understood
yet, thus, there is no reason to favor any particular cosmological theory over
others. As a result, various ways have been proposed to interpret the dynamics
of the dark sector in terms of dark matter and dark energy. The simplest
assumption is the consideration of independent evolution of these dark fluids.
The generalization of the above consideration is the assumption of a non-
gravitational interaction between these dark sectors. On the other hand, a
heuristic approach is to consider a unified dark fluid that can explain the
dynamics of dark energy and dark matter at cosmological scales. The attempt to
unify the dark sector of the Universe began long back ago. The most simplest
unified dark sector models can be constructed in the context of Einstein
gravity with the introduction of a generalized equation of state
$p=\mathcal{F}(\rho)$, where $p$ and $\rho$ are respectively the pressure and
energy density of the unified dark sector and $\mathcal{F}$ is an analytic
function of the energy density, $\rho$. The well known unified cosmological
models, such as the Chaplygin gas model Chaplygin (1904) and its successive
generalizations, namely, the generalized Chaplygin gas, modified Chaplygin
gas, see Refs. Kamenshchik  _et al._ (2001); Bilic  _et al._ (2002); Gorini
_et al._ (2003, 2005); Bento  _et al._ (2002); Xu and Lu (2010); Lu  _et al._
(2009); Benaoum (2022); Debnath  _et al._ (2004); Lu  _et al._ (2008); Xu  _et
al._ (2012); Escamilla-Rivera  _et al._ (2020) and some other unified
cosmological scenarios as well Hova and Yang (2017); Hernandez-Almada  _et
al._ (2019); Yang  _et al._ (2019a) belong to this classification. While it is
essential to mention that a subset of the unified models has been diagnosed
with exponential blowup in the matter power spectrum which is not consistent
with the observations Sandvik  _et al._ (2004), however, this does not rule
out the possibility of unified models aiming to cover a wide region of the
universe evolution because a new kind of unified fluid may avoid such
unphysical activities. The unified cosmological models can also be developed
by considering a relation like $p=\mathcal{G}(H)$ where $\mathcal{G}$ is an
analytic function of $H$, the Hubble function of the Friedmann-Lemaître-
Robertson-Walker (FLRW) line element. Apparently, theories with
$p=\mathcal{F}(\rho)$ and $p=\mathcal{G}(H)$ seem identical, however, this is
only true in spatially flat FLRW universe. For a curved universe, the two
approaches are not the same.

In the present work we are interested to study a particular class of unified
models endowed with bulk viscosity. The cosmological fluids allowing bulk
viscosity as an extra ingredient can explain the accelerating expansion of the
universe, and hence they are also enlisted as possible alternatives to the
standard $\Lambda$CDM cosmology in the literature Yadav  _et al._ (2021);
Normann and Brevik (2021). Following an earlier work Ref. Yang  _et al._
(2019b) where an evidence of non-zero bulk viscosity was preferred by the
current cosmological probes, in the present article, we use the simulated
Gravitational Waves Standard Sirens (GWSS) measurements from the Einstein
Telescope Maggiore  _et al._ (2020)111https://www.einsteintelescope.nl/en/ in
order to quantify the improvements of the cosmological parameters, if any,
from the future GWSS measurements. As the gravitational waves (GW) have opened
a new window for astrophysics and cosmology, therefore, it will be interesting
to investigate the contribution from the simulated GWSS data, once combined
with the current cosmological probes. This motivated many investigators to use
the mock GWSS data matching the expected sensitivity of the Einstein Telescope
to constrain a class of cosmological models, see for instance, Zhao  _et al._
(2011); Cai and Yang (2017); Zhang  _et al._ (2017); Cai  _et al._ (2018);
Yang  _et al._ (2020a, 2019c); Bachega  _et al._ (2020); Zhou  _et al._
(2019); Yang  _et al._ (2020b); Mitra  _et al._ (2021); Escamilla-Rivera and
Quevedo (2021); Pan  _et al._ (2021). In particular, the combined analysis of
simulated GWSS measurements from Einstein Telescope and the standard
cosmological probes has proven to be very effective for a class of
cosmological models, in the sense that the error bars in the key cosmological
parameters of these cosmological models are significantly reduced thanks to
the mock GWSS dataset Zhang  _et al._ (2019); Wang  _et al._ (2018); Di
Valentino  _et al._ (2018); Du  _et al._ (2019); Yang  _et al._ (2020a, 2019c,
b), however, in some specific $f(R)$ theories of gravity, the generated mock
GWSS from the Einstein Telescope may not be very much helpful to give
stringent constraints on them during its first phase of running Matos  _et
al._ (2021). Thus, one may expect that the constraining power of the Einstein
Telescope may depend on the underlying cosmological model. Aside from the
future GWSS measurements from the Einstein Telescope, one can also use the
simulated GWSS measurements from other GW observatories, such as, Laser
Interferometer Space Antenna (LISA) Allahyari  _et al._ (2022); Corman  _et
al._ (2021); Belgacem  _et al._ (2019); Cai  _et al._ (2017) and DECi-heltz
Interferometer Gravitational wave Observatory (DECIGO) Hou  _et al._ (2021);
Zheng  _et al._ (2021), TianQin Luo  _et al._ (2016). In this article, we
focus only on the simulated GWSS data from Einstein Telescope to constrain the
bulk viscous unified scenario.

The paper has been organized as follows: in Sec. II we discuss the
gravitational equations for the bulk viscous scenario. Sec. III describes the
observational data that we have considered for the analysis in this work. Sec.
IV presents the observational constraints on the bulk viscous models, and
mainly we discuss how the inclusion of gravitational waves data from the
Einstein Telescope improves the constraints. Finally, in Sec. V we present the
conclusions.

##  II Revisiting the bulk viscous scenarios: Background and perturbations

As usual, we consider the homogeneous and isotropic space time described by
the Friedmann-Lemaître-Robertson-Walker (FLRW) line element

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

where $a(t)$ is the expansion scale factor and $k$ denotes the spatial
curvature of the universe. For $k=0,-1,+1$, we have three different geometries
of the universe, namely, spatially flat, open and closed, respectively. In
this paper we restrict ourselves to the spatially flat scenario where we
assume that (i) the gravitational sector is described by the Einstein’s
gravity, (ii) the matter sector of the universe consists of the relativistic
radiation, non-relativistic baryons and a unified bulk viscous fluid which
combines the effects of dark matter and dark energy, (iii) all the fluids are
non-interacting with each other. Within this framework, we can write down the
gravitational field equations as follows (in the units where $8\pi G=1$)

$\displaystyle H^{2}=\frac{1}{3}\rho_{\rm tot},$ (2) $\displaystyle
2\dot{H}+3H^{2}=-\;p_{\rm tot},$ (3)

where an overhead dot indicates the derivative with respect to the cosmic time
$t$; $H\equiv\dot{a}/a$ is the Hubble expansion rate; $(\rho_{\rm tot},p_{\rm
tot})=(\rho_{r}+\rho_{b}+\rho_{u},p_{r}+p_{b}+p_{u})$ are the total energy
density and total pressure of the cosmic components in which
$(\rho_{r},p_{r})$, $(\rho_{b},p_{b})$, $(\rho_{u},p_{u})$ are the energy
density and pressure of radiation, baryons and the unified fluid,
respectively. The conservation equation for each fluid follows the usual law
$\dot{\rho}_{i}+3H(1+w_{i})\rho_{i}=0$, where the subscript $i$ refers to
radiation $(i=r)$, baryons $(i=b)$ and the unified fluid $(i=u)$ and $w_{i}$
are the standard barotropic state parameters: $w_{r}=p_{r}/\rho_{r}=1/3$,
$w_{b}=p_{b}/\rho_{b}=0$ and $w_{u}=p_{u}/\rho_{u}=(\gamma-1)$, where $\gamma$
is a constant parameter. In general for different values of $\gamma$, say for
instance, $\gamma=0$, we realize a cosmological constant-like fluid endowed
with the bulk viscosity and similarly $\gamma=1$ results in a dust-like fluid
endowed with the bulk viscosity. As the nature of the fluid is not clearly
understood and as the observational data play an effective role to understand
this nature, thus, in order to be more transparent in this direction we
consider $\gamma$ lying in the interval $[-3,3]$ which includes both exotic
($p_{u}/\rho_{u}=(\gamma-1)<-1/3$) and non-exotic
($p_{u}/\rho_{u}=(\gamma-1)>-1/3$ ) fluids. As already mentioned, since the
unified fluid has a bulk viscosity, therefore, it enjoys an effective pressure
Barrow (1988): $p_{\rm eff}=p_{u}-u^{\nu}_{;\nu}\eta(\rho_{u})$, where
$u^{\mu}_{;\mu}$ is the expansion scalar of this fluid and $\eta(\rho_{u})>0$
is the coefficient of the bulk viscosity. Thus, in the FLRW background, the
effective pressure of the bulk viscous fluid reduces to

$\displaystyle p_{\rm eff}=p_{u}-3H\eta(\rho_{u}).$ (4)

Since there is no unique selection for the bulk viscous coefficient,
$\eta(\rho_{u})$, therefore, we consider a well known choice for it in which
the bulk viscous coefficient has a power law evolution of the form Barrow
(1986, 1988, 1990):

$\displaystyle\eta(\rho_{u})=\alpha\rho_{u}^{m},$ (5)

where $\alpha$ is a positive constant and $m$ is any real number. Notice that
for the case $m=0$ we recover the scenario with a constant bulk viscous
coefficient. Now, with the consideration of the bulk viscous coefficient in
(5), the effective pressure of the unified fluid can be expressed as

$\displaystyle p_{\rm eff}=(\gamma-1)\rho_{u}-\sqrt{3}\alpha\rho_{\rm
tot}^{1/2}\rho_{u}^{m},$ (6)

and consequently, one can define the effective equation of state of the
viscous dark fluid as

$\displaystyle w_{\rm eff}=\frac{p_{\rm
eff}}{\rho_{u}}=(\gamma-1)-\sqrt{3}\alpha\rho_{\rm tot}^{1/2}\rho_{u}^{m-1}.$
(7)

The adiabatic sound speed for the viscous fluid is given by

$c_{a,\rm eff}^{2}=\frac{p_{\rm eff}^{\prime}}{\rho_{u}^{\prime}}=w_{\rm
eff}+\frac{w_{\rm eff}^{\prime}}{3\mathcal{H}(1+w_{\rm eff})}.$ (8)

where the prime denotes the derivative with respect to the conformal time
$\tau$ and $\mathcal{H}$ is the conformal Hubble parameter, $\mathcal{H}=aH$.
Note that depending on the nature of $w_{\rm eff}$, $c_{a,\rm eff}^{2}$ could
be negative, and hence $c_{a,\rm eff}$ could be an imaginary quantity. This
may invite instabilities in the perturbations. Thus, in order to avoid this
possible unphysical situation, we consider the entropy perturbations (non-
adiabatic perturbations) in the unified dark fluid following the analysis of
generalized dark matter Hu (1998).

Now we focus on the evolution of the unified bulk viscous fluid at the level
of perturbations. In the entropy perturbation mode, the true pressure
perturbation comes from the effective pressure given by

$\displaystyle\delta p_{\rm eff}$ $\displaystyle=$ $\displaystyle\delta
p_{u}-\delta\eta(\nabla_{\sigma}u^{\sigma})-\eta(\delta\nabla_{\sigma}u^{\sigma})$
(9) $\displaystyle=$ $\displaystyle\delta
p_{u}-3H\delta\eta-\frac{\eta}{a}\left(\theta+\frac{h^{\prime}}{2}\right).$

The effective sound speed of viscous dark fluid for the bulk viscous
coefficient (5) can be defined as

$\displaystyle c_{s,\rm eff}^{2}\equiv\left(\frac{\delta p_{\rm
eff}}{\delta\rho_{u}}\right)_{rf}$ $\displaystyle=c_{s}^{2}-\sqrt{3}\alpha
m\rho_{\rm
tot}^{1/2}\rho_{u}^{m-1}-\frac{\alpha\rho_{u}^{m-1}}{a\delta_{u}}\left(\theta+\frac{h^{\prime}}{2}\right),$
(10)

where ’$|_{rf}$’ denotes the rest frame. Following the analysis in Hu (1998),
the sound speed in the rest frame is assumed to be zero, i.e. $c^{2}_{s}=0$.

The density perturbation and the velocity perturbation can also be written as
Hu (1998)

$\displaystyle\delta_{u}^{\prime}$ $\displaystyle=$ $\displaystyle-(1+w_{\rm
eff})(\theta_{u}+\frac{h^{\prime}}{2})+\frac{w_{\rm eff}^{\prime}}{1+w_{\rm
eff}}\delta_{u}$ (11) $\displaystyle-$ $\displaystyle 3\mathcal{H}(c_{s,\rm
eff}^{2}-c_{a,\rm eff}^{2})\left[\delta_{u}+3\mathcal{H}(1+w_{\rm
eff})\frac{\theta_{D}}{k^{2}}\right],~{}\quad$
$\displaystyle\theta_{u}^{\prime}$ $\displaystyle=$
$\displaystyle-\mathcal{H}(1-3c_{s,\rm eff}^{2})\theta_{u}+\frac{c_{s,\rm
eff}^{2}}{1+w_{\rm eff}}k^{2}\delta_{u},$ (12)

Thus, following the evolution at the background and perturbation level
prescribed above, one can now be able to understand the dynamics of the bulk
viscous fluid. In this work we consider two different bulk viscous scenarios
characterized as follows: the bulk viscous model 1 (labeled as BVF1) where we
consider $\gamma=1$, and the bulk viscous model 2 (labeled as BVF2) where we
keep $\gamma$ as a free parameter. The common parameters in both BVF1 and BVF2
are $\alpha$ and $m$.

##  III Standard Cosmological probes, simulated GWSS data, and methodology

In this section we describe the cosmological data sets employed to perform the
statistical analyses of the present bulk viscous scenarios.

  * •

Cosmic Microwave Background (CMB): We use the CMB data from the Planck 2018
data release. Precisely, we use the CMB temperature and polarization angular
power spectra plikTTTEEE+lowl+lowE Aghanim  _et al._ (2020b, c).

  * •

Pantheon sample from Type Ia Supernovae (SNIa) data: Type Ia Supernovae are
the first astronomical data that probed the accelerating expansion of the
universe and hence indicated the existence of an exotic fluid with negative
pressure (dark energy). Here we use the Pantheon compilation of the SNIa data
comprising 1048 data points spanned in the redshift interval $[0.01,2.3]$
Scolnic  _et al._ (2018).

  * •

Gravitational Waves Standard Sirens (GWSS): We take the mock Gravitational
Waves Standard Sirens (GWSS) data generated by matching the expected
sensitivity of Einstein Telescope in order to understand the constraining
power of the future GWSS data from the Einstein Telescope. The Einstein
Telescope is a proposed ground based third-generation (3G) GW detector. The
telescope will take a triangular shape and its each arm length will be
increased to 10 km, compared to 3 km arm length VIRGO and 4 km arm length LIGO
Sathyaprakash  _et al._ (2012); Maggiore  _et al._ (2020). Thus, due to such
increased arm length, the Einstein Telescope will be a potential GW detector
by reducing all displacement noises Sathyaprakash  _et al._ (2012); Maggiore
_et al._ (2020). It is expected that after 10 years of operation, Einstein
Telescope will detect $\mathcal{O}(10^{3})$ GWSS events. Although the
detection of $\mathcal{O}(10^{3})$ GWSS events is very optimistic while the
number of detections could be low in reality Maggiore  _et al._ (2020). As
argued in Maggiore  _et al._ (2020), the Einstein Telescope will likely to
detect 20 $-$ 50 events per year, i.e. 200 $-$ 500 events in 10 years.
However, following the earlier works Zhao  _et al._ (2011); Cai and Yang
(2017); Wang  _et al._ (2018); Du  _et al._ (2019); Yang  _et al._ (2019c,
2020a, 2020b); Matos  _et al._ (2021); Pan  _et al._ (2021), in this article,
we restrict ourselves to the detection of $\mathcal{O}(10^{3})$ GWSS events by
the Einstein Telescope to constrain the bulk viscous scenarios. For more
features of the Einstein Telescope we refer the readers to Maggiore  _et al._
(2020).

We originally generate the mock GWSS luminosity distance measurements matching
the expected sensitivity of the Einstein Telescope after 10 years of full
operation. Specifically we create 1000 triples ($z_{i}$, $d_{L}(z_{i})$,
$\sigma_{i}$) where $z_{i}$ is the redshift of a GW source, $d_{L}(z_{i})$ is
the measured luminosity distance at redshift $z_{i}$ and $\sigma_{i}$ is the
uncertainty associated with the luminosity distance $d_{L}(z_{i})$. Let us
briefly summarize the procedure of generating the mock GWSS dataset and we
refer to Refs. Du  _et al._ (2019); Yang  _et al._ (2020a, 2019c) for more
technical details. The initial step for generating the mock GWSS dataset is to
identify the expected GW sources. We consider the GW events originating from
two distinct binary systems, namely, (i) a combination of a Black Hole (BH)
and a Neutron Star (NS) merger, identified as BHNS and (ii) the binary neutron
star (BNS) merger. Following the mass distributions as described in Ref. Du
_et al._ (2019), the ratio of the number of GW events for the BHNS merger
versus BNS merger is taken to be $0.03$ as predicted for the Advanced LIGO-
VIRGO network Abadie  _et al._ (2010). We then determine the merger rate
$R(z)$ of sources and from the merger rate of the sources, we determine the
redshift distribution of the sources, $P(z)$ given by Sathyaprakash  _et al._
(2010); Zhao  _et al._ (2011); Cai and Yang (2017); Wang  _et al._ (2018); Du
_et al._ (2019); Yang  _et al._ (2020a)

$\displaystyle P(z)\propto\frac{4\pi d_{C}^{2}(z)R(z)}{H(z)(1+z)},$ (13)

where $d_{C}(z)\equiv\int_{0}^{z}H^{-1}(z^{\prime})dz^{\prime}$ is the co-
moving distance and for $R(z)$ we take the following piece-wise linear
function estimated in Schneider  _et al._ (2001) (also see Cutler and Holz
(2009); Zhao  _et al._ (2011); Cai and Yang (2017); Wang  _et al._ (2018); Du
_et al._ (2019); Yang  _et al._ (2020a)): $R(z)=1+2z$ for $z\leq 1$,
$R(z)=\frac{3}{4}(5-z)$, for $1<z<5$ and $R(z)=0$ for $z>5$. After having
$P(z)$, we sample 1000 values of redshifts from this distribution which
represent the redshifts $z_{i}$ of our 1000 mock GWSS data.

The next step is to choose a fiducial model because while going from the
merger rate to the redshift distributions, a fiducial cosmological model is
needed since the expression for $P(z)$ includes both the comoving distance and
expansion rate at redshift $z$, i.e. $d_{L}(z)$ and $H(z)$ respectively. This
$H(z)$ corresponds to the fiducial model. As in this article we are interested
to investigate how the inclusion of GWSS data improves the constraints on the
BVF models, therefore, we generate two different mock GWSS datasets choosing
BVF1 and BVF2 as the fiducial models. We take the fiducial values of the
cosmological parameters given by the best fit values of the same cosmological
parameters of the BVF1 and BVF2 models obtained from the CMB+Pantheon data
analysis. Now, for the chosen fiducial model(s), one can now estimate the
luminosity distance at the redshift $z_{i}$ using the relation

$\displaystyle
d_{L}(z_{i})=(1+z_{i})\int_{0}^{z_{i}}\frac{dz^{\prime}}{H(z^{\prime})}\,.$
(14)

Thus, after having the luminosity distance $d_{L}(z_{i})$ of the GW source,
our last job is now to determine the uncertainty $\sigma_{i}$ associated with
this luminosity distance. The determination of the uncertainty $\sigma_{i}$
directly connects to the waveform of GW because the GW amplitude depends on
the luminosity distance (also on the so-called chirp mass Zhao  _et al._
(2011); Cai and Yang (2017); Wang  _et al._ (2018)) and hence one can extract
the information about $d_{L}(z_{i})$ and $\sigma_{i}$. We refer to Refs. Zhao
_et al._ (2011); Cai and Yang (2017); Wang  _et al._ (2018); Du  _et al._
(2019); Yang  _et al._ (2020a) for the technical details to calculate the
uncertainties on the luminosity distance measurements. The luminosity distance
measurement $d_{L}(z_{i})$ has two kind of uncertainties, one is the
instrumental uncertainty $\sigma_{i}^{\rm inst}$ and the other one is the weak
lensing uncertainty $\sigma_{i}^{\rm lens}$. The instrumental error can be
derived to be $\sigma_{i}^{\rm inst}~{}(\simeq 2d_{L}(z_{i})/\mathcal{S}$
where $\mathcal{S}$ is the combined signal-to-noise ratio of the Einstein
Telescope) using the Fisher matrix approach and assuming that the uncertainty
on $d_{L}(z_{i})$ is not correlated with the uncertainties on the remaining GW
parameters (see Zhao  _et al._ (2011); Cai and Yang (2017); Wang  _et al._
(2018); Du  _et al._ (2019); Yang  _et al._ (2020a)) and the lensing error is
$\sigma_{i}^{\rm lens}\simeq 0.05z_{i}d_{L}(z_{i})$ Zhao  _et al._ (2011).
Thus, the total uncertainty due to the instrumental and the weak lensing
uncertainties on $d_{L}(z_{i})$ is $\sigma_{i}=\sqrt{(\sigma_{i}^{\rm
inst})^{2}+(\sigma_{i}^{\rm lens})^{2}}$. Finally, let us note that the
combined signal-to-noise ratio of the GW detector is a very crucial quantity
in this context since for the Einstein Telescope, the combined signal-to-noise
ratio should be at least $8$ for a GW detection Sathyaprakash  _et al._
(2010). Thus, in summary, we generate 1000 GW sources up to redshift $z=2$
with $\mathcal{S}>8$. For more technical details we refer the readers to Refs.
Sathyaprakash  _et al._ (2010); Zhao  _et al._ (2011); Cai and Yang (2017);
Wang  _et al._ (2018); Du  _et al._ (2019); Yang  _et al._ (2020a).

To constrain the BVF scenarios we modify the publicly available CosmoMC
package Lewis and Bridle (2002) which is an excellent cosmological code
supporting the Planck 2018 likelihood Aghanim  _et al._ (2020c) and it has a
convergence diagnostic following the Gelman-Rubin statistic $R-1$ Gelman and
Rubin (1992). It is essential to mention that for both BVF1 and BVF2
scenarios, we have used the dimensionless quantity $\beta=\alpha
H_{0}\rho_{\rm tot,0}^{m-1}$ where $\rho_{\rm tot,0}$ is the present value of
$\rho_{\rm tot}$. We further mention here that $\beta=0$ (equivalently,
$\alpha=0$) implies no viscosity and hence the overall picture behaves like a
unified cosmic fluid without bulk viscosity. Thus, in summary, the parameter
space of BVF1 and BVF2 are as below:

$\displaystyle\mathcal{P}_{\rm BVF1}\equiv\\{\Omega_{b}h^{2},100\theta_{\rm
MC},\tau,n_{s},{\rm{ln}}(10^{10}A_{s}),\beta,m\\}$
$\displaystyle\mathcal{P}_{\rm BVF2}=\\{\Omega_{b}h^{2},100\theta_{\rm
MC},\tau,n_{s},{\rm{ln}}(10^{10}A_{s}),\beta,m,\gamma\\}$

where the description of the free parameters are as follows: $\Omega_{b}h^{2}$
is the baryons density, $100\theta_{MC}$ is the ratio of the sound horizon to
the angular diameter distance; $\tau$ is the optical depth, $n_{s}$ is the
scalar spectral index, $A_{s}$ is the amplitude of the initial power spectrum.
The flat priors on both cosmological scenarios are shown in Table 1.

Parameter | Priors (BVF1) | Priors (BVF2)
---|---|---
$\Omega_{b}h^{2}$ | $[0.005,0.1]$ | $[0.005,0.1]$
$\tau$ | $[0.01,0.8]$ | $[0.01,0.8]$
$n_{s}$ | $[0.5,1.5]$ | $[0.5,1.5]$
${\rm{ln}}(10^{10}A_{s})$ | $[2.4,4]$ | $[2.4,4]$
$100\theta_{\rm MC}$ | $[0.5,10]$ | $[0.5,10]$
$\beta$ | $[0,1]$ | $[0,1]$
$m$ | $[-2,0.5]$ | $[-2,0.5]$
$\gamma$ | $-$ | $[-3,3]$
Table 1: We show the flat priors on all the free parameters associated with the bulk viscous models.  Parameters | CMB+Pantheon | CMB+Pantheon+GWSS
---|---|---
$\Omega_{b}h^{2}$ | $0.02232_{-0.00015-0.00028}^{+0.00015+0.00029}$ | $0.02253_{-0.00014-0.00026}^{+0.00014+0.00028}$
$100\theta_{\rm MC}$ | $1.02780_{-0.00055-0.0011}^{+0.00058+0.0011}$ | $1.02808_{-0.00038-0.00073}^{+0.00037+0.00073}$
$\tau$ | $0.0537_{-0.0075-0.015}^{+0.0074+0.016}$ | $0.0567_{-0.0078-0.015}^{+0.0079+0.016}$
$n_{s}$ | $0.9641_{-0.0043-0.0084}^{+0.0043+0.0086}$ | $0.9686_{-0.0040-0.0080}^{+0.0041+0.0080}$
${\rm{ln}}(10^{10}A_{s})$ | $3.046_{-0.015-0.031}^{+0.016+0.031}$ | $3.048_{-0.016-0.033}^{+0.016+0.033}$
$\beta$ | $0.430_{-0.016-0.034}^{+0.017+0.033}$ | $0.4262_{-0.0078-0.015}^{+0.0079+0.016}$
$m$ | $-0.557_{-0.059-0.13}^{+0.068+0.12}$ | $-0.519_{-0.035-0.075}^{+0.038+0.074}$
$H_{0}$ | $68.1_{-1.1-2.3}^{+1.2+2.2}$ | $68.30_{-0.45-0.85}^{+0.46+0.91}$
Table 2: We report the observational constraints on the BVF1 scenario at 68%
and 95% CL for CMB+Pantheon and CMB+Pantheon+GWSS datasets.  Figure 1: For the
BVF1 scenario we show the 1-dimensional posterior distribution of some model
parameters and the 2-dimensional joint contours of the model parameters at 68%
and 95% CL for CMB+Pantheon and CMB+Pantheon+GWSS.

Figure 2:  The CMB $C^{TT}_{l}$ power spectrum versus multipole moment $l$
using the best fits values obtained for the BVF1 model using the join data
sets described, with three arbitrary $\beta$ and $m$ values.

##  IV Observational constraints: Results and analysis

In this section we present the constraints on the bulk viscous scenarios
considering CMB+Pantheon and CMB+Pantheon+GWSS. As we are interested to
estimate the improvement of the cosmological parameters in presence of the
GWSS measurements, and as the combined standard cosmological probes offer the
most stringent constraints on the cosmological parameters, therefore, the
inclusion of GWSS with the combined standard cosmological probes is
reasonable. As mentioned earlier, the key common free parameters of BVF1 and
BVF2 are $\beta$ and $m$, since $\beta\neq 0$ indicates the preference for a
non-zero bulk viscosity and $m\neq 0$ indicates that the coefficient of the
bulk viscosity is not constant in the redshift range considered. In the
following subsections we discuss the constraints on these two scenarios in
detail.

###  IV.1 Constraints on the BVF1 scenario

In Table 2 we have presented the constraints on the BVF1 scenario for
CMB+Pantheon and CMB+Pantheon+GWSS. The latter dataset is aimed to understand
the improvement expected from GWSS on the constraints from CMB+Pantheon. In
Fig. 1 we have compared these datasets graphically by showing the one
dimensional marginalized distribution of some model parameters and the two
dimensional joint contours. As discussed, this scenario has two main key
parameters, namely, $\beta$, quantifying the existence of bulk viscosity in
the cosmic sector, and $m$, which tells us whether the bulk viscosity will
have a dynamical nature (corresponding to $m\neq 0$) or not.

Since for CMB+Pantheon, we find an evidence for a non-zero bulk viscosity in
the cosmic sector at many standard deviations, i.e.
$\beta=0.430_{-0.016}^{+0.017}$ at 68% CL, this is further strengthen for
CMB+Pantheon+GWSS, where $\beta=0.4262_{-0.0078}^{+0.0079}$ at 68% CL.222It is
worthwhile to note here that the mean value of $\beta$ is not significantly
changed when the GWSS data are included, because we built the simulated data
using the best fit obtained from CMB+Pantheon.  One can clearly see that the
inclusion of GWSS to CMB+Pantheon improves the error bars on $\beta$ by a
factor of at least 2. This reflects the constraining power of GWSS. On the
other hand, focusing on the parameter $m$, which quantifies the time evolution
of the bulk viscosity, we see that $m$ remains non-zero at several standard
deviations for CMB+Pantheon, where the 68% CL constraint on $m$ is
$m=-0.557_{-0.059}^{+0.068}$, and becomes $m=-0.519_{-0.035}^{+0.038}$ for
CMB+Pantheon+GWSS. From the constraints on $m$, one can clearly see that the
uncertainty in $m$ is reduced by a factor of $\sim 1.7-1.8$ when the GWSS data
are included with the combined dataset CMB+Pantheon. Concerning the Hubble
constant, we find that $H_{0}$ assumes slightly higher values compared to the
$\Lambda$CDM based Planck. Actually, we have $H_{0}=68.1_{-1.1}^{+1.2}$ at 68%
CL for CMB+Pantheon, while $H_{0}=68.30_{-0.45}^{+0.46}$ at 68% CL for
CMB+Pantheon+GWSS, again improving the uncertainty in $H_{0}$ by a factor of
$2.5$. This shows that the effects of GWSS are clearly visible through these
parameters. In Fig. 1, one can compare the constraints on the model parameters
obtained from CMB+Pantheon and CMB+Pantheon+GWSS.

Finally, through Fig. 2 we examine how the model affects the CMB TT power
spectrum for different values of the model parameters, $\beta$ and $m$ with
respect to the standard $\Lambda$CDM scenario. In the upper panel of Fig. 2 we
depict the evolution in the CMB TT power spectrum for different values of
$\beta$ while in the lower panel of Fig. 2 we depict the evolution in the CMB
TT power spectrum for different values of $m$. From both the graphs, we notice
that as long as $\beta$ or $m$ increases, the model exhibits significant
differences in the lower multipoles ($\ell\leq 10$). For $\ell\geq 10$, we
observe that with the increasing values of $\beta$ and $m$, the peaks of the
CMB TT power spectrum increase significantly, particularly changing their
mutual ratio.

Parameters | CMB+Pantheon | CMB+Pantheon+GWSS
---|---|---
$\Omega_{b}h^{2}$ | $0.02241_{-0.00016-0.00032}^{+0.00016+0.00032}$ | $0.02238_{-0.00016-0.00031}^{+0.00015+0.00030}$
$100\theta_{\rm MC}$ | $1.02907_{-0.00082-0.00198}^{+0.00111+0.00180}$ | $1.02921_{-0.00040-0.00079}^{+0.00041+0.00080}$
$\tau$ | $0.0516_{-0.0072-0.015}^{+0.0074+0.015}$ | $0.0521_{-0.0078-0.014}^{+0.0071+0.015}$
$n_{s}$ | $0.9575_{-0.0066-0.012}^{+0.0053+0.012}$ | $0.9583_{-0.0039-0.0077}^{+0.0038+0.0075}$
${\rm{ln}}(10^{10}A_{s})$ | $3.038_{-0.017-0.032}^{+0.016+0.032}$ | $3.040_{-0.015-0.031}^{+0.015+0.032}$
$\beta$ | $0.447_{-0.022-0.042}^{+0.022+0.042}$ | $0.425_{-0.016-0.034}^{+0.018+0.032}$
$m$ | $-0.85_{-0.19-0.50}^{+0.30+0.46}$ | $-0.683_{-0.089-0.19}^{+0.099+0.18}$
$\gamma$ | $0.9970_{-0.0024-0.0036}^{+0.0015+0.0042}$ | $0.99757_{-0.00058-0.0011}^{+0.00049+0.0011}$
$H_{0}$ | $65.2_{-2.6-3.9}^{+1.7+4.4}$ | $64.91_{-0.60-1.2}^{+0.59+1.1}$
Table 3: We report the observational constraints on the BVF2 scenario at 68%
and 95% CL for CMB+Pantheon and CMB+Pantheon+GWSS. Figure 3: For the BVF2
scenario we show the 1-dimensional posterior distributions of some model
parameters and the 2-dimensional joint contours of the model parameters at 68%
and 95% CL for CMB+Pantheon and CMB+Pantheon+GWSS.

Figure 4: The CMB $C^{TT}_{l}$ power spectrum versus multipole moment $l$
using the best fits values obtained for each BVF2 models with three $\gamma$
values, respectively using the join data sets described.

###  IV.2 Constraints on the BVF2 scenario

In Table 3 we present the constraints on the BVF2 scenario for both
CMB+Pantheon and CMB+Pantheon+GWSS. And in Fig. 3, we compare the constraints
from these datasets explicitly showing the one dimensional marginalized
distribution of some model parameters and the two dimensional joint contours.
As already discussed, this scenario has three main key parameters, namely,
$\beta$, which quantifies the existence of bulk viscosity in the cosmic
sector, $m$, which tells us whether the bulk viscosity enjoys a dynamical
nature (corresponding to $m\neq 0$) or not, and finally, the parameter
$\gamma$, which indicates the fluid which endows the bulk viscosity. We note
that $\gamma=1$ refers to the dust fluid endowing the bulk viscosity in which
we are interested in, for which we recover the previous scenario BVF1.

For CMB+Pantheon, we find that $\beta\neq 0$ at several standard deviations
yielding $\beta=0.447\pm 0.022$ at 68% CL which gives a clear indication of a
non-zero bulk viscosity in the cosmic sector. When the GWSS are added to this
combination, i.e. CMB+Pantheon+GWSS, the conclusion about $\beta$ does not
change significantly ($\beta=0.425_{-0.016}^{+0.018}$ at 68% CL), indicating
that for this scenario GWSS do not provide any additional constraining power
on $\beta$. Looking at the dynamical nature of the bulk viscosity, we see that
for CMB+Pantheon, $m$ remains nonzero at more than 2 standard deviations
leading to $m=-0.85_{-0.19}^{+0.30}$ at 68% CL. However, this evidence could
be further strengthened by the inclusion of the GWSS data, that we forecast to
be $m=-0.683_{-0.089}^{+0.099}$ at 68% CL for CMB+Pantheon+GWSS, improving the
error bars up to a factor of 3. Finally, focusing on the parameter $\gamma$
which directly connects with the nature of the cosmic fluid endowing the bulk
viscosity, we can see that it is consistent with $1$, which corresponds to a
dust-like fluid, within 2 standard deviations for CMB+Pantheon
($\gamma=0.9970^{+0.0015}_{-0.0024}$ at 68% CL). Also for this parameter, the
addition of the GWSS further improves the constraining power of a factor
larger than $3$ to $4$, that we forecast to be
$\gamma=0.99757^{+0.00049}_{-0.00058}$ at 68% CL for CMB+Pantheon+GWSS.
Therefore, with respect to the BVF1 case, where the inclusion of the
forecasted GWSS was able to improve both $\beta$ and $m$, in this BVF2
scenario, the improvement of the constraining power is displayed only on $m$
and $\gamma$ but does not affect anymore $\beta$ significantly.

Furthermore, we find that for this scenario, the Hubble constant attains a
very low value for CMB+Pantheon compared to Planck’s estimation within the
$\Lambda$CDM paradigm. We also note that $H_{0}$ is correlated to all three
free parameters of this scenario, namely, $\beta$, $m$ and $\gamma$. With
$\beta$ and $\gamma$, $H_{0}$ is positively correlated while with $m$, it has
a strong anti-correlation. For CMB+Pantheon, $H_{0}=65.2_{-2.6}^{+1.7}$
km/s/Mpc at 68% CL and after the inclusion of GWSS it becomes
$H_{0}=64.91_{-0.60}^{+0.59}$ km/s/Mpc at 68% CL, reducing the uncertainty in
$H_{0}$ by a factor of $\sim 4$.

Finally, in Fig. 4, we examine the CMB TT power spectrum for this bulk viscous
scenario BVF2 considering different values of the free parameter $\gamma$ with
respect to the standard $\Lambda$CDM scenario. As $\gamma$ lies in the region
$[-3,3]$ and the nature of the cosmic fluid characterized by its equation of
state $p_{u}=(\gamma-1)\rho_{u}$ depends on the sign of $\gamma$, therefore,
we have considered two separate plots, one where $\gamma$ is non-negative
(i.e. $\gamma\geq 0$) and another plot where $\gamma$ allows both positive and
negative values including $\gamma=0$. From both the panels of Fig. 4, we
clearly see that any deviation from $\gamma=1$ makes significant changes in
the amplitude of the CMB TT power spectrum. In particular, we see that the
peaks of the CMB TT spectrum significantly increases and shift towards higher
multipoles for any value different from $\gamma=1$ at small scales, as well as
the Integrated Sachs Wolfe (ISW) plateau at large scales. As $\gamma=1$
indicates a cosmological constant-like fluid endowed with the bulk viscosity,
therefore, for $\gamma=1$, we replicate an equivalent behaviour of the
$\Lambda$CDM scenario.

##  V Conclusions

Although the $\Lambda$CDM cosmological model is extremely successful in
describing a large span of astronomical observations, it cannot explain
several theoretical and observational issues. This motivated the scientific
community to construct a variety of cosmological proposals and testing them
with the available astronomical data. Among these cosmological models, in this
article, we focus on the unified cosmological models allowing bulk viscosity
in the background. However, since these models do not recover the $\Lambda$CDM
scenario as a special case, our only ability in distinguishing them, once the
GWSS data will be available, will rely only on a Bayesian model comparison for
a better fit of the cosmological observations, as done in Ref. Yang  _et al._
(2019b). The unified cosmological scenarios endowed with bulk viscosity are
appealing from two different perspectives: the first one is the concept of a
unified picture of dark matter and dark energy, and the second is the
inclusion of bulk viscosity into that unified picture. Effectively, the
unified bulk viscous scenario is a generalized cosmic picture combining two
distinct cosmological directions. According to a recent paper on the unified
bulk viscous scenarios Yang  _et al._ (2019b), current cosmological probes
prefer a non-zero dynamical bulk viscosity in the dark sector at many standard
deviations. So, in light of the current cosmological probes, unified bulk
viscous cosmological scenarios are attractive. In this line of thought, what
about the future of such unified bulk viscous scenarios? In this article we
have focused on it with an answer.

Following Ref. Yang  _et al._ (2019b), in this work we have explored these
scenarios with the GWSS aiming to understand how the future distance
measurements from GWSS may improve the constraints on such scenarios. In order
to proceed toward this confrontation, we have generated ${\cal O}(10^{3})$
mock GWSS luminosity distance measurements matching the expected sensitivity
of the Einstein Telescope and added these mock data to the current
cosmological probes, namely CMB from Planck 2018 release333We mention that in
the earlier work Yang  _et al._ (2019b), CMB data from Planck 2015 were used
to constrain the bulk viscous scenarios.  and SNIa Pantheon sample. We find
that the inclusion of GWSS luminosity distance measurements together with the
current cosmological probes makes the possible future evidence for new physics
stronger, by reducing the uncertainty in the parameters in a significant way.
This is a potential behaviour of the GWSS luminosity distance measurements
since this makes the parameter much deterministic. Overall for both BVF1 and
BVF2 scenarios, we find a very strong preference of a non-zero time dependent
bulk viscous coefficient (alternatively, the viscous nature of the unified
dark fluid) at many standard deviations.

In conclusion, in the present paper we demonstrate that future GWSS distance
measurements from the Einstein Telescope might be powerful to extract more
information about the physics of the dark sector. Therefore, based on the
present results, we feel that it might just be a matter of time before we
convincingly detect the viscosity in the dark sector, if any.

##  VI Acknowledgments

The authors thank the referee for some important comments which helped us to
improve the article considerably. WY was supported by the National Natural
Science Foundation of China under Grants No. 12175096 and No. 11705079, and
Liaoning Revitalization Talents Program under Grant no. XLYC1907098. SP
acknowledges the financial support from the Department of Science and
Technology (DST), Govt. of India, under the Scheme “Fund for Improvement of
S&T Infrastructure (FIST)” [File No. SR/FST/MS-I/2019/41]. EDV is supported by
a Royal Society Dorothy Hodgkin Research Fellowship. CE-R is supported by the
Royal Astronomical Society as FRAS 10147 and by PAPIIT UNAM Project TA100122.
This article/publication is based upon work from COST Action CA21136
Addressing observational tensions in cosmology with systematics and
fundamental physics (CosmoVerse) supported by COST (European Cooperation in
Science and Technology). AP was supported in part by the National Research
Foundation of South Africa (Grant Number 131604). Also AP thanks the support
of Vicerrectoría de Investigación y Desarrollo Tecnológico (Vridt) at
Universidad Católica del Norte through Núcleo de Investigación Geometría
Diferencial y Aplicaciones, Resolución Vridt No - 096/2022.

## References

  * Weinberg (1989) S. Weinberg, Rev. Mod. Phys. 61, 1 (1989).
  * Zlatev  _et al._ (1999) I. Zlatev, L.-M. Wang, and P. J. Steinhardt, Phys. Rev. Lett. 82, 896 (1999), arXiv:astro-ph/9807002 .
  * Copeland  _et al._ (2006) E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006), arXiv:hep-th/0603057 .
  * Nojiri and Odintsov (2006) S. Nojiri and S. D. Odintsov, eConf C0602061, 06 (2006), arXiv:hep-th/0601213 .
  * Sotiriou and Faraoni (2010) T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82, 451 (2010), arXiv:0805.1726 [gr-qc] .
  * De Felice and Tsujikawa (2010) A. De Felice and S. Tsujikawa, Living Rev. Rel. 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
  * Capozziello and De Laurentis (2011) S. Capozziello and M. De Laurentis, Phys. Rept. 509, 167 (2011), arXiv:1108.6266 [gr-qc] .
  * Clifton  _et al._ (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rept. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Bamba  _et al._ (2012) K. Bamba, S. Capozziello, S. Nojiri, and S. D. Odintsov, Astrophys. Space Sci. 342, 155 (2012), arXiv:1205.3421 [gr-qc] .
  * Cai  _et al._ (2016) Y.-F. Cai, S. Capozziello, M. De Laurentis, and E. N. Saridakis, Rept. Prog. Phys. 79, 106901 (2016), arXiv:1511.07586 [gr-qc] .
  * Nojiri  _et al._ (2017) S. Nojiri, S. Odintsov, and V. Oikonomou, Phys. Rept. 692, 1 (2017), arXiv:1705.11098 [gr-qc] .
  * Bahamonde  _et al._ (2021) S. Bahamonde, K. F. Dialektopoulos, C. Escamilla-Rivera, G. Farrugia, V. Gakis, M. Hendry, M. Hohmann, J. L. Said, J. Mifsud, and E. Di Valentino, arXiv:2106.13793 (2021), arXiv:2106.13793 [gr-qc] .
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020a), arXiv:1807.06209 [astro-ph.CO] .
  * Wong  _et al._ (2020) K. C. Wong _et al._ , Mon. Not. Roy. Astron. Soc. 498, 1420 (2020), arXiv:1907.04869 [astro-ph.CO] .
  * Riess  _et al._ (2021) A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn, and D. Scolnic, Astrophys. J. Lett. 908, L6 (2021), arXiv:2012.08534 [astro-ph.CO] .
  * Shajib  _et al._ (2020) A. J. Shajib _et al._ (DES), Mon. Not. Roy. Astron. Soc. 494, 6072 (2020), arXiv:1910.06306 [astro-ph.CO] .
  * Birrer  _et al._ (2020) S. Birrer _et al._ , Astron. Astrophys. 643, A165 (2020), arXiv:2007.02941 [astro-ph.CO] .
  * Freedman (2021) W. L. Freedman, Astrophys. J. 919, 16 (2021), arXiv:2106.15656 [astro-ph.CO] .
  * Anand  _et al._ (2022) G. S. Anand, R. B. Tully, L. Rizzi, A. G. Riess, and W. Yuan, Astrophys. J. 932, 15 (2022), arXiv:2108.00007 [astro-ph.CO] .
  * Blakeslee  _et al._ (2021) J. P. Blakeslee, J. B. Jensen, C.-P. Ma, P. A. Milne, and J. E. Greene, Astrophys. J. 911, 65 (2021), arXiv:2101.02221 [astro-ph.CO] .
  * Pesce  _et al._ (2020) D. W. Pesce _et al._ , Astrophys. J. Lett. 891, L1 (2020), arXiv:2001.09213 [astro-ph.CO] .
  * Shah  _et al._ (2021) P. Shah, P. Lemos, and O. Lahav, Astron. Astrophys. Rev. 29, 9 (2021), arXiv:2109.01161 [astro-ph.CO] .
  * de Jaeger  _et al._ (2022) T. de Jaeger, L. Galbany, A. G. Riess, B. E. Stahl, B. J. Shappee, A. V. Filippenko, and W. Zheng, Mon. Not. Roy. Astron. Soc. 514, 4620 (2022), arXiv:2203.08974 [astro-ph.CO] .
  * Riess  _et al._ (2022) A. G. Riess _et al._ , Astrophys. J. Lett. 934, L7 (2022), arXiv:2112.04510 [astro-ph.CO] .
  * Heymans  _et al._ (2012) C. Heymans _et al._ , Mon. Not. Roy. Astron. Soc. 427, 146 (2012), arXiv:1210.0032 [astro-ph.CO] .
  * Hildebrandt  _et al._ (2017) H. Hildebrandt _et al._ , Mon. Not. Roy. Astron. Soc. 465, 1454 (2017), arXiv:1606.05338 [astro-ph.CO] .
  * Asgari  _et al._ (2021) M. Asgari _et al._ (KiDS), Astron. Astrophys. 645, A104 (2021), arXiv:2007.15633 [astro-ph.CO] .
  * Abbott  _et al._ (2018) T. M. C. Abbott _et al._ (DES), Phys. Rev. D 98, 043526 (2018), arXiv:1708.01530 [astro-ph.CO] .
  * Abbott  _et al._ (2022) T. M. C. Abbott _et al._ (DES), Phys. Rev. D 105, 023520 (2022), arXiv:2105.13549 [astro-ph.CO] .
  * Ade  _et al._ (2016) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A24 (2016), arXiv:1502.01597 [astro-ph.CO] .
  * Bocquet  _et al._ (2019) S. Bocquet _et al._ (SPT), Astrophys. J. 878, 55 (2019), arXiv:1812.01679 [astro-ph.CO] .
  * Hasselfield  _et al._ (2013) M. Hasselfield _et al._ , JCAP 07, 008 (2013), arXiv:1301.0816 [astro-ph.CO] .
  * Riess (2019) A. G. Riess, Nature Rev. Phys. 2, 10 (2019), arXiv:2001.03624 [astro-ph.CO] .
  * Verde  _et al._ (2019) L. Verde, T. Treu, and A. G. Riess, in  _Nature Astronomy 2019_ (2019) arXiv:1907.10625 [astro-ph.CO] .
  * Di Valentino (2021) E. Di Valentino, Mon. Not. Roy. Astron. Soc. 502, 2065 (2021), arXiv:2011.00246 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021a) E. Di Valentino _et al._ , Astropart. Phys. 131, 102604 (2021a), arXiv:2008.11285 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021b) E. Di Valentino _et al._ , Astropart. Phys. 131, 102605 (2021b), arXiv:2008.11284 [astro-ph.CO] .
  * Knox and Millea (2020) L. Knox and M. Millea, Phys. Rev. D 101, 043533 (2020), arXiv:1908.03663 [astro-ph.CO] .
  * Jedamzik  _et al._ (2021) K. Jedamzik, L. Pogosian, and G.-B. Zhao, Commun. in Phys. 4, 123 (2021), arXiv:2010.04158 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021c) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Class. Quant. Grav. 38, 153001 (2021c), arXiv:2103.01183 [astro-ph.CO] .
  * Perivolaropoulos and Skara (2022) L. Perivolaropoulos and F. Skara, New Astron. Rev. 95, 101659 (2022), arXiv:2105.05208 [astro-ph.CO] .
  * Schöneberg  _et al._ (2022) N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, Phys. Rept. 984, 1 (2022), arXiv:2107.10291 [astro-ph.CO] .
  * Abdalla  _et al._ (2022) E. Abdalla _et al._ , JHEAp 34, 49 (2022), arXiv:2203.06142 [astro-ph.CO] .
  * Kamionkowski and Riess (2022) M. Kamionkowski and A. G. Riess, arXiv:2211.04492 (2022), arXiv:2211.04492 [astro-ph.CO] .
  * Chaplygin (1904) S. Chaplygin, Sci. Mem. Moscow Univ. Math. Phys. 21, 1 (1904).
  * Kamenshchik  _et al._ (2001) A. Y. Kamenshchik, U. Moschella, and V. Pasquier, Phys. Lett. B 511, 265 (2001), arXiv:gr-qc/0103004 .
  * Bilic  _et al._ (2002) N. Bilic, G. B. Tupper, and R. D. Viollier, Phys. Lett. B 535, 17 (2002), arXiv:astro-ph/0111325 .
  * Gorini  _et al._ (2003) V. Gorini, A. Kamenshchik, and U. Moschella, Phys. Rev. D 67, 063509 (2003), arXiv:astro-ph/0209395 .
  * Gorini  _et al._ (2005) V. Gorini, A. Kamenshchik, U. Moschella, V. Pasquier, and A. Starobinsky, Phys. Rev. D 72, 103518 (2005), arXiv:astro-ph/0504576 .
  * Bento  _et al._ (2002) M. C. Bento, O. Bertolami, and A. A. Sen, Phys. Rev. D 66, 043507 (2002), arXiv:gr-qc/0202064 .
  * Xu and Lu (2010) L. Xu and J. Lu, JCAP 03, 025 (2010), arXiv:1004.3344 [astro-ph.CO] .
  * Lu  _et al._ (2009) J. Lu, Y. Gui, and L.-X. Xu, Eur. Phys. J. C 63, 349 (2009), arXiv:1004.3365 [astro-ph.CO] .
  * Benaoum (2022) H. Benaoum, Universe 8, 340 (2022), arXiv:hep-th/0205140 .
  * Debnath  _et al._ (2004) U. Debnath, A. Banerjee, and S. Chakraborty, Class. Quant. Grav. 21, 5609 (2004), arXiv:gr-qc/0411015 .
  * Lu  _et al._ (2008) J. Lu, L. Xu, J. Li, B. Chang, Y. Gui, and H. Liu, Phys. Lett. B 662, 87 (2008), arXiv:1004.3364 [astro-ph.CO] .
  * Xu  _et al._ (2012) L. Xu, Y. Wang, and H. Noh, Eur. Phys. J. C 72, 1931 (2012), arXiv:1204.5571 [astro-ph.CO] .
  * Escamilla-Rivera  _et al._ (2020) C. Escamilla-Rivera, M. A. C. Quintero, and S. Capozziello, JCAP 03, 008 (2020), arXiv:1910.02788 [astro-ph.CO] .
  * Hova and Yang (2017) H. Hova and H. Yang, Int. J. Mod. Phys. D 27, 1750178 (2017), arXiv:1011.4788 [astro-ph.CO] .
  * Hernandez-Almada  _et al._ (2019) A. Hernandez-Almada, J. Magana, M. A. Garcia-Aspeitia, and V. Motta, Eur. Phys. J. C 79, 12 (2019), arXiv:1805.07895 [astro-ph.CO] .
  * Yang  _et al._ (2019a) W. Yang, S. Pan, A. Paliathanasis, S. Ghosh, and Y. Wu, Mon. Not. Roy. Astron. Soc. 490, 2071 (2019a), arXiv:1904.10436 [gr-qc] .
  * Sandvik  _et al._ (2004) H. Sandvik, M. Tegmark, M. Zaldarriaga, and I. Waga, Phys. Rev. D 69, 123524 (2004), arXiv:astro-ph/0212114 .
  * Yadav  _et al._ (2021) A. K. Yadav, A. K. Yadav, M. Singh, R. Prasad, N. Ahmad, and K. P. Singh, Phys. Rev. D 104, 064044 (2021).
  * Normann and Brevik (2021) B. D. Normann and I. H. Brevik, Mod. Phys. Lett. A 36, 2150198 (2021), arXiv:2107.13533 [gr-qc] .
  * Yang  _et al._ (2019b) W. Yang, S. Pan, E. Di Valentino, A. Paliathanasis, and J. Lu, Phys. Rev. D100, 103518 (2019b), arXiv:1906.04162 [astro-ph.CO] .
  * Maggiore  _et al._ (2020) M. Maggiore _et al._ , JCAP 03, 050 (2020), arXiv:1912.02622 [astro-ph.CO] .
  * Zhao  _et al._ (2011) W. Zhao, C. Van Den Broeck, D. Baskaran, and T. G. F. Li, Phys. Rev. D 83, 023005 (2011), arXiv:1009.0206 [astro-ph.CO] .
  * Cai and Yang (2017) R.-G. Cai and T. Yang, Phys. Rev. D 95, 044024 (2017), arXiv:1608.08008 [astro-ph.CO] .
  * Zhang  _et al._ (2017) X. Zhang, J. Yu, T. Liu, W. Zhao, and A. Wang, Phys. Rev. D 95, 124008 (2017), arXiv:1703.09853 [gr-qc] .
  * Cai  _et al._ (2018) R.-G. Cai, T.-B. Liu, X.-W. Liu, S.-J. Wang, and T. Yang, Phys. Rev. D 97, 103005 (2018), arXiv:1712.00952 [astro-ph.CO] .
  * Yang  _et al._ (2020a) W. Yang, S. Pan, E. Di Valentino, B. Wang, and A. Wang, JCAP 05, 050 (2020a), arXiv:1904.11980 [astro-ph.CO] .
  * Yang  _et al._ (2019c) W. Yang, S. Vagnozzi, E. Di Valentino, R. C. Nunes, S. Pan, and D. F. Mota, JCAP 07, 037 (2019c), arXiv:1905.08286 [astro-ph.CO] .
  * Bachega  _et al._ (2020) R. R. A. Bachega, A. A. Costa, E. Abdalla, and K. S. F. Fornazier, JCAP 05, 021 (2020), arXiv:1906.08909 [astro-ph.CO] .
  * Zhou  _et al._ (2019) L. Zhou, X. Fu, Z. Peng, and J. Chen, Phys. Rev. D 100, 123539 (2019), arXiv:1912.02327 [gr-qc] .
  * Yang  _et al._ (2020b) W. Yang, S. Pan, D. F. Mota, and M. Du, Mon. Not. Roy. Astron. Soc. 497, 879 (2020b), arXiv:2001.02180 [astro-ph.CO] .
  * Mitra  _et al._ (2021) A. Mitra, J. Mifsud, D. F. Mota, and D. Parkinson, Mon. Not. Roy. Astron. Soc. 502, 5563 (2021), arXiv:2010.00189 [astro-ph.CO] .
  * Escamilla-Rivera and Quevedo (2021) C. Escamilla-Rivera and H. Quevedo, Class. Quant. Grav. 38, 115009 (2021), arXiv:2005.08784 [gr-qc] .
  * Pan  _et al._ (2021) Y. Pan, Y. He, J. Qi, J. Li, S. Cao, T. Liu, and J. Wang, Astrophys. J. 911, 135 (2021), arXiv:2103.05212 [astro-ph.CO] .
  * Zhang  _et al._ (2019) X.-N. Zhang, L.-F. Wang, J.-F. Zhang, and X. Zhang, Phys. Rev. D 99, 063510 (2019), arXiv:1804.08379 [astro-ph.CO] .
  * Wang  _et al._ (2018) L.-F. Wang, X.-N. Zhang, J.-F. Zhang, and X. Zhang, Phys. Lett. B 782, 87 (2018), arXiv:1802.04720 [astro-ph.CO] .
  * Di Valentino  _et al._ (2018) E. Di Valentino, D. E. Holz, A. Melchiorri, and F. Renzi, Phys. Rev. D 98, 083523 (2018), arXiv:1806.07463 [astro-ph.CO] .
  * Du  _et al._ (2019) M. Du, W. Yang, L. Xu, S. Pan, and D. F. Mota, Phys. Rev. D 100, 043535 (2019), arXiv:1812.01440 [astro-ph.CO] .
  * Matos  _et al._ (2021) I. S. Matos, M. O. Calvão, and I. Waga, Phys. Rev. D 103, 104059 (2021), arXiv:2104.10305 [gr-qc] .
  * Allahyari  _et al._ (2022) A. Allahyari, R. C. Nunes, and D. F. Mota, Mon. Not. Roy. Astron. Soc. 514, 1274 (2022), arXiv:2110.07634 [astro-ph.CO] .
  * Corman  _et al._ (2021) M. Corman, C. Escamilla-Rivera, and M. A. Hendry, JCAP 02, 005 (2021), arXiv:2004.04009 [astro-ph.CO] .
  * Belgacem  _et al._ (2019) E. Belgacem _et al._ (LISA Cosmology Working Group), JCAP 07, 024 (2019), arXiv:1906.01593 [astro-ph.CO] .
  * Cai  _et al._ (2017) R.-G. Cai, N. Tamanini, and T. Yang, JCAP 05, 031 (2017), arXiv:1703.07323 [astro-ph.CO] .
  * Hou  _et al._ (2021) S. Hou, X.-L. Fan, and Z.-H. Zhu, Mon. Not. Roy. Astron. Soc. 507, 761 (2021), arXiv:2106.01765 [astro-ph.CO] .
  * Zheng  _et al._ (2021) X. Zheng, S. Cao, Y. Liu, M. Biesiada, T. Liu, S. Geng, Y. Lian, and W. Guo, Eur. Phys. J. C 81, 14 (2021), arXiv:2012.14607 [astro-ph.CO] .
  * Luo  _et al._ (2016) J. Luo _et al._ (TianQin), Class. Quant. Grav. 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Barrow (1988) J. D. Barrow, Nucl. Phys. B 310, 743 (1988).
  * Barrow (1986) J. D. Barrow, Phys. Lett. B 180, 335 (1986).
  * Barrow (1990) J. D. Barrow, Phys. Lett. B 235, 40 (1990).
  * Hu (1998) W. Hu, Astrophys. J. 506, 485 (1998), arXiv:astro-ph/9801234 .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A8 (2020b), arXiv:1807.06210 [astro-ph.CO] .
  * Aghanim  _et al._ (2020c) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A5 (2020c), arXiv:1907.12875 [astro-ph.CO] .
  * Scolnic  _et al._ (2018) D. Scolnic _et al._ , Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * Sathyaprakash  _et al._ (2012) B. Sathyaprakash _et al._ , Class. Quant. Grav. 29, 124013 (2012), [Erratum: Class.Quant.Grav. 30, 079501 (2013)], arXiv:1206.0331 [gr-qc] .
  * Abadie  _et al._ (2010) J. Abadie _et al._ (LIGO Scientific), Nucl. Instrum. Meth. A 624, 223 (2010), arXiv:1007.3973 [gr-qc] .
  * Sathyaprakash  _et al._ (2010) B. S. Sathyaprakash, B. F. Schutz, and C. Van Den Broeck, Class. Quant. Grav. 27, 215006 (2010), arXiv:0906.4151 [astro-ph.CO] .
  * Schneider  _et al._ (2001) R. Schneider, V. Ferrari, S. Matarrese, and S. F. Portegies Zwart, Mon. Not. Roy. Astron. Soc. 324, 797 (2001), arXiv:astro-ph/0002055 .
  * Cutler and Holz (2009) C. Cutler and D. E. Holz, Phys. Rev. D 80, 104009 (2009), arXiv:0906.3752 [astro-ph.CO] .
  * Lewis and Bridle (2002) A. Lewis and S. Bridle, Phys. Rev. D 66, 103511 (2002), arXiv:astro-ph/0205436 .
  * Gelman and Rubin (1992) A. Gelman and D. B. Rubin, Statist. Sci. 7, 457 (1992).