Curvature perturbations from preheating with scale dependence

We extend the formalism to calculate non-Gaussianity of primordial curvature perturbations produced by preheating in the presence of a light scalar field. The calculation is carried out in the separate universe approximation using the non-perturbative delta N formalism and lattice field theory simulations. Initial conditions for simulations are drawn from a statistical ensemble determined by modes that left the horizon during inflation, with the time-dependence of Hubble rate during inflation taken into account. Our results show that cosmic variance, i.e., the contribution from modes with wavelength longer than the size of the observable universe today, plays a key role in determining the dominant contribution. We illustrate our formalism by applying it to an observationally-viable preheating model motivated by non-minimal coupling to gravity, and study its full parameter dependence.


Introduction
Observations of the cosmic microwave background radiation and large scale structure provide strong evidence for a period of inflation in the early universe, but they do not distinguish well between different specific models of inflation.One reason for this is that in typical models the inflationary era itself is very simple, consisting of slowly rolling scalar fields, and therefore the observational predictions can also be parameterised by a small number of slowroll parameters.Because of this it is interesting to consider the physics of reheating at the end of inflation, which can be very different in different models.
In particular, in many inflationary models reheating involves rapid particle production caused by a parametric resonance between the inflaton and other fields, commonly referred to as preheating.Preheating can be described as the stage after the end of inflation when field inhomogenities grow exponentially leading to large occupation numbers, which then back-react due to non-linear terms in the equation of motion causing the growth to stop at a point where universe enters radiation domination.Large occupation numbers mean this stage can be treated semi-classically and is equivalent to solving the full equations of motion numerically with initial conditions that are obtained from the end of inflation [1,2].
In the presence of light scalar fields, preheating produces curvature perturbations on observable scales [3].Lattice field theory simulations have shown that they can have observable levels of non-Gaussianity, especially in the massless preheating model [3][4][5][6].However, this model is incompatible with tensor-to-scalar ratio measurements from Planck satellite [7].In this paper, we address this issue by considering an observationally viable model of preheating with the inflaton ϕ non-minimally coupled to gravity, that decays into a massless spectator field χ.
The main highlight of our article is that we perform a full calculation by including time-dependence of Hubble rate until the end of inflation in order to obtain initial conditions for numerical simulations.Such a dependence was previously considered in Ref. [3] for the massless preheating model.However, as shown by Ref. [5] the dependence of scalar curvature perturbations on the spectator field value constitutes a spiky pattern, attributable to chaotic behaviour during preheating.This non-linear behaviour requires use of a non-perturbative treatment for obtaining the non-Gaussianity from scalar curvature perturbations [8,9].We find that including the time-dependence of Hubble rate places tight constraints on the "cosmic variance", the values which the mean spectator field can take.This is demonstrated for the preheating model we consider, but would in general be true for any other model of inflation and reheating.A negligible cosmic variance implies the mean spectator field value lies close to zero.As many inflationary potentials are symmetric around zero, it then implies the leading non-perturbative result would be the one at higher order in the non-perturbative delta N formalism.We calculate this result keeping in mind the scale dependence, showing that a simple scale-invariant way of solving momentum integrals [10] is insufficient to yield the correct non-Gaussianity from preheating.
The article is divided as follows: In Section 2, we begin by introducing curvature perturbations arising from preheating and the delta N formalism used to extract non-Gaussianity from them.Then we introduce and extend the non-perturbative delta N formalism to include scale dependence.These are the main results of our article.Next, in Section 3, we build an observationally viable model of preheating motivated by non-minimal coupling to gravity.In Section 4, we consider the time-dependence of Hubble rate to obtain initial conditions for simulations, and find that cosmic variance plays an important role in determining non-Gaussianity.Lastly, in Section 5, we present the f NL calculation from lattice simulations for our preheating model including its full parameter dependence.
2 Curvature perturbations from preheating

