Stochastic Inflation in General Relativity

We provide a formulation of Stochastic Inflation in full general relativity that goes beyond the slow-roll and separate universe approximations. We show how gauge invariant Langevin source terms can be obtained for the complete set of Einstein equations in their ADM formulation by providing a recipe for coarse-graining the spacetime in any small gauge. These stochastic source terms are defined in terms of the only dynamical scalar degree of freedom in single-field inflation and all depend simply on the first two time derivatives of the coarse-graining window function, on the gauge-invariant mode functions that satisfy the Mukhanov-Sasaki evolution equation, and on the slow-roll parameters. It is shown that this reasoning can also be applied to include gravitons as stochastic sources, thus enabling the study of all relevant degrees of freedom of general relativity for inflation. We validate the efficacy of these Langevin dynamics directly using an example in uniform field gauge, obtaining the stochastic e-fold number in the long wavelength limit without the need for a first-passage-time analysis. As well as investigating the most commonly used gauges in cosmological perturbation theory, we also derive stochastic source terms for the coarse-grained BSSN formulation of Einstein's equations, which enables a well-posed implementation for 3+1 numerical relativity simulations.


I. INTRODUCTION
Inflation theory was postulated more than 40 years ago as an explanation for the apparently fine-tuned initial conditions of the Hot Big Bang [1][2][3].The proposal gained traction as it also offers a natural mechanism for generating the initial density inhomogeneities [4][5][6][7][8][9] which in later stages of cosmic history led to the formation of cosmic structure via gravitational instability.These density fluctuations are directly observable in the Cosmic microwave background and their two-point statistics have been measured to very high precision [10].
Inflation also predicts the existence of a cosmological gravitational wave background [11], yet to be detected, as well as the possible existence of non-zero higher-order spatial correlation functions in the cosmological fluid [12].The latter's amplitude, and the amount of nongaussianity more generally, depend more heavily on the specifics of the inflationary model.They could be detectable in the cosmic microwave background (CMB) and cosmological structures and even lead to the formation of primordial black holes.If detected, any such features would open up remarkable windows into the early universe.
Central to the inflationary origin scenario is the assumption that our universe originated from quantum processes.This underscores the necessity to combine quantum mechanics and gravity to make precise predictions for inflationary models of the early universe.Although a theory of quantum gravity is lacking, despite decades of efforts [13][14][15][16][17][18], there are regimes where predictions can still be made by using techniques from Quantum Field Theory on Curved Spacetime (QFTCS) [15] or by constructing Effective Field Theories (EFTs), which have made continuous advancements in cosmology, inspired by the latter's success in flat space [19].Abandoning pretenses of completeness, an EFT establishes a region of validity, normally bounded by ultraviolet (UV) and/or infrared (IR) cutoffs, and the narrative of theoretical physics is implicitly about pushing these cutoffs to their limits.
Not long after the original concept of inflation was introduced, Stochastic Inflation (SI) was formulated in [20], was used as a non-perturbative methodology for light scalars in de Sitter in [21], and was later adopted as a special kind of EFT for inflationary physics, as part of a semiclassical stochastic gravity approach [22].Stochastic Inflation therefore, can be thought of as an EFT for light scalar fields in quasi-De Sitter spacetimes that sets its UV cutoff at the point where quantum operators governing the field perturbations can be described as classical stochastic variables and super-Hubble correlators can automatically exhibit decoherence, a distinctive feature of such systems.In most inflation models, this transition occurs when Fourier modes cross the Hubble radius and their transition into the EFT's range of validity is modelled by the action of a small stochastic noise perturbing the long wavelength fields described by General Relativity -in a sense, Stochastic Inflation is a combination of deterministic evolution with a continuous reset of initial conditions by small stochastic amounts coming from the influx of new modes which cross the UV threshold.The stochastic equations can then be used to derive non-perturbative results such as resummations of classical loops which are IR divergent in QFTCS [23][24][25].
In this work, as in most others on the topic, the sub-Hubble scales will be described by linear cosmological perturbation theory (CPT); although higher order perturbation theory can in principle be used for these sub-Hubble scales [26], we will not undertake such an endeavour here.Stochastic Inflation therefore, in its current formulation, sacrifices non-linearities arising in the sub-Hubble regime in favor of a fully non-perturbative theory in the super-Hubble regime, where QFTCS might fail as it has already been reported [24].It is therefore particularly suited for circumstances where non-linear evolution on super-Hubble scales, and possibly non-perturbative correlators or even full probability density functions (PDFs), are the relevant quantities to compute.
Although able to provide non-perturbative quantities, Stochastic Inflation as usually formulated comes with several approximations that simplify the complex dynamics of General Relativity, most notably relying on the Separate Universe Approximation (SUA) [27][28][29] and therefore altogether dropping the dynamics of some degrees of freedom of the gravitational field.To find the spectrum of the stochastic perturbations, it has also been claimed to require a specific choice for the time-slicing to ensure consistency with QFTCS, namely using a uniform e-fold time slicing [30][31][32][33].However, there are formulations where such restrictions are not required [34,35], which appears to be linked to an ongoing discussion about the appropriate inclusion of the momentum constraint [35,36].Furthermore, obtaining the observationally relevant curvature perturbation then requires a "first passage time" analysis (FPTA) to determine a stochastic number of efolds for a specified point of the scalar potential to be reached [37].
In this article, we endeavour to demonstrate that it is possible to generalize the stochastic inflation equations by making the linear treatment on sub-Hubble scales the sole assumption, thereby removing the need for all the aforementioned approximations.More specifically, we formulate equations for Stochastic Inflation within full General Relativity in its ADM formulation, retaining all the variables describing the gravitational field and the freedom to choose the time-slicing and the spatial coordinates within the 3D time-slices.
To achieve this, we will proceed as follows: After recalling the ADM formulation and its linearized version in Sec.II, we survey a range of commonly used gauges and gauge invariant variables in Sec.III, where we also provide explicit expressions for determining the spacetime foliation and the dynamical variables from knowledge of the gauge invariant comoving curvature perturbation R. Our procedure for coarse-graining and obtaining stochastic source terms for the full set of ADM equations is explained in Sec.IV.These stochastic sources are expressed in terms of R and are found to be identical in all the small gauges we examine.We therefore postulate that they represent the complete small gauge invariant stochastic continuation of the Einstein equations with linearized source terms, noting that they satisfy lin-ear perturbation theory by construction in any choice of gauge/foliation.Sec.IV E extends the result to include tensorial modes sources.As a short application, we present and solve stochastic equations for ∆N in Sec.V in a toy model exhibiting a short phase of ultra slow-roll, directly obtaining the PDF of the curvature perturbation without the need for a first-passage-time analysis.***** Notations: Units are such that ℏ = c = 1.The Planck mass is written M P l = (8πG) − 1 2 .Quantum operators will carry () while random variables will be bold symbols.We adopt the following Fourier transform convention:

ADM EQUATIONS AND SCALAR LINEAR PERTURBATIONS
In this section, we review the ADM formulation of the equations of General Relativity and their linearized perturbations around a homogeneous and isotropic background.For further purposes that will become clear further ahead, all equations are written as left-hand-sided.

A. ADM formulation of Einstein's equations
For an unknown 4-dimensional metric tensor with components g µν describing a spacetime with Ricci scalar R and 'matter' Lagrangian density L m , the Einstein-Hilbert action with vanishing Gibbons-Hawking-York boundary terms [38,39] is When varied with respect to the metric g µν it gives Einstein's equations where δLm δgµν such that ∇ µ T µν = 0 by construction.In this work, we will study a single ϕ scalar field fluid described by In what follows we will use the common notations for the 3+1 splitting of ADM formalism [40], letting n and P be the projectors along the normal of the 3-space submanifold and the submanifold itself respectively.The associated covariant derivative component on the 3-space submanifold will be denoted by a vertical bar () |i .
We start from the ADM parameterization of the metric, assuming the existence of a foliation of spacetime The Einstein-Hilbert action can be re-expressed to reflect this 3 + 1 splitting where 3 R the Ricci scalar of the spatial metric γ ij , and Π is the conjugate momentum and K is the extrinsic curvature which is usually split into its trace and traceless components K and Kij respectively.Variation of ( 5) with respect to ϕ yields the field evolution equation while variation with respect to γ ij yields the gravitational field's evolution equations where the 3 + 1 components of the energy-momentum tensor are the energy and momentum densities and the stress tensor S ij = P µi P νj T µν = T ij , decomposed as Finally, variation with respect to α and β i leads to two constraints usually referred to as the Hamiltonian and Momentum constraints respectively.

B. Scalar Linear equations
This section reviews the linearised version of the ADM formalism of scalar perturbation theory without making any gauge choice [41].The linear equations presented below match exactly the non-linear left-hand sides (LHS) of the previous subsection.Since the standard Einstein and ADM systems of equations are equivalent up to substitutions and factors, the equations below are also equivalent to the more commonly adopted formulations of cosmological perturbation theory.
We can incorporate general scalar perturbations through the line element where the four scalar functions Φ, Ψ, B and E define the 3 + 1 perturbations of the full ADM metric and a the usual background scale factor.The full dictionary from ADM to the 1 st order scalar decomposition follows from the identification of the metrics which leads to the linear extrinsic curvature where ∆ = ∂ i ∂ i = ∇ 2 /a 2 and the energy-momentum densities where σ = 0 (no anisotropic shear in this work, which is gauge-invariant to say), which means that Sij = 0 at first order.Note that the 3-Ricci tensor is also very simple and is given by at first order.Substituting these expressions into the full ADM equations yields a set of 0 th and 1 st order perturbation equations which do not yet incorporate a choice of gauge.
from the Hamiltonian constraint, the extrinsic curvature, and the field evolution equations respectively, while other equations are vanishing.These are of course the equations of homogeneous and isotropic FLRW cosmology.At first order, the gravitational evolution equations become: while for the field equation Finally, the first-order constraint equations are: Hence we have five equations for the five unknowns δϕ, Ψ, Φ, E and B. In general, we can use the gauge freedom to fix two of these variables.The two constraints then further reduce the dynamical scalar degrees of freedom down to one, usually encapsulated in a variable that is also gauge invariant.

III. SPACETIME FOLIATIONS AND SCALAR FIELDS FROM THE COMOVING CURVATURE PERTURBATION
In this section, we will explicitly show how a choice of gauge combined with knowledge of the gauge invariant variable R determine, to first order, the 3D spacetime hypersurfaces corresponding to the ADM foliation and all perturbation variables.We will use the relations of this section later when we derive our stochastic equations.We start by recalling a set of well-known gauge invariant variables used in the literature.
In the following and unless stated differently, α b = 1, which implies that () is a cosmological time derivative and H is the common Hubble rate.

A. Scalar perturbations, all-in-one
Gauge-invariant quantities present many advantages beyond their invariance under coordinates transformations: they naturally reduce the dynamics to the single dynamical scalar degree of freedom and their dynamical equations guarantee that the constraints of General Relativity are automatically satisfied.Many linear combinations of δϕ, Ψ, E, B and Φ are known to be gaugeinvariant at first order, in particular, Bardeen potentials Φ B and Ψ B [42], the gauge-invariant field perturbation δϕ gi , the curvature perturbation on uniform density hypersurfaces ζ gi or the curvature perturbation on comoving hypersurfaces R.They are given by Those are convenient for the community as they appear to have simple gauge-invariant evolution equations.For example, for a single field [P (X, ϕ) = X − V (ϕ)]-theory, the Fourier modes of R satisfy [4,43] Rk where ε 2 = − ε1 /Hε 1 and ε 1 = − Ḣ/H 2 are the slow-roll parameters.A remarkable feature of the above equation is the perfect cancellation of R's effective mass term.As we mentioned above, by knowing the background and a universal description of the perturbations such as R, one can then write all fields in any gauge, including ones employed in Numerical Relativity.Examples are given in the next subsections.
In this work, we will use R as our master gaugeinvariant variable.It is possible to express all other gauge invariant variables in terms of R as These relations are most easily obtained in the comoving gauge (see below) but since they are relations between gauge invariant variables they hold in all gauges.Eqn.(24) will be more than useful in some of the gauges we will examine below when the constraint equations will not be directly solvable.Note also that by using eqn.(23) one easily shows that Φ B = Ψ B which is known to hold for a scalar field as there is no anisotropic stress at linear order.We will however keep these two variables distinct as they will differ by a stochastic term when we consider coarse-graining -see eqn.(52).

