# The primordial black hole from running curvaton

Lei-Hua Liu1<EMAIL_ADDRESS>1Department of Physics, College of
Physics, Mechanical and Electrical Engineering, Jishou University, Jishou
416000, China

###### Abstract

In light of our previous work Liu:2019xhn , we investigate the possibility of
the formation of a primordial black hole in the second inflationary process
induced by the oscillation of curvaton. By adopting the instability of the
Mathieu equation, one could utilize the $\delta$ function to fully describe
the power spectrum. Due to the running of curvaton mass, we can simulate the
value of abundance of primordial black holes nearly covering all of the mass
ranges, in which we have given three special cases. One case could account for
the dark matter in some sense since the abundance of a primordial black hole
is about $75\%$. At late times, the relic of exponential potential could be
approximated to a constant of the order of cosmological constant dubbed as a
role of dark energy. Thus, our model could unify dark energy and dark matter
from the perspective of phenomenology. Finally, it sheds new light on
exploring Higgs physics.

##  I Introduction

The formation of a primordial black hole (PBH) in an inflationary period
consisting of preheating period offers an intriguing window for exploring the
early universe Hawking:1971ei ; Carr:1974nx ; Grindlay:1975eb . From another
perspective, the PBH could account for the formation of dark matter (DM)
according to the abundance of PBH Ivanov:1994pa ; Carr:2016drx ;
Gaggero:2016dpq ; Inomata:2017okj ; Kovetz:2017rvv ; Georg:2017mqk . The
explicit observation of gravitational waves is the most essential achievement
by the joint LIGO/Virgo collaboration, especially for the merge of black holes
(BHs) LIGOScientific:2016aoc ; LIGOScientific:2016emj ; LIGOScientific:2017bnn
, in which there are some merged BHs whose mass range is around $30M_{\odot}$
that cannot be achieved by the stellar evolution. However, the mass range of
PBH could cover this mass, thus the PBH is also of significance to exploring
the formation of BHs Bird:2016dcv ; Sasaki:2016jop ; Carr:2016drx . The range
of mass from $10^{-17}-10^{-15}M_{\odot}$ and $10^{-15}-10^{-13}M_{\odot}$
could almost explain the origin of DM since these scales cannot be constrained
by the observations. PBH is usually formed due to the gravitational collapse
of the overdense region. As a consequence, it will lead to the large amplitude
of curvature perturbation at some certain scales, which can be realized by
tuning the background dynamics of quantum field in inflationary universe
Garcia-Bellido:1996mdl ; Garcia-Bellido:2017mdw ; Domcke:2017fix ;
Kannike:2017bxn ; Carr:2017edp ; Ballesteros:2017fsr ; Hertzberg:2017dkh ;
Franciolini:2018vbk ; Kohri:2018qtx ; Ozsoy:2018flq ; Biagetti:2018pjj . Due
to its specification, PBH could leave the imprints for the observation, $\it
e.g.$ the seeds for the formation of galaxy, the evaporation of PBH could
interestingly interpret the observations of point-like gamma-ray sources
Belotsky:2014kca ; Khlopov:2008qy .

There are several ways for realizing the mechanism of enhancing the power
spectrum at small scales. One effective way named by ultra-slow-roll inflation
Martin:2012pe ; Motohashi:2014ppa , especially for the inflection point of the
potential Germani:2017bcs  dubbed as a very economical way of achieving the
ultra-slow-roll inflation for generating the enhanced power spectrum
Motohashi:2017kbs ; Ezquiaga:2017fvi ; Ballesteros:2018wlw . Under this
framework, it is not explicitly realized a viable mechanism for keeping the
e-folding number is around $50-60$ Passaglia:2018ixg ; Sasaki:2018dmp .
Another way of enhancing the power spectrum is implemented by the non-minimal
coupling and noncanonical kinetic term Fu:2019ttf ; Fu:2019vqc ;
Dalianis:2019vit ; Lin:2020goi ; Braglia:2020eai ; Gundhi:2020zvb ;
Cheong:2019vzl , in which the potential and the non-minimally coupling
function should have some special form with fine-tuning parameters of various
models. In light of $k$ inflation Armendariz-Picon:1999hyi ; Garriga:1999vw
and $G$ inflation Kobayashi:2010cm ; Kobayashi:2011nu . In particular, there
is a new mechanism for generating the large power spectrum under the framework
of noncanonical kinetic terms at small scales is called the resonant sound
speed during inflation, being different from the standard procedure for the
preheating period, in which its new place is the resonant sound speed appears
during the early inflationary period ($\Delta N$ is around $10$) Cai:2018tuh .
Due to its simplicity and richness from the phenomenology, it can be applied
to the stochastic gravitational waves Cai:2019jah  and this mechanism can be
embedded into DBI inflation Chen:2020uhe  and curvaton-inflaton mixed
inflation Chen:2019zza . A similar mechanism was proposed for which the sound
speed approaches zero at some stage of inflation in single field inflation
Ballesteros:2018wlw ; Kamenshchik:2018sig . The amplification of curvature
perturbation could also arise from the oscillation potential Cai:2019bmk .
Even the PBH can be formed due to the peak theory Wang:2021kbh .

To investigate the formation of PBH, we need the amplification of curvature
perturbation resulting from the inflationary perturbation, in which these
broad kinds of single-field inflationary theory are highly relevant to the
shape of inflation. In order to relax this stringent condition, such kind of
curvaton model was proposed, in which the curvature perturbation has arisen
from the curvaton Enqvist:2001zp ; Moroi:2001ct ; Lyth:2001nq  and it usually
is considered as an independent field. Consequently, it could account for the
various particle, $\it i.e.$ axion field Gong:2016yyb , further it could from
the axionic PBH due to the role of axion Kawasaki:2012wr ; Ando:2018nge ;
Chen:2019zza . In light of preheating process, we consider that the curvaton
field comes via inflaton decay Traschen:1990sw ; Kofman:1994rk ;
Shtanov:1994ce ; Prokopec:1996rr ; Greene:1997ge ; Kofman:1997yn ;
Greene:1997fu . If considering the multifield framework for curvaton
Liu:2020zzv , we could also investigate the impact from the noncanonical
kinetic term. Moreover, our paper also aims to unify the DM and dark energy
(DE) in light of Liu:2019xhn . From a thorough investigation of the curvaton
field, it will experience from its generation up to the DE epoch, we will give
a full analysis of the curvaton evolution.

This paper is organized as follows. In section II we revise the running
curvaton mechanism. In section III, we will investigate the formation of PBH
in light of the Mathieu equation. In section IV, we give a simple analysis of
the late time evolution of curvaton. Finally, section V gives the conclusion
and discussions.

We work in natural units in which $c=1=\hbar$, but retain the Newton constant
$G$.

##  II The revisit of running curvaton model

According to Ref. Liu:2019xhn , we investigate the formation of PBH during
preheating period. The main feature for realizing this process is that there
is a coupling between the curvaton field and the potential of inflation.
Subsequently, one can utilize the Mathieu equation to explore the formation of
PBH during the preheating period, recently a similar mechanism for realizing
the formation of PBH has been done in light of the stability of Mathieu
equation Cai:2021wzd .

###  II.1 Action

Before discussing the formation of PBH, we first revisit the running curvaton
model. At the very beginning, there is only one field called inflaton, then
the curvaton field will appear via the decay of inflaton. Meanwhile, the
explicit coupling is between the inflationary potential and curvaton field
being different from traditional coupling (Yukawa coupling for scalar fields).
Moreover, we add an exponential potential for the curvaton field to depict the
dynamical behavior of dark energy. Therefore, a total action can be written by

$S=\int
d^{4}x\sqrt{-g}\bigg{\\{}\frac{M_{P}^{2}}{2}R-\frac{1}{2}g^{\mu\nu}\nabla_{\text{$\mu$}}\phi\nabla_{\nu}\phi-\frac{1}{2}g^{\mu\nu}\nabla_{\text{$\mu$}}\chi\nabla_{\nu}\chi-V(\phi)-\frac{g_{0}}{M_{P}^{2}}\chi^{2}V(\phi)-\lambda_{0}\exp[-\lambda_{1}\frac{\chi}{M_{P}}]\bigg{\\}}$
(1)

where $\chi$ and $\phi$ are curvaton and inflaton, respectively, $R$ presents
the Ricci scalar and $g$ is the determinant of $g_{\mu\nu}$. and $g_{0}$,
$\lambda_{0}$, and $\lambda_{1}$ are the dimensionless parameters determined
by observations in Lagrangian, one point should be emphasized that
$\lambda_{0}$ is implemented to mimic the dynamical behavior of dark energy,
thus it is of order $10^{-120}$ in Planck units. And our curvaton model is
realized in preheating period for which the inflaton will oscillate its
potential, then the energy can be transformed into the curvaton field. As
Carrion:2021yeh  showed, as the inflaton oscillates near its minimal point, it
will lead to that the various inflationary potentials can be approximated by

$V(\phi)\approx\frac{1}{2}m^{2}\phi^{2},$ (2)

where $m$ is the effective mass of inflaton determined by
$\frac{d^{2}V(\phi)}{d\phi^{2}}$. Thus, our curvaton mechanism can be embedded
into many inflationary models, in which action (1) clearly shows that the
explicit coupling term $\frac{g_{0}\chi^{2}}{M_{P}^{2}}V(\phi)$ is consisting
of many forms. Due to this explicit coupling, the effective mass of curvaton
can be easily obtained by

$m_{\rm\chi}^{2}\approx\frac{g_{0}}{M_{P}^{2}}V(\phi)$ (3)

where we have neglected the contribution of the exponential part of curvaton’s
potential since it is too tiny before the inflaton finishes the decay. This
mass parameter will play a significant role in simulating the formation of PBH
in preheating period.

As we know, the wavelength is expanding slowly compared to the Hubble radius,
thus the perturbation mode of curvaton will cease at the subhorizon. To re-
enter the horizon, the second inflationary process is necessary for fulfilling
the re-entering of the horizon. In the next subsection, we will illustrate
this process via the interaction between the curvaton and other particles.

###  II.2 The second inflation

The second inflation is highly relevant with the evolution of background (for
inflaton and curvaton). In order to obtain a viable process of enhancing the
power spectrum, we first need the interaction that is similar to action (1),

$\mathcal{L}_{\rm int}=\frac{1}{2}m_{\chi}^{2}\chi^{2}+\frac{g_{1}}{M_{\rm
Pl}}\varphi^{2}V(\chi),$ (4)

where $\varphi$ is one kind of scalar particles ($\it e.g.$ ultral light
scalar particle) and $g_{1}$ is the dimensionless paramteter.

Following the Moroi:2005np , we will divide the preheating period into two
parts: (a) The first preheating comes via the oscillation of inflaton after
inflation which is named by $\phi$D or RD1 era. (b) The second inflation is
driven by the curvaton potential energy after $\chi$D dubbed as $\chi$D or RD2
era. The condition for the occurrence of the second condition is that the
initial amplitude of curvaton is large enough. The second thing we need to
figure out for which background will be alternated by the second inflation
process. As Moroi:2005np  pointed out, the curvaton model can be classified
into three types: the large field type ($m_{\chi}>M_{\rm Pl}$), the small
field type ($m_{\chi}<M_{\rm Pl}$) and the hybrid field type (the coupling
between curvaton and inflaton). In some sense, they worked out the most
general cases of curvaton scenario. In this paper, we focus on the hybrid
field type. According to their explicit calculation, we could obtain that
$N_{e}$ (e-folding number) will be enlarged to $70$ more or less, which is
very important for our later analysis of the formation scale of PBH during the
second inflation process. As for its observables, it will decrease the
spectral index $n_{\chi}$ in the hybrid type of curvaton model. Here, we will
focus on its curvature perturbation in the second inflation.

After introducing the second inflation, the curvature perturbation induced by
the curvaton will re-enter the Hubble radius. Consequently, the enhanced
density perturbation could have longer time for collapsing into the PBH after
re-enter the horizon.

###  II.3 Instability of Mathieu equation

In this subsection, we will obtain the EOM of $\delta\varphi$ field (quantum
fluctuation of some spectator field) from action (1). Here, we define that
$\varphi(x,t)=\bar{\varphi}(t)+\delta\varphi(x,t)$ where $\bar{\varphi}(t)$
denotes the background of spectator scalar particle only depending on time and
$\delta\varphi(x,t)$ is the quantum fluctuations. After this definition,
varying with action (1) and transferring into momentum space, one could obtain

