[**Seokcheon Lee$^{\,\dagger}$**]{}

$^\dagger$[*Research institute of natural science, Gyeongsang national university, 501 Jinju Daero, Jinju city, 52828, Korea*]{}

0.6in

**Abstract**

Starobinsky has suggested an inflation model which is obtained from the vacuum Einstein’s equations modified by the one-loop corrections due to quantized matter fields. Although the one-loop gravitational action is not known for a general FRW background, it can be obtained in a de Sitter space to give $\Mp^2 R + \alpha R^2 + \beta R^2 \ln (R/M^2)$. Thus, one needs to investigate the inflationary behavior of this model compared to the Starobinsky model ([*i.e.*]{} $\beta = 0$). The coefficient $\alpha$ can be changed by varying the renormalization scale $M^2$ and $\beta$ is obtained from the quantum anomaly which is related to the numbers of quantum fields. It has been assumed that $\alpha \gg \beta$. We investigate the viable values of $\alpha$ and $\beta$ based on the CMB observation. We also scrutinize the reheating process in this model.

Introduction
============

In semi-classical gravity, one treats the background as classical, but with taking the back-reaction of the matter into account. In a curved space-time, even in the absence of classical matter or radiation, quantum fluctuations of matter fields give nontrivial contribution to the expectation value of the energy-momentum tensor being quantized in some state [@BD; @9509057; @9707062; @BF; @PT]. $$\tensor{R}{_\mu_\nu} -\frac{1}{2} \tensor{g}{_\mu_\nu} R = 8 \pi G \langle \tensor{\hat{T}}{_\mu_\nu} \rangle \, ,  \label{1loop}$$ where $\langle \tensor{\hat{T}}{_\mu_\nu} \rangle$ includes $\infty^{4}$, $\infty^{2}$, and $\ln \infty$ divergences. The first and the second divergences can be eliminated by the renormalization of the cosmological constant and the Newtonian coupling constant, respectively. However, the logarithmic divergence can be eliminated only if one introduces an additional Lagrangian density of the form $\sqrt{-g} \left(c_1 R^2 + c_2 \tensor{R}{_\mu_\nu} \tensor{R}{^\mu^\nu} \right)$ [@Nariai73].

Among these quantum corrections to general relativity, if one includes only the quadratic terms in the Riemann tensors, then these terms are expressed as [@Mijic86] $$\alpha R^2 + \beta \tensor{R}{_\mu_\nu} \tensor{R}{^\mu^\nu} + \gamma \tensor{R}{_\mu_\nu_\lambda_\sigma} \tensor{R}{^\mu^\nu^\lambda^\sigma} = \epsilon R^2 + \xi \tensor{C}{_\mu_\nu_\lambda_\sigma} \tensor{C}{^\mu^\nu^\lambda^\sigma} + \eta \chi_{\text{E}} \label{quadraticR} \, ,$$ where $C$ is the Weyl tensor and $\chi_{\text{E}}$ is the density of the Euler number for the manifold. If one considers the Robertson-Walker metric, then $C^2$-term vanishes because the metric is conformally flat. Thus, Starobinsky includes the curvature-squared correction to the Einstein-Hilbert action and this term is important for the early universe [@Star80]. At the early universe ([*i.e.*]{} when the curvatures are large), the curvature squared term is dominant than the Ricci scalar term and this correction term leads to an effective cosmological constant. Therefore, the early universe went through an inflationary de Sitter era without introducing an inflaton field.

These quantum corrections take a simple form in the case of free, massless, conformally invariant fields. The contributions which all fields make to the geometrical part of the vacuum to vacuum amplitude. These contributions originate in the vacuum polarization, which the background geometry induces, and give rise to nonobservable renormalizations as well as physically real radiative corrections [@DeWitt67] $${\cal L}_{\text{vp}} = \beta R^2 \ln \left(\frac{R}{M^2} \right) \, . \label{Lvacuumpol}$$ Thus, when one regards the origin of the Starobinsky model which includes a quadratic Ricci scalar term in addition to the Hilbert-Einstein action, it is also natural to include the terms quadratic and logarithmic in the Ricci scalar. There have been investigation of the $f(R)$-gravity including the logarithmic correction for several aspects. The viable f(R)-gravity models replacing the dark energy have been considered [@0308176; @0309062; @0512109; @07053199; @07081482; @08031081; @09125474; @13057290]. The logarithmic f(R) models are also considered to investigate the relativistic stars [@13077977]. Applying the logarithmic f(R)-gravity models for the early inflation have been scrutinized [@Vilenkin85; @Baibosunov90; @14047349; @14061096; @14100631; @14116010; @151209324].

Any regular model Universe must be non-singular for an infinitesimal perturbation of the metric $${\cal L}_{\text{reg}} = - \frac{M^{2}}{6} \left( \frac{R}{M^2} \right)^n \, . \label{Lregular}$$ The possible regular solutions of the gravitational equations in the presence of a nonlinear increment of the four curvature to the Lagrangian density of the gravitational field can be written as [@BGS71] $$S_{\text{BGS}} = \frac{\Mp^2}{16 \pi} \int \sqrt{-g} \, d^4 x \left( R + M^2 f_{\text{c}}\left[ \frac{R}{M^2} \right] \right) \, , \label{BGSfR}$$ where $f_{\text{c}}(R/M^2)$ is a certain dimensionless function of the scalar curvature and $M$ is the characteristic mass. This work is motivated to obtain the bouncing Universe at $\tau = a t = 0$.

As the minimal extension of the Starobinsky model, the one with the logarithmic correction term is given by [@Starobinsky79; @Shore80; @Vilenkin85; @Baibosunov90; @14047349] $$\begin{aligned}
 S &= \frac{\Mp^2}{2} \int \sqrt{-g} \, d^4 x \left( R + \tilde{\alpha} R^2 + \tilde{\beta} R^2 \ln R \right) \nonumber \\ 
&=  \frac{\Mp^2}{2} \int \sqrt{-g} \, d^4 x R \left(1 + \alpha \frac{R}{M^2} + \beta \frac{R}{M^2} \ln \frac{R}{M^2} \right) \, , \label{SR2lnR} \end{aligned}$$ where $\alpha = \tilde{\alpha} + 2 \tilde{\beta} \ln M$ and $\beta = \tilde{\beta}$. This action is also naturally obtained when one considers the small correction in the Starobinsky model $ R + \alpha R^2 + R^{2+ \gamma}$ with $\gamma \ll 1$ as shown in [@14116010]. The inflationary behavior of this action was nicely investigated analytically by using the approximation [@14047349]. We improve the previous work numerically without using any approximation. We also investigate the reheating process in this model. Similar forms of the logarithmic model based on the stability and the viability of observation are also considered [@151209324; @160508016]. Constraints on the model parameters given in this model should be obtained from the CMB observation [@150202114].

Higgs Starobinsky inflation models have been investigated to solve the large Higgs values in the early universe[@150908882; @160502236; @160906887; @161104932; @170106636; @170107665; @170505638; @170700984]. In these models, non-minimal coupling of the Higgs field to the Ricci scalar induces the large quantum corrections. These approaches might explain both the dark matter and the inflation.

In the usual inflation models, soon after the end of the inflation, the inflaton fields begin to oscillate around the minimum of their effective potentials, producing particles, which interact with each other to reach the thermal equilibrium at the reheating temperature, $T_{\text{rh}}$ [@10012600]. The reheating dynamics after the inflation induced by Starobinsky model have been investigated both in Einstein frame and in Jordan frame [@09091737; @12041472; @12124466; @160604346].