B. Spacetime foliations in different gauges
1. Pure gravity -comoving gauge A gauge will be called pure gravity if it sets δϕ and E to small arbitrary spacetime functions δϕ * and E * , their smallness being of the order of other perturbative quantities such as R. The null case δϕ * = E * = 0 is commonly used in cosmology and is referred to as the comoving gauge.Such a pure gravity gauge is probably among the simplest ones to deduce the whole 3D hypersurface.Indeed, Φ is given directly via the definition of R together with the chosen δϕ * .Substituting it in the linearized energy and momentum constraints (eqn.(21) yields the result for Ψ and Alternatively, we could have used eqn.(24) directly to obtain the same relations.Either way, all quantities have been expressed in terms of R, and the constraints are satisfied by construction.
For the comoving gauge we have, in particular and both the foliation and the spatial geometry are expressed in terms of R.

Fixed spatial curvature -spatially flat gauge
A gauge will have fixed spatial curvature if it sets Φ and E to small arbitrary spacetime functions Φ * and E * .
The null case Φ * = E * = 0 corresponds to the wellknown spatially flat gauge of cosmological perturbation theory.
This family of gauges is as easy to handle as the pure gravity family above.Noting that R is proportional to one of the perturbation quantities, here δϕ, the Momentum constraint can be used to obtain Ψ while the Hamiltonian constraint gives B Again, we could have used eqn.(24) directly to obtain the same relations.We have therefore obtained the foliation and the field perturbation in terms of the only degree of freedom for the evolution of the perturbations.The hypersurface geometry is of course explicitly given by the gauge choice here.

Newtonian gauge
A gauge will be called Newtonian if it sets B and E to small arbitrary spacetime functions B * and E * .Unlike the two previous gauge families, in this type of gauge using the the constraints alone is not enough to express all quantities in terms of R. We can however resort to eqn.(24).Starting from χ * = −a 2 (B * − Ė * ), one gets Φ from Φ B and then Ψ is obtained from either the Mo- The standard case in cosmology is of course obtained for B * = E * = 0.

Uniform N
If one defines N to be the number of e-folds elapsed from an arbitrary time t 0 via then a gauge will be of uniform (or fixed ) N if δN is fixed to a small arbitrary spacetime function δN * while B = 0.As explained in sec.IV B, this type of gauge with δN * = 0 is the one mainly used in the formulations of Stochastic Inflation in the current literature.
As N ≡ − 1 3 t0 Kαdt, one gets N b = t0 Hdt at zeroth order and δN = Φ+ 1  3 k 2 E at first order [32].To calculate the whole hypersurface, one first needs to solve for E by using Φ B to obtain the solution of which provides E as with E(t 0 , ⃗ x) ≡ E 0 (⃗ x).Note that there is still some gauge freedom left in the choice of the initial value of E in space.Once E is known, one gets

Generalised synchronous gauges
A gauge will be called generalised synchronous if it sets Ψ and B to small arbitrary spacetime functions Ψ * and B * to stay within perturbation theory.The null case corresponds to the synchronous gauge in cosmology.Similarly to the Newtonian case, in this gauge family constraints are not enough to solve in terms of R, which is why we use eqn.(24) on top of it.This inconvenience is an indication of the gap between numerical relativity and cosmology gauges.
Using the definitions of the gauge-invariant Ψ B and Φ B and the latter's relation to R, see (24), we can obtain the following cascade which are all functions of R and its 1 st and 2 nd order time derivatives and their integrals.
In this gauge both χ and E need to be initialized at a given time by space-only functions χ 0 (⃗ x) and E 0 (⃗ x).It is indeed a well-known problem that the synchronous gauge (and so more generally any small extended synchronous gauge) does not fix all gauge degrees of freedom and that there is still spatial dependence [44], even though this can sometimes be confused with physical choices in the literature, see Appendix B for proof.In this work, we will choose gauges such that χ 0 = E 0 = 0 at an arbitrary time, hidden in an implicit boundary of the time integrals.

IV. STOCHASTIC EQUATIONS FOR GENERAL RELATIVITY A. Schematic formulation of effective stochastic IR dynamics
As mentioned in the introduction, Stochastic Inflation can be thought of as an Effective Field Theory (EFT) [45] of QFT on curved spacetime (QFTCS), valid on scales where quantum correlation functions can be well approximated by classical, stochastic ones.Its utility lies in treating fluctuations beyond perturbation theory and associated truncations in orders of non-linearity on super-Hubble scales, compared to QFTCS [46].
By definition, the EFT refers to variables that are coarse-grained beyond a certain length scale, normally commensurate with the Hubble radius, but there is no general method for this coarse-graining in inflationary spacetimes.The greatest difficulty of applying known EFT techniques with SI is probably the nature of its IR-UV split: the cutoff is spacetime-dependent.While building an EFT from a path integral [47,48] is a possible approach in simple cases of test fields in de Sitter (dS), no such approach that also includes the fluctuations of the metric in inflationary spacetimes has been achieved.
In this work, we use a coarse-graining of the equations of motion (EOM), which is probably the most common approach for reasons that will become clearer later.Given the complexity of the complete set of equations, we first schematically review its philosophy.Let's assume that we possess an IR classical theory giving access to some second-order tensorial partial differential equations and their linearization at first order in CPT of the form where X i , Ω i (X), are tensors of identical rank with Ω i being functions of the X i .δX i and δΩ i are their linearized perturbations, δΩ i being a linear combination of δX i and their partial derivatives.Those equations could be the previous section's or approximated versions of them or any others.No assumption is made on Λ, Ω, λ, and δΩ except their smoothness and that they depend on {X j , ∇X j } j , and spacetime coordinates.No straightforward link between λ and Λ or δΩ and Ω can be written down, although one has to bear in mind that δΩ i can carry perturbations of the {X j } j and of g along with their first partial derivatives.Note that the background equation has been kept implicit but answers the relation The background is most of the time assumed homogeneous, which will be the case in this work.
The next step is to find a solution of the first-order equations in Fourier space for {δX j ⃗ k } j , in at least one gauge.In the case of Inflation, the spatially flat gauge is the most common one to use where the dynamical field is equal to the curvature perturbation on comoving hypersurfaces following the Mukhanov-Sasaki equation, eqn.(23).After that, a gauge transform can bring the calculated spectrum to the desired gauge, although those extra terms are usually neglected under the right assumptions [32].The {δX j } j are assumed to remain dynamical quantities in this new gauge, i.e. with a non-zero conjugate momentum in the sense of Field Theory (FT).This assumption is necessary to make sure the equations are dynamical but also to make sure that canonical quantization is possible.At this stage, quantization is usually written down explicitly as an expansion in the modes where ) and any other commutator being 0. This thus gives us a Gaussian spectrum.
With the solutions of the linearized equations at hand, the principle of the literature is simple and is called IR-UV splitting [20].The quantised perturbations are indeed split as δ Xi = δ Xi> + δ Xi< , i.e. into long and short wavelength parts respectively1 , for instance with a spacetime-dependent window function in Fourier space For SI, the window function would typically let short modes become long ones when their wavelegth equals or exceeds the Hubble horizon, namely k ≲ aH, to ensure classicalisation.
The UV-IR splitting via some window function implies that the linearized equation for the long wavelength δX i> , simplified by the use of X i b 's and δX i 's equations, will have a non-zero right-hand side (RHS) term (compare with eq. ( 34)) where the RHS Σi is a Fourier expansion of functions of δ Xi> and its first derivatives, but also of the window function derivatives. 2Eq. ( 37) is a quantum Langevin equation [49] obtained after coarse-graining and quantisation.To get a classical equation, i.e. to get rid of the operators, we need the Stochastic approximation to write where α ⃗ k are now random variables such that To match the Gaussian spectrum of the linearized quantum operators, one should take and Arg[α ⃗ k ] ∼ U(0, 2π).The validity of this classicalisation can be tested by showing that correlators of δ X> and associated conjugate momenta receive negligible contributions from non-commutative terms on long wavelengths.For SI, the validity has already been studied and the following analysis will be restricted to the applicable cases [20,[50][51][52].
Finally, under the assumption that linear theory is enough at horizon-crossing for computing the stochastic Σ i , the LHS of (37) can be promoted to the full nonlinear equation in (34) with the understanding that it refers to the coarse-grained quantity X i> : The above equation constitutes the long wavelength, stochastic version of the original and is postulated to furnish an adequate approximation to the long wavelength fluctuation dynamics.The line of reasoning leading to (40), applied to a truncated subset of the Einstein equations, underlies most existing expositions of stochastic inflation equations in the literature.This linear source approximation, together with its assumed Gaussianity, is a common feature of the stochastic inflation literature [20,47,53], but its justification is rarely made explicit.The validity domain of this approximation lies in scenarios where, before horizoncrossing, UV scales have suppressed higher-order statistics compared to the UV tree-level or compared to their IR counterparts generated from the UV tree-level, after horizon-crossing.A quasi-dS universe is probably the safest case in this regard, which is convenient for its believed physical relevance.Some previous work [26] has even managed to account for suppressed next-order statistics in the noise.More generally, this stochastic evolution should be reliable for studying any non-linearity and non-gaussianity generated by the evolution that is larger than the first order in CPT.However, if no nonperturbative effect is obtained from SI, then one should stick to QFTCS, valid on IR scales and beyond tree-level before horizon-crossing, which is more precise than SI, as the disagreement between the two shows in perturbative regimes [35].
In this section and following the literature, the stochastic backreaction (different from the classical Einstein backreaction) has also been neglected because we have linearised the theory around the homogeneous background and not around the IR quantities by invoking the Starobinsky approximation X i> = X i b + O(δX i ) [20].This includes the calculus of the spectrum of {δX j ⃗ k } j .This assumption is crucial to get closer to a Markovian system and at least an additive noise, thus facilitating analytical solutions.
However, different heuristics have been used in the past for the stochastic update of those RHS and of the background quantities on which the modes evolve.Usually, background quantities X i b in the RHS of eqn.(40) and in the equation of {δX j ⃗ k } j are taken as their local IR versions X i> (t, ⃗ x).In this case, solving the equations requires a numerical approach [35,54,55].These studies have all reported a very small impact so far.In this work, we will keep it arbitrary and keep any RHS amplitude factor such as H or ε i undefined in this respect, until further comment.In the linear limit, one always retrieves eqn.(37).
To summarize, stochastic inflation is about reducing QFTCS approximations to UV scales only, allowing for the possibility to study fully non-linear and nonperturbative phenomena on IR scales, i.e. above Hubble scales in our case.The Langevin-type EOM of stochastic inflation can also be mapped to an associated Fokker-Planck (FP) equation and both can be solved either analytically in some simple cases or numerically.From these solutions, one can then derive non-perturbative correlators or even the full probability functional of the fields.It is in this ability to provide results beyond perturbation theory or even fixed n-point correlators where Stochastic inflation's appeal lies.

B. Current approaches to Stochastic Inflationary dynamics
So far we presented a schematic picture of how one might obtain equations describing the dynamics of Stochastic Inflation from coarse-graining the equations of motion.However, applying the above scheme to the full equations of General Relativity has not been implemented due to the latter's relative complexity, the a priori large number of variables, which can be both dynamical and constrained, and the necessity of making coordinate/gauge choices.As a result, various approximations have been made, mostly taking the long wavelength limit and reducing the number of dynamical variables.
Early work focused on the coarse-graining of the Klein-Gordon field equation only, assuming an unperturbed (initially dS) background and coarse-graining the field and sometimes its conjugate momentum [20,49,[56][57][58].It was only in [27,28] that the first stochastic equations were formulated fully including metric perturbations and so backreaction.This approach is still widely used.It consists in decreasing the interdependence of GR equations by using the long wavelength approximation and judicious gauge choices.In the full ADM formalism, even before linearization and coarse-graining, the lowest order gradient expansion of eqn.( 9) and (12) in the convenient β i = 0 slicing becomes [27,36,59] together with the constraints Furthermore, Si j is usually set to zero in the absence of anisotropic fluid sources, or because any possible contributions to it are considered higher order in the gradient expansion, leading to an exponential decay of Ki j .This results in one equation less with only K and Π dynamics remaining, that is, These equations form what is called the separate universe approach (SUA), widely used as a basis for the formulation of stochastic inflation.In many later works stemming largely from e.g.[29,60], the momentum constraint is not considered part of the long wavelength approximation, an approach recently referred to as (k=0)-SUA in [61].The relevant literature then interprets literally the similarity of the Π, K and the Hamiltonian constraint in eqn.(43) with Friedman's background equations of cosmology eqn.(18).The local quantity 1  3 K(x) replaces the homogeneous Hubble parameter H, so that the whole inhomogeneous universe is described as made up of patches, each of which evolves independently.The EOM coarse-graining has largely been applied to those two equations only: linearising provides the windowed RHS terms which can then be evaluated using UV modes solutions.
Unlike scalars on a fixed background, the inclusion of metric perturbations brings about the issue of the appropriate slicing to be used for both the linearisation and the noise calculation, especially beyond slow-roll scenarios.In this respect, the uniform-N gauge slicing (see sec.III B 4) has emerged as the most natural, with [30,31] arguing early on that in this slicing one recovers CPT equations and thus the long wavelength QFTCS limit when linearising the overdamped field equations.In general, the uniform-N gauge has emerged as the more commonly used because it allows direct access to the statistics of the fields in terms of N b slicing, a necessary step for the first-passage-time analysis (FTPA) and the associated stochastic ∆N formalism which provides information about the non-linear curvature perturbation on uniform field hypersurfaces [37].
When performing the EOM coarse-graining in this gauge, one gets the SI equations [32], which we will refer to as the (k=0)-SUA, adopting the term from [61], or the 'usual' equations where N b is the background number of e-folds, ε > 1 = −∂ N b ln H > the first slow-roll IR parameter, and In particular, this means that the coarse-graining has been assumed to apply as δϕ > = W δϕ and δπ > = W δπ in Fourier space [32]. 3By defining π > with a stochastic source, the literature does not use the canonical momentum of ϕ > .One can also re-write the previous equations as a unique equation in the ADM Π > , with , where H is the background quantity or not depending on the scheme used, see sec.IV A. In that sense, it is perfectly equivalent to coarse-grain the momentum and its definition (phasespace coarse-graining), and to coarse-grain the momentum via the field only, which will be our approach, following the coarse-graining in Wilsonian EFTs [23].It is common to keep the phase-space equations to apply the common overdamping assumption |δΠ| 2 ≪ |δϕ| 2 and thus reduce the problem to a first-order equation of the field only [20].This assumption will not be made in the following.
The role of the momentum constraint has been noted (refer e.g. to [27,35,36,59,61,62]) as being the only remaining link between the separate universes in the SUA.Neglecting it by invoking the long-wavelength approximation could lead to inconsistencies in CPT [35,62,63].This issue has been addressed in recent years by tackling the regime of validity of the (k=0)-SUA by comparing to CPT [32,33]: it seems that working far enough from the horizon-crossing scale with the uniform-N gauge is a safe choice to match CPT.Nevertheless, the notable results achieved by this approach [37,64,65] require transformations between gauges (flat gauge for the noise, uniform-N for the Langevin equations, and effectively uniform-ϕ in the FTPA) when giving up approximations such as slowroll [32].
Other SI approaches have emerged in parallel with somewhat different assumptions.Starting from the Hamilton-Jacobi formalism developed in [27] from eqn. (43), SI was formulated with the inclusion of the momentum constraint in [28], and the FPT approach of [37] was applied to this formalism in [36,66].Starting from the full set of (43), a set of stochastic equations for nonlinear variables like those given in (67) were developed in [34] which by construction reduce to coarse-grained, gauge invariant CPT in all time-slicings with β i = 0.As shown in [67,68], these equations produce the same perturbative results for the bispectrum as the ∆N .The more recent study [35] also writes stochastic equations in a gauge other than uniform N , the uniform Hubble gauge, by rehabilitating the momentum constraint in the gradient expansion.
The mini-literature review discussed above attempts to provide a flavour of the current state of play regarding the status of SI formalism and the approximations that have been deemed necessary to develop it.In the next section, we demonstrate how one can do away with the long wavelength approximation and the SUA, providing a set of SI equations that retain both the scalar field and all the variables of the gravitational field, the only approximation being scalar linear CPT for the computation of the noise terms by coarse-graining.

C. Stochastic equations for ADM
We now incorporate the continuous influx of modes that cross a smoothing scale commensurate with the Hubble scale, i.e. we apply the schematic procedure leading up to eqn.(40) to the full set of Einstein equations.This section contains the key results of this work, namely a computation of the stochastic sources associated with the ADM formulation of General Relativity.These stochastic equations are presented in (56).We perform this computation in all the gauges discussed above, always finding the same result for the stochastic source terms, presented in eqns.( 47) and ( 57); as can be seen there, the source terms are always given by the same functional of the gauge invariant R and the chosen window function.The stochastic source terms appear on the RHS of the dynamical equations but not the constraint equations.We stress that we do not impose any gradient or slow-roll expansion.

Coarse-graining linear theory
To coarse-grain GR we will apply the principle behind formula (40) to determine the stochastic source terms from linear theory.As explained above, the major problem is choosing a gauge, its associated dynamical quantities, and the window.As we will see, in this section and the next one we perform a gauge-invariant coarsegraining, i.e. we coarse-grain the gauge-invariant linear theory encoded in the usual gauge invariant quantities of eqn.(22).In particular, the coarse-graining is made by using a time-dependent window function in Fourier space W k (t) on a gauge-invariant quantity, here R. For the case of inflationary evolution, W k (t) would be activated after Hubble radius crossing when quantum modes can be assumed to behave classically.However, we stress that the coarse-graining method remains valid beyond this choice if one accepts operator-valued source terms that cannot be interpreted fully as classical random variables.
In practice, and as explained schematically through eqn.(37), coarse-graining means that we first search for the equations of the long-wavelength variable R > k = W k R k where R obeys (23).R > is still gauge-invariant of course but follows a slightly different equation of motion which can be obtained from eqn. ( 23) where the source term is We emphasize that in (47) the > operator has priority over time derivatives in all equations, i.e. for function ḟ > ≡ d (W f ) /dt.Note also that the source term S R vanishes when the window is constant, i.e. usually when it is super-Hubble or sub-Hubble enough, the latter yielding R > ≃ 0. This also shows that the existence of the source term is solely attributed to the time-dependent nature of the UV-IR split.Next, we derive the coarse-grained version of all the metric variables and the scalar field as given in sec.III B by using the replacement before taking any time derivative.This is done for any of the gauge families discussed in that section.We stress that this prescription ensures that the constraint equations will stay perfectly satisfied at first order, which we have verified for all the gauges we examined.

Coarse-graining General Relativity
The expressions obtained after making the replacement (49) in the formulae of sec.III B can be substituted into the linearized evolution and constraint equations of section II B for any of the aforementioned gauges.As explained above, this procedure will provide the source terms for each of the ADM equations.Although, as we will see below, the final result is very simple and identical in all gauges, the full computation involves rather long expressions.We will therefore provide here a summary of the full derivation of the coarse-grained field equation in a small generalised synchronous for illustration.All the other stochastic ADM equations are obtained similarly.
We need to calculate the RHS term of the following linearized equation for the field perturbation: ) which is just the re-writing of eqn.(20), using the dictionary of eqns.( 14), (15) in Fourier space.In this gauge, all relevant coarse-grained perturbation quantities are given in terms of R > as together with the coarse-grained Bardeen potentials where the operator > has priority over time derivatives.This non-zero difference of the long wavelength Φ > B and Ψ > B is not due to anisotropic stress but simply due to the appearance of a time derivative in the definition of Ψ B in (22); when coarse-grained Ψ B acquires an extra stochastic source term compared to Φ B .It is a transient horizoncrossing effect.This comes back to the usual equality for each mode if super-Hubble or deep sub-Hubble because the correction is negligible when the window function is constant.
Using the previous decomposition, further derivatives are needed where V ,ϕϕ is the background quantity, function of ε 1 , ε 2 and ε 3 , see Appendix A. We now have all the coarsegrained variables needed to plug into the LHS of the coarse-grained field equation (50).After a lengthy computation, a perfect cancellation occurs and only the S R term of δϕ > survives, giving which is completely independent of Ψ * and B * , i.e. of the specific functions corresponding to the choice of small gauge in which the computation was performed.Following the EOM approach we described in IV A, we can now assume that stochastic source terms derived for the linear spectrum amplitude equations remain valid for their non-linear parent equation so that the final, non-linear super-Hubble Langevin equation in this gauge is and where any LHS variable is implicitly understood as long wavelength, i.e. we've suppressed the the notation > .This whole procedure can now be performed to compute the new right-hand side (RHS) for any dynamical quantity in the ADM formulation and in any of the gauges discussed.Note that the computations are cumbersome and have been checked with Mathematica [69] which was also used to confirm that our previous R-decompositions were perfectly satisfying the original first-order equation eqs.( 19), ( 20) and ( 21) before applying the window function W k .
Applying the procedure outlined above to all the rest of the ADM equations, we find that they get augmented by Brownian terms as shown below where the RHS > superscripts are implicit, the source terms given by and with S R from eqn. (48).Importantly, any terms that might contribute to the source terms on the RHS of the constraints cancel completely and hence This a physically appealing result and a consequence of our coarse-graining philosophy: the stochastic noise terms can be thought of as a continuous readjustment of the "initial value data" at each timestep of the dynamical evolution.When setting up initial conditions in the ADM formalism, all relevant fields (determined by the choice of gauge) must be specified such that the constraints are satisfied (here up to O(S 2 R )).The above equations therefore ensure that this also remains true for each such stochastic readjustment per time step.
We have performed the computation in all families of gauges previously defined, always finding the same result: the Fourier transforms of the RHS coincide for all gauges.We therefore postulate that the above equations must hold for any arbitrary gauge choice even beyond those discussed here.Finally, following the above approach we have also coarse-grained both the more common formulation of Einstein's equations as well as their BSSN incarnation [44], see Appendix C. The results are consistent with eqn.(56).Writing down the BSSN equations is a necessary step to study well-posed GR in numerical relativity.

Discussion
Let's first recall that those equations' RHS are only valid for perturbations around a homogeneous background.In particular, without a separate universe approach, we don't provide here any heuristic to account for stochastic backreaction.
The most striking observation is probably the simplicity and similarity of the RHS terms, which contrasts with the NL LHS.This is completely due to the linear framework in the UV and the fact that all dynamics are encoded in one variable.This is also supported by the numerous null RHS in equations that encode either field definitions or constraints.In particular, the perfect satisfaction of the constraints after coarse-graining is a good sign that our spacetime is physical, i.e. here the horizon crossing is done coherently on the whole time hypersurface but also that we have addressed the insertion of the window.When setting Ψ > B = Φ > B one gets a violation of the constraints and a strong gauge-dependence of the RHS, which supports the choice of eqn.( 52) for an appropriate coarse-graining.
Another interesting term is that of the anisotropic evolution equation: first-order scalar perturbations source tensorial quantities at higher orders.This is in agreement with previous work [70,71] and is discussed further in section IV E.
Furthermore, the previous equations are valid for four major families of small gauges and any background time slicing (the latter being only a matter of variable change for straight time derivatives to get to α ′ b (t ′ ) ̸ = 1).This suggests that those equations could be true for any small gauge.At a linear level first, taking, for instance, eqn.(50): a 1st order gauge transformation would leave the LHS unchanged and similarly for the RHS as it is written in terms of a gauge-invariant quantity and background time derivatives.This works even if the transformation is a function of the window or R. The gaugeinvariance of the fully non-linear equations is less obvious.In particular, a gauge transformation will leave any non-linear LHS unchanged but the associated RHS will be consistent only if both the initial and the final gauges are close enough to the homogeneous background one so that we stay within Perturbation Theory when linearising the RHS.This is why there is evidence of small-gauge invariance of our equations at first order.
One needs to emphasize that the long-wavelength approximation has been removed from the equations.This is what allowed us to formulate equations in any gauge by bypassing the gauge-mapping issue in this regime [33].Note that the long wavelength approximation could still be applied to the choice of the window function if one wants to ensure the complete classicalisation of the crossing modes [52].However, it seems plausible that the window function can now be turned on much closer to the Hubble radius than in the (k = 0)-SUA, the latter being restricted by the quasi-isotropy assumption [33].In particular, one can now study safely regime transitions in SI where gradients are critical, as opposed to the main approach [72].Switching on the window closer to Hubble crossing also gives less interaction time (and so less higher-order effects) to the UV modes and so strengthens the linear and Gaussian source approximation.Of course, more study is needed before verifying this assertion, something we leave for future work.
Related to this matter, we finally want to discourage any attempt to go too far away from a dS spacetime.Although it is true that our derivations do not make any slow-roll assumption, using Gaussian sources, linear CPT, and no stochastic backreaction (i.e.full accountancy of UV-IR interactions) might not encode all necessary contributions from UV modes and thus questions the whole story of SI.To go in other regimes, one could study the UV within QFTCS and show the presence of a hierarchy in the correlation functions, order by order.

D. Usual Stochastic Inflation limit
In this section, we compare our results with eqn.(46) of SI [32].To achieve this, we need to change the LHS of the Π equation to its SUA limit eqn.(41) and specify our gauge.
By fixing ∂ t N = 1, i.e. choosing t as the background efolding N b , and β i = 0, we can write α > , and hence α > is given by the Hamiltonian constraint in the long wavelength limit.Substituting this in the Π equation's LHS yields which is equivalent to the usual eqn.( 46) but with our own noise term F −1 {S Π } at first order.Note that since the S R has no apparent Laplacian, we chose to leave it unchanged under the long-wavelength approximation.
Our formalism is now comparable to the literature's SI by looking at Σ and S R only.On the one hand, S R can be re-written in terms of background efolds slicing (α b = H −1 and H∂t = ∂N b ) as (60) On the other hand, we can make Σ explicit by calculating δϕ and ∂ N b δϕ from the long wavelength limit of our calculus in sec.III B 4 where δN * = 0.By using the long-wavelength limit4 which confirms to be the same as the spatially flat gauge decomposition.Finally, by substituting this in Σ ϕ and Σ π of eqn.(45), one exactly gets We just confirmed that our equations give the eqn.(46) when using the same assumptions.

E. Stochastic gravitons
As already pointed out, the scalar perturbations have an influence on the non-scalar degrees of freedom through the Kij evolution equation of eqn.(56).At second order and later non-perturbative orders, one thus expects scalar-induced and scalar-coupled gravitational waves [70,71].In that sense, eqn.( 56) is the first of its kind to provide a stochastic framework including both scalar and tensorial evolutions.This is not surprising as previous studies worked in the long-wavelength limit to drop tensorial dynamics.
In previous stochastic inflation work [73], it is suggested that we should also include the decohered first order gravitons in the stochastic sources.This can be simply added if we stay at linear order for the sources.Indeed, we only need to do the same work as for the scalar but for gravitonic perturbations.This can be done independently and just added to the final scalar result thanks to the scalar-decoupled limit.In practice, this appears to be straightforward as we can build a linear gaugeinvariant quantity h by writing the linear tensorial part of the metric in cosmic time as with h i i = ∂ i h ij = 0 (traceless and transverse).At linear order, K, Kij , 3 R and 3  Rij are the only quantities in the EOM inheriting contributions from h.This leads to the only equation in (56) having h terms at linear order, the linearised Kij equation which in Fourier space appear to be the known Mukhanov-Sasaki equations for the two linear spin-2 components of the gravitons which can also be initialised by a Bunch-Davies vacuum.
We now have everything we need to coarse-grain at linear order with the previous window method using ), assuming we want the same window for both polarisations.By coarse-graining eqn.(64) to get the spectrum of sources for eqn.(63), and promoting the sources to those of the Kij equation, we get to update S Kij as where The coarse-graining is now complete and appeared to be much easier because of the gauge invariance of the perturbation and the trivial satisfaction of the constraints at linear order.This extension is not without utility because eqn.(65) shows a competition between scalar and tensor sources, which seems to be in favor of the latter in slow-roll regimes.Tensorial degrees of freedom can now be studied in a non-perturbative framework, most likely numerically in the future.

