# Bulk viscous fluid in symmetric teleparallel cosmology: theory versus
experiment

Raja Solanki 0000-0001-8849-7688<EMAIL_ADDRESS>Department of
Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. Simran Arora 0000-0003-0326-8945
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute of
Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. P.K. Sahoo 0000-0003-2130-8832
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute
of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. P.H.R.S. Moraes
<EMAIL_ADDRESS>Universidade Federal do ABC (UFABC)-Centro de Ciências
Naturais e Humans (CCNH)-Avenida dos Estados 5001, 09210-580, Santo André, SP,
Brazil.

###### Abstract

The standard formulation of General Relativity Theory, in the absence of a
cosmological constant, is unable to explain the responsible mechanism for the
observed late-time cosmic acceleration. On the other hand, by inserting the
cosmological constant in Einstein’s field equations it is possible to describe
the cosmic acceleration, but the cosmological constant suffers from an
unprecedented fine-tunning problem. This motivates one to modify Einstein’s
space-time geometry of General Relativity. The $f(Q)$ modified theory of
gravity is an alternative theory to General Relativity, where the non-
metricity scalar $Q$ is the responsible candidate for gravitational
interactions. In the present work we consider a Friedmann-Lemâitre-Robertson-
Walker cosmological model dominated by bulk viscous cosmic fluid in $f(Q)$
gravity with the functional form $f(Q)=\alpha Q^{n}$, where $\alpha$ and $n$
are free parameters of the model. We constrain our model with the recent
Pantheon supernovae data set of 1048 data points, Hubble data set of 31 data
points and baryon acoustic oscillations data set consisting of six points. For
higher values of redshift, it is clear that the $f(Q)$ cosmology better fits
data than standard cosmology. We present the evolution of our deceleration
parameter with redshift and it properly predicts a transition from decelerated
to accelerated phases of the universe expansion. Also, we present the
evolution of density, bulk viscous pressure and the effective equation of
state parameter with redshift. Those show that bulk viscosity in a cosmic
fluid is a valid candidate to acquire the negative pressure to drive the
cosmic expansion efficiently.We also examine the behavior of different energy
conditions to test the viability of our cosmological $f(Q)$ model.
Furthermore, the statefinder diagnostics are also investigated in order to
distinguish among different dark energy models.

##  I Introduction

The acceleration of the universe expansion is one of the most active
discoveries of modern cosmology. Observational studies include Type Ia
supernovae [1, 2], large scale structure [3, 4], baryon acoustic oscillations
[5, 6] and cosmic microwave background radiation [7, 8]. The reason behind the
late-time acceleration is a mystery. Several models hypothesize the existence
of a component called dark energy, which makes up around 70% of the entire
universe and could possess the feature of speeding up the universe expansion.

The cosmological constant $\Lambda$ in General Relativity (GR) field equations
plays the role of dark energy, i.e., a fluid with constant energy density and
high negative pressure. There are some issues with the cosmological constant
model, like the cosmic coincidence problem [9], which is the fact that the
density of non-relativistic matter and dark energy are the same order today. A
more delicate issue surrounding the cosmological constant is the so-called
cosmological constant problem, which is the high discrepancy between the
astronomically observed value of $\Lambda$ [1, 2] and the Particle Physics
theoretically predicted value of the quantum vacuum energy [10].

With the main purpose of solving the above cosmological issues, dynamical
(time-varying) dark energy models such as Chaplygin gas model [11, 12],
k-essence [13, 14], quintessence [15, 16] and decaying vacuum models [17, 18,
19, 20] have been proposed for some time in the literature.

Modified theories of gravity have also been intensively investigated to
understand the origin of the cosmic acceleration as well as to address the
cosmological constant model problems. It is possible to predict late-time
cosmic acceleration by modifying GR action. Some possibilities can be seen
within $f(R)$ [21, 22, 23], $f(G)$ [24, 25], $f(R,\mathcal{T})$ [26] and
$f(T)$ [27, 28, 29] theories of gravitation, with $R$, $G$, $\mathcal{T}$ and
$T$ being respectively the Ricci, Gauss-Bonnet, energy-momentum and torsion
scalars.

In the present article we will work with the recently introduced $f(Q)$ theory
of gravity [30], for which $Q$ is the non-metricity scalar, to be presented
below. The non-metricity formulation has been discussed earlier by Hehl and
Ne’eman (see References [31, 32, 33, 34]). The symmetric teleparallel gravity
is proven in the so-called coincidence gauge by imposing that the connection
is symmetric [35]. Weyl geometry is also observed to be a particular example
of Weyl Cartan geometry in which torsion disappears. The non-metricity is
interpreted as a massless spin 3-field in the case of symmetric connections
[36, 37]. Also, it is noted in the literature that due to the appearance of
non-metricity, the light cone structure is not preserved during parallel
transport [41]. Further, fermions are an issue in TEGR because they couple to
the axial contorsion of the Weitzenbock connection. This difficulty is
eliminated in STEGR since Dirac fermions only couple to the completely
antisymmetric component of the affine connection and is unaffected by any
disformation piece.

Although recently proposed, the $f(Q)$ gravity theory already presents some
interesting and valuable applications in the literature [38, 39, 40]. The
first cosmological solutions in $f(Q)$ gravity appear in References [42, 43],
while $f(Q)$ cosmography and energy conditions can respectively be seen in
[44, 45].

Here we are going to consider $f(Q)$ cosmology in the presence of a viscous
fluid. When a cosmic fluid expands too fast, the recovering of thermodynamic
equilibrium generates an effective pressure. The high viscosity in a cosmic
fluid is the manifestation of such an effective pressure [46, 47].

Basically, there are two viscosity coefficients, namely shear viscosity and
bulk viscosity. Shear viscosity is related to velocity gradients in the fluid,
and by considering the universe as described by homogeneous and isotropic
Friedmann-Lemâitre-Robertson-Walker (FLRW) background, it can be omitted.
Anyhow, by dropping the FLRW background assumption, several cosmological
models with shear viscosity fluid have been constructed, as one can check, for
instance, [48, 49, 50, 51]. On the other hand, bulk viscosity, which we are
going to consider here, introduces damping associated with volumetric
straining. To get in touch with bulk viscous fluid cosmological models, one
can check References [52, 53, 54, 55, 56]. Moreover, some interesting
applications of bulk viscous cosmology in black holes presented in [57, 58].

Researchers examine dark energy reconstruction with numerous observations as
data increases. The majority of studies has been concentrated on observable
evidences from Type Ia supernovae, cosmic microwave background and baryon
acoustic oscillations (BAO), which are known to be helpful in constraining
cosmological models. The Hubble parameter dataset shows the intricate
structure of the expansion of the universe. The ages of the most massive and
slowly evolving galaxies offer direct measurements of the Hubble parameter
$H(z)$ at various redshifts $z$, resulting in the development of a new form of
standard cosmological probe [59].