$(\delta\ddot{\varphi}_{k}+3\frac{\dot{a}}{a}\delta\dot{\varphi}_{k})+\frac{k^{2}}{a^{2}}\delta\varphi_{k}+m^{2}_{\chi}\frac{g_{1}}{M_{P}^{2}}\chi^{2}\delta\varphi_{k}=0,$
(5)

there is an extra term $\frac{k^{2}}{a^{2}}\delta\varphi_{k}^{2}$ compared
with EOM of $\bar{\varphi}(t)$, in which we will analyze the quantum
fluctuations of the spectator field. Subsequently, one could utilize the EOM
of the background of curvaton field $\chi\approx\chi_{0}a^{-2/3}\sin(mt)$ in
preheating period and rearrange of eq. (5). Then, we could obtain the
following equation,

$\delta\ddot{\tilde{\varphi}}_{k}+(\frac{k^{2}}{a^{2}})\delta\tilde{\varphi}_{k}+\frac{g_{0}m^{2}_{\chi}}{2M_{P}^{2}}\chi_{0}^{2}\delta\tilde{\varphi}_{k}-\frac{m^{2}_{\chi}g_{1}}{2M_{P}^{2}}\delta\tilde{\varphi}_{k}\chi_{0}^{2}\cos(2m_{\chi}t)=0.$
(6)

Finally, we set $z=m_{\rm\chi}t$ and
$\delta\tilde{\varphi}_{k}=a^{3/2}\delta\varphi_{k}$, eq. (6) becomes

$\delta\tilde{\varphi}_{k}^{\prime\prime}+(A_{k}-2q_{k}\cos[2z])\delta\tilde{\varphi}_{k}=0$
(7)

where the correspondence can be found as follows,

$\displaystyle A_{k}$ $\displaystyle=$
$\displaystyle\frac{k^{2}}{m_{\chi}^{2}a^{2}}+\frac{g_{1}\chi_{0}^{2}}{2M_{P}^{2}},$
(8) $\displaystyle q_{k}$ $\displaystyle=$
$\displaystyle\frac{\chi_{0}^{2}g_{1}}{4M_{P}^{2}}.$ (9)

Thus, eq. (9) can become the standard form of the Mathieu equation. There is a
key feature of Mathieu equation, in which the solution of eq. (7) is
proportional to

$\varphi_{k}\propto\exp[\mu_{k}^{(n)}z]=\exp[\mu_{k}^{(n)}m_{\rm\chi}t],$ (10)

where $\mu_{k}^{(n)}$ is Lyapunov index depicting the $n$ th instability band
of Mathieu equation, in which the first band is sufficient for investigating
the production of PBH and its corresponding formula at the first band is
denoted by $\mu_{k}=\sqrt{\frac{q_{k}^{2}}{4}-(\frac{2k}{m_{\chi}}-1)^{2}}$
following the notation of Ref. Kofman:1997yn . The resonance occurs at the
$k=\frac{m_{\chi}}{2}$ and the first band could take the maximal value of
$\mu_{k}=\frac{q_{k}}{2}$ considered in the following calculation (namely
$m_{\chi}=2k$), where $m$ is the mass of inflaton in our paper, then solution
$\chi_{k}$ will become

$\varphi_{k}\propto\exp[\frac{q_{k}m_{\chi}t}{2}]=\exp[q_{k}kt],$ (11)

which will be analyzed for the PBH formation, meanwhile we will observe that
the formation of PBH is highly relevant with the curvaton mass $m_{\chi}$.
Here, we have connected $k$ to curvaton mass $m$. Due to this instability
(11), it is a key for ensuring the formation of PBH at some certain scales
which will be considered as an enhanced part of the power spectrum of
curvaton, in which it will not impact the CMB observation Akrami:2018odb .

In this section, we have revisited a large class of curvaton scenarios called
running curvaton, in which there will be second inflationary process for
forming the PBH as the initial amplitude of curvaton is large enough. During
the second inflation, the mass of curvaton will play a significant role for
the scale of PBH formation.

##  III Formation of PBH

In this section, we will proceed to the formation of PBH in light of the
instability of the Mathieu equation. Before the detailed calculation, we first
give a simple physical picture of the formation of PBH. Recall that the
formation occurs during the preheating period, the energy of curvaton will be
transferred into other spectator fields that derive the second inflation. The
mass of curvaton will play a role in the instability band.

To depict the super Hubble scale or the sub Hubble scale, it is determined by
the comoving Hubble radius $\mathcal{H}=aH$ with
$\mathcal{H}=\frac{da}{d\tau}$ ($\tau$ is the conformal time) and
$H=\frac{da}{dt}$ ($t$ is the physical time). Combined with our analysis for
the formation of PBH during the second inflation, the modes of quantum
fluctuations of the curvaton field will re-enter the horizon
($\mathcal{H}^{-1}$ as the horizon). In our calculation, we utilize the
physical time $t$. Thus, the evolution of power spectrum will be varying with
respect to $\tilde{k}=\frac{k}{a}$ not $k$. The previous analysis has shown
that the background of inflation will alternative, especially for which the
e-folding number will be changed into $70$ more or less. Then, we will use
this number as a reasonable input to analyze our Hubble radius.

###  III.1 Power spectrum within the instability of Mathieu equation

One of the most essential ingredients for producing the PBH is enhancing the
value of the power spectrum at certain scales beyond CMB constraints. In this
paper, we will utilize the instability of the Mathieu equation to obtain the
satisfying power spectrum.

To obtain the full power spectrum in light of this instability, we should
recall the content of the power spectrum utilizing $\delta N$ formalism
Liu:2019xhn , in which its corresponding formula is
$P_{\zeta}=\frac{H_{*}^{2}}{9\pi^{2}}\frac{r_{\rm decay}^{2}}{\chi^{2}}$
($H_{*}^{2}$ denotes the value at freezing time for Hubble radius and here
$\chi$ denotes the value of curvaton starts to oscillate), in which the range
of $P_{\zeta}$ is consistent with observations as taking $0.12<r_{\rm
decay}<1$. In a single field inflationary model, the power spectrum can be
obtained by $P_{\zeta}=\frac{H^{2}_{*}}{8\pi^{2}\epsilon}$, in which
$\epsilon$ is of order unity after inflation. Comparing these formulas for the
power spectrum, it was easily concluded that they are consistent with each
other once taking $0.12<r_{\rm decay}<1$. Consequently, we can utilize
$P_{\zeta}=\frac{H^{2}_{*}}{8\pi^{2}\epsilon}$ as the first part of the full
power spectrum, which could nicely recover the observations at large scales.
Subsequently, we can follow the standard procedure to rewrite the
$P_{\zeta}=\frac{H^{2}_{*}}{8\pi^{2}\epsilon}$ in terms of
$A_{s}\big{[}\frac{\tilde{k}}{\tilde{k}_{p}}\big{]}^{n_{s}-1}$ in light of
Baumann:2009ds , where $\tilde{k}_{p}$ is a fiducial comoving momentum whose
value is around $\tilde{k}_{p}=0.05~{}\rm Mpc^{-1}$.

Here, we need to emphasize that the power spectrum of our model is related to
the energy scale $k$, namely $k=\frac{m_{\chi}}{2}$. As for the observable
power spectrum, it is highly relevant to the comoving momentum. Consequently,
we need $k=a\tilde{k}$ relate the power spectrum to the observations.

From another perspective, we need the enhanced power spectrum at some certain
scales, to simulate the generation of PBH, similar mathematical property with
Cai:2018tuh  (mechanism is different since their enhanced amplitude of
curvature perturbation occurs at the inflation), the power spectrum (mainly
from curvaton) will be experienced a exponential factor via instability (11)
at some certain scales $k_{*}$,
$P_{\zeta}=\frac{H^{2}_{*}}{8\pi^{2}\epsilon}\exp[q_{k}mt]=\frac{H^{2}_{*}}{8\pi^{2}\epsilon}\exp[2q_{k}k_{\rm
scale}t]=\frac{H^{2}_{*}}{8\pi^{2}\epsilon}\exp[2q_{k}a\tilde{k}_{\rm
scale}t]$, where we have adopted the definition of
$P_{\zeta}=k^{3}|\zeta_{k}|^{2}/(2\pi^{2})$ and $|\zeta_{k}|\propto v_{k}$.
Here, one can see that $k_{\rm scale}$ 111Here, we adopt $k_{\rm scale}$ to
distinguish the $k$ in eq. (12).  is some so-called certain scales explicitly
related to the running mass $m$ via (10), the mass only needs to be smaller
compared with the upper limits of inflationary potential from COBE
normalization (analysis will be given later). Following Kofman:1997yn , it
clearly indicates that $\Delta k\propto q^{l}$ ($l$ is the l-th band of
instability of Mathieu equation) and our case is $q\ll 1$. Consequently, the
enhanced part of the power spectrum can be parametrized by the $\delta$
functions for some $k_{\rm scale}$. Taking these two factors into account, one
can get the full formula of the power spectrum as

$\displaystyle P_{\zeta}$ $\displaystyle=$ $\displaystyle
A_{s}\bigg{[}\frac{k}{k_{p}}\bigg{]}^{n_{s}-1}\bigg{[}1+\frac{q_{k}}{2}\exp\big{(}2q_{k}k_{\rm
scale}t\big{)}\delta(k-2k_{\rm scale})\bigg{]}$ (12) $\displaystyle=$
$\displaystyle
A_{s}\bigg{[}\frac{k}{k_{p}}\bigg{]}^{n_{s}-1}\bigg{[}1+\frac{q_{\rm
k}}{2}\exp(q_{\rm k}mt)\delta(k-m_{\chi})\bigg{]}$ $\displaystyle=$
$\displaystyle
A_{s}\bigg{[}\frac{k}{k_{p}}\bigg{]}^{n_{s}-1}\bigg{[}1+\frac{q_{\rm
k}}{2}\exp(q_{\rm k}mt)\delta(a\tilde{k}-m_{\chi})\bigg{]},$

where $A_{s}=\frac{H^{2}_{*}}{8\pi^{2}\epsilon}$, $q_{k}$ is defined in eq.
(9) dubbed as the amplitude of enhanced power spectrum in some sense and $t$
is the physical time and the coefficient of an exponential factor of
$\exp[q_{k}mt]$ comes via a triangle approximation. In the following
calculation, we will use the running mass as an input manifesting the main
feature of our model. Once obtaining this full power spectrum, one can
numerically simulate its range within the observational constraints.

Figure 1: The plot of power spectrum (12): The horizontal line corresponds to
energy scale whose range is $10^{-7}\rm~{}Mpc^{-1}\leqslant k\leqslant
10^{20}~{}\rm Mpc^{-1}$. The vertical line denotes the value of $P_{\zeta}$
whose range locates from $10^{-10}$ to $1$. The black solid line is
observational constraints come via Akrami:2018odb . The blue, green, red line
correspond to the enhanced part of power spectrum as taking various mass of
inflaton. $\tilde{k}_{1}=1.57\times 10^{7}~{}\rm Mpc^{-1}$ is for the blue
point, $\tilde{k}_{2}=1.6\times 10^{10}~{}\rm Mpc^{-1}$ and
$\tilde{k}_{3}=1.7\times 10^{12}~{}\rm Mpc^{-1}$ correspond to the green
point, red point, respectively. The spectral index is set by $n_{s}=0.965$.
$q_{k}\ll 1$.

In figure 1, it indicates that the varying trend of the full power spectrum
(12), in which the black solid line corresponds to the observational value
from COBE normalization Akrami:2018odb  whose order is of $10^{-9}$. In this
figure, we give three values of $\tilde{k}$ as an illustration of forming PBH
whose values are given in figure 1. To relate the realistic energy scale, we
unify the $\rm Kpc^{-1}$ into $\rm GeV$. We take $\tilde{k}_{3}$ as an
example, in which it is straightforward for obtain $\tilde{k}_{3}\approx
10^{-26}\rm GeV$. Keeping in mind that our Universe has experienced
exponential inflation and it is continuously expanding until the present, thus
scale factor is a monotonically increasing function. Compared with reheating
and preheating periods, the expansion of inflation is much larger. Thus, we
could use $a(t_{*})=\exp(70)$ (the inflation will last nearly $N=70$) as an
approximation for obtaining the energy scale
$k_{*}=a(t_{*})\tilde{k}_{3}\approx 10^{13}~{}\rm GeV$ whose value is
consistent with the energy scale of reheating or preheating period (depending
on various models). The realistic value of $k_{*}$ will be larger than this
value since we have set this approximation.

