Landscape tomography through primordial non-Gaussianity

In this paper, we show how the structure of the landscape potential of the primordial Universe may be probed through the properties of the primordial density perturbations responsible for the origin of the cosmic microwave background anisotropies and the large-scale structure of our Universe. Isocurvature fields -fields orthogonal to the inflationary trajectory- may have fluctuated across the barriers separating local minima of the landscape potential during inflation. We analyze how this process could have impacted the evolution of the primordial curvature perturbations. If the typical distance separating consecutive minima of the landscape potential and the height of the potential barriers are smaller than the Hubble expansion rate parametrizing inflation, the probability distribution function of isocurvature fields becomes non-Gaussian due to the appearance of bumps and dips associated to the structure of the potential. We show that this non-Gaussianity can be transferred to the statistics of primordial curvature perturbations if the isocurvature fields are coupled to the curvature perturbations. The type of non-Gaussian structure that emerges in the distribution of curvature perturbations cannot be fully probed with the standard methods of polyspectra; instead, the probability distribution function is needed. The latter is obtained by summing all the $n$-point correlation functions. To substantiate our claims, we offer a concrete model consisting of an axionlike isocurvature perturbation with a sinusoidal potential and a linear derivative coupling between the isocurvature and curvature fields. In this model, the probability distribution function of the curvature perturbations consists of a Gaussian function with small superimposed oscillations reflecting the isocurvature axion potential.


