Beyond \delta N formalism

We develop a theory of nonlinear cosmological perturbations on superhorizon scales for a multi-component scalar field with a general kinetic term and a general form of the potential in the context of inflationary cosmology. We employ the ADM formalism and the spatial gradient expansion approach, characterised by O(\epsilon^2), where \epsilon=1/(HL) is a small parameter representing the ratio of the Hubble radius to the characteristic length scale L of perturbations. We provide a formalism to obtain the solution in the multi-field case. This formalism can be applied to the superhorizon evolution of a primordial non-Gaussianity beyond the so-called \delta N formalism which is equivalent to O(\epsilon^0) of the gradient expansion. In doing so, we also derive fully nonlinear gauge transformation rules valid through O(\epsilon^2). These fully nonlinear gauge transformation rules can be used to derive the solution in a desired gauge from the one in a gauge where computations are much simpler. As a demonstration, we consider an analytically solvable model and construct the solution explicitly.


I. INTRODUCTION
Recent observations of the cosmic microwave background anisotropy show a very good agreement of the observational data with the predictions of conventional, single-field slow-roll models of inflation, that is, adiabatic Gaussian random primordial fluctuations with an almost scale-invariant spectrum [1]. Nevertheless, possible non-Gaussianities from inflation has been a focus of much attention in recent years, mainly driven by recent advances in cosmological observations. In particular, the PLANCK satellite [2] is expected to bring us preciser data and it is hoped that a small but finite primordial non-Gaussianity may actually be detected.
To study possible origins of non-Gaussianity, one must go beyond the linear perturbation theory. An observationally detectable level of non-Gaussianity cannot be produced in the conventional, single-field slow-roll models of inflation, since the predicted magnitude is extremely small, suppressed by the slow-roll parameters. Then a variety of ways to generate a large non-Gaussianity have been proposed. (See e.g. a focus section in CQG [3] and references therein for recent developments.) They may be roughly classified into two; multi-field models where non-Gaussianity can be produced classically on superhorizon scales, and non-canonical kinetic term models where non-Gaussianity can be produced quantum mechanically on subhorizon scales. In particular, in the former case, the δN formalism turned out to be a powerful tool for computing non-Gaussianities thanks to its full non-linear nature.
On the superhorizon scales, one can employ the spatial gradient expansion approach . It is characterised by an expansion parameter, ǫ = 1/(HL), representing the ratio of the Hubble radius to the characteristic length scale L of the perturbation. In the context of inflation, based on the leading order in gradient expansion, the δN formalism [6,11,12] or the separate universe approach [16] was developed. It is valid when local values of the inflaton field at each local point (averaged over each horizon-size region) determine the evolution of the universe at each point. This leading order in the gradient expansion provides a general conclusion for the evolution on superhorizon scales that the adiabatic growing mode is conserved on the comoving hypersurface [17].
In this paper, we consider the curvature perturbation on superhorizon scales up through next-to-leading order in gradient expansion, that is, to O(ǫ 2 ). To make our analysis as general as possible, we extend the δN formalism in the following two aspects: One is to go beyond the single-field assumption, and the other is to go beyond the slow-roll condition. While in the case of single-field inflation, the curvature perturbation remains constant as mentioned above, the superhorizon curvature perturbation can change in time in the case of multi-field inflation. Furthermore, even for single-field inflation, the time evolution can be non-negligible due to a temporal violation of the slow-roll condition. In order to study such a case, the δN formalism is not sufficient since the decaying mode cannot be neglected any longer, which usually appears at O(ǫ 2 ) of gradient expansion and is known to play a crucial role already at the level of linear perturbation theory [31,32], 1 not to mention the case of nonlinear perturbation theory [20,21,28,29].
Multi-field inflation may be motivated in the context of supergravity since it suggests the existence of many flat directions in the scalar field potential. In multi-field inflation, a non-slow-roll stage may appear when there is a change in the dominating component of the scalar field. For example, one can consider a double inflation model in which a heavier component dominates the first stage of inflation but damps out when the Hubble parameter becomes smaller than the mass, while a lighter component is negligible at the first stage but dominates the second stage of inflation after the heavier component has decayed out [24,34,35]. However, these previous analyses are essentially based on the δN formalism and it is in general necessary to extend it to O(ǫ 2 ), that is, to the beyond δN formalism. We focus on the case of a multi-component scalar field. As for a single scalar field, it has been developed in [28].
We mention that a multi-scalar case in the gradient expansion approach was studied previously [25,26]. However, it turns out to be valid only for a restricted situation (discussed later). Here we develop a general framework for fully nonlinear perturbations and present a formalism for obtaining the solution to O(ǫ 2 ). Then as an example we consider a specific model which allows an analytical treatment of the equations of motion.
This paper is organised as follows. In Sec. II, we introduce a multi-component scalar field and derive basic equations. We compare several typical time-slicing conditions and mention the differences of them from the single-field case. In Sec. III, we develop a theory of nonlinear cosmological perturbations on superhorizon scales. We formulate it on the uniform e-folding slicing (which is defined later). In Sec. IV, as a demonstration of our formalism, we consider an analytically solvable model and give the solution explicitly. Sec. V is devoted to a summary and discussions. Some details are deferred to Appendices. In Appendix A, the coincidence between some of time slicing conditions is discussed by using the Einstein equations. In Appendix B, we write down the basic equations on the uniform expansion slicing, and study the behaviour of the curvature perturbation in this slicing. In Appendix C, we give general nonlinear gauge transformation rules valid to next-to-leading order in gradient expansion. In Appendix D, we verify our formalism in a single-field model. Finally, in Appendix E, we discuss the structure of the Hamiltonian and momentum constraint equations in the gradient expansion.