Non-Gaussian curvature perturbations
The primordial curvature perturbation ζ is one of the central observables in cosmology.On large scales, it can be measured directly as the relative temperature anisotropy of the cosmic microwave background radiation, and on smaller scales it acts as the seed for structure formation.Assuming that the universe is statistically homogeneous and isotropic, the curvature perturbation can be characterised by its correlation functions in momentum space.The two-point and three-point correlation functions define the curvature power spectrum P ζ and bispectrum B ζ , respectively, ) It is also convenient to define the power spectrum per logarithmic scale This has been measured to be at the comoving scale k * corresponding to physical scale k phys = 0.05 MPc −1 [11].
If the curvature perturbation is exactly Gaussian, the bispectrum B ζ vanishes.Therefore the level of non-Gaussianity can be parameterised by the ratio [12] In general, this is a function of the momenta, but in many cases it is approximately constant.This constant value is referred to as "local" f NL , and it is constrained by the Planck observations as f local NL = −0.9± 5.1 [13].Simplest inflation models based on a single slowly-rolling inflaton field predict very highly Gaussian curvature perturbation, typically f NL < 0.1, which is well below observational bounds.Therefore a detection of a definite non-zero f NL would rule out those models.Conversely, any theory that predicts significant f NL can be tested and constrained by current and future observations.However, it is not necessarily the case that only the inflaton field generates curvature perturbations.They can equally be generated by non-adiabatic fluctuations in another scalar field, which get converted to adiabatic fluctuations during super-horizon evolution [14].Such a field is known as a curvaton.The curvaton can be thought of as taking away the burden from the inflaton to both generate inflation and curvature.This in turn allows for a wide variety of inflaton origin theories such as axion inflation [15], string inflation [16] and others to become viable.
It was noted quite early on [17] that the period when energy from inflaton is transferred to Standard Model particles known as reheating [18] occurs in a non-linear manner.Massless preheating [19] was the toy model with potential, used to study reheating between inflaton ϕ and a scalar field χ [20].Here, the second field χ can also be thought of as a massless curvaton.During preheating, energy is rapidly transferred from the inflaton to the spectator field via parametric resonance.
In terms of perturbation theory, both fields can be split into a background part and a fluctuation.During inflation the background χ field is negligible.Inflation ends when the background ϕ field reaches near its minimum and slow-roll condition fails.However, due to the build up of kinetic energy, it overshoots the minimum and begins oscillating around it.The linear order equations of motion for the fluctuations of both fields now have an oscillating potential arising from the background ϕ.It is a generic feature of second order differential equations with an oscillating term that their solutions can develop an instability for particular parameter values, and grow exponentially.For the massless preheating case, parametric resonance occurs between g 2 /λ = 1 and g 2 /λ = 3 for the zero momentum mode [19].
In reality, fluctuations cannot grow indefinitely and stop when non-linear terms in their equations of motion start to become significant.After a sufficient time has passed, the fluctuations die down and the background inflaton settles into its minimum.With massless preheating (2.6), since both fields are massless, the universe then exits into a radiation dominated equation of state.The entire evolution from end of inflation to onset of radiation domination can be captured fully by using lattice simulations [3][4][5] which predict large non-Gaussianity for the massless preheating model.However, this model is not realistic as it predicts a tensor-to-scalar ratio that is too high [6] when compared to current observations [7].

Delta N formalism
In typical inflationary models, the primordial curvature perturbation arises predominantly from the inflaton field and is nearly Gaussian, but as was shown in Refs.[3][4][5], preheating can produce a potentially observable non-Gaussian contribution.
Because we are interested in the curvature perturbations on very large, superhorizon scales, we can use the separate universe approximation to describe its evolution.This approximation is based on the observation that if two spatial volumes are sufficiently far apart, so that light cannot travel from one to the other during the time interval we are interested in, they evolve independently of each other.In the calculation, they can therefore be thought of as two separate universes.In practice, the size of these separate spatial volumes can be taken to be of the order of the Hubble length, and therefore we will refer to them as Hubble volumes.
Furthermore, it is often the case that the universe can be assumed to be nearly homogeneous and isotropic on the scale of one Hubble length, and in that case we can approximate each separate universe with the Friedmann-Robertson-Walker (FRW) metric, (2.7) In this approximation, each Hubble volume therefore has its own scale factor a(t) which evolves according to the Friedmann equation where ρ is the local energy density in that Hubble volume.The curvature perturbation ζ can be defined as the logarithm of the scale factor on a uniform energy density slice.For practical calculations it is convenient to define an initial flat slice with a uniform scale factor a ini .The curvature perturbation on a slice with uniform energy density ρ ref is then given by the number of e-foldings between the two slices [21], and N is the mean value of N in the observable universe.
If the equation of state is the same everywhere, every Hubble volume expands in the same way, and therefore no curvature perturbation is generated.However, if there is a scalar field χ that has a different local value χ(⃗ x) in each Hubble volume, it can affect the amount of expansion, generating a position-dependent curvature perturbation [22], (2.10) In this paper, we assume that there are two scalar fields: the inflaton ϕ which dominates the energy density during inflation and satisfies the slow roll conditions, and another field χ, which is light compared to the Hubble rate during inflation so that it has significant superhorizon correlations, and which interacts sufficiently weakly with the inflaton so that its effect can be neglected in the early stages of inflation.We also assume that its selfinteractions are weak so that it can be assumed to be Gaussian.
The total curvature perturbation ζ = ζ ϕ + ζ χ also includes a contribution from the inflaton field ϕ.This additive form is valid if during inflation, the inflaton field dominates the expansion and is slowly rolling and that the back-reaction from χ is negligible, as we assume.Equation (2.10) shows that to obtain the statistics of the curvature perturbation generated by the field χ, one needs to determine δN (χ) which is a function of the local value of the field χ, and the statistics of the field χ.
Assuming statistical homogeneity and isotropy, the statistics of χ are fully described by the two-point correlation function or its Fourier transform Σ(k) which satisfies where δχ = χ − χ, and χ is mean value of χ in the observable universe.In analogy with Eq. ( 2.3) we also define If δN is not a linear function of χ, then it gives rise to a non-Gaussian contribution to the curvature perturbation.In many simple scenarios, δN can be well approximated by a quadratic Taylor expansion, Assuming that the dominant contribution to the curvature perturbation is given by ζ ϕ and that it is Gaussian, the leading term in the bispectrum is [10] Considering an equilateral configuration, However, if δN ′ is sufficiently small, which happens if the theory has symmetry under the sign change χ → −χ and χ is small, the dominant term is given by whose Fourier transform gives In this case f NL is therefore given by [10] However, numerical simulations [3][4][5][6]23] have demonstrated that in preheating scenarios, δN can be a very complicated function and therefore this simple Taylor expansion is not applicable.Hence an improvement to this formalism in the form of non-perturbative delta N formalism [8,9] is needed.

