# Towards a unified interpretation of the early Universe in $R^{2}$-corrected
dark energy model of $F(R)$ gravity

Hua Chen<EMAIL_ADDRESS>Institute of Astrophysics, Central China
Normal University, Wuhan 430079, China  Taishi Katsuragawa
<EMAIL_ADDRESS>Institute of Astrophysics, Central China Normal
University, Wuhan 430079, China  Shinya Matsuzaki<EMAIL_ADDRESS>Center for
Theoretical Physics and College of Physics, Jilin University, Changchun,
130012, China

###### Abstract

$R^{2}$-corrected dark energy (DE) models in $F(R)$ gravity have been widely
investigated in recent years, which not only removes the weak singularity
potentially present in DE models but also provide us with a unified picture of
the cosmic history, including the inflationary and DE epochs. Towards the
unified interpretation of dynamical DE all over the cosmic history in the
class of $R^{2}$-corrected DE models, we explore the universal features of the
scalaron dynamics in the radiation-dominated epoch, along with the chameleon
mechanism, by keeping our eyes on the inflationary and DE epochs. We show that
the scalaron evolution does not follow a surfing solution and is mostly
adiabatic before big bang nucleosynthesis (BBN), even properly including the
kick by the nonperturbative QCD phase transition, hence a catastrophic
consequence claimed in the literature is not applied to this class of DE
models. This is due to the presence of the gigantic scale hierarchy between
$R^{2}$ correction and DE, so is the universal feature for the class of
$R^{2}$-corrected DE models. The prospects for the post- or onset-inflationary
epoch would be pretty different from what the standard $R^{2}$ inflationary
scenario undergoes due to the presence of the chameleon mechanism.

##  I Introduction

Growing observational evidence indicates two accelerations of cosmic
expansion: the early-time inflation expected to occur before radiation
domination Akrami  _et al._ (2020); Martin (2016), and the late-time
acceleration sourced by dark energy (DE) Riess  _et al._ (1998); Perlmutter
_et al._ (1999) (See Zyla  _et al._ (2020) for reviews). There are two ways to
realize an accelerated phase: to modify the matter content by introducing a
new component, with the equation of state $w<-1/3$, into the matter
Lagrangian, or to modify the gravitational theory Nojiri  _et al._ (2017);
Heisenberg (2019). For instance, $F(R)$ gravity is one of the simplest
extensions of general relativity (GR) by improving the Ricci scalar $R$ to be
an arbitrary function $F(R)$ De Felice and Tsujikawa (2010). As a result,
$F(R)$ gravity introduces one dynamical-field degree of freedom, dubbed
scalaron, which has the gravitational origin and could source the inflation
and DE without requiring exotic components into the matter Lagrangian.

Although phenomenologically successful, it is known that viable DE models in
$F(R)$ gravity generally suffer from a problem of the weak curvature
singularity with $\left|R\right|\to\infty$ at finite field value Appleby and
Battye (2008); Frolov (2008). A promising way to cure this problem is to add
higher-curvature corrections, for instance, the $\alpha R^{2}$ term Dev  _et
al._ (2008); Kobayashi and Maeda (2009); Appleby  _et al._ (2010); Lee  _et
al._ (2012). Moreover, $R^{2}$ gravity is also well known as the inflation
model Starobinsky (1980), which has proven very successful in light of the
cosmic microwave background observation Akrami  _et al._ (2020). The two
successes above may imply $R^{2}$ gravity as a viable unified theory. In
addition to removing the curvature singularity, $R^{2}$-corrected DE models
could naturally realize the two accelerations by the same scalaron field
Odintsov and Oikonomou (2016); Odintsov  _et al._ (2017); Nojiri  _et al._
(2020); Oikonomou (2021a); Odintsov and Oikonomou (2019, 2020, 2021);
Oikonomou and Giannakoudi (2022a). Other setups have also been considered, for
instance, $R+\alpha R^{n}-\beta R^{2-n}$ Artymowski and Lalak (2014) and its
constraints Geng  _et al._ (2015); Yashiki  _et al._ (2020), and a single
scalar field extension to $F(R)$ gravity Oikonomou (2021b); Oikonomou and
Giannakoudi (2022b).

One may then wonder if such a unified model of $F(R)$ gravity can reproduce
the overall cosmic history as much as we currently know. Although the fact of
a gigantic hierarchy between the inflation and DE scales has allowed us to
ignore one of two epochs in studying one, $R^{2}$-corrected DE models in
$F(R)$ gravity should be able to explain the intermediate epochs in the early
Universe, such as the radiation-dominated epoch. The modification for DE might
have nontrivial effects on the potential structure of the scalaron, even in
the early Universe. In that case, the complete picture of the whole cosmic
history should be examined along with the nontrivial contributions to a link
bridging over the early and later epochs. The research on such a link will
shed light on the particle cosmology inspired by modified gravity in the early
Universe.

This work makes the first step along this new avenue: we investigate the
dynamic scalaron evolution in the radiation-dominated epoch, including
dynamics of matter fields, in correlation with the inflation and DE epochs.
The exponential coupling with matter fields makes the scalaron possess an
environment-dependent mass, so that it can pass gravitational constraints in
the local and dense solar systems and act as DE at the sparse cosmic scale De
Felice and Tsujikawa (2010). Due to such a chameleon mechanism, the decoupling
of SM particles from the thermal equilibrium in the early Universe has a
significant impact on the dynamics of the chameleon field in terms of scalar-
tensor theories or the scalaron in $F(R)$ gravity. In particular, it was shown
that there exists a novel solution, called the surfing solution, so that
chameleons would be kicked by matter fields through the chameleon mechanism,
and be surfed in a constant and over the $\mathrm{GeV^{2}}$ velocity,
resulting in nonadiabatic changes in the chameleon’s effective mass and the
explosive production of scalarons with trans-Planckian momenta Erickcek  _et
al._ (2013, 2014).

Besides the SM-particle ensemble kicks (SMPE kicks) due to the particle
decoupling, we show that there is another powerful kick due to the QCD phase
transition, which we call the QCDPT kick in this work. The QCDPT kick would
exacerbate the catastrophic consequence reported in Erickcek  _et al._ (2013,
2014). Using data on the temperature dependence of the equation of state
observed as the crossover reported from the nonperturbative lattice QCD
simulation Bali  _et al._ (2014), we convert it into the change in the energy-
momentum tensor around the QCDPT epoch and show that the QCDPT kick is more
significant on the scalaron evolution than those estimated from the
conventional thermal decoupling of the SM particles based on the perfect fluid
approximation. Nevertheless, in terms of the scalaron in $F(R)$ gravity, we
find that the surfing solution is always absent, owing to the unique and weak
matter coupling. Thus, the kicks will be harmless and merely stabilize the
scalaron at its effective potential minimum.

We further find that the scalaron evolution at the radiation-dominated epoch
is mostly adiabatic before BBN, so the catastrophic consequence claimed in the
literature is not applied to this class at all. The null of catastrophe is
universal in this class of $R^{2}$-corrected DE models because it essentially
stems from the hierarchy between the $R^{2}$ correction and DE scales. Thus,
the $R^{2}$ correction is also a solution to the catastrophic consequence, as
well as the curvature singularity in general DE models of $F(R)$ gravity. The
prospects for the post- or onset-inflationary epoch would differ from the
standard $R^{2}$ inflationary scenario, which will also be briefly addressed.

This paper is organized as follows. In Section II, we briefly review the
$R^{2}$-corrected DE model of $F(R)$ gravity and formulate the scalaron field
equation in the flat Friedmann–Lemaître–Robertson–Walker (FLRW) spacetime. In
Section III, we investigate the chameleon mechanism in the radiation-dominant
epoch and demonstrate the kick solutions by the SMPE and QCDPT kicks.
Moreover, we evaluate the time evolution of the scalaron mass to show that the
$R^{2}$ term suppresses the nonadiabaticity. In Section IV, we briefly discuss
the prospects of the post- or onset-inflationary epoch characteristic of the
$R^{2}$-corrected DE models.

Throughout this paper, a “dot” or “prime” denotes a derivative with respect to
cosmic time or the argument of a function, respectively. Especially
$\varphi^{\prime}=d\varphi/d\tilde{N}$, where $\tilde{N}$ is the number of
e-folds with “tilde” denoting a physical quantity in the Einstein frame.

##  II $F(R)$ Gravity Theory

###  II.1 Rudiments of $F(R)$ Cosmology

We begin by introducing the action of $F(R)$ gravity, which is given by

$\displaystyle S=\int d^{4}x\ \sqrt{-g}\left[\frac{1}{2\kappa^{2}}F(R)+{\cal
L}_{\text{M}}\left(g^{\mu\nu},\Phi_{\text{M}}\right)\right]\,,$ (1)

where $\kappa^{2}\equiv 8\pi G\equiv 1/M_{\text{Pl}}^{2}$, $g$ is the
determinant of Jordan-frame metric $g_{\mu\nu}$, and $\mathcal{L}_{\text{M}}$
is the matter Lagrangian as a function of $g_{\mu\nu}$ and matter fields
$\Phi_{\text{M}}$. Varying the action with respect to $g_{\mu\nu}$ gives the
field equation

$\displaystyle
F^{\prime}(R)R_{\mu\nu}-\frac{1}{2}F(R)g_{\mu\nu}-\nabla_{\mu}\nabla_{\nu}F^{\prime}(R)+g_{\mu\nu}\square
F^{\prime}(R)=\kappa^{2}T_{\mu\nu}^{(\text{M})}\,,$ (2)

where $\square\equiv\nabla^{\mu}\nabla_{\mu}$. One finds that a dynamical
field degree of freedom shows up when taking the trace

$\displaystyle 3\square F^{\prime}(R)+F^{\prime}(R)R-2F(R)=\kappa^{2}T_{\
\mu}^{\mu(\text{M})}\,.$ (3)

Under the Weyl transformation,

$\displaystyle\tilde{g}_{\mu\nu}$
$\displaystyle=e^{2\sqrt{1/6}\kappa\varphi}g_{\mu\nu}\equiv
F^{\prime}(R)g_{\mu\nu}\,,$ (4)

this field $\varphi$, dubbed scalaron, appears explicitly in the Einstein-
frame action,

$\displaystyle\tilde{S}=\int d^{4}x\
\sqrt{-\tilde{g}}\left[\frac{1}{2\kappa^{2}}\tilde{R}-\frac{1}{2}\tilde{g}^{\mu\nu}(\partial_{\mu}\varphi)(\partial_{\nu}\varphi)-V(\varphi)+e^{-4\sqrt{1/6}\kappa\varphi}\mathcal{L}_{\text{M}}\right]\,,$
(5)

where the potential $V(\varphi)$ is defined as

$\displaystyle V(\varphi)$
$\displaystyle\equiv\frac{1}{2\kappa^{2}}\frac{RF^{\prime}(R)-F(R)}{F^{\prime
2}(R)}\,,$ (6)

with $R=R(\varphi)$ given by the Weyl transformation in Eq. (4). Varying Eq.
(5) with respect to the Einstein-frame metric $\tilde{g}_{\mu\nu}$ yields

$\displaystyle\tilde{R}_{\mu\nu}-\frac{1}{2}\tilde{R}\tilde{g}_{\mu\nu}-\kappa^{2}\left[-\frac{1}{2}\tilde{g}_{\mu\nu}\tilde{g}^{\rho\sigma}(\partial_{\rho}\varphi)(\partial_{\sigma}\varphi)+(\partial_{\mu}\varphi)(\partial_{\nu}\varphi)-\tilde{g}_{\mu\nu}V(\varphi)\right]$
$\displaystyle=\kappa^{2}\tilde{T}_{\mu\nu}^{(\text{M})}\,.$ (7)

In flat FLRW spacetime with perfect fluid approximation, $T_{\
\mu}^{\mu(\text{M})}=-\rho_{\text{M}}+3P_{\text{M}}$, Eq. (7) gives the
Friedmann equation,

$\displaystyle
3\tilde{H}^{2}=\kappa^{2}\left[\frac{1}{2}\left(\frac{d\varphi}{d\tilde{t}}\right)^{2}+V(\varphi)+\tilde{\rho}_{\text{M}}\right]\,,$
(8)

and the acceleration equation,

$\displaystyle\frac{\ddot{\tilde{a}}}{\tilde{a}}=$
$\displaystyle-\frac{\kappa^{2}}{6}\left[2\left(\frac{d\varphi}{d\tilde{t}}\right)^{2}-2V(\varphi)+\tilde{\rho}_{\text{M}}+3\tilde{P}_{\text{M}}\right]\,.$
(9)

Here, $a$ is the scale factor and $H\equiv\dot{a}/a$ is the Hubble parameter.
Varing Eq. (5) with respect to $\varphi$, the equation of motion of the
scalaron reads

$\displaystyle\frac{d^{2}\varphi}{d\tilde{t}^{2}}+3\tilde{H}\frac{d\varphi}{d\tilde{t}}=-\frac{dV}{d\varphi}-\frac{\kappa}{\sqrt{6}}e^{-4\sqrt{1/6}\kappa\varphi}T_{\
\mu}^{\mu(\text{M})}\,,$ (10)

from which we define the effective potential ignoring the $\varphi$-dependence
in the matter Lagrangian

$\displaystyle V_{\text{eff}}(\varphi)$
$\displaystyle=V(\varphi)-\frac{1}{4}e^{-4\sqrt{1/6}\kappa\varphi}T_{\
\mu}^{\mu(\text{M})}\,.$ (11)

Then the mass of the scalaron is defined as

$\displaystyle m_{\varphi}^{2}\equiv
V_{\text{eff}}^{\prime\prime}(\varphi)=V^{\prime\prime}(\varphi)-\frac{2\kappa^{2}}{3}e^{-4\sqrt{1/6}\kappa\varphi}T_{\
\mu}^{\mu(\text{M})}\,.$ (12)