A. The Einstein equations
We develop a theory of nonlinear cosmological perturbations on superhorizon scales. For this purpose we employ the ADM formalism and the gradient expansion approach. In the ADM decomposition, the metric is expressed as where α is the lapse function, β i is the shift vector and Latin indices run over 1,2 and 3. We introduce the extrinsic curvature K ij defined by whereD is the covariant derivative with respect to the spatial metricγ ij . In addition to the standard ADM decomposition, the spatial metric and the extrinsic curvature are further decomposed so as to separate trace and trace-free partsγ where a(t) is the scale factor of a fiducial homogeneous Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime and the determinant of γ ij is normalised to be unity and A ij is trace free. The explicit form of K is given by where H is the Hubble parameter defined by H(t) ≡ da(t) dt a(t). As for a matter field, let us focus on a minimally-coupled multi-component scalar field, where I, J and K run over 1, 2, · · · , M with φ K denoting the K-th component of the scalar field. Note that we do not assume a specific form of both the kinetic term and potential, which are arbitrary functions of X IJ and φ K . This type of Lagrangian can be applied to, for example, multi-field DBI inflationary models. For the calculation of their non-Gaussianities, see, e.g. [36,37] and also [38,39] for recent developments. The equation of motion for the scalar field is given by where the subscript I in P I represents a derivative with respect to φ I and P (IJ) is defined as The energy-momentum tensor is Notice that this energy-momentum tensor cannot be written in the perfect fluid form any more, which is one of main differences from the single-field case.
All the independent components of the energy-momentum tensor are conveniently expressed in terms of E and J i as (2. 10) and T ij , where n µ is the unit vector normal to the time constant surfaces and is given by For convenience, we further decompose T ij in the same way as Eq. (2.4), Now we write down the Einstein equations. In the ADM decomposition, the Einstein equations are separated into four constraints, the Hamiltonian constraint and three momentum constraints, and six dynamical equations for the spatial metric. The constraints are where R ≡ R[γ] is the Ricci scalar of the normalised spatial metric γ ij , D i is the covariant derivative with respect to γ ij , D 2 ≡ γ ij D i D j , γ ij is the inverse of γ ij , and the spatial indices are raised or lowered by γ ij and γ ij , respectively. As for the dynamical equations for the spatial metric, we rewrite Eq. (2.2) as The equations for the extrinsic curvature (K,A ij ) are given by where ∂ ⊥ ≡ n µ ∂ µ , and we have introduced the trace-free projection operator [...] T F defined for a tensor Q ij as Finally, the equations of motion for the scalar field (2.7) are (2.20)

B. Gradient expansion and assumption
In the gradient expansion approach we suppose that the characteristic length scale L of a perturbation is longer than the Hubble length scale 1/H of the background, i.e. HL ≫ 1. Therefore, ǫ ≡ 1/(HL) is regarded as a small parameter and we can systematically expand equations in the order of ǫ, identifying a spatial derivative is of order ǫ, To clarify the order of gradient expansion, we introduce the superscript (n). For example, (2) α means the lapse function at second order in gradient expansion.
As a background spacetime, we consider a FLRW universe. At O(ǫ 0 ) of the gradient expansion, there is apparently no spatial gradient and the universe is locally homogeneous and isotropic. This leads to the following condition on the spatial metric: (2.21) Since we adopt this assumption, the spatial metric at leading order is given by an arbitrary spatial function of the spatial coordinates, 22) under the condition that the eigenvalues of f ij are all positive definite everywhere. From the definition of A ij , Eq. (2.16), the above assumption implies (2.23) Throughout this paper, in order to simplify the equations, we set the shift vector to zero up to second order in gradient expansion, (2.24) Let us call this choice of the spatial coordinates as the time-slice-orthogonal threading. Here we mention that the above condition does not completely fix the spatial coordinates. As discussed later, one can actually make an arbitrary coordinate transformation of the form, x i →x i = f i (x j ).