In our present work we include 31 measurements of Hubble expansion spanned
using differential age method [60] and BAO data consisting of six points [61].
Scolnic et al. recently published a large Type Ia supernovae sample named
Pantheon, with 1048 points and covering the redshift range $0.01<z<2.3$ [62].
Our analysis uses the $H(z)$, BAO and Pantheon samples to constrain the
cosmological model.

This work aims to describe the recently observed late-time acceleration with
the help of bulk viscosity of cosmic fluid (without including any dark energy
component) in the framework of $f(Q)$ theory of gravity. The manuscript is
organized as follows: in Sec. II we discuss the $f(Q)$ gravity formalism. In
Sec. III we describe the FLRW universe dominated by bulk viscous non-
relativistic matter and derive the expression for the Hubble parameter and the
deceleration parameter. Further, in Sec. IV, we analyze the observational data
to find the best-fit ranges for the parameters using the Hubble data set
containing 31 points, BAO sample and the Pantheon data set of 1048 samples.
Moreover, we analyze the behavior of different cosmological parameters such as
Hubble, density, effective pressure, deceleration parameter and effective
equation of state (EoS) parameter. In Sec. V, we investigate the consistency
of our bulk viscous fluid model by analyzing the different energy conditions.
In Sec. VI, we analyze the behavior of statefinder parameters on the values
constrained by the observational data to differentiate between dark energy
models. Lastly, we discuss our results in Sec. VII.

##  II Fundamental formulations in $f(Q)$ gravity

In $f(Q)$ gravity theory, the spacetime is established with the help of non-
metricity and symmetric teleparallelism condition, i.e.
$\nabla_{\alpha}g_{\mu\nu}\neq 0$ and $R^{\rho}_{\sigma\mu\nu}=0$. The
associated affine connection is given by

$\Upsilon^{\alpha}_{\ \mu\nu}=\Gamma^{\alpha}_{\ \mu\nu}+L^{\alpha}_{\
\mu\nu},$ (1)

with

$\displaystyle\Gamma^{\alpha}_{\
\mu\nu}\equiv\frac{1}{2}g^{\alpha\lambda}(g_{\mu\lambda,\nu}+g_{\lambda\nu,\mu}-g_{\mu\nu,\lambda}),$
(2) $\displaystyle L^{\alpha}_{\ \mu\nu}\equiv\frac{1}{2}(Q^{\alpha}_{\
\mu\nu}-Q_{\mu\ \nu}^{\ \alpha}-Q_{\nu\ \mu}^{\ \alpha}),$ (3)

being the Christoffel symbols and the distortion tensor, respectively, where

$\displaystyle Q_{\alpha\mu\nu}\equiv\nabla_{\alpha}g_{\mu\nu},$ (4)

is the non-metricity tensors.

The non-metricity tensor given by Eq.(4) has following two traces,

$Q_{\alpha}=Q_{\alpha}\>^{\mu}\>{}_{\mu},$ (5)

$\tilde{Q}_{\alpha}=Q^{\mu}\>_{\alpha\mu}.$ (6)

In addition, the superpotential tensor or non-metricity conjugate is given by

$\displaystyle 4P^{\lambda}\>_{\mu\nu}$ $\displaystyle=$
$\displaystyle-Q^{\lambda}_{\ \mu\nu}+Q_{\mu\ \nu}^{\ \lambda}+Q_{\nu\ \mu}^{\
\lambda}+(Q^{\lambda}-\tilde{Q}^{\lambda})g_{\mu\nu}$ (7) $\displaystyle-$
$\displaystyle\frac{1}{2}(\delta^{\lambda}_{\ \mu}Q_{\nu}+\delta^{\lambda}_{\
\nu}Q_{\mu}).$

Then the trace of the non-metricity tensor can be acquired as [63]

$Q=-Q_{\lambda\mu\nu}P^{\lambda\mu\nu}.$ (8)

The definition of the energy-momentum tensor for matter is

$\mathcal{T}_{\mu\nu}=\frac{-2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}L_{m})}{\delta
g^{\mu\nu}}.$ (9)

Furthermore, one can obtain the following relation between the curvature
tensor $R^{\rho}_{\sigma\mu\nu}$ and $\mathring{R}^{\rho}_{\sigma\mu\nu}$
corresponding to the connection $\Upsilon$ and $\Gamma$ as

$R^{\rho}_{\sigma\mu\nu}=\mathring{R}^{\rho}_{\sigma\mu\nu}+\mathring{\nabla}_{\mu}L^{\rho}_{\nu\sigma}-\mathring{\nabla}_{\nu}L^{\rho}_{\mu\sigma}+L^{\rho}_{\mu\lambda}L^{\lambda}_{\nu\sigma}-L^{\rho}_{\nu\lambda}L^{\lambda}_{\mu\sigma}$
(10)

and so

$R_{\sigma\nu}=\mathring{R}_{\sigma\nu}+\frac{1}{2}\mathring{\nabla}_{\nu}Q_{\sigma}+\mathring{\nabla}_{\rho}L^{\rho}_{\nu\sigma}-\frac{1}{2}Q_{\lambda}L^{\lambda}_{\nu\sigma}-L^{\rho}_{\sigma\lambda}L^{\lambda}_{\rho\sigma}$
(11)

$R=\mathring{R}+\mathring{\nabla}_{\lambda}Q^{\lambda}-\mathring{\nabla}_{\lambda}\tilde{Q}^{\lambda}-\frac{1}{4}Q_{\lambda}Q^{\lambda}+\frac{1}{2}Q_{\lambda}\tilde{Q}^{\lambda}-L_{\rho\nu\lambda}L^{\lambda\rho\nu}$
(12)

Now, the connection (1) can be parameterized as [30]

$\Upsilon^{\alpha}\,_{\mu\beta}=\frac{\partial
x^{\alpha}}{\partial\xi^{\rho}}\partial_{\mu}\partial_{\beta}\xi^{\rho}.$ (13)

Here, $\xi^{\alpha}=\xi^{\alpha}(x^{\mu})$ is an invertible relation. Hence,
it is always possible to find a coordinate system so that the connection
$\Upsilon^{\alpha}_{\ \mu\nu}$ vanishes. This situation is called coincident
gauge and the covariant derivative $\nabla_{\alpha}$ reduces to the partial
one $\partial_{\alpha}$. But in any other coordinate system in which this
affine connection does not vanish, the metric evolution will be affected and
result in a completely different theory [65, 64]. Thus in the coincident gauge
coordinate , we have

$Q_{\alpha\mu\nu}=\partial_{\alpha}g_{\mu\nu}$ (14)