For later convenience, we define some dimensionless quantities: the density
ratio of cold matter $\tilde{\rho}_{\text{m}}$ and radiation
$\tilde{\rho}_{\text{r}}$, i.e.
$f_{\text{m}}\equiv\tilde{\rho}_{\text{m}}/\tilde{\rho}_{\text{r}}$, the kick
function
$\Sigma\equiv(\tilde{\rho}_{\text{r}}-3\tilde{P}_{\text{r}})/\tilde{\rho}_{\text{r}}$,
the scalaron $\phi\equiv\kappa\varphi$, and the number of e-folds
$\tilde{N}\equiv\ln(\tilde{a}/\tilde{a}_{i})$. Then, Eqs. (8) and (10) can be
recast as

$\displaystyle\tilde{H}$
$\displaystyle=\kappa\sqrt{\frac{V+\left(1+f_{\text{m}}\right)\tilde{\rho}_{\text{r}}}{3\left(1-\frac{1}{6}\phi^{\prime
2}\right)}}\,,$ (13)

and

$\displaystyle\frac{V+\left(1+f_{\text{m}}\right)\tilde{\rho}_{\text{r}}}{3\left(1-\frac{1}{6}\phi^{\prime
2}\right)}\phi^{\prime\prime}+\left(V+\frac{2+3f_{\text{m}}+\Sigma}{6}\tilde{\rho}_{\text{r}}\right)\phi^{\prime}$
$\displaystyle=-\frac{dV}{d\phi}+\frac{1}{\sqrt{6}}\left(\Sigma+f_{\text{m}}\right)\tilde{\rho}_{\text{r}}\,,$
(14)

where $\phi^{\prime}=d\phi/d\tilde{N}$. Since in the radiation dominated
epoch, $f_{\text{m}}\lesssim 10^{-6}$ Erickcek  _et al._ (2014), we shall
hereafter ignore $f_{\text{m}}$. Finally, the scalaron velocity with respect
to the cosmic time can be easily read off as

$\displaystyle\frac{d\varphi}{d\tilde{t}}=\phi^{\prime}\sqrt{\frac{V+\tilde{\rho}_{\text{r}}}{1-\frac{1}{6}\phi^{\prime
2}}}\,.$ (15)

###  II.2 $R^{2}$-corrected DE models

We now consider $R^{2}$-corrected DE models with

$\displaystyle F(R)=F_{\text{DE}}(R)+\alpha R^{2}\,,$ (16)

and take the Starobinsky’s DE model as a benchmark:

$\displaystyle F(R)=R-\beta
R_{\text{c}}\left[1-\left(1+\frac{R^{2}}{R_{\text{c}}^{2}}\right)^{-n}\right]+\alpha
R^{2}\,,$ (17)

where $\beta>0$, $n>0$, and $R_{\text{c}}\sim\Lambda\simeq 4\times 10^{-84}\
[\mathrm{GeV}^{2}]$ are model parameters. The potential is a multi-valued
function (see, e.g., Refs. Frolov (2008); Katsuragawa  _et al._ (2019)). In
the present study we focus on the case where $+\infty>R>R_{\text{b}}$ with the
branch located at $R_{\text{b}}\simeq 0.93R_{\text{c}}$ for $\beta=2$ and
$n=1$. The fifth force experiment constrains $\alpha<10^{22}\
[\mathrm{GeV}^{-2}]$ Cembranos (2009), and $\alpha=1/(6m_{0}^{2})\sim
10^{-27}\ [\mathrm{GeV}^{-2}]$ if one improves $\alpha R^{2}$ to be
responsible for the inflation with mass scale $m_{0}\sim 10^{13}\
\mathrm{[GeV]}$.

As is in the case of the $R^{2}$ inflation, the scalaron potential has a flat
plateau, which realizes the slow-roll inflation, as shown in Fig. 1.

Figure 1:  Black solid curve shows the potential of the scalaron in the class
of $R^{2}$-corrected model for $F(R)$ gravity, where we have chosen
$\alpha=1/6\times 10^{-26}\ \mathrm{[GeV^{-2}]}$, $\beta=2$ and $n=1$. The
scalaron field has been normalized by the Planck scale in the plot
($\phi=\kappa\varphi$). The potential minimum is set so as to realize the DE
scale $V_{0}\equiv\frac{R_{\text{c}}}{2\kappa^{2}}\sim\rho_{\Lambda}\approx
10^{-47}\ [\mathrm{GeV}^{4}]$. Though being visibly shallow at the potential
minimum, the potential asymptotically reaches infinity, $V(\varphi)\to\infty$,
when $\varphi\to-\infty$. The red dashed curve shows $R^{2}$ inflation for
comparison.

On the contrary to $R^{2}$ inflation, the potential possesses a shallow plain
with a global minimum at $\varphi\simeq 0.11M_{\text{Pl}}$, responsible for
the DE vacuum in the current Universe. Therefore, one can ignore the DE
modification only when discussing the inflationary epoch triggered by the
$R^{2}$ term; as will be clarified later, however, the DE modification can
alter the standard cosmic history after the inflation, even in the early
Universe.

We note that the structure of the potential above is common among the unified
models of the $R^{2}$ inflation and DE. There should be two (almost) flat
parts of potential to realize the two accelerations in the early and late-time
Universe, because cosmological observations indicate that inflation is
described by quasi-de Sitter spacetime in the early Universe and that DE
almost acts as the cosmological constant in the late-time Universe. Since the
viable DE model should reproduce the $\Lambda$-CDM model in the large-
curvature limit, the model-dependence is negligible in the regime where
$\phi>0$. Moreover, considering the matter contribution to the effective
potential, the flat potential for the DE background is overwhelmed by the
trace of the energy-momentum tensor in the region where $\phi<0$. In other
words, we can discuss the intermediate stage between the above two
accelerations in a model-independent way, although our analysis in the present
work specifies the DE model to the Starobinsky model.

Furthermore, the potential structure connecting two flat potentials by a
gigantic scale difference is related to the hierarchy problem of DE. Thus,
even though the $R^{2}$ correction will not be identified as $R^{2}$
inflation, satisfying the constraint on the parameter $\alpha$, the flat
potential for the DE background is always smaller than that for the $R^{2}$
correction. The details of potential analysis and argument about the model
independence are summarized in Appendix. A.

##  III Scalaron in Radiation domination

###  III.1 Absence of Surfing Solution

After the production of SM particles via the reheating epoch, the Universe is
almost full of radiation, with the trace of the energy-momentum tensor
vanishing. Due to the cooling of the Universe, however, SM particles decouple
successively from the thermal plasma, resulting in nonzero contributions to
the trace of the energy-momentum tensor, and thereby affect the scalaron
dynamics due to the matter coupling. Therefore, first of all, we estimate the
energy density at the radiation-dominated epoch, $\rho_{\text{r}}(T)$, and
$\Sigma(T)$.

Assuming the conventional ensemble of SM particles in the perfect fluid to
well describe the early Universe, the trace of the energy-momentum tensor for
each relativistic particle is evaluated

$\displaystyle\rho_{i}-3P_{i}$
$\displaystyle=\frac{g_{i}T^{4}}{2\pi^{2}}x^{2}\int_{0}^{\infty}dy\
\frac{y^{2}}{\sqrt{x^{2}+y^{2}}}\frac{1}{\exp\left(\sqrt{x^{2}+y^{2}}\right)\pm
1}\,,$ (18)

