Stochastic Inflation at NNLO

Stochastic Inflation is an important framework for understanding the physics of de Sitter space and the phenomenology of inflation. In the leading approximation, this approach results in a Fokker-Planck equation that calculates the probability distribution for a light scalar field as a function of time. Despite its successes, the quantum field theoretic origins and the range of validity for this equation have remained elusive, and establishing a formalism to systematically incorporate higher order effects has been an area of active study. In this paper, we calculate the next-to-next-to-leading order (NNLO) corrections to Stochastic Inflation using Soft de Sitter Effective Theory (SdSET). In this effective description, Stochastic Inflation manifests as the renormalization group evolution of composite operators. The leading impact of non-Gaussian quantum fluctuations appears at NNLO, which is presented here for the first time; we derive the coefficient of this term from a two-loop anomalous dimension calculation within SdSET. We solve the resulting equation to determine the NNLO equilibrium distribution and the low-lying relaxation eigenvalues. In the process, we must match the UV theory onto SdSET at one-loop order, which provides a non-trivial confirmation that the separation into Wilson-coefficient corrections and contributions to initial conditions persists beyond tree level. Furthermore, these results illustrate how the naive factorization of time and momentum integrals in SdSET no longer holds in the presence of logarithmic divergences. It is these effects that ultimately give rise to the renormalization group flow that yields Stochastic Inflation.

Stochastic Inflation [21][22][23] is a framework for treating the IR dynamics of a massless scalar field in a dS background, and has long been suspected to provide the nonperturbative resolution to the IR issues associated with massless fields in dS [24][25][26][27][28][29][30][31][32][33][34][35]. The idea is to reframe the problem in terms of the probability distribution for the scalar field as a function of time, resulting in a Fokker-Planck equation that depends on the scalar field potential. There are two contributions to the evolution, resulting from quantum noise induced by fluctuations of the field as it crosses the dS horizon and classical drift due to the potential. This framework forms the conceptual basis for slow-roll eternal inflation [36][37][38] and can be used to describe the onset of eternal inflation quantitatively [39]. Moreover, these results hint at the physical meaning of the dS entropy [40][41][42], which remains a significant unsolved problem [43].
Notwithstanding the conceptual and technical appeal of Stochastic Inflation, it is necessarily approximate. For example, we expect that there are non-Gaussian contributions to the quantum noise that result from the UV interaction, which are not modeled by Stochastic Inflation. Furthermore, the Fokker-Planck formalism obscures the connection to cosmological correlators, and it is not a prioi obvious how to incorporate higher-order corrections. It would be ideal if we could understand how the success of Stochastic Inflation relates to other results regarding the IR behavior of fields in dS, such as the freeze-out of superhorizon metric fluctuations which has been shown to all orders in perturbation theory [44][45][46], or the loop generated anomalous scaling for the time-evolution of massive fields [47,48]. One of our goals in exploring the corrections to Stochastic Inflation is to understand how they fit into the broader context of quantum field theory in dS.
A framework that accomplishes this ambitious goal is the Soft de Sitter Effective Theory (SdSET) [33]. By following the standard Effective Field Theory (EFT) playbook, this approach isolates the dynamics that persist in the superhorizon limit, yielding more efficient calculations of loop corrections to long wavelength cosmological correlators. Taking a real scalar field in dS as the UV description, SdSET describes the IR limit of this model by relying on two degrees of freedom that correspond to the growing and decaying modes which are familiar from the solving the Klein-Gordon equation classically in a dS background. This representation admits a power counting prescription that systematically expands about the long wavelength limit in terms of a local Lagrangian. Loop dependence on time and space is manifestly factorized throughout the calculation, allowing an efficient isolation of the time dependent IR divergent logs. Such logs lead to secular growth for both massive and light fields, and appear in SdSET as contributions to the anomalous dimensions of local operators. In the case of light fields, an infinite number of operators become degenerate and Starobinsky's model of Stochastic Inflation is equivalent to the leading order (LO) dynamical renormalization group (RG) that governs their mixing as a function of time. 1 This implies that corrections to Stochastic Inflation can be computing by simply extending the RG analysis to higher orders.
Taking the UV description to be massless λφ 4 theory, the endpoint of this RG flow is a non-trivial fixed point where the field values are φ ∼ Hλ −1/4 . Corrections to this description around this fixed point must account for this non-perturbative scaling with λ. In this paper, we will calculate the evolution of operators to next-to-next-to leading order (NNLO) in this power counting. At NLO, our results reproduce previous calculations [30,34]; as we will show, these contributions can be attributed to field definitions within SdSET. In contrast, at NNLO we find a universal correction in the form of a two-loop anomalous dimension that introduces the first higher derivative correction to Stochastic Inflation. In the process, we perform the full one-loop matching in SdSET, which further elucidates the relationship between the EFT and the UV descriptions.
One novel feature of SdSET is that consistently matching a UV theory onto the EFT requires specifying both Wilson coefficients and (time-independent) initial conditions. Deriving the RG that yields Stochastic Inflation at NNLO requires performing this matching explicitly at one-loop order. This provides a highly non-trivial check of the SdSET formalism, and these results can be utilized for a wide variety of correlator calculations. We will also use this calculation as an opportunity to demonstrate the power of the symmetry preserving "dynamical dimensional regularization" technique introduced in [33]. This paper is organized as follows. We begin with a review of Stochastic Inflation in Sec. 2, with an emphasis on its origins as a Markovian process, which provides a framework with which we can organize corrections. Then Sec. 3 reviews the most salient aspects of the SdSET formalism. The new calculations begin in Sec. 4, where we present the oneloop matching results that are relevant for our applications here. These are then applied in Sec. 5, where we compute the composite operator anomalous dimensions that feed into Stochastic Inflation up to NNLO, and leads to the main result of this work in Eq. (5.33). We then explore the implications of this formula in Sec. 6, and finally conclude in Sec. 7. An appendix on the relevant, but somewhat technical, six-point function matching is provided in App. A, and the hard cutoff version of the main calculations are given in App. B.
Guide for the reader: We have attempted to make this paper accessible to a wide ranging audience. For the casual reader, we suggest (i) reading Sec. 2 and Sec. 3.5, (ii) studying the matching diagrams in Sec. 4, and (iii) studying the composite operator mixing diagrams in Sec. 5. This brings the reader to (iv) the nearly final result in Eq. (5.25), where the explicit corrections are given in Eq. (5.26). Finally, (v) Sec. 5.3 provides the derivation that results in the simplified form of the NNLO formula for Stochastic Inflation given in Eq. (5.33). One might also be interested in the NNLO equilibrium probability distribution and relaxation eigenvalues presented in Sec. 6, which brings the reader to the conclusions in Sec. 7.

Stochastic Inflation
In a theory of a massless scalar field φ in a dS background with Hubble constant H, the field's value will fluctuate by O(H) as each momentum mode crosses the dS horizon. An equivalent point of view is that these stochastic fluctuations are the result of the non-zero temperature within dS. This effect has a natural interpretation as a random walk, an idea that was made precise by Starobinsky [21] and led the formalism known as Stochastic Inflation. In this section, we discuss this approach, and emphasize the structure of higher order corrections.
For concreteness, we will assume the canonical example of λφ 4 theory in a dS background, whose UV description in terms of a scalar field φ is where g µν is the dS metric and g ≡ det g µν as usual. The essential formalism developed here holds for general models. However, we will not be able to derive corrections to Stochastic Inflation generically, and so we will work with this simple and well studied example when we calculate explicit corrections.
In the process of discussing the general structure of corrections to Stochastic Inflation, we will arrive at a natural interpretation for higher-order corrections in terms of the transition amplitudes for the field φ. Unfortunately, how to determine the corrections directly is not transparent in this description. In Sec. 3.5 we will show how the corrections discussed in this Section arise from operator mixing, and the remainder of the paper is devoted to deriving these corrections and their implications using SdSET.
Before moving on, we emphasize that the formalism we develop in this section will rely on the assumption that the late time evolution of φ can be modeled as a Markovian system (as described in Sec. 2.2 below). This will be justified by the concrete calculation of the dynamical renormalization group using the SdSET that is developed later in this paper. It ultimately is due to the fact that the dynamics of the SdSET degrees of freedom are governed by a first order equation, which is a consequence of the EFT power counting.