As for the various value of $k$, it could be determined by the original
definition according to
$\mu_{k}=\sqrt{\frac{q_{k}^{2}}{4}-(\frac{2k}{m}-1)^{2}}$ as taking the
maximal value for $\mu_{k}$. Especially for the enhanced part of the power
spectrum, it reveals that its corresponding value could reach the order of
$10^{-1}$ which is sufficient for the generation of PBH. All of these
numerical simulations are done with $q_{k}\ll 1$ belonging to the narrow
resonance. As for the choice of physical time, we set $t=10^{2}$ sec as a
reasonable input since the formation of PBH occurs at the deep preheating
period, in which it is after the matter-radiation equality with $10~{}\rm
s<t<3~{}\rm min$.

Additionally, this enhanced part of the power spectrum will not impact the
observational constraints since the observational scale of CMB is
approximately located from $10^{-5}~{}\rm Mpc^{-1}$ to $10~{}\rm Mpc^{-1}$,
which leads to being consistent with observations.

###  III.2 Formation of PBH in light of instability of Mathieu equation

We will utilize the enhancement of the primordial power spectrum to
investigate the formation of PBH in our theoretical framework. It can be seen
clearly that we have large parameter spaces to construct the large amplitude
of power spectrum to collapse into PBH, in which we give some specific values
of curvaton mass and various $q_{k}$ whose value is much smaller than one.
This formation process occurs during the radiation period induced by the
oscillation of curvaton (the seoncond inflation). Meanwhile, the PBH will be
generated after the perturbations of curvaton re-entry the horizon. The mass
of PBH related to the horizon mass at the horizon re-entry with various co-
moving wavenumber $k$ is

$M(k)=\gamma\frac{4\pi}{\kappa^{2}H}\approx
M_{\odot}\bigg{(}\frac{\gamma}{0.2}\bigg{)}\bigg{(}\frac{g_{*}}{10.75}\bigg{)}^{-\frac{1}{6}}\bigg{(}\frac{k}{1.9\times
10^{6}~{}\rm Mpc^{-1}}\bigg{)},$ (13)

where $\kappa^{-1}=M_{\rm pl}=2.4\times 10^{18}~{}\rm GeV$ is the reduced
Planck mass, $H$ is evaluated at $k=aH$ ($H$ is Hubble parameter and $a$ is
scale factor), $\gamma$ is defined to the ratio of PBH mass to the horizon
mass indicating the efficiency of collapse, in which its value is
approximately set by $0.2$ from Ref. Carr:1975qj , $g_{*}$ is the degrees
freedom of energy densities at the formation of PBH. Since its formation
occurs during the preheating period, the process happens during the deep
radiation period whose value can be determined by $g_{*}=106.75$. Here, we
should emphasize the difference with the traditional curvaton mechanism, in
which curvaton dubs as an independent field leading to the decay of curvaton
after the actual preheating. However, our model is different since the
curvaton generates and practically disappears as the inflaton decays, to be
more precise, the curvaton formed as the inflaton starts to decay, and
meanwhile, the curvaton will disappear as the inflaton finishes the decay
process due to the coupling between the curvaton and inflaton. During the
whole process of preheating for inflaton, the PBH will be formed due to the
instability of curvaton. That is the reason we mentioned this process happens
during the deep radiation period.

To investigate the abundance of PBH with mass $M$, we need the mass fraction
$\beta(M)$ against the total energy at the formation of PBH. Under the
assumption of the distribution is Gaussian, it can be expressed by
Tada:2019amh ; Young:2014ana

$\beta(M)=\frac{\rho_{\rm PBH}}{\rho_{\rm total}}=\frac{1}{2}\rm
erfc\bigg{(}\frac{\delta_{c}}{\sqrt{2\sigma^{2}(M)}}\bigg{)},$ (14)

where $\rm erfc$ is the complementary error function and $\delta_{c}$ denotes
the threshold of perturbation at the formation of PBH whose value is around
$0.4$ in light of Musco:2012au ; Harada:2013epa , in which Ref. Martin:2019nuw
; Hidalgo:2017dfp  also proposed the value of $\delta_{c}$, even Ref.
Escriva:2019phb ; Escriva:2020tak  discussed the value of $\delta_{c}$ is also
relavant with compaction function leading to a tiny deviation with $0.4$.
$\delta(M)$ is the variance of the density perturbation with mass $M$ for PBH
which can be associated with power spectrum,

$\sigma^{2}(M(k))=\int d\ln
qW^{2}(qk^{-1})\frac{16}{81}(qk^{-1})^{4}P_{\zeta}(q),$ (15)

where $W(x)=\exp(-x^{2}/2)$ (Gaussian window function). Once having these two
physical quantities, we can define the abundance of PBH, namely the fraction
of PBH to the total DM,

$f_{\rm PBH}=\frac{\Omega_{\rm PBH}}{\Omega_{\rm DM}}=2.0\times
10^{8}\bigg{(}\frac{\gamma}{0.2}\bigg{)}^{1/2}\bigg{(}\frac{10.75}{g_{*}}\bigg{)}^{1/4}\bigg{(}\frac{M_{\odot}}{M}\bigg{)}^{1/2}\beta(M),$
(16)

where $\Omega_{\rm DM}$ is the current energy density of DM from Planck 2018
results giving its value is around $\Omega_{\rm DM}h^{2}\approx 0.12$. In
light of basic estimation for the abundance of PBH, the amount of power
spectrum of $P_{\zeta}$ should reach the order of $10^{-2}$ to produce the
sizable PBH abundance on small scales, in which we have given our numerical
simulation from figure 1.

Figure 2: The plot of the abundance of PBH (16): The horizontal line
corresponds to the ratio of $M_{\rm PBH}$ to $M_{\rm\odot}$ whose range is
from $10^{-19}$ to $10^{4}$, which could cover the whole mass range of PBH.
The vertical line is the abundance of PBH, in which its corresponding range is
from $10^{-5}$ to $1$. The red, green, and blue dashed lines correspond to the
cases as depicted in figure 1 with various masses of the inflaton. The blue
allowed region comes via the ultrashort-timescale microlensing events of OGLE
data Niikura:2019kqi . The other shadow part denote the current allowed
observations, the extragalactic gamma-rays of PBH evaporation (EG$\gamma$)
Carr:2009jm , the white dwarfs explosion (WD) Graham:2015apa , the galactic
center $511~{}keV$ gamma-ray line (INTEGRAL) Laha:2019ssq  (Ref.
Dasgupta:2019cae  discussed the allowing maximum rotation can significantly
improve and extend the constraints from $511~{}keV$ in higher mass windows),
the microlensing events with Subaru HSC (Subaru HSC) Acharyya:2019nwy , with
the Kepler satellite (Kepler) Griest:2013esa , the EROS/MACHO (EROS/MACHO)
EROS-2:2006ryy , and the constraints from CMB (CMB) Ali-Haimoud:2016mbv .

Figure 2 indicates the abundance of PBH to the total energy density of DM in
terms of the running curvaton scenario. There are three cases corresponding to
the various masses as showing in figure 1, in which their mass are
$10^{-12}M_{\odot}$, $10^{-8}M_{\odot}$ and $10^{0}M_{\odot}$, respectively,
in which their the corresponding mass is obtained by the definition of
$\mu_{k}$ as taking its maximal value that is namely $2k=m$ as setting
$m_{1}$, $m_{2}$, $m_{3}$, $\rm e.t.c$. One these three cases, case one (the
red point) could result in the main component of DM whose percentage is around
$76.5\%$, namely, it could account for the DM in some sense. As for case 2
(green point) and case 3 (blue point), they are far from the observational
constraints, especially compared with Kepler Griest:2013esa  and Subaru HSC
Acharyya:2019nwy . Case 3 has a similar situation compared with CMB Ali-
Haimoud:2016mbv . However, Ref. Hutsi:2020sol  claimed that the mass range of
PBH is from $2$ to $400$ solar mass, which could only account for the $0.2\%$
DM, meanwhile Carr:2020gox  also discussed that the less contribution of DM
will play an important role and supply a significant probe for the early
Universe.

Here, we only vary with mass of inflaton, in which the potential of inflaton
can be constrained as follows Liu:2018hno ,

$\frac{V(\phi)}{M_{\rm P}^{4}}\approx 3.0\times
10^{-10}\bigg{(}\frac{r_{*}}{5\times
10^{-3}}\bigg{)}\bigg{(}\frac{A_{s}}{2.1\times 10^{-9}}\bigg{)},$ (17)

where $r_{*}$ is the tensor-to-scalar ratio whose value is less than $0.06$
and $A_{s}$ is the amplitude of the power spectrum of curvature perturbation
whose value is around $2.1\times 10^{-9}$. The constraint comes from COBE
normalization, which tells the range of inflationary potential during
inflation. After inflation, the Universe undergoes the preheating period, and
the value of inflationary potential is smaller than $3.0\times 10^{-9}$ since
the value of the inflaton field will be decreased by transferring the energy
into other fundamental particles (i.e. Higgs particles, and other Fermions).
Thus, it gives us lots of freedom to simulate the range of value for inflaton
mass, which we have shown in figure 2 and 1. Consequently, our mechanism will
nicely recover the whole mass range of PBH. In some sense, it could account
for the origin of DM. In our numerical simulation, the parameters we adapted
belong to the narrow resonance, namely $q_{k}\ll 1$. In figure 3, we have
plotted the range $q_{k}$, in which we have set $M_{\rm P}=1$ and $a=1$ (scale
factor), and the other parameters can be explicitly shown in this figure. From
figure 3, it indicates that the range of $g_{0}$ cannot be large in order to
obtain $q_{k}\ll 1$, whose range is compatible with Torres-Lomas:2014bua . As
for the broad narrow resonance, it will also lead to the sufficient production
of PBH. However, the over-production of PBH is not a generic feature even in
broad resonance, namely $q\gg 1$ Torres-Lomas:2014bua , which is determined by
the coupling constant between the inflaton field and target field. By
considering our model, the curvaton comes via the decay of inflaton, thus the
formation of PBH is an inevitable process whose value of PBH abundance is
sufficient analysis in Ref. Suyama:2004mz , in which the criteria of PBH
formation is the duration time of preheating, thus the power spectrum is
mainly characterized by the scale, not the key feature of instability of
Mathieu equation describing the preheating period. To depict this feature, we
characterize the power spectrum by the $\delta$ function highly relevant to
the scale on small scales. For the large scales, it is almost the Gaussian for
power spectrum (nearly scale-invariant).

Figure 3: Contour Plot of $q_{k}$ via eq. (9): The horizontal line corresponds
to the amplitude of inflaton whose range is $0\leqslant phi_{0}\leqslant 20$
determined by the e-folding number is around $60$. The vertical line denotes
the value of $g_{0}$ whose range locates from $0$ to $0.0001$, we have set
$M_{\rm P}=1$ and $a=1$. The right panel shows that the value of $q_{k}$
matching its corresponding color.

In this section, we have thoroughly investigated the formation of PBH during
the preheating period. The key ingredient for generating the PBH is utilizing
the instability of the Mathieu equation characterized by the $\delta$
functions in the power spectrum at small scales. In the next section, we will
proceed with the study of the equation of state $EoS$ to investigate the dark
energy epoch.

##  IV The late time evolution of curvaton field

The curvaton mechanism occurs during the preheating period, meanwhile, it
comes via the decay of inflaton. From action (1), it indicates that the main
contribution will be disappeared as the inflaton finish the decay process. To
clarify this situation, we will make a plot of potential consisting of the
inflaton and curvaton, however, the curvaton comes via the transferring of
energy of inflaton, meanwhile, there are lots of uncertainties about
preheating mechanism. We will only show varying values of inflaton and
curvaton, in which the inflaton field will be decreased and the curvaton will
be enhanced.