where $g_{i}$ is the degree of freedom, $x\equiv m_{i/T}$, $y\equiv P_{i}/T$,
with “$i$” labelling different SM species, and “$+/-$” are applied to fermion
and boson, respectively. The energy density of thermal bath $\rho_{\text{r}}$
is obtained by

$\displaystyle\rho_{\text{r}}\equiv\frac{\pi^{2}}{30}g_{*}(T)T^{4}\,,$ (19)

where $g_{*}(T)$ is relativistic effective degrees of freedom.

Then the kick function $\Sigma(T)$ reads

$\displaystyle\Sigma_{i}(T)=\frac{\rho_{i}-3P_{i}}{\rho_{\text{r}}}=\frac{15}{\pi^{4}}\frac{g_{i}}{g_{*}(T)}x^{2}\int_{0}^{\infty}dy\frac{y^{2}}{\sqrt{x^{2}+y^{2}}}\frac{1}{\exp\left(\sqrt{x^{2}+y^{2}}\right)\pm
1}\,.$ (20)

Summing up all contributions from the thermal decoupling of SM particles, we
find four peaks corresponding to what we call the SMPE kicks, which would
affect the scalaron evolution significantly, as will be seen later.

Using the conservation of the entropy density per comoving volume in the
Jordan frame, $sa^{3}=\frac{2\pi^{2}}{45}g_{*s}(T)T^{3}a^{3}=\text{const.}$,
one can extract $\Sigma(\tilde{N})$ by making a replacement

$\displaystyle
Tg_{*s}^{1/3}(T)=g_{*s}^{1/3}(T_{i})T_{i}e^{\sqrt{1/6}\left(\phi-\phi_{i}\right)-\tilde{N}}$
(21)

to obtain $T(\tilde{N})$.

It has been shown that there exists a novel solution, dubbed surfing solution,
in scalar-tensor theories, which would cause nonadiabaticity and explosive
particle production Erickcek  _et al._ (2013, 2014). Here, however, we show
that in the case of $F(R)$ gravity, there exists no surfing solution.

In scalar-tensor theories, the conformal factor, related by the Weyl
transformation, as in Eq.(4), includes free parameters dependent of particle
spices, $Q_{i}$, like

$\displaystyle\tilde{g}_{\mu\nu}=e^{2Q_{i}\kappa\varphi}g_{\mu\nu}\,.$ (22)

In the case of $F(R)$ gravity, we have the universal coupling
$Q=Q_{i}=1/\sqrt{6}$. Since the DE scale is very small, one can assume that
$V(\phi)\ll\rho_{\text{r}}$ and $dV/d\phi$ is negligible compared to the
driving term in $\Sigma$ in the early Universe. Then Eq. (14) is reduced to

$\displaystyle\phi^{\prime\prime}=\left(1-\frac{1}{6}\phi^{\prime
2}\right)\left[-\left(1+\frac{\Sigma}{2}\right)\phi^{\prime}+3Q\Sigma\right]\,.$
(23)

Eq. (21) allows for the existence of a surfing solution, i.e.
$\phi^{\prime\prime}=0$ and $\phi^{\prime}=1/Q$. This solution is valid if and
only if the kick function satisfies

$\displaystyle\Sigma(T_{\text{s}})\geq\frac{2}{6Q^{2}-1}\,.$ (24)

One finds that this condition does not hold in $F(R)$ gravity with
$Q=1/\sqrt{6}$. Thus, the surfing solution is absent in $F(R)$ gravity no
matter how the kick function varies. Note that our condition is different from
that in Erickcek  _et al._ (2013, 2014) since we do not assume that
$1+\Sigma\simeq 1$ in this work. We also emphasize that the absence of the
surfing solution is independent of specific $F(R)$ models because the
consequence relies only on the coupling constant $Q_{i}$ uniquely determined
by the Weyl transformation, although it is an arbitrary parameter in the
general scalar-tensor theory.

###  III.2 SMPE and QCDPT kicks

It turns out that even without a surfing solution, the SMPE and QCDPT kicks
have a nontrivial impact on the dynamics of scalaron. Utilizing data about
effective degrees of freedom from Ref. Husdal (2016) together with the result
on the equation of state from the lattice QCD simulation Bali  _et al._
(2014), ranging from $T=114-305$ [MeV], we plot the effective degrees of
freedom for the energy density and entropy density in Fig. 2, compare the
trace of the energy-momentum tensor with the energy density in Fig. 3, and
thereby obtain the kick function in Fig. 4. We find that the strongest kick
among the SMPE kicks comes from the electron-positron annihilation, which
yields $\Sigma\simeq 0.1$, while the QCDPT kick is by six times greater:
$\Sigma\simeq 0.6$. The detail of data handling is shown in Appendix. B.

Figure 2:  The effective degrees of freedom as a function of the Jordan-frame
temperature for entropy (solid blue line) and energy density (orange dashed
line).  Figure 3:  Comparison of the trace of the energy-momentum tensor and
the energy density of radiation as a function of the Jordan-frame temperature.
Figure 4:  The kicks as a function of the Jordan-frame temperature during the
radiation-dominant epoch, including the QCD phase transition and the perfect-
fluid description ensemble SM particles in the thermal equilibrium. The solid
curve around the center of the plot denotes the contribution from the QCDPT
kick. The dashed curves have been drawn by accumulating the ensemble
contributions from the thermal decoupling of SM particles.

We discuss the implications of the QCDPT kick to the scalaron evolution, in
comparison with the four conventional SMPE kicks. First of all, we note that
due to the big hierarchy of the inflation and DE scales, the scalaron cannot
climb up the potential barrier at $\varphi\geq 0$, so that it rebounds after
hitting the potential barrier, as seen from Fig. 5. Meanwhile, we observe that
the four SMPE kicks and full five kicks push the scalaron at most $0.64\
M_{\text{Pl}}$ and $1.41\ M_{\text{Pl}}$, respectively, which means that the
QCDPT kick is more effective than the SMPE kicks. In terms of scalar-tensor
theories, the surfing solution takes effect if
$Q\geq\sqrt{\frac{2+\Sigma(T_{\text{s}})}{6\Sigma(T_{\text{s}})}}\gtrsim
0.85$, and the velocity can be as large as $\phi^{\prime}=1/Q\gtrsim 1.17$,
which exacerbates the existing catastrophic results claimed in Erickcek  _et
al._ (2013, 2014).

Figure 5:  The scalaron evolution under the SMPE kicks (blue solid line) and
total kicks (orange dashed line), where just as sample points, we have chosen
initial conditions $\varphi_{i}/M_{\text{Pl}}=-(1.2,0.6,0.11,0.01)$ at
$T_{i}=10000$ [GeV], with $g_{*s}\left(T_{i}\right)=106.75$.

On the other hand, as seen from Fig. 5, the matter contribution related to the
chameleon mechanism in $F(R)$ gravity gets the opposite effect compared with
that in scalar-tensor theories. The chameleon mechanism indeed protects the
scalaron from leaving the potential minimum. In other words, the balance
between the kicks, bare potential, and Hubble friction will eventually
stabilize the scalaron at the potential minimum. This brings two advantages to
$F(R)$ gravity compared with other models: 1\. surfing solution is absent in
$F(R)$ gravity; 2. No specific initial conditions are required.