Non-perturbative delta N formalism
Non-perturbative delta N formalism [8,9] performs expansion of the curvature correlator in terms of field covariance Σ(x) instead of expanding δN in powers of the field δχ as is done in perturbative delta N approach.In general, correlation functions of the curvature perturbation ζ χ generated by the field χ are given by the joint probability distribution p(χ 1 , . . ., χ n ) of field values Because we are assuming that the field χ is a Gaussian random field, its probability distribution can be expressed in terms of its two-point correlator Σ(x) and its Gaussian one-point probability distribution p(χ), which is determined by its mean χ and variance ⟨δχ 2 ⟩ as For example [9], the two point correlator can be expanded as where

and the coefficients are given by
(2.23) The non-perturbative expansion is valid if Σ ij << ⟨δχ 2 ⟩.
To calculate f NL , we also need the three-point correlator which can be expanded in terms of the covariance as where "perms" indicates the sum of different permutations of 1, 2 and 3. To obtain the power spectrum and bispectrum, we need to respectively Fourier transform the two-point and three-point coordinate space correlators above.

Scale dependence
In Ref. [9] the non-perturbative delta N expansion Eq. (2.22) and Eq.(2.24) were truncated at first order in covariance to yield power spectrum and the bispectrum which can be seen to be analogous to Eq. (2.15), with the substitutions δN ′ → Ñχ and δN ′′ → Ñχχ .Therefore the value of the non-Gaussianity parameter f NL also has the same form, (2.27) However, as in the perturbative case, this is not necessarily the dominant term, especially if Ñχ is small.Therefore we go to the next order of the expansion in powers of Σ(k), Considering the case Ñχ → 0, we therefore obtain for the equilateral case which is again analogous to the earlier result (2.19) but is valid even when the Taylor expansion (2.14) is not.In Ref. [10], the authors assumed a scale-invariant power spectrum for the field χ, where H * is the Hubble rate at the time when the current observational scale left the horizon.With that assumption, the non-Gaussianity parameter f NL given by Eq. (2.30) at the observational scale k * becomes approximately with an infrared cut-off L −1 corresponding to the size of the observable universe.Therefore its k-dependence appears logarithmic.However, in practice the field χ will not be exactly scale invariant because of its mass and the time-dependence of the Hubble rate during inflation.To capture the effect of the scale dependence, we approximate the field variance with a power law, where the spectral index n χ (> 0) and the amplitude A χ , which has energy dimensions 2−n χ , depend on the parameters of the specific inflation model as we will discuss in more detail in Section 5.2.This form implies the real space correlator behaves on asymptotically long distances as Σ ij ∼ |⃗ x i − ⃗ x j | −nχ and is therefore small on observable scales, justifying the non-perturbative delta N expansion.Then the integral in Eq. (2.30) gives for the equilateral case, This shows that f NL has power-law dependence on the observational scale k * .Because cosmological observations are carried out on comoving scales that are many orders of magnitude larger than the Hubble scale during inflation, this leads to a significant suppression of f NL if n χ > 0.