Figure 4: Contour Plot of curvaton. It shows that the potential of curvaton
field, in which $V(\chi)=\frac{g_{0}\chi^{2}}{M_{\rm
P}}V(\phi)+\lambda_{0}\exp(-\lambda_{1}\frac{\chi}{M_{\rm P}})$ with
$V(\phi)=\frac{1}{2}m^{2}\phi^{2}$ as an instance. we have set $M_{\rm P}=1$,
$g_{0}=0.001$ and $\lambda_{1}$ has been obsorbed in the curvaton field.

In figure 4, we have shown the potential of curvaton during preheating which
just is considered an illustration of the range of curvaton. The field value
of curvaton comes via the inflaton decay. Consequently, its contour plot is
highly probably different from figure 4, however, it will not change the range
of field for curvaton. We could see that the range of curvaton is around
$10^{-15}$ in Planck units as taking proper parameters, which means that its
energy density will be dominant after inflaton decay but much more tiny
compared with the energy scope of inflationary potential (matching with
curvaton assumption). As for the curvaton decay, it will nearly disappear as
the completion of inflaton decay only keeping the relic of the exponential
potential of $curvaton$, namely, it is almost impossible to decay into other
particles determined by its effective mass $m_{\rm
eff}^{2}=\frac{d^{2}V(\chi)}{d\chi^{2}}$, since Ref. Liu:2019xhn  has
indicated that the order of effective mass for curvaton is of order
cosmological constant that means that the effective mass of curvaton is much
smaller compared with all kinds of fundamental particles including Higgs
particles, fermions, and gauge field particles, $\rm e.t.c.$, in which the new
place of our model is that the lifetime of curvaton is almost the same with
inflaton since the main mass part of curvaton is proportional to the
inflationary potential. Simple and explicit analysis for the curvaton
potential is taken into account, in which it will be naturally approached to a
constant being of the order of cosmological constant dubbed as a role of dark
energy.

From this section, we give a simple analysis that shows that the range of
curvaton is less than $10^{-15}$ in Planck units. In this scope, the potential
of curvaton is dominant during preheating but much more tiny compared with the
inflationary energy scale matching the assumption of curvaton. At late times,
the curvaton field will also decay into other fundamental particles. However,
there is a relic of the exponential potential of curvaton playing a role in
cosmological constant, which can be dubbed as the dark energy dominating the
current epoch. In some sense, this curvaton mechanism adding the formation of
PBH during preheating could unify the dark energy and dark matter.

##  V Conclusion and discussions

In this paper, we have investigated the formation of PBH during the preheating
period in light of the instability of the Mathieu equation. Being different
from previous work Martin:2019nuw ; Torres-Lomas:2014bua ; Suyama:2004mz , we
will implement the running mass of curvaton which is proportional to (2) to
investigate the formation of PBH during the deeply preheating period (induced
by the second inflation). Thus, the mass of curvaton will be varying from COBE
normalization to the DE scale determined by $V_{\rm
eff}=\frac{d^{2}V}{d\chi^{2}}$ as the inflaton finishes the decay. Due to its
huge range, we can relate the exponential growth of curvature perturbation to
the mass scale which is the instability band corresponding to the $\tilde{k}$.
Then, we can adapt the property in light of the property of $\delta$
functions, the other energy scale is almost scale-invariant as shown in eq.
(12). Thus, the simple analytical power spectrum perfectly agreed with
observations whose detail can be seen in figure 1.

Once obtained this key result of this paper, we can numerically simulate the
abundance of PBH among the DM. Figure 2 clearly indicates the value of $f_{\rm
PBH}$ with some specific values of mass corresponding to the certain
$\tilde{k}$. Here, we only use the property of the Lyapunov index to find the
corresponding various values of $k$ could simulate different values of $f_{\rm
PBH}$ as choosing proper parameters. Especially for the case $1$ of figure 2,
it shows that the value of abundance of PBH could reach the $75\%$ of DM which
could account for the DM in some sense. During the preheating period, we have
shown that the potential is of an order of $10^{-15}$ dominating the main
content of the Universe and meanwhile being consistent with the assumption of
curvaton. At late times, as the inflaton almost completes the process of decay
(curvaton disappears), there is a relic of curvaton exponential potential that
will be dubbed as the role of dark energy as illustrated in section IV. Thus,
our model could unify dark energy and dark matter.

Our model is highly relevant to the coupling structure. Consequently, we could
use this similar mechanism to explore the possibility of the formation of the
Higgs field during preheating, and the mass of the Higgs field is stringently
constrained by the observation and also its rich decay channels. It sheds a
light on exploring Higgs physics. Finally, we should emphasize that our non-
Guaassianity parameter can be extended into a new one associated with $w$
(EoS) and $r_{\rm decay}$ Liu:2020zlr .

## Acknowledgements

LH is funded by the Hunan Provincial Department of Education, NO. 19B464, and
the National Natural Science Foundation of China Grants (NSFC) NO. 12165009.
We are grateful for the discussions with Diego Cruces for the value of
$\delta$ and the numerical simulation of Wu-Long Xu.

## References