To examine the second point, we determine the scalaron field value in the BBN
epoch. First, note that the Weyl transformation, in Eq. (4), implies that a
constant mass scale $m_{i}$ in the Jordan frame will be modified due to the
$\varphi$ dependence like
$\tilde{m}_{i}=m_{i}\exp\left(-\sqrt{1/6}\kappa\varphi\right)$ in the Einstein
frame. Thus the proton-neutron mass difference could also be altered by the
presence of the scalaron field with somewhat a large field value. As was
addressed in Chen  _et al._ (2020), the size of the scalaron filed, normalized
to the Planck scale, is constrained so as not to spoil the successful BBN.
When a conservative limit of 10% in size (in the unit of the Planck scale) is
taken, we find that the allowed interval of the scalaron field value during
the BBN epoch is $-0.33\lesssim\kappa\varphi_{\text{BBN}}\lesssim 0.15$, which
requires the initial value as $\varphi_{i}\gtrsim-0.97\ M_{\text{Pl}}$ or
$\varphi_{i}\gtrsim-1.74\ M_{\text{Pl}}$ under the SMPE kicks or the total
five kicks. More stringent bound as evaluated in Chen  _et al._ (2020) might
be placed, where $-0.05\lesssim\kappa\varphi_{\text{BBN}}\lesssim 0.03$. In
that case, we would have $\varphi_{i}\gtrsim-0.69\ M_{\text{Pl}}$ with only
the four SMPE kicks taken into account, or $\varphi_{i}\gtrsim-1.46\
M_{\text{Pl}}$ with the total kicks considered111 We have assumed the field
value to be negative at this point since we just regard the unified model as a
dark energy model and have not related it with inflation yet.. Thus, the QCDPT
kick is also significant in relaxing the initial value of the scalaron field.

##  IV Nonadiabaticity in Dynamics of Scalaron

The DE scale is extremely small compared to other characteristic scales in the
early Universe, in which scale distance the scalaron travels. Hence the motion
of the scalaron could get too fast to be adiabatic. The adiabaticity is
violated when

$\displaystyle\left|\frac{d}{d\tilde{t}}\frac{1}{m_{\varphi}\left(\tilde{t}\right)}\right|\gtrsim
1\,.$ (25)

From Eqs. (12) and (15), we have

$\displaystyle\left|\frac{d}{d\tilde{t}}\frac{1}{m_{\varphi}\left(\tilde{t}\right)}\right|$
$\displaystyle=\left|-\frac{\phi^{\prime}}{2}\sqrt{\frac{V(\phi)+\tilde{\rho}_{\text{r}}}{3\left(1-\frac{1}{6}\phi^{\prime
2}\right)}}\frac{V^{\prime\prime\prime}(\phi)-\frac{8}{3\sqrt{6}}\Sigma\tilde{\rho}_{\text{r}}}{\left(V^{\prime\prime}(\phi)+\frac{2}{3}\Sigma\tilde{\rho}_{\text{r}}\right)^{3/2}}\right|\,.$
(26)

In Fig. 6 we plot this discriminator of the adiabaticity as a function of the
Jorndan-frame temperature in the radiation dominated epoch.

As shown in Fig. 6, the strong signals of the nonadiabaticity appear just
after the electron-positron annihilation at $T\sim 10^{-5}~{}[\mathrm{GeV}]$.

Figure 6:  The nonadiabaticity as a function of the Jordan-frame temperature,
where the solid-blue line and orange dashed line denote that the scalaron is
under SMPE kicks and total kicks, respectively. The initial conditions are
chosen as the same as in Fig. 5.

Although the results depend on the initial condition we have chosen, we
numerically confirmed that the large violation of adiabaticity below
$T=10^{-5}~{}[\mathrm{GeV}]$ is insensitive to the DE potential and thus the
energy density of matter. Neglecting the DE scale in Eq. (26), the
discriminator is approximated to be

$\displaystyle\left|\frac{d}{d\tilde{t}}\frac{1}{m_{\varphi}\left(\tilde{t}\right)}\right|\approx\left|\phi^{\prime}\sqrt{\frac{1}{3\Sigma\left(1-\frac{1}{6}\phi^{\prime
2}\right)}}\right|\,.$ (27)

From this form, we see that since in the conventional perfect fluid
description, the kick function $\Sigma(T)$ goes vanishing after the electron-
positron annihilation as shown in Fig. 4, it is unavoidable that the
adiabaticity is badly violated. The rapid drop of the kick function is mainly
due to the fast decrease of the energy density of radiation, compared with the
trace of the energy-momentum tensor, as shown in Fig. 3. One also notes that
the adiabatic condition is slightly violated at around $T=10^{-2}$ [GeV], as
shown in the upper panel in Fig. 6. Again, it is due to the faster drop of the
trace of the energy-momentum tensor than the energy density. Those rapid drops
might be smeared to be not so much sharper than what we observe in Fig. 4, by
taking into account not only the conventional ensemble of single free quanta
but also interacting multi-body processes (as encoded in collision terms in
the Boltzmann equation) in evaluating $\Sigma$. Therefore, we claim that
emergence of the nonadiabaticity could be an artifact coming from the
invalidity of the conventional ensemble description.

We could qualitatively argue the validity of the conventional ensemble
description applied in the existing works. First of all, note that the matter
contents can be categorized into two, i.e., relativistic radiation and
nonrelativistic matter. When the temperature $T$ is lower than the mass $m$,
Eq. (18) suggests $\rho_{i}-3P_{i}\approx\rho\propto T^{3/2}e^{-m/T}$. Since
$\rho_{\text{r}}\propto T^{4}$, the kick function Eq. (20) in the thermal
equilibrium evolves as

$\displaystyle\Sigma(T)\propto T^{-5/2}e^{-m/T}\,.$ (28)

Thus, as the temperature decreases, the kick function exponentially decreases
to zero, which causes the strong nonadiabaticity in Eq. (27).

On the other hand, the above picture is invalid after the thermal decoupling
of SM particles, and the number density $n$ satisfies $n\propto a^{-3}$ to
produce the thermal relic abundance. Because the scale factor and temperature
are related to each other as $a\propto T^{-1}$, the energy density should
follow $\rho=mn\propto T^{3}$ after the decoupling. Thus, taking into account
the effect of decoupling, we find

$\displaystyle\Sigma(T)\propto T^{-1}\,.$ (29)

Therefore, the kick function no longer decreases monotonically, which can
prevent the scalaron from being nonadiabatic. It is remarkable that the
decoupling significantly affects the scalaron dynamics in the thermal history,
and it shows an adaptive limit of the conventional ensemble description in the
study of the chameleon mechanism. Note that the chameleon mechanism works to
stabilize the scalaron with the thermal relics of SM particles and to avoid
the nonadiabaticity.

##  V Conclusion and Discussions

The unification of DE and inflation in $F(R)$ gravity has been widely
investigated. This work has explored the scalaron evolution with a specific
potential structure combining the DE and inflationary epochs. The scalaron
should have been frozen around its de Sitter vacuum in the radiation-dominated
epoch due to Hubble friction. However, as shown in this work, the decoupling
of SM particles from the thermal bath would drive the scalaron to evolve and
climb its effective potential. In particular, the contribution from the QCD PT
has a nontrivial impact on the scalaron dynamics. We found that the kicks
with/without the QCD PT would pull the scalaron up to $0.64\ M_{\text{Pl}}$
and $1.41\ M_{\text{Pl}}$, respectively.