Leading Order
The framework of Stochastic Inflation results in a probability distribution P (φ, t) for the field φ at a time t. To leading order, P (φ, t) obeys a Fokker-Planck equation where the first term captures the stochastic noise from the inherent quantum variance of φ, while the second term is due to the classical drift induced by the potential where V (φ) ≡ ∂V /∂φ. One interesting application of this equation is to solve for the fixed point that the scalar field would reach if it lived in an eternal dS background. To find the fixed point, we enforce that ∂P eq /∂t = 0 for the equilibrium solution P eq , which implies Integrating both sides of this equation twice leads to the solution We can use this leading order solution to organize corrections to Stochastic Inflation as a perturbative series in the UV coupling.

Beyond Leading Order
Having reviewed the leading order formalism and its consequence for pure dS, we now turn to exploring the form we can expect corrections to take. In order to generalize this Fokker-Planck equation, we return to its origins. The underlying assumption is that the system is Markovian, in that the time slice of interest is entirely determined by the information contained in the previous step. In other words, a Markovian system has no "memory." Since this assumption holds for the spectrum of scalar field fluctuations at horizon crossing, the resulting formalism will tell us what kinds of corrections to Stochastic Inflation we can expect. The Markovian assumption leads directly to the Chapman-Kolmogorov equations, which describe a probability distribution P (φ, t + dt) that is fully determined by P (φ, t): where W (φ|φ ) is a "transition rate" in that it sets the rate for transitioning to φ from ∆φ ∆φ Figure 1: Visualization of the Kramers-Moyal expansion. The top panel shows the probability of "hopping" from φ to φ , W (φ |φ) or from φ to φ, W (φ|φ ), or equivalently any other point such as φ . On the bottom, we see the process in terms of the probability W (∆φ|φ) to hop from a specific starting point φ by a distance ∆φ.
another value φ in a differential amount of time. This equation simply expresses that the probability distribution for φ at t + dt is fully determined by the weighted sum of the possible transitions that yield φ minus the sum of all the weighted transitions for φ to change value. Next, we will reorganize Eq. (2.5) using what is known as the Kramers-Moyal expansion, visualized in Fig. 1. This is effectively the assumption that the transitions are dominated by "local" jumps [52]. The first step is to make a substitution of ∆φ = φ − φ in the first term and ∆φ = φ − φ in the second term. In other words, when φ is in the final state, then φ = φ − ∆φ and when φ is in the final state φ = φ + ∆φ. This yields where we have included two compensating relative minus signs, one from the different changes of variables for the two terms, and another due to needing to flip the limits of integration after switched the integration variable ∆φ → −∆φ in the second term. Since we are assuming the main support of W comes from local jumps, we want to Taylor expand the first term for fixed φ. To make performing this expansion transparent, it is then useful to redefine W using so that Eq. (2.6) becomes Then Taylor expanding the first term about a fixed value of φ yields where in the last line, we used the fact that ∆φ is independent of φ to pull the derivatives outside of the integral. Plugging this expansion into Eq. (2.8), we see that the n = 0 term cancels so that where and we have used Eq. (2.7) to write this expression in terms of the original transition rate W . Note that all terms in this expansion are total derivatives, as required for the conservation of the total probability. We also see that Ω n (φ) encodes the φ-dependence of the n th moment of the distribution. For the theories of interest here, Ω n (φ) will admit a polynomial expansion in φ: Hence, for n = 1, m tracks the polynomial interactions that could appear in a generic potential. Moving to the n = 2 "noise" term, we again can compare to Eq. (2.2) to find In this case, the m > 0 terms correspond to higher order corrections.
To summarize, if we assume that the UV theory has only a λφ 4 interaction, we conclude that the generalized evolution equation that describes Stochastic Inflation takes the form where V is the φ derivative of the potential, which includes the matching corrections required to obtain the accuracy of interest.
To develop some intuition for what the Ω n corrections are capturing, we can interpret W (φ + ∆φ|∆φ) as the probability distribution of transitions of size ∆φ. Then Ω n (φ) is simply the n th moment of this distribution; the first and second moments Ω 1 (φ) and Ω 2 (φ) are the complete set of inputs for a Gaussian distribution. Furthermore, if Ω 2 has non-trivial φ dependence, i.e., Ω (m =0) 2 = 0, then the variance of the noise depends on the starting location of the jump. Finally, if Ω n=3 (φ) = 0, then we know our distribution is non-Gaussian. More generally, we interpret the n > 2 terms in the generalized equation governing Stochastic Inflation in Eq. (2.16) as encoding contributions from non-Gaussian noise generated by the UV λφ 4 interaction. Therefore, we conclude that corrections to Eq. (2.2) are of three types: • Higher order "noise" terms captured by ∂ n /∂φ n with n > 2, see Eq. (2.11).
• Higher order terms in the potential V via corrections to coefficients and the generation of higher polynomial φ terms.
Next, we will argue for how to relate the expansion in each of these quantities to the expansion in the UV quartic coupling λ as it corrects the equilibrium solution in Eq. (2.4).

Organizing Corrections Systematically
Due to the underlying λφ 4 potential, we expect that the statistics of φ are neither Gaussian nor independent of the background value of the field. We therefore expect corrections to the equation for Stochastic Inflation of the form discussed previously. In this section, we will take the equilibrium solution in Eq. (2.4) and apply it to a UV theory with where φ 4 eq = 9H 4 /(π 2 λ). This distribution has support over a field range |φ| φ eq such that Therefore, we will organize the possible corrections by assuming the equilibrium scaling We will further assume that the corrections are generated as an expansion in perturbation theory, so that which we will see agrees with the explicit calculations presented below. Putting all of this together allows us to determine the order in λ for each term that appears in Eq. (2.16).
Note that we will assume the φ → −φ UV symmetry is preserved, which explains the absence of many terms. These contributions are as follows.
Leading Order (LO): Stochastic Inflation at leading order is given by Eq. (2.2): where we have used the known leading order results Ω (0) 2 = H 3 /(4π 2 ) and V (φ) = λφ 3 /3!. Both terms on the right hand side of this equation are O λ 1/2 × P (φ, t), which defines what we mean by "LO." Next-to-leading Order (NLO): Accounting for both Ω (m) n and corrections to the potential, we can determine that the next-to-leading corrections to Stochastic Inflation should take the form where the correction to the noise term Ω 2 ∼ λ and the correction to the potential c 6 ∼ λ will both be determined below, and Ω (1) 2 = 0 due to the φ → −φ symmetry. We see that these NLO terms are O(λ) × P (φ, t). These corrections have been previously calculated in [30,34].
Next-to-next-to-leading Order (NNLO): Following the same logic, we can find the form that the next order terms take: where Ω 2 ∼ λ 2 and Ω 3 ∼ λ will be determined by operator mixing in the next section and c 8 ∼ λ 3 , and Ω (0) 3 = 0 due to the φ → −φ symmetry. These NNLO terms are Note that in addition to these corrections, we must also include subleading corrections to the parameters that already appear at lower order; these do not change the structure of the equation, but will of course be accounted for as we perform the calculation. The rest of this paper is devoted to determining these coefficients systematically using the framework of SdSET, with a brief discussion of the physical implications of working with NNLO Stochastic Inflation.
Soft de Sitter Effective Theory offers a method to compute the equations of Stochastic Inflation systematically to any order. The key advantage offered by SdSET is that power counting is manifest, thereby making corrections easy to identify. Furthermore, loop integrals are scaleless and regulated by (dynamical) dim reg and thus preserve the power counting that is manifest in the action. In addition to fixing the values of the SdSET Wilson coefficients, the UV theory sets the initial conditions for the effective theory fields. We will be specifically interested in understanding the corrections to Stochastic Inflation for the concrete example of λφ 4 theory in dS, where the UV action is given above in Eq. (2.1).
In this section, we will review the machinery of SdSET and how it arises from a given UV theory. In the subsequent sections, we will use this technology to match SdSET to λφ 4 theory and then use it derive the equations of Stochastic Inflation at NNLO.