## References

  * (1) L. H. Liu and W. L. Xu, Chin. Phys. C 44, no.8, 085103 (2020) doi:10.1088/1674-1137/44/8/085103 [arXiv:1911.10542 [astro-ph.CO]].
  * (2) S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971)
  * (3) B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399-415 (1974)
  * (4) J. E. Grindlay, H. F. Helmken, R. H. Brown, J. Davis and L. R. Allen, Astrophys. J. 201, 82-89 (1975) doi:10.1086/153861
  * (5) P. Ivanov, P. Naselsky and I. Novikov, Phys. Rev. D 50, 7173-7178 (1994) doi:10.1103/PhysRevD.50.7173
  * (6) B. Carr, F. Kuhnel and M. Sandstad, Phys. Rev. D 94, no.8, 083504 (2016) doi:10.1103/PhysRevD.94.083504 [arXiv:1607.06077 [astro-ph.CO]].
  * (7) D. Gaggero, G. Bertone, F. Calore, R. M. T. Connors, M. Lovell, S. Markoff and E. Storm, Phys. Rev. Lett. 118, no.24, 241101 (2017) doi:10.1103/PhysRevLett.118.241101 [arXiv:1612.00457 [astro-ph.HE]].
  * (8) K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada and T. T. Yanagida, Phys. Rev. D 96, no.4, 043504 (2017) doi:10.1103/PhysRevD.96.043504 [arXiv:1701.02544 [astro-ph.CO]].
  * (9) E. D. Kovetz, Phys. Rev. Lett. 119, no.13, 131301 (2017) doi:10.1103/PhysRevLett.119.131301 [arXiv:1705.09182 [astro-ph.CO]].
  * (10) J. Georg and S. Watson, JHEP 09, 138 (2017) doi:10.1007/JHEP09(2017)138 [arXiv:1703.04825 [astro-ph.CO]].
  * (11) B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116, no.6, 061102 (2016) doi:10.1103/PhysRevLett.116.061102 [arXiv:1602.03837 [gr-qc]].
  * (12) B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116, no.13, 131103 (2016) doi:10.1103/PhysRevLett.116.131103 [arXiv:1602.03838 [gr-qc]].
  * (13) B. P. Abbott et al. [LIGO Scientific and VIRGO], Phys. Rev. Lett. 118, no.22, 221101 (2017) [erratum: Phys. Rev. Lett. 121, no.12, 129901 (2018)] doi:10.1103/PhysRevLett.118.221101 [arXiv:1706.01812 [gr-qc]].
  * (14) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli and A. G. Riess, Phys. Rev. Lett. 116, no.20, 201301 (2016) doi:10.1103/PhysRevLett.116.201301 [arXiv:1603.00464 [astro-ph.CO]].
  * (15) M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Phys. Rev. Lett. 117, no.6, 061101 (2016) [erratum: Phys. Rev. Lett. 121, no.5, 059901 (2018)] doi:10.1103/PhysRevLett.117.061101 [arXiv:1603.08338 [astro-ph.CO]].
  * (16) J. Garcia-Bellido, A. D. Linde and D. Wands, Phys. Rev. D 54, 6040-6058 (1996) doi:10.1103/PhysRevD.54.6040 [arXiv:astro-ph/9605094 [astro-ph]].
  * (17) J. Garcia-Bellido and E. Ruiz Morales, Phys. Dark Univ. 18, 47-54 (2017) doi:10.1016/j.dark.2017.09.007 [arXiv:1702.03901 [astro-ph.CO]].
  * (18) V. Domcke, F. Muia, M. Pieroni and L. T. Witkowski, JCAP 07, 048 (2017) doi:10.1088/1475-7516/2017/07/048 [arXiv:1704.03464 [astro-ph.CO]].
  * (19) K. Kannike, L. Marzola, M. Raidal and H. Veermäe, JCAP 09, 020 (2017) doi:10.1088/1475-7516/2017/09/020 [arXiv:1705.06225 [astro-ph.CO]].
  * (20) B. Carr, T. Tenkanen and V. Vaskonen, Phys. Rev. D 96, no.6, 063507 (2017) doi:10.1103/PhysRevD.96.063507 [arXiv:1706.03746 [astro-ph.CO]].
  * (21) G. Ballesteros and M. Taoso, Phys. Rev. D 97, no.2, 023501 (2018) doi:10.1103/PhysRevD.97.023501 [arXiv:1709.05565 [hep-ph]].
  * (22) M. P. Hertzberg and M. Yamada, Phys. Rev. D 97, no.8, 083509 (2018) doi:10.1103/PhysRevD.97.083509 [arXiv:1712.09750 [astro-ph.CO]].
  * (23) G. Franciolini, A. Kehagias, S. Matarrese and A. Riotto, JCAP 03, 016 (2018) doi:10.1088/1475-7516/2018/03/016 [arXiv:1801.09415 [astro-ph.CO]].
  * (24) K. Kohri and T. Terada, Class. Quant. Grav. 35, no.23, 235017 (2018) doi:10.1088/1361-6382/aaea18 [arXiv:1802.06785 [astro-ph.CO]].
  * (25) O. Özsoy, S. Parameswaran, G. Tasinato and I. Zavala, JCAP 07, 005 (2018) doi:10.1088/1475-7516/2018/07/005 [arXiv:1803.07626 [hep-th]].
  * (26) M. Biagetti, G. Franciolini, A. Kehagias and A. Riotto, JCAP 07, 032 (2018) doi:10.1088/1475-7516/2018/07/032 [arXiv:1804.07124 [astro-ph.CO]].
  * (27) K. M. Belotsky, A. D. Dmitriev, E. A. Esipova, V. A. Gani, A. V. Grobov, M. Y. Khlopov, A. A. Kirillov, S. G. Rubin and I. V. Svadkovsky, Mod. Phys. Lett. A 29, no.37, 1440005 (2014) doi:10.1142/S0217732314400057 [arXiv:1410.0203 [astro-ph.CO]].
  * (28) M. Y. Khlopov, Res. Astron. Astrophys. 10, 495-528 (2010) doi:10.1088/1674-4527/10/6/001 [arXiv:0801.0116 [astro-ph]].
  * (29) J. Martin, H. Motohashi and T. Suyama, Phys. Rev. D 87, no.2, 023514 (2013) doi:10.1103/PhysRevD.87.023514 [arXiv:1211.0083 [astro-ph.CO]].
  * (30) H. Motohashi, A. A. Starobinsky and J. Yokoyama, JCAP 09, 018 (2015) doi:10.1088/1475-7516/2015/09/018 [arXiv:1411.5021 [astro-ph.CO]].
  * (31) C. Germani and T. Prokopec, Phys. Dark Univ. 18, 6-10 (2017) doi:10.1016/j.dark.2017.09.001 [arXiv:1706.04226 [astro-ph.CO]].
  * (32) H. Motohashi and W. Hu, Phys. Rev. D 96, no.6, 063503 (2017) doi:10.1103/PhysRevD.96.063503 [arXiv:1706.06784 [astro-ph.CO]].
  * (33) J. M. Ezquiaga, J. Garcia-Bellido and E. Ruiz Morales, Phys. Lett. B 776, 345-349 (2018) doi:10.1016/j.physletb.2017.11.039 [arXiv:1705.04861 [astro-ph.CO]].
  * (34) G. Ballesteros, J. Beltran Jimenez and M. Pieroni, JCAP 06, 016 (2019) doi:10.1088/1475-7516/2019/06/016 [arXiv:1811.03065 [astro-ph.CO]].
  * (35) S. Passaglia, W. Hu and H. Motohashi, Phys. Rev. D 99, no.4, 043536 (2019) doi:10.1103/PhysRevD.99.043536 [arXiv:1812.08243 [astro-ph.CO]].
  * (36) M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Class. Quant. Grav. 35, no.6, 063001 (2018) doi:10.1088/1361-6382/aaa7b4 [arXiv:1801.05235 [astro-ph.CO]].
  * (37) C. Fu, P. Wu and H. Yu, Phys. Rev. D 100, no.6, 063532 (2019) doi:10.1103/PhysRevD.100.063532 [arXiv:1907.05042 [astro-ph.CO]].
  * (38) C. Fu, P. Wu and H. Yu, Phys. Rev. D 101, no.2, 023529 (2020) doi:10.1103/PhysRevD.101.023529 [arXiv:1912.05927 [astro-ph.CO]].
  * (39) I. Dalianis, S. Karydas and E. Papantonopoulos, JCAP 06, 040 (2020) doi:10.1088/1475-7516/2020/06/040 [arXiv:1910.00622 [astro-ph.CO]].
  * (40) J. Lin, Q. Gao, Y. Gong, Y. Lu, C. Zhang and F. Zhang, Phys. Rev. D 101, no.10, 103515 (2020) doi:10.1103/PhysRevD.101.103515 [arXiv:2001.05909 [gr-qc]].
  * (41) M. Braglia, D. K. Hazra, F. Finelli, G. F. Smoot, L. Sriramkumar and A. A. Starobinsky, JCAP 08, 001 (2020) doi:10.1088/1475-7516/2020/08/001 [arXiv:2005.02895 [astro-ph.CO]].
  * (42) A. Gundhi and C. F. Steinwachs, Eur. Phys. J. C 81, no.5, 460 (2021) doi:10.1140/epjc/s10052-021-09225-2 [arXiv:2011.09485 [hep-th]].
  * (43) D. Y. Cheong, S. M. Lee and S. C. Park, JCAP 01, 032 (2021) doi:10.1088/1475-7516/2021/01/032 [arXiv:1912.12032 [hep-ph]].
  * (44) C. Armendariz-Picon, T. Damour and V. F. Mukhanov, Phys. Lett. B 458, 209-218 (1999) doi:10.1016/S0370-2693(99)00603-6 [arXiv:hep-th/9904075 [hep-th]].
  * (45) J. Garriga and V. F. Mukhanov, Phys. Lett. B 458, 219-225 (1999) doi:10.1016/S0370-2693(99)00602-4 [arXiv:hep-th/9904176 [hep-th]].
  * (46) T. Kobayashi, M. Yamaguchi and J. Yokoyama, Phys. Rev. Lett. 105, 231302 (2010) doi:10.1103/PhysRevLett.105.231302 [arXiv:1008.0603 [hep-th]].
  * (47) T. Kobayashi, M. Yamaguchi and J. Yokoyama, Prog. Theor. Phys. 126, 511-529 (2011) doi:10.1143/PTP.126.511 [arXiv:1105.5723 [hep-th]].
  * (48) Y. F. Cai, X. Tong, D. G. Wang and S. F. Yan, Phys. Rev. Lett. 121, no.8, 081306 (2018) doi:10.1103/PhysRevLett.121.081306 [arXiv:1805.03639 [astro-ph.CO]].
  * (49) Y. F. Cai, C. Chen, X. Tong, D. G. Wang and S. F. Yan, Phys. Rev. D 100, no.4, 043518 (2019) doi:10.1103/PhysRevD.100.043518 [arXiv:1902.08187 [astro-ph.CO]].
  * (50) C. Chen, X. H. Ma and Y. F. Cai, Phys. Rev. D 102, no.6, 063526 (2020) doi:10.1103/PhysRevD.102.063526 [arXiv:2003.03821 [astro-ph.CO]].
  * (51) C. Chen and Y. F. Cai, JCAP 10, 068 (2019) doi:10.1088/1475-7516/2019/10/068 [arXiv:1908.03942 [astro-ph.CO]].
  * (52) A. Y. Kamenshchik, A. Tronconi, T. Vardanyan and G. Venturi, Phys. Lett. B 791, 201-205 (2019) doi:10.1016/j.physletb.2019.02.036 [arXiv:1812.02547 [gr-qc]].
  * (53) R. G. Cai, Z. K. Guo, J. Liu, L. Liu and X. Y. Yang, JCAP 06, 013 (2020) doi:10.1088/1475-7516/2020/06/013 [arXiv:1912.10437 [astro-ph.CO]].
  * (54) Q. Wang, Y. C. Liu, B. Y. Su and N. Li, Phys. Rev. D 104, no.8, 083546 (2021) doi:10.1103/PhysRevD.104.083546 [arXiv:2111.10028 [astro-ph.CO]].
  * (55) K. Enqvist and M. S. Sloth, “Adiabatic CMB perturbations in pre - big bang string cosmology,” Nucl. Phys. B 626 (2002) 395 doi:10.1016/S0550-3213(02)00043-3 [hep-ph/0109214].
  * (56) D. H. Lyth and D. Wands, “Generating the curvature perturbation without an inflaton,” Phys. Lett. B 524 (2002) 5 doi:10.1016/S0370-2693(01)01366-1 [hep-ph/0110002].
  * (57) T. Moroi and T. Takahashi, “Effects of cosmological moduli fields on cosmic microwave background,” Phys. Lett. B 522 (2001) 215 Erratum: [Phys. Lett. B 539 (2002) 303] doi:10.1016/S0370-2693(02)02070-1, 10.1016/S0370-2693(01)01295-3 [hep-ph/0110096].
  * (58) J. O. Gong, N. Kitajima and T. Terada, JCAP 1703 (2017) 053 doi:10.1088/1475-7516/2017/03/053 [arXiv:1611.08975 [hep-ph]].
  * (59) M. Kawasaki, N. Kitajima and T. T. Yanagida, Phys. Rev. D 87 (2013) no.6, 063519 doi:10.1103/PhysRevD.87.063519 [arXiv:1207.2550 [hep-ph]].
  * (60) K. Ando, M. Kawasaki and H. Nakatsuka, Phys. Rev. D 98 (2018) no.8, 083508 doi:10.1103/PhysRevD.98.083508 [arXiv:1805.07757 [astro-ph.CO]].
  * (61) J. H. Traschen and R. H. Brandenberger, “Particle Production During Out-of-equilibrium Phase Transitions,” Phys. Rev. D 42 (1990) 2491. doi:10.1103/PhysRevD.42.2491
  * (62) L. Kofman, A. D. Linde and A. A. Starobinsky, “Reheating after inflation,” Phys. Rev. Lett. 73 (1994) 3195 doi:10.1103/PhysRevLett.73.3195 [hep-th/9405187].
  * (63) Y. Shtanov, J. H. Traschen and R. H. Brandenberger, “Universe reheating after inflation,” Phys. Rev. D 51 (1995) 5438 doi:10.1103/PhysRevD.51.5438 [hep-ph/9407247].
  * (64) T. Prokopec and T. G. Roos, “Lattice study of classical inflaton decay,” Phys. Rev. D 55 (1997) 3768 doi:10.1103/PhysRevD.55.3768 [hep-ph/9610400].
  * (65) B. R. Greene, T. Prokopec and T. G. Roos, “Inflaton decay and heavy particle production with negative coupling,” Phys. Rev. D 56 (1997) 6484 doi:10.1103/PhysRevD.56.6484 [hep-ph/9705357].
  * (66) L. Kofman, A. D. Linde and A. A. Starobinsky, “Towards the theory of reheating after inflation,” Phys. Rev. D 56 (1997) 3258 doi:10.1103/PhysRevD.56.3258 [hep-ph/9704452].
  * (67) P. B. Greene, L. Kofman, A. D. Linde and A. A. Starobinsky, “Structure of resonance in preheating after inflation,” Phys. Rev. D 56 (1997) 6175 doi:10.1103/PhysRevD.56.6175 [hep-ph/9705347].
  * (68) L. H. Liu and T. Prokopec, JCAP 06, 033 (2021) doi:10.1088/1475-7516/2021/06/033 [arXiv:2005.11069 [astro-ph.CO]].
  * (69) R. G. Cai, C. Chen and C. Fu, Phys. Rev. D 104, no.8, 083537 (2021) doi:10.1103/PhysRevD.104.083537 [arXiv:2108.03422 [astro-ph.CO]].
  * (70) K. Carrion, J. C. Hidalgo, A. Montiel and L. E. Padilla, JCAP 07, 001 (2021) doi:10.1088/1475-7516/2021/07/001 [arXiv:2101.02156 [astro-ph.CO]].
  * (71) T. Moroi and T. Takahashi, Phys. Rev. D 72, 023505 (2005) doi:10.1103/PhysRevD.72.023505 [arXiv:astro-ph/0505339 [astro-ph]].
  * (72) Y. Akrami et al. [Planck Collaboration], “Planck 2018 results. X. Constraints on inflation,” arXiv:1807.06211 [astro-ph.CO].
  * (73) D. Baumann, doi:10.1142/9789814327183_0010 [arXiv:0907.5424 [hep-th]].
  * (74) B. J. Carr, Astrophys. J. 201, 1-19 (1975) doi:10.1086/153853
  * (75) Y. Tada and S. Yokoyama, Phys. Rev. D 100, no.2, 023537 (2019) doi:10.1103/PhysRevD.100.023537 [arXiv:1904.10298 [astro-ph.CO]].
  * (76) S. Young, C. T. Byrnes and M. Sasaki, JCAP 07, 045 (2014) doi:10.1088/1475-7516/2014/07/045 [arXiv:1405.7023 [gr-qc]].
  * (77) I. Musco and J. C. Miller, Class. Quant. Grav. 30, 145009 (2013) doi:10.1088/0264-9381/30/14/145009 [arXiv:1201.2379 [gr-qc]].
  * (78) T. Harada, C. M. Yoo and K. Kohri, Phys. Rev. D 88, no.8, 084051 (2013) [erratum: Phys. Rev. D 89, no.2, 029903 (2014)] doi:10.1103/PhysRevD.88.084051 [arXiv:1309.4201 [astro-ph.CO]].
  * (79) J. Martin, T. Papanikolaou and V. Vennin, JCAP 01, 024 (2020) doi:10.1088/1475-7516/2020/01/024 [arXiv:1907.04236 [astro-ph.CO]].
  * (80) J. C. Hidalgo, J. De Santiago, G. German, N. Barbosa-Cendejas and W. Ruiz-Luna, Phys. Rev. D 96, no.6, 063504 (2017) doi:10.1103/PhysRevD.96.063504 [arXiv:1705.02308 [astro-ph.CO]].
  * (81) A. Escrivà, C. Germani and R. K. Sheth, Phys. Rev. D 101, no.4, 044022 (2020) doi:10.1103/PhysRevD.101.044022 [arXiv:1907.13311 [gr-qc]].
  * (82) A. Escrivà, C. Germani and R. K. Sheth, JCAP 01, 030 (2021) doi:10.1088/1475-7516/2021/01/030 [arXiv:2007.05564 [gr-qc]].
  * (83) H. Niikura, M. Takada, S. Yokoyama, T. Sumi and S. Masaki, Phys. Rev. D 99, no.8, 083503 (2019) doi:10.1103/PhysRevD.99.083503 [arXiv:1901.07120 [astro-ph.CO]].
  * (84) B. J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, Phys. Rev. D 81, 104019 (2010) doi:10.1103/PhysRevD.81.104019 [arXiv:0912.5297 [astro-ph.CO]].
  * (85) P. W. Graham, S. Rajendran and J. Varela, novae,” Phys. Rev. D 92, no.6, 063007 (2015) doi:10.1103/PhysRevD.92.063007 [arXiv:1505.04444 [hep-ph]].
  * (86) R. Laha, Phys. Rev. Lett. 123, no.25, 251101 (2019) doi:10.1103/PhysRevLett.123.251101 [arXiv:1906.09994 [astro-ph.HE]].
  * (87) B. Dasgupta, R. Laha and A. Ray, Phys. Rev. Lett. 125, no.10, 101101 (2020) doi:10.1103/PhysRevLett.125.101101 [arXiv:1912.01014 [hep-ph]].
  * (88) A. Acharyya, I. Agudo, E. O. Angüner, R. Alfaro, J. Alfaro, C. Alispach, R. Aloisio, R. Alves Batista, J. P. Amans and L. Amati, et al. Astropart. Phys. 111, 35-53 (2019) doi:10.1016/j.astropartphys.2019.04.001 [arXiv:1904.01426 [astro-ph.IM]].
  * (89) K. Griest, A. M. Cieplak and M. J. Lehner, Phys. Rev. Lett. 111, no.18, 181302 (2013) doi:10.1103/PhysRevLett.111.181302
  * (90) P. Tisserand et al. [EROS-2], Astron. Astrophys. 469, 387-404 (2007) doi:10.1051/0004-6361:20066017 [arXiv:astro-ph/0607207 [astro-ph]].
  * (91) Y. Ali-Haïmoud and M. Kamionkowski, Phys. Rev. D 95, no.4, 043534 (2017) doi:10.1103/PhysRevD.95.043534 [arXiv:1612.05644 [astro-ph.CO]].
  * (92) G. Hütsi, M. Raidal, V. Vaskonen and H. Veermäe, JCAP 03, 068 (2021) doi:10.1088/1475-7516/2021/03/068 [arXiv:2012.02786 [astro-ph.CO]].
  * (93) B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, [arXiv:2002.12778 [astro-ph.CO]].
  * (94) L. H. Liu, T. Prokopec and A. A. Starobinsky, Phys. Rev. D 98, no.4, 043505 (2018) doi:10.1103/PhysRevD.98.043505 [arXiv:1806.05407 [gr-qc]].
  * (95) E. Torres-Lomas, J. C. Hidalgo, K. A. Malik and L. A. Ureña-López, Phys. Rev. D 89, no.8, 083008 (2014) doi:10.1103/PhysRevD.89.083008 [arXiv:1401.6960 [astro-ph.CO]].
  * (96) T. Suyama, T. Tanaka, B. Bassett and H. Kudoh, Phys. Rev. D 71, 063507 (2005) doi:10.1103/PhysRevD.71.063507 [arXiv:hep-ph/0410247 [hep-ph]].
  * (97) L. H. Liu, B. Liang, Y. C. Zhou, X. D. Liu, W. L. Xu and A. C. Li, Phys. Rev. D 103, no.6, 063515 (2021) doi:10.1103/PhysRevD.103.063515 [arXiv:2007.08278 [astro-ph.CO]].


