More on loops in reheating: Non-gaussianities and tensor power spectrum

We consider the single field chaotic m^2\phi^2 inflationary model with a period of preheating, where the inflaton decays to another scalar field \chi\ in the parametric resonance regime. In a recent work, one of us has shown that the \chi\ modes circulating in the loops during preheating notably modify the<\zeta\zeta>correlation function. We first rederive this result using a different gauge condition hence reconfirm that superhorizon \zeta\ modes are affected by the loops in preheating. Further, we examine how \chi\ loops give rise to non-gaussianity and affect the tensor perturbations. For that, all cubic and some higher order interactions involving two \chi\ fields are determined and their contribution to the non-gaussianity parameter f_{NL} and the tensor power spectrum are calculated at one loop. Our estimates for these corrections show that while a large amount of non-gaussianity can be produced during reheating, the tensor power spectrum receive moderate corrections. We observe that the loop quantum effects increase with more \chi\ fields circulating in the loops indicating that the perturbation theory might be broken down. These findings demonstrate that the loop corrections during reheating are significant and they must be taken into account for precision inflationary cosmology.

The constancy of the superhorizon curvature perturbation ζ is very crucial for the inflationary predictions to hold. This helps to determine the cosmic microwave background (CMB) fluctuations from the correlation functions evaluated at the horizon crossing time. Technically, the conservation of ζ sets an upper bound for the time integrals that appear in the in-in perturbation theory [6].
On the other hand, it is well known that the entropy perturbations can cause superhorizon evolution of the curvature perturbation (see e.g. [7]). In reheating, those fields that are unimportant during the exponential expansion are exited by the inflaton decay and they start to dominate the universe. Moreover, in some models the decay can occur violently in a preheating stage [8][9][10][11]. Although these entropy perturbations are not produced at cosmologically interesting scales, reheating stage ends with highly nonlinear processes (see e.g. [12][13][14]). While these nonlinearities can be effectively described by fluid dynamics that only affect local quantities [15], some of them are known to have important consequences (see e.g. [16,17]).
Classical, and similarly quantum, nonlinearities imply that Fourier modes do not evolve independently. As a result of this mode-mode coupling, short distance fluctuations are expected to affect the long wavelength modes. For example, a cubic interaction term in a Lagrangian would allow two modes with nearly equal large momenta to change the amplitude of a mode with small momentum. In quantum theory, there are also virtual modes circulating in the loops that affect the correlation functions. Evidently, it is crucial to determine the size of such effects. In a recent work, one of us has shown that in the chaotic m 2 φ 2 model with a period of preheating where the inlaton φ decays to another scalar χ, the parametrically amplified χ modes appearing in the loops would meaningfully modify the curvature power spectrum [18]. This is an example of the entropy perturbations affecting the superhorizon curvature variable, however not by the real physical fluctuations but because of the virtual entropy modes appearing in the loops (entropy modes are known to give power loop infrared divergences during inflation [20,21]).
Since ζ becomes an ill defined dynamical variable during reheating, the calculations in [18] have been carried out in the ζ = 0 gauge. In that case, one may first calculate the χ loop corrections to the inflaton fluctuation power spectrum until the coherency of the inflaton oscillations is lost.
After that moment, the possible effects on the superhorizon evolution are expected to be averaged out and become negligible. One may then apply a gauge transformation to read the ζ power spectrum. Since in the first stage of the preheating the background inflaton oscillates coherently, the superhorizon modes are affected without violating causality [22]. Note that as long as the relativistic equations are treated properly, there should not arise any issue with causality.
In this paper, our first aim is to carry out the calculation of [18] in the ϕ = 0 gauge, i.e. we will use ζ directly as the main dynamical variable. Because ζ is only ill defined at isolated times when the inflaton velocity vanishes, the propagator has "spikes" and it diverges at these moments. We smooth out these spikes by using the time averaged background quantities in the ζ action. In [18], only the loops arising from the interaction potential have been considered. Here, we determine all cubic interactions involving χ and ζ, and estimate the total one loop correction to the ζ power spectrum. Not surprisingly, our computations confirm the findings of [18] and show a significant contribution to the ζζ correlation function.
Our second aim in this work is to determine the χ loop contributions to the non-gaussianity and to the tensor power spectrum. Notable modifications to the scalar power spectrum found in [18] indicate the existence of similar significant corrections for these observables. We estimate the amount of non-gaussianity from the ζ-three point function by calculating the one loop graphs arising from the cubic interactions. It turns out that these corrections to the three point function can be expressed in terms of the two point function and it is possible to read the shape independent nongaussianity order parameter f N L . Similarly, the χ field coupled to the tensor fluctuations yield loop corrections to the tensor power spectrum. Since the tensor field behaves like a test field propagating on the background, the tensor calculation is not affected by different time slices of spacetime. We find that the tensor power spectrum is moderately corrected by the loops in reheating.
The organization of the paper is as follows. In the next section, we consider the chaotic m 2 φ 2 model with the extra preheating scalar field χ to which the inflaton decays in the parametric resonance regime. In section III, we determine the cubic interaction terms involving the curvature perturbation ζ, the tensor mode γ ij and the preheating scalar χ. We then calculate one loop corrections to the scalar power spectrum ζζ , the three-point function ζζζ , the tensor power spectrum γγ and make order of magnitude estimates of these corrections by using the theory of preheating. In section IV, we further consider some higher order interactions involving two χ fields and determine their loop effects. In V we conclude with remarks and future directions.