C. Leading order in gradient expansion
In this subsection, we study the leading order gradient expansion and make clear the correspondence between the leading order equations and background equations. This correspondence can be used to construct the solution at leading order in gradient expansion in terms of the background solution.
At leading order in gradient expansion, the Einstein equations are (2.25) 27) and the scalar field equation is where we have introduced the proper time τ by In terms of τ , the expression of K in Eq. (2.5) is simplified under the time-slice-orthogonal threading condition, Under the identifications, ae ψ ⇔ a , and τ ⇔ t , (2. 31) one also has the correspondence, K ⇔ 3H. This means that the basic equations at leading order, Eqs. (2.27) and (2.28), take exactly the same form as those in the background modulo above identifications. Namely, given a background solution, one can construct the solution at leading order in gradient expansion as All the information of inhomogeneities is contained in the initial condition as well as in the proper time τ through Eq. (2.29). Thus it is sufficient to solve the background equations to obtain the solution at leading order in gradient expansion.
In passing, we note that the e-folding number is often used as the time coordinate to describe the background evolution. For convenience, we define it as the number of e-folds counted backward in time from a fixed final time. That is, (2.34) Accordingly, the scale factor is expressed as By replacing t with τ and H with K/3 we can generalise the e-fold number to the one defined locally in space as . One needs to specify the gauge condition to study perturbations in perturbation theory or in gradient expansion. Since spatial coordinates have been already fixed by the time-slice-orthogonal threading, one has to determine the time-slicing condition. Here, let us list various slicings and their definitions, Uniform e-folding number ; N (t, x i ) = N (t) . (2.41) Hereinafter, we call the uniform expansion, uniform energy and uniform e-folding number slicings as the uniform K, uniform E and uniform N slicings, respectively. We mention that there is a remaining gauge degree of freedom in the synchronous or uniform N slicing, while the time slices are completely fixed in the uniform K and uniform E slicings. As for the uniform N slicing, the gauge condition demands ∂ t ψ to vanish from Eqs. (2.30) and (2.36). This means one can freely choose the initial value of ψ (and hence its spatial configuration at any later time because ψ is conserved). This corresponds to the freedom in the choice of the initial time-slice as we see later. Utilising this freedom, we can make a scalar quantity, one of scalar fields φ I or K for example, homogeneous on the initial slice.
E. Towards the next-to-leading order in gradient expansion As we have seen in subsection II C, the leading order solutions are given by functions of τ in terms of the background solutions. At next-to-leading order in gradient expansion, terms with spatial derivatives of the leading order solution appear in the evolution equations. To evaluate those terms, one needs to calculate the spatial derivative of the lapse function, for example in ∂ i φ, (2.42) However the leading order (0) α is in general given explicitly only after solving the following equation for α: As a demonstration, the analysis on the uniform K slices is performed in Appendix B, and it is clearly shown that it is almost impossible to solve this equation, at least in an analytical way. This problem did not appear in the single-field case. It is because one can show that various different slicings become identical at leading order in gradient expansion. In particular, all the slicings listed in subsection II D coincide with each other as shown in Appendix A: where → means the left ones imply the right ones and means it holds when one chooses the initial slice to be the comoving, uniform K or uniform E slicing by using the remaining gauge degree of freedom. Thus the lapse function is homogeneous in all the slicings in the above, and we may set α = 1 if desired.
On the other hand, one has to face this problem in the case of multi-field inflation. We overcome this problem by choosing the synchronous slicing or uniform N slicing, which gives us a homogeneous time coordinate. On these slicings, one can evaluate the spatial derivatives of the leading order solution which appear as source terms and construct a solution to next-to-leading order in gradient expansion by integrating those terms.
There are two necessary steps before reaching the goal. Once we have a solution, it is necessary to construct a conserved quantity out of it that can be directly related to observable quantities. It is widely known that the comoving curvature perturbation eventually become conserved in a single-field model in linear theory. In non-linear theory, there exists a corresponding quantity, ψ on the comoving, uniform K or uniform E slicing, which is conserved at leading order in gradient expansion [17]. Even in the multi-field case, the system effectively reduces to a single-field system after the so-called non-adiabatic pressure has died out, that is, when the adiabatic limit is reached. Therefore it is necessary to perform a nonlinear gauge transformation from the uniform N slicing to one of those three slicings. This is one of the steps. Since the comoving slicing is not well defined in general in the multi-field case [12], we choose the uniform K slicing as the target gauge.
The other step to be taken is related to the definition of the curvature perturbation at next-to-leading order. In linear theory the curvature perturbation is named so because it determines the three-dimensional Ricci scalar, and ψ can be called so to full nonlinear order in the context of the leading order in gradient expansion. At next-to-leading order, however, ψ itself is no longer adequate to be called the curvature perturbation [28]. One needs to add the contribution from part of γ ij , which we call χ, to obtain a properly defined curvature perturbation conserved through O(ǫ 2 ). Therefore, after transforming from the uniform N slicing to the uniform K slicing, one has to evaluate the combination, R K ≡ ψ K + χ K /3. This is the Beyond δN formalism.
Before concluding this section, we mention the difference between our work and that of Weinberg [25,26]. There it was assumed that the lapse function can be chosen to be equal to unity at leading order in gradient expansion, hence all the scalar fields are homogeneous. This severely constrains the class of scalar field models as well as the initial condition because the curvature perturbation must be always conserved at leading order in gradient expansion. Here we do not impose such assumptions and perform a completely general analysis.

III. BEYOND δN FORMALISM
Let us first summarise the five steps in the Beyond δN formalism.

Write down the basic equations (the Einstein equations and scalar field equation) in the uniform N slicing
with the time-slice-orthogonal threading. For convenience let us call the choice of the coordinates in which one adopts the uniform X slicing with the time-slice-orthogonal threading the X gauge. So the above choice is the N gauge. In this gauge the metric components at leading order are trivial since both ψ and γ ij are independent of time.
2. First solve the leading order scalar field equation under an appropriate initial condition and then the next-toleading order scalar field equation which involves spatial gradients of the leading order solution.
3. Solve the next-to-leading order Einstein equations for the metric components and their derivatives.
4. Determine the gauge transformation from the N gauge to the K gauge and apply the gauge transformation rules to obtain the solution in the K gauge.

Evaluate the curvature perturbation
In what follows, we describe these steps in detail but only formally. An example in which these steps can be computed analytically will be discussed in Sec. IV.
Step 1: First, we rewrite the uniform N slicing condition from Eqs. (2.30) and (2.36) as Hence ψ is constant in time and is given by a function of the spatial coordinates alone, In the N gauge, the Einstein equations are reduced to the following equations. The constraints are The evolution equations for K, A ij and γ ij are The scalar field equation is We rewrite Eq. (3.8) by eliminating ∂ N K with Eq. (3.5) as Thus once K 2 is expressed in terms of the scalar field and its derivatives, the above equation gives a closed scalar field equation. An explicit derivation of the closed scalar field equation is possible only after we specify the explicit form of P (X IJ , φ) as a function of X IJ and φ. Here we describe generally but formally the procedure to obtain K 2 as a function of the scalar field. First, we separate the term with time-derivatives in X IJ and denote it by (3.10) At leading order, we can neglect the spatial-derivative term in Eq. (3.10). Then, P and P (IJ) are given by P (K 2 Y IJ , φ I ) and P (IJ) (K 2 Y IJ , φ I ), respectively. To next-to-leading order, P and P (IJ) can be expanded as Inserting these expressions into Eq. (3.3), one obtains an algebraic equation for K 2 . Solving it gives an expression of K 2 in terms of the scalar field. Then a closed equation for the scalar field is obtained by plugging it into Eq. (3.9).
Step 2: Although we can keep our discussion completely general, below we focus on the case of a multi-component canonical scalar field, This choice is taken purely for the sake of simplicity and clarity, because the expansion K can be explicitly expressed in terms of the scalar field in this case. In general, one cannot obtain an explicit expression of K in terms of the scalar field unless the form of P is explicitly specified. Nevertheless, the discussion below also applies to the general case perfectly.
From Eq. (3.3) we find (3.14) Inserting this into Eq. (3.9), one obtains the following closed equation: Since each term in the right-hand side of the equation involves two spatial derivatives, we can neglect them at leading order. At next-to-leading order, they can be understood as source terms whose time evolution have been already determined from the leading order solution. As noted in the above, although one cannot obtain an equation like the above explicitly for general P , one can always derive a closed equation for the scalar field once a specific form of P is given.
Solving the closed equation for the scalar field, the solution is formally obtained as where the subscript N indicates the N gauge. In the arguments on the right-hand side, the subscript 0 denotes the initial value, and D(· · · ) the spatial derivatives of the quantities inside the parentheses.
Step 3: Once the solution for the scalar field is obtained, one can solve Eqs. (3.5), (3.6) and (3.7) to obtain the metric quantities. As for K, however, it is simpler and in fact better to use Eq. (3.14), which is essentially the Hamiltonian constraint, since the integral constants appearing from integrating Eq. (3.5) are not freely specifiable but must satisfy the Hamiltonian constraint.
Step 4: Given the solution in the N gauge, the next step is to find the gauge transformation from the N gauge to the K gauge. It can be determined as follows. As noted above, the expression for K N in terms of the scalar field is obtained from Eq. (3.14). Since the leading order and next-to-leading order scalar field solutions are expressed as Eqs. (3.16), the same is true for K N , Let the transformation from the uniform N slicing to the uniform K slicing be given by N →Ñ = N + n(N, x i ) or converselyÑ +ñ(Ñ ,x i ) = N , where the uniform K slice is given byK =const.. This nonlinear gauge transformation is discussed in detail in Appendix C. The nonlinear gauge transformation generatorñ(Ñ , x i ) from the N gauge to K gauge is then determined by the condition that K is spatially homogeneous on the uniform K slice by definition: Then the other quantities in the K gauge are obtained by applying the above gauge transformation. In particular, the determinant of the spatial metric ψ K in the K gauge is obtained from Eqs. (C28) and (C33), and the unimodular part of the spatial metric γ ij K from Eqs. (C29) and (C34).
Step 5: We are to construct the nonlinear curvature perturbation, R K = ψ K + χ K /3, where χ K is the scalar part of the metric γ ij K . The scalar part of γ ij is defined in the same way as the single-scalar case [27], where △ −1 is the inverse Laplacian operator on the flat background. Extracting χ K from γ ij K , and combining ψ K and χ K , we finally obtain the nonlinear curvature perturbation R K = ψ K + χ K /3 in the K gauge.

IV. SOLVABLE EXAMPLE
In this section, we demonstrate how to obtain the solution up to next-to-leading order in gradient expansion by applying our formalism to a specific, analytically solvable model.

A. Model and equations
For simplicity, we consider a canonical scalar field with exponential potential [23], where W is a constant. The leading order scalar field equation is given by where we have omitted a summation symbol over the field component indices, J. Hereafter summation is implied over repeated component indices. Further we assume the two slow-roll conditions on the leading order equation. The first one is that we can neglect the "kinetic energy" in the energy density of the scalar field, The second one is that we can neglect the "acceleration", It is important that we apply these assumptions only to (0) φ. We do not impose the slow-roll conditions on (2) φ. So we can rewrite Eq. (4.2) as At next-to-leading order, we have where Inserting the leading order equation (4.5) into (4.6), it gives This is the basic equation at next-to-leading order.