# Accelerating expansion of the universe in modified symmetric teleparallel
gravity

Raja Solanki 0000-0001-8849-7688<EMAIL_ADDRESS>Department of
Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. Avik De 0000-0001-6475-3085
<EMAIL_ADDRESS>Department of Mathematical and Actuarial Sciences,
Universiti Tunku Abdul Rahman, Jalan Sungai Long, 43000 Cheras, Malaysia
Sanjay Mandal 0000-0003-2570-2335<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.

###### Abstract

The fundamental nature and origin of dark energy are one of the premier
mysteries of theoretical physics. In General Relativity Theory, the
cosmological constant $\Lambda$ is the simplest explanation for dark energy.
On the other hand, the cosmological constant $\Lambda$ suffers from a delicate
issue so-called fine-tunning problem. This motivates one to modify the
spacetime geometry of Einstein’s GR. The $f(Q)$ gravity is a recently proposed
modified theory of gravity in which the non-metricity scalar $Q$ drives the
gravitational interaction. In this article, we consider a linear $f(Q)$ model,
specifically $f(Q)=\alpha Q+\beta$, where $\alpha$ and $\beta$ are free
parameters. Then we estimate the best fit values of model parameters that
would be in agreement with the recent observational data sets. We use 31
points of the CC 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.
Further, we analyze the physical behavior of cosmological parameters such as
density, deceleration, and the EoS parameters corresponding to the constrained
values of the model parameters. The evolution of deceleration parameter
predicts a transition from decelerated to accelerated phases of the universe.
Further, the evolution of equation of state parameter depicts quintessence
type behavior of the dark energy fluid part. We find that our $f(Q)$
cosmological model can effectively describe the late time cosmic acceleration
without invoking any dark energy component in the matter part.

##  I Introduction

In the last two decades, evidences from Supernova searches [1, 2], WMAP
experiment [3, 4], CMBR [5, 6], LSS [7, 8] and the BAO measurements [9, 10]
indicate an accelerating cosmological expansion. The final fate of the
universe is a topic of great concern. The root cause triggering this
cosmological expansion is attributed to some sort of negative pressure dark
energy (DE). The ultimate fate of our universe strongly depends on the
fundamental nature of dark energy. The DE is usually characterized by an
equation of state parameter $\omega_{DE}\equiv\frac{p_{DE}}{\rho_{DE}}$ that
is a ratio of spatially homogeneous pressure $p_{DE}$ to the energy density
$\rho_{DE}$ of dark energy. According to recent cosmological observations, the
ambiguities are too large to differentiate among the cases : $\omega<-1$,
$\omega=-1$, and $\omega>-1$. The value of the equation of state parameter for
dark energy obtained by the WMAP [11] which combined data from the $H_{0}$
measurements, supernovae, CMB, and BAO shows that $\omega=-1.084\pm 0.063$
while in the year 2015, the Planck collaboration indicates that
$\omega=-1.006\pm 0.0451$ [12] and further in 2018 it reported that
$\omega=-1.028\pm 0.032$ [13]. The cosmological constant $\Lambda$ in GR is
the simplest explanation for the dark energy and it is characterized by
$\omega=-1$. However, there is a high discrepancy between the observed value
of the cosmological constant $\Lambda$ and its expected value from quantum
gravity [14]. This inconsistency in the value of $\Lambda$ is referred to as a
cosmological constant problem. Another widely explored time varying DE model
is the model with quintessence dark energy that is characterized by an
equation of state $-1<\omega<-\frac{1}{3}$ [15, 16]. In such models, density
of dark energy decreases with time [17]. Further, the least theoretically
understood dark energy characterized by $\omega<-1$ is called phantom energy.
The phantom energy case has gained much attention among theorists due to its
strange properties. The phantom model represents growing dark energy that
results in an extreme future expansion which leads to finite-time future
singularity. For the classification of singularities, see the references [18,
19, 20]. It also violates all the four energy conditions that help to
constrain wormholes [21]. In this article, we follow a different mechanism in
which dark energy evolves from the gravitational sector instead of the matter
part. Such an approach has been widely used in the literature so-called
modified theories of gravity [22, 23, 24]. Recently, $f(Q)$ gravity has been
proposed by J. B. Jiménez et al [25] and it has gained much attention among
cosmologists. The $f(R)$ theory of gravity is a generalization of GR in which
the space-time is described by the non-vanishing curvature with vanishing
torsion and non-metricity [26]. The $f(\mathcal{T})$ gravity theory is a
generalization of the teleparallel equivalent of GR in which the space-time is
described by the non-zero torsion with vanishing curvature and non-metricity
[27]. Finally, the $f(Q)$ theory of gravity, which will be described in
Section II, is a generalization of the symmetric teleparallel equivalent of GR
in which the non-metricity scalar describes the gravitational interactions
with zero curvature and torsion. Furthermore, the cosmological realisation of
$f(R)$ theories remain close to those of GR, while models based on $f(T)$
theories suffer from strong coupling issues on generic FLRW backgrounds [28].
This coupling problems observed in $f(T)$ theories are absent in $f(Q)$
models. Also, at the small scale quasi-static limit the predictions of the
models based $f(T)$ and $f(Q)$ theory coincide whereas at higher scales model
of $f(Q)$ theories generically propagate two scalar degrees of freedom that
are absent in the case of $f(T)$. These two degrees of freedom disappear
around maximally symmetric backgrounds and thus, cause the discussed strong
coupling problem [29].

Although recently proposed, the $f(Q)$ gravity theory already presents some
interesting and valuable applications in the literature. The first
cosmological solutions in $f(Q)$ gravity appear in References [29, 30], while
$f(Q)$ cosmography and energy conditions can respectively be seen in [31, 32].
The geodesic deviation equation in $f(Q)$ theory from its covariant
formulation was studied [33]. Quantum cosmology have been studied for a power-
law model [34]. Cosmological solutions and growth index of matter
perturbations have been investigated for a polynomial functional form of
$f(Q)$ [35]. Isotropization process in an anisotropic Bianchi universe was
investigated in a polynomial $f(Q)$ model [36]. Harko et al. analyzed the
coupling matter in $f(Q)$ gravity by assuming a power-law function [37]. The
first evidence that non-metricity $f(Q)$ gravity can challenge $\Lambda$CDM
has been presented in [38]. A set of constraints for $f(Q)$ cosmology with a
$\Lambda$CDM background is presented in [39] using ground and space based
gravitational waves observatories. Propagation of gravitational waves and
redshift space distortions have been analyzed under the framework of $f(Q)$
gravity [40, 41]. A designer approach to $f(Q)$ gravity with signatures and
its cosmological implications presented in [42, 43]. Furthermore, there are
plenty of interesting works have been done in this theory such as black hole
solution, wormhole solution, and general covariant symmetric teleparallel
cosmology presented in [44, 45, 46]. Recently, some extensions of $f(Q)$
theory have been introduced in the literature namely $f(Q,T)$ theory [47] and
Weyl $f(Q,T)$ theory [48].

During last two decades the wealth of observational data increases. The
majority of studies has been concentrated on evidences from Type Ia
supernovae, baryon acoustic oscillations (BAO) and cosmic microwave background
(CMB). The Hubble parameter dataset shows the intricate structure of the
expansion of the universe. The ages of the most massive galaxies offer direct
measurements of the Hubble parameter $H(z)$ at different redshifts $z$,
resulting in the development of a new form of standard cosmological probe
[49]. In the present manuscript we include 31 measurements of $H(z)$ spanned
using differential age method which is treated as Cosmic chronometer (CC) data
sets [50] in further study, BAO data consisting of six points [51]. Recently,
Scolnic et al. published a large Type Ia supernovae sample named Pantheon data
sets consisting 1048 points that covers the redshift range $0.01<z<2.3$ [52].
Our analysis uses the CC, BAO and Pantheon samples to constrain the
cosmological model.

This manuscript is organized as follows. In Sec II, we present the fundamental
formulations in $F(Q)$ gravity. In Sec III, Friedmann’s equation for a flat
cosmology along with a dark energy fluid part is discussed. In Sec IV, we
consider a $F(Q)$ gravity model and derive the expressions for density,
equation of state (EoS), and the deceleration parameter. In Sec V we
constraint the model parameters by using 31 points of the CC data sets, 6
points of the BAO data sets, and 1048 points from the Pantheon supernovae
samples. Further, we investigate the physical behavior of cosmological
parameters such as density, deceleration, and the EoS parameters corresponding
to the constraint values of the model parameters. Finally, we discuss our
conclusions in Sec VI.

##  II Fundamental Formulations in $F(Q)$ Gravity Theory

The spacetime curvature demonstration of gravity is one of the most
fundamental notions in natural sciences. Yet apart from the curvature, there
are apparently two geometrical objects, namely, the torsion and the non-
metricity of the spacetime’s geodesic structure which can establish the
gravity. The standard theory of gravity governed by the General Relativity
(GR) makes use of the spacetime curvature to determine gravity, with a
torsion-free and metric-compatible connection, the very special Levi-Civita
connection. As a byproduct, it assures that the geodesics are also
autoparallel in this case. The second notion benefits from a metric-
compatible, curvature-free connection with torsion, and is called the
Teleparallel Equivalent of GR (TEGR). Whereas the last one determines gravity
in a curvature and torsion-free environment, in which the non-metricity takes
the complete charge of gravity, and is called the Symmetric Teleparallel
Equivalent of GR (STEGR). As discussed in [53], the geometrical framework of
this last setting is perhaps the simplest among the three equivalent theories
of gravity because the connection can be globally completely cut off from the
discussion by an appropriate choice of coordinates, called the coincident
gauge coordinates. Moreover, in this particular coordinates system, it entails
only the first derivatives of the metric tensor, unlike the standard theory of
gravity, and thus generates a well-posed variational principle without any
Gibbons-Hawking-York (GHY) boundary terms.