Potential
Calculations of non-Gaussianity in the massless preheating model have been made by simulating full non-linear behaviour in the perturbative delta N formalism [3].Similar calculations using the non-perturbative delta N formalism [6] yield significant non-Gaussianity f NL but the tensor to scalar ratio is very high for inflaton power spectrum within observational bounds.It has been shown that introduction of non-minimal coupling of inflaton to gravity is able to satisfy observational bounds on spectral index and tensor to scalar ratio [7,24].In order to maintain the theoretical relaxations made possible by the curvaton scenario described above, along with satisfying non-Gaussianity, spectral index and tensor to scalar ratio bounds, we need to study massless preheating with non-minimal coupling to gravity.
Even if no non-minimal coupling is included in the classical action, one loop calculations made using QFT in curved spacetime automatically generate the non-minimal coupling [25].Therefore, it is necessary to include the non-minimal coupling.Indeed, a minimal coupling would require there to be some extra symmetry which forbids the non-minimal coupling term from appearing in the action.The coupling takes a scale invariant form: ξϕ 2 R/2 where ϕ is the inflaton field, R is the scalar curvature and ξ is the non-minimal coupling parameter which is dimensionless.The QFT calculations further show that the non-minimal coupling parameter runs logarithmically with energy scale but does not possess any fixed point [26].Thus, there is no a priori constraint on the value of the non-minimal coupling parameter.In this work, we will assume that it has a small positive value 0 < ξ ≪ 1/6.
The Jordan frame action for massless preheating with only the inflaton field coupled non-minimally to gravity is taken to be where This Jordan frame action can be converted to the Einstein frame by conformally rescaling the metric to remove the non-minimal coupling and then redefining the fields to bring their kinetic terms back to the canonical form.
In our case, using a conformal rescaling gµν = (1 + ξ ϕ 2 where and G ij is the metric in field space with components We can check that the Ricci scalar of this field space metric is not zero at order ξ.Hence even with only the inflaton non-minimally coupled to gravity it is not possible to bring back canonical kinetic terms for both fields as the field space metric is not conformally flat.
Furthermore it has been shown [27] that bringing back canonical kinetic terms is generically not possible if more than one field is non-minimally coupled.However if only a single field ϕ was present in the Jordan frame action, then the field space metric has only one component G ϕϕ and a redefined field variable φ which satisfies the differential equation where the final form is valid for ξ ≪ 1/6, brings back the canonical kinetic term.This gives us the relation and the Einstein frame potential Such a potential is also motivated by α−attractor T models of inflation [28].Therefore, in order to simplify numerical simulations by not having to include noncanonical kinetic terms, we choose a potential for the two-field case as (3.9) A similar approach to preheating simulations in α-attractor motivated potentials has been used in Ref. [29].Here, we have kept the identification ϕ → (M P / √ ξ) tanh( √ ξϕ/M P ) in the interaction between the two fields to preserve the flatness of the potential at large inflaton values, which maintains the desirable features of α-attractor models.Henceforth we shall drop the tilde on the fields for simplicity, bearing in mind that they belong to the Einstein frame.

Parameter values
Our non-minimally coupled (NMC) preheating model has three parameters: λ, g 2 , and ξ.Of these, we shall see that λ is constrained by the power spectrum observation and ξ is constrained by tensor-to-scalar ratio observation.
During inflation we assume the field χ is negligible, so we have effectively a single-field inflationary model with the inflaton potential The first slow-roll parameter is Setting ϵ = 1, we find that inflation ends when ϕ reaches the value ϕ end given by cosh We normalise the scale factor a to be unity at this time, a end = 1.
Using the slow-roll field and Friedmann equations, we can then relate the number of e-foldings until the end of inflation N = − ln a to the inflaton field value, and the Hubble rate to N as The curvature power spectrum due to the inflaton field is therefore [30] where we have used Eqs.(3.11) and (3.14), and ϕ * is the field value N * e-foldings before the end of inflation, when the observational scale k * exited the horizon, The precise value of N * depends on post-inflationary evolution of the universe, but to be specific we will assume N * = 55.Setting P * equal to its observed value, P * = 2.1 × 10 −9 , we obtain As we can see, λ is fully determined once we fix ξ.
Parameter ξ is constrained from below by the tensor to scalar ratio observation [7].At the observational scale N * , our model gives [30]

Initial Conditions
After introducing our model in the previous section, the next step is to calculate non-Gaussianity parameter f NL given by Eq. (2.30).
In order to apply the non-perturbative delta N formalism, we need to know the statistics of the χ field at the initial time slice a ini defined in Eq. (2.9).Because the field is weakly coupled, its evolution during inflation can be described using linear theory in momentum space.It is convenient to divide the comoving momenta k in three separate ranges: Modes with k < a 0 H 0 , where a 0 and H 0 are the scale factor and Hubble rate today, have wavelength longer than the size of the currently observable universe.Therefore from the perspective of any observation, they appear as a uniform, non-zero mean value, which we denote as χ.Because of this random origin, the precise value of χ cannot be computed, but if we know the whole inflationary history, we can compute its variance ⟨χ 2 ⟩, which is often referred to as the "cosmic variance".
Modes with a 0 H 0 < k < a ini H ini , where H ini is the Hubble rate at the initial slice, have a wavelength that is shorter than the size of the observable universe today but longer than the size of a single "separate universe".Therefore these modes give rise to a different initial value χ ini in each separate universe and are responsible for the produced curvature perturbation.More precisely, in the context of the non-perturbative delta N calculation, they determine the field correlation function Σ(x) and the one-point probability distribution (2.21), which depends on the mean χ and the variance ⟨δχ 2 ⟩.
Modes with k > a ini H ini , where H ini is the Hubble rate at the initial slice, have a wavelength shorter than the size of a single "separate universe".Therefore they are statistically homogeneous on the observational scales and do not contribute directly to the curvature perturbation.However, they play an important role in the calculation and are discussed in Section 5.1.
At linear order, a mode χ k with comoving momentum k satisfies the mode equation where the effective ϕ-dependent mass of the χ field is we can see that assuming g 2 ξ/λ < 1, the χ field is light early on during inflation, and its mass increases gradually relative to the Hubble rate during inflation.For convenience, we choose the initial time slice a ini in Eq. (2.9) to correspond to the moment when which gives Because the field is light, we can approximate that when the comoving mode χ k exists the horizon at time t k , its amplitude is Afterwards, its evolution is overdamped and its amplitude decays with rate At the time slice t ini , the amplitude is therefore where N k and N ini denote the number of e-foldings until the end of inflation at times t k and t ini , respectively, and which for N ini in Eq. (4.5) is given by This implies the power spectrum (4.8) is scale dependent through the function F (N k ).Therefore we need to go beyond the scale-invariant approximation used in other works Refs.[6,10].
In order to compute f NL using Eq.(2.30), we need the two-point correlator Σ(k), which is given by Eq. (2.13), as well as the one-point probability distribution p(χ), which is given by Eq. (2.21) and depends on the variance ⟨δχ 2 ⟩ and the mean value χ.As discussed earlier, the latter is a random variable whose variance ⟨χ 2 ⟩ we can compute.
The two variances are both given by a similar integral, with a different integration range, ) where the latter forms are obtained by changing the integration variable to N k , the value of N at the time when mode k left the horizon, with H(N k ) given by Eq. (3.14).
Throughout this article we consistently choose the number of efolds at which the currently observable scale left the horizon to be N * = 55.Cosmic variance is two orders of magnitude less than the variance of the initial χ field in the range 1 < g 2 /λ < 3 which is relevant for parametric resonance [19].For g 2 /λ > 3, it is more than four orders less.It only becomes comparable to the initial variance at g 2 /λ < 0.5 which is well outside the parametric resonance band.
The cosmic variance integral Eq. ( 4.11) at large N k goes as From which we may deduce that for ξ = 0 the integral has no exponential suppression and hence cosmic variance for the minimal coupling case is infinite [4].However, the introduction of even a small positive non-minimal coupling ξ > 0 gives rise to an exponential suppression and hence we expect cosmic variance to be negligible for the non-minimal coupling case.This is indeed what we find by numerically integrating for the cosmic and initial variances in Eqs.(4.11) and (4.12), as summarised in Figure 1.
Whereas an infinite cosmic variance would imply that the mean χ that appears in the probability density p(χ) can have any value and can therefore be considered to be a free parameter, a negligible cosmic variance forces us to choose χ ≈ 0. Symmetry of the potential around zero dictates Ñχ will be small which, as we saw in Section 2.3, makes the dominant term in the non-perturbative delta N formalism to be Eq.(2.30).Thus we illustrate the important role played by cosmic variance in calculating the non-Gaussianity from preheating.If cosmic variance is not negligible, as happens for the massless preheating toy model, then one needs to use Eq.(2.27) as done in Ref. [6].However if it is negligible, as happens in our NMC preheating model, then we need to use Eq.(2.30) to calculate f NL .The only qualification is for the potential to be symmetric in χ, which is true in most preheating scenarios.