The layout of this manuscript proceeds as follows. In the next section, we briefly review the equation of motions of the $R + \alpha R^2 + \beta R^2 \ln R$ inflation model both in Jordan frame and in Einstein frame. One needs to be careful when $\beta$ is negative. The approximate analysis of the inflationary behaviors was investigated in the reference [@14047349]. We scrutinize the exact inflation in this model for the different values of $\alpha$ and $\beta$ to be consistent with the Planck results. We investigate the oscillatory epoch and the reheating period of this model in section 3. Finally, we conclude in section 4.

$f(R)$ Inflation
================

$f(R)$ inflation models are conventionally analyzed in the Einstein frame after the conformal transformation of the original action in the Jordan frame. In this section, we briefly review this process and provide the analytic relations formulae between the general form of $f(R)$ and the cosmological parameters obtained from the observations. We also investigate the specific form of $f(R)$ model obtained from the quantum correction of the matter to constrain its parameters.

General formalism
-----------------

The action for the general $f(R)$-gravity theories in 4 dimensional space-time can be rewritten as $$S \equiv \frac{1}{16 \pi G} \int \sqrt{-g} \, d^4 x f(R) \equiv \frac{\Mp^2}{2} \int \sqrt{-g} \, d^4 x M^2 \bar{f}(y) \, , \label{Sfx}$$ where $\Mp$ is the reduced Planck mass, $M$ is the characteristic mass scale, $y = R/M^2$ is the normalized dimensionless variable, and $\bar{f}$ is a dimensionless function. The above Jordan frame (JF) action is conformally transformed into the Einstein frame (EF) action by doing the Weyl transformation of the metric and the field redefinition $$\begin{aligned}
 S_{\text{E}} &= \int  \sqrt{-g^{\text{E}} } \, d^4 x \, \left( \frac{\Mp^2 }{2} R_{\text{E}} - \frac{1}{2} \partial^{\mu} \phi \partial_{\mu} \phi - V \left( \phi \right) \right) \label{ESR2lnR} \, , \\
& \text{where}  \,\,\,\,\,\, g^{(\text{E})}_{\mu \nu} \equiv \frac{\partial f}{\partial R} g_{\mu\nu} \, , \label{gE} \\ 
\phi  &\equiv \Mp \bphi = \sqrt{\frac{3}{2}} \Mp \ln \left[ \, \frac{\partial f }{\partial R} \, \right]  = \sqrt{\frac{3}{2}} \Mp \ln \left[ \, \frac{\partial \bar{f} }{\partial y} \, \right] \equiv \sqrt{\frac{3}{2}} \Mp \ln \bar{f}_{,y} \left( y \right) \, , \label{bphi} \\ 
V \left( \phi \right) &\equiv \frac{\Mp^2}{2} \frac{R \left( \partial f / \partial R \right) - f}{ \left( \partial f / \partial R \right)^2}  = \frac{\Mp^2 M^2}{2} \frac{y \bar{f}_{,y} - \bar{f}}{ \bar{f}_{,y}^2} \, . \label{V} \end{aligned}$$ The potential of the inflaton field, $\phi$ is obtained from the specific form of $f(R)$ from the theory. Thus, one can analyze the usual inflationary behavior of this model in the EF.

One can obtain the slow-roll parameters as a function of $\bar{f}$ by using the Eq.(\[V\]) $$\begin{aligned}
 \epsilon_{\text{V}} &\equiv \frac{\Mp^2}{2} \left( \frac{V_{,\phi}}{V} \right)^{2} = \frac{\left(-2 \bff+y \bff_{,y} \right)^2}{3 \left(\bff - y \bff_{,y}\right)^2}\label{epsilon} \, ,  \\
\eta_{\text{V}} &\equiv \Mp^2 \frac{V_{,\phi\,\phi}}{V} = -\frac{2 \left(\bff_{,y}^2-4 \bff \bff_{,yy} +y \bff_{,y} \bff_{,yy} \right)}{3 \left( \bff -y \bff_{,y} \right) \bff_{,yy}}\label{eta} \, , \\
\xi_{\text{V}}^{2} &\equiv 2 \Mp^4 \frac{V_{,\phi} V_{,\phi\phi\phi}}{V^2} = \frac{8 \left(2 \bff -y \bff_{,y}\right) \left(3\bff_{,y}^2 \bff_{,yy}^2-8 \bff \bff_{,yy}^3+y \bff_{,y} \bff_{,yy}^3+\bff_{,y}^3 \bff_{,yyy}\right)}{9 \left(\bff-y \bff_{,y}\right)^2 \bff_{,yy}^3}\label{xi} \, , \end{aligned}$$ where $V_{,\phi}$ means the derivative of $V$ with respect to the inflaton field, $\phi$. The second expressions in Eqs.(\[epsilon\])-(\[xi\]) are described by the dimensionless quantities and one can use these formulae without any mistake in units. The scalar spectral index $n_{s}$, the tensor-to-scalar ratio $r$, the running of the spectral index $d n_{s} /d \ln k$ and the scalar primordial amplitude $A_{s}$ are obtained from the above slow-roll parameters $$\begin{aligned}
 n_{s} &= 1 - 6 \epsilon_{\text{V}} + 2 \eta_{\text{V}}  \, , \label{ns} \\
r &= 16 \epsilon_{\text{V}} \label{r} \, , \\
\frac{d n_{s}}{d \ln k} &= 16 \epsilon_{\text{V}} \eta_{\text{V}} -24 \epsilon_{\text{V}}^2 -2 \xi_{\text{V}}^{2} \label{runns} \, , \\
A_{s} &= \frac{32}{75 \Mp^4} \frac{V}{\epsilon_{\text{V}}} \label{As} \, . \end{aligned}$$