In-In Correlators
As we will see below, Stochastic Inflation is equivalent to the renormalization group equations that govern how composite operators mix. One approach to determining the operator mixing is to compute the divergences of in-in correlation functions involving composite operators. This section is devoted to setting up the relevant framework. We work in the interaction picture, where fields are quantized using the solutions to their quadratic equations of motion. A free scalar fields in dS can be expressed as a mode expansion: where τ = −1/[aH] is the conformal time, and a † k and a k are the canonical creation and annihilation operators respectively that satisfy In the Bunch-Davies vacuum, one finds the positive frequency modes are given bȳ so that ν = 3/2 corresponds to a massless field (see e.g. [53] for review). The observables of this theory are equal-time in-in correlation functions, which are computed via We will be interested in multi-field correlators in the long wavelength limit, so that Q(t) = φ k 1 , t ... φ k n , t . In general, this expression must be generalized to allow for the i prescription that projects the initial state onto the interacting vacuum, but this form is the most useful starting point for understanding the superhorizon evolution. For illustration, we can compute the tree-level four point correlation function as where the φ k fields are all evaluated at the same time t, and we have introduced the notation At tree-level, we can simply use the massless mode functions, 2 8) to evaluate this expression for a real operator Q(t).

Taking the Long Wavelength Limit
In this section, we review how to determine in-in correlators in the soft limit using SdSET. The starting point is to decompose the UV fields according to is split into soft (superhorizon) and hard (subhorizon) modes, and we have introduced a dimensionless time variable t ≡ Ht so that the mass dimension of operators tracks the EFT power counting. The parameters α and β are subject to the constraint α + β = 3. This is straightforward to derive from the top down, as these parameters are determined by the mass via α = 3/2 − ν and β = 3/2 + ν, where 2 We will use 3 2 − ν ≡ α = 0 as a regulator when we encounter loops.
ν is defined in Eq. (3.3). The decomposition of φ S into ϕ + and ϕ − is exact in the free theory. Taking the limit k τ 1 and using τ = −1/[aH] in Eq. (3.3), we find and (3.14) The operatorsã k andb k are given in terms of the UV creation and annihilation operators of the formã From the UV theory, we determine that the operators commute with themselves Nevertheless, these operators still have non-zero correlation function where .. ≡ 0|..|0 , and |0 is the vacuum that is annihilated by a k . This gives rise to classical statistical power spectra where C α and D β are defined in Eq. (3.14) above. Note that in the massless limit α → 0, we reproduce the famous scale invariant power spectrum. The corrections to this mapping can be systematically accounted for by matching between the UV theory and the EFT, see Sec. 4. The fields ϕ + and ϕ − have well defined power counting; they carry operator dimension α and β respectively. After utilizing field redefinitions, on-shell conditions, and power counting to remove redundant operators, the low energy effective action is given by where the c n,1 are dimensionless Wilson coefficients. Note that t carries dimension zero by SdSET power counting, so marginal operators are dimension three. This explains why we have only included operators with a single factor of ϕ − since these are the only terms that become marginal in the massless limit (α → 0). Additionally, we have not included any terms with ∂, which start at dimension five and are therefore power suppressed by at least In addition to an action, SdSET requires specifying initial conditions for the fields ϕ + and ϕ − that are acquired from the time evolution prior to horizon crossing. These initial conditions are random such that, to leading order, ϕ + behaves as a classical stochastic variable with correlations fixed by matching where K is a reference momentum scale, F (n) { q i } encode the dependence on the rescaled momenta q i = k i /K, and the (n) subscripts track the order in the λ perturbative expansion for each contribution; the two point correlators in Eq. (3.18) should be viewed as ... IC (0) . Because ϕ + is time independent to leading order in the EFT, the initial conditions are determined by matching the time independent terms. We evaluate time integrals in the EFT using where we assume that this holds even when γ < 0. This analytic continuation enforces that the contributions from early times vanish, thereby ensuring that power law divergences associated with physics at horizon crossing are automatically absorbed into the initial conditions (in close analogy with how dim reg treats power law divergences).
Cosmological correlators are determined in SdSET using the same in-in formalism as applied to the UV theory, see Sec. 3.1. For illustration, we can compute the tree-level trispectrum using Eq. (3.4) and the canonical commutator Performing the time integrals using Eq. (3.21), we find where c 3,1 is the Wilson coefficient for the ϕ 3 + ϕ − operator, see Eq. (3.19), and we used Eq. (3.18a) to evaluate the field contractions. In addition, we must include any trispectrum associated with the initial conditions.

(Dynamical) Renormalization
Loops corrections are calculated in the SdSET using dynamical dimensional regularization (dynamical dim reg). Rather than varying the spacetime dimension, we instead float the dynamical exponents α, and evaluate loop integrals by analytic continuation in α. Then when we encounter divergences as α → 0, they will be accompanied by log corrections to the time evolution, in exact analogy with conventional dim reg. To keep the units fixed as we vary α, we will introduce the necessary powers of [aH] such that ϕ + stays dimensionless. Then we take α = 0 at the end of the calculation, so that: so that ϕ + corresponds to a massless mode. Since we are working within the EFT, we will typically encounter vanishing scaleless integrals. Then we can isolate the UV divergence in the usual way by regulating the IR with a dimensionful parameter K: where we have taken the limit α → 0 in the third line. Having isolated these UV divergent contributions, we can use them to determine the (dynamical) RG flow, which resums a series in log[aH].
Having regulated the divergence, we can then absorb it into the renormalization of the operator where [aH] and K are energy and momentum scales we have invented to make the logs small, i.e., subtraction points. By definition the bare operator O is independent of these arbitrary parameters, and thus where the anomalous dimension γ is independent of [aH]. We are working with a scheme where Z O is diagonal at the scale [aH] . In general, γ → γ ij is a matrix that acts on the space of operators, and which encodes both the anomalous scaling and mixing of these operators.