I. INTRODUCTION
The search for primordial non-Gaussianity (NG) has been guided by our ability to make predictions within the inflationary paradigm [1][2][3][4][5]. The simplest models of inflation predict that the main departures from Gaussianity are to be found in the form of small, but nonvanishing, three-point correlation functions of the primordial comoving curvature perturbation ζ [6][7][8][9]. However, interactions involving the inflaton and other fields could enhance the amplitude of the three-point or higher point correlation functions (see [10][11][12][13] for reviews). These interactions could be self-interactions of the inflaton or interactions of the inflaton with other degrees of freedom. While current cosmic microwave background (CMB) constraints on the bispectrum are consistent with Gaussian statistics [14], it is possible that the method of threepoint or higher-point correlation functions do not constitute the most efficient parametrization of primordial NG hidden in the data [15,16].
Multifield  and quasi-single-field  models of inflation constitute a particularly useful framework to study the generation of large primordial NG. The interaction between curvature perturbations and isocurvature fields are known to offer NG departures by enhancing the amplitude of the three-point correlation function of ζ. For example, a light  or massive  isocurvature field ψ may affect the dynamics of curvature perturbations during horizon crossing and/or at superhorizon scales due to a special type of derivative coupling of the form L int ∝ζψ. ( This interaction appears in the Lagrangian describing the dynamics of fluctuations when the inflationary path experiences a turn in the multiple-field target space [32][33][34]. This derivative coupling, or other types of couplings appearing at higher order, may communicate different kinds of nonlinearities present in the multiple-field space to the curvature perturbation ζ. In this work we extend this mechanism [involving the derivative coupling of Eq. (1)] and show that the probability distribution function (PDF) of primordial curvature perturbations may inherit a novel class of non-Gaussianity. It relies on the existence of an isocurvature field ψ that acquires non-Gaussian statistics through its own self-interactions [78], which are transferred to curvature perturbations on superhorizon scales. The selfinteractions of the isocurvature field are determined by a given potential ∆V (ψ) characterized for having a rich structure within a field range smaller than H, the Hubble expansion rate of the Universe during inflation. More precisely, the potential for the isocurvature field is assumed to have several local minima and barriers separated by a characteristic distance ∆ψ. We also assume that the scale of the barrier height ∆V 1/4 is smaller than H (see Fig. 1). The resulting picture is reminiscent of the string landscape [79] (see Ref. [80] for a recent discussion on certain implications of the landscape picture for cosmology). In this situation the isocurvature field will be able to vigorously jump and diffuse across the FIG. 1. Illustration of a situation where the isocurvature mode would be able to fluctuate with an amplitude that would traverse several minima of the landscape. We are interested in situations where both the characteristic distance between close minima ∆ψ and the characteristic height of the potential barrier ∆V 1/4 are smaller than H, which is the typical amplitude of fluctuations around horizon crossing.
barriers. As a result, its PDF will be such that, after horizon crossing, it becomes more probable to measure a given amplitude of ψ that coincides with a local minimum of the landscape potential. As a consequence, the PDF of ζ will inherit a similar non-Gaussian profile.
The type of non-Gaussianity transferred to the curvature fluctuations cannot be fully parametrized with low n-point correlation functions such as three-and fourpoint spectra. Instead, to describe the type of NG that we encounter here, one needs to take into account a larger set of n-point correlation functions, revealing information about the nontrivial structure of the ζ PDF. This means that, in order to constrain or unveil the effects discussed in this article, one needs methods different from those often employed to test NG. More to the point, to test this class of NG one needs to constrain the entire shape of the PDF for ζ. Methods to constrain the primordial PDF have already been introduced in the past [15,16] but understandably have received much less attention than those useful to constrain the three-and four-point functions. In this work, we will not focus on such methods. Instead, we will study the mechanism underlying the emergence of this type of non-Gaussianity in a specific well-motivated example, involving an axion field.
The main result that we derive here is the PDF for ζ in the particular case that ∆V (ψ) = Λ 4 [1 − cos(ψ/f )]. This PDF is given by Eq. (106) and plotted in Figs. 5-7 for various choices of parameters. In a companion Letter [81], we argue that this result can be extended to more general potentials ∆V (ψ) and show that the shape of the PDF is entirely determined by the potential ∆V (ψ) for ψ. Together with [81], the results of this article imply that, in principle, this type of NG is able to provide information about the shape of a section of the landscape at the time of horizon crossing: a landscape tomography. Whereas, by restricting oneself to the moments of this PDF (i.e., the n-point correlation functions) one completely misses this structure.
If observations ever confirm the existence of this type of non-Gaussianity, we would have a concrete way of identifying a new energy scale parametrizing the field range of the landscape potential involving degrees of freedom additional to the inflaton. Indeed, the shapes in momentum space of the n-point functions leading to the PDF are of the local type. These shapes cannot be obtained through self-interactions arising in single-field inflation since they would violate Maldacena's consistency relation [9] (or more generally, the soft limit theorems in single-field inflation [82][83][84][85][86][87][88][89][90][91]). Therefore, tomographic non-Gaussianity can only be attributed to degrees of freedom (fluctuating across the landscape) whose effect on the spectra cannot be integrated out.
The present article has been organized as follows: We begin in Sec. II by describing the general mechanism by which non-Gaussianity may arise in the distribution of curvature perturbation due to isocurvature selfinteractions. In Sec. III we introduce the specific perturbation theory (and describe the in-in formalism) that will be used to compute n-point correlation functions. In Sec. IV we discuss the linear theory that arises in the absence of the potential ∆V (ψ), but in the presence of the derivative interaction, coupling together the curvature and isocurvature perturbations, responsible for the transfer of the NG statistics. Then, in Sec. V we proceed to compute all n-point correlation functions due to the self-interactions of the isocurvature field in a nonperturbative manner. In Sec. VI we use the general expression obtained in the previous section to derive the probability distribution function incorporating the non-Gaussian structure reflecting the axionlike potential. We discuss our results and provide our concluding remarks in Sec. VII. In order to alleviate the exposition, we have left some (important) material for the appendices: In Appendix A, we offer a concrete example of a specific UV complete model from where the action for perturbations used in Section III descents. Last but not least, in Appendices B and C we provide details of some steps omitted in the computations of Secs. V and VI, respectively.

II. DESCRIPTION OF THE MECHANISM
It is well known that, in two-field models of inflation with turning trajectories, to quadratic order in the fluctuations, the evolution of the comoving curvature perturbation ζ coupled to a single isocurvature field ψ is described by the following Lagrangian [32][33][34]: where ≡ −Ḣ/H 2 is the usual first slow-roll parameter. Here, µ corresponds to the so-called entropy mass of ψ. In the long wavelength limit, ψ satisfies the following equation of motion (obtained after integrating the equation of motion for ζ once) from where it is possible to read that µ determines the mass of ψ on superhorizon scales. Notice that ζ and ψ interact through the coupling α appearing in the special combination determining the kinetic term of ζ. In general the coupling α depends on time, and its appearance may be understood as the consequence of bends of the inflationary trajectory in the multifield target space (or more precisely, nongeodesic motions in target space) [33,34,44,45,74,75]. If the entropy mass vanishes (µ = 0), the field ψ becomes "ultralight", and the system gains a symmetry given by where C is an arbitrary constant. The consequences of this symmetry were investigated in Ref. [43]. We summarize the findings of [43] as follows: First, the symmetry of the Lagrangian (2) under the transformation (5) ensures the existence of a constant solution for ψ. This can be seen directly in Eq. (3). This solution, say ψ * , spontaneously breaks the symmetry, and dominates on superhorizon scales. Second, the symmetry of the Lagrangian under the transformation (6) implies that the constant solution ψ * will source the evolution of ζ on superhorizon scales. Concretely, if for simplicity we assume that α is nearly constant, on superhorizon scales one finds: If we conveniently identify ψ * as the value of the field at horizon crossing, then ∆N corresponds to the number of e-folds after that event. A given n-point correlation function is then given by: In the particular case of n = 2, we obtain a relation between the power spectrum of ζ and the power spectrum of ψ: Moreover, it is possible to show that ζ has a negligible influence on the evolution of ψ, which behaves as a massless field before and after horizon crossing [recall our comments about µ after Eq. (2)]. This implies that the power spectrum for ψ is given by where H * is the Hubble parameter evaluated at horizon crossing. Now the key issue to stress about Eqs. (8) and (9) is that the statistics of the field ζ is completely determined by the statistics of ψ. In this case, given that we are only considering a quadratic Lagrangian without higher order self-interactions for ψ, the statistics of ψ is found to be Gaussian. Then, the statistics inherited by ζ is also found to be Gaussian, with non-Gaussian deviations suppressed by slow-roll parameters as usual [9]. Now, in more general situations we expect a selfinteraction affecting the dynamics of the isocurvature field ψ. Thus, instead of the Lagrangian (2), we may consider the following extension: Notice that the potential ∆V (ψ) is replacing the initial mass term of Eq. (2). Without any concrete knowledge of ∆V (ψ) we would expect that it could be expanded in a power series of the form: However, such an expansion assumes that the coefficients µ 2 , g, etc... are such that, for amplitudes of ψ characteristic of horizon crossing, the higher order terms of the expansion remain suppressed. In this work we want to explore those situations where the fluctuations ψ are such that we cannot disregard the structure of ∆V (ψ) by assuming the hierarchical expansion of (12). To make this statement more concrete, we will consider the following specific axionlike potential: where f is the axion decay constant. In Appendix A we shall describe a concrete example wherein this Lagrangian emerges.
To continue, notice that the potential (13) breaks the shift symmetry (for µ = 0) of the Lagrangian (2) down to a discrete symmetry: This time, ψ may acquire constant solutions that minimize the sinusoidal potential. On superhorizon scales, any of these solutions will dominate the behavior of ψ. Just as before, ζ will be sourced by ψ, but this time the enhancement will happen for those values of ψ that minimize the potential. This result suggests that the statistics transferred from ψ to ζ will continue to be operative in this new context, but in a manner that it will be enhanced at those values in which ψ coincides with a minimum of the potential, and suppressed for those values in which ψ coincides with a maximum. Therefore, the structure of the potential will be necessarily inherited by the PDF of the curvature perturbation ζ, which becomes non-Gaussian. Presumably, the Lagrangian (11) is the result of perturbing a more fundamental multifield theory, with a scalar field potential of the form V = V 0 + ∆V [as already mentioned, in Appendix A we show that this is one way to derive (11)]. We will be interested in the regime Λ 4 /3H 2 M 2 Pl 1, so the potential ∆V has little to say about the background dynamics of the full system, and inflation is driven by a piece V 0 . Then, the background equations of motion require V 0 ∼ 3H 2 M 2 Pl . In addition, if Λ 4 /3M 2 Pl H 2 1 then the symmetry breaking is mild and, for all practical purposes, before and during horizon crossing the field ψ will behave as an ultralight field. This means that at horizon crossing ψ will freeze, and Eq. (7) will describe how ψ transfers its statistics to ζ. As times passes, the nonlinearities due to ∆V (ψ) will start to become accentuated, and one expects a nonlinear contribution to Eq. (7) coming from the nonlinear evolution of ψ that does not freeze. To leading order in α we expect that any level of nonlinearity will be communicated to ζ through a non-Gaussian contribution to the n-point correlation functions of the form The reason behind this guess is the following: First, in the absence of interactions between the two fields (α = 0), the fluctuation ψ will acquire a non-Gaussian distribution due to its potential ∆V (ψ) The non-Gaussian contribution to the n-point correlation functions were computed in Ref. [78] using the in-in formalism, and are found to be given by where σ 2 0 is the variance of ψ appearing from loop resummations. In the previous expression the subscript c indicates that we are only taking into account the diagrammatic contributions due to the potential ∆V (ψ) that are fully connected [which is why there is a single overall Dirac-delta function on the right-hand side of (17) enforcing momentum conservation]. This set of n-point correlation functions are generated during horizon crossing.
Second, let us turn back on the coupling α = 0. Then, because we are assuming that ∆V (ψ)/3H 3 M 2 Pl 1 the field ψ is essentially massless and the linear relation (7) will remain valid on superhorizon scales, independent of the nonlinear dynamics. This implies that any non-Gaussianity gained by ψ during horizon crossing will be transferred to ζ via Eq. (16) after horizon crossing. As we shall see, a detailed computation leads to where the factor 1/2 comes from the interaction structure implied by certain nested integrals appearing in the computation of the n-point correlation functions using the in-in formalism.
The importance of Eqs. (17) and (18) is that they allow us to infer a probability distribution function for ζ. This PDF is characterized by a class of non-Gaussianity that cannot be fully captured with three-or four-point correlation functions, as opposed to the case springing from the ansatz (12). This PDF is found to be given by (see Sec. VI E for the derivation) where σ ζ is the variance of ζ parametrizing the Gaussian part of the distribution. In the previous expression K(x), f ζ (x) and A are given functions and quantities determined by parameters related to ∆V (ψ) that will be deduced in the next sections. The main characteristic of the distribution function ρ(ζ) is that, in spite of the x integral, it inherits the structure of the potential ∆V (ψ). That is, the probability of measuring ζ increases (decreases) if the field ψ sourcing its amplitude is at a local minimum (maximum) of ∆V (ψ). The mechanism described here is certainly not exclusive to the potential ∆V (ψ) given in Eq. (13). It should be safe to suspect that any potential ∆V (ψ) with a rich structure (i.e., characterized by field distances ∆ψ smaller than H) will imply the existence of some level of non-Gaussianity for ζ revealing the structure of ∆V (ψ). In fact, we will prove this statement in the companion Letter [81]. Thus, we see that the type of non-Gaussian departures discussed here in principle gives us nontrivial information about the landscape, offering us tomographic information about the shape of the multifield potential.
Before finishing this section, let us mention that the field ψ is not expected to be a true axion as realized in QCD or string theory [92,93] for the range of parameters that we are interested in. The reason is that large fluctuations of ψ traversing many minima of the potential would destabilize the radial field fixing the value of the axion decay constant f [94]. For this reason, we take the potential of Eq. (13) to be representative of systems with potentials with a rich structure, as expected in the string landscape. See Ref. [95] for a previous work that has studied the system (11) with an axionlike potential (for the decoupled case α = 0) analyzing issues related to the landscape. For a recent review on the role of axions in cosmology and inflation see Ref. [96] and references therein.
In the rest of this article, we set ourselves to derive Eqs. (18) and (19) using the in-in formalism, giving details about every step of the computations.

III. SETTING THE STAGE
One of the main technical goals of this article is to compute n-point correlation functions of ζ resulting from the potential ∆V (ψ) given in Eq. (13). To proceed, we introduce canonical fields u and v defined as: To simplify any computation involving u and v we will assume a purely de Sitter background, with a(t) = e Ht . This means that our computations will have corrections of order not accounted for. It is also convenient to introduce conformal time τ (via dτ = dt/a = dte −Ht ) and write a(τ ) = −1/Hτ, with −∞ < τ < 0. Then, putting together Eqs. (20)- (22) back into the action (2), we obtain where we have introduced the dimensionless coupling which is taken to be nonvanishing in the de Sitter limit. From Eq. (23), we infer that the canonical momenta associated with u and v are, respectively, given by These momenta satisfy the equal time commutation relations, given by with every other commutator vanishing. From (25) and (26) we see that the Hamiltonian of the system is given by where we have introduced the notation x = d 3 x.

A. Splitting the theory
We may now split the Hamiltonian into three contributions as H = H 0 + H λ + H Λ , where H 0 corresponds to the free Hamiltonian of the system (obtained in the limit λ = Λ = 0). Notice that H 0 describes a system with two decoupled massless scalar perturbations On the other hand, H λ contains the interaction term proportional to λ, and H Λ contains the self-interactions for v We may now quantize the system by adopting the interacting picture framework. That is, the quantum fields u and v are expressed as where u I (x, τ ) and v I (x, τ ) are the interaction picture fields, which evolve as quantum fields of the free theory. Explicitly, they are given by with k = (2π) −3 d 3 k, and wherê Here, the pairs a ± (k) and a † ± (k) correspond to creation and annihilation operators satisfying the commutation relations: with b, c ∈ {+, −}. The mode functions u k (τ ) and v k (τ ) are both given by which corresponds to the standard expression for a massless mode on a de Sitter spacetime with Bunch-Davies initial conditions. On the other hand, U (τ ) is the propagator in the interaction picture, which is given by where T stands for the time ordering symbol. In a given product of operators, T instructs us to put operators evaluated at later times on the left and operators evaluated at earlier times on the right. In addition, we have where is a small positive number introduced to select the correct interaction picture vacuum.
Finally, H I in Eq. (41) is given by where Notice that in the previous expressions the canonical momenta Π I u and Π I v in the interaction picture are simply given by In order to deal with H Λ I , we will consider the Taylor expansion of the cosine function as: This expansion gives us an infinite number of vertices to deal with. As we shall see in Sec. V, it will be possible to resum this expansion back into an exponential contribution, leading to nonperturbative results in terms of the ratio H/f . The computation of n-point correlation functions of the following sections may be organized diagrammatically. These computations involve contractions of the Hamiltonians H λ I and H Λ I and the fields u I and v I . In the present context, a contraction is the result of a commutation between creation and annihilation operators introduced in Eqs. (37) and (38) that results from their normal ordering. A commutation involving the pair a † − and a − is represented by a dashed line joining vertices (or external legs) labeled by the u fields that contain these operators. Similarly, a commutation involving the pair a † + and a + is represented by a solid line labeled by the v fields that contain these operators. If the field participating in the commutation comes from H λ I , then the line joins a vertex with an empty circle. Otherwise, if the field participating in the commutation comes from H Λ I , then the line joins a vertex with a solid circle. Figure 2 shows the various classes of diagrams appearing in the computation of n-point correlation functions.

B. Perturbativity conditions
In the next sections, we compute the n-point correlation functions perturbatively. Given that we will resum the expansion in H/f appearing in (46), we will not impose any condition on the size of f . On the other hand, it is worth counting with some criteria on how large λ and Λ can be in order to have a well-behaved perturbation theory. A naive estimation (that does not take into account renormalization) may be obtained by rewriting the Lagrangian (11) in a dimensionless form, weighting spacetime variables and fields by their characteristic values. In de Sitter, a characteristic length scale is given by H −1 . Moreover, the amplitude of massless scalar fields around horizon crossing is of order H. Thus, redefining spacetime and field variables as and after writingL = L/H 4 , the Lagrangian (11) becomes where λ = √ 2 α/H is the dimensionless coupling already defined in Eq. (24). Here, derivatives are with respect to the dimensionless variables, and we have further defined the ratio g ≡ H/f . By asking that the dimensionless couplings remain small, we obtain the following perturbativity conditions for the potential (13). Note that the first condition is stronger than Λ 4 /3M 2 Pl H 2 1 already required for the background not to be affected by ψ. Also, recall that we are not restricting the value of g = H/f , as the results of the next sections are nonperturbative with respect to this parameter.
As we shall see in Sec. V, the loop corrections due to the resummation of the sinusoidal potential (46) will renormalize the bare coupling Λ. The consequence of this is that the correct perturbative parameter will turn out to be where σ 2 S is a short wavelength contribution to the variance of the field ψ (we will compute this quantity in Sec. VI A). In summary, our results will be perturbative in the couplings Λ 4 ren /H 4 and λ, but nonperturbative in the parameter H/f because (46) will be resummed eventually.

IV. LINEAR THEORY
Before computing the n-point correlation functions using the full nonlinear theory, let us have a look into the linear theory obtained in the limit Λ 4 → 0. In the absence of nonlinear interactions, the statistics will be Gaussian, and the only meaningful quantity to compute is the power spectrum for ζ. This system was investigated in Ref. [43], and here we show the main steps allowing one to deduce the value of the power spectrum. This result will be important later on.
To start with, notice that if Λ = 0 the system consists of two canonical massless fields u and v coupled through the interaction Hamiltonian (43). The power spectrum for ζ may be obtained by computing the two-point correlation function u(x, τ )u(y, τ ) . We will perform this computation up to order λ 2 . Given a fluctuation ϕ, we define its power spectrum P ϕ (k) as To proceed with the computation of P ζ , we write u(x, τ ) = U † (τ )u I (x, τ )U (τ ), where U (τ ) is given by Eq. (41), with H I as in (43). Up to second order in λ this quantity is given by The two integrals may be solved in the superhorizon limit k|τ | 1. The first integral is found to be given by whereas the second integral reads These expressions, together with (54), allow one to compute the two-point correlation function u(x, τ )u(y, τ ) .
In momentum space, one obtains where A 1 and A 2 are numbers given by Their numerical values are A 1 −2.11 and A 2 1.46. Note that in putting together (54)- (56) to compute the two-point correlation function, the divergent logarithms ln(−kτ ) cancel out. The computation of the two-point correlation function of Eq. (57) may be thought of as the result of the diagrammatic expansion of Fig. 3. The zeroth order contribution corresponds to the first diagram, whereas the contribution of order λ 2 corresponds to the The two diagrams contributing to the computation of the two-point function. The first diagram gives the standard power spectrum for ζ, whereas the second diagram gives the correction due to the λ derivative interaction.
second diagram, where the two external legs are mediated by a v propagator. Now, horizon crossing happens when k|τ | 1. Thus, the number of e-folds after horizon crossing is given by It may be seen that after several e-folds the contribution to the power spectrum (57) quadratic in ∆N dominates, and we obtain For this to happen, we need λ 2 ∆N 2 1. The power spectrum for v may be found through a similar computation, which gives valid up to order O(λ 4 ). Combining these two results, we then derive the following relation between the two power spectra where P ψ (k) is found to be given by This result is consistent with the behavior shown in Eq. (7) based on symmetry arguments. It shows that the power spectrum for ζ is proportional to the power spectrum of ψ, with a factor that grows with the number of e-folds. Let us briefly comment on the validity of the result shown in Eq. (61). Our perturbative method consisted of separating the theory between a zeroth order Hamiltonian H 0 [given in Eq. (30)] and an interaction Hamiltonian (given in Eq. (43)), proportional to λ. On the one hand, we have argued that, in our final result for the power spectrum (57), we are allowed to retain as the dominant piece the term proportional to λ 2 . On the other hand, notice that our perturbative method is valid as long as λ 2 1. These two statements are not in contradiction: the computation admits a cumulative effect that grows with the number of e-folds as λ 2 ∆N 2 , which may be larger than 1 (after ∆N 60). This effect was discussed in detail in Ref. [43], and it will play an important role in Sec. VI.
Also, the example of the derivative coupling we used here has a special property that, at superhorizon scales, the linear equation of motion for the isocurvaton field ψ has no source term from the curvature mode ζ. Therefore, ψ does not grow once it exits the horizon, while ζ does. This means that, if we were to solve the coupled linear equation iteratively to all orders, the enhancement factor from ∆N would stay at the order ∆N 2 . Therefore, the requirement of the perturbation theory is only that λ 1, and λ∆N can be greater than 1. Last but not least, even if the two conditions λ 2 ∆N 2 1 and λ 2 1 may seem to be fine-tuned, the condition λ 2 1 has only been adopted in order to be able to perform analytic computations. The requirement λ 2 ∆N 2 1 is valid independent of the perturbativity condition λ 2 1, and can already be inferred from the symmetry arguments around Eq. (9). In principle, Eq. (9) [or Eq. (63)] should be valid independently of the value of λ.

V. COMPUTATION OF CORRELATION FUNCTIONS
In this section, we describe how to compute the npoint correlation functions of ζ at the end of inflation (details of this computation are shown in Appendix B). The quantity of interest corresponds tõ which, in terms of the canonically normalized fields introduced in Sec. III, is given bỹ Our main goal is to obtain an expression for this function up to order Λ 4 . Given that the interaction Hamiltonian (44) determined by the potential ∆V (ψ) does not depend on u I , the n-point correlation functions will acquire a dependence on Λ 4 only through the mixing Hamiltonian (43) involving the coupling λ. This means that the fully connected contribution to (66) will necessarily involve at least one factor λ per field u(x n , τ ). In other words, the lowest order contribution to (66) represented by fully connected diagrams, will be of order Λ 4 λ n . The diagrammatic representation of this computation in momentum space is shown in Fig. 4. The ⊗ vertex denotes the exact vertex, up to order Λ 4 , connecting the n external legs participating in (66). Because the expansion (46) contains an infinite number of vertices (with an even number of legs), the exact vertex consists of the sum of an infinite number of diagrams involving loops that start and finish on the same vertices. Due to overall momentum conservation, these loops do not carry external momenta.
To proceed with the computation ofG (n) ζ , we start by recalling that a given u(k j , τ ) entering the n-point function û(k 1 , τ ) · · ·û(k n , τ ) of Eq. (66) has the form u(k j , τ ) = U † (τ )u I (k j , τ )U (τ ). Let us for a moment disregard the prescription determining the integration limits ∞ ± . Then, by expanding the propagator U (τ ) in this expression, we obtain We only need to keep terms up to order Λ 4 . Then, because H Λ I commutes with u I , but not with v I , the previous equation may be further reduced to where the ellipses of the last line denote terms that will not contribute to the piece that we want to compute. For instance, by inserting another H λ I Hamiltonian (through a commutator) between τ and τ j we would be computing a correction to the ζ propagator and not to the fully connected part of order Λ 4 λ n . Now, Eq. (68) tells us that the structure of u(k j , τ ) in terms of creation and annihilation operators a ± is of the following form: The computation of n-point correlation functions requires us to plug this form of u back into (66) and perform every possible contraction between creation and annihilation operators of the various terms appearing in (69).
The final result that we are pursuing is an expression containing only terms of order Λ 4 λ n , and thus many of the contractions correspond to loops involving pairs of a + operators. The diagrammatic expansion of this computation is shown in Fig. 4. Since we are computing the fully connected contribution to (69), in every contraction involving a + operators, at least one of them must come from a term of order Λ 4 in (69).
The details of this computation are shown in Appendix B. The final result is found to be given bỹ for even n, as it vanishes for odd n because the potential is even under ψ → −ψ. Here σ 0 is the variance of ψ, defined through the relation The factor exp(−σ 2 0 /2f 2 ) of Eq. (70) appears as the consequence of the resummation shown in Fig. 4. This result may be compared to that shown in Eq. (17) and obtained in Ref. [78]. This comparison proves the result quoted in Eq. (18) of the Introduction.
In Eq. (70), ∆N corresponds to the number of e-folds elapsed since a reference time τ 0 around which the set of modes k l crossed the horizon (k l τ 0 ∼ 1) This definition coincides with that of Eq. (60) in the case of the two-point function. Notice that the shape in momentum space of the non-Gaussianity parametrized by these n-point correlation functions is of the local type.
In obtaining this result we have assumed that the condition ∆N 1 holds. If instead one has a relatively small ∆N , then other terms that were neglected might need to be included back. Given that our perturbativity conditions demand λ 1, this is in agreement with the power spectrum shown in (63), valid for λ 2 ∆N 2 1.

A. About horizon exit
Here let us make a remark on the value of ∆N . Since the ψ field in our model example is exactly massless, ∆N should be evaluated from when modes exit the horizon until the end of inflation. Because longer modes exit the horizon earlier, ∆N is not exactly a constant and has a logarithmic dependence on the magnitude of relevant momenta. Here we ignore this weak momentum dependence and approximate ∆N ∼ 60 as a constant. On the other hand, for our purpose we do not have to regard the mass of the scalar field ψ to be exactly zero (but still require µ H, so our model conditions are satisfied). Such a field would have decayed before it stays at the superhorizon for the entire 60 e-folds, as it happens in the lower mass range case of the quasi-single-field inflation models [44,45]. In this case, all modes of ψ will stay for the same amount of e-folds, ∆N , after the horizon exit and before the decay, and 1 < ∆N < 60 is now exactly a constant.
All the connected diagrams contributing toG (n) at order Λ 4 . The ⊗ vertex represents the resummation of all the loop contributions coming from the expansion of the cosine function shown in Eq. (46). In other words, for a given number of legs, the ⊗ vertex contains all the relevant effects due to the cosine. Because of the combinatorial factors of each diagram, and given that there are no external momenta running through the loops, the resummation reduces to a constant factor given by exp(−σ 2 /2f 2 ).

B. Regularization: IR and UV cutoffs
Before we finish this section, let us briefly come back to Eq. (71). If we replace the mode solutions u k (τ ) of Eq. (40) back into (71), and define the dimensionless integration variable q = k|τ |, we obtain Observe that σ 0 is independent of time τ . However, it contains divergences coming from the integration limits q → 0 and q → +∞. We may therefore introduce infrared and ultraviolet cutoff scales q IR and q UV , respectively, and obtain Notice that the variable q = k|τ | = k/aH is the physical momenta per unit of the Hubble scale. The UV cutoff refers to a scale that is deep inside the horizon and is the scale of new physics and the limit of low energy effective theory. This cutoff contributes to the renormalization of the coupling Λ 4 as in what happens in flat spacetime. The logarithmic IR divergence is due to the random walk of the massless field in the dS space. Actual observations do not have access to all the scales, and so the IR cutoff should be set by the size of the observable Universe. We will come back to this issue in Sec. VI A, where we consider the need of defining the variance of modes available to cosmological observers.

VI. TOMOGRAPHIC NON-GAUSSIANITY
The expression for the n-point correlation functions given by (70) may be Fourier transformed back into coordinate space as This expression may be used to deduce the probability distribution function of measuring an amplitude ζ(x) at a given position x. To this end, we need to compute the moments ζ n that are given by evaluating all the coordinates in (75) at a common value x where the subscript c denotes that ζ n c comes from fully connected diagrams, hence it is proportional to a Diracdelta that conserves momentum. Due to this momentum conservation, ζ n c defined in Eq. (76) is independent of x. In the following subsections, we first obtain a concrete expression for the n-point functions of Eq. (76) valid for long wavelength modes, and then we proceed to derive the PDF from where these n-point functions are computed.
A. n-point functions for long wavelength modes The quantityG (n) ζ (τ, k 1 , · · · , k n ) of Eq. (70) only shows the leading IR contribution to the full n-point correlation function. For the same reason, in (75) we cannot integrate along the entire momentum space, and we are forced to introduce a cutoff momentum k L . This is not a technical limitation, but all the contrary. We are interested in making predictions of inflation valid for superhorizon perturbations (that will later on reenter the horizon after inflation), and so we want to compute correlation functions of long wavelength ζ modes. This is normally done by introducing window functions selecting the relevant scales for the computation of correlation functions in coordinates space. For simplicity, here we consider a window function with a hard cutoff k L . With this purpose in mind, we introduce the cutoff in terms of physical momentum q phys ≡ k|τ | (per unit of H) instead of comoving momentum k. That is, we choose a hard cutoff momentum q L and split the curvature perturbation as where ζ L only includes modes of wavelengths larger than some fixed value 2π/q L . Horizon crossing happens at q phys = 1, and so we must impose In other words ζ L contains superhorizon contributions (at the end of inflation) between the physical cutoff scales q L and q IR . Explicitly, ζ L is given by where k L = q L /|τ |. Thus, we will compute a more restricted version of (76) given by where G (n) ζ,L (τ, x 1 , · · · , x n ) reads as in (75), but now with the momenta integrated up to k L . Explicitly, we have The division of scales (77) forces us to split σ 2 0 = ψ 2 , introduced in Eq. (73), into short and long wavelength contributions as σ 2 0 = σ 2 S + σ 2 L , in such a way that σ 2 S and σ 2 L receive contributions larger and smaller than k L respectively. From Eq. (73) we see that σ 2 S and σ 2 L are given by where we have introduced the ratio which is a measure of the range of scales spanned by the long mode contributions ζ L . The logarithmic dependence of (83) suggests that σ 2 L H 2 . For instance, if we take q L = 1, then ξ corresponds to the ratio between the largest wavelength and the Hubble radius at the end of inflation. Then ln ξ 60, and one obtains σ 2 L H 2 . Notice, however, that, in general, ξ parametrizes the window function selecting the scales, hence its value should be determined by the range of momenta available to cosmological observations. In the particular case of the CMB, this ratio is approximately given by ln ξ ∼ 8.
Next, to obtain an expression for ζ n L c , we evaluate the arguments of (81) at a single coordinate value x. Because of momentum conservation, the argument of the exponential vanishes, and we are left with the following expression where Here, Λ 4 ren = e −σ 2 S /2f 2 Λ 4 is the renormalized coupling introduced in Eq. (52) resulting from the loop resummation introduced in Sec. V. Because this resummation is always induced by Λ, the combination Λ 4 ren will be present at all orders in perturbation theory (disregarding higher order loop corrections carrying external momenta that start appearing at order Λ 8 ). In Ref. [81] we discuss the renormalization of ∆V (ψ) more generally, paying special attention to the running of the parameters defining ∆V (ψ) in order to make observables independent of the cutoff scales. To continue, the coefficient g n in (85) is defined as for even n and zero otherwise becauseG vanishes if n is odd. Here I n corresponds to the following integral: Equation (85) is written in terms of the variance σ 2 L associated with the probability distribution function of ψ. It will be more useful to write ζ n L c in terms of the variance σ 2 ζ instead of σ 2 L . Recall that in Sec. IV we derived the power spectrum of ζ in terms of the power spectrum of ψ, given in Eq. (63). When λ 2 ∆N 2 1, this result implies Then, by defining f ζ as it is direct to find This is the general form of the n-point correlation function that we need in order to reconstruct the tomographic PDF for ζ. It is important to emphasize here that we can also consider the regime λ 2 ∆N 2 < 1, where one finds σ 2 ζ = σ 2 L /2 . This case would give us a slightly different expression for (91) but would not change the form of the reconstructed PDF (except for the way in which some parameters appear). For simplicity, we stick to the regime λ 2 ∆N 2 1.

B. Dependence of the n-point functions on the cutoff scales
The integral I n determining the form of the factor g n [through Eq. (87)] is a function of the order n and the ratio ξ introduced in Eq. (84). Indeed, in Appendix C, we show that I n can be written in terms of a single integration variable as where G(ξ, x) and F (ξ, x) are given by Here Ci(x) is the cosine integral function. In Appendix C, we also show that in the formal limit ξ → ∞, the integral asymptotes to a simple function I 0 n (ξ) given by Then, given that ln ξ = 4π 2 σ 2 L /H 2 [which, can be seen from Eq. (83)], one obtains g n = n, implying a very simple general expression for the n-point functions ζ n L c . However, given that k L is at most the horizon exit scale, and that k IR is the largest scale available to present observers, we have the bound ln ξ ≤ 60, which implies that ln ξ is too small to allow us to take I n as I 0 n . The reason for this is that, with ln ξ 60, the correction ∆I n = I n − I 0 n is already one-tenth of I 0 n for n ∼ 35. This in turn implies that the PDF derived with I 0 n starts to deviate significantly from the one derived with I n when f is smaller than σ L ∼ H, which is precisely the interesting region of parameters that we wish to explore. A proof of this statement is given in Appendix C.
In what follows, we devote ourselves to reconstruct the PDF out of the n-point function ζ n L c given in Eq. (85). We will first do this in Sec. VI D for the case in which I n is taken to be as I 0 n (ξ) shown in Eq. (95). Then, in Sec. VI E, we will show how to obtain the PDF for the full expression for I n (ξ) shown in Eq. (92). Before deriving these two PDFs we describe the general idea behind its reconstruction.

C. PDF reconstruction: general idea
Recall that we have focused our interest on the computation of the fully connected contributions to the n-point correlation functions. Had we focused instead on the full n-point correlation functions, including disconnected diagrams, we would have arrived at the more general expression Here, the factors σ 2m ζ come from propagators connecting pairs of external lines. The combinatorial factor n!/m!(n − 2m)!2 m consists of the total number of ways to connect the n external legs in such a way that 2m of them are connected by propagators, and the rest are connected to the Λ 4 vertex.
The probability distribution function ρ(ζ) that we are searching for must be such that To find ρ(ζ), it is useful to notice that the term m = n/2 in Eq. (97) is given by which corresponds to the n-point correlation function of a Gaussian distribution. This means that ρ(ζ) is given by a leading Gaussian distribution with a non-Gaussian correction proportional to Λ 4 . Thus, to find ρ(ζ) we may try the following ansatz where is the Gaussian part giving rise to the subset of n-point functions given in Eq. (99). The piece ∆ρ(ζ) corresponds to the correction resulting from the nonlinear interactions proportional to Λ 4 (or equivalently, A 2 ).
In what follows, we determine the form of ∆ρ(ζ) due to ζ n L c shown in Eq. (85). The procedure crucially depends on knowing how ζ n L c depends on n, which requires us to deal with I n (ξ) of Eq. (92). To proceed, we find it useful to first show how to deduce ∆ρ(ζ) in the case where I n (ξ) is given by its asymptotic form I 0 n (ξ). This will then allow us to deal easily with the more general situation in which I n (ξ) is given by the full expression given in (92).

D. Asymptotic reconstruction
If we take Eq. (85) with I n (ξ) replaced by its asymptotic form I 0 n (ξ) given in Eq. (95), then g n = n when n is even, and in this case we simply have (and zero if n is odd). This equation may be used to derive the PDF ρ(ζ) determining the probability of measuring a certain value of the curvature perturbation at an arbitrary position. To find ∆ρ(ζ) we may try the following ansatz where we have used the fact that ρ(ζ) must be even under the change ζ → −ζ, for only the even moments are nonvanishing. It is direct to find that B 0 and C 0 are the only nonvanishing coefficients and, therefore, that the full PDF is given by [97] (104) This PDF satisfies Eq. (97), and given that it corresponds to a small, absolutely continuous deformation of a Gaussian distribution, it is unique [that is, it is the only possible reconstruction from the moments of Eq. (102)].
The probability distribution function (104) is valid in the formal limit ξ → ∞. If σ L H we could trust this result for the case f σ L , in which case the PDF shows nontrivial structures in the form of superimposed oscillations. However, given that σ L ∼ H (because ln ξ ≤ 60), we cannot trust the regime f σ L (see Appendix C), and we are forced to consider the more general case in which I n is given by its full form shown in (92). In spite of this limitation, Eq. (104) constitutes one of our main results. It gives a simple non-Gaussian probability distribution function for ζ in terms of various parameters related to the landscape shape. It may be verified that the PDF is already normalized as ρ(ζ)dζ = 1. This probability distribution function is plotted in Fig. 5 for specific values of f ζ /σ ζ and A 2 . Notice that the second term inside the squared parenthesis in Eq. (104) accounts for the increase in probability of finding values of ζ that are sourced by those values of ψ which minimize the cosine potential of Eq. (13). On the other hand, the third term, linear in ζ may be interpreted as a contribution accounting for the diffusion of ζ (one could in fact absorb the third term into the second term by slightly shifting f ζ ).

E. Full reconstruction
We now consider the task of deriving the full PDF, valid for any value of ln ξ > 0. To proceed, it is helpful to realize that the most important aspect about the reconstruction performed in the previous section was the dependence of ζ n L c on n as shown in Eq. (102). In the general case, if we consider the x integral of Eq. (92) explicitly, we see that the dependence of ζ n L c on n is exactly the same, except that this time it happens for each value of x. Then, a simple comparison with (85) shows that this time the reconstruction amounts to identifying an x dependent decay constant that satisfies f ζ (0) = f ζ . Hence, keeping track of all the numerical factors, we find where the kernel K ξ (x) is given by (107) The result shown in Eq. (106) is our main result. It gives us the PDF for any value of the ratio ξ = k L /k IR . In particular, we can trust this result well inside the regime f ζ < σ ζ for values ln ξ 8, which corresponds to the range of scales available to CMB observations [as opposed to the case of the PDF of Eq. (104)]. An outstanding property of (106) is that it preserves the oscillatory structure of the potential in a strikingly similar manner as (104). The main difference, is that now there is a filtering function that accounts for the effects that arise when one considers only the bounded region of k space which we are able to probe. The consequences of this filtering can be appreciated by looking at Fig. 6 (plotted for ln ξ = 8). There we see by comparison to Fig. 5 that the amplitude of the oscillations in the full PDF is suppressed, since the value chosen for A 2 in this last plot is 100 larger than the previous one. Moreover, as illustrated in Fig. 7 (also plotted for ln ξ = 8), decreasing the value of f ζ from this point does not enhance the amplitude as we might have thought while looking at Eq. (104) but the opposite: the amplitude actually gets smaller, mostly because of the exponential factor in the kernel K ξ (x). Note that in the formal limit ln ξ → ∞ the probability distribution function (106) becomes independent of ξ, and we recover (104). This corresponds to the ideal situation whereby the entire range of momenta is available to observers. In that case, one can directly infer the parameters Λ 4 ren and f of the landscape potential by performing statistics with observations of primordial curvature perturbations. Otherwise, as long as observers can only have access to a limited amount of modes, as parametrized by ln ξ, the filtering function appearing in the kernel K ξ (x) will wash out the structure of the potential. This is simply because ln ξ restricts the number of modes in momentum space that can add up to increase the effect of nonlinearities due to ∆V on the PDF in coordinate space: momentum conservation through the vertex implies that while certain modes are probing large (small) momenta limited by q L (q IR ), other momenta will probe more restricted regions in momentum space. As a result, to extract information about ∆V one has to take 7. An example of the PDF of Eq. (106) for the choice of parameters f ζ /σ ζ = 2 × 10 −2 , A 2 = 2.5 × 10 −2 and ln ξ = 8 (solid curve). For comparison, we have plotted a Gaussian distribution of variance σ ζ (dashed curve). Here A 2 has the same value as that of Fig. 6, but f ζ /σ ζ is smaller.
into account the role of ln ξ. We will come back to this issue in Ref. [81].

VII. DISCUSSION AND CONCLUSIONS
We have examined a regime of multifield inflation where the shape of the landscape potential in the isocurvature direction can be probed using non-Gaussianity of primordial density perturbations. At the level of npoint correlation functions (or polyspectra), these non-Gaussianities take the local form, as in all multifield models  and quasi-single-field models  with sufficiently light isocurvatons. However, a novel point of this paper is that the information about the shape of the potential is not manifest in the individual n-point correlation functions, but rather in the resummed probability distribution function given in Eq. (106). This shows that local non-Gaussianity may have a rich structure inherited by the self-interactions of isocurvature fields, together with a derivative coupling common to multifield models.
Although the mechanism of statistical transfer examined in this article is based on the derivative coupling L int ∝ζψ, our results are likely to be more general. We therefore expect that other classes of interactions between the curvature perturbation and other scalar fields scanning the landscape lead to similar conclusions.
Also, as mentioned in the Introduction, the particular form of the potential ∆V (ψ) should not be so crucial. While it is true that the cosine function used in this work comes with the right properties making the loop resummations possible, more general potentials are in fact not intractable. In a companion Letter [81], we will extend Eq. (106) to the case of more general potentials ∆V (ψ) and analyze the possibility of reconstructing such potentials out of observable quantities.
In what follows, we discuss various relevant aspects related to our main result.

A. Aspects of nonperturbativity
The probability distribution function of Eq. (106) was the result of two resummations. The first resummation corresponded to the loop expansion coming from the infinite number of vertices following from the Taylor expansion of Eq. (46). This makes our PDF nonperturbative in terms of the parameter H/f (or, equivalently, σ L /f and thus σ ζ /f ζ ). The second resummation corresponded to the derivation of the PDF itself, from the infinite set of n-point functions given by (91). This second resummation gave us back the appearance of the cosine function, weighted by the kernel K ξ (x) of Eq. (107) and with a periodicity determined by f ζ (x) as in Eq. (105). These two functions are such that the non-Gaussian correction to the PDF preserves the structure of the landscape potential ∆V .
These two steps (the resummations) work in tandem and inform us about some crucial aspects of our result. While each n-point function depends nonperturbatively on σ ζ /f ζ , by themselves they cannot give an account of the oscillatory structure of the PDF. In particular, the lowest non-Gaussian correlation function (the four-point function) can hardly inform us about the nature of the class of non-Gaussianity from where it is derived. This is the most salient point of this work: in order to correctly describe the nonlinear effects of the isocurvature field on the curvature perturbation ζ, we cannot just limit ourselves to the lowest non-Gaussian n-point correlation function.
As discussed in Sec. VI, observables depend on the renormalized parameter Λ 4 ren = e −σ 2 S /2f 2 Λ 4 , which is a result of the loop resummation. Given that this resummation will reappear at any order of the expansion with respect to Λ 4 , we should take the perturbative expansion of the hight of the potential as being controlled by the renormalized parameter Λ 4 ren /H 4 instead of Λ 4 /H 4 . In other words, higher order corrections to the PDF of Eq. (106) are expected to be of order A 2 .
It is important to stress that, while our result is perturbative in Λ 4 ren /H 4 and λ, the main features of the PDF (106) should persist to other regimes. For example, we have stayed in the region λ 1 in order to be able to treat the mixing term (31) as part of the interaction Hamiltonian. This regime has the limitation that, in order for the effects coming from λ to be relevant, one needs a large number of e-folds ∆N after horizon crossing. On the other hand, we expect that the qualitative characteristics found in the PDF of Eq. (106) springing from λ will be enhanced in the case λ 1.

B. Relation to previous works
Our analysis has some similarities (but also important differences) with previous works studying the implications of isocurvature fields on the production of primordial non-Gaussianity. In quasi-single-field inflation models, the isocurvature field is assumed to be massive and to have some interactions, such as the cubic self-interaction as in Eq. (12) In these models, the mass of ψ in principle can be a free parameter and plays an important role: it controls the extent to which the fluctuations of ψ can survive at superhorizon scales and interact with ζ. This is because the amplitude of ψ decays after horizon crossing as where R is the real part of 1 − 4µ 2 /9H 2 . Although many results apply for general µ, the most interesting cases in the quasi-single-field literature  are those with µ ∼ H. Since the potential (108) is assumed to hold within a field range much larger than the amplitude of ψ, O(H), the isocurvaton field ψ is confined within this potential and does not fluctuate outside to explore the fuller structure of the landscape. In fact, the PDF of the density perturbation of the quasi-single-field inflation model may be worked out in a similar fashion and it should encode the shape of the potential, although it has much less rich structure. The main predictions for non-Gaussianities coming from (108) are some nontrivial polyspectra, such as bispectra and trispectra.
In the case studied in this article, the field ψ does not decay. Note that, classically, ψ has a mass coming from the cosine potential of Eq. (13), and this quantity may be even larger than H 2 . But given that the potential barriers are small (∆V 1/4 ∼ Λ H at tree level) the ψ fluctuations will still be able to traverse vigorously the barriers of the potential and explore the potential landscape. In other words, the ψ field is effectively massless at the leading order. In this case, the classical mass term is only part of the rich structure in the small perturbation ∆V and appears as the first term in the series expansion of ∆V . Therefore, this series expansion needs to be resummed in the final result. In the case of the cosine potential studied here, these aspects are summarized in the fact that all the vertices depend on just two parameters (Λ 4 and 1/f ), and so every vertex contributes decisively in the computation of the n-point correlation functions.

C. After inflation
So far we have established how a scalar field ψ with a non-Gaussian distribution function can transfer its statistics to the curvature perturbation ζ during inflation. The mechanism by which the statistics is transferred is cumulative: ψ transfers its statistics (both Gaussian and non-Gaussian) to ζ as long as λ = 0, and the non-Gaussianity of ζ becomes more accentuated as time passes. After a long enough period, the statistics of ζ (which in the absence of the λ coupling would be nearly Gaussian) becomes completely dominated by that of the curvature perturbation ψ.
Three main things could happen after such a period that bring this mechanism to an end: (1) As mentioned, if ψ is not exactly massless, after some e-folds at superhorizon it will naturally decay, as in quasi-single-field inflation models. (2) The potential ∆V (ψ) changes drastically. In more realistic scenarios we would expect that the potential ∆V (ψ) depends explicitly on time due to its dependence on the background. Before the end of inflation ∆V (ψ) could introduce a new scale that makes ψ a very massive field within the relevant amplitude range σ L H. In that case, the amplitude of ψ would quickly decay (due to the kinematics of a massive field in an expanding Universe) and ψ would not be able to source ζ any more. (3) The third possibility is that λ effectively vanishes before the end of inflation [before even ∆V (ψ) changes]. Here, even though the amplitude of ψ has not decayed, the sourcing offered by ψ ends.
In all of the previous cases, the non-Gaussian statistics of ζ persists, simply because on superhorizon scales ζ remains constant (after ψ has done its job of sourcing its statistics). In other words, the statistics transferred by ψ while λ = 0 and ψ = 0 serves as the initial condition for a posterior phase where λ = 0 and/or ψ = 0. Thus, because the statistical transfer is cumulative, the new phase with λ = 0 and/or ψ = 0 would not imply that the non-Gaussian statistics of ζ is erased. All the contrary, if λ = 0 or ψ becomes massive, then ζ would kinematically decouple from ψ and would continue to evolve independently, with a frozen amplitude, preserving its non-Gaussian statistics. Of course, the statistics of ζ would then survive reheating until horizon entry, fixing the initial conditions for perturbations in the hot Big-Bang era.

D. Current constraints on the PDF
We can constrain the level of non-Gaussianity in the probability distribution function (106) by looking into current bounds on the trispectrum [21,[98][99][100][101][102] set by Planck, as this model has an identically vanishing bispectrum [see Eq. (70)], and consequently, we cannot use it to constrain the resulting PDF.
Specifically, Planck is able to constrain the parameter g local NL that appears in the following relation involving the four-point function for ζ, and its power spectrum: This expression may be compared with our general expression (70) for the specific case n = 4, which is given by To compare both expressions, it is necessary to recall, from the discussion around Eq. (9), that the power spectrum for ζ is given by Then, it follows that g local NL is given by We can turn this expression into a more useful result by recalling, from Eqs. (89) and (90), that σ 2 ζ = σ 2 L λ 2 ∆N 2 /2 and f ζ /σ ζ = f /σ L . We obtain With the help of Eq. (83) we see that σ 2 ζ is related to the power spectrum (113) through (recall that ξ = q L /q IR ). Planck observations [103] currently constrain the amplitude of the power spectrum as P ζ (k)k 3 /2π 2 = (2.196 ± 0.158) × 10 −9 . Then, by setting ln ξ = 8, the range of scales corresponding to the CMB, we may write σ 2 ζ 1.3 × 10 −7 . Then g NL becomes Furthermore, current constraints on the primordial trispectrum by Planck [14] imply g NL = (−9.0 ± 15.4) × 10 4 at 95% C.L. It then follows that A 2 and the ratio f 2 ζ /σ 2 ζ must satisfy the following restriction: (118) Figure 8 shows the allowed values for the parameter space spanned by A 2 and f ζ /σ ζ . It may be seen that A 2 becomes less constrained for both, large and small values of f ζ /σ ζ . It is interesting to note that the values used to plot both Fig. 6 (f ζ /σ ζ = 5 × 10 −2 and A 2 σ 2 ζ /f 2 ζ = 10) and Fig. 7 (f ζ /σ ζ = 2 × 10 −2 and A 2 σ 2 ζ /f 2 ζ = 62.5) are well within the allowed region. However, it is hard to conceive that peaks in the non-Gaussian PDF larger than those shown in Fig. 6 are not excluded. This suggests that even strong constraints on the four-point function would not compete with other methods aiming to constrain the shape of the probability distribution function.

E. Concluding remarks
To conclude, we have studied the generation of a novel class of non-Gaussianity, in multifield inflation with the help of a well-motivated model involving an axionlike field coupled to the comoving curvature perturbation. We showed that the isocurvature fluctuations can traverse the periodic potential by fluctuating across the barriers at the time of horizon crossing. As a result, the isocurvature field acquires rich statistics that reflect the landscape structure. These statistics may be transferred to the adiabatic mode via a derivative coupling which is the lowest possible coupling that arises from an effective field theory point of view (a concrete UV completion of this effective field theory can be found in Appendix A).
The outcome of nonperturbative resummations is a probability distribution function for the curvature perturbation which consists of a Gaussian part corrected by a term involving information about the landscape potential. In such a case, the traditional perturbative f NL parametrization cannot probe the full non-Gaussian structure and new estimators probing the whole PDF of temperature anisotropies should be developed.
Our study gives an example where a stringy landscape may have calculable effects on the CMB and large-scale structure that have not yet been thoroughly searched for in the current data. Last but not least, our result could have important implications for understanding the generation of primordial black holes [104].
−2/R 0 . With this form of the metric, the background equations of motion are found to be given by Let us now specify the potential V (X , Y) by splitting it as The term V 0 (X , Y) corresponds to the main contribution driving inflation, whereas ∆V (Y) will give rise to the selfinteraction of the fluctuation ψ. We need to set V 0 (X , Y) appropriately so that the action for the perturbations reduces to (11). We choose V 0 (X , Y) by demanding that it has the structure shown in (A2), with a function W = W (X ) (a function of X only). That is, we write Now, it is a simple matter to show that the system admits the following background solution: as long as the constant value of Y coincides with a local minimum of ∆V (Y), such that ∆V (Y) = 0, and That is, the field X evolves and drives inflation while Y stabilizes to a local minimum of ∆V (Y). This background solution is characterized by an inflationary trajectory in field space with tangent and normal vectors T a and N a given by Then, one may verify that the covariant time variation of T a is linked to N a through the following equation where the covariant time derivative D/dt is defined to act on an arbitrary vector V a as DV a /dt =V a + Γ a bc V bφc . The quantity Ω corresponds to the rate of turn (angular velocity) characterizing the bend of the inflationary trajectory. In the present case, it reads Thus, we see that in this model the background dynamics consists of an inflationary background with a nonvanishing rate of turn Ω. To have the appropriate amount of inflation, one needs to choose W (X ) appropriately, but the details of this choice are not relevant for this discussion.
One may now perturb the system as φ a (x, t) = φ a 0 + δφ a (x, t). In other words, we may write where X 0 (t) and Y 0 are the background fields satisfying (A11)-(A14).
To find the perturbed system we choose the comoving gauge, whereby the fluctuation along the inflationary trajectory is set to vanish T a δφ a = 0. This is equivalent to δX (x, t) = 0. In this gauge, the relevant fluctuations are the field fluctuation normal to the trajectory N a δφ a = ψ and the comoving curvature perturbation ζ, which is introduced through the following perturbed metric written with the help of the ADM decomposition: (A20) where N and N i are the usual lapse and shift functions. It is now straightforward to deduce the action for the fluctuations ζ and ψ. If one keeps ζ and ψ quadratic everywhere in the action, except for ∆V , where one keeps ψ to all orders, then one obtains the Lagrangian of Eq. (11). The parameter α of Eq. (11) is related to Ω by As a final step, the potential (13) is obtained for the particular choice The reason why we have kept ψ up to quadratic order everywhere in the action except for ∆V is that we are interested in treating the range of field fluctuations in ∆V nonperturbatively. In the case of (A22), this is equivalent to treating 1/f nonperturbatively.
The subscript c reminds us that we are keeping fully connected contributions only. As discussed after Eq. (69), to obtain these contributions, we must keep only contractions including at least one operator a + (or a † + ) appearing in H Λ I . To track these contractions, we have introduced the sum Ip over the sets I p = 1 p , 2 p , · · · , consisting of all possible arrangements (of dimension p) of momenta labels. For example, we could write: 1 2 = (1, 2), 2 2 = (1, 3), 3 2 = (2, 3), etc.
To perform the contractions appearing in (B1), it is useful to define the following two quantities: In terms of these quantities, one can write the following field commutators: These relations allow us to further infer the form of the commutators involving H λ I appearing in (B1). These are found to be given by Then, using these back into (B1), and performing the relevant contractions, we find if n is even, and zero otherwise because the expectation value of an odd number of fields in the interaction picture vanishes identically. This is a consequence of the potential's being an even function of the isocurvature field. In Eq. (B8), the sum over "perm q l " refers to the sum over all possible permutations l → q l of labels belonging to the set I p . These permutations affect the arguments τ q l appearing in some functions in the second and fourth lines of this equation.
The previous result can be simplified with the help of the following two steps: first, the sum appearing in the third line comes from the infinite loop contributions shown in Fig. (4). They are the consequence of contractions between pairs of creation and annihilation operators a † + and a + appearing in H Λ I . These terms can be resummed back into the following exponential: where σ 0 is nothing but the variance of ψ already defined in Eq. (71). Notice that σ 0 is independent of time, and so we may factorize exp(−σ 2 0 /2f 2 ) out of the τ integral. Secondly, in Eq. (B8) we may relabel every integration variable of the form τ i by a new variable τ l (with l ∈ I p ) in such a way that, in the functions' arguments, the k l 's are always accompanied by a τ l . This relabeling allows us to recognize that the sum over all possible permutations q l becomes a sum over all domains of integration for the variables {τ l } l∈Ip . As a result, the nested integrals of Eq. (B8) unravel, and we are led to u(k 1 , τ ) · · · u(k n , τ ) c = i n+1 (2π) 3 δ (3) The previous expression can be simplified further by performing the summation over the index p (including the sum over the sets I p introduced earlier). Notice that these sums allow one to rewrite the second and fourth lines of Eq. (B10) as the multiplication of pairs of terms where Θ(τ ) is the usual Heaviside function. The functions G ± (k l , τ , τ ) are nothing but the mixed propagators defined in Ref. [106], with the external -or boundary, in the language of [106] -time, denoted by τ here, left explicit as an argument. Moreover, subtracting them we find Replacing these definitions back into Eq. (B11) and noting that consecutive terms with alternating signs cancel out in the sum over j, we end up with the following simple expression for the n-point correlation function It is rewarding to notice that this result is exactly what we would have obtained had we used the diagrammatic rules of [106], after adequately treating the loop contributions.
Let us now perform the integrals of Eq. (B16) [or, equivalently, Eq. (B11)] to obtain a simple and useful expression for the n-point correlation function in momentum space. As discussed in Sec. IV, recall that the effect of the H λ I is to source the evolution of the amplitude of ζ (or, equivalently, u) on superhorizon scales. Thus, for a given fixed set of k l 's, the integration domain of each integral appearing in Eq. (B16) can be split into two parts. Before horizon crossing, where |τ l |k l > 1 (and |τ |k l > 1) and after horizon crossing, where the opposite is true: |τ l |k l < 1 (and |τ |k l < 1). The integrants are highly oscillatory in the domain |τ l |k l > 1 (and |τ |k l > 1). These oscillatory contributions would have vanished had we kept track of the prescription. And so, we may simply disregard the contributions to (B11) coming from these domains.
On the other hand, the integration of these functions over the domain |τ l |k l < 1 (and |τ |k l < 1) gives us expressions that dominate if the upper limit τ is such that |τ |k l 1 (in fact, the integration over the domain |τ l |k l < 1 diverges as τ → 0). Thus, instead of obtaining exact expressions for these integrals, we may seek the infrared divergent contributions that dominate on superhorizon scales. To do this explicitly for each integral, we introduce an arbitrary time τ 0 such that |τ 0 |k l < 1 for all k l 's. We will use τ 0 to cut off the lower limit of every time integral. The upper limit may be either the end of inflation or some other value, depending on the exact mass of the scalar field. In either case, we assume that the amount of e-folds that different modes spend outside the horizon is approximately the same, as discussed in the main text. ln(τ 0 /τ ) as the number of e-folds after all the modes left the horizon: ∆N ln(τ 0 /τ ). (B26) Then, it is possible to verify that the infrared contributions to the integrals leading to (B24) dominate when the condition is satisfied. Recall that, in order to deal with the system perturbatively, our computation assumed λ 1. This means that only after a time ∆N ∼ 1/λ the effects (both linear and nonlinear) induced by ψ start to take over, as evidenced by the computations of Sec. IV.
neglect the correction to the PDF, at the very least we must have This result severely limits the applicability of the PDF of Eq. (104), and forces us to consider the more general reconstruction carried out in Sec. VI E.

Derivation of ∆In
To show Eqs. (C11)-(C13), we may split the integral (C3) into two parts, one containing the term sin(x) − x cos(x) and the other containing the term − sin(x/ξ) + (x/ξ) cos(x/ξ). Then by changing the integration variable of the second part as x → x/ξ, we end up with the following expression I n (ξ) = nπ 2(2π 2 ) n+1 Ī n (ξ) + (−1) nĪ n (1/ξ) , where we have definedĪ n (ξ) as (C16) Notice that for large ξ the second contribution to (C15), given byĪ n (1/ξ), is very suppressed compared to the first one, given byĪ n (ξ). To see this, it is enough to see that F (1/ξ, x) is a function that becomes quickly independent of ξ for values x > 1/ξ, and so one finds that in the relevant integration domain the function is essentially given by This implies that the contribution to (C18) coming from I n (1/ξ) is at most of order 1. Then, we are left with I n (ξ) = nπ 2(2π 2 ) n+1Ī n (ξ) + O(1). (C18) Next, notice that the function F (ξ, x) that appears inside the integralĪ n (ξ) satisfies F (ξ, 0) = ln 1 ξ .
In addition, in the range 1 < x ≤ ξ, it is well approximated by F ξ (x) 1 − γ + ln ξ x , whereas, for values x > ξ, the function F ξ (x) vanishes quickly. We can therefore write where γ is the Euler-Mascheroni constant and (y) is a slowly varying function that remains small in the relevant interval 1 < x ≤ ξ. In fact, this function is suppressed everywhere 1 < x ≤ ξ and its largest value is of order 0.1 when x/ξ = 1. It is enough to take (y) = y 2 /12. Then, we can writē where x * ξ has been introduced to cut off the integral, since for values larger than x * the function F (ξ, x) is highly suppressed. The introduction of x * has the benefit of allowing us to use approximation (C20) inside the integral. Now, taking a derivative with respect to ξ, we findĪ n (ξ) = (n − 1) Because ξx (ξx) is of order in the entire domain 1 < xξ, the second term is clearly subleading with respect to the first one. Then, we can simply disregard the second term, and writeĪ n (ξ) = (n − 1) ξĪ n−1 (ξ). (C24)