Lattice simulations
We use the program HLattice [31] to perform fully non-linear simulations of preheating dynamics.For a specific value of g 2 /λ and ξ, we draw N runs = 10, 000 initial field values χ i , i ∈ {1, . . ., N runs }, from the probability distribution (2.21) with zero mean, i.e., χ = 0, and with ⟨δχ 2 ⟩ given by Eq. (4.12).A lattice simulation is run and the value of δN is extracted for each of the N runs values separately.This amounts to importance sampling the probability distribution, and therefore the coefficients Ñχ and Ñχχ defined in Eqs.(2.23) can be obtained as averages over the runs, The fields are evolved in a linear fashion till the condition that the potential energy is equal to the kinetic energy is met.Thereafter non-linear simulations begin.Using the principles first outlined in [2], we can study the semi-classical dynamics occurring during preheating by classically simulating an equivalent non-equilibrium, non-linear problem.In HLattice, this is achieved through capturing the quantum nature by initialising the fields with random Gaussian fluctuations for inhomogeneous modes.
Full equations of motion for the fields and Friedmann equation are solved on a comoving lattice with 32 3 grid points and length of a side equal to 20/H, where H is the Hubble rate at the start of the non-linear evolution.The lattice size has been chosen as a compromise between two factors: Ideally, we would want the lattice to be large enough so that light cannot travel across the lattice during the nonlinear non-equilibrium stage of the evolution, but if it is larger than the comoving Hubble length, then we would need to include metric perturbations and the use of the FRW metric would not be justified.For our lattice size, the comoving Hubble length is initially somewhat smaller than the lattice size but always lies above the lattice spacing.We choose the discretisation scheme LATTICEEASY for evolution in HLattice as we are not including any metric perturbations.
We have also checked that changing the lattice spacing from 32 3 to 16 3 grid points while keeping lattice size 20/H fixed as well as changing the lattice spacing from 32 3 to 64 3 grid points while changing lattice size from 20/H to 40/H does not alter our simulation results significantly.Furthermore, the comoving mass scale (am χ ) −1 where m 2 χ = g 2 ⟨tanh 2 √ ξϕ⟩/ξ that is a representative of the field dynamics, lies between the lattice spacing and lattice size throughout the non-linear evolution.
Plotting the raw a versus H data obtained from HLattice in the form of deviation from radiation domination, Figure 2 shows that all curves become flat at late times, indicating onset of radiation domination.It also shows the change in asymptote for two different values of χ ini .Using Gaussian filtering in conformal time with standard deviation of 100 time steps, we can precisely pick the value of N = ln(a) at H ref = 5 × 10 −15 M P which would correspond to a constant energy density ρ ref to yield the N function Eq. (2.9).The two curves deviate due to field excursion in transverse direction [5], giving rise to a spike in the delta N function at χ ini = 1.581 × 10 −6 M P as compared to χ ini = 0.
Figure 3 compares the δN functions obtained for three different values of the nonminimal coupling ξ.Particularly relevant are the values ξ = 0 which corresponds to the minimal coupling case and ξ = 0.004 which is the lower bound from observations Eq. (3.19).The same inhomogeneous modes are used to initialise simulation for all χ ini .We expect the uncertainty arising from different random realisations of the inhomogeneous modes to be largely uncorrelated between different Hubble volumes, making its effect subdominant to   1.The plot is not exactly symmetric under χ → −χ due to same random realisation of inhomogeneous modes being used for each run, which breaks the symmetry.the statistical error in sampling of χ ini .As a quantitative check we carry out 32 runs with different random realisations of the inhomogeneous modes for g 2 /λ = 2 and ξ = 0.004 at χ ini ≈ 1.581 × 10 −6 M P , 1.203 × 10 −6 M P and 0M P .We find the standard deviation in δN due to this effect is approximately 40% of the spike heights in Figure 3. Therefore the spikes are a genuine feature of our simulations and this effect is subdominant.
Even though the cosmic variance ⟨χ 2 ⟩ is small, it is not exactly zero, and therefore it is important to consider the full range of possible values of χ.Because the relevant values are  1.Numerical values of quantities to be specified in simulations for different ξ but all with the same g 2 /λ = 2. small, we can achieve this without carrying out new simulations by re-weighting the data: For a general small value of χ, the coefficients (2.23) are given by Ñχ = 1 Figure 4 shows the non-perturbative coefficient Ñχχ over the range of χ allowed by the cosmic variance.The final value of Ñχχ and error bar for each χ are obtained by using the bootstrap (resampling with replacement) method.We use the symmetry χ → −χ to double the number of sample points δN (−χ) used in calculating Ñχχ .

Calculating f NL
In section 2.4, we described how a scale dependent way of solving momentum integrals can significantly change the non-Gaussianity calculated from preheating.Here, we illustrate it using our NMC preheating model.
From lattice simulations, Ñχχ ∼ 10 6 M −2 P for g 2 /λ = 2 and ξ = 0.004.Putting these values together gives in the scale invariant approximation.However, the correct f NL shall be obtained by taking a scale-dependent estimate of the field variance.For our model of NMC preheating, this essentially implies including the over-damped oscillator contribution F (N ) from Eq. (4.8) in the variance.To obtain N in terms of k, we need to invert the horizon crossing relation k = a(N )H(N ), This is difficult to invert analytically.However, for large N we notice that the expression is dominated by the exponential factor.Therefore in order to have an analytical estimate, we use the approximation H = const that allows us to capture the large N behaviour.Fixing H from Eq. (3.14) to H * λ 12 and identifying it in the scale-dependent power spectrum at initial time slice (4.8) with the exponential damping coefficient (4.10) becoming the correlator takes the form where (5.10) This form separates the momentum q dependence explicitly and is precisely the form (2.33) we used in Section 2.4 to calculate an analytical estimate (2.34) for f NL .
For typical value of parameters g 2 /λ = 2 and ξ = 0.004 in our preheating model, we have n χ ∼ 0.1, A χ ∼ 10 −12 M 2−nχ P .Ñχχ ∼ 10 6 M −2 P for g 2 /λ = 2 and ξ = 0.004 from numerical simulations remains the same.However there is a small subtlety in specifying the value of comoving scale k * .We have fixed scale factor to be one at the end of inflation, a end = 1.However the physical scale k phys = 0.05MPc −1 at which Planck satellite measurements take place assumes scale factor is one at present times, a 0 = 1.Therefore k * = (a 0 /a end )k phys .Assuming instantaneous reheating [32] The order of magnitude estimate for f NL can be made more precise by calculating the momentum integrals involved in Eq. (2.30) numerically.We split the domain of integration for the momentum q into q < k * and q > k * .For q < k * , the constant H * approximation is valid and momentum integral can be obtained analytically as in Eq. (2.34).For q > k * , we invert q(N ) numerically to give Σ num (q) and momentum integral as (5.14) Combining the integrals for the above two domains yields the full integral from q = 0 to q = q(N ini ).The remaining ingredient to be used in Eq. (2.30) is the non-perturbative coefficient Ñχχ .From the lattice simulation results at different χ, figure 4, we take the average Ñχχ and average error ∆ Ñχχ to estimate f NL and its error through standard propagation.For g 2 /λ = 2, ξ = 0.004 we find which matches well our order of magnitude estimate above.We can see that there is a significant decrease in the non-Gaussianity estimate if we perform a complete linearised calculation, rather than using the scale invariant approximation (2.31).Essentially, this occurs because the momentum integral in Eq. (2.30) is power law suppressed in k * ∼ 10 −30 M P which results in a very low overall factor in f NL .Thus we illustrate that for non-Gaussianity calculation in any model of preheating, it is important to include scale-dependence during inflation to yield an accurate result.