while in an arbitrary coordinate system,

$Q_{\alpha\mu\nu}=\partial_{\alpha}g_{\mu\nu}-2\Upsilon^{\lambda}_{\alpha(\mu}g_{\nu)\lambda}.$
(15)

Locally, GR does not distinguish between gravitational and inertial effects,
however, by invoking frame fields, it is possible to covariantly define
gravitational energy in the teleparallel approach [66]. The canonical frame is
identified in the absence of both curvature and torsion and the canonical
coordinates in the absence of inertial effects. New physics can emerge from
such a formalism. Symmetric teleparallel gravity is broadly derived in three
gravity theories based on the coordinate transformations [67]. In this case,
we use the spatially flat case, i.e. $f(Q)$ gravity, whose field equations are
much easier to understand.

The $f(Q)$ gravity is described by the action [30]:

$S=\int{\frac{1}{2}f(Q)\sqrt{-g}d^{4}x}+\int{L_{m}\sqrt{-g}d^{4}x},$ (16)

in which $f(Q)$ is an arbitrary function of the non-metricity, $g=\text{det}\
g_{\mu\nu}$ and $L_{m}$ is the lagrangian density of matter. One can check in
Reference [30] that the functional form $f(Q)=-Q$ corresponds to the STEGR
(symmetric teleparallel equivalent to General Relativity) limit.

The gravitational field equations obtained by varying action (16) with respect
to the metric are

$\frac{2}{\sqrt{-g}}\nabla_{\lambda}(\sqrt{-g}f_{Q}P^{\lambda}\>_{\mu\nu})+\frac{1}{2}g_{\mu\nu}f+f_{Q}(P_{\mu\lambda\beta}Q_{\nu}\>^{\lambda\beta}-2Q_{\lambda\beta\mu}P^{\lambda\beta}\>_{\nu})=-\mathcal{T}_{\mu\nu},$
(17)

in which we defined $f_{Q}=df/dQ$.

Moreover, by varying the action with respect to the connection, one obtains

$\nabla_{\mu}\nabla_{\nu}(\sqrt{-g}f_{Q}P^{\mu\nu}\>_{\lambda})=0.$ (18)

##  III The cosmological model

We consider the flat FLRW metric for our analysis [68], such that:

$ds^{2}=-dt^{2}+a^{2}(t)(dx^{2}+dy^{2}+dz^{2}).$ (19)

In Equation (19) above, $a(t)$ is the cosmic scale factor, as usually.

From now, we will fix the coincident gauge so that connection becomes trivial
and the metric is only a fundamental variable.

We are going to assume a bulk viscous fluid and below we present some
considerations favoring such an assumption.

Firstly, to consider bulk viscosity in a fluid can be seen as an attempt to
refine its description, minimizing its ideal properties. This can be checked,
for instance, in the stellar astrophysics realistic models in References [69,
70, 71].

Under conditions of spatial homogeneity and isotropy (which refers to the
cosmological principle, as one can check, for instance, Reference [72]), the
bulk viscous pressure is the unique admissible dissipative phenomenon. In a
gas dynamical model, the existence of an effective bulk pressure can be traced
back to a non-standard self-interacting force on the particles of the gas
[55]. The bulk viscosity contributes negatively to the total pressure, as one
can chegk, for instance References [73, 74, 75].

Due to spatial isotropy, the bulk viscous pressure is the same in all spatial
directions and hence proportional to the volume expansion $\theta=3H$, with
$H=\dot{a}/a$ being the Hubble parameter and a dot represents time derivative.

The effective pressure of the cosmic fluid becomes [76, 77, 78]

$\bar{p}=p-\zeta\theta=p-3\zeta H,$ (20)

in which $p$ is the usual pressure and $\zeta>0$ is the bulk viscosity
coefficient, which we will assume as a free parameter of the model.

The corresponding energy-momentum tensor is given by

$\mathcal{T}_{\mu\nu}=(\rho+\bar{p})u_{\mu}u_{\nu}+\bar{p}g_{\mu\nu},$ (21)

in which $\rho$ is the matter-energy density and the four-velocity $u^{\mu}$
is such that its components are $u^{\mu}=(1,0,0,0)$.

The relation between normal pressure and matter-energy density follows [79]
$p=(\gamma-1)\rho$, with $\gamma$ being a constant lying in the range
$0\leq\gamma\leq 2$. Then the effective equation of state for the bulk viscous
fluid is given by the following:

$\bar{p}=(\gamma-1)\rho-3\zeta H.$ (22)

The Friedmann-like equations for our $f(Q)$ gravitational model are obtained
from the substitution of Equations (19)-(22) into Equation (11) and read as
follows (check, for instance, References [63, 80])

$3H^{2}=\frac{1}{2f_{Q}}\left(-\rho+\frac{f}{2}\right),$ (23)
$\dot{H}+\left(3H+\frac{\dot{f_{Q}}}{f_{Q}}\right)H=\frac{1}{2f_{Q}}\left(\bar{p}+\frac{f}{2}\right).$
(24)

In particular, for $f(Q)=-Q$ we retrieve the usual Friedmann equations [80],
as expected, since as we have mentioned above, this particular choice for the
functional form of the function $f(Q)$ is the STEGR limit of the theory.

For our investigation of bulk viscosity fluid cosmological model, we consider
the following $f(Q)$ functional form:

$f(Q)=\alpha Q^{n},$ (25)

with $\alpha\neq 0$ and constant $n$. This particular functional form for
$f(Q)$ was motivated by a polinomial form applied, for instance, in Reference
[45].

For the above choice of the $f(Q)$ function (Equation (25)), we rewrite
Equations (23)-(24) as follows

$\rho=\alpha 6^{n}\left(\frac{1}{2}-n\right)H^{2n},$ (26)
$\dot{H}+\frac{3}{2n}H^{2}=\frac{6^{1-n}\bar{p}}{2\alpha n(2n-1)}H^{2(1-n)},$
(27)

in which for the former we have isolated $\rho$.

From Equation (27) and Equation (22), we have the following:

$\dot{H}+\frac{3\gamma}{2n}H^{2}=-\frac{6^{2-n}\zeta}{4\alpha
n(2n-1)}H^{3-2n}.$ (28)

Now, we replace the term $d/dt$ by $d/dlna$ via the expression $d/dt=Hd/dlna$,
such that Equation (28) becomes

$\frac{dH}{dlna}+\frac{3\gamma}{2n}H=-\frac{6^{2-n}\zeta}{4\alpha
n(2n-1)}H^{2(1-n)}.$ (29)

The integration of Equation (29) yields the following solution

