Towards a unified interpretation of the early Universe in R²-corrected dark energy model of F(R) gravity

Growing observational evidence indicates two accelerations of cosmic expansion: the early-time inflation expected to occur before radiation domination [1, 2], and the late-time acceleration sourced by dark energy (DE) [3,4] (See [5] 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 [6, 7]. 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) $ [8]. 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 [9, 10]. A promising way to cure this problem is to add higher-curvature corrections, for instance, the $ \alpha R^2 $ term [11–14]. Moreover, $ R^2 $ gravity is also well known as the inflation model [15], which has proven very successful in light of the cosmic microwave background observation [1]. 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 [16–23]. Other setups have also been considered, for instance, $ R+\alpha R^{n} -\beta R^{2-n} $ [24] and its constraints [25, 26], and a single scalar field extension to $ F(R) $ gravity [27, 28].

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 [8]. 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 [29, 30].

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 [29, 30]. Using data on the temperature dependence of the equation of state observed as the crossover reported from the nonperturbative lattice QCD simulation [31], 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}={\rm d}\varphi/{\rm d}\tilde{N}$, where $ \tilde{N} $ is the number of e-folds with "tilde" denoting a physical quantity in the Einstein frame.

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

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

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

Under the Weyl transformation,

this field φ, dubbed scalaron, appears explicitly in the Einstein-frame action,

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

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

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

and the acceleration equation,

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

from which we define the effective potential ignoring the φ-dependence in the matter Lagrangian

Then the mass of the scalaron is defined as

