Production of gravitational waves during preheating in the Starobinsky inflationary model

The production of GWs during preheating in the Starobinsky model with a nonminimally coupled auxiliary scalar field is studied through the lattice simulation in this paper. We find that the GW spectrum $\Omega_{\rm gw}$ grows fast with the increase of the absolute value of coupling parameter $\xi$. This is because the resonant bands become broad with the increase of $|\xi|$. When $\xi<0$, $\Omega_{\rm gw}$ begins to grow once the inflation ends and grows faster than the case of $\xi>0$. $\Omega_{\rm gw}$ reaches the maximum at $\xi=-20$ ($\xi=42$ for the case $\xi>0$) and then decreases with slight oscillation. Furthermore we find that the GWs produced in the era of preheating satisfy the limits from the Planck and next-generation CMB experiments.


Introduction
Inflation, a phase of quasi-exponential expansion in the early universe, is an elegant idea proposed to resolve most of theoretical problems in the big bang standard cosmology [1][2][3][4]. Furthermore, the scalar quantum fluctuations during inflation provide the seed for the formation of large scale structure [5][6][7][8][9]. Inflation usually is assumed to be driven by a slow-roll single scalar field, which predicts a scale-invariant spectrum of scalar fluctuations on the super-horizon scale. This prediction is consistent with the Cosmic Microwave Background (CMB) observations [10,11]. The latest CMB limits the spectral index of the power spectrum to be n s = 0.9649 ± 0.0042 at the 68% confidence level (CL) [12]. Except for the scalar perturbations, there also exists the tensor fluctuations during inflation, which is also scale-invariant on large scales. It is worthy to note that tensor perturbations in the scope of general relativity were first and quantitatively correctly calculated in [13], and the quantitatively correct expressions for both scalar and tensor perturbations were first presented a e-mail: pxwu@hunnu.edu.cn (corresponding author) b e-mail: hwyu@hunnu.edu.cn in [14] for the Starobinsky model based on modified f (R) gravity [1]. Since the tensor perturbations lead to the B-mode polarization for the CMB photon, its amplitude can be constrained by the CMB observations. The allowed region of the ratio of tensor to scalar fluctuations r is r 0.002 < 0.10 at the 95% CL [12]. Thus, the combination of n s and r can discriminate a host of inflationary models. In this regard, it is well-known that the Starobinsky inflationary model is in excellent agreement with the latest CMB observations [12].
After inflation, the inflaton usually will undergo a periodical oscillation around the minimum of its potential. During this oscillation the inflation field will decay into some light particles to thermalize the universe, which is called reheating [15]. In the first stage of reheating, i.e., preheating, the periodical oscillation of the inflaton may lead to an explosive production of the inflaton quanta or other light particles coupled to the inflation field through the parametric resonance [16][17][18][19][20][21][22]. When the mode momenta of the inflaton quanta or the light fields are in the resonance bands, these Fourier modes will grow exponentially. Since only a part of modes have the resonant momenta, the matter distribution has large and time-dependent inhomogeneities in the position space, which results in that the matter possesses substantial quadruple moments, and becomes an effective source of significant gravitational waves (GWs) [23]. The production of GWs during preheating has been studied widely in [24][25][26][27][28][29][30]. Different from vacuum fluctuations of tensor perturbations during inflation, the amplitude of the GW spectrum generated during preheating is independent of the energy scale of inflation which only determines the present peak frequency [31,32]. For low energy scale inflationary models, the peak frequency of GWs produced after inflation may well occur in the range which in principle can be detected by future direct detection experiments (like LIGO/VIRGO) [29,[33][34][35]. This opens a unique observational window for us to test inflation and the subsequent dynamics of the very early universe.
It has been found that in the Starobinsky model the creation of matter can be realized through the decay of scalarons [1,36,37]. However, the inflaton field in the Starobinsky mode is absent of self-resonance, and in order to realize the preheating process an auxiliary scalar field χ coupled non-minimally with the gravity is required [38]. Recently, the rescattering between the χ particles and the inflaton condensate during preheating has been studied in [39] by using the lattice simulation. The result shows that the rescattering is an efficient mechanism promoting the growth of the χ field variance, and knocks copious inflaton particles out of the inflaton condensate. Thus, both the auxiliary field and the inflaton field become the effective gravitational wave sources. However, Ref.
[39] does not investigate the production of GWs during preheating, which motivates us to finish the present work.
In this paper, we investigate the production of GWs during preheating in the Starobinsky inflation model. The energy density and spectrum of GWs will be analyzed detailedly by using the lattice simulation. Throughout this paper, we adopt the metric signature (−, +, +, +). Latin indices run from 0 to 3, Greek letters do from 1 to 3, and the Einstein convention is assumed for repeated indices.

Basics equations
For the Starobinsky model, when an auxiliary scalar field χ coupled non-minimally with gravity is considered, the Lagrangian has the form [38] where g is the determinant of metric tensor g μν , R is the Ricci scalar, κ −1 = M p with M p being the reduced Planck mass, μ is a constant, which is fixed at 1.3 × 10 −5 M p by the magnitude of the primordial density perturbations [14,40], and ξ is an arbitrary coupling parameter. In our following analyses, besides the weak coupling, the case of strong coupling |ξ | 1 is also considered. The part of the reason is that with the increase of |ξ | (ξ < 0) the value of r decreases and enters the 68% region allowed by the Planck observations in the nonminimally coupled inflation model with a self-coupling potential [41]. Usually, it is convenient to discuss the cosmic dynamics in the Einstein frame. By taking a conformal transformation:ĝ μν = 2 g μν with 2 being 2 = 1 − ξκ 2 χ 2 + R/(3μ 2 ), and introducing a scalar field φ ≡ √ 3/2κ −1 ln 2 as the inflaton, one can obtain the Lagrangian in the Einstein frame where In a flat Friedmann-Robertson-Walker (FRW) background, we obtain the following field equations from the Lagrangian and Here a dot denotes a derivative with respective to the cosmic time t, a is the scale factor, and H is the Hubble parameter. During the inflation which is driven by the inflaton field φ, the χ field has a negligible effect and thus can be neglected. After inflation, the field φ behaves as an inflaton with the quadratic potential and oscillates coherently around φ = 0. This oscillations lead to that the χ field may have a tachyonic mass since the square of its effective mass, which has the form can be negative. The tachyonic instability causes the parametric resonance of the χ particles, which excites the production of copious χ particles of small-momentum modes.
Since there are rescatterings between the χ particles and the inflaton condensate, the produced χ particles can knock inflaton quanta out of the condensate into low-momentum modes. So, the rapid growth of the inflaton and χ fluctuations can be expected, which means that they can become significant GW sources. The best way to obtain the production of GW during preheating in the Starobinsky model is to perform numerical lattice simulations. Here, we use the publicly available package HLattice [42] to do simulations. In original HLattice the symplectic integrator is used to handle scalar fields with canonical kinetic terms. Since the scalar fields considered in this paper have non-canonical kinetic terms, we modify the HLattice by adopting the fourth-order Runge-Kutta integrator instead of the symplectic one. The results are simulated with N 3 = (128) 3 points. The lattice length of side L is chosen to to satisfy L H i < 2π , which means that all field modes are always within the horizon, where H i ( 3.87 × 10 −6 M p ) is the Hubble parameter when the preheating begins. Different values of L H i are chosen for the different values of ξ . For example, when ξ = 5 we set L H i to be 5, which indicates that the minimum of the comoving wave number k is k min = 1.25H i . Since all field modes are sub-horizon, we will initialize the fluctuations of φ and χ fields to be those in the Bunch-Davies vacuum. In addition, we set a UV cutoff k UV on k when initializing the fluctuations, i.e., k UV = 30k min when ξ = 5, which is larger than the maximum resonance frequency in order to prevent it from affecting the results. Furthermore, we have checked that the results are independent of the values of k UV . In the original HLattice, the produced results are also limited to be in the regions k ∈ [k min , k UV ].
Here, we relax this limitation to include some momenta larger than k UV . This operation is similar to what was done in [43]. The initial conditions of the homogeneous part of φ field during preheating are determined by the time when the inflation ends, while for the χ field the homogeneous part is initialized to zero. We stop the simulation when the density spectrum of the inflaton quanta and the χ field do not change appreciably.

Production of gravitational waves
Now we study the production of GWs during preheating in the Starobinsky inflationary model. GW can be represented by the transverse traceless part of the spatial metric disturbance of the FRW background The perturbation h i j corresponds to two independent tensor degrees of freedom and satisfies the equation of motion where the source term T T i j is the transverse-traceless part of the anisotropic stress i j , which is defined as i j ≡ T i j − p g i j (10) with T i j and p being the energy-momentum tensor and the pressure of the system, respectively. For the model considered in this paper the energy-momentum tensor has the form and the anisotropic stress can be expressed as From the above expression, one can obtain the transversetraceless part of the anisotropic stress In our numerical calculation, it is more convenient to work in the Fourier space. Using the convention one can obtain that the GW equation (9) in Fourier space has the form where k = |k|. After introducing the transverse-traceless projection operator [25]: is the spatial projection operators, we achieve the transverse-traceless part of i j in the momentum space It is easy to see that k i T T i j (t, k) = T T ii (t, k) = 0 in any time t. Using this projection operator, one can define a new tensor u i j which satisfies the following relation in the momentum space Then the GW equation (9) can be re-written as We assume no GWs at the beginning of preheating, so u i j and its derivative are initialized to be zero. The GW energy density is given by Here · · · is the spatial average. According to the Parseval's theorem, the GW energy density can be expressed as an integral in momentum space where L is the length of one side of the lattice. The corresponding GW spectra, normalized to the critical energy density ρ c , can be obtained through Since the produced GWs during preheating are determined by the spectra of scalar fields, we first investigate the evolutions of χ and φ and show their spectra in Fig. 1 where only ξ = −5 and 5 are considered as examples. It is easy to see that due to the parametric resonance χ particles are produced rapidly. When ξ = −5 the production is quicker than when ξ = 5, which results from that all modes of the χ field with k 2 /a 2 < |m 2 χ,eff | are already tachyonic when the preheating begins [39]. After the increase of the χ spectrum, the spectrum of the φ field begins to increase since inflaton particles are knocked out of the inflaton condensate through the rescattering between χ particles and the inflaton condensate. Finally, both the inflaton field and the auxiliary field become the effective GW sources. Figure 2 shows the evolutions of the GW density spectra from the lattice simulation with ξ = 5, 10, 20 and 30. The outputs are plotted from bottom to top per same time step t = 0.0625H −1 i , where H i denotes the value of Hubble rate when inflation ends. One can see clearly that the low frequency modes increase rapidly as parametric resonance begins, and then the growth of the higher frequency modes is very fast due to the nonlinear effect. The amplitude increases faster and faster with the increase of ξ , which results from the fact that the resonant bands become broader and broader with the increase of ξ . The blue curves, which show the final GW spectra, indicate that the maximum value of gw (k), gw,max , is not a monotonous function of ξ , since when ξ = 10 the value is larger than the ones of ξ = 5 and ξ = 30, which is different from the case of non-minimally coupled inflation model in which the GW spectra produced in preheating increase with the increase of the coupling parameter [35]. Figure 3 give the evolutions of the GW density spectra for the case of negative ξ . We take four different values: ξ = −5, −10, −20 and −30. The results are very similar with the ones shown in Fig. 2. This is because the resonant bands become broad with the increase of |ξ |. But in the case of ξ < 0, gw begins to grow once the inflation ends and grows faster than the case of ξ > 0. The reason is that when ξ < 0, all modes of the χ field with k 2 /a 2 < |m 2 χ,eff | are already tachyonic when the preheating begins [39]. In addition, we find when Although Figs. 2 and 3 do not indicate a simple relation between gw,max and the coupling constant ξ , there, however, indeed exists a subtle relation as we will now demonstrate. For this purpose, we plot the evolution of gw,max as a function of ξ in Fig. 4. We find that when ξ is very small the production of GWs is very weak, which is agreement with the result obtained in [45] where it was found that the direct generation of GWs after the end of inflation is suppressed in the Starobinsky model. One can see that the maximum of gw is about 6.67 × 10 −3 , which occurs at ξ −20 . In the case of ξ > 0, the maximum of gw is about 5.46 × 10 −3 at ξ 42. It is easy to see that gw,max first increases fastly with the increase of |ξ |, and then reaches the maximum after several oscillations. Finally it decreases with slight oscillation with the further increase of |ξ |.
Since we are interested in the energy density of GWs today, the GW spectra generated in the preheating era need to be transformed to the present-day case. The present scale factor compared to that at the time when GW production stops can be expressed as where subscripts * , j and 0 denote the time when GW production stops, thermal equilibrium is established and today, respectively. Here, w is the mean equation of state of the cosmic energy from t * to t j , ρ r is the radiation energy density, g is the number of effectively relativistic degrees of freedom, and ρ * is the total energy density of scalar fields. So the present-day GW physical frequency is given by Here, we take the reheating temperature to be about 3.1 × 10 9 Gev [44], which means that one can takeḡ j /ḡ 0 = 106.75/3.36 31 in analysis. The transformed function to obtain today's GW amplitude is [32] gw,0 ( f )h 2 = gw ( f ) During reheating the mean equation of state w varies with time, but its form is unknown. In our analysis we assume w = 1/3 between t * and t j for simplicity. This assumption is reasonable since in [39] it has been found that after preheating the value of w is close to 0.3. Since the GW energy density scales like radiation, it may be constrained by CMB and BBN measurements of the total radiation density in species beyond the standard model [46]. If we assume that beyond the standard model all extra radiation density today during the formation of the CMB is comprised of GWs, the total GW energy density is constrained from a bound on N eff [47] gw,0 h 2 γ,0 h 2 = where γ,0 h 2 = 2.47 × 10 −5 is the present energy density of photons, N eff is the effective extra relativistic degrees of freedom and N eff = N eff − 3.046. The Planck 2018 results [48] limit | N eff | 0.33, which requires that the GW energy density must satisfy gw,0 h 2 1.85 × 10 −6 . The next-generation CMB experiments, such as CMB-S4, will probe N eff ≤ 0.03 at 1σ CL and N eff ≤ 0.06 at 2σ Fig. 5 The total GW energy density today as a function of ξ CL [49], which potentially constrains the energy density to be gw,0 h 2 1.68 − 3.36 × 10 −7 (27) In Fig. 5, we show the total GW energy density today as a function of ξ . One can see that the maximum of energy density is 1.09×10 −7 , which occurs at ξ −21. Thus, the GWs produced during preheating in the Starobinsky inflationary model satisfy the Planck limit and the next-generation CMB experiment constraints.

Conclusion
The Starobinsky inflationary model satisfies the latest CMB observations very well. In this inflationary model, to achieve the preheating process requires an auxiliary scalar field which couples non-minimally with the gravity. This auxiliary scalar field becomes an effective GW source due to the parametric resonance. Furthermore, the rescattering between auxiliary scalar field and inflation field can knock copious inflaton particles out of the inflaton condensate, which results in that the inflaton field also becomes an effective GW source. In this paper, using the lattice simulation, we study the production of GWs during preheating in the Starobinsky model. We find that the GW spectrum gw grows fast with the increase of the absolute value of the coupling parameter ξ . This is because the resonant bands become broad with the increase of |ξ |. When ξ < 0, gw begins to grow once the inflation ends and grows faster than the case of ξ > 0. gw reaches the maximum at ξ −20 (ξ 42 for the case ξ > 0) and then decreases with slight oscillation. By assuming that all extra radiation density during the CMB formation is comprised of GWs, we investigate the constraint on the total GW energy density from observations and find that the GWs producing during preheating in the Starobinsky mode satisfy the limits from the Planck and next-generation CMB experiments. and Technology Innovation Plan of Hunan province under Grant No. 2017XK2019.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The published work is theoretical and so no data should be deposited.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .