Non-Gaussianity from resonant curvaton decay
We calculate curvature perturbations in the scenario in which the curvaton field decays into another scalar field via parametric resonance. As a result of a nonlinear stage at the end of the resonance, standard perturbative calculation techniques fail in this case. Instead, we use lattice field theory simulations and the separate universe approximation to calculate the curvature perturbation as a nonlinear function of the curvaton field. For the parameters tested, the generated perturbations are highly non-Gaussian and not well approximated by the usual parameterisation. Resonant decay plays an important role in the curvaton scenario and can have a substantial effect on the resulting perturbations.
Observations of the cosmic microwave background radiation are consistent with Gaussian perturbations , but there are tantalising hints of non-Gaussianity  at a level that would be clearly observable with the Planck satellite and other future experiments . This would rule out the simplest inflationary models with slowly rolling scalar fields, which can only produce very-nearly-Gaussian perturbations [4, 5, 6, 7].
Models which can generate large non-Gaussianities either during, or at the end of, inflation include those with non-canonical Lagrangians [8, 9, 10, 11, 12, 13, 14, 15, 16], those where there is breakdown of slow-roll dynamics during inflation due to sharp features in the potential [17, 18], those with multiple fields with specific inflationary trajectories in the field space [19, 20, 21, 22, 23, 24, 25, 26, 27] and those where additional light scalars affect the dynamics of perturbative or non-perturbative inflaton decay [28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39].
A well-known example of a multi-field model where large non-Gaussianities can be generated after the end of inflation is the curvaton scenario [40, 41, 42, 43, 44]. In this model, the primordial perturbations arise from the curvaton field, a scalar field which is light relative to the Hubble rate during inflation but, unlike the inflaton, gives a subdominant contribution to the energy density.
After inflation the relative curvaton contribution to the total energy density increases. This affects the expansion of space, and the curvaton perturbations become imprinted on the metric fluctuations. The standard adiabatic hot big bang era is recovered when the curvaton eventually decays and thermalises with the existing radiation. The mechanism can be seen as a conversion of initial isocurvature perturbations into adiabatic curvature perturbations during the post-inflationary epoch. If the curvaton remains subdominant even at the time of its decay, its perturbations need to be relatively large to yield the observed amplitude of primordial perturbations. Consequently, in this limit the curvaton scenario typically generates significant non-Gaussianities [45, 46, 47, 48, 49, 50, 51, 52].
Almost all the analysis of the curvaton scenario is based on the assumption of a perturbative curvaton decay. As pointed out in refs. [53, 54], it is, however, possible that the curvaton decays through a non-perturbative process analogous to inflationary preheating [55, 56, 57]. This is a natural outcome in models where the curvaton is coupled to other scalar fields which acquire effective masses proportional to the value of the curvaton field.
During this short period of non-equilibrium physics, part of the curvaton condensate decays rapidly into quanta of other light scalar fields. In analogy to inflationary preheating, the decay is not complete and the remaining curvaton particles have to decay perturbatively. The properties of primordial perturbations generated in the curvaton scenario depend sensitively on the dynamics after the end of inflation until the time of the curvaton decay [58, 59, 60]. As the preheating dynamics are highly nonlinear, it is natural to ask if this stage could generate large non-Gaussianities.
Our aim in this work is to calculate the amplitude and non-Gaussianity of perturbations from curvaton preheating. To do this we use the separate universe approximation [61, 62, 63, 64] and classical lattice field theory simulations [65, 66]. This method was applied to inflationary preheating in refs. [37, 38, 39]. We solve coupled field and Friedmann equations to determine the expansion of the universe during the non-equilibrium dynamics and the fraction of the curvaton energy which is converted into radiation. From these, we calculate the curvature perturbation as a function of the local value of the curvaton field. We find that, at least for our choice of parameters, this dependence is highly nonlinear, implying very high levels of non-Gaussianity.
The structure of the paper is as follows. In section 2 we review the curvaton model and derive useful general results. In section 3 we present an analytic calculation of production of curvature perturbations using linearised approximation of preheating dynamics, and demonstrate that the answers depend on quantities that are not calculable in linear theory. A full nonlinear computation is, therefore, needed. In section 4 we show how this can be done, and compute the amplitude and non-Gaussianity of the perturbations for three sets of parameters. Finally, we present our conclusions in section 5.
2 General theory
We study the curvature perturbations generated in a simple curvaton model with canonical kinetic terms and the potential
where is the inflaton, is the curvaton and is a scalar field which the curvaton will decay into. During inflation the inflaton potential dominates the energy density. We assume standard slow roll inflation, so the inflaton field rolls down its potential until it reaches a critical value at which the slow roll conditions fail and inflation ends. As usual, the perturbations of the inflaton field generate curvature perturbations, but we assume that their amplitude is negligible. This requires , where is the Hubble rate and is the slow roll parameter, both of which are determined by the inflaton potential .
Instead, we assume that the observed primordial perturbations are generated solely by the curvaton field , which should be light during inflation so that it develops nearly scale-invariant quantum fluctuations similar to those of the inflaton field. The curvaton field remains nearly frozen at the value it had at the end of inflation, until the time at which the Hubble parameter has decreased to and the curvaton starts to oscillate around the minimum of its potential. These fluctuations contribute to the energy density, and, therefore, they affect the expansion of the universe and generate curvature perturbations.
Eventually, the curvaton decays into lighter degrees of freedom. In the standard curvaton scenario this decay is assumed to take place perturbatively. The interaction term in (1) only allows pair annihilations, the rate of which falls quickly when the universe expands. Some of the curvaton particles would, therefore, survive until today, behaving as dark matter. While this might be an interesting scenario to study, we follow most of the literature and assume that the curvaton is also coupled to other fields. In particular, a Yukawa coupling to a light (possibly Standard Model) fermion field of the form
would give the decay rate 
Without further knowledge of the curvaton couplings, the decay rate appears as a free phenomenological parameter, which is bounded from below because the curvaton cannot decay arbitrarily late without spoiling the success of the hot big bang cosmology. The most conservative lower bound for the decay time is given by nucleosynthesis: .
In this work, we assume that the other scalar field in (1) is heavy during inflation, , so that its value remains close to zero and it develops no long-range perturbations. It was recently shown in ref.  that in this case some of the curvaton particles can decay via a parametric resonance into quanta of the field. The resonance does not destroy all the curvaton particles, so a Yukawa coupling is still needed to allow those which remain to decay perturbatively. However, it reduces their density significantly.
2.2 Curvature perturbations
Using the approach, the curvature perturbation on superhorizon scales is given by the expression 
where is the scale factor at some fixed energy density , and denotes the difference from the mean value over the whole currently observable universe. The scale factor is normalised to be constant at some earlier flat hypersurface, which we choose to be at the end of inflation. This measures the difference in the integrated expansion between Friedmann–Robertson–Walker (FRW) solutions with different initial conditions which are evolved until the same final energy density . In our case, the variations of initial conditions arise from the super-horizon fluctuations of the curvaton field produced during inflation. This means that the curvature perturbation is a local function of the curvaton field fluctuations . As we know the statistics of , this determines the statistics of the curvature perturbation completely. In this paper we compute this function.
Observations show that the primordial perturbations were nearly Gaussian, which implies that should be approximately linear. Therefore, it is natural to Taylor expand (4) as
where the prime denotes a partial derivative with respect to the curvaton value at the end of inflation evaluated at constant final energy density . This implies that, to leading order in , the power spectrum of the curvature perturbations is
Neglecting the inherent non-Gaussianity of the curvaton perturbations yields an error proportional to slow-roll parameters. This is irrelevant in the limit , which is the main focus of this work. The WMAP data and other recent observations give the constraint [1, 2, 3, 70].
In this paper we focus on the case in which the curvaton field decays into particles through a parametric resonance, at the end of which the fields undergo a period of potentially very complicated, nonlinear, non-equilibrium dynamics. Afterwards, the fields equilibrate, and we assume that eventually the universe behaves as a mixture of non-interacting matter and radiation. This assumption is checked using simulations in section 4 and we find it to be accurate enough for our current purposes. Therefore, we parameterise the total energy density as a sum of matter and radiation components
where , and are the scale factor, energy density and matter fraction at some arbitrary reference time well after the resonance is over.111 We define following ref. . This differs from the variable used in ref.  by a factor of in the limit .
We assume that the inflaton field has either decayed into ultrarelativistic degrees of freedom or is itself ultrarelativistic, so that it contributes to the radiation component. The field is also ultrarelativistic; the only degree of freedom that contributes to the matter component is the curvaton.
If the resonance destroyed all of the curvaton particles, would be zero, but, in practice, some curvatons are left over. We assume which corresponds to the curvaton being subdominant at the end of the resonance. We can then rearrange (9) to give, at leading order in , the scale factor at energy density
where we have assumed that the curvaton is still subdominant at this time. That is
The curvature perturbation is given by combining (4) with (10). In general, , and all vary between one separate universe and another. In our case, the variation is entirely due to fluctuations of , the value of the curvaton field during inflation. In order to calculate the curvature perturbation using (5) we differentiate (10) with respect to keeping the final energy density fixed. This gives
The remaining curvatons decay within a fixed time where the decay rate depends on the curvaton interactions. We assume that all matter in the universe is ultrarelativistic after the curvaton decay, so that the universe becomes radiation dominated and evolves adiabatically. Following the sudden decay approximation , we assume that this decay takes place instantaneously at a fixed energy density , which is determined by the curvaton decay rate . The final curvature perturbation is therefore given by setting in (12) and (13), or equivalently , where
The decay rate is unknown, so we can treat as a free parameter.
2.3 Standard perturbative decay
As a simple example of the calculation of perturbations, and to aid comparison with our results for resonant curvaton decay, we first summarise the standard perturbative calculation of perturbations generated in the curvaton model [43, 45].
This calculation ignores the resonance, so (9) is assumed to be valid from the start of the curvaton oscillations, which we can therefore choose as our reference time, . This time is determined by the condition , which implies that independently of . Furthermore, the scale factor is also independent of to good approximation, because the curvaton contribution to the energy density is negligible during inflation. Perturbations are therefore generated solely by derivatives of , the matter fraction at the start of the oscillations, and (12) and (13) simplify to
giving the simple result 
3 Linearised calculations
Our general result in (12) and (13) depends on the quantities , and and their derivatives. Ideally, we would like to be able to calculate them analytically for general parameter values. We first attempt to carry out this calculation in linear theory of resonance , i.e. neglecting the backreaction of particles produced during the resonance. This approximation gives a clear physical picture of the resonance, and also an accurate quantitative description of many aspects of it. However, we find that it is not suitable for calculating the curvature perturbation, and, therefore, this calculation acts, ultimately, as a motivation for the nonlinear approach in section 4.
As we assume the field is heavy during inflation and has a vanishing vacuum expectation value, and that the universe becomes radiation dominated after the end of inflation, the equation of motion for the curvaton is given by
which can easily be put into the canonical form of a Bessel equation. The general solution of (19) which remains bounded as is given by
where is a Bessel function of the first kind and denotes the initial curvaton value set by the dynamics during inflation. In the model (1) that we consider here, the value of the curvaton field is constrained by
where the upper bound comes from the subdominance of the curvaton and the lower bound is required to keep the field massive during inflation and to enable broad parametric resonance .
The curvaton remains nearly frozen at the value until the Hubble parameter decreases to and the field starts oscillate around the minimum of its potential. We assume the curvaton is still subdominant at this time and obeys the solution (20). After the onset of oscillations, , (20) can be approximated by the asymptotic expression
Following the notation of , we choose to normalise the scale factor to unity at
Up to corrections , the energy density of the oscillating curvaton is given by
where is the envelope of the oscillatory solution (22) at
The time appears as an unphysical reference point in our analysis but as it formally corresponds to a quarter of the first oscillation cycle it can also be thought of as the beginning of curvaton oscillations, hence the subscript.
The coherent curvaton oscillations induce a time-varying effective mass for the field which can lead to copious production of particles as discussed in . The process is analogous to the standard inflationary preheating scenario [56, 57, 71], except that the universe is dominated by radiation instead of non-relativistic matter during the resonance. The resonant curvaton decay is efficient if
which corresponds to broad resonance bands in momentum space. For curvaton values in the range (21) this condition is always satisfied, as one can immediately see by taking into account the lightness of the curvaton during inflation, .
Using the solution (22) and applying the analytic methods developed to analyse particle production during the parametric resonance , the comoving number density of particles at late times () can be estimated as
where and is an effective growth index. Due to the stochastic nature of parametric resonance in expanding space, the actual growth index depends very sensitively on other parameters and varies rapidly between and between one cycle of oscillation and another.
As the number density of particles grows exponentially, the contribution to the effective curvaton mass eventually comes to dominate over the bare mass . The time when the two contributions are equal—after which the backreaction of the particles can no longer be neglected—is found by equating the number density with . Using the result (27), this yields an equation for
whose solution is given by the branch of Lambert W-function
Here denotes the value of the resonance parameter (26) at the beginning of oscillations. The result is consistent with the assumption for all provided that we place the mild constraint , which corresponds to . In this limit (29) can be approximated to reasonable accuracy by
The derivatives of with respect to , which are needed for computing the curvature perturbation, are
where and the relation following from (25) has been used. The terms involving derivatives of the effective growth index would be very difficult to estimate analytically and we therefore leave them in an implicit form.
Approximate energy conservation gives an estimate
for the resonance parameter (26) at the onset of backreaction. If , the resonance does not terminate immediately at . Instead, the effective curvaton mass starts to grow exponentially, which soon brings the resonance to end. The nonlinear stage of the resonance after therefore makes only a small contribution to the total duration of the resonance  which can be reasonably well estimated by the duration of the linear stage given by (29).
3.2 Curvature perturbation
It is not possible to describe the backreaction and the equilibration of the fields using linear theory. To carry out the linearised calculation, we therefore simply assume that the resonance ends instantaneously at time , and that some fraction of the curvaton energy is transferred into radiation. The matter and radiation energy densities immediately after backreaction are therefore given by
In order to calculate the curvature perturbation, we choose this backreaction time as our reference time in (9). Therefore we have , and
where and we have assumed that , which is equivalent to having throughout the resonance. Substituting these into (12) and (13), we can work out the expression for the curvature perturbation. Using the leading terms in (32), (33), we find
By setting in the above expressions, the standard curvaton results (18) for a perturbative decay are recovered.
The terms scaling as in (39) and (40) represent the contribution due to radiation inhomogeneities created by the curvaton decay. If the curvaton particles left over after the resonance have a long perturbative lifetime, this component effectively redshifts away and the expressions (39) and (40) reduce to
To make use of these results, one would have to know the first and second derivatives of and , which, unfortunately, are difficult to calculate. In principle, the exponential growth rate is calculable in linearised theory, but obtaining its derivatives reliably is hard because it is an averaged quantity that describes a stochastic process. In contrast, is fully determined by the nonlinear dynamics at the end of the resonance, and, therefore, it is not possible to calculate it using linearised equations.
We conclude that the perturbations generated by resonant curvaton decay are not calculable using linear theory. Instead, we have to carry out a fully nonlinear calculation using lattice field theory methods.
4.1 Simulation method
In section 3 we saw that the curvature perturbations produced by the resonance cannot be calculated using linear theory. Therefore, we adopt a completely different approach. We use nonlinear three-dimensional classical lattice field theory simulations222We use a modification of the corrected version of the code used in refs.  and  described in the erratum of . to determine how the scale factor evolves as the fields fall out of equilibrium at the end of preheating and to compute the required quantities , and directly. This is a standard method of solving the field dynamics in such systems [65, 66, 72, 73, 74, 75, 76, 77], and recently [37, 38, 39] it was shown how it can be used, with the separate universe approximation, to compute curvature perturbations.
The and fields are taken to be position-dependent, and the scale factor is homogeneous over the whole lattice. The latter assumption is justified when the lattice is smaller than the Hubble volume, and allows us to describe the expansion of the universe using the Friedmann equation. We have a coupled system of field equations,
and the Friedmann equation
where the energy density is calculated as the average energy density in the simulation box,
The coupled system of equations (44), (45) and (46) are solved on a comoving lattice with periodic boundary conditions. The inflaton component is included as idealised radiation where is the density in the field at the beginning of the simulation. We solve the system in conformal time () using a fourth order Runge-Kutta algorithm. The system of equations is then
The second order Friedmann equation is
where the averaged density and pressure are
The initial inflaton energy density depends on the inflaton potential . For simplicity, we assume it has a quartic form,
We choose the initial value of to be , where we use the subscript to denote the initial values for the simulation. We still have the freedom to choose the value of the conformal time at the start of our simulations. We take
so that the scale factor is simply proportional to conformal time, even at .
The initial values for the and fields consist of a homogeneous background component or , which represents fluctuations with wavelength longer than the lattice size, and Gaussian inhomogeneous fluctuations which represent short-distance quantum fluctuations.
The field is massive during inflation; its long-wavelength fluctuations are heavily suppressed. In contrast, the curvaton field is still frozen at the start of the simulation. As a result of superhorizon fluctuations, its local value at the end of inflation varies between one Hubble patch and another; this sets the initial value for the curvaton field in our simulation, . By running simulations with different initial curvaton values , we can, therefore, determine how , and depend on and calculate the curvature perturbation using (5).
The relevant range of is determined by the power spectrum of the curvaton field. Assuming that the curvaton mass is very small, the power spectrum is given by
where and are the Hubble rate and the number of -foldings before the end of inflation evaluated at the time when the mode left the horizon. A more accurate calculation for non-zero is given in appendix 1 of .
To completely determine the observable curvature perturbation, we need to run the simulation for all values of that were realised in our current Hubble volume. The curvaton field is a Gaussian random field, so the mean value over a Hubble volume today, , has a Gaussian distribution. The variance of this distribution is
Note that this depends on , the total number of -foldings of inflation, so is essentially a free variable. This distribution is centred around zero, assuming no homogeneous classical curvaton background.
The width of the range of that we need to cover is given by the variance of within our current Hubble volume, between different Hubble volumes at the end of inflation. This is given by
where is the number of e-foldings after the largest currently observable scales left the horizon. This gives the width of the range of that we need to consider. Note that the only free parameter this depends on is .
For the initial values of the inhomogeneous field modes, we follow the standard approach [65, 66, 72, 73, 75]. The field is given random initial conditions from a Gaussian distribution whose two-point functions are the same as those in the tree-level quantum vacuum state,
where All other two-point correlators vanish. The inhomogeneous modes of are populated similarly to those of .
In order to estimate the error bars in, for example, figure 3 we repeat each simulation with the same parameters multiple times using different realisations of the above random initial conditions. In earlier work on preheating [37, 38, 39], the whole evolution was calculated using full nonlinear equations. In the current case that would be computationally very expensive, because the model is not conformally invariant and the universe expands by a large factor during the evolution. To fit the relevant wavelengths inside the lattice throughout the whole simulation, the lattice size would have to be much larger than , which is not possible given the number of simulations which need to be run.
Instead, we take a shortcut, and make use of the fact that until the magnitude of the fluctuations in become large enough to backreact on (when ) the field dynamics are linear to a very good approximation. Therefore, we can evolve the whole probability distributions of the initial conditions in -space using linearised versions of (48)-(52). More precisely, we need to solve the linear equations for the inhomogeneous modes,
in the background provided by the solutions of the nonlinear equations for the homogeneous modes,
Crucially, the general solution of the linear mode (59) and (60) is obtained by solving the equations for two different initial values. This is because, as a consequence of linearity, the general solution can be represented as
where is a two-by-two matrix. Consider now the solutions with initial conditions of position 1 and velocity 0 and vice versa:
Linearity implies that the solution of arbitrary initial conditions and is given simply by
We use this approach to evolve the fields until some scale factor , which is well before the nonlinear terms become important. Only then, we Fourier transform the fields to coordinate space and start to evolve them using the full nonlinear equations (48)-(52). This is shown in figure 1. We checked that our results do not depend on the chosen value of .
This method has two advantages which make the calculation presented here possible. Firstly, we evolve the whole distribution, so the linear evolution only needs to be calculated once, and only the relatively short nonlinear stage has to be repeated for each realisation of the initial conditions. Secondly, and more significantly, it reduces the factor by which the universe grows during the evolution of (48)-(52) from 100 to 10. This allow us to use lattices with much lower resolution to achieve comparably accurate results, thus significantly reducing computation time.
4.2 Simulation results
Our model has three free parameters, , , . It is not possible for us to probe this whole parameter space using current computational technology. Therefore, we simply pick one set of parameters, , and . We take .
We also have to choose the initial value of the curvaton field . In order to cover all values of in our currently observable universe, we carried out simulations for a range of
whose width is given by (57),
The probability distribution of the central value is given by (56), but as a results of its dependence on the total number of e-foldings , which is unknown, we can essentially choose it freely. Nevertheless, the value of is restricted by several constraints. The resonance being brought to an end by backreaction demands that (see (26)). The characteristic frequency of the field oscillations is and the time the system takes to reach backreaction is proportional to . Therefore, the simulation time is proportional to . In our simulations we used three values:
These corresponds to 50, 200 and 750 respectively. With these parameters one point lattice simulation takes 3 hours. We repeat each simulation several (50) times with different random realisations of initial conditions to obtain statistical errors. Although, this computational method is very resource intensive, the independence of each simulation (or Hubble patch) means that the problem scales perfectly on a multiprocessor machine. A time-step of and lattice spacing of were used.
Each simulation gives us time-streams of the scale factor , the density , given by (51), and the pressure , given by (52). We measure the matter fraction in two different ways: from the ratio of the pressure and energy density,
and by fitting the dependence of the scale factor on the energy density,
in a narrow range about the point of interest. In continuous time these two expressions would agree, and therefore any discrepancy between them would be a sign of time discretisation errors.
The measured matter fraction is shown in figure 2 as a function of the energy density. The two curves correspond to the two definitions (73) and (74), which clearly agree very well. Initially, it grows because the matter component decreases more slowly than the radiation component. At , it drops rapidly as the resonance destroys most of the curvatons. Afterwards again grows as the universe expands.
To calculate the curvature perturbation using (12) and (13), we need to choose a reference point at which we measure the scale factor , energy density and matter fraction . We chose this be at fixed energy density . We interpolated our measurements using a simple least-squares fit in and about to find the scale factor and the matter fraction . The results are shown in figure 3 for three choices of .
We found that the choice of has no effect on the results as long as it is well after the end of the resonance. This is demonstrated by figure 2, in which the (blue) dashed curve shows the matter fraction calculated from the assumption (9) of non-interacting matter and radiation components with the measured values of and . This agrees very well with the measured at late times.
We calculate the first and second derivatives of and with respect to . Assuming that the expansion (5) works well over the whole range of