B. Solution
We can easily solve the leading order and next-to-leading order scalar field equations. The solution of Eq. (4.5) is given by 9) and the solution of Eq. (4.8) is obtained as where N 0 is an initial time and C φ I and D φ I represent the initial values of the scalar field and its time derivative. Note that the solution (0) φ I satisfies the slow-roll conditions (4.3) and (4.4) if all the masses are small, Here let us show the time-independence of S φ I , which is given by Eq. (4.7). From Eq. (3.14) we have The leading order potential (0) V is given by where C V is the initial value of the potential, (4.14) Substituting the above solution into Eq. (4.7) gives The time-independence of S φ I is now clear since it is expressed solely in terms of C ψ , C V , C φ I and (0) γ ij which are all time-independent functions. Therefore we can rewrite the second order solution (4.10) in a simpler manner as where Given the solution for the scalar field up to second order in gradient expansion, we can now obtain those for K, A ij and γ ij . At leading order, K is expressed in terms of the scalar field through Eq. (4.12) as Because of the assumption (2.21), (0) γ ij and (0) A ij are trivial, At next-to-leading order, Eq. (3.14) becomes 20) and Eq. (3.6) reduces to where S A ij is independent of time. Specifically, Substituting the solution of the scalar field into Eq. (4.20), we obtain (2) K as where As for A ij , we obtain it by integrating Eq. (4.21), Finally, (2) γ ij is given by solving Eq. (3.7), Before concluding this subsection, let us fix the remaining gauge degree of freedom on the uniform N slicing mentioned in subsection II D. We can make (0) K spatially homogeneous on the initial slice by using this gauge degree of freedom, where (0) V 0 is a pure constant independent of both space and time. With this choice, we can substantially simplify the above expressions of the solution because the spatial derivative of C V vanishes. To summarise, thus obtained solution is

30)
where (0) V = (0) V 0 exp M 2 (N − N 0 ) and C φ I , D φ I , C γ ij and C A ij are initial values of φ, ∂ N φ, γ ij and A ij , respectively. The source functions S φ I and S A ij are time-independent whose explicit forms are 33) and the functions I φ , S K , I A , I γ and J γ are given by C. Solution on the uniform K slice In this subsection, we derive the solution on the K gauge by applying a gauge transformation to the solution on the N gauge. To do so, we first need to determine the generator of the gauge transformation between the two slices, N →Ñ = N + n(N, x i ) or converselyÑ +ñ(Ñ ,x i ) = N . The nonlinear gauge transformation is discussed in detail in Appendix C.
As is clear from Eq. (4.18) with the initial condition C V = (0) V 0 =const., Eq. (4.27), the uniformity of K on the uniform N slices is kept in this model at leading order if it is chosen so on the initial slice. Thus no gauge transformation is necessary at leading order.
At next-to-leading order, the gauge transformation of K is from Eq. (C35) given by Since (2)K vanishes on the uniform K slice, (2)ñ is determined from Eq. (4.20) as The solution for ψ is obtained from Eq. (C33) as where C ψ is given in (3.2) and that for γ ij is same as γ ij on the N gauge from Eq. (C34), At next-to-leading order in gradient expansion, we have to extract the scalar component χ from γ ij to obtain the curvature perturbation R. In linear theory, the variable χ reduces to the traceless scalar-type component H Lin T , and the curvature perturbation is given by R Lin = (H Lin L + H Lin T /3)Y , where ψ → H Lin L in the linear limit, and we have followed the notation of [40].
Neglecting the contribution of gravitational waves and focusing on that of the scalar-type perturbations, we can assume that (0) γ ij in the K gauge approaches the flat metric at late times when the adiabatic limit is reached, This condition completely fixes the remaining spatial gauge degrees of freedom, say, within the class of the time-sliceorthogonal threading. Under this condition, we rewrite γ ij on the K gauge given in Eq. (4.30) to manifest its time and spatial dependence, where C γ ij = δ ij from Eq. (4.43), and C χ ij is given by The time-dependent functions f A and f χ are given by where the functions I γ and J γ are defined in Eqs. (4.36) and (4.37), respectively. Now we extract the scalar part from C A ij and C χ ij by using the definition of χ given by Eq. (3.20). We note that this definition is unique in the sense that it reduces to the standard scalar part in the limit of linear theory and gives the O(ǫ 2 ) correction unambiguously.
The contribution of C A ij to χ may be determined by evaluating Eq. (3.4) on the initial slice with the help of Eqs. (4.28) and (4.31), Applying the definition of χ in Eq. (3.20) to the above, we find the contribution from C A ij is As for the contribution of C χ ij to χ, we can make use of the relation, where the left-hand side is the Ricci scalar of e 2C ψ δ ij . Thus the contribution from C χ ij is Adding both contributions together, we obtain where χ A and χ χ are given, respectively, by Eqs. (4.49) and (4.51).
Finally summing all the contributions to the curvature perturbation R = ψ K + χ K /3, where ψ K and χ K are given by Eqs. (4.41) and (4.52), respectively, we obtain What we need to know is the final value of R at sufficiently late times, N → 0 (a → a 0 e N0 ). We take N 0 to be a time around which the scales relevant to cosmological observations crossed the Hubble horizon, hence N 0 50. In this case, at N = 0, the curvature perturbation reduces to The first term, (0) C ψ represents the leading order curvature perturbation obtainable in the usual δN formalism, and the remaining terms represent O(ǫ 2 ) contributions, the calculation of which is the main purpose of the beyond δN formalism. As a check, we have confirmed that the above result is consistent with the one obtained in linear perturbation theory with the same background solution.
It should be noted that there is no contribution from χ K at N = 0 by definition. However, this does not mean the evaluation of χ K was meaningless. In general it can make an important contribution to R around the horizon crossing time N ∼ N 0 , and by matching our R with that evolved from inside the horizon at N ∼ N 0 , the final value R(N = 0) is determined both by the values of ψ K and χ K and by their derivatives at N ∼ N 0 .
It may be noted that the last term proportional to D φ I comes the time derivative of . Although this is generically small by construction, the contribution of the term itself can become large since it is proportional to m I /M 2 ∼ 1/M where M 2 ≪ 1 as assumed in Eq. (4.11).