A metric-affine spacetime equipped with a metric $g_{\mu\nu}$ that encodes
angles and distances and a general affine connection $X^{\alpha}_{\mu\nu}$
that defines the notion of parallel transport and covariant derivatives. From
differential geometry, it is well-known that the generic affine connection
$X^{\alpha}_{\mu\nu}$ admits a splitting into three parts [53],

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

with the Levi-Civita connection of the metric tensor $g_{\mu\nu}$,

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

and the contortion tensor,

$K^{\alpha}_{\ \mu\nu}\equiv\frac{1}{2}(T^{\alpha}_{\ \mu\nu}+T_{\mu\ \nu}^{\
\alpha}+T_{\nu\ \mu}^{\ \alpha})$ (3)

and the distortion tensor,

$L^{\alpha}_{\ \mu\nu}\equiv\frac{1}{2}(Q^{\alpha}_{\ \mu\nu}-Q_{\mu\ \nu}^{\
\alpha}-Q_{\nu\ \mu}^{\ \alpha})$ (4)

The last two terms are called torsion tensor and non-metricity tensor
respectively and defined as

$T^{\alpha}_{\ \mu\nu}\equiv X^{\alpha}_{\ \mu\nu}-X^{\alpha}_{\ \nu\mu}$ (5)

and

$Q_{\alpha\mu\nu}\equiv\nabla_{\alpha}g_{\mu\nu}$ (6)

The geometrical framework we use has a flat and torsionless connection so that
it corresponds to pure coordinate transformation from the trivial connection
as described in [25]. The connection can be parameterized as

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

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
$X^{\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 [54, 33]. Thus in the coincident gauge
coordinate , we have

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

while in an arbitrary coordinate system,

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

It is clear from the previous discussion that in STEGR under a coincident
gauge coordinates, $X^{\alpha}_{\ \mu\nu}$ and $K^{\alpha}_{\ \mu\nu}$ vanish,
and thus from equation (1) we can conclude that

$\Gamma^{\alpha}_{\ \mu\nu}=-L^{\alpha}_{\ \mu\nu}$ (10)

using which we can calculate the required tensors and scalars in that specific
coordinates.

The gravitational interactions in modified symmetric teleparallel geometry or
$F(Q)$ gravity is governed by the following action:

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

Here $F(Q)$ is an arbitrary function of the non-metricity scalar $Q$,
$g=det(g_{\mu\nu})$ and $L_{m}$ is the Lagrangian density of matter.

Due to the symmetry of the metric tensor $g_{\mu\nu}$, We can have only two
independent traces from the non-metricity tensor $Q_{\alpha\mu\nu}$, namely,

$Q_{\alpha}=Q_{\alpha}\>^{\mu}\>{}_{\mu}\>and\>\tilde{Q}_{\alpha}=Q^{\mu}\>_{\alpha\mu}$
(12)

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

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

The non-metricity scalar is acquired by [55]

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

Moreover, the stress-energy momentum tensor for the cosmic matter content is
given by

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

We denote $F_{Q}=\frac{dF}{dQ}$ for the sake of convenience.

The field equation describing the gravitational interactions in $F(Q)$ gravity
is obtained by varying the action (11) with respect to the metric tensor as

$\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})=-T_{\mu\nu}.$
(16)

Moreover, one can obtain the following result by varying the action with
respect to the connection,

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

##  III Flat FLRW Universe in $F(Q)$ Cosmology

Taking into account the homogeneity and isotropy of the universe, we describe
our universe by the spatially flat FLRW line element [56] in Cartesian
coordinates, which is, as a matter of fact also a coincident gauge
coordinates, therefore from now connection becomes trivial and metric is the
only fundamental variable,

$ds^{2}=-dt^{2}+a^{2}(t)[dx^{2}+dy^{2}+dz^{2}]$ (18)

Here, $a(t)$ is a measure of the cosmic expansion at a cosmic time $t$, known
as the scale factor . Now one can obtain the non-metricity scalar for the line
element (18) as

$Q=6H^{2}$ (19)

The stress-energy momentum tensor characterizing the matter-content of the
universe by its matter-energy density $\rho$ and isotropic pressure $p$ for
the line element (18) is

$\mathcal{T}_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu}$ (20)

Here $u^{\mu}=(1,0,0,0)$ are components of the four velocities of the perfect
cosmic fluid.

Then the Friedmann like equations ruling the dynamics of the universe for the
function $F(Q)=-Q+f(Q)$, are [55]

$f+Q-2Qf_{Q}=2\rho$ (21)

and

$\dot{H}=\frac{\rho+p}{2\left(-1+f_{Q}+2Qf_{QQ}\right)}$ (22)

In particular, for the function $F(Q)=-Q$, we can recover the usual Friedmann
equations of GR.

Now, the trace of the field equations leads to the following matter
conservation equation

$\dot{\rho}+3H\left(\rho+p\right)=0$ (23)

The equation of state that relates the usual pressure and matter-energy
density of the cosmic fluid is given as [57]

$p=\omega\rho$ (24)

Here $\omega$ is constant called equation of state (EoS) parameter.

The equation (23) can be solved with the help of equation (24), to give

$\rho=\rho_{0}a(t)^{-3\left(1+\omega\right)}$ (25)

These equations (21) and (22) can be interpreted as symmetric teleparallel
equivalent to GR (STG) cosmology with an additional component coming due to
non-metricity of space-time that behaves like dark energy fluid part. These
dark energy components coming due to non-metricity are defined by

$\rho_{DE}=-\frac{f}{2}+Qf_{Q}$ (26)

and

$p_{DE}=-\rho_{DE}-2\dot{H}\left(f_{Q}+2Qf_{QQ}\right)$ (27)

Now, the dimensionless density parameter for the dark energy fluid part is
defined as

$\Omega_{DE}=\frac{\rho_{DE}}{3H^{2}}$ (28)

Moreover, the equation of state (EoS) parameter that relates the energy
density and pressure of the dark energy component is

$\omega_{DE}=\frac{p_{DE}}{\rho_{DE}}=-1+4\dot{H}\left(\frac{f_{Q}+2Qf_{QQ}}{f-2Qf_{Q}}\right)$
(29)

Using (21) and (22), we have

$\omega_{DE}=-1+\left(1+\omega\right)\frac{\left(f+Q-2Qf_{Q}\right)\left(f_{Q}+2Qf_{QQ}\right)}{\left(-1+f_{Q}+2Qf_{QQ}\right)\left(f-2Qf_{Q}\right)}$
(30)

Thus, the effective Friedmann equations along with a dark energy fluid part
coming due to non-metricity reads

$H^{2}=\frac{1}{3}\left[\rho+\rho_{DE}\right]$ (31)

$\dot{H}=-\frac{1}{2}\left[\rho+p+\rho_{DE}+p_{DE}\right]$ (32)

Furthermore, these dark energy components satisfy the standard continuity
equation

$\dot{\rho}_{DE}+3H\left(\rho_{DE}+p_{DE}\right)=0$ (33)

##  IV Cosmological $F(Q)$ Model

Despite aforementioned issues, cosmological constant $\Lambda$ in GR is most
successful model so far, therefore this motivate us to consider following
linear $f(Q)$ model [58],

$f(Q)=\alpha Q+\beta$ (34)

Here $\alpha$ and $\beta$ are free model parameters.

Then for this particular $f(Q)$ model, we have a first-order differential
equation for a universe consisting of non-relativistic pressureless matter
reads as

$\dot{H}\left[\alpha-1\right]+\frac{3}{2}H^{2}\left[\alpha-1+\frac{\beta}{6H^{2}}\right]=0$
(35)

By using equations (26) and (34), we obtain

$\rho_{DE}=3\alpha H^{2}-\frac{\beta}{2}$ (36)

Again, by using equations (28) and (30), we obtain

$\Omega_{DE}=\frac{\rho_{DE}}{3H^{2}}=\alpha-\frac{\beta}{6H^{2}}$ (37)

and

$\omega_{DE}=-1+\frac{\left(\alpha-1-\frac{\beta}{6H^{2}}\right)\alpha}{\left(\alpha-1\right)\left(\alpha-\frac{\beta}{6H^{2}}\right)}$
(38)

Thus, the effective EoS parameter for our model is

$w_{eff}=\frac{p_{eff}}{\rho_{eff}}=\frac{p_{DE}}{\rho+\rho_{DE}}$ (39)

where, $p_{eff}$ and $\rho_{eff}$ correspond to the total pressure and energy
density of the universe. Then we have

$w_{eff}=-\left(\alpha-\frac{\beta}{6H^{2}}\right)+\frac{\left(\alpha-1-\frac{\beta}{6H^{2}}\right)\alpha}{\left(\alpha-1\right)}$
(40)

The deceleration parameter is a key component to describe the expansion phase
of the universe, it is defined as [59]

$q=\frac{1}{2}\left(1+3\Omega_{DE}\omega_{DE}\right)$ (41)

Using (37) and (38), we get

$q=\frac{1}{2}+\frac{3}{2}\bigl{\\{}-\left(\alpha-\frac{\beta}{6H^{2}}\right)+\frac{\left(\alpha-1-\frac{\beta}{6H^{2}}\right)\alpha}{\left(\alpha-1\right)}\bigr{\\}}$
(42)

Now by solving equation (35), we obtained expression for Hubble parameter in
terms of redshift as

$H(z)=\bigl{\\{}H_{0}^{2}(1+z)^{3}+\frac{\beta}{6(\alpha-1)}\left[1-(1+z)^{3}\right]\bigr{\\}}^{\frac{1}{2}}$
(43)

where $H_{0}=67.9\pm 0.5$ km/s/Mpc [62] is the present value of the Hubble
parameter. Now, our aim is to constraint the free parameters using
observational data.

##  V Observational Constraints

To constrain the model parameters of our cosmological model, we use the most
recent Cosmic chronometer data set, BAO data set, and Supernovae observations.
We use 31 points of the CC 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 [60].

###  V.0.1 Cosmic Chronometer 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 dt. Here we use the collection of
$H(z)$ data points that consists of 31 points measured from differential age
technique (listed in [61]). To calculate mean values of the model parameters
$\alpha$ and $\beta$, we have taken the chi-square function as:

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

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

###  V.0.2 Pantheon datasets

Initially, the observational studies on sample of 50 points of type Ia
supernovae revealed that our universe is accelerating. In the last two
decades, studies on more sample of type Ia supernovae data sets has been
increased. In this manuscript, we have used a sample of 1048 spectroscopically
confirmed type Ia supernovae known as Pantheon data sets. Scolnic et al. [52]
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. For a flat universe [62],
the luminosity distance is given by

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

Here $c$ is the speed of light.

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

$\mu(z)=5log_{10}D_{L}(z)+\mu_{0},$ (46)

with

$\mu_{0}=5log(1/H_{0}Mpc)+25,$ (47)

such that

$\chi^{2}_{SN}(p_{1},....)=\sum_{i,j=1}^{1048}\bigtriangledown\mu_{i}\left(C^{-1}_{SN}\right)_{ij}\bigtriangledown\mu_{j},$
(48)

where $p_{j}$ represents the assumed model’s free parameters and $C_{SN}$ is
the covariance metric [52], and

$\displaystyle\quad\bigtriangledown\mu_{i}=\mu^{th}(z_{i},p_{1},...)-\mu_{i}^{obs}.$

Here $\mu_{th}$ represents theoretical value of the distance modulus,
$\mu_{obs}$ represents its observed.

###  V.0.3 BAO datasets

The Baryonic Acoustic Oscillation (BAO) data set consists of 6dFGS, SDSS and
WiggleZ surveys that comprise BAO measurements at six different redshifts in
table-1. The characteristic scale of BAO is governed by the sound horizon
$r_{s}$ at the photon decoupling epoch $z_{\ast}$ using 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_{0b}/4\Omega_{0\gamma})a}}$
(49)

Here $\Omega_{0b}$ and $\Omega_{0\gamma}$ represent present densities of
baryons and photons respectively.

In this work, BAO datasets of six points for $d_{A}(z_{\ast})/D_{V}(z_{BAO})$
is taken from the references [51, 63, 64, 65, 66, 67]. Here the redshift at
the photon decoupling epoch is taken as $z_{\ast}\approx 1091$ and
$d_{A}(z)=\int_{0}^{z}\frac{dz^{\prime}}{H(z^{\prime})}$ 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 distance datasets is used as [67]

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