For later convenience, we define some dimensionless quantities: the density ratio of cold matter $ \tilde{\rho}_{\rm{m}} $ and radiation $ \tilde{\rho}_{\rm{r}} $, i.e. $ f_{\rm{m}}\equiv\tilde{\rho}_{\rm{m}}/\tilde{\rho}_{\rm{r}} $, the kick function $ \Sigma\equiv(\tilde{\rho}_{{\rm{r}}}- 3\tilde{P}_{{\rm{r}}})/\tilde{\rho}_{{\rm{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

and

where $\phi^\prime={\rm d}\phi/{\rm d}\tilde{N}$. Since in the radiation dominated epoch, $ f_{\rm{m}}\lesssim 10^{-6} $ [30], we shall hereafter ignore $ f_{\rm{m}} $. Finally, the scalaron velocity with respect to the cosmic time can be easily read off as

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

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

where $ \beta > 0 $, $ n > 0 $, and $ R_{\rm{c}}\sim\Lambda\simeq4\times10^{-84}\ [\mathrm{GeV}^{2}] $ are model parameters. The potential is a multi-valued function (see, e.g., Refs. [10, 32]). In the present study we focus on the case where $ +\infty>R>R_{\rm{b}} $ with the branch located at $ R_{\rm{b}} \simeq 0.93 R_{\rm{c}} $ for $ \beta=2 $ and $ n=1 $. The fifth force experiment constrains $ \alpha<10^{22}\ [\mathrm{GeV}^{-2}] $ [33], and $\alpha = 1/(6 m_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. On the contrary to $ R^{2} $ inflation, the potential possesses a shallow plain with a global minimum at $ \varphi \simeq 0.11 M_{\rm{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.

Figure 1.  (color online) 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_{{\rm{c}}}}{2\kappa^{2}}\sim\rho_{\Lambda}\approx10^{-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.

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 Λ-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 α, 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.

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_{{\rm{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

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_{\rm{r}} $ is obtained by

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

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

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}=\dfrac{2\pi^{2}}{45}g_{*s}(T)T^{3}a^{3}= {\rm{const.}} $, one can extract $ \Sigma (\tilde{N}) $ by making a replacement

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 [29, 30]. 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

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_{\rm{r}} $ and ${\rm d} V/{\rm d}\phi$ is negligible compared to the driving term in Σ in the early Universe. Then Eq. (14) is reduced to

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

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 [29, 30] 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.

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. [34] together with the result on the equation of state from the lattice QCD simulation [31], 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.  (color online) 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.  (color online) 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_{\rm{Pl}} $ and $ 1.41\ M_{\rm{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{\dfrac{2+\Sigma(T_{{\rm{s}}})}{6 \Sigma(T_{{\rm{s}}})}}\gtrsim 0.85 $, and the velocity can be as large as $ \phi'=1/Q\gtrsim 1.17 $, which exacerbates the existing catastrophic results claimed in [29, 30].

Figure 5.  (color online) 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_{\rm{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 φ 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 [5], 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_{{\rm{BBN}}}\lesssim 0.15 $, which requires the initial value as $ \varphi_i\gtrsim -0.97\ M_{\rm{Pl}} $ or $ \varphi_{i}\gtrsim -1.74\ M_{\rm{Pl}} $ under the SMPE kicks or the total five kicks. More stringent bound as evaluated in [35] might be placed, where $-0.05\lesssim \kappa\varphi_{{\rm{BBN}}}\lesssim0.03$. In that case, we would have $ \varphi_i\gtrsim -0.69\ M_{\rm{Pl}} $ with only the four SMPE kicks taken into account, or $ \varphi_{i}\gtrsim -1.46\ M_{\rm{Pl}} $ with the total kicks considered^①. Thus, the QCDPT kick is also significant in relaxing the initial value of the scalaron field.

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

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

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

Figure 6.  (color online) 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.

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

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 Σ. 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} {\rm e}^{-m/T}$. Since $ \rho_{\rm{r}} \propto T^{4} $, the kick function Eq. (20) in the thermal equilibrium evolves as

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 = m n \propto T^{3} $ after the decoupling. Thus, taking into account the effect of decoupling, we find

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.

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_{\rm{Pl}} $ and $ 1.41\ M_{\rm{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 [23], 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

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_{\rm{Pl}} $, experiences one single rebound, and eventually oscillates at around the effective potential minimum. 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. Note that $ V^{\prime \prime}(\varphi) $ is negative during the slow-roll inflation, and we cannot define the mass in Eq. (26).

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_{\rm{i}}=5.5 $.

Figure 8.  (color online) 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.

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 [29, 30]. 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 [36, 37].

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 [29, 30]. 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 [38–42] (for reviews, see e.g., [42–44]). 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 [45] 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.

We are grateful to Seishi Enomoto for the valuable discussion. H.C. appreciates Taotao Qiu's insightful comments.

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 φ in Eq. (4) and scalaron potential in Eq. (6) in terms of Ricci scalar R:

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

one finds

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

The above is common even in other viable DE models of $ F(R) $ gravity so that DE models of $ F(R) $ gravity reproduce the Λ-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 (A5), the scalaron field and potential in terms of the Ricci scalar are given as

Figures. A1 and A2 shows the scalaron field φ and potential $ V(\varphi) $ in terms of the Ricci scalar R. 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 {\rm{const.}} $ when $ R\rightarrow \infty $, which shows the singularity problem [10].

Figure A1.  (color online) Comparison among the $ R^2 $-corrected Starobinsky model (gray solid) in Eq. (17), the Starobinsky model (blue dashed) in Eq. (A3), approximated model for intermediate scale (green dot-dashed) in Eq. (A5), and the $ R^{2} $ inflation model (red dotted) with respect to the scalaron fields φ as the function of the Ricci scalar R. The parameters are chosen as $ n=1 $, $ \beta=2 $, $ \alpha R_{c}=10^{-10} $.

Figure A2.  (color online) The same as Fig. A1 with respect to the scalaron potential $ V(\varphi) $. Eq. (A7) 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.

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

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. (A8) (see, Fig. 1). Although the potential for $ \varphi<0 $ should be written by the Starobinsky model (A3), 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. (A5) and the trace of the energy-momentum tensor determines the scalaron dynamics in the $ R^{2} $-corrected DE models of $ F(R) $ gravity.

In this section, we show how to handle data from Refs. [31, 34]. Data regarding the degrees of freedom for energy density ρ and entropy density s together with temperature have been listed in Table A1 in Ref. [34]. 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 [31], ranging from $ T=114-305 $ MeV, where data regarding entropy density s, energy density ρ, and pressure P of quarks and mesons plasma have been provided^②. We extract the degrees of freedom for entropy density from

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.