Matching and Initial Conditions
SdSET provides an effective description for the time evolution of scalar modes that have crossed the Hubble horizon. Their state at the time of horizon crossing cannot be computed within the EFT, and instead has to be provided as an additional input. This is why SdSET requires matching for both the initial conditions and the EFT Wilson coefficients. This is not unique to SdSET, but is necessarily part of any EFT description of the post inflationary universe as well, see e.g. [54]. When defined using dynamical dim reg, SdSET is a so-called continuum EFT [55]. Concretely, the time integrals include arbitrary early times, even though the EFT does not provide a model of the subhorizon physics, since it relies on the use of the long-wavelength mode functions at all times. Importantly, as with all continuum EFT, these early time integrals only make scaleless (and therefore vanishing) contributions. Thus we can integrate over all times, such that the regulated integrals respect the low energy symmetries and the EFT power counting. Underlying the validity of this procedure is the fact that the subhorizon physics only alters the initial conditions and thus we can fully account for all subhorizon evolution by matching a UV theory onto SdSET.
One can demonstrate how matching separates into Wilson coefficient and initial condition corrections more directly using a hard cutoff to evaluate time and momentum integrals, i.e., treating the theory as a Wilsonian EFT [55]. Let k i denote the magnitudes of the momenta appearing in the cosmological correlator, and let the cutoff of the momentum integrals Λ be much greater than any of the k i . In addition, we denote the time cutoff by t Λ ; this is the time when all the EFT modes are in the superhorizon limit, and thus it specifies when we set the initial conditions for the EFT.
Before t Λ , a subset of EFT modes are subhorizon, and so one must use the UV theory to describe them. To match onto the EFT, it is useful to split the full theory time evolution into pieces before and after t Λ . Let U I (t, t ) represent the interaction picture propagator, and |Ω be the UV vacuum state. This decomposition can then be written as where φ I is the interaction picture field. One can trivially re-write Eq. (3.29) as an expectation value of |ψ = U I (t Λ , −∞) |Ω , the state of the full theory fields at t Λ . One can then integrate out the modes whose wave vector magnitudes satisfy k > Λ, so that the remaining modes are superhorizon after t Λ . It was shown in [33] that, after integrating out these so-called hard modes, the resulting action for the superhorizon modes ϕ + and ϕ − is local. We can then evolve these modes from t Λ to t using the unitary time evolution operators defined within the SdSET itself: where U I,EFT is the interaction picture propagator obtained from Eq. (3.19), and φ S is given by Eq. (3.24). The state |ψ EFT (t Λ ) is the EFT state inherited from |ψ that results from integrating out the hard modes; this state encodes the initial conditions for ϕ + and ϕ − . One then fixes the EFT parameters in Eq. (3.19) and the initial state |ψ EFT (t Λ ) by matching to full theory correlators. In practice, the full form of |ψ EFT (t Λ ) is more infor-mation than is needed to derive the late time behavior of n-point correlation functions of φ. Instead, it is sufficient to determine a finite number of n-point functions, from matching. Furthermore, this shows that all of the contributions from t < t Λ are encode in the state |ψ EFT (t Λ ) or, equivalently, in the initial conditions set at t Λ . Finally, since t Λ is an unphysical cutoff parameter, no physical results can depend on it. In particular, the initial conditions can be identified as the time-independent contribution from the UV correlators. This implies that we do not need to rely on hard cutoff to derive the initial conditions. In what follows, we will regulate the theory using dynamical dim reg, and will derive the initial conditions by identifying the time-independent contributions to correlators that appear when matching.
To illustrate the matching procedure, it is useful to first consider a simple example. By expanding Eq. (3.4) to leading order in c 3,1 and setting |in = |ψ EFT we find that the tree-level EFT prediction for the four point function of φ is As in all perturbative SdSET calculations, the time dependence and field contractions factorize at the integrand level. The EFT prediction is then parameterized by the coupling c 3,1 and two expectations values of |ψ EFT .
The expectation values parameterize the subhorizon evolution of the EFT modes before t Λ . If the full theory is perturbative, then the subhorizon evolution is approximately Gaussian. Assuming the UV theory is given by λφ 4 , then to O(λ) we can write where the final term on the RHS encodes the non-Gaussian contribution to the subhorizon evolution of the modes. There is a similar formula for the expectation value on the second line of Eq. (3.32). However, since it is already multiplied by c 3,1 , we only need the ... IC (0) contribution to the expectation value at this order.
The SdSET two point functions derived from the free theory are given in Eq. (3.18). To include the impact of the UV interaction on the EFT, we compute ϕ I,+ k 1 ... ϕ I,+ k 4 IC and c 3,1 by matching to the full theory. The superhorizon evolution between t Λ and t generates the term proportional to c 3,1 , and so we can isolate the initial conditions contribution by evaluating both sides of Eq. (3.32) at t Λ , giving where the "connected" subscript refers to the fact that this does not include the contributions from products of lower point contractions. The RHS of Eq. (3.34) can be computed using Eq. (3.4), where H I is given by the full theory interaction Hamiltonian and the time integrals extend from −∞ to t Λ . Since the integral's main region of support occurs when the modes are subhorizon, we cannot replace the mode functions with their late time behavior, and instead have to use their full UV form given in Eq. (3.3). The UV mode functions simplify for scalars whose mass is much lighter than the Hubble constant and can be approximated by Eq. (3.8).
One can then fix c 3,1 by demanding that Eq. (3.32) reproduces the full theory prediction for the correlator in the regime t > t Λ . While the split between subhorizon (t < t Λ ) and superhorizon (t > t Λ ) evolution is manifest in the EFT, this split has to be inputted by hand in the full theory, as was done in Eq. (3.29); this is effectively making a choice of scheme. Practically, one can decompose the time integrals in Eq. (3.4) into regions before and after t Λ . In our four point example, subhorizon contribution is already taken care of by the initial conditions, i.e., Eq. (3.34), while the second term in Eq. (3.32) must reproduce the superhorizon evolution. The Wilson coefficient c 3,1 is fixed by this condition.
Since the split time t Λ is arbitrary from the perspective of the full theory, all fulltheory and EFT predictions of the φ correlators must be independent of it. This provides an additional check on the matching calculation in the hard cutoff scheme. Fortunately, when we calculate with dynamical dim reg, the time integrals are manifestly independent of t Λ , while still maintaining the split between initial conditions and time-evolution. In this sense, the hard-cutoff scheme proves the validity of dynamical dim reg, while dynamical dim reg makes it manifest that we can implement this procedure without breaking symmetries. This matches the more intuitive argument that our treatment of initial conditions and time evolution in SdSET is identical to the continuum EFT approach.

Stochastic Inflation from SdSET
From the point of view of our EFT, Stochastic Inflation can be understood as a consequence of operator mixing. Specifically, for light fields for which α → 0, the composite operators ϕ n + are degenerate to leading order (in that they have the same dimension as determined by the EFT power counting). Assuming the correlations of these fields are only due to the Gaussian contribution given in Eq. (3.18), one encounters a UV divergence from a one-loop The dynamical RG associated with this operator mixing can be written as where in the second line we plugged in Eq. (2.2), Ht = t, the α = 0 relation φ| ϕ − =0 → Hϕ + , and This shows that Eq. (3.36) is equivalent to the leading order equation for Stochastic Inflation Eq. (2.2). Therefore, calculating corrections to Stochastic Inflation has been reduced to the straightforward task of computing the higher order dynamical RG equation using SdSET. Concretely, we would expect to find mixing between composite operators ϕ n + and all possible ϕ n + such that Note the role of the binomial coefficient which will originate from the number of fields inside ϕ n + whose contractions are responsible for mixing with a given operator, leading to a single log divergence. Repeating the above argument, we see that this dynamical RG is equivalent to where we see that number of derivatives is related to the binomial coefficient of the associated mixing term. Comparing with Sec. 2.2, we see that the NLO corrections are determined by c 5,1 and b 1 while the NNLO coefficients are c 7,1 , b 2 and d 0 . Achieving NNLO accuracy requires matching λφ 4 theory onto the SdSET at one loop, the subject of the next section. Note, however, that we have described Stochastic Inflation in terms of ϕ + rather than the UV field φ. This distinction will be important because the equations of Stochastic Inflation are not invariant under field redefintions. We will address these issues in detail in Sec. 5.3. Finally, we note that this result demonstrates the Markovian assumption which led to Eq. (2.16) does in fact hold as a consequence of SdSET power counting. Specifically, the fact that the ϕ − dynamics are irrelevant to the evolution of the ϕ + correlators implies that the evolution of the system is indeed linear, and thus it has no "memory."

Matching λφ Onto SdSET at One-Loop
In this section we will show how to match correlators of φ to correlators of ϕ + in the SdSET to determine the EFT parameters in terms of UV data. This then serves as input for any calculation of cosmological correlators for λφ 4 theory in the long wavelength limit, including the corrections to Stochastic Inflation we will discuss in subsequent sections. Furthermore, by extending this program to one-loop order for the first time, this calculation will serve as a non-trivial check of the SdSET framework. Unlike conventional EFTs, we match both the couplings of SdSET and the stochastic initial conditions, the former results from studying the time-dependent terms in the EFT while the latter is fixed at the time of horizon crossing. As a result, consistent matching of time-dependence of the UV correlators is non-trivial and requires that the SdSET is a complete representation of the long wavelength dynamics. In contrast, time-independent contributions to a given correlator can always be absorbed into the initial conditions (up to composite operators as discussed in Sec. 4.3). As a result, for tree and one-loop matching we will be particularly focused on time-dependent UV contributions.