$H(a)=\left\\{(H_{0}a^{-\frac{3\gamma}{2n}})^{2n-1}+\frac{6^{1-n}\zeta}{\gamma\alpha(2n-1)^{2}}[a^{-\frac{3\gamma(2n-1)}{2n}}-1]\right\\}^{\frac{1}{2n-1}},$
(30)

with $H_{0}$ being a constant of integration to be found below.

We obtain the Hubble parameter in terms of redshift by using relation [72]
$a(t)=1/(1+z)$ in Equation (31). By making $z=0$ in (31) we find that
$H(0)=H_{0}$. The deceleration parameter is defined as
$q=-\ddot{a}a/\dot{a}^{2}=-\ddot{a}/(H^{2}a)$. Henceforth from Equation (30)
we have

$\displaystyle H(z)$ $\displaystyle=$
$\displaystyle\left\\{[H_{0}(1+z)^{\frac{3\gamma}{2n}}]^{2n-1}+\frac{6^{1-n}\zeta}{\gamma\alpha(2n-1)^{2}}[(1+z)^{\frac{3\gamma(2n-1)}{2n}}-1]\right\\}^{\frac{1}{2n-1}},$
(31) $\displaystyle q(z)$ $\displaystyle=$
$\displaystyle\frac{3}{2n}\left\\{\frac{\zeta}{\alpha
6^{n-1}(2n-1)\left\\{[H_{0}(1+z)^{\frac{3\gamma}{2n}}]^{2n-1}+\frac{6^{1-n}\zeta}{\gamma\alpha(2n-1)^{2}}[(1+z)^{\frac{3\gamma(2n-1)}{2n}}-1]\right\\}}+\gamma\right\\}-1.$
(32)

##  IV Observational constraints

To examine the observational features of our cosmological model, we use the
most recent cosmic Hubble and Supernovae observations. We use 31 points of the
Hubble data sets, 6 points of the BAO data sets and 1048 points from the
Pantheon supernovae samples. We apply the Bayesian analysis and likelihood
function along with the Markov Chain Monte Carlo (MCMC) method in emcee python
library [81].

###  IV.1 Hubble datasets

The Hubble parameter can be expressed as $H(z)=-dz/[dt(1+z)]$. As $dz$ is
derived from a spectroscopic survey, the model-independent value of the Hubble
parameter may be calculated by measuring the quantity $dt$.

We incorporate the set of 31 data points that are measured from the
differential age approach [82] to avoid extra correlation with BAO data. The
mean values of the model parameters $\zeta$, $\alpha$, $\gamma$ and $n$ are
calculated using the chi-square function as follows:

$\chi_{H}^{2}(\zeta,\alpha,\gamma,n)=\sum\limits_{k=1}^{31}\frac{[H_{th}(\zeta,\alpha,\gamma,n,z_{k})-H_{obs}(z_{k})]^{2}}{\sigma_{H(z_{k})}^{2}}.$
(33)

Here, $H_{th}$ is the Hubble parameter value predicted by the model, $H_{obs}$
represents its observed value and the standard error in the observed value of
$H$ is $\sigma_{H(z_{k})}$.

From the Hubble dataset, we obtain the best fit values for $\zeta$, $\alpha$,
$\gamma$, $n$ as the $1-\sigma$ and $2-\sigma$ contour plots in Fig.1.

Figure 1: The $1-\sigma$ and $2-\sigma$ likelihood contours for the model
parameters using the Hubble datasets. Figure 2: The error bar plot of $H$
versus $z$ for the considered $f(Q)$ model. The solid red line is the curve
for $f(Q)$ model whereas the black dotted line represents the $\Lambda$CDM
model. The blue dots depict the 31 points of the Hubble data.

The best fit values of the model parameters are $\zeta=0.65^{+0.10}_{-0.10}$,
$\alpha=-0.0166^{+0.0079}_{-0.0053}$, $\gamma=1.36^{+0.14}_{-0.12}$ and
$n=0.974^{+0.019}_{-0.027}$.

Fig.2 shows the error bar plot of the considered model and $\Lambda$CDM or
standard cosmological model, with cosmological constant density parameter
$\Omega_{\Lambda_{0}}=0.7$, matter density parameter $\Omega_{m_{0}}=0.3$ and
$H_{0}=69$ km/s/Mpc.

###  IV.2 BAO datasets

The BAO distance dataset, which includes the 6dFGS, SDSS and WiggleZ surveys,
comprise BAO measurements at six different redshifts in Table 1. The
characteristic scale of BAO is ruled by the sound horizon $r_{s}$ at the epoch
of photon decoupling $z_{\ast}$ that is given by the following relation:

$r_{s}(z_{\ast})=\frac{c}{\sqrt{3}}\int_{0}^{\frac{1}{1+z_{\ast}}}\frac{da}{a^{2}H(a)\sqrt{1+(3\Omega_{b0}/4\Omega_{\gamma
0})a}}.$ (34)

Here, $\Omega_{b0}$ and $\Omega_{\gamma 0}$ correspond to the present
densities of baryons and photons respectively.

The following relations are used in BAO measurements

$\displaystyle\triangle\theta$ $\displaystyle=$
$\displaystyle\frac{r_{s}}{d_{A}(z)},$ (35) $\displaystyle d_{A}(z)$
$\displaystyle=$ $\displaystyle\int_{0}^{z}\frac{dz^{\prime}}{H(z^{\prime})},$
(36) $\displaystyle\triangle z$ $\displaystyle=$ $\displaystyle H(z)r_{s},$
(37)

where $\triangle\theta$ represents the measured angular separation, $d_{A}$ is
the angular diameter distance and $\triangle z$ represents the measured
redshift separation of the BAO feature in the 2 point correlation function of
the galaxy distribution on the sky along the line of sight.

In this work, BAO datasets of six points for $d_{A}(z_{\ast})/D_{V}(z_{BAO})$
is taken from the References [83, 84, 85, 86, 87, 88], where the redshift at
the epoch of photon decoupling is taken as $z_{\ast}\approx 1091$ and
$d_{A}(z)$ is the co-moving angular diameter distance together with the
dilation scale $D_{V}(z)=\left[d_{A}(z)^{2}z/H(z)\right]^{1/3}$. The chi-
square function for the BAO datasets is taken to be [88]

$\chi_{BAO}^{2}=X^{T}C^{-1}X\,,$ (38)

$X=\left(\begin{array}[]{c}\frac{d_{A}(z_{\star})}{D_{V}(0.106)}-30.95\\\
\frac{d_{A}(z_{\star})}{D_{V}(0.2)}-17.55\\\
\frac{d_{A}(z_{\star})}{D_{V}(0.35)}-10.11\\\
\frac{d_{A}(z_{\star})}{D_{V}(0.44)}-8.44\\\
\frac{d_{A}(z_{\star})}{D_{V}(0.6)}-6.69\\\
\frac{d_{A}(z_{\star})}{D_{V}(0.73)}-5.45\end{array}\right)\,,$

