Moment transport equations for nonGaussianity
Abstract
We present a novel method for calculating the primordial nonGaussianity produced by superhorizon evolution during inflation. Our method evolves the distribution of coarsegrained inflationary field values using a transport equation. We present simple evolution equations for the moments of this distribution, such as the variance and skewness. This method possesses some advantages over existing techniques. Among them, it cleanly separates multiple sources of primordial nonGaussianity, and is computationally efficient when compared with popular alternatives, such as the framework. We adduce numerical calculations demonstrating that our new method offers good agreement with those already in the literature. We focus on two fields and the parameter, but we expect our method will generalize to multiple scalar fields and to moments of arbitrarily high order. We present our expressions in a fieldspace covariant form which we postulate to be valid for any number of fields.
Keywords: Inflation, Cosmological perturbation theory, Physics of the early universe, Quantum field theory in curved spacetime.
pacs:
98.80.k, 98.80.Cq, 11.10.Hi, ,
1 Introduction
Inflation generically predicts a primordial spectrum of density perturbations which is almost precisely Gaussian [1, 2, 3, 4, 5]. In recent years the small nonGaussian component [6, 7, 8, 9, 10, 11, 12] has emerged as an important observable [13], and will be measured with good precision by the Planck Surveyor satellite [14]. In the near future, as observational data become more plentiful, it will be important to understand the nonGaussian signal expected in a wide variety of models, and to anticipate what conclusions can be drawn about earlyuniverse physics from a prospective detection of primordial nonGaussianity.
In this paper, we present a novel method for calculating the primordial nonGaussianity produced by superhorizon evolution in twofield models of inflation. Our method is based on the realspace distribution of inflationary field values on a flat hypersurface, which can be thought of as a probability density function whose evolution is determined by a form of the collisionless Boltzmann equation. Using a cumulant representation [15, 16, 17, 18] to expand our density function around an exact Gaussian, we derive ordinary differential equations which evolve the moments of this distribution. Further, we show how these moments are related to observable quantities, such as the dimensionless bispectrum measured by [14, 10]. We present numerical results which show that this method gives good agreement with other techniques. It is not necessary to make any assumptions about the inflationary model beyond requiring a canonical kinetic term and applying the slowroll approximation. While there are already numerous methods for computing the superhorizon contribution to , including the widely used formalism, we believe the one reported here has a number of advantages.
First, this new technique is ideally suited to unraveling the various contributions to . This is because we follow the moments of the inflaton distribution directly, which makes it straightforward to identify large contributions to the skewness or other moments. The evolution equation for each moment is simple and possesses clearly identifiable source terms, which are related to the properties of the inflationary flow on field space. This offers a clear separation between two key sources of primordial nonGaussianity. One of these is the intrinsic nonlinearity associated with evolution of the probability density function between successive flat hypersurfaces; the other is a gauge transformation from field fluctuations to the curvature peturbation, . It would be difficult or impossible to observe this split within the context of other calculational schemes, such as the conventional formalism.
A second advantage of our method is connected with the computational cost of numerical implementation. Analytic formulas for are known in certain cases, mostly in the context of the framework, but only for very specific choices of the potential [19, 20, 21, 22] or Hubble rate [23, 24]. These formulas become increasingly cumbersome as the number of fields increases, or if one studies higher moments [25, 26]. In the future, it seems clear that studies of complex models with many fields will increasingly rely on numerical methods. The numerical framework requires the solution to a number of ordinary differential equations which scales exponentially with the number of fields. Since some models include hundreds of fields, this may present a significant obstacle [27]. Moreover, the formalism depends crucially on a numerical integration algorithm with low noise properties, since finite differences must be extracted between perturbatively different initial conditions after efolds of evolution. Thus, the background equations must be solved to great accuracy, since any small error has considerable scope to propagate.
In this paper we ultimately solve our equations numerically to determine the evolution of moments in specific models. Our method requires the solution to a number of differential equations which scales at most polynomially (or in certain cases perhaps even linearly) with the number of fields. It does not rely on extracting finite differences, and therefore is much less susceptible to numerical noise. These advantages may be shared with other schemes, such as the numerical method recently employed by Lehners & RenauxPetel [28].
A third advantage, to which we hope to return in a future publication, is that our formalism yields explicit evolution equations with source terms. From an analysis of these source terms, we hope that it will be possible to identify those physical features of specific models which lead to the production of large nonGaussianities.
This paper is organized as follows. In §2, we show how the nonGaussian parameter can be computed in our framework. The calculation remains in real space throughout (as opposed to Fourier space), which modifies the relationship between and the multipoint functions of the inflaton field. Our expression for shows a clean separation between different contributions to nonGaussianity, especially between the intrinsic nonlinearity of the field evolution and the gauge transformation between comoving and flat hypersurfaces. In §3, we introduce our model for the distribution of inflaton field values, which is a “moment expansion” around a purely Gaussian distribution. We derive the equations which govern the evolution of the moments of this distribution in the one and twofield cases. In §4, we present a comparison of our new technique and those already in the literature. We compute numerically in several twofield models, and find excellent agreement between techniques. We conclude in §5.
Throughout this paper, we use units in which , and the reduced Planck mass is set to unity.
2 Frameworks for computing
In this section, we introduce our new method for computing the nonGaussianity parameter . This method requires three main ingredients: a formula for the curvature perturbation, , in terms of the field values on a spatially flat hypersurface; expressions for the derivatives of the number of efoldings, , as a function of field values at horizon exit; and a prescription for evolving the field distribution from horizon exit to the time when we require the statistical properties of . The first two ingredients are given in Eqs. (9)–(10) and (27)–(2.3), found at the end of §2.2 and §2.3 respectively. The final ingredient is discussed in §3.
2.1 Calculations beyond linear order
Once it became clear that nonlinearities of the microwave background anisotropies could be detected by the WMAP and Planck survey satellites [14], many authors studied higherorder correlations of the curvature perturbation. In early work, direct calculations of a correlation function were matched to the known limit of local nongaussianity [9, 10, 29, 30, 31, 32]. This method works well if isocurvature modes are absent, so that the curvature perturbation is constant after horizon exit. In the more realistic situation that isocurvature modes cause evolution on superhorizon scales, all correlation functions become time dependent. Various formalisms have been employed to describe this evolution. Lyth & Rodríguez [11] extended the method [33, 34] beyond linear order. This method is simple and wellsuited to analytical calculation. Rigopoulos, Shellard and van Tent [35, 36] worked with a gradient expansion, rewriting the field equations in Langevin form. The noise term was used as a proxy for setting initial conditions at horizon crossing. A similar ‘exact’ gradient formalism was written down by Langlois & Vernizzi [37, 38, 39]. In its perturbative form, this approach has been used by Lehners & RenauxPetel to obtain numerical results [28]. Another numerical scheme has been introduced by Huston & Malik [40].
What properties do we require of a successful prediction? Consider a typical observer, drawn at random from an ensemble of realizations of inflation. In any of the formalisms discussed above, we aim to estimate the statistical properties of the curvature perturbation which would be measured by such an observer. Some realizations may yield statistical properties which are quite different from the ensemble average, but these large excursions are uninteresting unless anthropic arguments are in play.
Next we introduce a collection of comparablysized spacetime volumes whose mutual scatter is destined to dominate the microwave background anisotropy on a given scale. Neglecting spatial gradients, each spacetime volume will follow a trajectory in field space which is slightly displaced from its neighbors. The scatter between trajectories is determined by initial conditions set at horizon exit, which are determined by promoting the vacuum fluctuation to a classical perturbation.^{1}^{1}1The scatter among trajectories may grow with the overall volume, owing to backreaction effects in de Sitter space [41, 42, 43, 44]. If inflation can end in more than one vacuum, leading to different predictions for observable quantities, it may be helpful to evaluate this enhanced dispersion for the purpose of deciding into which vacuum the field will fall. Once this minimum has been decided, however, largevolume effects needlessly complicate the calculation. In this paper, we will always assume that predictions are being made by following a sufficient number of trajectories to determine the statistical properties with reasonable precision, but no more. A correct prediction is a function of the trajectories followed by every volume in the collection, taken as a whole. One never makes a prediction for a single trajectory.
Each spacetime volume follows a trajectory, which we label with its position at some fixed time, to be made precise below. Throughout this paper, superscript ‘’ denotes evaluation on a spatially flat hypersurface. Consider the evolution of some quantity of interest, , which is a function of trajectory. If we know the distribution we can study statistical properties of such as the moment ,
(1) 
where we have introduced the ensemble average of ,
(2) 
In Eqs. (1)–(2), stands for a collection of any number of fields. It is the which are observable quantities.
Eq. (1) defines what we will call the exact separate universe picture. It is often convenient to expand as a power series in the field values centered on a fiducial trajectory, labelled ‘fid,’
(3) 
When Eq. (3) is used to evaluate the , we refer to the ‘perturbative’ separate universe picture. If all terms in the power series are retained, these two versions of the calculation are formally equivalent. In unfavorable cases, however, convergence may occur slowly or not at all. This possibility was discussed in Refs. [45, 46]. Although our calculation is formally perturbative, it is not directly equivalent to Eq. (3). We briefly discuss the relation of our calculation to conventional perturbation theory in §5.
2.2 Calculating in the separate universe picture
By definition, the curvature perturbation measures local fluctuations in expansion history (expressed in efolds ), calculated on a comoving hypersurface. In many models, the curvature perturbation is synthesized by superhorizon physics, which reprocesses a set of Gaussian fluctuations generated at horizon exit. In a singlefield model, only one Gaussian fluctuation can be present, which we label . Neglecting spatial gradients, the total curvature perturbation must then be a function of alone. For , this can be wellapproximated by
(4) 
where is independent of spatial position. Eq. (4) defines the socalled “local” form of nongaussianity. It applies only when quantum interference effects can be neglected, making a welldefined object rather than a superposition of operators [47]. If this condition is satisfied, spatial correlations of may be extracted and it follows that can be estimated using the rule
(5) 
where we have recalled that is nearly Gaussian, or equivalently that .
With spatially independent, Eq. (4) strictly applies only in singlefield inflation. In this case one can accurately determine by applying Eq. (4) to a single trajectory with fixed initial conditions, as in the method of Lehners & RenauxPetel [28]. Where more than one field is present, may vary in space because it depends on the isocurvature modes. In this case one must determine statistically on a bundle of adjacent trajectories which sample the local distribution of isocurvature modes. Eq. (5) is then indispensible. Following Maldacena [10], and later Lyth & Rodríguez [11], we adopt Eq. (5) as our definition of , whatever its origin. In real space, the coefficient in Eq. (5) depends on the convention . More generally, this follows from the definition of , Eq. (1). In Fourier space, either prescription is automatically enforced after dropping disconnected contributions, again leading to Eq. (5).
To proceed, we require estimates of the correlation functions and . We first describe the conventional approach, in which ‘’ denotes a flat hypersurface at a fixed initial time. The quantity denotes the number of efoldings between this initial slice and a final comoving hypersurface, where indexes the species of light scalar fields. The local variation in expansion can be written in the fiducial picture as
(6) 
where .
Subject to the condition that the relevant scales are all outside the horizon, we are free to choose the initial time—set by the hypersurface ‘’—at our convenience. In the conventional approach, ‘’ is taken to lie a few efolds after our collection of spacetime volumes passes outside the causal horizon [11, 12]. This choice has many virtues. First, we need to know statistical properties of the field fluctuations only around the time of horizon crossing, where they can be computed without the appearance of large logarithms [44, 48]. Second, as a consequence of the slowroll approximation, the are uncorrelated at this time, leading to algebraic simplifications. Finally, the formula subsumes a gauge transformation from the field variables to the observational variable . Using Eqs. (1)–(2), (5) and (6), one finds that can be written to a good approximation [11]
(7) 
where and for simplicity we have dropped the ‘’ which indicates time of evaluation. A similar definition applies for .
Eq. (7) is accurate up to small intrinsic nonGaussianities present in the field fluctuations at horizon exit. As a means of predicting it is pleasingly compact, and straightforward to evaluate in many models. Unfortunately, it also obscures the physics which determines . For this reason it is hard to infer, from Eq. (7) alone, those classes of models in which is always large or small.^{2}^{2}2Even in simple models it can be quite subtle to determine what range of is dynamically allowed. See, for example, Refs. [21, 22].
Our strategy is quite different. We choose ‘’ to lie around the time when we require the statistical properties of . The role of the formula, Eq. (6), is then to encode only the gauge transformation between the and . In §2.3 below, we show how the appropriate gauge transformation is computed using the formula. In the present section we restrict our attention to determining a formula for under the assumption that the distribution of field values on ‘’ is known. In §3, we will supply the required prescription to evolve the distribution of field values between horizon exit and ‘’.
Although the are independent random variables at horizon exit, correlations can be induced by subsequent evolution. One must therefore allow for offdiagonal terms in the twopoint function. Remembering that we are working with a collection of spacetime volumes in real space, smoothed on some characteristic scale, we write
(8) 
does not vary in space, but it may be a function of the scale which characterizes our ensemble of spacetime volumes. In all but the simplest models it will vary in time. It is also necessary to account for intrinsic nonlinearities among the , which are small at horizon crossing but may grow. We define
(9) 
Likewise, should be regarded as a function of time and scale. The permutation symmetries of an expectation value such as (9) guarantee that, for example, .^{3}^{3}3We are assuming that these expectation values are ensemble averages over classical stochastic fields, and are therefore invariant under reordering of the fields. When written explicitly, we place the indices of symbols such as in numerical order. Neglecting a small () intrinsic fourpoint correlation, it follows that
(10) 
Now we specialize to a twofield model, parametrized by fields and . Using Eqs. (1)–(2), (6) and (8), it follows that the twopoint function of satisfies
(11) 
The threepoint function can be written
(12) 
where we have identified two separate contributions, labelled ‘1’ and ‘2’. The ‘1’ term includes all contributions involving intrinsic nonlinearities, those which arise from nonGaussian correlations among the field fluctuations,
(13) 
The ‘2’ term encodes nonlinearities arising directly from the gauge transformation to
(14) 
After use of Eq. (5), Eqs. (12)–(14) can be used to extract the nonlinearity parameter . This decomposes likewise into two contributions , which we shall discuss in more detail in §4.
2.3 The derivatives of
To compute in concrete models, we require expressions for the derivatives and . For generic initial and final times, these are difficult to obtain. Lyth & Rodríguez [11] used direct integration, which is effective for quadratic potentials and constant slowroll parameters. Vernizzi & Wands [19] obtained expressions in a twofield model with an arbitrary sumseparable potential by introducing Gaussian normal coordinates on the space of trajectories. Their approach was generalized to many fields by Battefeld & Easther [20]. Productseparable potentials can be accommodated using the same technique [49]. An alternative technique has been proposed by Yokoyama et al. [50].
A considerable simplification occurs in the present case, because we only require the derivative evaluated between flat and comoving hypersurfaces which coincide in the unperturbed universe. For any species , and to leading order in the slowroll approximation, the number of efolds between the flat hypersurface ‘’ and a comoving hypersurface ‘’ satisfies
(15) 
where and are the field values evaluated on ‘’ and ‘,’ respectively. Under an infinitesimal shift of , we deduce that obeys
(16) 
Note that this applies for an arbitrary , which need not factorize into a sum or product of potentials for the individual species . In principle a contribution from variation of the integrand is present, which spoils a naïve attempt to generalize the method of Refs. [19, 20, 49] to an arbitrary potential. This contribution vanishes in virtue of our supposition that ‘’ and ‘’ are infinitesimally separated.
To compute it is helpful to introduce a quantity , which in the sumseparable case coincides with the conserved quantity of Vernizzi & Wands [19, 51]. For our specific choice of a twofield model, this takes the form
(17) 
where the integrals are evaluated on a single spatial hypersurface. In an field model, one would obtain conserved quantities which label the isocurvature fields. The construction of these quantities is discussed in Refs. [20, 25]. For sumseparable potentials one can show using the equations of motion that is conserved under time evolution to leading order in slowroll. It is not conserved for general potentials, but the variation can be neglected for infinitesimally separated hypersurfaces.
Under a change of trajectory, varies according to the rules
(18) 
and
(19) 
The comoving hypersurface ‘’ is defined by
(20) 
We are assuming that the slowroll approximation applies, so that the kinetic energy may be neglected in comparison with the potential . Therefore on ‘’ we have
(21) 
Combining Eqs. (18), (19) and (21) we obtain expressions for , namely
(22)  
(23)  
(24)  
(25) 
where we have defined
(26) 
Eqs. (22)–(25) can alternatively be derived without use of by comparing Eq. (16) with the formulas of Ref. [52], which were derived using conventional perturbation theory. Applying (16), we obtain
(27) 
To proceed, we require the second derivatives of . These can be obtained directly from (27), after use of Eqs. (22)–(25). We find
(28) 
An analogous expression for can be obtained after the simultaneous exchange . The mixed derivative satisfies
(29) 
Now that the calculation is complete, we can drop the superscripts ‘’ and ‘,’ since any background quantity is the same on either hypersurface. Once this is done it can be verified that (despite appearances) Eq. (2.3) is invariant under the exchange .
3 Transport equations
In this section we return to the problem of evolution between horizon exit and the time of observation, and supply the prescription which connects the distribution of field values at these two times.
3.1 Nongaussian distribution in the singlefield case
We begin by discussing the singlefield system, which lacks the technical complexity of the twofield case, yet still exhibits certain interesting features which recur there. Among these features are the subtle difference between motion of the statistical mean and the background field value, and the hierarchy of moment evolution equations. Moreover, the structure of the moment mixing equations is similar to that which obtains in the twofield case. For this reason, the onefield scenario provides an instructive example of the techniques we wish to employ.
Recall that we work in real space with a collection of comparably sized spacetime volumes, each with a slightly different expansion history, and the scatter in these histories determines the microwave background anisotropy on a given angular scale. Within each volume the smoothed background field takes a uniform value described by a density function , where in this section we are dropping the superscript ‘’ denoting evaluation of spatially flat hypersurfaces. Our ultimate goal is to calculate the reduced bispectrum, , which describes the third moment of . In the language of probability this is the skewness, which we denote . A Gaussian distribution has skewness zero, and inflation usually predicts that the skew is small. For this reason, rather than seek a distribution with nonzero third moment, as proposed in Ref. [18], we will introduce higher moments as perturbative corrections to the Gaussian. Such a procedure is known as a cumulant expansion.
The construction of cumulant expansions is a classical problem in probability theory. We seek a distribution with centroid , variance , and skew , with all higher moments determined by and alone. A distribution with suitable properties is
(30) 
where
(31) 
is a pure Gaussian and denotes the Hermite polynomial, for which there are multiple normalization conventions. We choose to normalize so that
(32) 
which implies that the leading term of is . This is sometimes called the “Probabilist’s convention.” We define expectation values by the usual rule,
(33) 
The probability density function in Eq. (30) has the properties^{1}^{1}1These formulas apply for arbitrary values of , and do not depend on the approximation that is small. However, for large the density function (30) may become negative for some values of . It then ceases to be a probability density in the strict sense. This does not present a problem in practice, since we are interested in distributions which are approximately Gaussian, and for which will typically be small. Moreover, our principal use of Eq. (30) is as a formal tool to extract evolution equations for each moment. For this reason we will not worry whether defines an honest probability density function in the strict mathematical sense.
(34) 
The moments , , and may be timedependent, so evolution of the probability density in time can be accommodated by finding evolution equations for these quantities.
The density function given in Eq. (30) is wellknown and has been applied in many situations. It is a solution to the problem of approximating a nearlyGaussian distribution whose moments are known. Eq. (30) is in fact the first two terms of the Gram–Charlier ‘A’ series, also sometimes called the Gauss–Hermite series.^{2}^{2}2In the physics literature, this series has sometimes erroneously been called the Edgeworth expansion. In recent years it has found multiple applications to cosmology, of which our method is closest to that of Taylor & Watts [53]. Other applications are discussed in Refs. [16, 17, 15, 54, 18, 53, 55, 56, 57, 58, 59, 60]. For a review of the ‘A’ series and related nearlyGaussian probability distributions from an astrophysical perspective, see [61]. In this paper, we will refer to Eq. (30) and its natural generalization to higher moments as the “moment expansion.”
In the slowroll approximation, the field in each spacetime volume obeys a simple equation of motion
(35) 
where records the number of efoldings of expansion. We refer to as the velocity field. Expanding about the instantaneous centroid gives
(36) 
where
(37) 
The value of evolves with time, so each expansion coefficient is timedependent. Hence, we do not assume that the velocity field is globally welldescribed by a quadratic Taylor expansion, but merely that it is welldescribed as such in the neighborhood of the instantaneous centroid. We expand the velocity field to second order, although in principle this expansion could be carried to arbitrary order.
It remains to specify how the probability density evolves in time. Conservation of probability leads to the transport equation
(38) 
Eq. (38) can also be understood as the limit of a Chapman–Kolmogorov process as the size of each hop goes to zero. It is well known—for example, from the study of Starobinsky’s diffusion equation which forms the basis of the stochastic approach to inflation [62]—that the choice of time variable in this equation is significant, with different choices corresponding to the selection of a temporal gauge. We have chosen to use the efolding time, , which means that we are evolving the distribution on hypersurfaces of uniform expansion. These are the spatially flat hypersurfaces whose field perturbations enter the formulas described in §2.
In principle, Eq. (38) can be solved directly. In practice it is simpler to extract equations for the moments of , giving evolution equations for , and . To achieve this, one need only resolve Eq. (38) into a Hermite series of the form
(39) 
The Hermite polynomials are linearly independent, and application of the orthogonality condition (32) shows that the must all vanish. This leads to a hierarchy of equations , which we refer to as the moment hierarchy. At the top of the hierarchy, the equation is empty and expresses conservation of probability.
The first nontrivial equation requires and yields an evolution equation for the centroid ,
(40) 
The first term on the righthand side drives the centroid along the velocity field, as one would anticipate based on the background equation of motion, Eq. (35). However, the second term shows that the centroid is also influenced as the wings of the probability distribution probe the nearby velocity field. This influence is not captured by the background equation of motion. If we are in a situation with , then the wings of the density function will be moving faster than the center. Hence, the velocity of the centroid will be larger than one might expect by restricting attention to . Accordingly, the mean fluctuation value is not following a solution to the background equations of motion.
Evolution equations for the variance and skew are obtained after enforcing , yielding
(41)  
(42) 
In both equations, the first term on the righthand sides describes how and scale as the density function expands or contracts in response to the velocity field. These terms force and to scale in proportion to the velocity field. Specifically, if we temporarily drop the second terms in each equation above, one finds that and . This precisely matches our expectation for the scaling of these quantities. Hence, these terms account for the Jacobians associated with infinitesimal transformations induced by the flow .
For applications to inflationary nonGaussianity, the second terms in (41) and (42) are more relevant. These terms describe how each moment is sourced by higher moments and the interaction of the density function with the velocity field. In the example above, if we are in a situation where , the tails of the density function are moving faster than the core. This means that one tail is shrinking and the other is extending, skewing the probability density. The opposite occurs when . These effects are measured by the second term in (42). Hence, by expanding our PDF to the third moment, and our velocity field to quadratic order, we are able to construct a set of evolution equations which include the leadingorder source terms for each moment.
3.2 The twofield case
There is little conceptually new as we move from one field to two. The new features are mostly technical in nature. Our primary challenge is a generalization of the moment expansion to two fields, allowing for the possibility of correlation between the fields. With this done, we can write down evolution equations whose structure is very similar to those found in the singlefield case.
The twofield system is described by a twodimensional velocity field , defined by
(43) 
where again we are using the number of efolds as the time variable. The index takes values in . While we think it is likely that our equations generalize to any number of fields, we have only explicitly constructed them for a twofield system. As will become clear below, certain steps in this construction apply only for two fields, and hence we make no claims at present concerning examples with three or more fields.
The twodimensional transport equation is
(44) 
Here and in the following we have returned to our convention that repeated species indices are summed. As in the singlefield case, we construct a probability distribution which is nearly Gaussian, but has a small nonzero skewness. That gives
(45) 
where is a pure Gaussian distribution, defined by
(46) 
In this equation, defines the center of the distribution and describes the covariance between the fields. We adopt a conventional parametrization in terms of variances and a correlation coefficient ,
(47) 
The matrix defines twopoint correlations of the fields,
(48) 
All skewnesses are encoded in . Before defining this explicitly, it is helpful to pause and notice a complication inherent in Eqs. (46)–(47) which was not present in the singlefield case. To extract a hierarchy of moment evolution equations from the transport equation, Eq. (38), we made the expansion given in (39) and argued that orthogonality of the Hermite polynomials implied the hierarchy . However, Hermite polynomials of the form are not orthogonal under the Gaussian measure of Eq. (46). Following an expansion analogous to Eq. (39) the moment hierarchy would comprise linear combinations of the coefficients. The problem is essentially an algebraic question of Gram–Schmidt orthogonalization.
To avoid this problem it is convenient to diagonalize the covariance matrix , introducing new variables and for which Eq. (46) factorizes into the product of two measures under which the polynomials and are separately orthogonal. The necessary redefinitions are
(49) 
and
(50) 
A simple expression for can be given in terms of and ,
(51) 
We now define the nonGaussian factor, which encodes the skewnesses, to be
(52) 
In these variables we find , but . In addition, we have
(53) 
In order for Eq. (52) to be useful, it is necessary to express the skewnesses associated with the physical variables in terms of and . By definition, these satisfy
(54) 
After substituting for the definition of these quantities inside the expectation values in Eq. (53) we arrive at the relations
(55)  
(56)  
(57)  
(58) 
The moments , and are timedependent, but for clarity we will usually suppress this in our notation.
Next we must extract the moment hierarchy, which governs evolution of , , and . We expand the velocity field in a neighborhood of the instantaneous centroid according to
(59) 
where we have defined
(60) 
As in the singlefield case, these coefficients are functions of time and vary with the motion of the centroid. The expansion can be pursued to higher order if desired.
Our construction of and implies that the twofield transport equation can be arranged as a double Gauss–Hermite expansion,
(61) 
Because the Hermite polynomials are orthogonal in the measure defined by , we deduce the moment hierarchy
(62) 
We define the “rank” of each coefficient by . We terminated the velocity field expansion at quadratic order, and our probability distribution included only the first three moments. It follows that only with rank five or less are nonzero. If we followed the velocity field to higher order, or included higher terms in the moment expansion, we would obtain nontrivial higherrank coefficients. Inclusion of additional coefficients requires no qualitative modification of our analysis and can be incorporated in the scheme we describe below.
A useful feature of the expansion in Eq. (61) is that the rank coefficients give evolution equations for the order moments. Written explicitly in components, the expressions that result from (61) are quite cumbersome. However, when written as fieldspace covariant expressions they can be expressed in a surprisingly compact form.
 Rank 0