As a result, successful BBN is preserved as long as one chooses the proper
initial conditions of the scalaron field at the beginning of the radiation-
dominated epoch. Unlike usual chameleon models, due to the scale hierarchy of
inflation and DE, we clarified that the $R^{2}$ term protects the scalaron
mass from receiving nonadiabatic change in time. On the other hand, the
nonadiabaticity relies on the choice of the initial value, which is determined
by the inflationary epoch. Therefore, to estimate the impact of the plain part
of the potential on the scalaron dynamics, one has to study the inflation,
reheating, and radiation-dominated epochs systematically and inclusively.
Studying such intermediate stages is crucial in assessing the validity of
$F(R)$ gravity, and we will leave the complete analysis in future work.

In closing, we make possible prospects for the post/onset-inflationary epoch
and particle production to supply the reheating of the Universe. In the
current work, we have discussed the impact of the $R^{2}$ term at the
radiation-dominated epoch, whose potential is so steep that the scalaron
hardly climbs up even under the most potent QCDPT kick. Although we
demonstrated the numerical analysis for four different initial conditions
chosen by hand, the ambiguity of the initial condition can be embedded into
the inflation dynamics if we consider the whole cosmic history. Thus, it is
inevitable to investigate the scalaron dynamics between the inflation and
radiation-dominated epoch. We shall briefly discuss the connection of $R^{2}$
inflation to DE. Although the DE scale has little impact on the dynamics of
inflation Geng  _et al._ (2015), the post inflationary dynamics, i.e.
(p)reheating, would be significantly changed due to the presence of plain at
$\varphi<0$.

Assuming the cold inflation, the reheating after the inflation induces the
radiation-dominated epoch. In absence of matter fields during the inflation,
one observes that the equation of motion is reduced to

$\displaystyle\phi^{\prime\prime}=-3\left(1-\frac{1}{6}\phi^{\prime
2}\right)\left(\phi^{\prime}+\frac{dV}{Vd\phi}\right)\,.$ (30)

By solving the equation of motion, as shown in Fig 7, we find the scalaron
would quickly pass its potential minimum at $\varphi\simeq-0.11M_{\text{Pl}}$,
experiences one single rebound, and eventually oscillates at around the
effective potential minimum.

Figure 7:  The scalaron evolution in absence of matter fields, where
$\tilde{N}$ is the number of e-folds in the Einstein frame. In order for
enough e-folds during inflation, we have chosen $\phi_{\text{i}}=5.5$.

We also reevaluate the nonadiabaticity in the scalaron mass in Eq. (26) with
the matter contribution ignored. As in Fig. 8, we find that the
nonadiabaticity is triggered when the scalaron passes through its bare
potential minimum.

Figure 8:  The nonadiabaticity of the scalaron in the absence of matter
fields, depicted by the orange-dashed line. The blue-solid line corresponds to
the field evolution (the same as Fig. 7), while the black-dotted line marks
the corresponding e-folding number when the scalaron reaches the bare
potential minimum.

Note that $V^{\prime\prime}(\varphi)$ is negative during the slow-roll
inflation, and we cannot define the mass in Eq. (26).

The strong nonadiabaticity around the bare potential minimum originates from
the hierarchy between the inflation and DE scales: that is, the role-shift of
the scalaron from the inflaton to dynamical DE triggers the huge change in
time for the scalaron mass from the inflationary epoch to the post-
inflationary epoch. The nonadiabaticity indicates the nonperturbative self-
production of the scalaron quanta and the breakdown of classical analysis
Erickcek  _et al._ (2013, 2014). Consequently, compared with the original
$R^{2}$ inflation, the perturbative oscillation picture around $\varphi=0$
disappears, thus, in the unified $R^{2}$-corrected class of models in $F(R)$
gravity, we need to alter the post-inflationary dynamics of the standard
cosmology. This consequence relies only on the hierarchy between the inflation
and DE scales, irrespective of the details of DE models Motohashi and
Nishizawa (2012); Nishizawa and Motohashi (2014).

Therefore, refining the early cosmology in $R^{2}$-corrected DE models of
$F(R)$ gravity is mandatory. With the results based on the current study, the
following approaches would be plausible:

  * •

In light of the nonadiabaticity of the scalaron, the effective scalaron mass
would presumably include correlation functions of the scalaron quanta Erickcek
_et al._ (2013, 2014). As in Fig. 8, the scalaron self-production due to the
large nonadiabaticity could happen when the scalaron rolls the potential, but
before it passes through the DE vacuum and rebounds ($\tilde{N}\lesssim
67.5$). Then, the produced and accumulated scalaron quanta can highly dominate
in the effective potential over the bare potential of the DE scale. Then it
acts as the potential barrier to push the scalaron back to $\varphi=0$. In
other words, the energy density of the scalaron quanta serves as the chameleon
mechanism. It lifts the potential, which prevents the scalaron from going to
the bare potential minimum. This barrier would instantaneously be present and
gone when the nonperturbative-scalaron quanta-production ends. Then the
background scalaron could get escaped back to the standard DE vacuum. Such a
self-chameleon mechanism could trigger the nontrivial scalaron dynamics after
the inflation. Moreover, when the nonperturbative production is included, the
e-folding $\tilde{N}$ could be also affected. It is intriguing to study the
possible scalaron oscillation around $\varphi=0$ on the effective potential
and consecutive particle production and investigate how to connect to the
radiation-dominant epoch.

  * •

The motion of SM particles (number densities) coupled to scalaron can be
nonadiabatic when $\varphi\rightarrow+\infty$ due to the exponential coupling
form, $\exp\left(-\sqrt{1/6}\kappa\varphi\right)$. This implies that
nonperturbative production of the SM particles could happen during or
immediately after the $R^{2}$ inflation, in the manner called the preheating
Dolgov and Kirilova (1990); Traschen and Brandenberger (1990); Kofman  _et
al._ (1994); Shtanov  _et al._ (1995); Kofman  _et al._ (1997) (for reviews,
see e.g., Kofman  _et al._ (1997); Amin  _et al._ (2014); Lozanov (2019)).
This also tells us a possibility that the chameleon mechanism works even right
after the inflation is over, hence the scalaron might not smoothly move down
to the DE vacuum. This scenario might be similar to the situation known as the
warm inflation Berera (1995) instead of the standard cold inflation, the model
parameters can thereby severely be constrained by cosmological observations.
It is worth pursuing this scalaron preheating along with the SM particle
production in detail.

At any rate, the prospected scenarios above shall fill the missing link in the
$F(R)$ gravity cosmology, which deserves in other publications.

###### Acknowledgements.

We are grateful to Seishi Enomoto for the valuable discussion. H.C.
appreciates Taotao Qiu’s insightful comments. T.K. is supported by the
National Key R&D Program of China (2021YFA0718500). S.M. is supported by the
National Science Foundation of China (NSFC) under Grant No.11975108, 12047569,
12147217 and the Seeds Funding of Jilin University.