Tree-level Matching and Field Redefinitions
Tree-level matching in the interacting theory is non-trivial due to the impact on the initial conditions. We will need to introduce non-Gaussian initial conditions in order to match higher point correlation functions as calculated by the UV theory.
We can understand many important aspects of matching by Taylor expanding the UV calculation in the long wavelength limit. As a simple demonstration, we can explore the superhorizon behavior of the operator φ. At first order in the coupling, we can apply the definition of the in-in correlator given in Eq. (3.4) with Q(t) = φ, which gives where φ int are the interaction picture fields. Because the time integral runs over all times, this includes both the regime where the modes are hard (UV) and the long wavelength limit where the EFT applies. Nevertheless, if we expand in the long wavelength limit and evaluate the integrals with dynamical dim reg, we will only get contributions from late times. Since φ in are the free field operators, we can use the map given in Eq. (3.11) to determine the long wavelength behavior of the full theory: 3 where in the last line we expanded in α 1. Note that the scaling dimension of ϕ + is still α and thus will provide the necessary distance scale to make the log dimensionless inside of a correlator, as is familiar from conventional dim reg. Now we turn to exploring the same effect within the EFT direction by calculating the time evolution of ϕ + using SdSET. Using Eq. (3.11) with Q = ϕ + and H int = c 3,1 ϕ 3 where it is trivial to match the tree-level UV interaction to the EFT interaction, such that The equality between c n,0 and c n−m,m found in matching is also fixed by the reparametrization invariance of SdSET. We see that the EFT is capturing the first term in the Taylor expansion of the UV theory given in Eq. (4.3), but not the second. The origin of the missing term is two-fold. First, we are only considering correlations of ϕ + instead of the full UV field, which also includes ϕ − . This alone would not matter, since ϕ − is suppressed by [aH] −3 . However, in order to organize the interactions within the EFT, we removed the c 4,0 ϕ 4 + /4! term by a field redefinition. Specifically, to remove the c n,0 ϕ n + /n! operator, we take Therefore, keeping track of the field redefinition implies that we should use with α → 0 and β → 3. Now the quantities on the RHS live purely in the EFT. As a result ϕ − will not contribute to correlation functions of ϕ because they are suppressed by powers of [aH] −3 . Combining this with Eq. (4.4), we find This result also provides the map between SdSET and Refs. [31,35], which derived the soft behavior by explicitly expanding the UV in-in correlator in the superhorizon limit. The tree-like structure they observe is a consequence of our power counting as only c n,1 is marginal and hence interactions only include a single factor of ϕ − . The nested set of commutators in Eq. (3.4) ensures that the marginal operators always have a tree-like structure. In SdSET, this is manifest from dynamical RG, and all the additional finite terms arise from the field redefinitions.
By similar considerations, we must apply the field redefinition to match the UV potential onto V (ϕ + , ϕ − ) in the EFT. Although λϕ 4 + has been removed, this procedure introduces higher order terms, such as Removing the first term will introduce a ϕ 7 + ϕ − interaction at order λ 3 and so on. As a result, our field redefinition requires that Finally, using c 2,2 = c 3,1 = c 4,0 = λ from matching, we arrive at (3!) . The correction to c 5,1 is equivalent to the NLO corrections to the effective potential calculated in Refs. [30] and [34] using complementary techniques. While these two references approach this problem from different perspectives, the wavefunction of the universe and the dS static patch respectively, both effectively integrate out the decaying mode ϕ − which leads to an additional term in the potential. Instead, when ϕ − is included, our corrections arise from insuring ϕ − does not mix with ϕ + at higher orders in perturbation theory. As the dimensions of ϕ + and ϕ − are well separate for α → 0, removing mixing can always be achieved by such a field redefinition. Additionally, it is easy to determine the correction to c 2n+1,1 ∝ λ n by repeated application of Eq. (4.6). Most importantly, we do not integrate out ϕ − to ensure we have a local action, rather than an open EFT for the growing mode alone [29].
What we have accomplished thus far is to determine the correct basis of operators to match the EFT and UV descriptions. We have ensured that the superhorizon limit of φ and ϕ agree as operators at higher orders in λ. However, in order to match the correlators of the UV theory, which include the subhorizon evolution, we will need to determine the stochastic initial conditions beyond the Gaussian limit.
In order to correctly match the four-point function, we write the EFT trispectrum to order λ as 4 where c 3,1 = c 4,0 = λ as before, and the subscript "C" denotes that this is only the connected contributions. Matching this EFT expression to the UV result in Eq. (3.9) fixes the non-Gaussian contribution to the initial conditions: where k t = k 1 + k 2 + k 3 + k 4 . Most significantly, all the time-dependence of the full UV trispectrum is already captured by the EFT and, as expected, the initial conditions are only required for matching the time-independence contributions. This result is extended to the six-point function in App. A as expected from general arguments.
Matching Derivative Operators: Before moving on the loop-level matching, let us briefly comment on the case where the UV theory is itself an effective theory. Specifically, we are only considering the case of a λφ 4 interaction in the UV, while in principle there could be a variety of higher derivative (irrelevant) interactions as well. The first such operator we can write down is (∇ µ φ∇ µ φ) 2 /M 4 . The full tree level trispectrum was calculated in Ref. [56] and is given by (4.14) What we notice right away is that this gives us a completely time-independent result. It can therefore only be absorbed into the initial conditions of ϕ + : This result illustrates a broader feature of physics in dS: higher derivatives decouple at long wavelength and thus contribute, at most, time-independent correlation functions. This property is made manifest in SdSET and must hold in matching. As a result, a variety of possible UV theories only impact the initial conditions and not the long wavelength dynamics. In this sense, our choice to focus λφ 4 is not missing more complicated superhorizon evolution. Instead, higher derivative terms are trivially matched in SdSET and thus do not further illuminate the structure of the EFT.

One-loop Matching
When calculating loops in the UV theory, rather than the EFT, our EFT regulator (dynamical dim reg) is not effective. 5 Loop calculations in the UV are not scaleless and therefore we will need to regulate them differently than the EFT approach. This difference will be absorbed into the matching calculating. In App. B, we match using a hard cutoff in both theories for an example with consistent regulators in both and find identical results, up to scheme dependent coefficients. To avoid the usual challenges of working with a hard cutoff, we will use dimensional regularization in the UV theory via where ν = d 2 /4 − m 2 /H 2 . To simplify calculations, we can fix ν = 3/2 for any dimension, and then can regulate integrals that appear via an analytic continuation in d. We note that while this will regulate the one-loop divergences that appear in this section, dim reg alone is not sufficient in general, as we will see below.
We will begin with the one-loop power spectrum of the growing mode, illustrated in the left side of Fig. 2. In the UV theory, a standard in-in calculations gives us Regulating the IR by substituting p 2 → p 2 + K 2 in the denominator, we get This is the complete one-loop power spectrum of the UV theory.
In SdSET, the one-loop power spectrum is given by several terms where c 4,0 ∼ λ, and ϕ is given in Eq. (4.7). The subscript (n) labels the order in λ in which correlator is calculated and δα (m) is the contribution to the correlator from a shift in the value of α at order λ m . The first term is the one-loop power spectrum in the EFT, . (4.20) We must also allow for the possibility that, due to matching, we will need to correct the value of α from its free-field value. Substituting α → α+δα (1) in Eq. (3.18a) and expanding to linear order, we have We emphasize that since the coefficient of ϕ 3 + in the definition of ϕ is fixed by matching the superhorizon six-point function, there is no additional freedom within this term that can be used to match the 1-loop power spectrum.
The final term δ ϕ k ϕ k IC is time independent and is determined by matching the time-independent part of the UV calculation. Combining these results and using c 4,0 → λ, we find

(4.23)
Comparing this result to Eq. (4.18), we see that the time-dependence of the UV and EFT agree after matching the single log coefficient with This is a non-trivial result as we only have a single free parameter to match the full timedependence on the UV result. All the time-independent contributions can be matched by adjusting the initial conditions, IC (1) , which in this case is just a renormalization of the amplitude of the two point function, namely a shift in C α .

One-loop Trispectrum
Now we move on to the one-loop trispectrum, as illustrated in the right side of Fig. 2. For our purposes here, it suffices to simply match the UV divergent terms which result in time dependent log[aH] factors. 6 This will result in a correction to the c 3,1 Wilson coefficient in the EFT. While the initial conditions also receive corrections from matching the trispectrum, these would contribute to Stochastic Inflation beyond NNLO, and so we will not compute them here. For convenience, we break the UV calculation into two terms where H int is given in Eq. (3.5). In what follows, we will show that I 4 is UV finite, while K 4 is UV divergent. Then we will match this divergent contribution between the UV theory and the EFT to derive the one-loop correction to c 3,1 . We begin by evaluating I 4 . Expanding in the limit p k i , we have (4.27) 6 We remind the reader that the "UV divergences" that we are isolating here take their origins from the IR divergences of the full theory that are resummed by Stochastic Inflation. The expansions being done to simplify the integrals that appear on the UV theory side of this matching calculation amount to decomposing the correlator via the method of regions [58,59].
The integrals over τ 1 and τ 2 are independent and we can evaluate them as usual. Expanding in k i τ i 1, we have (4.28) Putting this together yields (4.29) This integral will not lead to any UV divergences. Furthermore, this result is exact in pτ since we have only expanded in k/p. Therefore all the subleading terms will be more UV convergent, and correspond to higher power corrections in the EFT. This shows that all we need to include when matching c 3,1 is the correction due to K 4 . The second term, K 4 , gives rise to more interesting UV behavior. Again expanding for large loop momentum p k i , we have (4.30) We have expressed the commutator acting on the loop momenta in φ (3) in terms of the imaginary part to simplify the calculation. In this form, the time integrals can be evaluated exactly, giving We can see the appearance of higher powers of log from the perturbative EFT contribution to the one-loop trispectrum where the compensating dimensionful factors that appear inside the log[aH] term come from expanding the prefactor in the small α limit. One can check that this term matches the coefficient of the log 3 divergence of the UV calculation. In order to match the linear log term, we need to keep track of the EFT field redefinition to order λ 2 . Specifically, we need (4.33) to remove the ϕ 6 + operator from the EFT potential. This O λ 2 term contributes to the trispectrum at one loop: which matches the leading UV term in Eq. (4.31), namely the log term proportional to a factor of 10/81. After matching this term, we see a fairly complicated expression remains for the linear log. This can be absorbed into the c 3,1 Wilson coefficient in the potential using Eq. (3.23), such that This is the one-loop matching correction to an EFT Wilson coefficient, in analogy with Eq. (4.24) above. A priori, one might think that we have to consider the divergence from the mixing of the operators ϕ + and ϕ 3 + . In later sections, we will see that such a mixing is equivalent to a shift in the potential of the form in Eq. (4.35), and therefore these two interpretations of the logarithmic growth are related to each other by a field redefinition.