Thus, if one solves the evolution of $f(R)$, then one can obtain the model predictions for the inflationary evolution. In order to obtain the numerical analysis for these predictions, one needs to solve the field equation obtained from Eq.(\[ESR2lnR\]) $$\ddot{\phi} + 3H \dot{\phi} + V_{,\phi} = 0 \label{phiEq} \, ,$$ where dot denotes the derivative with respect to the time, $t$. If one adopts the normalized field definition $\bMp \bphi = \phi$ and changes variable from $t$ to the e-folding number $N = \ln a$, then the above field equation becomes $$\bphi'' +  \left( 6 - \bphi^{'2} \right) \frac{\bphi'}{2} + \left( 6 - \bphi^{'2} \right) \frac{V_{,\bphi}}{2V} = 0 \label{bphiEq} \, .$$ The evolution of the given $f(R)$ as a function of $N$ is obtained when one inserts Eqs.(\[bphi\])-(\[V\]) into the Eq.(\[bphiEq\]). We investigate the specific model of $f(R)$ in the following subsection.

$R + \alpha \frac{R^2}{M^2} + \beta \frac{R^2}{M^2} \ln \frac{R}{M^2}$ model
----------------------------------------------------------------------------

As we mentioned in the introduction, the one loop corrections involving the matter fields in the curved space-time generally provide the correction terms with the quadratic and logarithmic in the Ricci scalar [@Nariai73; @DeWitt67]. This action is given in Eq.(\[SR2lnR\]), one can rewrite this as $$\begin{aligned}
 S &= \frac{\Mp^2}{2} \int \sqrt{-g} \, d^4 x R \left(1 + \alpha \frac{R}{M^2} + \beta \frac{R}{M^2} \ln \frac{R}{M^2} \right) \equiv \frac{\Mp^2}{2} \int \sqrt{-g} \, d^4 x M^2 \bar{f}(y) \nonumber \\ 
&= \frac{\bMp^2 M^2}{2} \int \sqrt{-g} \, d^4 x \, y \left(1 + \alpha y + \beta y \ln y \right) \, , \label{Sx2lnx} \end{aligned}$$ where $M$ is the characteristic mass scale can be tuned in order to obtain the proper inflationary behavior, $y = R/M^2$ is the normalized dimensionless variable (curvature), and $\bar{f}$ is a dimensionless function. For the given $f(R)$, one can obtain the analytic form of the Ricci scalar, and the scalar potential, $V(\phi)$ as nicely shown in the reference [@14047349] $$\begin{aligned}
 R &= \frac{e^{\sqrt{\frac{2}{3}} \bphi} -1 }{2 \beta W_{k}(X)} = \frac{\bff_{,y} - 1}{2 \beta W_{k}(X)}  \label{RforR2} \,\,  , \\
X &\equiv \frac{e^{(2\alpha+\beta)/2\beta}}{2\beta} \left( e^{\sqrt{\frac{2}{3}} \bphi} -1 \right) = \frac{e^{(2\alpha+\beta)/2\beta}}{2\beta} \left( \bff_{,y} - 1 \right)   \label{XforR2} \,\, , \\
V(\bphi) &= \left(1 - e^{-\sqrt{\frac{2}{3}} \bphi} \right)^{2} \frac{1+2W_{k}(X)}{16 \beta W_{k}(X)^{2}} = \left( \frac{1 - \bff_{,y}}{\bff_{,y}} \right)^2  \frac{1+2W_{k}(X)}{16 \beta W_{k}(X)^{2}}  \label{VforR2} \, , \end{aligned}$$ where $W_{k}$ is the ProductLog (Lambert function) of branch $k=0$ for the positive value of $\beta$ and $k=-1$ for the negative one. Even though there exist these exact analytic solutions, the numerical results for the dynamics of the system are more clear than the analytic ones due to the complexity of the Lambert function. Thus, we numerically solve the evolutions of the inflaton field. The behaviors of the inflaton potential for the different values of $\beta$ are shown in the Fig.\[fig-1\]. $M = 6 \times 10^{12}$ GeV is used in this figure. The solid, dotted, and dashed lines depict the evolution of the inflaton potential when $\beta = 0, -0.02$, and 0.02, respectively. The requirement for the attractive gravity and ghost free graviton is given by $\partial f/ \partial R \, (\partial \bff/ \partial y) > 0$. In order to avoid a curvature singularity and Dolgov-Kaasaki instability, one also needs the condition $\partial^2 f/ \partial R^2 \, (\partial^2 \bff/ \partial y^2) > 0$ [@09091737]. These conditions provide the upper bound on the values of $\beta$ for the given value of $\alpha$ at each $R$. Thus, for the negative value of $\beta$, the model has the upper bound on the $\phi$ ([*i.e. UV incomplete*]{}). When $\beta = 0$, the model is identical to the Starobinsky model and show the plateau region to have the enough number of e-foldings. When the $\beta$ is positive, the potential becomes unstable with a runaway direction for the large values of $\phi$. In this case, the hilltop type model is realized when $\phi$ rolls towards the origin and $\beta$ value will be constrained in order to get the enough e-folding number.

![image](Fig1.pdf){width="1\linewidth"}

The detail evolution of the inflation model given in the action Eq.(\[Sx2lnx\]) is investigated for the different values of coupling constant $\beta$. These results are obtained from Eqs.(\[bphi\]), (\[V\]), and (\[bphiEq\]). We use $M = 6 \times 10^{12}$ GeV and $\alpha = 1$. In order to obtain the enough number of e-foldings, one needs to tune the initial values of the scalar field, $\bphi_{i}$. For an illustration, we show a typical evolutionary behavior of $\phi$ and $V(\phi)$ in Fig.\[fig-2\]. In this figure, $\beta = 0$ and N-efolings are fixed at 60. In this model, $\bphi_{i}$ is chosen to be 41.36 in order to obtain $N_{\text{e-folding}} = 60$ ([*i.e.*]{} $\epsilon(N=60) =1$). In this case, the scalar field slowly rolls down through the flat plateau region and reaches to the minimum to oscillate. The oscillatory behavior of the $\bphi$ is shown in the left panel of the Fig.\[fig-2\] and the slow-roll behaviors is represented on the right panel of the Fig.\[fig-2\].

  -- --
     
  -- --

In tables.\[tab-1\] and \[tab-2\], we summarize the values of the cosmological parameters obtained from the different values of $\beta$ when we fix the number of e-foldings. In table.\[tab-1\], we choose the e-folding number as 60 and vary the value of $\beta$ from -0.02 to 0.05. This range of $\beta$ provides the allowed regions for the cosmological parameters given in the reference [@150202114]. The smaller the values of $\beta$, the steeper the potential in the slow-roll plateau. Thus, one needs the larger initial field values for the smaller values of $\beta$ to obtain the same value of the e-folding number. The $\bphi_{i}$ varies from 75.2 to 16.4 for $-0.02 \leq \beta \leq 0.05$. One of slow-roll parameters, $\epsilon_{\text{V}}$ decreases as $\beta$ increases. However, as $\beta$ decreases, so does the other slow-roll parameter, $\eta_{\text{V}}$. The scalar spectral index, $n_{s}$ decreases as $\beta$ decreases. $n_{s}$ varies from 0.98 to 0.93 in the given range of $\beta$. The tensor-to-scalar ratio decreases from $7.4 \times 10^{-3}$ to $4.4 \times 10^{-4}$ when $\beta$ changes from -0.02 to 0.05. The magnitude of the running of the spectral index is maximum around Starobinsky model. The amplitude of the scalar primordial curvature perturbation increases as $\beta$ increases.

  -------- ------- ------- -------------- ----------------------- ------------------- ------- -------------
                                $r$        $\epsilon_{\text{V}}$   $\eta_{\text{V}}$              $A_s$
                            $(10^{-3}) $        $(10^{-4})$          $(10^{-4}) $       $$     $(10^{-9})$
   $ -20$   75.19   0.980       7.41               4.63                  0.76          -5.22      0.78
   $-10$    54.37   0.973       4.76               2.98                  1.53          -5.45      1.14
   $ -1$    42.42   0.968       3.25               2.03                  2.42          -5.51      1.61
     0      41.36   0.967       3.12               1.95                  2.54          -5.51      1.67
    $1$     40.33   0.966       2.99               1.87                  2.65          -5.51      1.73
    $10$    32.70   0.960       2.07               1.29                  3.79          -5.45      2.42
   $20 $    26.65   0.953       1.39               0.87                  5.31          -5.30      3.48
   $ 50$    16.43   0.932       0.44               0.28                  11.4          -4.49      10.0
  -------- ------- ------- -------------- ----------------------- ------------------- ------- -------------

  : These are the cosmological parameters predicted from the models when one chooses $N_{\text{e-foldings}} = 60$.[]{data-label="tab-1"}

The table.\[tab-2\] shows the viable cosmological parameters when the e-folding number is 50. The $\bphi_{i}$ varies from 57.5 to 15.3 when $\beta$ changes from -0.02 to 0.05. $\epsilon_{\text{V}}$ decreases from $5.8 \times 10^{-4}$ to $5.5 \times 10^{-5}$ for the given range of $\beta$. $\eta_{\text{V}}$ increases from $1.33 \times 10^{-4}$ to $1.33 \times 10^{-3}$ for the same interval of $\beta$. The scalar spectral index, $n_{s}$ decreases as $\beta$ decreases. $n_{s}$ varies from 0.97 to 0.93 in the given range of $\beta$. The tensor-to-scalar ratio decreases from $9.2 \times 10^{-3}$ to $8.9 \times 10^{-4}$ when $\beta$ changes from -0.02 to 0.05. The magnitude of the running of the spectral index is maximum around Starobinsky model. The amplitude of the scalar primordial curvature perturbation increases from $6.2 \times 10^{-10}$ to $5.0 \times 10^{-9}$ as $\beta$ increases.

  -------- ------- ------- -------------- ----------------------- ------------------- ------- -------------
                                $r$        $\epsilon_{\text{V}}$   $\eta_{\text{V}}$              $A_s$
                            $(10^{-3}) $        $(10^{-4})$          $(10^{-4}) $       $$     $(10^{-9})$
   $ -20$   57.45   0.973       9.23               5.77                  1.33          -7.67      0.62
   $-10$    43.76   0.967       6.35               3.97                  2.35          -7.88      0.85
    $-1$    35.41   0.961       4.60               2.88                  3.46          -7.95      1.13
     0      34.64   0.960       4.44               2.78                  3.60          -7.95      1.17
    $1$     33.90   0.960       4.29               2.68                  3.73          -7.95      1.20
    $10$    28.27   0.954       3.15               1.97                  5.08          -7.89      1.58
    $20$    23.63   0.947       2.27               1.42                  6.79          -7.73      2.13
    $50$    15.31   0.927       0.89               0.55                  13.3          -6.90      5.00
  -------- ------- ------- -------------- ----------------------- ------------------- ------- -------------

  : These are the cosmological parameters predicted from the models for $N_{\text{e-foldings}} = 50$.[]{data-label="tab-2"}

Thus, even though the logarithmic corrections to the Starobinsky model is natural if one accepts the fact that $f(R)$-gravity is originated from the quantum corrections of the matter in the curved space-time, the correction is strongly constrained in order to satisfy the observation. In semi-classical theory, one can obtain the $R^2$-term from the local action and $R^2 \ln R$-term from the non-local action. The coupling constant $\alpha$ can be changed by varying the renormalization scale $M$, while the coefficient $\beta$ is related to the trace anomaly. $\beta$-value can be determined by the numbers of quantum fields for the different spins. And usually $\alpha \gg \beta$ is satisfied [@Vilenkin85]. Thus, if one can obtain the constrain of $\beta$ accurately, then one might be able to obtain the number of fields during the early universe. Even though,$\beta$ should be less than 1% compared to $\alpha$, its effect on the cosmological parameters can be distinguishable if the accuracies of observations reach to percentage level.

Reheating
=========

The basic reheating process of the Universe at the end of the inflation is that the oscillating inflaton fields produce radiations via the tree-level decay of inflaton particles into relativistic species. Reheating is terminated when the rate of expansion of the Universe becomes smaller than the total decay rate of the inflaton field into new fields. Even though this perturbative reheating process is not the full story, it is good enough to estimate the reheating epoch. Thus, we show the oscillatory epoch and reheating period based on the given model. The reheating process and reheating temperature based on the Starobinsky model have been investigated in the literature [@09091737; @12041472; @12124466; @160604346]. We probe the reheating process for the $R^2 \ln R$-correction model and investigate any difference compared to those of the Starobinsky model.

Oscillation
-----------

At the end of the inflation, the scalar fields reach to the minimum of the potential and start to oscillate around it. In order to describe this period properly, one needs to describe this epoch semi-analytically. The field equation of the general $f(R)$-gravity theories is obtained from the action (\[Sfx\]) $$f'(R) \tensor{R}{_\mu_\nu} - \frac{1}{2} f(R) \tensor{g}{_\mu_\nu} + \left[ \tensor{g}{_\mu_\nu} - \nabla_{\mu} \nabla_{\nu} \right] f'(R) = \frac{1}{\Mp^2} \tensor{T}{_\mu_\nu} \, , \label{fEq}$$ where prime means the derivatives with respect to $R$. For the homogeneous and isotropic background, one can obtain two field equations from the above Eq. (\[fEq\]) $$\begin{aligned}
 \dot{R} &= - \frac{f'}{f''} H + \frac{1}{6H} \frac{f' R - f}{f''} + \frac{1}{f''H} \frac{\rho}{3 \Mp^2} \, , \label{dotR} \\
\ddot{R} &= - 3H \dot{R} - \frac{f'''}{f''} \dot{R}^2 + \frac{ R f' - 2f}{3 f''} + \frac{1}{f''} \frac{\rho - 3P}{3 \Mp^2} \, , \label{ddotR} \end{aligned}$$ where $\rho$ and $P$ are the energy density and the pressure of some species, respectively. If one adopts the specific form of $f(R)$ given by Eq.(\[SR2lnR\]), then one obtains $$\begin{aligned}
 & \dot{R} + \left(\frac{M^2 }{2\alpha+3\beta+2\beta \ln x} \right) H + \left(\frac{2\alpha+\beta+2\beta \ln x }{2\alpha+3\beta+2\beta \ln x} \right) H R \nonumber \\ 
&- \left(\frac{ \alpha+\beta+\beta \ln x }{2\alpha+3\beta+2\beta \ln x} \right) \frac{R^2}{6H} = \left(\frac{M^2}{2\alpha+3\beta+2\beta \ln x} \right)  \frac{\rho}{3 \bMp^2 H} \, , \label{dotR2} \\
& \ddot{R} + 3H \dot{R} + \left(\frac{2\beta}{2\alpha+3\beta+2\beta \ln x} \right)\frac{\dot{R}^2}{R} + \frac{M^2}{3} \left(\frac{1-\beta x}{2\alpha+3\beta+2\beta \ln x} \right) R \nonumber \\ 
&= \left(\frac{M^2}{2\alpha+3\beta+2\beta \ln x} \right)\frac{\rho - 3P}{3 \bMp^2} \, . \label{ddotR2} \end{aligned}$$ One can use the relation between the Ricci scalar and the Hubble parameter, $R = 6 \left( \dot{H} + 2 H^2 \right)$ to investigate the evolution of $H$. Then, the equation (\[dotR\]) is rewritten as $$\begin{aligned}
 & \ddot{H} - \left( \frac{\alpha+\beta+\beta \ln x}{2\alpha+3\beta+2\beta \ln x} \right) \frac{\dot{H}^2}{H} + 3 H \dot{H} + \frac{M^2}{6 \left( 2\alpha+3\beta+2\beta \ln x \right)} H \ \nonumber \\
& - \frac{2\beta}{2\alpha+3\beta+2\beta \ln x } H^3 =  \left(\frac{M^2}{2\alpha+3\beta+2\beta \ln x} \right) \frac{\rho}{3 \bMp^2 H} \, . \label{dotR3} \end{aligned}$$ Now one can scrutinize the behaviors of $R$ and $H$ from Eqs.(\[dotR2\])-(\[dotR3\]). First, one can puts $\dot{R} = 0$ and $\rho = 0$ in the slow-roll region. Then one obtains $$\begin{aligned}
 & R_{i} \simeq 3 b H_{i}^2 \left( 1 + \sqrt{1 + \frac{2c}{3b^2 H_{i}^2}} \, \right) \label{Ri} \, , \\
& \text{where} \,\, b \equiv \frac{2\alpha+\beta+2\beta \ln x}{\alpha+\beta+\beta\ln x} \,\, , \,\, c \equiv \frac{M^2}{\alpha+\beta+\beta\ln x} \, . \label{bc} 
 \end{aligned}$$ As one expects $R_{i} \sim 12 H_{i}^2$ in this region, because $b \sim 2$ when $\alpha \gg \beta$. Even though $H$ is almost constant, it decreases almost linearly in time with a small slope. This can be obtained from Eq.(\[dotR3\]) $$\begin{aligned}
 \dot{H} &\simeq - \frac{M^2}{18 \left(2\alpha+3\beta+2\beta \ln x \right)} + \frac{2\beta}{3 \left(2\alpha+3\beta+2\beta \ln x \right) } H^2 \nonumber \\ 
& \simeq - \frac{M^2}{18 \left(2\alpha+3\beta+2\beta \ln x \right)} \label{SRdotH} \, , \end{aligned}$$ where we use the approximation $\beta \ll 1$ in the second equality. In order to obtain the finite period of inflation, $\left(2\alpha+3\beta+2\beta \ln x \right)$ should be positive. One also obtains the Hubble parameter from the Eq.(\[SRdotH\]) $$H(t) \simeq H_{i} - \frac{M^2}{18 \left(2\alpha+3\beta+2\beta \ln x \right)} \left( t - t_{i} \right) \, . \label{Hit}$$ From this, one can estimate the expansion of the scale factor of the Universe during this region $$a(t_{\text{end}}) \simeq a_{i} \exp \left[ \frac{9 \left( 2\alpha + 3\beta + 2\beta \ln x_{i} \right) H_{i}^2}{M^2} \right] \label{SRa} \, ,$$ where we use the fact that $\ln x$ is almost constant in this region. Thus, one can constrain the $M^2$ to satisfy the e-foling number for the given values of $\alpha$ and $\beta$. We already do this numerically and do not need to probe any detail in this section. Second, we consider the oscillation period. During the slow-roll period, $\Bigl|\frac{1}{2} \frac{\dot{H}^2}{H} \Bigr| \ll \Bigl|3 \dot{H} H \Bigr|$ is satisfied. However, the magnitudes of these two terms become comparable as $H$ decreases and becomes small. And oscillatory phase is followed with $\Bigl|3 \dot{H} H \Bigr| \simeq 0$. $$\ddot{H} - \left( \frac{\alpha+\beta+\beta \ln x}{2\alpha+3\beta+2\beta \ln x} \right) \frac{\dot{H}^2}{H} + \frac{M^2}{6 \left( 2\alpha+3\beta+2\beta \ln x \right)} H + \frac{2\beta}{2\alpha+3\beta+2\beta \ln x } H^3 =  0 \, . \label{dotR3osc}$$ In order to obtain the oscillation period, $H^3$-term should be negligible ([*i.e.*]{} $\beta \simeq 0$). With this condition, the solution for Eq.(\[dotR3osc\]) is given by $$\begin{aligned}
 & H(t) = c_{1} \times \cos^{\frac{M^2}{6 \omega^2}} \left[ \omega (t - t_{\text{osc}}) \right] \, , \nonumber \\ 
 & \text{where} \,\, \omega = \frac{M \sqrt{\alpha+2\beta+\beta \ln x}}{\sqrt{6} (2\alpha+3\beta+2\beta \ln x)} \label{Hosc} \,\,  \text{and} \,\, t_{\text{osc}} =  6 c_{2} \left( 2\alpha+3\beta+2\beta \ln x \right) \, . \end{aligned}$$ $c_1$ and $c_2$ are integral constants. $t_{\text{osc}}$ is the time when the oscillation period begins with $\Bigl|\frac{1}{2} \frac{\dot{H}^2}{H} \Bigr| = \Bigl|3 \dot{H} H \Bigr|$. One can estimate the order of the frequency, $\omega$ at this epoch. ${\cal O} (\omega) =  {\cal O} (M/\hbar) \sim 10^{37}$ Hz. We point out the oscillatory behavior before derive the approximate solution. The oscillation is powered by $\frac{M^2}{6\omega^2}$ and thus it should be an integer. This means $\beta \simeq 0$. Also, the oscillation should be damped and thus, one needs to replace $c_{1}$ with some function $g(t)$ in Eq.(\[Hosc\]). Now, one needs to specify the $g(t)$ in the $\alpha \gg \beta$ limit. For this we compare the $H(t)$s in both before and during the oscillation region. $$\label{Htiosc}
        H(t) = \begin{dcases} H_{i} - \dfrac{2\omega^2}{3} \left(t-t_{i} \right) = \dfrac{\omega}{3} \left( 1 - 2 \omega \left( t - \tosc \right) \right)  &, \,\,  t_{i} < t \leq  \tosc \\
                                    \dfrac{\cos^{2} \left[\omega \left(t - \tosc \right) \right] }{\left(\frac{15}{2\omega }+6(t-\tosc)-\frac{9}{2\omega } \cos [2\omega (t-\tosc)]\right)}  & , \,\, \tosc \leq t \, ,
                    \end{dcases}$$ where $\tosc = t_{i} + (3/2\omega^2) H_{i} - (1/2\omega)$. We show the details to obtain these solutions in the appendix. Thus, after the long inflationary plateau, $H(t)$ decreases linearly and reaches to the oscillation period. The approximate solutions for the scale factor $a$ are obtained from $H$s in Eqs.(\[Htiosc\]) $$\label{attiosc} a(t) = \begin{dcases} a_{i} e^{H_{i}\left(t - t_{i}\right) - (\omega^2/3) \left( t -t_{i} \right)^2} &, \,\, t_{i} < t \leq t_{\text{osc}}  \\
             a_{\osc} \left( 1 + \dfrac{\omega \left(t - \tosc \right)}{4} \right)^{2/3}  &, \,\, \tosc \leq t \, , 
             \end{dcases}$$ where one can obtain the $a(t)$ in the oscillatory phase by integrating the $H$ averaged over a few cycles. During the corresponding periods, the Ricci scalar evolves as $$\label{Rtiosc} 
    R(t) = \begin{dcases} -4 \omega (\omega -1) - 8\omega ^2 \left( t-t_{i} \right) &, \,\, t_{i} < t \leq t_{\text{osc}} \, , \\ 
                                              \frac{-4 \omega  \cos \left[\left(t-t_{\text{osc}}\right) \omega \right]}{\left(5+4 \omega \left( t - t_{\text{osc}} \right) -3 \cos \left[2 \left(t-t_{\text{osc}}\right) \omega \right]\right){}^2} \times & \\ 
                                              \Biggl( \Bigl(-7+\left(4-8 t+8 t_{\text{osc}}\right) \omega \Bigr) \cos \left[\left(t-t_{\text{osc}}\right) \omega \right] & \\ 
                                              + 3 \cos \left[3 \left(t-t_{\text{osc}}\right) \omega \right]+8 \omega  \Bigl( 2+\left(t -t_{\text{osc}}\right) \omega \Bigr) \sin [(t-\tosc) \omega ] \Biggr) &, \,\, \tosc \leq t  \, . \end{dcases}$$

  -- --
     
  -- --

We show the evolution of $H$ and $R$ in Fig.\[fig-2\]. One should note that the shape and size of the oscillatory phase depend only on the $M$ and the $\alpha$ because one should use $\alpha \gg \beta$ in order to obtain the oscillatory behavior on the Hubble parameter as shown in the above.

Reheating
---------

The scalaron reaches to its minimum potential at $\phi = 0$ ([*i.e.*]{} $R=0$) and starts to oscillate. These oscillations excite the fields and reheat the Universe. To estimate the reheating, one can add the simple case of a scalar field $\chi$ to the original action in Eq.(\[Sfx\]) $$\begin{aligned}
 S &\equiv \int d^{4} x \sqrt{-g} \left( \frac{\bMp^2}{2} f \left(R \right) + {\cal L}_{\chi} + {\cal L}_{r} \right) \label{SfRchi}  \\
&= \frac{1}{2} \int d^{4} x \sqrt{-g} \left( \bMp^2 f \left(R \right) -  \tensor{g}{^\mu^\nu} \partial_{\mu} \chi \partial_{\nu} \chi - m_{\chi}^2 \chi^2 \right) + S_{r} \, , \nonumber \end{aligned}$$ where $S_{r}$ is the action of the radiation which is produced by the decay of $\chi$ field after the inflation. In this scenario, we do not consider the coupling between $R$ and $\chi$ ([*i.e.*]{} no decay of the scalaron). One can use the same field equation given in Eq.(\[fEq\]) by specifying the energy momentum tensor as $$\tensor{T}{^\mu^\nu} \equiv {T_{\chi}}\tensor{\vphantom{T}}{^\mu^\nu} + {T_{r}}\tensor{\vphantom{T}}{^\mu^\nu}  = \partial^{\mu} \chi \partial^{\nu} \chi - \tensor{g}{^\mu^\nu} {\cal L}_{\chi} + \left( \rho_{r} + P_{r} \right) U^{\mu} U^{\nu} + P_{r} \tensor{g}{^\mu^\nu} \, , \label{Tmunu}$$ where $U^{\mu}$ is the four-velocity and $\rho_{r}$ and $P_{r}$ are the energy density and pressure of the radiation, respectively. In the FRW metric, one obtains the field equations for the $\chi$ and radiation, $$\begin{aligned}
 \ddot{\chi} + 3H \dot{\chi} + m_{\chi}^2 \chi &= -\Gamma_{\chi} \dot{\chi} \,\, , \label{chifeq} \\
\dot{\rho}_{r} + 4 H \rho_{r} &= \Gamma_{\chi} \dot{\chi}^{2} \,\, . \label{dotrho} \end{aligned}$$ Because one is interested in the oscillatory region in the background produced by the scalaron, the background evolutions are given by Eqs.(\[Htiosc\])-(\[Rtiosc\]). It means that one can ignore any backreaction on the background evolutions from $\chi$ and $\rho_{r}$ during this period. Thus, one can solve the field equations (\[chifeq\]) and (\[dotrho\]) both analytically and numerically. For this purpose, one can rewrite these equations as a function of e-folding numbers, $N$ $$\begin{aligned}
  \frac{d^2 \chi}{dN^2} + \left( \frac{3}{2} + \frac{\Gamma_{\chi}}{H\left(N\right)} \right) \frac{d \chi}{d N} + \frac{m_{\chi}^2}{H(N)^2} \chi &= 0 \,\, , \label{chifeqN} \\
\frac{\rho_{r}}{dN} + 4 \rho_{r} - \Gamma_{\chi} H(N) \left( \frac{d \chi}{d N}  \right)^{2}&= 0 \,\, , \label{Nrhor} \\
 \text{where} \,\, N \equiv \ln \left[ \frac{a}{a_{\osc}} \right] = \frac{2}{3} \ln \left[ 1 + \frac{\omega \left(t - \tosc \right)}{4} \right] \,\, , \,\,  H(N) &= \frac{\omega}{6} e^{-\frac{3}{2} N} \,\, .  \label{HN} \end{aligned}$$ As can be seen in Eq.(\[Nrhor\]), $\rho_{r}$ is created from the decaying of the $\chi$ field. At the end of the inflation, $\Gamma$ is typically much smaller than the Hubble parameter. Thus, at the beginning of the phase of the inflationary oscillations, the energy loss into particles is negligible compared to the energy loss due to the expansion of the space time. $\chi$ particle production becomes effective only when the Hubble expansion rate decreases to a value comparable to $\Gamma$. Thus, $\rho_{r}$ becomes the constant value of $\Gamma_{\chi} H(N_{\osc}) (d\chi/dN |_{N=N_{\osc}})^2$ at the not too long after the Universe has come into the oscillation phase. One can define the end of the reheating when the 90% of the energy density of the $\chi$ field is converted to that of the radiation ([*i.e.*]{} $\rho_{r} = 9 \rho_{\chi}$) where $\rho_{\chi}$ is given by $$\rho_{\chi} = \frac{1}{2} \left( H^2 \left( \frac{d \chi}{d N} \right)^2 + m_{\chi}^2 \chi^2 \right) = \frac{H^2 M^2}{2} \left( \left( \frac{d \bar{\chi}}{d N} \right)^2 + \frac{m_{\chi}^2}{H^2} \bar{\chi}^2 \right) \, , \label{rhochi}$$ where we define the dimensionless field variable $\bar{\chi} = \chi /M$. The energy density of the radiation at the temperature $T$ is given by $$\rho_{r}(N) = \frac{g_{\ast}(N) \pi^2}{30} T(N)^{4} \, , \label{rhor}$$ where $g_{\ast}(N)$ is the number of relativistic species degrees of freedom and it is around ${\cal O}(100)$ at an early epoch. From the definition of the reheating epoch, one can obtain the reheating temperature $$T_{\rh} = \left( \frac{9}{50 \pi^2} \left( \frac{100}{g_{\ast}} \right) \rho_{\chi}(N_{\rh}) \right)^{1/4} \, . \label{Trh}$$ The reheating temperature is determined by numerically solving Eqs.(\[chifeq\]) and(\[dotrho\]). This depends on model parameters, $\alpha$, $\beta$, $M$, $\Gamma_{\chi}$, and $m_{\chi}$. We demonstrate the evolutions of $\rho_{\chi}$ and $\rho_{r}$ for the different models in Figs.\[fig-4\] and \[fig-5\].

  -- --
     
  -- --

In the figure.\[fig-4\], we depict the evolutions of $\rho_{r}$ and $\rho_{\chi}$ when $m_{\chi} = M$ and $\beta = 10^{-2} \alpha$. In the left panel of Fig.\[fig-4\], we choose the small decay rate, $\Gamma = 10^{-3} m_{\chi}$. In this case, it takes longer time ([*i.e.*]{} larger $N_{\text{rh}}$) to reheat the Universe compared to the larger decay rate models. Thus, one obtains the lower reheating, $T_{\text{rh}} = 2.3 \times 10^{10}$ GeV. In this model, the $\rho_{r}$ reaches to its maximum values at $N_{\text{max}} = 0.37$ after the beginning of the oscillatory epoch and terminates the reheating at $N_{\text{rh}} = 3.77$. In the right panel of the figure.\[fig-4\], the decay rate is now $3 m_{\chi}$. Thus, reheating process becomes faster and after $\rho_{r}$ reaches its maximum at $N_{\text{max}} = 0.20$, it becomes thermalized at $N_{\text{th}} =  0.31$. In this model, the reheating temperature is higher than the first model and becomes $T_{\text{rh}} = 9.1 \times 10^{10}$ GeV.

  -- --
     
  -- --

In the figure.\[fig-5\], the evolutions of $\rho_{r}$ and $\rho_{\chi}$ are depicted for different values of $m_{\chi}$ with $\Gamma = 3 m_{\chi}$ and $\beta = 10^{-2} \alpha$. In the left panel of Fig.\[fig-5\], the large mass scalar field, $m_{\chi} = 1.5 M$ is chosen. In this case, the decay process is more efficient than the smaller mass model and it takes less time ([*i.e.*]{} $N_{\text{rh}}$) to reach both the maximum $\rho_{r}$ and the reheating of the Universe compared to the smaller mass of the scalar field models. Thus, one obtains the higher reheating, $T_{\text{rh}} = 1.2 \times 10^{11}$ GeV. In this model, the $\rho_{r}$ reaches to its maximum values at $N_{\text{max}} = 0.15$ after the beginning of the oscillatory epoch and terminates the reheating at $N_{\text{rh}} = 0.21$. In the right panel of the figure.\[fig-5\], the mass of the scalar field, $m_{\chi} = 0.5 M$. Thus, reheating process becomes slower and after $\rho_{r}$ reaches its maximum at $N_{\text{max}} = 0.32$, it becomes thermalized at $N_{\text{th}} =  0.64$. In this model, the reheating temperature is lower than the model of the left panel, $T_{\text{rh}} = 5.9 \times 10^{10}$ GeV.

  ---------------- -------------- ------------ --------------------------------- -------- ------------------------------------------- -------- ----------------------
                                                $\rho _r\left(N_{\max }\right)$            $\rho _{\chi }\left(N_{\text{rh}}\right)$              $T_{\text{rh}}$
                                                 $\left[(\text{GeV})^4\right]$                   $\left[(\text{GeV})^4\right]$                     $[\text{GeV}]$
      \[1ex\]                      $-10^{-2}$        $1.33\times 10^{44}$         $0.37$             $1.49\times 10^{43}$              $3.77$   $2.28\times 10^{10}$
      \[1ex\]                      $10^{-2}$         $1.30\times 10^{44}$         $0.37$             $1.48\times 10^{43}$              $3.77$   $2.28\times 10^{10}$
    \[0.5ex\] M     $3m_{\chi }$   $10^{-2}$         $3.33\times 10^{46}$          0.20              $3.70\times 10^{45}$               0.31    $9.06\times 10^{10}$
   \[0.5ex\] 1.5M   $3m_{\chi }$   $10^{-2}$         $8.68\times 10^{46}$          0.15              $9.64\times 10^{45}$               0.21    $1.15\times 10^{11}$
   \[0.5ex\] 0.5M   $3m_{\chi }$   $10^{-2}$         $5.98\times 10^{45}$          0.32              $6.64\times 10^{44}$               0.64    $5.90\times 10^{10}$
     \[0.5ex\]                                                                                                                                 
  ---------------- -------------- ------------ --------------------------------- -------- ------------------------------------------- -------- ----------------------

  : The reheating temperature for the different models when $M = 6 \times 10^{12}$ GeV. We assume that the reheating phase is obtained when $\rho_{\chi}$ reaches to 90 % of the maximum $\rho_{r}$.[]{data-label="tab-3"}

We summarize the results in the table.\[tab-3\]. As shown in the figures.\[fig-4\] and \[fig-5\], the larger the decay rate, the earlier the reheating epoch. Thus, for the same values of $m_{\chi}$ and $\beta$, one obtains the higher reheating temperature for the larger decay rate. Also, if $m_{\chi}$ increases, so does the $T_{\text{rh}}$ due to the increasing of efficiency of the decay of the scalar field. We are also interested in the effect of the $\beta$ on the reheating temperature in order to investigate the possible constraint on its values from the reheating temperature. Even though, there exists small deviations on the reheating temperatures for the different values of $\beta$, the difference is less than 1 %.

Conclusions
===========

We have investigated the general modification of the Starobinsky inflation model by including the logarithmic correction in addition to the quadratic in the Ricci scalar. We show that the values of coefficients $\alpha$ and $\beta$ are strongly constrained from the observation. The maximum value of the coefficient $\beta$ is about 1 % of that of the coefficient $\alpha$ in order to be satisfied with the CMB results. $\beta \neq 0$ means the deviation of the model from the Starobinsky’s one. When we vary the value of $\beta$ from 0 to |0.02| ([*i.e.*]{} 2 % change in $\beta$), the scalar spectral index is changed by 1% only. We also investigate the reheating process in order to obtain any deviation from the Starobinsky model. The changes in the reheating temperatures are less than 1 % when we compare all the viable $\beta$ values. Even though, the effects of the deviation from the Starobinsky model on the known observational quantities are quite small, one might still be able to confirm the models with upcoming more accurate observations and (or) from other cosmological observables.

In this appendix, we show the detail derivation of the damped oscillation behavior given in the Eq.(\[Hosc\]).

Oscillation
===========

If one uses the Eq.(\[bphi\]), then Eq.(\[bphiEq\]) is given by $$\label{ffieldEq}$$

One can obtain the approximate damped oscillation of the Hubble parameter by replacing $c_{1}$ by function $g(t)$ in the equation (\[Hosc\]) $$H_{\text{osc}} (t) \simeq g(t) \times \cos^{2} \left[\omega \left(t - \tosc \right) \right] \,\, , \,\, \tosc \leq t \, , \label{Hoscdampapp}$$ where $\omega^2 = M^2/(24 \alpha)$. In order to obtain the proper damped oscillation, the form of $g(t)$ should be given by $$g(t) = \left(  c_{1(\osc)} + c_{2(\osc)} \left( t - \tosc \right) + c_{3(\osc)} \cos \left[ 2 \omega \left( t - \tosc \right) \right] \right)^{-1}  \, , \label{gtapp}$$ where $c_{i(\osc)}$ are constant to be determined from the boundary values of $H(\tosc)$ before and during the oscillation period. The boundary conditions are given by $$\begin{aligned}
 H(\tosc) &= H_{i} - \frac{2}{3} \omega^2 \left( \tosc - t_{i} \right) = H_{\osc} \left(\tosc \right) = \frac{1}{c_{1(\osc)} + c_{3(\osc)}} \, , \label{BC1} \\ 
\dot{H}(\tosc) &= -\frac{2}{3} \omega^2 = \dot{H}_{\osc} \left(\tosc \right) = - \frac{c_{2(\osc)}}{\left( c_{1(\osc)} + c_{3(\osc)}\right)^2} \, , \label{BC2} \\\ddot{H}(\tosc) &= 0 = \ddot{H}_{\osc} \left(\tosc \right) = \frac{2\left( c_{2(\osc)}^{2} +\left( c_{3(\osc)}^2 - c_{1(\osc)}^2 \right) \omega^2 \right)}{\left( c_{1(\osc)}  + c_{3(\osc)} \right)^3} \, , \label{BC3}  \end{aligned}$$ From these equations.(\[BC1\])-(\[BC3\]), one obtains $$c_{1(\osc)}  = \frac{15}{2\omega} \, , \, c_{2(\osc)}  = 6 \, , \, c_{3(\osc)} = - \frac{9}{2\omega} \, . \label{c1c2c3}$$ Thus, one obtains the analytic solution of $H(t)$ during the oscillation period as $$H_{\osc} (t) = \left(\frac{15}{2\omega }+6(t-\text{tosc})-\frac{9}{2\omega } \cos [2\omega (t-\tosc)]\right)^{-1} \cos^{2} \left[\omega \left(t - \tosc \right) \right] \, . \label{Hosc3}$$

energy density of the radiation $$\begin{aligned}
 \rho_{r}(N_{\osc}) &\simeq \Gamma_{\chi} H(N_{\osc}) M^2 \chi'\left(N_{\osc}\right)^2 \, , \label{rhoNosc} \\
\rho_{\chi}(N_{\osc}) &\simeq \frac{H(N_{\osc})^2 M^2}{2} \left( \bar{\chi}'(N_{\osc})^{2} + \frac{m_{\chi}^2}{H(N_{\osc})^2} \bar{\chi}(N_{\osc})^2 \right) \, ,\label{rhochiNosc} \\
\rho_{\chi}(N_{\rh}) &= \frac{H(N_{\rh})^2 M^2}{2} \left( \bar{\chi}'(N_{\rh})^{2} + \frac{m_{\chi}^2}{H(N_{\rh})^2} \bar{\chi}(N_{\rh})^2 \right) \, , \label{rhochiNrh} \\
\rho_{r}(N_{\rh}) &= \frac{\pi^2}{2} g_{\ast} T_{\rh}^4 \equiv 9 \rho_{\chi}(N_{\rh}) \, . \label{Trh} \end{aligned}$$ Thus, one can obtain reheating temperature $T_{\rh}$ as $$T_{\rh} = \left( \frac{9}{50 \pi^2} \left( \frac{100}{g_{\ast}} \right) \rho_{\chi}(N_{\rh}) \right)^{1/4} \, . \label{Trh2}$$

Acknowledgments
===============

SL is supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and Future Planning (Grant No. NRF-2017R1A2B4011168).

[99]{}

N. D. Birrell and P. C. W. Davies. Quantum Fields in Curved Space. Cambridge University Press, April 1984.

R. M. Wald. Quantum field theory in curved spacetime. September 1995, \[gr-qc/9509057\].

L. H. Ford. Quantum field theory in curved spacetime. Jul 1997, \[gr-qc/9707062\].

C. Baer and K. Fredenhagen, editors. Quantum Field Theory on Curved Spacetimes. Springer, 2009.

L. Parker, D. Toms, Quantum Field Theory in Curved Spacetime, Cambridge University Press, 2009.

H. Nariai, Prog. of Theore. Phys [**49**]{}, 165 (1973). M. B. Mijic, M. S. Morris and W.-M. Suen, Phys. Rev. D [**34**]{}, 2934 (1986). A. A. Starobinsky, Phys. Lett. B [**91**]{}, 99 (1980). B. S. DeWitt, Phys. Rev [**162**]{}, 1239 (1967) S. Nojiri and S. D. Odintsov, Gen. Rel. Grav. [**36**]{}, 1765 (2004) \[arXiv:hep-th/0308176\].

X.-H. Meng and P. Wang, Phys. Lett. B [**584**]{}, 1 (2004) \[arXiv:hep-th/0309062\]. I. Navarro and K. V. Acoleyen, JCAP [**0609**]{}, 006 (2006) \[arXiv:gr-qc/0512109\]. S. A. Appleby and R. A. Battye, Phys. Lett. B [**654**]{}, 7 (2007) \[arXiv:0705.3199\]. R. Saffari and S. Rahvar, Phys. Rev. D [**77**]{}, 104028 (2008) \[arXiv:0708.1482\]. S. A.‘Appleby and R. A. Battye, JCAP [**0805**]{}, 019 (2008) \[arXiv:0803.1081\]. Z. Girones, A. Marchetti, O. Mena, C. Pena-Garay and N. Rius, JCAP [**1011**]{}, 004 (2010) \[arXiv:0912.5474\]. J.-Q. Guo and A. V. Frolov, Phys. Rev. D [**88**]{}, 124036 (2013) \[arXiv:1305.7290\]. H. Alavirad and J. M. Weller, Phys. Rev. D [**88**]{}, 124034 (2013) \[arXiv:1307.7977\]. T. Gurovich and A. A. Starobinsky, Sov. Phys. JETP [**50**]{}, 844 (1979). G. M. Shore, Ann. Phys. (N.Y.) [**128**]{}, 376 (1980). A. Vilenkin, Phys. Rev. D [**32**]{}, 2511 (1985). M. B. Baibosunov, V. Ts. Gurovich, and U. M. Imanaliev, Sov. Phys. JETP [**71**]{}, 636 (1990). I. Ben-Dayan, S. Jing, M. Torabian, A. Westphal, and L. Zarate, JCAP [**09**]{}, 005 (2014) \[arXiv:1404.7349\]. M. Rinaldi, G. Cognola, L. Vanzo and S. Zerbini, JCAP [**08**]{}, 015 (2014) \[arXiv:1406.1096\]. M. Rinaldi, G. Cognola, L. Vanzo and S. Zerbini, Phys. Rev. D [**91**]{}, 123527 (2015) \[arXiv:1410.0631\]. B. J. Broy, F. G. Pedro and A. Westphal, JCAP [**03**]{}, 029 (2015) \[arXiv:1411.6010\]. M. Amin, S. Khalil, and M. Salah, JCAP [**08**]{}, 043 (2016) \[arXiv:1512.09324\]. J. Sadeghi and H. Farahani, Phys. Lett. B [**751**]{}, 89 (2015) \[arXiv:1605.08016\]. P. A. R. Ade [*et al*]{}. \[Planck Collaboration\], Astron. Astrophys. [**594**]{}, 20 (2016) \[arXiv:1502.02114\]. B. N. Breizman, V. Ts. Gurovich, and V. P. Sokolov, Sov. Phys. JETP [**32**]{}, 155 (1971). T. d. P. Netto, A. M. Pelinson, I. L. Shapiro, and A. A. Starobinsky, Eur. Phys. J. C [**76**]{}, 544 (2016) \[arXiv:1509.08882\]. X. Calmet and I. Kuntz, Eur. Phys. J. C [**76**]{}, 289 (2016) \[arXiv:1605.02236\]. Z. Lalak and L. Nakonieczny, Phys. Dark. Univ. [**15**]{}, 125 (2017) \[arXiv:1609.06887\]. T. Alanne, F. Sannino, T. Tenkanen, and K. Tuominen, Phys. Rev. D [**95**]{}, 035004 (2017) \[arXiv:1611.04932\]. Y.-C. Wang and T. Wang, Phys. Rev. D [**96**]{}, 123506 (2017) \[arXiv:1701.06636\]. Y. Ema, Phys. Lett. B [**770**]{}, 403 (2017) \[arXiv:1701.07665\]. T. Mori, K. Kohri, and J. White, JCAP [**10**]{}, 044 2017 \[arXiv:1705.05638\]. A. Karam, T. Pappas, and K. Tamvakis, Phys. Rev. D [**96**]{}, 064036 (2017) \[arXiv:1707.00984\]. R. Allahverdi, R. Brandenberger, F.-Y. Cyr-Racine and A. Mazumdar, Annu. Rev. Nucl. Part. Sci. [**60**]{}, 27 (2010) \[arXiv:1001.2600\]. S. Appleby, R. Battye, and A. Starobinsky, JCAP [**06**]{}, 005 (2010) \[arXiv:0909.1737\]. H. Motohashi and A. Nishizawa, Phys. Rev. D [**86**]{}, 083514 (2012) \[arXiv:1204.1472\]. D. Gorbunov and A. Tokareva, JCAP [**12**]{}, 021 (2013) \[arXiv:1212.4466\]. C. van de Bruck, P. Dunsby, L.E. Paduraru, , Int. J. Mod. Phys. D [**26**]{}, 1750152 (2017) \[arXiv:1606.04346\].