$z_{BAO}$ | $0.106$ | $0.2$ | $0.35$ | $0.44$ | $0.6$ | $0.73$
---|---|---|---|---|---|---
$\frac{d_{A}(z_{\ast})}{D_{V}(z_{BAO})}$ | $30.95\pm 1.46$ | $17.55\pm 0.60$ | $10.11\pm 0.37$ | $8.44\pm 0.67$ | $6.69\pm 0.33$ | $5.45\pm 0.31$
Table 1: Values of $d_{A}(z_{\ast})/D_{V}(z_{BAO})$ for distinct values of
$z_{BAO}$.

The inverse covariance matrix $C^{-1}$ is defined in [88]

$C^{-1}=\left(\begin{array}[]{cccccc}0.48435&-0.101383&-0.164945&-0.0305703&-0.097874&-0.106738\\\
-0.101383&3.2882&-2.45497&-0.0787898&-0.252254&-0.2751\\\
-0.164945&-2.454987&9.55916&-0.128187&-0.410404&-0.447574\\\
-0.0305703&-0.0787898&-0.128187&2.78728&-2.75632&1.16437\\\
-0.097874&-0.252254&-0.410404&-2.75632&14.9245&-7.32441\\\
-0.106738&-0.2751&-0.447574&1.16437&-7.32441&14.5022\end{array}\right)\,.$

Figure 3: The $1-\sigma$ and $2-\sigma$ likelihood contours for the model
parameters using the BAO datasets.

The values that fit observations are $\zeta=0.66^{+0.12}_{-0.12}$,
$\alpha=-0.0125^{+0.0047}_{-0.0023}$, $\gamma=1.191^{+0.070}_{-0.070}$ and
$n=0.996^{+0.010}_{-0.012}$ in Fig.3.

###  IV.3 Pantheon datasets

Scolnic et al. [62] put together the Pantheon samples consisting of 1048 type
Ia supernovae in the redshift range $0.01<z<2.3$. The PanSTARSS1 Medium Deep
Survey, SDSS, SNLS and numerous low-z and HST samples contribute to it. The
empirical relation used to calculate the distance modulus of SNeIa from the
observation of light curves is given by $\mu=m_{B}^{*}+\alpha X_{1}-\beta
C-M_{B}+\Delta_{M}+\Delta_{B}$, where $X_{1}$ and $C$ denote the stretch and
color correction parameters, respectively [62, 89], $m_{B}^{*}$ represents the
observed apparent magnitude and $M_{B}$ is the absolute magnitude in the
B-band for SNeIa. The parameters $\alpha$ and $\beta$ are the two nuisance
parameters describing the luminosity stretch and luminosity color relations,
respectively. Further, the distance correction factor is $\Delta_{M}$ and
$\Delta_{B}$ is a distance correction based on predicted biases from
simulations.

The nuisance parameters in the Tripp formula [90] were reconstructed using a
novel technique called BEAMS with Bias Corrections [91, 92] and the observed
distance modulus was reduced to the difference between the corrected apparent
magnitude $m_{B}$ and the absolute magnitude $M_{B}$, which is
$\mu=m_{B}-M_{B}$. We shall avoid marginalizing the over nuisance parameters
$\alpha$ and $\beta$ but marginalize over the Pantheon data for $M_{B}$.
Hence, we ignore the values of $\alpha$ and $\beta$ for the present
investigation of the model.

The luminosity distance read as

$D_{L}(z)=(1+z)\int_{0}^{z}\frac{cdz^{\prime}}{H(z^{\prime})},$ (39)

with $c$ being the speed of light.

The $\chi^{2}$ function for type Ia supernovae is obtained by correlating the
theoretical distance modulus

$\displaystyle\mu(z)$ $\displaystyle=$ $\displaystyle
5log_{10}D_{L}(z)+\mu_{0},$ (40) $\displaystyle\mu_{0}$ $\displaystyle=$
$\displaystyle 5log(1/H_{0}Mpc)+25,$ (41)

such that

$\chi^{2}_{SN}(\zeta,\alpha,\gamma,n)=\sum_{k=1}^{1048}\dfrac{\left[\mu_{obs}(z_{k})-\mu_{th}(\zeta,\alpha,\gamma,n,z_{k})\right]^{2}}{\sigma^{2}(z_{k})},$
(42)

where $\mu_{th}$ is the theoretical value of distance modulus, $\mu_{obs}$ is
the observed value whereas $\sigma^{2}(z_{k})$ is the standard error in the
observed value.

Using the Pantheon supernovae datasets, we obtain the best fit values for
$\zeta$, $\alpha$, $\gamma$ and $n$ as the $1-\sigma$ and $2-\sigma$ contour
plots in Fig.4. The values that fit the model are
$\zeta=0.67^{+0.12}_{-0.12}$, $\alpha=-0.00999^{+0.0047}_{-0.0024}$,
$\gamma=1.34^{+0.15}_{-0.12}$ and $n=1.001^{+0.024}_{-0.024}$.

Figure 4: The $1-\sigma$ and $2-\sigma$ likelihood contours for the model
parameters using the Pantheon datasets.

###  IV.4 Cosmological Parameters

The evolution of the Hubble parameter, deceleration parameter, energy density,
pressure with bulk viscosity and the effective EoS parameter for the redshift
range $-1<z<8$ are presented below, in order to test the late time cosmic
expansion history and the future of expanding universe [93]. In order to do so
we use the set of values constrained by Hubble, BAO and the Pantheon data sets
for the model parameters.

Figure 5: Profile of the Hubble parameter for the given model corresponding to
the values of the parameters constrained by the Hubble, BAO and the Pantheon
data point sets. Figure 6: Profile of the deceleration parameter for the given
model corresponding to the values of the parameters constrained by the Hubble,
BAO and the Pantheon data point sets. Figure 7: Profile of the density
parameter for the given model corresponding to the values of the parameters
constrained by the Hubble, BAO and the Pantheon data point sets. Figure 8:
Profile of the pressure for the given model corresponding to the values of the
parameters constrained by the Hubble, BAO and the Pantheon data point sets.
Figure 9: Profile of the EoS parameter for the given model corresponding to
the values of the parameters constrained by the Hubble, BAO and the Pantheon
data point sets.

From Fig.6, it is clear that the deceleration parameter shows the transition
from a decelerated ($q>0$) to an accelerated ($q<0$) phase of the universe
expansion for the constrained values of the model parameters. The transition
redshift is $z_{t}\thickapprox 0.142$, $z_{t}\thickapprox 0.776$ and
$z_{t}\thickapprox 0.622$ corresponding to the Hubble, BAO and the Pantheon
data sets, respectively. The present value of the deceleration parameter is,
respectively, $q_{0}=-0.127$ , $q_{0}=-0.436$ and $q_{0}=-0.454$.