Initial Conditions for Composite Operators
The above matching procedure is sufficient to regulate the correlation function of φ and match ϕ + correlators at separated points. Composite operators are defined when some of these operators are at coincident points. In Fourier space, this involves a convolution integral which can produce divergences that require renormalizing the composite operator itself.
Composite operators can be defined in this way purely within the EFT. Since we are regulating loops in the EFT with dynamical dim reg, this implies that we will need to know initial conditions for general α. Fortunately, we will be interested in limits where the momenta are hierarchical (UV divergences of the momentum integrals), which simplifies the matching considerably.
We will start by taking the tree-level initial conditions for the trispectrum, and tying two of the legs together to form a ϕ 2 + [ x = 0] composite operator: where we are interested in determining the integrand of the right hand side in the limit p k i . For α = 0, we determined the initial conditions exactly in Eq. (4.13). In order to regulate the UV divergence that comes from tying the two legs together, we want to evaluate this for general α, while isolating the term of interest, which is proportional to P + (k 1 )P + (k 2 ), where P + is defined by (4.37) so that P + (k) tree = (2k 3 ) −1 . Since k p, the initial conditions will arise at the horizon crossing of the modes carrying momentum p when the k i -modes are superhorizon. As discussed in Sec. 4.2, we therefore match using massive mode functions for only the k i fields, where we can also Taylor expand in kτ 1. As a result, the initial conditions are given by where [aH] −2α = (−τ ) 2α . The RHS of this expression is the calculation in the UV theory with the appropriate choice of masses for the mode functions. We are implicitly evaluating the correlation function at a time τ 0 and extracting the τ 0 -independent piece in the τ 0 → 0 limit. Evaluating the integral, we find . (4.39) Corrections to this result are suppressed by powers of k i /p which will not contribute to the one-loop divergences that we will use to determine the RG for composite operator mixing below in Sec. 5.1. For two-loop divergences, one must determine the initial conditions (see Sec. 5.2 below): where we will take the limit p i k, such that − k − p 1 − p 2 − p 1 − p 2 . We calculate this contribution for general α i by matching to the full theory: is the positive frequency mode for a field with α = 3 2 − ν. Again, since the mode functions at horizon crossing behave as if they are effectively massless, this contribution can be determine using massless mode function of pτ 1 , but we must use the massive mode functions for pτ . The result of integrating over τ 1 is  43) and p t = p 1 + p 2 + p 3 . Finally, we must also determined the tree-level six-point initial conditions, illustrated Figure 3: Tree level pentaspectrum in Fig. 3. We are specifically interested in the correlator As in the case of the trispectrum, we will isolate the term such that (4.46) By direct calculation we find that Γ 2,4 (p) = 1 216p 3+4α 16 + 4γ E (−11 + 3γ E ) + 3π 2 + 4(−11 + 6γ E + 3 log 2) log 2 . (4.47)

Composite Operator Mixing
From our above discussion, we argued that corrections to Stochastic Inflation are uniquely determined by the correlation functions of composite operators; computing those that are relevant to correcting Stochastic Inflation up to NNLO is the topic of this section. We will start by setting up the problem of calculating composite operator renormalization. We are interested in operators of the form ϕ n + ( x). These are well-defined purely within the EFT, so we must be able to discuss their correlation functions and renormalizations given only the EFT data. Of course, the crucial information is the initial conditions, which therefore depend on first matching to the UV.
In general, given the EFT field operator ϕ + ( x, t), we can always define a composite operator In a free theory, we find the structure: where P + is defined in Eq. (4.37). Taking the Fourier transform, For simplicity, it will often be easiest to use ϕ n Given this definition, we see that we should be able to define all such correlators of composite operators as simply integrals over the correlators of ϕ + . For example, (5.5) A very important feature of this formula is that the integrand on the RHS is free of divergences, in the sense that we should have already regulated and renormalized the expression. As a result, all of the renormalization of the composite operator itself (and hence the matrix of anomalous dimensions) has to be associated with integrals over p i as opposed to the loop integrals that appear in the calculation of ϕ + p 1 ... ϕ + p n ϕ + k 1 ... ϕ + k m itself.

One-loop Corrections Trispectrum
We will start by computing b 1 , which is determined from the four-point function via time Figure 4: One loop correction to ϕ 2 + that looks like an anomalous dimension (ϕ 2 + → ϕ 2 + ). We start from the tree level trispectrum and integrate over two of the fields to form the composite operator ϕ 2 + . Left: The Witten diagram with a boundary at future infinity. Right: The Feynman diagram with the same momentum flow.
The relationship between this loop contribution and the tree-level trispectrum is illustrated in Fig. 4. The contribution to the four-point function from the time evolution already yields a log[aH]. Therefore, any anomalous scaling (which should also only be a single log) must arise from the initial conditions. The relevant contribution was calculated in Eq. (4.39).
Performing the integration over p using dynamical dim reg, we find Expanding for α → 0 we have We can repeat this calculation with ϕ n + [0] to determine the one-loop anomalous dimension for all n. Keeping track of combinatorics, one finds This NLO log[aH] dependence can be resummed by including the following correction in the dynamical RG: (5.10)

Six-point
Our next task is to compute b 2 , which we determine by evaluating the six point function, with p k i . The relationship between this one-loop contribution and the tree-level sixpoint function is illustrated in Fig. 5. As before, only the initial conditions arising from will contribute to the operator mixing, where Γ 2,4 (p) was calculated in Eq. (4.47). Performing the momentum integration using dynamical dim reg, we find Now we restore the factor (k 1 k 2 k 3 k 4 ) α [aH] −4α associated with P (k i ) for general α and take the limit α → 0 to get . (5.15) Repeating this calculation for ϕ n + we have .
This NLO log[aH] dependence can be resummed by including the following correction in the dynamical RG: time Figure 5: Diagram of contribution one loop contribution to Γ 2,4 (ϕ 2 + → ϕ 4 + ). We start from the tree level pentaspectrum (6 points function) and integrate over two of the fields to form the composite operator ϕ 2 + . Left: The Witten diagram with a boundary at future infinity. Right: The Feynman diagram with the same momentum flow.

Two-loop Corrections
Next, we move to the calculation of the two-loop anomalous dimension that generates the NNLO non-Gaussian noise term for Stochastic Inflation. This represents a novel contribution which is calculated here for the first time. In particular, our goal is to calculate: where p 1,2,3 k. The relationship between this two-loop contribution and the tree-level trispectrum is illustrated in Fig. 6. The correlation function on the RHS was calculated in Eq. (4.42) for general α i (for each p i ) such that we can regulate the integral with dynamical dim reg. Making the above substitution where κ is fixed by the full calculation for the reasons described above, see Eq. (4.43). We note, that this result will require us to calculate several integrals of the form Figure 6: Diagram of contribution two loop contribution to the mixing ϕ 3 + → ϕ + . We start from the tree level trispectrum and integrate over three of the fields to form the composite operator ϕ 3 + . We first drew this as a Witten diagram with a boundary at future infinity. At the bottom we have shown the Feynman diagram with the same momentum flow.
where a, b and c are half integers when α = 0. We can evaluate this integral as follows: where we introduced an IR regulator K as we did for our 1-loop divergence. This IR regulator is only needed when a + b + c 3 (the integral vanishes by dynamical dim reg otherwise). One should also notice that enforcing a + b + c → 3 restores the invariance under permutations of a, b and c.
We always have at least one log from Γ[a + b + c − 3]. Therefore, the only contributions that are not log 2 (or higher) are those where all the other Γ functions are finite. Up to permutations, there are three relevant cases (i) a = 3/2, b = 3/2, c = 0, (ii) a = 3/2, b = 1, c = 1/2, and (iii) a = b = c = 1. Only (iii) is a single log as expected from above. Isolating just the single log term, for example using a = 1 − α/2 with b = c = 1, we get Repeating this calculation for ϕ n + we have . (5.23) This NNLO log[aH] dependence can be resummed by including the following correction in the dynamical RG: (5.24)