##  Appendix A Potential Analysis for Numerical Analysis

Despite the successful connection between inflation and DE, we face the
hierarchy problem numerically. The smallness of the DE scale suggests the
gigantic separation from any other physical scale; however, we are interested
in physics at such an intermediate scale. For $F(R)$ model (16), we examine
the scalaron field $\varphi$ in Eq. (4) and scalaron potential in Eq. (6) in
terms of Ricci scalar $R$:

$\displaystyle e^{2\sqrt{1/6}\kappa\varphi}$
$\displaystyle=F^{\prime}_{\text{DE}}(R)+2\alpha R$ (31) $\displaystyle
V\left(\varphi(R)\right)$
$\displaystyle=\frac{1}{2\kappa^{2}}\frac{\left[RF^{\prime}_{\text{DE}}(R)-F_{\text{DE}}(R)\right]+\alpha
R^{2}}{\left[F^{\prime}_{\text{DE}}(R)+2\alpha R\right]^{2}}\,.$ (32)

The intermediate scale indicates $R_{c}\ll R$, and in Starobinsky’s DE model

$\displaystyle F_{\text{DE}}(R)=R-\beta
R_{\text{c}}\left[1-\left(1+\frac{R^{2}}{R_{\text{c}}^{2}}\right)^{-n}\right]\,,$
(33)

one finds

$\displaystyle-\beta
R_{\text{c}}\left[1-\left(1+\frac{R^{2}}{R_{\text{c}}^{2}}\right)^{-n}\right]\approx$
$\displaystyle-\beta R_{\text{c}}\,.$ (34)

And thus, the total $F(R)$ function is approximated as GR with the
cosmological constant and $R^{2}$ correction:

$\displaystyle F(R)=R-\beta R_{\text{c}}+\alpha R^{2}\,.$ (35)

The above is common even in other viable DE models of $F(R)$ gravity so that
DE models of $F(R)$ gravity reproduce the $\Lambda$-CDM model in the large-
curvature limit. Thus, using the above $F(R)$ model for the large-curvature
region $R_{c}\ll R$, we can investigate the $R^{2}$-corrected DE model or the
unification of the inflation and DE in a model-independent way.

For the model (35), the scalaron field and potential in terms of the Ricci
scalar are given as

$\displaystyle e^{2\sqrt{1/6}\kappa\varphi}$ $\displaystyle=1+2\alpha R\,,$
(36) $\displaystyle V\left(\varphi(R)\right)$
$\displaystyle=\frac{1}{2\kappa^{2}}\frac{\beta R_{\text{c}}+\alpha
R^{2}}{\left[1+2\alpha R\right]^{2}}\,.$ (37)

Figs. 9 and 10 shows the scalaron field $\varphi$ and potential $V(\varphi)$
in terms of the Ricci scalar $R$.

Figure 9:  Comparison among the $R^{2}$-corrected Starobinsky model (gray
solid) in Eq. (17), the Starobinsky model (blue dashed) in Eq. (33),
approximated model for intermediate scale (green dot-dashed) in Eq. (35), and
the $R^{2}$ inflation model (red dotted) with respect to the scalaron fields
$\varphi$ as the function of the Ricci scalar $R$. The parameters are chosen
as $n=1$, $\beta=2$, $\alpha R_{c}=10^{-10}$.  Figure 10:  The same as Fig. 9
with respect to the scalaron potential $V(\varphi)$. Eq. (37) suggests that DE
term and $R^{2}$ correction are comparable when $R/R_{c}=\sqrt{\beta/\alpha
R_{c}}\sim 10^{-5}$ in this setup.

One can find that the approximated model for the intermediate scale is
consistent with the full $R^{2}$-corrected one for the range $\varphi\geq 0$,
and it is not valid for the range $\varphi<0$ (see Fig. 1). Note that
$\varphi\rightarrow 0$ and $V(\varphi)\rightarrow\text{const.}$ when
$R\rightarrow\infty$, which shows the singularity problem Frolov (2008).

Furthermore, one finds an analytic form of the scalaron potential in terms of
the scalaron field,

$\displaystyle V(\varphi)$ $\displaystyle=\frac{\beta
R_{\text{c}}}{2\kappa^{2}}e^{-4\sqrt{1/6}\kappa\varphi}+\frac{1}{8\kappa^{2}\alpha}\left(1-e^{-2\sqrt{1/6}\kappa\varphi}\right)^{2}\,.$
(38)

Note that the second term in the potential $V(\phi)$ represents the $R^{2}$
inflation potential. Thus, the potential for $\varphi\geq 0$ is approximately
described by Eq. (38) (see, Fig. 1). Although the potential for $\varphi<0$
should be written by the Starobinsky model (33), the matter contribution
$T^{\mu}_{\ \mu}$ in the effective potential Eq. (11) overwhelms the bare
potential because of the hierarchy between the DE scale and the other physical
scale. Finally, one can conclude that in the existence of matter and chameleon
mechanism, Eq. (35) and the trace of the energy-momentum tensor determines the
scalaron dynamics in the $R^{2}$-corrected DE models of $F(R)$ gravity.

##  Appendix B Construction of kick function

In this section, we show how to handle data from Refs. Husdal (2016); Bali
_et al._ (2014). Data regarding the degrees of freedom for energy density
$\rho$ and entropy density $s$ together with temperature have been listed in
Table A1 in Ref. Husdal (2016). In order to study the effect of QCD PT,
however, we also use the result on the equation of state from the lattice QCD
simulation Bali  _et al._ (2014), ranging from $T=114-305$ [MeV], where data
regarding entropy density $s$, energy density $\rho$, and pressure $P$ of
quarks and mesons plasma have been provided 222 In the literature Bali  _et
al._ (2014) the lattice simulation has also incorporated nonzero magnetic
fields so that the system gets anisotropic, hence the pressure $P$ is actually
measured along the $z$-direction, $P_{z}$. The data plots that we quote in the
present study, include the cases with and without the applied magnetic field,
from which we have taken only the case with zero magnetic field. . We extract
the degrees of freedom for entropy density from

$\displaystyle g_{*s}(T)\equiv\frac{45}{2\pi^{2}}\frac{s(T)}{T^{3}}\,,$ (39)

Furthermore, we have assumed that $g_{*}=g_{*s}$ during QCD PT, whose
differences only appear after electron-positron annihilation.

Then, we set $T=305$ [MeV] and $T=114$ [MeV] as dividing points to separately
deal with quantities before, in, and after QCD PT. In particular, we add
$14.25$, the remaining degree of freedom after the annihilation of $\pi^{0}$,
into that during QCD PT by hand. With the degree of freedom for the energy
density in hand, one obtains the energy density in Eq. (19) and thereby
obtains the kick function from Eq. (20). The tabulated data regarding Fig. 2
is available in `degrees_of_freedom.dat` on arXiv.