From Fig. 5 and 7 it is clear that the Hubble and density parameter shows the
positive behavior for all the constrained values of the model parameters,
which is expected.

Fig 8 indicates that the bulk viscous cosmic fluid exhibits, for lower
redshifts, the negative pressure that make bulk viscosity to be a viable
candidate to drive the cosmic acceleration. Furthermore, the effective EoS
parameter presented in Fig.9 indicates that the cosmic viscous fluid behaves
like quintessence dark energy. The present values of EoS parameter
corresponding to the Hubble, BAO and the Pantheon samples are
$\omega_{0}=-0.433$ , $\omega_{0}=-0.625$, and $\omega_{0}=-0.635$.

##  V Energy Conditions

In the present section we are going to construct the energy conditions for the
solutions of the present model. The energy conditions are relations applied to
the matter energy-momentum tensor with the purpose of satisfying positive
energy. The energy conditions are derived from the Raychaudhuri equation and
are written as [94]

  * •

Null energy condition (NEC) : $\rho_{eff}+p_{eff}\geq 0$;

  * •

Weak energy condition (WEC) : $\rho_{eff}\geq 0$ and $\rho_{eff}+p_{eff}\geq
0$;

  * •

Dominant energy condition (DEC) : $\rho_{eff}\pm p_{eff}\geq 0$;

  * •

Strong energy condition (SEC) : $\rho_{eff}+3p_{eff}\geq 0$,

with $\rho_{eff}$ being the effective energy density.

In Figs.10 and 11 it is evident that the NEC and DEC exhibit positive behavior
for all the constrained values of the model parameters. As WEC is the
combination of energy density and NEC, we conclude that NEC, DEC and WEC are
all satisfied in the entire domain of redshift. Fig.12 indicates that the SEC
exhibits, for lower redshifts, the negative behavior that is related to cosmic
acceleration [95]. This is also reflected in the deceleration parameter
behavior in Fig.6.

Figure 10: Profile of the null energy condition for the given model
corresponding to the values of the parameters constrained by the Hubble, BAO
and the Pantheon data point sets. Figure 11: Profile of the dominant energy
condition for the given model corresponding to the values of the parameters
constrained by the Hubble, BAO and the Pantheon data point sets. Figure 12:
Profile of the strong energy condition for the given model corresponding to
the values of the parameters constrained by the Hubble, BAO and the Pantheon
data point sets.

##  VI Statefinder analysis

The cosmological constant $\Lambda$ suffers from two major drawbacks, namely
the aforementioned cosmological constant and cosmic coincidence problems. To
surmount these problems, dynamic models of dark energy have been introduced in
the literature as we have also mentioned in Introduction. To discriminate
between these time-varying dark energy models, an appropriate tool was
required. In this direction, V. Sahni et al. introduced a new pair of
geometrical parameters known as statefinder parameters $(r,s)$ [96]. The
statefinder parameters are defined as

$\displaystyle r$ $\displaystyle=$ $\displaystyle\frac{\dddot{a}}{aH^{3}},$
(43) $\displaystyle s$ $\displaystyle=$
$\displaystyle\frac{(r-1)}{3(q-\frac{1}{2})}.$ (44)

The parameter $r$ can be rewritten as $r=2q^{2}+q-\frac{\dot{q}}{H}$.

For different values of the statefinder pair $(r,s)$, it represents the
following dark energy models:

${r=1,s=0}$ represents $\Lambda$CDM model,
${r>1,s<0}$ represents Chaplygin gas model,
${r<1,s>0}$ represents quintessence model.

In Figures 13 and 14, we plot the $s-r$ and $q-r$ diagrams for our
cosmological model by taking the values of the parameters constrained by the
Hubble, BAO and the Pantheon data sets.

Figure 13: Plot of the trajectories in the $r-s$ plane for the given
cosmological model corresponding to the parameters values constrained by the
Hubble, BAO and the Pantheon data sets. Figure 14: Plot of the trajectories in
the $q-r$ plane for the given cosmological model corresponding to the
parameters values constrained by the Hubble, BAO, and the Pantheon data sets.

Figures 13 and 14 show that our bulk viscous model lies in the quintessence
region. Also, the evolutionary trajectories of our model departure from the
$\Lambda$CDM point. The present values of the statefinder parameters
corresponding to the values of the model parameters constrained by the Hubble,
BAO and the Pantheon samples are $r_{0}=0.837$ and $s_{0}=0.086$, $r_{0}=0.62$
and $s_{0}=0.135$, $r_{0}=0.762$ and $s_{0}=0.083$ respectively.

##  VII Discussions and conclusions

In the present section we will discuss the results obtained in Sections IV,V,
and VI for the bulk viscous symmetric teleparallel cosmological model here
developed and presented.

Cosmology has been on the agenda mainly for two reasons: dark energy and dark
matter. While dark energy has been deeply discussed throughout the paper, dark
matter is predicted within $\Lambda$CDM model as a sort of matter that does
not interact electromagnetically, so that it cannot be seen, but its
gravitational effects otherwise can well be detected. Still we have not yet
detected or even associated dark matter to a particle of standard model or
beyond [97, 98, 99]. Modified (or alternative) theories of gravity have also
been used to describe dark matter effects [100, 101]. In these cases, dark
matter is simply an effect of modification of gravity.

Returning to the dark energy question, it is highly counter-intuitive that the
expansion of the universe is actually accelerating. Although the vacuum
quantum energy can well explain this dynamical effect via the cosmological
constant in Einstein’s field equations of GR, the aforementioned important and
persistent problems related to $\Lambda$ supply the search for alternative
explanations.

In the present article, as an attempt to describe dark energy we have assumed
the symmetric teleparallel gravity as the underlying gravity theory.

The $f(Q)$ gravity was recently proposed by Jiménez et al. in [30] as a
ramification of the geometric trinity, that says that space-time manifold can
be described by curvature, torsion or non-metricity. Particularly, the
symmetric teleparallel gravity describes gravitational interactions via the
non-metricity scalar, with null curvature and torsion.

Our $f(Q)$ cosmological model was based on a spatially homogeneous and
isotropic flat metric and an energy-momentum tensor describing a bulk viscous
fluid. We let the $f(Q)$ function be a power of $n$ as $f(Q)\sim Q^{n}$, with
$n$ a free parameter.