Stochastic Inflation at NNLO
Now we have computed all the necessary pieces. To summarize our results, we have found that at NNLO, the Fokker-Planck equation for Stochastic Inflation becomes The rest of this section is devoted to expressing this result in a particularly simple basis. As we have emphasized above, the coefficient b 1 and b 2 are basis dependent, in the sense that they can be changed by taking a field redefinition. For the purposes of solving for the NNLO equilibrium distribution, we will find it useful to first perform a field redefinition that moves their effects into the potential. For contrast, we note that the coefficient d 0 is basis independent, and furthermore field redefinitions of this (∂/∂ϕ + ) 3 term will only induce higher order (NNNLO) terms that we will neglect.
Concretely, we want to redefine the fieldφ + = f (ϕ + ) such thatb 1 =b 2 = 0. Using P = (dφ + /dϕ + )P under such a field redefinition, we see that Next, in order to setb 1 =b 2 = 0, we defineφ + so that We can integrate this equation to determinẽ The remaining term is then determined to be d 2φ This basis change also impacts the effective potential that appears in the Fokker-Planck equation: The first term, 3 2 b 1φ+ , simply provides an O(λ) correction to α, i.e., the quadratic term in the potential c 1,1 can always be removed by redefining α, just as we did above when matching the 1-loop matching power spectrum, see Eq. (4.24). The new contribution to the second term, b 2 = O(λ 2 ), is simply a correction to the definition c 3,1 , again just as was computed above when matching to the 1-loop trispectrum, see Eq. (4.35). The same is true for higher powers of ϕ + that we have dropped. In short, we see that b 1 and b 2 simply shift the definition of the couplings within SdSET at higher order in λ but do not introduce any new terms. As a result, we may simply define λ eff = λ + 18b 2 + 3! δc 3,1 (5.32a) where δc 3,1 is the correction from one-loop matching in Eq. (4.35). For the other contributions to Stochastic Inflation, λ eff = λ is sufficient to achieve NNLO accuracy. In this sense, the coefficients of the NLO and NNLO corrections to the potential are independent of redefinitions of λ. This same argument explains the scheme-independence of β-functions up to two loops. Putting this all together and relabelingφ + → ϕ + , we arrive at a canonical form for the NNLO equation: Presumably this freedom to put these equations into such a canonical form can be recast in terms of a covariant description in field space [60,61], and multi-field generalizations thereof.

Implications
Having derived novel corrections to the equations that govern the Markovian evolution of the probability distribution for ϕ + , we now turn to solving them. In particular, this section will explore the physical implications of these new terms by calculating the NNLO equilibrium distribution for ϕ + , assuming a static dS background. Then we will extend this to account for the dynamics of the system as it relaxes to the equilibrium state. To this end, we will set up the formalism to calculate the "relaxation eigenvalues," and will numerically solve for these quantities to O(λ 3/2 ).

Equilibrium Probability Distribution
Starting with the canonical form of the NNLO equation governing Stochastic Inflation, we can understand the impact on the equilibrium probability distribution P eq (ϕ + ), which by definition satisfies ∂ ∂t P eq (ϕ + ) = 0. We can solve the problem non-perturbatively in the absence of the higher derivative term proportional to d 0 . Therefore, our strategy will be to find the solution in terms of general V eff , and then include the correction from d 0 as a perturbation.
The equilibrium contribution from V eff can be determined non-perturbatively following Sec. 2.1. The equilibrium solution satisfies and P V eq is defined as being a solution to this equation. 7 Integrating twice gives the solution Note that V eff is only a function of ϕ + and should not be confused with V (ϕ + , ϕ − ) in Eq. (4.11). Since this solution holds for any V eff , it gives the answer at both LO and NLO provided we include the NLO contributions to V eff . At NNLO, in addition to the correction to V eff , we must include d 0 = λ eff /(192π 2 ) which alters the equation for the equilibrium solution We can integrate this equation once to get This can be solved using separation of variables, so we will make the ansatz P eq (ϕ + ) = P V eq (ϕ + )Q(ϕ + ) , (6.6) 7 Note that P V eq (ϕ + ) includes any higher order correction in V eff but not does not including higher derivative terms that arise at NNLO and beyond. In this sense P V eq is not to be confused with the LO solution, as it contains some (but not all) contributions at every order.
where P V eq is given in Eq. (6.2), including the NNLO corrections to V eff . Plugging this ansatz into Eq. (6.5) gives We can solve this equation perturbatively in Q. The zeroth order term is Q = constant. At the next order, clearly Q is order d 0 so we can neglect derivatives of Q in the first terms such that such that the solution becomes Finally, we use d 0 = λ eff /(192π 2 ) and V eff λ eff ϕ 4 + /4!, which consistently captures effects up to NNLO accuracy. We then evaluate the integrals and simplify the expression to find Using the fact that λ eff ϕ 4 + = O(1) for the LO equilibrium solution, we see that both terms in Q are O(λ eff ), as expected for NNLO accuracy. Recall the NLO and one additional contribution at NNLO are encoded in V eff (φ + ) and are included in P V eq (φ). Combining these terms and writing P eq = CP LO (ϕ + )P NLO (ϕ + )P NNLO (ϕ + ), we have

Relaxation Eigenvalues
In this section, we will explore the implications for the time dependence of P (ϕ + , t). This can be characterized by computing the so-called "relaxation eigenvalues" as we explain below. In the previous section, we could calculate the analytic NNLO equilibrium probability distribution where we only had to treat the (∂/∂ϕ + ) 3 term perturbatively. Here, we must resort to a numerical evaluation of the perturbative expansion, which requires that we treat all higher order corrections as perturbations.
To begin, we return to the full NNLO equation that governs Stochastic Inflation given in Eq. (5.33a), and rewrite it as a Euclidean Schrödinger equation [23,31]: This equation can be solved using separation of variables where V eff (ϕ + ) is defined in Eq. (6.3), we have assumed t 0 = 0, and Φ n are the eigenfunctions of [23,31] where Λ n are the non-negative "relaxation eigenvalues," and the Schrödinger potential is The lowest eigenvalue is zero with the eigenfunction Φ 0 (ϕ + ) ∝ exp − 4π 2 3 V eff (ϕ + ) . Since all other Λ n are positive, at late times the distribution P (ϕ + , t) relaxes to the fixed point This reproduces the static result above in Eq. (6.2). However, for the numerical evaluations we take V eff (ϕ + ) = λ eff 3! ϕ 3 + , and treat the additional correction to the potential as perturbations. In addition, for the remainder of this section, we will drop subscript 'eff' on the coupling, λ eff → λ, for brevity.
To explore the effect of the NNLO correction we need to determine the eigenvalues Λ n of Eq. (5.33a) for n ≥ 1. We can find the eigenvalues and eigenfunctions of Eq. (5.33a) as perturbative expansions in powers of √ λ, n are the solutions of Eq. (6.14) with Schrödinger potential U = U (0) . This potential is obtained from Eq. (6.15) by setting V eff (ϕ + ) = λ 3! ϕ 3 + , the leading term in Eq. (5.33b). Explicitly, Since ϕ + ∼ λ −1/4 we see that the Schrödinger potential, the term in square brackets, indeed scales as O λ 1/2 . This equation can be solved numerically by the shooting method [62]. Higher order terms in Eq. (5.33b) as well as the NNLO term ∂ 3 ϕ + (ϕ + P ) can be treated as perturbative corrections to Eq. (6.18). That is, if we substitute Eq. (6.13) into Eq. (5.33a) and keep terms up to O λ 3/2 , the eigenfunctions Φ n must solve Finally, the perturbative corrections to the eigenvalues are computed numerically using The first few relaxation eigenvalues are (recall that λ = λ eff here) The contribution to the eigenvalues at O λ 3/2 includes both corrections to the equations of Stochastic Inflation at NNLO as well as perturbative corrections to the eigenvalues from the LO and NLO equations. Working with the full NNLO equations was crucial to obtaining the detailed numeric values. We note that the NNLO contribution to V eff in Eq. (5.33b) dominates; the small numerical coefficient of d 0 9 × 10 −5 λ suppresses its impact on these eigenvalues. As the form of V eff is determined by SdSET field redefinitions to all orders, it may prove useful in future studies to simply include higher order corrections to the potential as an approximation. We expect this minor impact of d 0 is due to the fact that we are assuming the UV theory is λφ 4 such that the non-Gaussian noise and corrections to the potential are determined by the same parameter. In contrast, if one were to consider inflationary models with primordial non-Gaussianity, these two effects are controlled by independent parameters, such that the non-Gaussian contribution to the noise could become important [4]. In either case, this investigation is only possible because we have framework in which all corrections to Stochastic Inflation, including contributions from non-Gaussian noise, can be systematically computed.