$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}$

where

$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)\,,$

and $C^{-1}$ is the inverse covariance matrix reads as [67].

$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)\,.$

We have calculated the best fit ranges for parameters $\alpha$ and $\beta$ of
our cosmological model by minimizing the chi-square function for the
combination $CC+BAO+Pantheon$. The obtained best fit values of the model
parameters are $\alpha=0.998760\pm 0.000048$ and $\beta=-26.01\pm 0.98$. The
$1-\sigma$ and $2-\sigma$ likelihood contour for the model parameters using
the combine CC+BAO+Pantheon data sets is presented in Fig 1.

Figure 1: The $1-\sigma$ and $2-\sigma$ likelihood contours for the model
parameters using the combination CC+BAO+Pantheon data sets.

###  V.0.4 Cosmological Parameters

The evolution of different cosmological parameters of our cosmological model
such as density, deceleration, and the EoS parameters corresponding to the
constrained values of the model parameters are presented below.

Figure 2: Profile of the density of the dark energy component vs redshift z .
Figure 3: Profile of the deceleration parameter vs redshift z . Figure 4:
Profile of the EoS parameter for the dark energy component vs redshift z .
Figure 5: Profile of the effective EoS parameter vs redshift z .

Form fig. 2 it is clear that the energy density of the dark energy component
of the universe decreases with cosmic time and it falls off to be zero in the
far future. Hence our model represents a decaying dark energy model. Fig. 3
show that a transition from decelerating to accelerating phase has been
experienced by the universe in the recent past with transition redshift
$z_{t}=0.844^{+0.0027}_{-0.0025}$. Further from fig. 4 snd 5, we observe that
the EoS parameter for the dark energy fluid evolving due to non-metricity and
effective EoS parameter shows quintessence like behavior. The present values
of the deceleration parameter and EoS parameter for the dark energy fluid part
correspond to the constrained values of the model parameters are
$q_{0}=-0.63^{+0.0011}_{-0.0012}$ and $w_{0}=-0.75^{+0.0007}_{-0.0008}$.

##  VI Conclusion

Understanding the evolution of dark energy is a great challenge for modern
cosmology. Dark energy occupies nearly 68.3% of the entire universe while dark
matter and baryonic matter occupy nearly 26.8% and 4.9% of the total energy
content of the universe. However, the vacuum quantum energy can well describe
the origin of dark energy by means of the cosmological constant $\Lambda$ in
the field equations of GR. Although the aforementioned issues related to
$\Lambda$ supply the search for alternative explanations for the dark energy
candidate.

In this article, we attempted to describe the evolution of dark energy from
the geometry of spacetime. We considered a linear $f(Q)$ model $f(Q)=\alpha
Q+\beta$, where $\alpha$ and $\beta$ are free parameters. Then we found the
expressions for density, deceleration, and the EoS parameters for our
cosmological model. Further, to constrain the model parameters we used CC data
sets consisting 31 data points, 6 points of the BAO data sets, and 1048 points
from the Pantheon supernovae samples. We have calculated the best fit ranges
of the model parameters $\alpha$ and $\beta$ for the combine $CC+BAO+Pantheon$
data sets. The obtained best fit values of the model parameters are
$\alpha=0.998760\pm 0.000048$ and $\beta=-26.01\pm 0.98$. Further, we have
studied the evolution of different cosmological parameters corresponding to
these best fit values of the model parameters. The evolution trajectory of the
deceleration parameter shows that our universe had experienced a transition
from deceleration to acceleration phase in the recent past with the transition
redshift $z_{t}=0.844^{+0.0027}_{-0.0025}$. The obtained present value of the
deceleration parameter is $q_{0}=-0.63^{+0.0011}_{-0.0012}$. However, in case
of standard $\Lambda$CDM model, the present day value of the deceleration
parameter is $q_{0}\approx-0.56$ with the transition redshift $z_{t}\approx
0.66$ [62]. Therefore, the behavior of our $f(Q)$ model indicates a deviation
from that of $\Lambda$CDM. Further, the present value of EoS parameter i.e.
$w_{0}=-0.75^{+0.0007}_{-0.0008}$ for the dark energy fluid part indicating a
quintessence type behavior of our model. Thus the acquired results indicate
that our cosmological $f(Q)$ model well establishes the requirements to
describe late-time cosmic acceleration without invoking any dark energy part
in the matter content. The cosmological $f(Q)$ model considered in this
article have great significance. A power law correction to the STEGR will give
rise to branches of solution applicable either to the early universe or to
late-time cosmic acceleration [29]. Our cosmological $f(Q)$ model which is
nothing but STEGR with constant term, provides a correction to the late-time
cosmology, where they can give rise to the dark energy. Also our $f(Q)$ model
agrees with the investigation done in [35], which indicated that the universe
evolves from a radiation dominated epoch towards a matter dominated epoch and
eventually settles to an accelerated dark energy dominated epoch.

## Acknowledgments

RS acknowledges University Grants Commission(UGC), New Delhi, India for
awarding Junior Research Fellowship (UGC-Ref. No.: 191620096030). AD is
supported by the grant FRGS/1/2021/STG06/UTAR/02/1. SM acknowledges Department
of Science and Technology (DST), Govt. of India, New Delhi, for awarding
Senior Research Fellowship (File No. DST/INSPIRE Fellowship/2018/IF18D676).
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]. We are very much
grateful to the honorable referees and to the editor for the illuminating
suggestions that have significantly improved our work in terms of research
quality, and presentation.

## References

  * [1] A.G. Riess et al., Astron. J. 116, 1009 (1998).
  * [2] S. Perlmutter et al., Astrophys. J. 517, 565 (1999).
  * [3] C.L. Bennett et al., Astrophys. J. Suppl. 148, 119-134 (2003).
  * [4] D.N. Spergel et al., [WMAP Collaboration], Astrophys. J. Suppl. 148, 175 (2003).
  * [5] R.R. Caldwell, M. Doran, Phys. Rev. D 69, 103517 (2004).
  * [6] Z.Y. Huang et al., JCAP 0605, 013 (2006).
  * [7] T. Koivisto, D.F. Mota, Phys. Rev. D 73, 083502 (2006).
  * [8] S.F. Daniel, Phys. Rev. D 77, 103513 (2008).
  * [9] D.J. Eisenstein et al., Astrophys. J. 633, 560 (2005).
  * [10] W.J. Percival at el., Mon. Not. R. Astron. Soc. 401, 2148 (2010).
  * [11] G. Hinshaw et al., Astrophys. J. Suppl. Ser. 208, 19 (2013).
  * [12] P.A.R. Ade et al., Astron. Astrophys. 594, A13 (2015).
  * [13] N. Aghanim et al., Astron. Astrophys. 641, A6 (2020).
  * [14] S.Weinberg, Rev. Mod. Phys. 61, 1 (1989).
  * [15] B. Ratra and P.J.E. Peebles, Phys. Rev. D 37, 3406 (1998).
  * [16] M. S. Turner and M. White, Phys. Rev. D 56, 4439-4443 (1997).
  * [17] L. Xu et al., Phys. Rev. D 84, 123004 (2011).
  * [18] R.R. Caldwell and M. Kamionkowski, and N.N. Weinberg, Phys. Rev. Lett. 91, 071301 (2003).
  * [19] J.D. Barrow, Class. Quant. Grav. 21, L79 (2004).
  * [20] J.D. Barrow, Class. Quant. Grav. 21, 5619 (2004).
  * [21] Z. Hassan, S. Mandal and P.K. Sahoo, Fortschr. Phys. 69, 2100023 (2021).
  * [22] L. Amendola et al., Phys. Rev. D 75, 083504 (2007).
  * [23] S. Appleby and R. Battye, Phys. Lett. B 654, 7 (2007).
  * [24] R. Ferraro, F. Fiorini, Phys. Rev. D 75, 084031 (2007).
  * [25] J. B. Jiménez et al., Phys. Rev. D 98, 044048 (2018).
  * [26] A. A. Starobinsky, Pisma Zh. Eksp. Teor. Fiz. 30, 719(1979).
  * [27] G.Bengochea, R. Ferraro, Phys. Rev. D 79, 124019(2009).
  * [28] A. Golovnev and T. Koivisto, JCAP 11, 012 (2018).
  * [29] J.B. Jiménez et al.,  Phys. Rev. D 101, 103507(2020).
  * [30] W. Khyllep et al.,  Phys. Rev. D 103, 103521(2021).
  * [31] S. Mandal et al.,  Phys. Rev. D 102, 124029(2020).
  * [32] S. Mandal et al.,  Phys. Rev. D 102, 024057(2020).
  * [33] J.T. Beh, T.H. Loo, and A. De, Chinese Journal of Physics, DOI: https://doi.org/10.1016/j.cjph.2021.11.026 (2021).
  * [34] N. Dimakis, A. Paliathanasis, and T. Christodoulakis, arXiv, arXiv:2108.01970 (2021).
  * [35] W. Khyllep, A. Paliathanasis, and J. Dutta, Phys. Rev. D 103, 103521 (2021).
  * [36] A. De, S. Mandal, J.T. Beh, T.H. Loo and P.K. Sahoo, Eur. Phy. J. C. 82, 72 (2022).
  * [37] T. Harko et al., Phys. Rev. D 98, 084043 (2018).
  * [38] F. K. Anagnostopoulos, S. Basilakos, and E. N.Saridakis, Phys. Lett. B 822, 136634 (2021).
  * [39] J. Ferreira et al., arXiv, arXiv:2203.13788 (2022).
  * [40] M. Hohmann et al, Phys. Rev. D 99, 024009 (2019).
  * [41] B. J. Barros et al., Phys.Dark Univ. 30, 100616 (2020).
  * [42] I. S. Albuquerque and N. Frusciante, arXiv, arXiv:2202.04637 (2022).
  * [43] N. Frusciante, Phys. Rev. D 103, 044021 (2021).
  * [44] F. D Ambrosio et al., Phys. Rev. D 105, 024042 (2022).
  * [45] G. Mustafa et al., Phys. Lett. B 821, 136612 (2021).
  * [46] M. Hohmann, Phys. Rev. D 104, 124077 (2021).
  * [47] Y. Xu et al., Eur. Phys. J. C 79, 708 (2019).
  * [48] Y. Xu et al., Eur. Phys. J. C 80, 449 (2020).
  * [49] R. Jimenez, A. Loeb, Astrophys. J. 573, 37(2002).
  * [50] G. S. Sharov et al., Mon. Not. R. Astron. Soc. 466, 3497(2017).
  * [51] C. Blake et al.,  Mon. Not. Roy. Astron. Soc. 418, 1707 (2011).
  * [52] D.M. Scolnic et al., Astrophys. J. 859, 101(2018).
  * [53] T. Ortin, Gravity and Strings, Cambridge Monographs on Mathematical Physics (Cambridge University Press (2015).
  * [54] N. Dimakis, A. Paliathanasis, and T. Christodoulakis,  arXiv, arXiv:2108.01970 (2021).
  * [55] R. Lazkoz et al., Phys. Rev. D 100, 104027 (2019).
  * [56] B. Ryden, Introduction to Cosmology (Addison Wesley, SanFrancisco, United States of America, 2003).
  * [57] S. Chakrabarti, S. L. Said, and K. Bamba,  Eur. Phys. J. C. 79, 454 (2019).
  * [58] G. Mustafa et al., Annals Phys. 437, 168751 (2022).
  * [59] M. Jamil et al.,  J. Phys. Soc. Japan. 81, 114004 (2012).
  * [60] D. F. Mackey et al., Publ. Astron. Soc. Pac. 125, 306(2013).
  * [61] A. Gómez-Valent1 and Luca Amendola, JCAP 04, 051 (2018).
  * [62] Planck Collaboration,  Astron. Astrophys. 641, A6(2020).
  * [63] W. J. Percival et al., Mon. Not. Roy. Astron. Soc. 401, 2148 (2010).
  * [64] F. Beutler et al.,  Mon. Not. Roy. Astron. Soc. 416, 3017 (2011).
  * [65] N. Jarosik et al.,  Astrophys. J. Suppl. 192, 14 (2011).
  * [66] D. J. Eisenstein et al.,  Astrophys. J. 633, 560 (2005).
  * [67] R. Giostri et al., J. Cosm. Astropart. Phys. 1203, 027 (2012).