V. SUMMARY AND DISCUSSION
In this paper, we developed a theory of nonlinear cosmological perturbations on superhorizon scales in the context of inflationary cosmology. We considered a multi-component scalar field with a general kinetic term and a general form of the potential. To discuss the superhorizon dynamics, we employed the ADM formalism and the spatial gradient expansion approach.
Different from the single-field case, there is a difficulty in solving the equations in the multi-field case. At leadingorder, the equations take the same form as those for the homogeneous and isotropic FLRW background with suitable identifications of variables. In particular, there are correspondences between the cosmic time in the background and the proper time along each comoving trajectory, τ ⇔ t, and the scale factor with the one defined locally ae ψ ⇔ a. This implies that given the background solution, the solution at leading-order in gradient expansion is automatically obtained.
In the single-field case, one can show the coincidence among the comoving, uniform expansion, and synchronous slicings at leading order. This allows us to set the lapse function to unity, and replace the proper time by the cosmic time in the solution for the scalar field. Then the metric is expressed in terms of the scalar field easily, and the nextto-leading order equations can be solved straightforwardly because the space-time dependence of the source terms consisting of the leading order quantities is explicitly known.
In cosmological perturbation theory, the most important quantity to be evaluated is the curvature perturbation on the comoving slices which is conserved on superhorizon scales after the universe has reached the adiabatic limit. This quantity accurate to next-to-leading order may be relatively easily obtained in the single-field case because of the above mentioned coincidence among the comoving, uniform expansion and synchronous slicings.
On the other hand, in the multi-field case, such a coincidence between different slicings does not hold. This implies the following. One can express the lapse function as a function of the scalar field, but the scalar field is also a function of the proper time. Thus one has the equation, To go beyond the leading order, one needs to solve this equation for α, but it seems almost impossible.
In this paper, we developed a formalism to go beyond the leading order which avoids the above problem. Namely, we first solve the field equations in a slicing in which the lapse function is trivial. The synchronous slicing is one of such slicings, but we adopt the uniform e-folding number slicing in which the time slices are chosen in such a way that the number of e-folds along each orbit orthogonal to the time slices, N , is spatially homogeneous on each time slice. In other words, the uniform N slicing is a synchronous slicing if N is used as the time coordinate.
In this slicing we can solve the equations to next-to-leading order without encountering the above mentioned problem. After the solution to next-to-leading order is obtained, we transform it to the one in the uniform expansion slicing which is known to be identical to the comoving slicing on superhorizon scales in the adiabatic limit. Thus the gauge transformation laws play an essential role in our formalism. We derived them which are accurate to next-toleading order. Note that they are fully nonlinear in nature in the language of the standard perturbation approach. They are summarised in Appendix C.
As a demonstration of our formalism, we considered an analytically solvable model and constructed the explicit form of the solution. Namely, we considered a multi-component canonical scalar field with exponential potential, Following the general procedure discussed above, we first solved for the scalar field and the metric in the uniform N slicing. Then using a remaining gauge degree of freedom of this slicing, we set the initial condition such that it coincides with the uniform expansion slicing at leading order. In this slicing with this initial condition, the next-to-leading order solution was straightforwardly found. Then we applied the derived nonlinear gauge transformation to it to obtain the solution in the uniform expansion slicing. Finally, from thus obtained spatial metric which takes the form, e ψ γ ij where γ ij is a unimodular metric, we constructed the generalised, nonlinear curvature perturbation R defined by By inspecting the form of the final value of the curvature perturbation, we argued that the decaying modes of the scalar field may give rise to non-negligible contribution to the final, conserved curvature perturbation. To quantify the effect, it is necessary to match our solution on super-horizon scales with the one solved from sub-horizon scales. The evaluation of these next-to-leading order corrections to the spectrum as well as to the bispectrum of the curvature perturbation is left for future study.
where τ is defined in terms of the lapse function as Here let us choose the uniform K slicing, Then Eq. (A2) immediately implies that the scalar field is uniform on this slice, This shows the uniform K slice coincides with the comoving slice. Now since the left-hand side of Eq. (A1) as well as the potential term on the right-hand side is homogeneous, it follows that the kinetic term is homogeneous. This means (A6) Therefore the uniform K slicing coincides with the synchronous slicing. Then the spatial homogeneity of both K and α implies that of the number of e-folds N by definition, Eq. (2.36). Thus we have shown the coincidence of these four slicings at leading order in gradient expansion. As we noticed in subsection II D, there is a remaining gauge degree of freedom in the synchronous and uniform N slicings. In the above proof, this freedom is tacitly fixed since we took the uniform K slicing as a starting point. In other words, the uniform K (or comoving) slicing implies it is synchronous and uniform N . However, if one starts from the synchronous or uniform N slicing, we have to set the initial data such that K or φ is uniform on the initial slice in order to show the coincidence. In fact, if we use this remaining degree of freedom to set a different condition on the initial data, for example ψ = 0, the synchronous or uniform N slicing will not coincide with the uniform K or comoving slicing.
In passing let us also comment on the conservation of the curvature perturbation. From the definition of K, Eq. (2.5), one has This shows how the inhomogeneous part on the right-hand side contributes to the evolution of the curvature perturbation. As we have discussed above, there is no such inhomogeneities in the uniform K or comoving slicing in the single-scalar case. This is a proof of the conservation of the nonlinear curvature perturbation on super-horizon scales.
On the other hand, in the case of multi-field, the corresponding scalar component on the right-hand side of the momentum constraint, Eq. (A2), becomes uniform on the uniform K slice. However this scalar component does not appear in Eq. (A1) in general, and hence one cannot show the homogeneity of the lapse function as in the single-field case.
We note that this does not mean that the Hamiltonian and momentum constraints are not related. Actually the leading order momentum constraint reduces to the spatial derivative of the Hamiltonian constraint in the slow-roll limit as shown in Appendix E 2 a. For a non-slow roll case, this relation between the Hamiltonian and momentum constraints holds only if we take into account the next-to-leading order corrections properly (see Appendix E 2 b).