Conclusions
Understanding the nature of quantum dS space is one of the most basic conceptual problems in cosmology [43]. Stochastic Inflation [21][22][23] informs much of the physical intuition for how we think about accelerating cosmologies, particularly as we approach the eternally inflating regime that is dominated by quantum fluctuations [36][37][38][39]. Yet, Stochastic Inflation is itself an approximation whose regime of validity, and corrections thereof, should follow from a more basic starting point. Ultimately, a complete description should include dynamical gravity, although the simpler case of quantum field theory in a fixed dS background studied here already provides a non-trivial challenge.
In this paper, we demonstrated precisely how corrections to Stochastic Inflation arise from quantum field theory in dS, namely as a natural consequence of dynamical renormalization group flow within the EFT that emerges in the superhorizon limit. By working with SdSET, the origin of the stochastic description is a direct consequence of EFT power counting, which also explains why this effect is only relevant for light (massless) scalars. (This same power counting scheme also explains the all orders conservation of the adiabatic mode [33].) By matching λφ 4 theory onto SdSET up to one loop, we could then calculate the log enhanced corrections to the mixing of EFT operators up to two loops. This allowed us to derive the corrections to the equations of Stochastic Inflation at NNLO accuracy. These results include the first higher derivative correction to the framework, which is the leading signature of the non-Gaussian contribution to the noise as modes cross the horizon.
This work extends derivations of Stochastic Inflation from quantum field theory in dS at LO [29,31,33,35] and NLO [30,32,34] to NNLO. Yet, even at NLO, we showed that the "universal" correction to the effective potential follows from a field redefinition and can be extended to all orders. This result agrees with Refs. [30,32,34], which arrive at this NLO correction by (effectively) integrating-out the decaying mode at tree-level. Furthermore, the first appearance of non-Gaussianity in the stochastic noise appears at NNLO and requires a genuine two-loop calculation. Higher-loop calculations of inflationary correlators are notoriously difficult, but they are made manageable by working with SdSET, which facilitates the use of the symmetry preserving dynamical dimensional regularization. Most importantly, SdSET reduces the problem of calculating any corrections to Stochastic Inflation to the determination of the matrix of anomalous dimensions. Rather than being a mysterious feature of dS space, we now see that the derivation of Stochastic Inflation and corrections thereof is conceptually and technically similar to calculating the scaling dimensions of operators at the Wilson-Fisher fixed point in d = 4 − dimensions. Finally, there is an intriguing connection that could be made with a thermodynamic interpretation of the equilibrium probability distributions, P LO ∼ exp(−βE), where the inverse temperature is β = 2π/H. It would be interesting to understand the meaning of the NNLO corrections derived here from this point of view.
Phenomenologically, Stochastic Inflation is an important tool for understanding the predicted non-Gaussianity in multi-field inflation [44,[63][64][65][66][67], where superhorizon evolution can give rise to non-trivial correlations. Previous work has included the non-Gaussian contributions for the non-linear superhorizon evolution, but thus far the effects of non-Gaussian noise has been missing. It will be interesting to explore models where both effects are simultaneously important. For example, one might hope this techniques would elucidate the physics of the small mass regime of quasi-single field inflation [68], which is known to produce large logs.
Conceptually, Stochastic Inflation serves as the basis for much of our understanding of slow-roll eternal inflation. This description requires coupling a light scalar field to gravity, yet much of the structure is determined by the quantum noise in the Fokker-Planck equation. Specifically, the regime of slow roll eternal inflation is the limit where the potential becomes flat and the quantum noise dominates the time evolution until inflation ends. While understanding this regime is often considered a conceptual problem, it may have important consequences for cosmological solutions to hierarchy problems, such as [69][70][71][72].
Finally, underlying our results on Stochastic Inflation is a demonstration that SdSET is a consistent description of dS quantum field theory at loop level. Calculations in a wide variety of cosmological settings are beset with challenges stemming from the underlying time evolution and lack of consistent regulator. The successful implementation of SdSET as an organizing principle for calculating quantum correlators in dS offers hope that more of these cosmological problems may be organized and simplified when described with the right degrees of freedom. The emergence of Stochastic Inflation as a simple consequence of EFT power counting is a non-trivial example of these principles in action. and Our goal is to match this expression onto the EFT, so we need to take the limit where all of the fields are superhorizon. The additional contributions that arise from the subhorizon region, k i τ j = O(1), will be absorbed into the initial conditions, which we do not need to calculate explicitly for our purposes in this work.
To match the superhorizon behavior, we can again expand the operator using Eq. It is easy to see that the superhorizon contribution to A 6 is determined by [φ] λ , This shows that by matching the trispectrum with Eq. (4.8), we can also match A 6 . Clearly we cannot determine B 6 using [φ] λ . The most straightforward way to determine the superhorizon contribution is to write where G k, t , t = i φ k, t , φ k , t Integrating this expression using dynamical dim reg, we find the superhorizon contribution is Next, we would like to see how this formula arises in SdSET. First, we calculate the contribution to the six-point function at second order in c 3,1 . We take the same commutator structure as B 6 , where one [ϕ + , ϕ − ] acts on the external line and the other on an internal commutator. The result is Combining these terms and using c 3,1 = c 4,0 = λ, we match the UV six-point function Note that there is some ambiguity in the constant term due to scheme dependence associated with regulating our (divergent) time integrals. Although our expression matches the constant as well, in some other schemes, the initial conditions may play a role in matching.
On the other hand, all powers of log aH must match in any scheme, as we find here.

B Hard Cutoff Calculations
In the main text, we used dynamical dim reg for the EFT loop calculations. Loops in the UV calculations were, in some cases, regulated with dim reg rather than dynamical dim reg. These regulators offer some technical advantages but one might worry about using different regulators in the matching calculation. We can therefore gain further conceptual insight by redoing these calculations with a hard cutoff. This regulator can be easily implemented in both the UV and EFT and also makes the origin of divergences more transparent. Furthermore, we can work directly with the massless mode functions, thereby avoiding the complications of working with massive modes. In this appendix we will repeat the calculations from the main text using a hard cutoff, reproducing all the above results up to differences in scheme dependent coefficients.

B.1 Matching
In this section, we compute the matching for λφ 4 onto the SdSET up to one-loop order. We will use a hard cutoff for both momentum and time integrals. Specifically, we regulate the momentum integral with a UV cutoff Λ = [aH] and an IR cutoff K. For time integrals, noting that the UV region of integration does not contribute due to our i prescription, we simply regulate the IR with a cutoff t , which corresponds to the time of horizon crossing for a mode k.
Combing these results we have This expression differs from our result with dynamical dim reg, Eq. (4.24). This is not entirely surprising as the precise definitions of the parameters in the UV are scheme dependent and thus this scheme dependence is also inherited through matching.

Trispectrum
Next, we will perform the calculation to match the trispectrum taking α = 0 on all legs and regulating all integrals with a hard cutoff. As was argued above, K 4 is the only term that can generate log divergences from the loop momentum integrals. Taking p k i we have  As we did in the main text, we are focused on the single log[aH]/K term because it cannot be absorbed into the initial conditions and the RG implies that higher powers of log should be products of logs already present in lower order diagrams. In order to compare this to the EFT, we need to keep track of our field redefinition to order λ 2 . Specifically, we need The sum over permutations is essential in this calculation as the split into p 1 , x 2 and x 3 breaks the manifest permutation invariance of the measure of integration, which is only valid if the integrand itself is permutation invariant. This calculation reproduces the result from the main text where this contribution is log 2 , although the need for two separate regulators makes this less transparent.