## References

  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Martin (2016) J. Martin, Astrophys. Space Sci. Proc. 45, 41 (2016), arXiv:1502.05733 [astro-ph.CO] .
  * Riess  _et al._ (1998) A. G. Riess _et al._ (Supernova Search Team), Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 .
  * Perlmutter  _et al._ (1999) S. Perlmutter _et al._ (Supernova Cosmology Project), Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 .
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Nojiri  _et al._ (2017) S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Phys. Rept. 692, 1 (2017), arXiv:1705.11098 [gr-qc] .
  * Heisenberg (2019) L. Heisenberg, Phys. Rept. 796, 1 (2019), arXiv:1807.01725 [gr-qc] .
  * De Felice and Tsujikawa (2010) A. De Felice and S. Tsujikawa, Living Rev. Rel. 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
  * Appleby and Battye (2008) S. A. Appleby and R. A. Battye, JCAP 05, 019 (2008), arXiv:0803.1081 [astro-ph] .
  * Frolov (2008) A. V. Frolov, Phys. Rev. Lett. 101, 061103 (2008), arXiv:0803.2500 [astro-ph] .
  * Dev  _et al._ (2008) A. Dev, D. Jain, S. Jhingan, S. Nojiri, M. Sami, and I. Thongkool, Phys. Rev. D 78, 083515 (2008), arXiv:0807.3445 [hep-th] .
  * Kobayashi and Maeda (2009) T. Kobayashi and K.-i. Maeda, Phys. Rev. D 79, 024009 (2009), arXiv:0810.5664 [astro-ph] .
  * Appleby  _et al._ (2010) S. A. Appleby, R. A. Battye, and A. A. Starobinsky, JCAP 06, 005 (2010), arXiv:0909.1737 [astro-ph.CO] .
  * Lee  _et al._ (2012) C.-C. Lee, C.-Q. Geng, and L. Yang, Prog. Theor. Phys. 128, 415 (2012), arXiv:1201.4546 [astro-ph.CO] .
  * Starobinsky (1980) A. A. Starobinsky, Phys. Lett. B 91, 99 (1980).
  * Odintsov and Oikonomou (2016) S. D. Odintsov and V. K. Oikonomou, Class. Quant. Grav. 33, 125029 (2016), arXiv:1602.03309 [gr-qc] .
  * Odintsov  _et al._ (2017) S. D. Odintsov, V. K. Oikonomou, and L. Sebastiani, Nucl. Phys. B 923, 608 (2017), arXiv:1708.08346 [gr-qc] .
  * Nojiri  _et al._ (2020) S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Phys. Dark Univ. 29, 100602 (2020), arXiv:1912.13128 [gr-qc] .
  * Oikonomou (2021a) V. K. Oikonomou, Phys. Rev. D 103, 044036 (2021a), arXiv:2012.00586 [astro-ph.CO] .
  * Odintsov and Oikonomou (2019) S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 99, 104070 (2019), arXiv:1905.03496 [gr-qc] .
  * Odintsov and Oikonomou (2020) S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 101, 044009 (2020), arXiv:2001.06830 [gr-qc] .
  * Odintsov and Oikonomou (2021) S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 104, 124065 (2021), arXiv:2112.06269 [gr-qc] .
  * Oikonomou and Giannakoudi (2022a) V. K. Oikonomou and I. Giannakoudi, (2022a), arXiv:2205.08599 [gr-qc] .
  * Artymowski and Lalak (2014) M. Artymowski and Z. Lalak, JCAP 09, 036 (2014), arXiv:1405.7818 [hep-th] .
  * Geng  _et al._ (2015) C.-Q. Geng, C.-C. Lee, and S. Lin, Astrophys. Space Sci. 360, 21 (2015), arXiv:1506.03253 [astro-ph.CO] .
  * Yashiki  _et al._ (2020) M. Yashiki, N. Sakai, and R. Saito, Phys. Rev. D 102, 043504 (2020), arXiv:2004.00864 [astro-ph.CO] .
  * Oikonomou (2021b) V. K. Oikonomou, Phys. Rev. D 103, 124028 (2021b), arXiv:2012.01312 [gr-qc] .
  * Oikonomou and Giannakoudi (2022b) V. K. Oikonomou and I. Giannakoudi, Nucl. Phys. B 978, 115779 (2022b), arXiv:2204.02454 [gr-qc] .
  * Erickcek  _et al._ (2013) A. L. Erickcek, N. Barnaby, C. Burrage, and Z. Huang, Phys. Rev. Lett. 110, 171101 (2013), arXiv:1304.0009 [astro-ph.CO] .
  * Erickcek  _et al._ (2014) A. L. Erickcek, N. Barnaby, C. Burrage, and Z. Huang, Phys. Rev. D 89, 084074 (2014), arXiv:1310.5149 [astro-ph.CO] .
  * Bali  _et al._ (2014) G. S. Bali, F. Bruckmann, G. Endrödi, S. D. Katz, and A. Schäfer, JHEP 08, 177 (2014), arXiv:1406.0269 [hep-lat] .
  * Katsuragawa  _et al._ (2019) T. Katsuragawa, T. Nakamura, T. Ikeda, and S. Capozziello, Phys. Rev. D 99, 124050 (2019), arXiv:1902.02494 [gr-qc] .
  * Cembranos (2009) J. A. R. Cembranos, Phys. Rev. Lett. 102, 141301 (2009), arXiv:0809.1653 [hep-ph] .
  * Husdal (2016) L. Husdal, Galaxies 4, 78 (2016), arXiv:1609.04979 [astro-ph.CO] .
  * Chen  _et al._ (2020) H. Chen, T. Katsuragawa, S. Matsuzaki, and T. Qiu, JHEP 02, 155 (2020), arXiv:1908.04146 [hep-ph] .
  * Motohashi and Nishizawa (2012) H. Motohashi and A. Nishizawa, Phys. Rev. D 86, 083514 (2012), arXiv:1204.1472 [astro-ph.CO] .
  * Nishizawa and Motohashi (2014) A. Nishizawa and H. Motohashi, Phys. Rev. D 89, 063541 (2014), arXiv:1401.1023 [astro-ph.CO] .
  * Dolgov and Kirilova (1990) A. D. Dolgov and D. P. Kirilova, Sov. J. Nucl. Phys. 51, 172 (1990).
  * Traschen and Brandenberger (1990) J. H. Traschen and R. H. Brandenberger, Phys. Rev. D 42, 2491 (1990).
  * Kofman  _et al._ (1994) L. Kofman, A. D. Linde, and A. A. Starobinsky, Phys. Rev. Lett. 73, 3195 (1994), arXiv:hep-th/9405187 .
  * Shtanov  _et al._ (1995) Y. Shtanov, J. H. Traschen, and R. H. Brandenberger, Phys. Rev. D 51, 5438 (1995), arXiv:hep-ph/9407247 .
  * Kofman  _et al._ (1997) L. Kofman, A. D. Linde, and A. A. Starobinsky, Phys. Rev. D 56, 3258 (1997), arXiv:hep-ph/9704452 .
  * Amin  _et al._ (2014) M. A. Amin, M. P. Hertzberg, D. I. Kaiser, and J. Karouby, Int. J. Mod. Phys. D 24, 1530003 (2014), arXiv:1410.3808 [hep-ph] .
  * Lozanov (2019) K. D. Lozanov, (2019), arXiv:1907.04402 [astro-ph.CO] .
  * Berera (1995) A. Berera, Phys. Rev. Lett. 75, 3218 (1995), arXiv:astro-ph/9509049 .