In Section IV we started testing our cosmological solutions. We started
confronting the Hubble parameter (31) with 31 data points that are measured
from differential age approach. In Fig.1 we have obtained the best fit values
for the model free parameters. We have then plotted Fig.2, in which our $H(z)$
is confronted with cosmological data and compared with $\Lambda$CDM
prediction. We can see the $f(Q)$ model describes observations with good
agreement and specially for higher values of redshift it is clear that it
provides a better fit when compared to $\Lambda$CDM model. Further, in
Figs.(3) and (4), we have obtained the best fit values for the model free
parameters.

From Fig. 5 and 7 we found that the Hubble and density parameter shows the
expected positive behavior for all the constrained values of the model
parameters. Fig.8 indicates that the bulk viscous cosmic fluid exhibit the
negative pressure that make bulk viscosity to be a viable candidate to drive
the cosmic acceleration. This is also reflected in the deceleration parameter
behavior in Fig.6, which shows a transition from decelerated to accelerated
phases of the universe expansion. Furthermore, the effective EoS parameter
presented in Fig.9 indicates that the cosmic viscous fluid behaves like
quintessence dark energy.

In Section V, we investigated the consistency of our model by analyzing the
different energy conditions. We found that NEC, DEC, and WEC all are satisfied
in the entire domain of redshift (presented in Fig10 and 11) while the SEC,
presented in Fig12, is violated for lower redshifts that implies the cosmic
acceleration and satisfied for higher redshifts that implies a decelerated
phase of the universe.

In Section VI, Figs.13 and 14 show that the evolutionary trajectories of our
model is departed from $\Lambda$CDM fixed point ${r=1,s=0}$. In the present
epoch they lie in the quintessence region ${r<1,s>0}$. The present model is,
therefore, a good alternative to explain the universe dynamics, particularly
with no necessity of invoking the cosmological constant.

Finally, it is worth remarking that $f(Q)$ gravity still needs to be applied
in several different regimes in order to be stablished as a viable
gravitational formalism. In [30] it was shown that the $f(Q)$ gravity can be
compatible with Solar System constraints. Recently, in [102] it was shown that
the $f(Q)$ gravity predicts a modified gravitational wave propagation and a
different gravitational wave luminosity distance when compared to the standard
electromagnetic one, although this investigation was done for a different
functional form of the $f(Q)$ function.

## ACKNOWLEDGMENTS

RS acknowledges University Grants Commission (UGC), New Delhi, India for
awarding Junior Research Fellowship (UGC-Ref. No.: 191620096030). SA
acknowledges CSIR, New Delhi, India for SRF. PKS acknowledges CSIR, New Delhi,
India for financial support to carry out the Research project
[No.03(1454)/19/EMR-II Dt.02/08/2019] and IUCAA, Pune, India for providing
support through visiting associateship program.