Appendix B: Uniform K slicing
In this Appendix, we derive the basic equations on the uniform K slicing, K(t, x i ) = 3H(t), and discuss a difficulty in constructing the solution at next-to-leading order in gradient expansion.
In the uniform K slicing and time-slice-orthogonal threading, the constraint equations are given by The evolution equation for K and A ij are where and τ is defined in terms of the lapse function as and the equation for the scalar field is First, let us consider the leading order. The Hamiltonian constraint (B1) reduces to hence (0) E is homogeneous on this slice. This means the solution for (0) E exactly coincides with the background solution E BG . As for P , the leading order solution is easily constructed in terms of the background solution once it is known as discussed in Sec. II C. To be specific, if the background pressure is expressed as a function of t, the leading order solution for P is given by The lapse function can be obtained from Eq. (B3) as Thus the lapse function becomes inhomogeneous if (0) P is so. From Eq. (2.5), the evolution of ψ is related with the inhomogeneity of (0) α in this slicing as which reflects the well-known fact that the curvature perturbation evolves if a non-adiabatic pressure exists. However there is a profound problem in Eq. (B11). Although it looks like Eq. (B11) gives the solution for (0) α, it is not so because the right-hand side depends also on (0) α. In reality, Eq. (B10) gives (0) P as a function of τ , which is given by the time integration of the lapse function. Therefore Eq (B11) can be schematically expressed as In general it is very difficult to solve this equation, at least analytically. As clear from Eq. (B3), one needs to evaluate the spatial derivatives of α for example at next-to-leading order. But this will be almost impossible. This is the main problem in constructing a solution at next-to-leading order in gradient expansion on the uniform K slicing. This kind of problem arises also on other slicings except for the synchronous or uniform N slicings. This is why the uniform N slicing is used in our formalism.
As a closing remark of this Appendix, we emphasise that the above difficulty becomes actually problematic only at next-to-leading order. It is not a problem at leading order. This is because the solution of the "curvature perturbation" (0) ψ is given by a time integration which depends only on the initial and final values of N . From Eq. (B12) we have Thus, at leading order in gradient expansion, as far as we are interested in the curvature perturbation there is no need to know the solution for (0) α. This is why the δN formalism works well despite the presence of the above problem.
This gives the leading order gauge transformation of ψ as As is clear from Eq. (C6), γ ij is invariant at this order, Under the gauge transformation, the scalar field is essentially invariant except for the change of the arguments, Then, the leading order gauge transformation is Here we note that what we call a gauge transformation should perhaps be called a coordinate transformation as far as the change of time slicing is concerned. For example, for a scalar quantity Q, the gauge transformation under the temporal shift (t → t + T ) is given by which is derived from the invariance of Q as a scalar under the time coordinate transformation, In the context of gradient expansion, we define the nonlinear gauge transformation of a scalar quantity by Eq. (C16) not by Eq. (C15), since the gradient expansion is a completely nonlinear approach. If we were to expand perturbatively the left-hand side of Eq. (C16), we would obtain an infinite number of terms, which is definitely something to be avoided. On the other hand, as we see shortly below, to the accuracy of our interest, i.e. to O(ǫ 2 ), the spatial coordinate transformation can be linearised. Hence we interpret it as the standard gauge transformation and compare the quantities at the same coordinate values. Once we have the gauge transformation rules for the metric components, we can readily obtain those for the extrinsic curvature components, K and A ij . Since A ij vanishes at leading order, we only have to consider the transformation of the expansion K. From its definition, Eq. (2.5), and using the transformation rules for the metric components derived above, we find it is invariant at leading order, Again, as discussed around Eq. (C16), we note that K is invariant only in the sense of the nonlinear gauge transformation we have defined.