A. The background
Let us consider the chaotic model that has the following potential where φ is the inflaton and χ is the reheating scalar, which are propagating in a flat FRW background ds 2 = −dt 2 + a 2 (dx 2 + dy 2 + dz 2 ). (2) This model can be seen to be the prototype of the chaotic inflationary paradigm and preheating. As it is well known, a period of inflation can be realized if initially φ >M p and the nearly exponential expansion ends roughly when φ ∼M p /20 (see e.g. [11]), whereM 2 p = G (we define M p to be the reduced Planck mass M 2 p = 8πG/3). During inflation and in the first stage of preheating where the backreaction effects are negligible, the χ background vanishes Following the exponential expansion, φ starts oscillating about its minimum φ = 0. Assuming which is generically satisfied in this model, the background field equations can be approximately solved as where Note that the amplitude obeysΦ + 3HΦ/2 ≃ 0, where the dot denotes the time derivative.
We define t R and t F as follows: t R : Beginning of reheating, t F : End of the first stage of preheating.
After the time t F , the χ particles created out of the vacuum start affecting the background and thus the backreaction effects are set in. Our aim is to calculate the χ loop corrections to the cosmological correlation functions, which are effective in the time interval (t R , t F ).
Some features of preheating depend on the parameters of the model and many cases are discussed numerically in [11]. For our estimates, we will use the following canonical set that gives the broad parametric resonance: In that case, the first stage of preheating ends after about 11 inflaton oscillations and one has [11] One may also note that which can be determined from m/H(t F ) = 3mt F /2 ≃ 100 and the fact that the first stage ends after 11 oscillations. The initial amplitude in (6) is given by Φ 0 ≃M p /20.
In the physical momentum space, the first resonance band is given by q phys ∈ (0, q * ) where In general, there are other resonance bands which can be important for preheating [11]. In the model we are studying, the first instability band gives the largest contribution and in the following we simply underestimate the loop corrections by neglecting the effects of other resonance bands.
The χ momentum modes sitting in the band (0, q * ) encounter exponential amplification. In determining q * , we will use the smallest value of Φ, i.e Φ(t F ) in (8). The first stage ends when the interaction potential energy density g 2 φ 2 χ 2 becomes comparable to the inflaton potential energy density m 2 φ 2 , since after that moment the frequency of the inflaton oscillations are affected by the χ particles. This implies [11] As we will see below, (11) is important for estimating the χ loop corrections.
Since the background value of χ is zero, we have θ = 0, σ = φ and δs = δχ, which shows that in this model φ is the adiabatic mode and χ is the entropy mode.