The rank0 coefficient is identically zero. This expresses the fact that the total probability is conserved as the distribution evolves.
 Rank 1

The rank1 coefficients and give evolution equations for the centroid . These equations can be written in the form
(63) We remind the reader that here and below, terms like , and represent the velocity field and its derivatives evaluated at the centroid . The first term in (63) expresses the nonanomalous motion of the centroid, which coincides with the background velocity field of Eq. (43). The second term describes how the wings of the probability distribution sample the velocity field at nearby points. Narrow probability distributions have small components of and hence are only sensitive to the local value of . Broad probability distributions have large components of and are therefore more sensitive to the velocity field far from the centroid.
 Rank 2

The rank2 coefficients , and give evolution equations for the variances and the correlation . These can conveniently be packaged as evolution equations for the matrix
(64) This equation describes the stretching and rotation of as it is transported by the velocity field. It includes a sensitivity to the wings of the probability distribution, in a manner analogous to the similar term appearing in (41). Hence the skew acts as a source for the correlation matrix.
 Rank 3

The rank3 coefficients , , and describe evolution of the moments . These are
(65) The first term describes how the moments flow into each other as the velocity field rotates and shears the coordinate frame relative to the coordinate frame. The second term describes sourcing of nonGaussianity from inhomogeneities in the velocity field and the overall spread of the probability distribution.
Some higherrank coefficients—in our case, those of ranks four and five—are also nonzero, but do not give any new evolution equations. These coefficients measure the “error” introduced by truncating the moment expansion. If we had included higher cumulants, these higherrank coefficients would have given evolution equations for the higher moments of the probability distribution. In general, all moments of the density function will mix so it is always necessary to terminate our expansion at a predetermined order—both in cumulants and powers of the field fluctuation. The order we have chosen is sufficient to generate evolution equations containing both the leadingorder behavior of the moments—namely, the first terms in Eqs. (63), (64) and (65)—and the leading corrections, given by the latter terms in these equations.
4 Numerical results
At this point we put our new method into practice. We study two models for which the nonGaussian signal is already known, using the standard formula. For each case we employ our method and compare it with results obtained using . To ensure a fair comparison, we solve numerically in both cases. Our new method employs the slowroll approximation, as described above. Therefore, when using the approach we produce results both with and without slowroll simplifications.
First consider double quadratic inflation, which was studied by Rigopoulos, Shellard & van Tent [36, 63] and later by Vernizzi & Wands [19]. The potential is
(66) 
We use the initial conditions chosen in Ref. [36], where , and the fiducial trajectory has coordinates and .
We plot the evolution of in Fig. 1, which also shows the prediction of the standard formula (with and without employing slow roll simplifications). We implement the algorithm using a finite difference method to calculate the derivatives of . A similar technique was used in Ref. [19]. This model yields a very modest nonGaussian signal, below unity even at its peak. If inflation ends away from the spike then is practically negligible.
Eq. (12) shows that the method of moment transport allows us to separate contributions to from the intrinsic nonGaussianity of the field fluctuations, and nonlinearities of the gauge transformation to . As explained in §2.2, we denote the former and the latter , and plot them separately in Fig. 2.
Inspection of this figure clearly shows that is determined by a cancellation between two much larger components. Its final shape and magnitude are exquisitely sensitive to their relative phase. Initially, the magnitudes of and grow, but their sum remains small. The peak in Fig. 1 arises from the peak of , which is incompletely cancelled by . It is remarkable that initially evolves in exact opposition to the gauge transformation, to which it is not obviously connected.
In the double quadratic model, is always small. However, it has recently been shown by Byrnes et al. that a large nonGaussian signal can be generated even when slowroll is a good approximation [21, 22]. The conditions for this to occur are incompletely understood, but apparently require a specific choice of potential and strong tuning of initial conditions. In Figs. 3–4 we show the evolution of in a model with the potential
(67) 
which corresponds to Example A of Ref. [21, §5] when we choose and initial conditions , . It is clear that the agreement is exact.
In this model, is overwhelmingly dominated by the contribution from the secondorder gauge transformation, , as shown in Fig. 4.
This conclusion applies equally to the other large examples discussed in Refs. [21, 22], although we make no claim that this is a general phenomenon.
In conclusion, Figs. 1 and 3 show excellent agreement between our new method and the outcome of the numerical formula. These figures also compare the moment transport method and without the slowroll approximation. We conclude that the slowroll estimate remains broadly accurate throughout the entire evolution.
5 Discussion
Nonlinearities are now routinely extracted from allsky observations of the microwave background anisotropy. Our purpose in this paper has been to propose a new technique with which to predict the observable signal. Present data already give interesting constraints on the skewness parameter , and over the next several years we expect that the Planck survey satellite will make these constraints very stringent. It is even possible that higherorder moments, such as the kurtosis parameter [64] will become better constrained [65]. To meet the need of the observational community for comparison with theory, reliable estimates of these nonlinear quantities will be necessary for various models of earlyuniverse physics.
A survey of the literature suggests that the ‘conventional’ method, originally introduced by Lyth & Rodríguez, remains the method of choice for analytical study of nonGaussianity. In comparison, our proposed moment transport method exhibits several clear differences. First, the conventional method functions best when we base the expansion on a flat hypersurface immediately after horizon exit. In our method, we make the opposite choice and move the flat hypersurface as close as possible to the time of observation. After this, the role of the formula is to provide no more than the nonlinear gauge transformation between field fluctuations and the curvature perturbation. We substitute the method of moment transport to evolve the distribution of field fluctuations between horizon exit and observation.
Second, in integrating the transport equation one uses an expansion of the velocity field such as the one given in Eqs. (59)–(60). This expansion is refreshed at each step of integration, so the result is related to conventional perturbative calculations in a very similar way to renormalizationgroup improved perturbation theory [66]. In this interpretation, derivatives of play the role of couplings. At a given order, , in the moment hierarchy, the equations for lowerorder moments function as renormalization group equations for the couplings at level, resumming potentially large terms before they spoil perturbation theory. This property is shared with any formalism such as which is nonperturbative in time evolution, but may be an advantage in comparison with perturbative methods. We also note that although is nonperturbative as a point of principle, practical implementations are frequently perturbative. For example, the method of Vernizzi & Wands [19] and Battefeld & Easther [20] depends on the existence of quantities which are conserved only to leading order in , and can lose accuracy after efoldings.
Numerical calculations confirm that our method gives results in excellent agreement with existing techniques. As a byproduct of our analysis, we note that the large nongaussianities which have recently been observed in sum and productseparable potentials [21, 22] are dominated by nonlinearities from the secondorder part of the gauge transformation from to . The contribution from intrinsic nonlinearities of the field fluctuations, measured by the skewnesses , is negligible. In such cases one can obtain a useful formula for by approximating the field distribution as an exact Gaussian. The nonGaussianity produced in such cases arises from a distortion of comoving hypersurfaces with respect to adjacent spatially flat hypersurfaces.
Our new method joins many wellestablished techniques for estimating nonGaussian properties of the curvature perturbation. In our experience, these techniques give comparable estimates of , but they do not exactly agree. Each method invokes different assumptions, such as the neglect of gradients or the degree to which time dependence can be accommodated. The mutual scatter between different methods can be attributed to the theory error inherent in any estimate of . The comparison presented in §4 shows that while all of these methods slightly disagree, the moment transport method gives good agreement with other established methods.
References
References
 [1] S. W. Hawking, The Development of Irregularities in a Single Bubble Inflationary Universe, Phys. Lett. B115 (1982) 295.
 [2] S. W. Hawking and I. G. Moss, Fluctuations in the inflationary universe, Nucl. Phys. B224 (1983) 180.
 [3] J. M. Bardeen, P. J. Steinhardt, and M. S. Turner, Spontaneous Creation of Almost Scale  Free Density Perturbations in an Inflationary Universe, Phys. Rev. D28 (1983) 679.
 [4] D. H. Lyth, Large Scale Energy Density Perturbations and Inflation, Phys. Rev. D31 (1985) 1792–1798.
 [5] A. H. Guth and S.Y. Pi, The Quantum Mechanics of the Scalar Field in the New Inflationary Universe, Phys. Rev. D32 (1985) 1899–1920.
 [6] T. Falk, R. Rangarajan, and M. Srednicki, The angular dependence of the threepoint correlation function of the cosmic microwave background radiation as predicted by inflationary cosmologies, Astrophys. J. 403 (1993) L1, [arXiv:astroph/9208001].
 [7] A. Gangui, F. Lucchin, S. Matarrese, and S. Mollerach, The threepoint correlation function of the cosmic microwave background in inflationary models, Astrophys. J. 430 (1994) 447–457, [arXiv:astroph/9312033].
 [8] T. Pyne and S. M. Carroll, HigherOrder Gravitational Perturbations of the Cosmic Microwave Background, Phys. Rev. D53 (1996) 2920–2929, [arXiv:astroph/9510041].
 [9] V. Acquaviva, N. Bartolo, S. Matarrese, and A. Riotto, Secondorder cosmological perturbations from inflation, Nucl. Phys. B667 (2003) 119–148, [arXiv:astroph/0209156].
 [10] J. M. Maldacena, NonGaussian features of primordial fluctuations in single field inflationary models, JHEP 05 (2003) 013, [arXiv:astroph/0210603].
 [11] D. H. Lyth and Y. Rodriguez, The inflationary prediction for primordial nongaussianity, Phys. Rev. Lett. 95 (2005) 121302, [arXiv:astroph/0504045].
 [12] D. Seery and J. E. Lidsey, Primordial nongaussianities from multiplefield inflation, JCAP 0509 (2005) 011, [arXiv:astroph/0506056].
 [13] E. Komatsu et al., NonGaussianity as a Probe of the Physics of the Primordial Universe and the Astrophysics of the Low Redshift Universe, arXiv:0902.4759.
 [14] E. Komatsu and D. N. Spergel, Acoustic signatures in the primary microwave background bispectrum, Phys. Rev. D63 (2001) 063002, [arXiv:astroph/0005036].
 [15] R. Juszkiewicz, D. H. Weinberg, P. Amsterdamski, M. Chodorowski, and F. Bouchet, Weakly nonlinear Gaussian fluctuations and the Edgeworth expansion, Astrophys. J. 442 (1995) 39.
 [16] F. R. Bouchet and R. Juszkiewicz, Perturbation theory confronts observations: Implications for the ‘initial’ conditions and , arXiv:astroph/9312007.
 [17] F. R. Bouchet, Introductory overview of Eulerian and Lagrangian perturbation theories, arXiv:astroph/9603013.
 [18] P. Fosalba, E. Gaztanaga, and E. Elizalde, Gravitational Evolution of the LargeScale Density Distribution: The Edgeworth & Gamma Expansions, arXiv:astroph/9910308.
 [19] F. Vernizzi and D. Wands, NonGaussianities in twofield inflation, JCAP 0605 (2006) 019, [arXiv:astroph/0603799].
 [20] T. Battefeld and R. Easther, Nongaussianities in multifield inflation, JCAP 0703 (2007) 020, [arXiv:astroph/0610296].
 [21] C. T. Byrnes, K.Y. Choi, and L. M. H. Hall, Conditions for large nonGaussianity in twofield slowroll inflation, JCAP 0810 (2008) 008, [arXiv:0807.1101].
 [22] C. T. Byrnes, K.Y. Choi, and L. M. H. Hall, Large nonGaussianity from twocomponent hybrid inflation, JCAP 0902 (2009) 017, [arXiv:0812.0807].
 [23] C. T. Byrnes and G. Tasinato, NonGaussianity beyond slow roll in multifield inflation, arXiv:0906.0767.
 [24] D. Battefeld and T. Battefeld, On NonGaussianities in MultiField Inflation ( fields): Bi and Trispectra beyond SlowRoll, arXiv:0908.4269.
 [25] D. Seery and J. E. Lidsey, Nongaussianity from the inflationary trispectrum, JCAP 0701 (2007) 008, [arXiv:astroph/0611034].
 [26] C. T. Byrnes, M. Sasaki, and D. Wands, The primordial trispectrum from inflation, Phys. Rev. D74 (2006) 123519, [arXiv:astroph/0611075].
 [27] S. Dimopoulos, S. Kachru, J. McGreevy, and J. G. Wacker, Nflation, JCAP 0808 (2008) 003, [arXiv:hepth/0507205].
 [28] J.L. Lehners and S. RenauxPetel, Multifield Cosmological Perturbations at Third Order and the Ekpyrotic Trispectrum, arXiv:0906.0530.
 [29] P. Creminelli, On nongaussianities in singlefield inflation, JCAP 0310 (2003) 003, [arXiv:astroph/0306122].
 [30] M. Zaldarriaga, NonGaussianities in models with a varying inflaton decay rate, Phys. Rev. D69 (2004) 043508, [arXiv:astroph/0306006].
 [31] N. ArkaniHamed, P. Creminelli, S. Mukohyama, and M. Zaldarriaga, Ghost Inflation, JCAP 0404 (2004) 001, [arXiv:hepth/0312100].
 [32] M. Alishahiha, E. Silverstein, and D. Tong, DBI in the sky, Phys. Rev. D70 (2004) 123505, [arXiv:hepth/0404084].
 [33] A. A. Starobinsky, Multicomponent de Sitter (Inflationary) Stages and the Generation of Perturbations, JETP Lett. 42 (1985) 152–155.
 [34] M. Sasaki and E. D. Stewart, A General analytic formula for the spectral index of the density perturbations produced during inflation, Prog. Theor. Phys. 95 (1996) 71–78, [arXiv:astroph/9507001].
 [35] G. I. Rigopoulos and E. P. S. Shellard, Nonlinear inflationary perturbations, JCAP 0510 (2005) 006, [arXiv:astroph/0405185].
 [36] G. I. Rigopoulos, E. P. S. Shellard, and B. J. W. van Tent, Nonlinear perturbations in multiplefield inflation, Phys. Rev. D73 (2006) 083521, [arXiv:astroph/0504508].
 [37] D. Langlois and F. Vernizzi, Evolution of nonlinear cosmological perturbations, Phys. Rev. Lett. 95 (2005) 091303, [arXiv:astroph/0503416].
 [38] D. Langlois and F. Vernizzi, Conserved nonlinear quantities in cosmology, Phys. Rev. D72 (2005) 103501, [arXiv:astroph/0509078].
 [39] D. Langlois and F. Vernizzi, Nonlinear perturbations of cosmological scalar fields, JCAP 0702 (2007) 017, [arXiv:astroph/0610064].
 [40] I. Huston and K. A. Malik, Numerical calculation of second order perturbations, arXiv:0907.2917.
 [41] D. H. Lyth, The curvature perturbation in a box, JCAP 0712 (2007) 016, [arXiv:0707.0361].
 [42] M. S. Sloth, On the one loop corrections to inflation and the CMB anisotropies, Nucl. Phys. B748 (2006) 149–169, [arXiv:astroph/0604488].
 [43] M. S. Sloth, On the one loop corrections to inflation. II: The consistency relation, Nucl. Phys. B775 (2007) 78–94, [arXiv:hepth/0612138].
 [44] D. Seery, Oneloop corrections to the curvature perturbation from inflation, JCAP 0802 (2008) 006, [arXiv:0707.3378].
 [45] H. R. S. Cogollo, Y. Rodríguez, and C. A. ValenzuelaToledo, On the Issue of the Series Convergence and Loop Corrections in the Generation of Observable Primordial NonGaussianity in SlowRoll Inflation. Part I: the Bispectrum, JCAP 0808 (2008) 029, [arXiv:0806.1546].
 [46] Y. Rodríguez and C. A. ValenzuelaToledo, On the Issue of the Series Convergence and Loop Corrections in the Generation of Observable Primordial NonGaussianity in SlowRoll Inflation. Part II: the Trispectrum, arXiv:0811.4092.
 [47] D. Seery, M. S. Sloth, and F. Vernizzi, Inflationary trispectrum from graviton exchange, JCAP 0903 (2009) 018, [arXiv:0811.3934].
 [48] D. Seery, K. A. Malik, and D. H. Lyth, Nongaussianity of inflationary field perturbations from the field equation, JCAP 0803 (2008) 014, [arXiv:0802.0588].
 [49] K.Y. Choi, L. M. H. Hall, and C. van de Bruck, Spectral running and nonGaussianity from slowroll inflation in generalised twofield models, JCAP 0702 (2007) 029, [arXiv:astroph/0701247].
 [50] S. Yokoyama, T. Suyama, and T. Tanaka, Primordial NonGaussianity in MultiScalar Inflation, Phys. Rev. D77 (2008) 083511, [arXiv:0711.2920].
 [51] J. GarciaBellido and D. Wands, Metric perturbations in twofield inflation, Phys. Rev. D53 (1996) 5437–5445, [arXiv:astroph/9511029].
 [52] C. Gordon, D. Wands, B. A. Bassett, and R. Maartens, Adiabatic and entropy perturbations from inflation, Phys. Rev. D63 (2001) 023506, [arXiv:astroph/0009131].
 [53] A. Taylor and P. Watts, Evolution of the cosmological density distribution function, arXiv:astroph/0001118.
 [54] L. Amendola, Non Gaussian likelihood function and COBE data, Mon. Not. Roy. Astron. Soc. 283 (1996) 983–989.
 [55] S. Matarrese, L. Verde, and R. Jimenez, The abundance of highredshift objects as a probe of non Gaussian initial conditions, Astrophys. J. 541 (2000) 10, [arXiv:astroph/0001366].
 [56] L. Amendola, The dependence of cosmological parameters estimated from the microwave background on nongaussianity, Astrophys. J. 569 (2002) 595–599, [arXiv:astroph/0107527].
 [57] P. Watts and P. Coles, Statistical Cosmology with Quadratic Density Fields, Mon. Not. Roy. Astron. Soc. 338 (2003) 806, [arXiv:astroph/0208295].
 [58] M. LoVerde, A. Miller, S. Shandera, and L. Verde, Effects of ScaleDependent NonGaussianity on Cosmological Structures, JCAP 0804 (2008) 014, [arXiv:0711.4126].
 [59] D. Seery and J. C. Hidalgo, NonGaussian corrections to the probability distribution of the curvature perturbation from inflation, JCAP 0607 (2006) 008, [arXiv:astroph/0604579].
 [60] T. Y. Lam and R. K. Sheth, Halo abundances in the model, arXiv:0905.1702.
 [61] S. Blinnikov and R. Moessner, Expansions for nearly Gaussian distributions, Astron. Astrophys. Suppl. Ser. 130 (1998) 193–205, [arXiv:astroph/9711239].
 [62] A. A. Starobinsky, Stochastic de Sitter (inflationary) stage in the early universe, . In De Vega, H.J. (Ed.), Sanchez, N. (Ed.): Field Theory, Quantum Gravity and Strings, 107126.
 [63] G. I. Rigopoulos, E. P. S. Shellard, and B. J. W. van Tent, Quantitative bispectra from multifield inflation, Phys. Rev. D76 (2007) 083512, [arXiv:astroph/0511041].
 [64] M. Sasaki, J. Valiviita, and D. Wands, Nongaussianity of the primordial perturbation in the curvaton model, Phys. Rev. D74 (2006) 103003, [arXiv:astroph/0607627].
 [65] V. Desjacques and U. Seljak, Signature of primordial nonGaussianity of type in the mass function and bias of dark matter haloes, arXiv:0907.2257.
 [66] M. GellMann and F. E. Low, Quantum electrodynamics at small distances, Phys. Rev. 95 (1954) 1300–1312.