V. LANGEVIN EQUATIONS IN THE UNIFORM FIELD GAUGE
The stochastic equations of the previous section can be applied in a variety of gauges.It is now time to extract useful (gauge invariant) quantities such as curvature scalars on certain matter hypersurfaces.In the linear theory, it refers to ζ gi and R defined earlier, curvature perturbations on spatial hypersurfaces of uniformdensity and uniform-field respectively.Beyond the first order in CPT, it is possible to construct such non-linear gauge invariant quantities although the literature provides different levels of assumptions [29,59,74].One can in particular define the following ones [59,74] where N was defined above in sec.III B 4. It is a gaugedependent quantity as the determinant is a density-2 tensor.According to [74], the usefulness of (67) holds even beyond the long wavelength approximation.
From these variables, it is clear that knowing N in uniform-density or uniform-field gauges provides a direct gauge-invariant extraction and this is why the FPTA is required in usual studies where the efolds are not stochastic [37].In this context and since a gauge-invariant formulation is being proposed in this paper, we propose to apply the uniform field gauge directly to our equations and avoid the FPTA.Note that this has been attempted recently but starting from the literature's usual equations and assumptions [75].
We start from our equations (56) and set the coordinates such that the field evolves uniformly according to the background dynamics, ϕ(t, ⃗ x) = ϕ b (t), where ϕ b (t) follows eqn.(18) with α b = 1.Hence ϕ can act as a clock labeling the 3D spatial hypersurfaces.The spatial coordinates on these spatial slices are fixed such that β i = 0. Obviously, this gauge choice can only be valid if the background field is non-static.In addition, ϕ does not now receive stochastic impulses and the stochastic dynamical variables are N and α for which Langevin equations can be derived.
Leaving a more complete numerical study for future work, we make in this section the long-wavelength approximation for the Hamiltonian constraint, re-writing eqn.(43) by inserting the definition of N .From this we immediately obtain by making use of the background Friedman equation.
No slow-roll approximation has been made so far.The second equation we need is given by the field equation which becomes in this gauge an equation for the lapse where S Π is the Langevin noise found previously in sec.IV C 2. Note again that any multiplier of S Π is left undecided concerning the stochastic backreaction.To get to ∂ t α = ∂ t f where f (Z t , t) = φb Zt for Z t = φb α , the Ito lemma [76][77][78] is used to get We have included the Ito correction to the derivative of a function of a random variable as the second term in the above equation.Note however that the effect of this term is minimal (higher order) and the results of the simulations we describe below are practically unaffected by it.Substituting the field background equation in and the efolds evolution from eqn. (69), changing the time variable by dividing by H and consequently updating the variance of the Wiener process yields where S = S Π /H φb is the final stochastic source.In the following, we will keep the coefficient α 3 multiplying S as a non-background, stochastic quantity, and won't set it to 1.Note that when linearized, these equations still match CPT at first order.Equations ( 69) and ( 72) together form a coupled system of stochastic PDEs.
In the SUA philosophy, it is common to restrain our case to one patch of the universe for the treatment of the noise, i.e. using the previous equation at one given point ⃗ x 0 and using only S(N b , ⃗ x 0 ).This is completely justified because the window function derivatives yield suppressed correlations beyond the Hubble scale (e.g. a Heaviside window in Fourier space would lead to a cardinal sinus in real space) [20].In this framework, it is easier to solve this system of two coupled Langevin equations or stochastic ODEs.
An analytical solution would require specifying the background dynamics and writing the Fokker-Planck equation to get the PDF.To our knowledge, there is no known way to do that analytically without further approximations such as overdamping, model-dependent simplifications, or higher-order correlations neglection.
As a proof of concept, we decided to provide numerical results instead.
The simulations were realised with the Stratonovich evolver of Mathematica [69].The following amplitude for S, computed in Appendix D, is valid for both slowroll and ultra-slow regimes (USR) The main difference with usual amplitudes is the √ ε 1 in the denominator.This comes from the fact that the lapse equation depends on the amplitude of R and not the usual Mukhanov-Sasaki variable (a √ 2ε 1 M P l × R) to which field equations are sensitive to.Of course, the ε 1 → 0 limit is problematic but this is purely a gauge artifact as already noted earlier when defining our gauge and as explained in [79].For the simulation, we choose to evaluate this quantity with the background evolution and so to neglect the stochastic backreaction for simplicity, as opposed to the α 3 factor acting on S.
Figures 1 and 2 show 50.000paths starting from uniform lapse and efolds.0.01-efolds steps were used in the numerical scheme.Without worrying about the realism of our model of Inflation concerning observational constraints, we considered a plateau Inflation followed continuously by a 3rd order polynomial slope of the form . Initial position and momentum together with the start of the slope are finetuned 5 for two reasons: we want a decent amount of dispersion and non-Gaussianity by having the field almost stopped on the plateau, but not too much as the gauge definition makes the simulation crash if the velocity is too low.Constraints also need to be satisfied.It appears that USR is reached from 2 efolds onwards (ε 1 ≪ 1, ε 2 ≃ 6) which is why the noise is only turned on at this time.When the slope is reached at ≃ 3.04 efolds, the velocity starts to increase extremely slowly, still within the USR regime.
It is well known that a flat potential or a transition can leave strong non-Gaussian imprints on perturbations such as exponential tails [65].This is actually what we confirm here in Figures 3 and 4: both the lapse and the efolds get an exponential tail on the plateau between 3 (end of plateau) and 4 efolds, which can be confirmed by diverse fittings.This non-Gaussianity can also be tracked in time by looking at the skewness and the kurtosis in Figure 2.
When it comes to the second phase -and any analytical attempts would probably fail to describe it fully, things eventually stabilise along the slope.From these figures, it becomes clearer that the lapse acts as an extremely non-Gaussian efolds' momentum.This implies that very strong non-Gaussian changes are given to the efolds until stopping and eventually the efolds distribution's non-Gaussianity stabilizes later to a lower remnant level6 when all realisations are in the same regime.We have checked that switching the noise terms off at 6 efolds makes the lapse come back to its attractor α b = 1 and that the efolds, which are meant to describe R N L , are indeed conserved.If the noise was frozen later on this potential, one would see that leaving the USR phase makes the lapse and the efolds back closer to Gaussianity.This is because PDF in real space is not a good estimator: adding many Gaussian contributions lowers the relative non-Gaussianity.However the Non-Gaussianity we produced is still imprinted by the end of our simulation, in particular on the statistics of the scales which crossed the Horizon before 6 efolds.This advocates for the necessity to look at quantities such has n-spectra or coarse-grained PDFs when looking at data.
It is important to mention that the simulation has also been run by setting α 3 to α 3 b = 1.Unusual left-skewed PDFs have been generated and highlight the importance of stochastic backreaction in certain cases.Here it is critical because of the non-perturbative behavior of the lapse and the efolds.In particular, this term is a good barrier to reaching an unphysical α = 0 with stochastic kicks, without adding a prior when solving the Fokker-Planck equation.
What is to be remembered from these equations and simulation is their potential: having coarse-grained but satisfied constraints and an idea of what the stochastic amplitude is, can help us skip the tedious and approximate stochastic ∆N and FPT formalisms by writing it all directly in the right gauge.
However, sticking to the literature, the SUA was used and so the validity is far from the crossing scale.This allows us to reduce to a simple GR framework where we don't need a full Numerical Relativity code which probably does not exist in such a gauge.Furthermore, the validity of the long-wavelength approximation is questionable when perturbations become smaller in magnitude than gradient corrections.This could be the case here for the lapse.For these reasons, traceless modes and other terms should be fully accounted for and could stop those perturbations from vanishing completely.We thus emphasize here that it is necessary to run a full NR code even if it takes to use an NR gauge.In that sense, this section mainly aims at illustrating how this new framework might compete with the FPT formalism using the same assumptions for the evolution.