Next-to-leading order gauge transformation
Now we derive the gauge transformation at next-to-leading order. To do so, we first have to fix the spatial coordinate transformation. It is determined by Eq. (C5), which results from the requirement that the shift vector should vanish to O(ǫ 3 ). In terms ofñ, the gauge transformation generatorL i is given as where l i is a time-independent spatial vector which represents the remaining gauge degrees of freedom in the spatial coordinates in the time-slice-orthogonal threading.
To consider the gauge transformation at next-to-leading order, we expand the nonlinear gauge transformation generatorñ in Eq. (C1) to the leading order term, the next-to-leading order term and so on, n ≡ (0)ñ + (2)ñ + O(ǫ 4 ) . (C19) As for the spatial shiftL i , it is unnecessary to expand it, since it is O(ǫ) already and the next-to-leading order term inL i is O(ǫ 3 ) which only affects the gauge transformation at O(ǫ 4 ). First we consider the metric components. From Eq. (C4), the next-to-leading order gauge transformation for the lapse function is determined as where ∂ĩ = ∂/∂x i and (0)Ñ (N, Note that all the quantities on the right-hand side are those evaluated at N = (0)Ñ and x i =x i . Taking the determinant of Eq. (C6), we obtain the gauge transformation for (2) ψ, Then using Eq. (C6) again, the gauge transformation for γ ij is obtained as Next we consider the extrinsic curvature. At next-to-leading order, Eq. (2.5) gives And also Eq. (2.16) determines the gauge transformation for A ij as Finally, the next-to-leading order gauge transformation for the scalar field is given from Eq. (C13) as To conclude this Appendix, let us summarise the nonlinear gauge transformation derived above. The leading order transformation rules areα The next-to-leading order transformation rules are (2)γ and whereL i is given byL 1. Solution in the uniform N gauge First we construct the solution in the N gauge by transforming it from the one in the comoving gauge. As we have seen in Appendix A, the leading order solution in the comoving gauge is where ∂ĩ = ∂/∂x i andL i is given explicitly bỹ After some manipulations, Eqs. (D9), (D11) and (D12) may be re-expressed as where Dĩ is a covariant derivative with respect to (0) γ ijN (or equivalently to C ij ). These expressions are convenient for later use.

Consistency check for Aij
Here we explicitly show that the solution for A ij obtained in the previous subsection satisfies the evolution equation in the N gauge. This verifies the consistency of the derived gauge transformation.
In the N gauge, the evolution equations for A ij is while in the comoving gauge it is We rewrite the second line inside the brackets of Eq. (D18) as where we have used the leading order Einstein equation, By subtracting Eq. (D19) from Eq. (D18), we find with the aide of Eq. (D20), where Q N c = Q N − Q c . On the other hand, taking the time derivative of Eq. (D17), which we have obtained by the gauge transformation, we obtain Comparing Eqs. (D22) and (D23), we see the precise coincidence between the two.

Consistency check for K
In the N gauge, the evolution equation for K is and that in the comoving gauge is Again we rewrite the second line inside the brackets of Eq. (D24) as where we have used Eq. (D21). By subtracting Eq. (D25) from Eq. (D24), we find with the aide of Eqs. (D15) and (D26), On the other hand, taking the time derivative of Eq. (D16), which has been obtained by the gauge transformation, we obtain ∂Ñ K We see that Eqs. (D27) and (D28) coincide precisely with each other.
where the suffix tmp indicates that it is obtained from the temporal evolution equation. Now let us carefully examine the solution (E3) in comparison with that obtained in the text, Eq. (4.29). First we compare the initial values. This determines the initial value of (2) K tmp as (2) K tmp (N 0 ) Next we subtract the solution (E3) together with thus obtained initial data (E4) from the solution (4.29), After integration parts and some manipulations, we obtain The right-hand side of the above equation must vanish at all times. This demands that the term inside the square brackets should vanish, Apparently this gives a non-trivial constraint among initial data, C ψ and C φ I . This is the "hidden" constraint corresponding to the Hamiltonian constraint.

Consistency with the momentum constraint
Below we show that the momentum constraint is automatically satisfied if the Hamiltonian constraint as well as the scalar field equation is satisfied. In the case of slow-roll inflation, the leading order Hamiltonian constraint is sufficient to show it. For general (non-slow-roll) inflation, one needs the next-leading order Hamiltonian constraint.

a. Slow-roll case
Under the slow-roll approximation, the Hamiltonian and momentum constraints reduce respectively to and the scalar field equation becomes Now we multiply both sides of the momentum constraint (E9) by 2K, and substitute the scalar field equation (E10) into the right-hand side of it. This gives (left − handside) = 4 3 One can easily see that these are merely a spatial derivative of the Hamiltonian constraint (E8). Thus we conclude that, at leading order in gradient expansion, the momentum constraint is automatically satisfied if the Hamiltonian constraint and the scalar field equations are satisfied.

b. Non slow-roll case
Here we show that in general the leading order momentum constraint is automatically satisfied once the Hamiltonian constraint and the energy conservation equations are satisfied to next-leading order in gradient expansion.
From Eqs. (2.13) and (2.14) with A ij = O(ǫ 2 ), the Hamiltonian and momentum constraint equations are 1 a 2 e 2ψ R − 4D 2 ψ + 2D i ψD i ψ + The energy conservation law, n ν ∇ ν T µν = 0, in the N gauge is Now let us multiply the evolution equation for K, Eq. (3.5), by 2K/3 and subtract Eq. (E15) from it. This gives Substituting the Hamiltonian constraint (E13) into the left-hand side of the above equation, we find ∂ N 1 3 K 2 − E = −∂ N 1 2a 2 e 2ψ R − 4D 2 ψ + 2D i ψD i ψ = − 1 a 2 e 2ψ R − 4D 2 ψ + 2D i ψD i ψ . (E17) Comparing the above two equations, we obtain Again one can see that the leading order momentum constraint (E14) holds automatically if the next-to-leading order Hamiltonian constraint (E13) holds. What does this mean? In the gradient expansion, we have assumed A ij does not contribute to the leading order dynamics. This corresponds to neglecting the adiabatic decaying mode in linear theory. In a single field model, it is well known this decaying mode appears only at the next-leading order both in linear theory [31] and in non-linear theory [28,29]. We also mention the work [41] in which they discussed the behavior of decaying modes in different choices of gauge. The absence of the decaying mode at leading order in gradient expansion should hold also in the case of multi-field inflation, at least as long as the background homogeneous solution is stable against a homogeneous but anisotropic perturbation.
To summarise, the point is that the initial data for the scalar field and its time derivatives are not freely specifiable in general, and in particular we need to take into account the next-leading order terms in the Hamiltonian constraint for the non-slow-roll case. Once we take into account the next-leading order Hamiltonian constraint, the leading order momentum constraint is automatically satisfied. As for the next-to-leading order momentum constraint, it constrains the initial value of A ij as in Eq. (4.47).