Parameter dependence
The ξ dependence of momentum integral occurs through the factor A 3 χ k 3nχ * /n χ in Eq. (2.34).Figure 5 shows that this factor drops down almost exponentially as ξ is increased from the lowest value of ξ = 0.004 allowed in our model.The non-perturbative coefficient Ñχχ Eq. (2.23) captures the departure from a constant δN within a spread of χ ini governed by the initial variance.Now as can be seen from Figure 3, the height of spikes reduces as ξ increases.Furthermore, the number of spikes decreases as the spread in χ ini shrinks.Figure 5 also shows that the initial variance does not change significantly with ξ when compared to the change in momentum integral.Combining both the factors, we surmise f NL will be maximum for lowest observational bound on ξ, ξ = 0.004 for a given value of g 2 /λ.
We can therefore search for the highest value of f NL in our model by keeping ξ = 0.004 fixed and scanning over different g 2 /λ. Figure 6 shows that the maximum value of f NL occurs around g 2 /λ = 1.625.At exactly g 2 /λ = 1.625 we find it to be f NL = −(6.3± 0.7) × 10 −4 . (5.16) This value still lies well-within the current observation bound of f local NL = −0.9± 5.1 [13], but it is several orders of magnitude larger than what we obtain for generic g 2 /λ, e.g., Eq. (5.15).
Based on a linearised analysis of the parametric resonance for the massless preheating case, one might expect the maximum value of f NL to lie at g 2 /λ = 1.875 [4].However, our simulations indicate this is shifted to g 2 /λ ≈ 1.625.
To investigate this shift we can compare the growth of zero mode with root mean square of the inhomogeneous fluctuations in Figure 7.We see that the zero mode for g 2 /λ = 1.875 grows with a slightly higher slope than g 2 /λ = 1.625 in accordance with parametric resonance considerations.However the major difference in evolution for the two coupling strengths arises from the non-linear behaviour of fluctuations.Figure 7 shows that after the evolution becomes non-linear, the fluctuations catch up significantly slower for the g 2 /λ = 1.625 value, thus allowing for exponential growth of zero mode to persist for a longer time and hence increasing the δN (χ) that increases f NL as compared to g 2 /λ = 1.875 value.This effect can only be captured by lattice simulations.Thus by finding many orders of magnitude change in f NL with the coupling strength, we have also demonstrated the need for full lattice simulations to study parameter dependence of preheating models.