## References

  * [1] A.G. Riess et al., Astron. J. 116, 1009(1998).
  * [2] S. Perlmutter et al., Astrophys. J. 517, 565(1999).
  * [3] T. Koivisto, D.F. Mota, Phys. Rev. D 73, 083502(2006).
  * [4] S.F. Daniel, Phys. Rev. D 77, 103513(2008).
  * [5] D.J. Eisenstein et al., Astrophys. J. 633, 560(2005).
  * [6] W.J. Percival at el., Mon. Not. R. Astron. Soc. 401, 2148(2010).
  * [7] R.R. Caldwell, M. Doran, Phys. Rev. D 69, 103517(2004).
  * [8] Z.Y. Huang et al., J. Cosm. Astrop. Phys. 0605, 013(2006).
  * [9] N. Dalal et al., Phys. Rev. Lett. 87, 141302(2001).
  * [10] S. Weinberg, Rev. Mod. Phys. 61, 1(1989).
  * [11] M. C. Bento et al., Phys. Rev. D 66, 043507(2002).
  * [12] A. Y. Kamenshchik et al., Phys. Lett. B 511, 265(2001).
  * [13] T. Chiba et al., Phys. Rev. D 62, 023511(2000).
  * [14] C. Armendariz-Picon et al., Phys. Rev. Lett. 85, 4438(2000).
  * [15] S. M. Carroll, Phys. Rev. Lett. 81, 3067(1998).
  * [16] Y. Fujii, Phys. Rev. D 26, 2580(1982).
  * [17] L. Xu et al., Phys. Rev. D 84, 123004(2011).
  * [18] M. Tong and H. Noh, Eur. Phys. J. C 71, 1586(2011).
  * [19] K. Freese et al., Nucl. Phys. B 287, 797(1987).
  * [20] A.-M.M. Abdel-Rahman, Phys. Rev. D 45, 3497(1992).
  * [21] S. Appleby and R. Battye, Phys. Lett. B 654, 7(2007).
  * [22] L. Amendola et al., Phys. Rev. D 75, 083504(2007).
  * [23] R. Saffari and S. Rahvar, Phys. Rev. D 77, 104028(2008).
  * [24] G. Cognola et al., Phys. Rev. D 73, 084007(2006).
  * [25] B. Li et al., Phys. Rev. D 76, 044027(2007).
  * [26] P.H.R.S. Moraes and J.R.L. Santos, Eur. Phys. J. C 76, 60(2016).
  * [27] X. Ren et al., Phys. Dark Univ. 32, 100812(2021).
  * [28] G.L. Nashed, Gen. Rel. Grav. 47, 75(2015).
  * [29] M.R. Setare and N. Mohammadipour, J. Cosm. Astrop. Phys. 01, 015(2013).
  * [30] J. B. Jiménez et al., Phys. Rev. D 98, 044048(2018).
  * [31] F. W. Hehl et al., Phys. Rep. 258, 1-171 (1995).
  * [32] F. W. Hehl et al., Phys. Lett. B 78, 102-106 (1978).
  * [33] F. W. Hehl and Y. N. Obukhov, arXiv, arXiv: 0001010 (2000).
  * [34] Y. Ne’eman and F. W. Hehl, Class. Quant. Grav. 14, 251-260 (1997).
  * [35] J.M. Nester and H.-J. Yo, Chin. J. Phys. 37, 113(1999).
  * [36] N. Boulanger and I. Kirsch, Phys. Rev. D 73, 124023 (2006).
  * [37] P. Baekler et al., Phys. Rev. D 74, 125009 (2006).
  * [38] N. Frusciante, Phys. Rev. D 103, 044021(2021).
  * [39] Bruno J. Barros et al., Phys. Dark Univ. 30, 100616(2020).
  * [40] J. Ferreira et al., Phys. Rev. D 105, 123531 (2022); I. S. Albuquerque, N. Frusciante, Phys. Dark Univ. 35, 100980 (2022); S. Capozziello, R. D’Agostino, Phys. Lett. B 832, 137229 (2022); L. Atayde, N. Frusciante, Phys. Rev. D 104, 064052 (2022); I. Ayuso, R. Lazkoz, V. Salzano, Phys. Rev. D 103, 063505 (2021); W. Khyllep et al., arxiv: 2207.02610; S. Sahlu, E. Tsegaye, arxiv:2206.02517.
  * [41] H. Weyl, Sitzungsber. Preuss. akad. Wiss. 465, 1 (1918).
  * [42] J.B. Jiménez et al., Phys. Rev. D 101, 103507(2020).
  * [43] W. Khyllep et al., Phys. Rev. D 103, 103521(2021).
  * [44] S. Mandal et al., Phys. Rev. D 102, 124029(2020).
  * [45] S. Mandal et al., Phys. Rev. D 102, 024057(2020).
  * [46] J. R. Wilson et al., Phys. Rev. D 75, 043521(2007).
  * [47] H. Okumura, F. Yonezawa, Physica A 321, 207(2003).
  * [48] R. Bali and D.R. Jain, Astrophys. Spa. Sci. 141, 207(1988).
  * [49] R. Bali and D.R. Jain, Astriphys. Spa. Sci. 139, 175(1987).
  * [50] Y. Deng and P.D. Mannheim, Phys. Rev. D 44, 1722(1991).
  * [51] W.-H. Huang, J. Math. Phys. 31, 659(1990).
  * [52] G.C. Samanta et al., Zeits. Naturfor. 72, 365(2017).
  * [53] J. Satish and R. Venkateswarlu, Chin. J. Phys. 54, 830(2016).
  * [54] A. Beesham, Phys. Rev. D 48, 3539(1993).
  * [55] R. Colistete Jr. et al., Phys. Rev. D 76, 103516(2007).
  * [56] S. Davood Sadatian, EPL, 126, 30004(2019); I. Brevik, Entropy, 14(11), 2302-2310(2012); C.P. Singh, P. Kumar,  Eur. Phys. J. C, 74, 3070(2014); M. Srivastava, C.P. Singh, Astrophys. Space Sci., 363, 117 (2018).
  * [57] F. D. Paolis, M. Jamil, and A. Qadir, Int. J. Theor. Phys. 49, 621-632 (2010).
  * [58] I. Brevik, M. Jamil, Int. J. Geom. Meth. Mod. Phys. 16, 1950030 (2019).
  * [59] R. Jimenez, A. Loeb, ApJ 573, 37(2002).
  * [60] G. S. Sharov et al., Mon. Not. R. Astron. Soc. 466, 3497(2017).
  * [61] C. Blake et al., Mon. Not. R. Astron. Soc. 418, 1707(2011).
  * [62] D.M. Scolnic et al., ApJ 859, 101(2018).
  * [63] R. Lazkoz et al.,  Phys. Rev. D 100, 104027 (2019).
  * [64] J.T. Beh, T.H. Loo, and A. De, arXiv, arXiv: 2107.04513 (2021).
  * [65] N. Dimakis, A. Paliathanasis, and T. Christodoulakis,  arXiv, arXiv:2108.01970 (2021).
  * [66] J.W. Maluf,  Ann. Phys. 525, 339(2013).
  * [67] Manuel Hohmann, Phys. Rev. D 104, 124077(2021).
  * [68] S. Vagnozzi, A. Loeb., and M. Moresco,  Astrophys. J. 908, 84 (2021).
  * [69] M.E. Gusakov,  Phys. Rev. D 76, 083001(2007).
  * [70] M.E. Gusakov and E.M. Kantor,  Phys. Rev. D 78, 083006(2008).
  * [71] P. Haensel et al.,  Astron. Astrophys. 381, 1080(2002).
  * [72] B. Ryden, Introduction to Cosmology (Addison Wesley, SanFrancisco, United States of America, 2003).
  * [73] S.D. Odintsov et al.,  Phys. Rev. D 101, 044010(2020).
  * [74] J.C. Fabris et al.,  Gen. Rel. Grav. 38, 495(2006).
  * [75] X.-H. Meng and X. Dou,  Comm. Theor. Phys. 52, 377(2009).
  * [76] I. Brevik and O. Gorbunova, Gen. Rel. Grav. 37, 2039(2005).
  * [77] O. Gron,  Astrophys.Space Sci. 173, 191-225 (1990).
  * [78] C. Eckart,  Phys. Rev. 58, 919 (1940).
  * [79] J. Ren and X.-H. Meng,  Physics Letters B 633, 1-8 (2006).
  * [80] T. Harko et al.,  Phys.Rev.D 98, 084043 (2018).
  * [81] D. F. Mackey et al., Publ. Astron. Soc. Pac. 125, 306(2013).
  * [82] G.S. Sharov and V.O. Vasilie, Mathematical Modelling and Geometry 6 1(2018).
  * [83] C. Blake et al.,  Mon. Not. Roy. Astron. Soc. 418, 1707 (2011).
  * [84] W. J. Percival et al., Mon. Not. Roy. Astron. Soc. 401, 2148 (2010).
  * [85] F. Beutler et al.,  Mon. Not. Roy. Astron. Soc. 416, 3017 (2011).
  * [86] N. Jarosik et al.,  Astrophys. J. Suppl. 192, 14 (2011).
  * [87] D. J. Eisenstein et al.,  Astrophys. J. 633, 560 (2005).
  * [88] R. Giostri et al., J. Cosm. Astropart. Phys. 1203, 027 (2012).
  * [89] P. Mukherjee, N. Banerjee, Eur. Phys. J. C 81, 36(2021).
  * [90] R. Tripp, Astron. Astrophys. 331, 815(1998).
  * [91] R. Kessler, D. Scolnic, Astrophys. J. 836, 56(2017).
  * [92] Fotios K. Anagnostopoulos et al., Phys. Rev. D 103, 104013(2021).
  * [93] S. Vagnozzi, F. Pacucci, and A. Loeb., JHEAp 36, 27-35 (2022).
  * [94] A. Raychaudhuri,  Phys. Rev. 98, 1123 (1955).
  * [95] Sanjay Mandal, P.K. Sahoo, and J.R.L. Santos, Phys. Rev. D, 102, 024057 (2020); S. Capozziello, F.S.N. Lobo, and J.O. Mimoso, Phys. Lett. B, 730, 280–283 (2014).
  * [96] V. Sahni et al., JETP Lett. 77, 201(2003).
  * [97] CDMS II Collaboration,  Science 327, 1619(2010).
  * [98] D.S. Akerib et al.,  Phys. Rev. Lett. 112, 091303(2014).
  * [99] R. Essig et al.,  Phys. Rev. Lett 109, 021301(2012).
  * [100] C.G. Bohmer et al.,  Astrop. Phys. 29, 386(2008).
  * [101] P.D. Mannheim and J.G. O’Brien,  Phys. Rev. D 85, 124020(2012).
  * [102] N. Frusciante, Physic. Rev. D  103, 044021(2021).