VI. CONCLUSIONS AND DISCUSSION
In this work we have examined the formulation of SI in full General Relativity, dropping most of the approximations made historically such as the long wavelength approximation and the corresponding reduction of gravitational field variables, such as anisotropic degrees of freedom or the momentum constraint, that have been integral to existing versions of the SUA.We have also addressed the issue of time-slicing choice; the choice of lapse and shift in the language of ADM.The only essential approximation kept here is the requirement that CPT is enough to compute the noise source terms, a requirement that could be possibly lifted by retaining higher orders in CPT.Although we have been invariably discussing 'stochastic' noise terms, the coarse-graining behind our derivations could have been made at any scale, provided one was content with operator-valued source terms.Of course, the classicalization of cosmological perturbations is highly convenient, and it is this requirement that sets the IR scale to be placed at some slightly above the Hubble radius.
Returning to the core of 3+1 general relativity, which is essentially the evolution of one temporal 3D-hypersurface to the next, we have proposed a method to coarse-grain the linear theory in a gauge-invariant way using its only degree of freedom, here chosen to be R the comoving curvature perturbation.We have validated our procedure for several of the most common gauges used in cosmology or numerical relativity, finding that the choice of any of them always results in the same source terms, at least at the linear level.Linking the resulting Langevin terms to their counterpart non-linear equations and adding the treatment of stochastic gravitonic sources, we have provided the first complete set of GR equations for Stochastic Inflation.Looking at their form and exploring alternative gauges has offered strong evidence for our postulate that they are indeed gauge-invariant in that any gauge choice would provide identical results.
From our equations, we were able to recover the limits where existing SI equations apply.We were also able to go beyond the usual approaches to demonstrate broader applications; our example showed that our formalism could obtain results for the stochastic dynamics of the e-folds N directly in the uniform field gauge, without shifting 3D-hypersurfaces as required in the stochastic ∆N formalism, for the same SUA assumptions normally made in the existing literature.
We want to highlight the potential of such results.The possible applications are numerous and most notably a key focus of our future work will be numerical.The present article provides all the tools needed for sourcing a full 3+1 numerical relativity code (most of which are running with BSSN equations) with stochastic perturbations.Note that this includes the Langevin terms but also the initial conditions (see sec.III B 4) which constitute the main challenge of numerical relativity.The present systematic treatment gives the opportunity to quantitatively study the nonlinear evolution of inhomogeneities, taking forward previous work which considered inflationary initial conditions in a variety of contexts [82][83][84][85][86]. Full GR simulations of super-Hubble dynamics during inflation should provide insights about the nonlinear generation of higher-order correlators.This is important because currently there is no alternative to QFTCS methods except Stochastic Inflation, which has to date traded greater scope on super-Hubble scales in exchange for numerous other approximations.
Finally, we note that there remain many further extensions to be considered for Stochastic Inflation.For instance, our methods should be compatible with a greater range of inflationary scenarios such as multiple fields [34], accounting for anisotropic sources from these, or modified and higher energy theories of gravity in the context of the EFT of Inflation [87].Note that these would require extra-work to find well-posed formulations, enabling numerical solving.It would also be interesting to investigate non-quasi-dS spacetimes, though this would either require rigorous justifications or higher-order perturbations and statistics.With or without quasi-dS scenarios, the incorporation of higher-order effects also lacks a full GR framework.As stated in [23,24,58,88] in fixed dS spacetimes, we believe that the most rigorous approach would begin from a path integral approach rather than the EOM, however, efforts should be made towards a full GR framework [48].In particular, the coarse-graining approximations would be under control and this could enable the incorporation of quantum loops, clarifying the validity of Starobinsky's approximation with all gravitational degrees of freedom.This is left for future work.

FIG. 3 .
FIG.3.PDFs (blue) of a scalar field undergoing USR dynamics in uniform field gauge, see Figure1for parameters of the simulation.50 000 paths were used for Kerner Density Estimation to probe the PDFs at five background efold times (N b = 3, 4, 5, 6 from top to bottom) for both the lapse (left) and the efold (right) differences to the background, against normal PDFs with the same first two moments (black dashed).95% confidence intervals (red) of these pdf estimators were calculated using bootstrapping.