Conclusion and Discussion
We have extended the non-perturbative formalism [3,6,9] used to calculate non-Gaussianity generated during preheating and illustrated it by applying to an observationally viable model  of massless preheating motivated by non-minimal coupling to gravity.
Interestingly, we find that calculating non-Gaussianity from preheating at currently observed scale, which occurs after the end of inflation, involves taking into account the evolution of modes at very early times, well-before they left the horizon.
In a separate universe picture, the evolution of each individual causal volume of our currently observable universe depends sensitively on initial conditions [5].However, the overall probability distribution of initial conditions is Gaussian and depends on the initial variance as well as the mean over our currently observable universe.This mean itself depends on the cosmic variance over the entire universe, of which our observable universe is just a small part.
If the spectrum of the spectator field is sufficiently blue-tilted during inflation, then the cosmic variance becomes negligible and hence the dependence on modes before the currently observed scale left the horizon disappears.As an example we have shown how this occurs in our preheating model.We are able to arrive at this finding precisely because we have considered the full time-dependence of Hubble rate during inflation.
When time-dependence of Hubble rate is included, it translates to the spectator field having a scale-dependent spectrum.Thus we may no longer use a simple scale-invariant approximation to calculate non-Gaussianity [10].We have also shown how considering scale dependence yields a f NL value that is many orders of magnitude different from just using a scale-invariant approximation.
As noted in the text, we drop the effects of non-canonical kinetic terms arising from conversion of Jordan to Einstein frame in the two-field case to simplify the simulations.Simulations with non-canonical kinetic terms have been performed [33] for non-minimal coupling parameter ξ > 1.Such simulations, which incorporate the curvature in field space, are able to capture the full effects of non-minimal coupling and can be used to calculate non-Gaussianity using the formalism outlined in this article.Lattice simulations of preheating with non-minimal coupling can also be performed directly in the Jordan frame [34].It would be interesting to perform these simulations for the non-minimally coupled massless preheating model we have considered and compare our results.
We have scanned for f NL over a parameter range g 2 /λ = 1 to 3 in our model, which is relevant from parametric resonance considerations [19].We find that f NL changes by many orders of magnitude within this range.Therefore both detection or non-detection of non-Gaussianity can be used to constrain the coupling strength in preheating models.The maximal f NL ∼ 10 −4 occurs around g 2 /λ = 1.625 with ξ = 0.004 fixed.This value of the non-Gaussianity parameter is consistent with current observational bounds.It is too small to be detected by cosmological experiments in the near future, but perhaps is large enough that a different measure of non-Gaussianity, for example using wavelet transform [35], can be used to detect it.Finding such a measure that captures the inherently non-linear nature of the fields during preheating is a possible direction for future work.
Nonetheless, the full calculation of f NL in our model enables a complete illustration of our formalism to obtain non-Gaussianity from preheating.This method can very well be applied to other interesting models of inflation and reheating, possibly yielding high, detectable non-Gaussianity.

Figure 1 .
Figure 1.Ratio of the cosmic variance to the initial variance for g 2 /λ ∈ [0.5, 3] at fixed ξ = 0.004.Cosmic variance is two orders of magnitude less than the variance of the initial χ field in the range 1 < g 2 /λ < 3 which is relevant for parametric resonance[19].For g 2 /λ > 3, it is more than four orders less.It only becomes comparable to the initial variance at g 2 /λ < 0.5 which is well outside the parametric resonance band.

Figure 2 .
Figure 2. The deviation from radiation domination a ∝ 1/ √ H at two χ ini values for g 2 /λ = 2 and ξ = 0.004.χ ini = 1.581 × 10 −6 M P corresponds to spike in the delta N function.Inset shows zoomed in data around H ref = 5 × 10 −15 M P where we pick the value of N = ln(a) after filtering.

Figure 3 .
Figure 3. δN at H ref = 5 × 10 −15M P for ξ = 0, 0.004, 0.007 and g 2 /λ = 2.The relevant initial quantities used to plot δN are summarised in Table1.The plot is not exactly symmetric under χ → −χ due to same random realisation of inhomogeneous modes being used for each run, which breaks the symmetry.

Figure 4 .
Figure 4. Ñχχ for 1000 χ drawn from a Gaussian distribution with mean equal to zero and variance equal to cosmic variance ⟨χ 2 ⟩.The model parameters used are g 2 /λ = 2 and ξ = 0.004.The final value and error bar for each χ are obtained by bootstrapping 100 resamples.