B. Quadratic actions and mode functions
The full action governing the dynamics of the system can be written in the ADM form as (we where N and N i are the standard lapse and shift functions of the metric We define the perturbations as where the gauge is completely fixed by imposing Here, ϕ denotes the inflaton fluctuation and in this gauge the inflaton takes its background value φ = φ(t) given in (5). Note that we use the same letter χ to denote the reheating scalar fluctuation in (20) since the background value of χ vanishes. As pointed out in [23], the lapse N can be solved exactly as N 2 = B/A. However, to determine the action up to cubic order it is enough to solve the constraints to linear order, which gives [6] Note that neither χ nor γ ij appear in the solutions of N and N i to this order. By expanding the action (15), one may obtain the following well known quadratic actions which are valid both during inflation and reheating. The ζ kinetic term vanishes at times wheṅ φ = 0 and the ζ propagator diverges at these times. This divergence must be cured to make the loop contributions well defined.
The free fields can be expanded as where s = 1, 2 and the ladder operators obey the usual commutator relations, e.g. [a k , a † k ′ ] = δ 3 (k − k ′ ). The polarization tensor ǫ s ij has the following properties To satisfy the canonical commutation relations, the mode functions must obey the Wronskian On the other hand, the linearized mode equations becomë Note that the equation for χ k gets a contribution from the potential (1), which is responsible for the parametric resonance.
We will be interested in the superhorizon ζ k and γ k modes. Neglecting the k 2 /a 2 terms in (28) one can easily obtain two linearly independent superhorizon solutions which can be written as where ζ k , c k and d k are constants and As usual, the modes (29) have the constant and the decaying pieces, and the normalization conditions in (27) imply One may note the mass dimensions 1 of the constants as [ζ To be able to calculate the χ loop effects, we need to determine the behavior of the χ modes, especially the ones in the resonance band, in detail. For that, one may write the mode function in the WKB form as follows where The Wronskian condition is satisfied by imposing |α q | 2 − |β q | 2 = 1. During inflation, χ becomes a very massive field with mass gΦ 0 . As a result, for the modes of interest the Bunch-Davies mode function in the beginning of reheating can be written up to an irrelevant phase as This shows that at the end of the exponential expansion these individual χ q modes are suppressed by a −3/2 and this is the main reason for the metric preheating scenario of [22,24,25] to break down, as it is discussed in [26][27][28] (it is possible to circumvent this suppression in some models, as it is shown in [29][30][31]). During preheating, χ q changes non-adiabatically as the inflaton passes through the potential minimum φ = 0. This process can be formulated as the particle creation by parabolic potentials which gives the exponential increase β q = e µqmt for the modes in the instability bands, where µ q is an index characterizing the exponential growth.
From (32) one may find that For |α q | ≃ |β q | ≫ 1, it is possible to see that |χ q | 2 oscillates between 1/(2a 3 ω q ) and 4|β q | 2 /(2a 3 ω q ) with the frequency ω q . To determine the phase of χ q , one may define θ q as Then, the Wronskian condition (27) gives i.e. up to an unimportant constant the phase is uniquely fixed by the amplitude |χ q |.
The growth of the modes in the first instability band can be described by introducing an effective index µ q ≃ µ, and for the parameters given in (7) one has [11] µ ≃ 0.13.
To estimate the magnitude of the amplitude |χ q | at the end of the first stage of preheating, one may look at the expectation value χ 2 , which is given by where in the last equality we restrict the momentum integral to the first (and the most important) instability band, which is supposed to give the dominant contribution to the vacuum expectation value; we switch to the physical momentum space and introduce |χ q * | to denote a mean value for the modes in this instability band. Note that the 4π factor in (39) comes from the angular directions in the momentum space. Comparing with (11) one may deduce that at the end of the first stage As pointed out above, the amplitude |χ q | is actually an oscillating function that has frequency ω q .
However, one has ω q * ≫ m and thus |χ q * | oscillates much faster than the background inflaton field.
As a result, (40) should be divided by 2 to give a time averaged value for the amplitude. We also use the index µ to obtain the amplitude in the middle of the period and define The phase corresponding to (41) can determined from (37) as These estimates will be crucial in determining the strength of a graph in the in-in perturbation theory.
We define the scalar and the tensor power spectra in the momentum space, i.e. P ζ k and P γ k , from the two point functions in the from where the polarization tensor Π ijkl , which is defined as obeys Π ijkl Π klmn = 2Π ijmn . The tree level standard results can be read from (29) as The constants ζ k can be determined from the mode functions of the free fields during inflation and as it is well known they depend on the horizon crossing time for a given k (see [32] for a study of loop corrections to the mode functions during inflation).

C. Smoothing out spikes of ζ
In finding f (t) from (30), an infinity arises when the limits of the integration contains a moment givingφ = 0. To avoid these singularities one may try to fix f (t) by an indefinite integral since one only needs a function whose derivative gives (30). However, the function obtained in this way is unavoidably singular at times whenφ = 0. Moreover, the loop corrections turn out to involve the time integrals of f (t) or df /dt, and these also diverge when f (t) obeys (30).
This pathologic behavior arises due to the bad choice of gauge. 2 Namely, ϕ = 0 gauge breaks down at times whenφ = 0 giving rise to the spikes of ζ. This has already been noted in some earlier work, see e.g. [17,33]. As discussed in [33], although ζ becomes an ill defined variable in reheating, (1 + w)ζ becomes well defined, where w is the equation of state parameter. In our model To smooth out the spikes of ζ, we first note that the Einstein's equations for the background Since we use (6) to approximate the Hubble parameter H, one may defineφ 2 av by using (6) in (47) that yieldsφ It is clear thatφ 2 av gives the "time" average of the oscillating functionφ 2 . To make ζ well defined, one may now replaceφ 2 byφ 2 av in the free action of ζ in (23). In the context of the discussion carried out in [33], this is equivalent to using an average equation of state parameter w av instead of the actual one. Consequently, one simply treats the ζ variable as if it evolves in a matter dominated universe. In that case, the new function obeys A simple integration then gives where we use (49) and (30) for f (t) and g(t), respectively. As we will see below, the loop contributions turn out to depend on the difference of two f (t) or the difference of two g(t) functions, and therefore there is no need to fix the integration constants in (50).

III. CUBIC INTERACTIONS AND LOOP CORRECTIONS
Using (20) and (22) in (15), a straightforward calculation gives the following cubic action involving two χ fields: Combining this cubic action with the quadratic ones given in (23) and switching to the Hamiltonian formulation, one may find the cubic interaction Hamiltonian containing two χ fields as where Although it is not indicated explicitly, all the fields appearing in (52) can be taken to be the interaction picture fields that enter in the in-in perturbation theory as it is formulated in [34]. In obtaining (52) we only spatially integrate by parts the last term in the first line of (51) to replace the shift N i by its potential ψ given in (22).
For any given operator O, the in-in formalism can be applied to obtain the following perturbative expansion for the vacuum expectation value [34] where the lower limit of the time integrals is set to t R rather than −∞ since we are interested in the loop effects during reheating. In general, the two terms in a given commutator in (56) (10). This yields an extra enlargement factor of a 3 .
Using (56) for the operator ζ(t, x)ζ(t, y) with N = 2 gives the following vacuum expectation values of the nested commutators: Using (29) in these commutators, we see that for superhorizon modes the first commutator gives Defining the function F (t 1 , t 2 , k) by and using (29), (31) and (43), one can determine the largest correction as where k denotes the comoving cosmological superhorizon scale of interest and t F marks the end of the first stage of preheating as defined above. It is remarkable that the one-loop correction P ζ k (t F ) (1) becomes a multiple of the the tree level function P ζ(0) k given in (46). From (53) and (54), where only the contribution of the first terms in (53) and (54) are written explicitly.
The momentum integral in (65) diverges and it must be regularized before making any order of magnitude estimates. To figure out the contribution of the modes in the resonance band and for regularization, we simply cutoff the integral in (65) with aq * , where q * is given by (10). Indeed, this corresponds to the adiabatic regularization where one uses the WKB mode function (32) and discard the pieces with α q that give infinities. Note that initially one has α q (t R ) = 1 and β q (t R ) = 0; and β q increases with time in the resonance band and stays vanishingly small for high energy modes since they propagate adiabatically. From (7) and (8) one sees that q * ≃ 10 −6M p ≪M p .
Consequently, in a legitimate renormalization procedure that is more systematic than the simple adiabatic regularization, the UV subtractions should not change our estimates appreciably since the cutoff scale q * corresponds to a relatively low energy scale.
The correction (64) modifies not only the amplitude but also the index of the power spectrum, but the change in the index is negligible for cosmologically interesting scales 3 since in that case k ≪ aq * . Therefore, the k dependence of (64) is negligible and to a very good approximation one may ignore it by setting k = 0.
As discussed above, a(t 1 ) 3 and a(t 2 ) 3 terms cancel out the scale factor suppressions of the four χ q modes. The 1/a 3 factor that appears in df /dt in (49)  In what follows we estimate (64) to determine the size of the loop effects in reheating. We first focus on the term that is explicitly shown in (65) and then confirm that others give similar contributions. Since the resonant χ modes encounter most of their growth near the end of the 3 The index is meaningfully modified for the modes entering the horizon during reheating that may change the primordial black hole formation, see [35].
first stage, one may focus on the last inflaton oscillation for the time integrals in (64), namely, the lower and the upper limits can be set to mt F − 2π and mt F , respectively. Using (36), the square brackets in (64) yields the following factor We see that the leading order contribution does not cancel out since the phase factors have different time arguments. In (65), there are four χ q modes integrated out in the first instability band, which can be estimated as q 3 * |χ q * | 4 , where |χ q * | is the mean value of the modes introduced in (41). The function df /dt can be read from (49). Treating the slowly changing factors like H and Φ as constants one finally finds that where the dimensionless constant C 1 is given by Recall that the phase θ q * is defined in (42).
For our set of parameters (7), the constant C 1 can be determined by a numerical integration that yields C 1 ≃ 0.078. Using then (7), (8), (10) and (41) in (67), we obtain which is consistent with the estimate given in [18].
Let us now consider the contributions of the other terms in (65), which can be determined from the definition (63). From (53) and (54), these consist of the products of four χ fields, on which certain time or spatial derivatives act (there is also a nonlocal term with 1/∂ 2 that involves the Green function of the Laplacian). In (64), only the imaginary part of F (k 1 , k 2 ) appears in the square brackets. One can easily see that after taking the imaginary part, each product yields a term similar to (66) and thus the leading order contributions do not cancel out. On the other hand, the time integrals are very similar to (68) and they can all be estimated to give C/m 2 . One may also note that a partial derivative ∂ i would produce q i in momentum space andχ q ≃ ω q χ q , where ω q is given in (33). Therefore, to estimate the size of a correction one may simply replace g 2 Φ 2 factor in the second square bracket in (67), which arises due to g 2 φ 2 terms in (53) and (54), by ω 2 q * corresponding toχ 2 or q 2 * corresponding to (∂χ) 2 (note that 1/a 2 factor, which multiplies (∂χ) 2 in (53) and (54), converts the comoving momentum scale arising from the spatial partial derivative to the physical momentum scale). Similarly, the magnitudes of the nonlocal terms can be estimated by using the Green function for the Laplacian and the correlation length corresponding to the χ fluctuations, which is roughly equal to 1/q * as shown in [36]. In all these different cases one may see that the contributions have the same order of magnitude with (69), since for our numerical choice of parameters (7) one has gΦ ∼ q * ∼ ω q * . The sign of each contribution depends on the phases through the expressions like (66), which is sensitive to the initial conditions [11].
Because the one loop correction (69) is larger than the tree level result, the in-in perturbation theory might be broken down in this model. Since the modes of the χ field is exponentially amplified during preheating, the quantum corrections are enlarged when more χ fields circulate in the loops.
As we will see, the results of the next section will support this expectation, i.e. the lower order loop corrections that are supposed to give larger contributions than (69) become smaller due to the less number of χ modes circulating in the loops. A similar situation also arises for f N L as we will discuss in the next section.
The function P (k 1 , k 2 , k 3 ) measures the size of the non-gaussianity involving the comoving superhorizon scales k 1 , k 2 and k 3 that obey k 1 + k 2 + k 3 = 0. To pin down the loop corrections one may use (56) for ζ(t, x)ζ(t, y)ζ(t, z) and since H I is linear in ζ the first nonzero contribution arises for N = 3, which gives the diagram in Fig. 2.
one may straightforwardly express the leading order one loop correction in terms of G(k 1 , k 2 , k 3 ) where the extra two terms, which can be obtained by cyclic interchange of momenta, are not written explicitly.
Using (53) and (54) in (71), it is possible to express G(k 1 , k 2 , k 3 ) as a loop momentum integral of the mode functions. Indeed a straightforward calculation gives where the contributions of the first terms in (53) and (54) are expressed explicitly. If q denotes the loop variable that is restricted to the instability band (0, aq * ), again one has q ≫ k 1 , k 2 , k 3 .
Since the modes in the loop integral in (73) become functions of q + k i , i.e. χ q+k i , the dependence of G(k 1 , k 2 , k 3 ) on its arguments is very weak and one may write G(k 1 , k 2 , k 3 ) ≃ G. Using (36) we obtain |G| ≃ d 3 q 36 (2π) 9 Since the largest contribution to this loop integral comes when q runs near aq * , one may set q = aq * and use d 3 q → 4πq 3 * to estimate the integral. It is now possible to use (74) in (72) to read the three point function. As before, the largest contribution to the time integrals come from the last oscillation period in which χ modes are amplified most. Keeping the slowly changing factors like Φ and H as constants in this last cycle, we obtain where the dimensionless constant C 2 is given by We would like to recall that in this expression the scale factors cancel out each other and the time dependent dimension-full quantities are evaluated at the end of the first stage of preheating.
The non-gaussianity parameter f N L can be defined as [6,37] where ζ g denotes the corresponding free quantum field. This definition introduces a shape independent parameter that gives an overall order of magnitude estimate for the scalar non-gaussianity.
Calculating the three point function by using (77) and comparing with (75) one finds (78) For our canonical set (7), C 2 can be found by a numerical integration that gives C 2 ≃ 0.00057 (recall that θ q * is fixed in (42)). Using the values of other dimension-full parameters in (78) we This is a very large amount of non-gaussianity that is solely produced in reheating and it is obviously inconsistent with observations. On the other hand, by comparing (69) and (79) we observe that although they measure different one loop corrections, the latter has more χ modes circulating in the loops and it produces a much bigger number. Therefore, the large amount obtained in (79) can be an artifact of perturbation theory, which might become invalid in this model. It is possible to produce large non-gaussianity in inflationary models (see e.g. [38]), but the single scalar field models generically give f N L = O(ǫ), where ǫ is the slow roll parameter. Although we are not capable of making non-perturbative estimates, our computations show that a large non-gaussianity can be produced during reheating.

C. The tensor power spectrum
The interaction Hamiltonian (52) also modifies the tensor power spectrum due to the last term involving the graviton coupling. One may first think that this interaction is suppressed by 1/a 2 , however this factor simply converts the two comoving momenta arising from the two partial derivatives to the physical scale. The tensor field γ ij is similar to a spectator field since its background value vanishes. As a result, the tensor power spectrum is not affected by the (infinitesimal) changes of the spacetime slicing and the gauge can be fixed in a natural way without giving rise to any complications. Moreover, unlike the ζ propagator, the tensor propagator does not contain any singularities. The correction corresponding to (52) can be pictured as in Fig. 3.
Using (56) for γ ij (t, x)γ kl (t, y) with N = 2, which gives the first nonzero contribution, and applying the identity [AB, C] = A[B, C] + [A, C]B, one finds terms with single or two graviton commutators. It is easy to see that the terms with two graviton commutators are suppressed by 1/a 3 and hence they become completely negligible. A straightforward calculation then gives the following one loop correction to the tensor power spectrum in momentum space where g(t) is defined in (30) and In (80) we reintroduce the Planck mass M p , which can be fixed either by dimensional analysis or by keeping track of its presence starting from the action (15). Once again, the one loop correction in momentum space becomes a multiple of the tree level power spectrum. This is mainly because of the fact that the expectation value [O ij (t 2 , z 2 ), O kl (t 1 , z 1 )] , which appears due to last term of the interaction Hamiltonian (52), produces δ ik δ jl + δ il δ kj and this index structure acting on the polarization tensor Π ijkl , which is introduced in (45), gives the same tensor.
Converting the comoving integration variable in (81) to the physical scale generates the power a 7 , and this factor together with a(t 1 )a(t 2 ) in (80) completely compensate the suppressions of the mode functions χ q and the 1/a 3 decay of the function g(t). As before, the change in the spectral index is negligible due to the large hierarchy between the superhorizon scale k and the scale q * characterizing the instability band. Therefore, in (80) one may ignore the k dependence, set q = q * and let d 3 q → 4πq 3 * . For the χ modes, one may use (36) and (42). Finally, to estimate the time integral, we introduce the time dependence of the background quantities using (6). As a result we where For our canonical set of parameters (7), we numerically integrate (83) that yields C 3 ≃ 0.029.
Using (8) for the Hubble parameter one finds The main reason for this correction to be tiny compared to the scalar power spectrum (69) is that the tensor corrections are suppressed by the Planck mass M p . Nevertheless, the modification (84) is much larger than the quantum corrections that arise during inflation, which are suppressed by the ratio H/M p [6].

IV. SOME HIGHER ORDER INTERACTIONS AND LOOPS
The results of the previous section show that cubic interactions involving two χ fields modify the scalar and the tensor power spectra and give rise to non-gaussianity. Although the interaction Hamiltonian (52) is cubic, the first nonzero contributions come from (56) with N = 2 for the scalar and the tensor power spectra, and with N = 3 for the three point function. The corresponding one loop corrections are sixth and ninth order in fluctuations, respectively.
In this section, we consider some higher order (e.g. fourth and fifth order) interactions, again involving two χ fields, and calculate the corresponding one loop effects. Our aim in considering such interactions is two fold. First, we would like to use (56) with N = 1. Therefore, by a naive counting in perturbation theory the effects are supposed to be more prominent than the ones we have studied in the previous section (although this turns out to be incorrect as we will see below).
Second, the loop effects calculated in the previous section involve the commutators of the χ fields and thus one must carefully treat the phase factors as we did in (66). The loop corrections we consider in this section demonstrate the modifications more directly.
A. The scalar power spectrum and non-gaussianity Starting from the action (15), one may obtain the following terms in the interaction Hamiltonian where ρ χ is the energy density of the χ field given by The first term in (85) contributes to the scalar power spectrum and the second one produces scalar non-gaussianity. Note that the linear ζ term in (85) agrees with the cubic hamiltonian in (52).
Let us first consider the one loop correction to the scalar power spectrum arising from (85) that can be pictured as in Fig. 4. Using (56) for the ζ(t, x)ζ(t, y) with N = 1, a straightforward calculation gives This equation clearly shows how the correction enlarges in time during preheating as the energy density ρ χ increases as a result of χ particle creation. Note that (87) only modifies the amplitude of the spectrum since the correction multiplying the tree level result does not depend on the external momentum k. At the end of the first stage of preheating the energy density of the created χ particles catches up the background energy density, which gives ρ χ (t 1 ) ≃ 3H 2 M 2 p . Reading f (t) from (50), it is easy to see that Indeed, using (6) for the background quantities one finds that where, as before, we restrict the time integral to the last inflaton oscillation cycle.
One may find other terms in the interaction Hamiltonian that modifies the scalar power spectrum. For instance, by introducing exp(3ζ) factor in (52), which arises from √ h, one obtains a fourth order term H (4) After using (90) in (56) with N = 1, one encounters terms either with ζζ [ζ, ζ] or ζζ [ζ, ζ]. It is easy to see that the latter is suppressed by 1/a 3 and the former yields From (53), one has O 1 ≃ ρ χ /H and using (49) we obtain The main conclusion here is that although the corrections (88) and (92) correspond to lower order in perturbation theory, they give smaller contributions compared to (69).
The fifth order term ζ 3 ρ χ in (85) corrects the three point function and thus it gives rise to non-gaussianity. The corresponding graph is pictured in Fig. 5. Using (56) for ζ(t, x)ζ(t, y)ζ(t, z) with N = 1 and using the definition of the three point function in momentum space given in (70), one finds that From (77), the corresponding f N L parameter can be calculated as As in (90), by introducing exp(3ζ) factor in (52) gives the following interaction Hamiltonian: It's contribution to f N L can be found as In both of these cases it is easy to estimate the integrals so that Therefore a small amount of non-gaussianity is produced by these interactions. As in the case of the power spectrum, the loop corrections to f N L coming from the interactions that can be pictured as in Fig. 5 become much smaller than the previous one (79).
B. Fourth order interactions that has the form γγχχ and the tensor power spectrum Till now in this section we have considered some higher order interactions that modify the scalar power spectrum and the f N L parameter. It is clear that in a systematic study one should work out the complete fourth order action to determine the corrections more accurately. In that case, the lapse N and the shift N i must be solved up to second order. This is a complicated calculation and the complete fourth order action is not very illuminating for the scalar field. However, the interactions studied above are generic enough to indicate that other corrections to the scalar power spectrum and f N L will be similar to the ones found above.
In this subsection, we determine the complete fourth order action involving the interactions of the tensor field γ ij and the reheating scalar χ. Our aim is again to compare the corresponding corrections with (84) to see how the perturbation theory is working. Since we solely concentrate on the tensor modes we set (recall that we have been working in the ϕ = 0 gauge). The quartic interactions involving γ ij and χ are necessarily in the from γγχχ since the background values of γ ij and χ are zero. Similarly, there is no linear term in N and N i after one sets ζ = 0. We define where ∂ i N (2)i T = 0. To determine these second order quantities, one may use the exact solution for the lapse where A and B are defined in (17) and (18), and work out the momentum constraint, which reads where K ij and K are defined above (17). Up to second order in fluctuations, the Ricci scalar R (3) of the constant time hypersurface can be found as After a relatively long but straightforward calculation we find Similarly, the transverse part of the shift reads In all these expressions the indices are contracted with the Kronecker delta and we set M p = 1.
Before discussing the loop corrections, it is interesting to check the validity of the perturbation theory from the quadratic expressions given for the lapse and the shift. As discussed in [17], the perturbation theory is applicable if one has While the first condition is needed for keeping the time coordinate to be proper, the second ensures that the original foliation of the spacetime that is presumed for perturbation theory is not destroyed by the fluctuations. It is obvious that the terms containing the χ field are dangerous for the conditions (106). From (104) we find FIG. 6. The 1-loop graph arising from the action (110) that contributes to the tensor spectrum γγ .
Since near the end of the first stage χ 2 ≃ m 2 /g 2 and ω q * ≃ gΦ, one has From this expression it is easy to see that the first condition in (106) is safe. On the other hand, using (104) the second condition in (106) demands It is clear that when the energy density of χ particles catches up the background energy density, i.e. ρ χ ≃ H 2 M 2 p , and this condition is invalidated. This result is independent of our loop considerations and separately indicates the failure of the perturbation theory in this model.
Returning to the interactions involving χ and γ ij , one can use the solutions for the lapse and the shift in (15) to find where a subindex on N or N i indicates that only the relevant terms must be kept in (104) and (105).
To fix the action completely, one should also determine the fourth order terms in K ij K ij − K 2 , however we will not need them for our analysis below. The corresponding corrections can be pictured as in Fig. 6.
It is clear that in (110) A straightforward calculation gives the following one loop correction to the tensor power spectrum: To estimate this correction, we first note that ∂ i χ∂ i χ ≃ a 2 q 2 * χ 2 . We then focus on the last inflaton oscillation cycle in which χ 2 reaches its maximum value. Treating χ 2 as a constant and using (6) for the background quantities one may estimate For our canonical case (7), the integral can be evaluated numerically to yield 0.28. Using (11) and the values of the other dimension-full parameters we obtain We see that this correction is two orders of magnitude smaller than (84). As before, a correction which is supposed to be larger according to the naive counting in perturbation theory turns out to be smaller. Note that both corrections (84) and (114) are still larger than the quantum effects produced during inflation, which are characterized by the ratio H/M p ≃ 10 −8 [6].

V. CONCLUSIONS
In a recent work [18], one of us has shown that the loop quantum effects during reheating significantly modify the scalar power spectrum. In this paper, in an attempt to extend the findings of [18] we consider how loops in reheating produce non-gaussianity and affect the tensor power spectrum in the chaotic m 2 φ 2 model. Based on the tree level results, this model is actually ruled out by 95% confidence level by the Planck data, however our findings show that quantum effects during reheating can change this conclusion since the corresponding corrections can alter the tree level results appreciably.
In most of the scalar field inflationary models, inflation is followed by a period of coherent inflaton oscillations where the background is still homogeneous and isotropic. This phase continues until the backreaction effects are set in. As pointed out in [22], in such a background causality does not preclude the emergence of the superhorizon effects because by coherency the same physical influence can appear at different positions at the same time. Therefore, the quantum effects can be important for cosmological variables in the first stage of reheating. On the other hand, it is known that the entropy perturbations can cause nontrivial superhorizon evolution of the curvature perturbation. Consequently, it is not surprising to see that the effects of entropy modes circulating in the loops become significant, especially in the parametric resonance regime. Indeed, we observe that the corrections get larger as the number of χ modes circulating in the loops increases, which indicates that the perturbation theory might become invalid.
It is well known that in the chaotic model we have studied, the curvature perturbation ζ becomes an ill defined variable during reheating. Because of that reason in [18], the calculations have been carried out in the ζ = 0 gauge till the end of the first stage of reheating and then a gauge transformation has been applied to read the ζζ correlation function. In this paper, we utilize a different strategy and smooth out the spikes of ζ by using the averaged out background variables in the quadratic ζ-action. As it is shown above, the results obtained in this way is consistent with [18] and thus our conclusions about the scalar power spectrum (and non-gaussianity) are firm. Note that the tensor calculation is free from the gauge fixing issues.
It is possible to develop the results of this paper in different directions. Due to the importance of the chaotic m 2 φ 2 model, it would be valuable to perform a full numerical check of the loop corrections that are estimated in this paper. It would also be crucial to see whether the loops in reheating modify the predictions of the models that are favored by Planck data, like the Starobinsky model [39]. Finally, it would be interesting to determine the loop effects when the inflaton decay occurs perturbatively. In that case while the reheating scalar modes cannot take large values, the decay process is completed in a long time that might enhance the quantum effects, since according to in-in formalism (56), the quantum corrections are proportional to the duration of the process.