Stochastic Inflationary Numerical Relativity

• Stochastic Inflationary Numerical Relativity is a framework that splits scalar and tensor fields into IR and UV components to incorporate gauge-invariant stochastic evolution in full 3+1 dimensions.
• It integrates nonlinear gravitational backreaction, gradient retention, and non-Gaussian statistics to simulate extreme inflationary phenomena like primordial black hole formation.
• The methodology employs high-order ADM/BSSN evolution with stochastic source terms on a discretized cubic grid ensuring robust constraint preservation and realistic initial conditions.

Stochastic inflationary numerical relativity combines the stochastic (quantum diffusion) approach to inflationary cosmology with the dynamical simulation of spacetime and matter fields using full 3+1-dimensional numerical relativity (NR). This framework realizes a gauge-invariant, nonperturbative theory of inflation in which all scalar and tensor metric degrees of freedom, as well as the inflaton and its stochastic fluctuations, are evolved with the complete Einstein equations. Stochastic inflation in numerical relativity extends standard stochastic inflation by incorporating nonlinear gravitational backreaction, retaining gradient information, transcending the separate-universe approximation, and enabling the study of extreme, nonperturbative inflationary phenomena such as primordial black hole formation, strong-resonance effects, and non-Gaussian statistics.

1. Formulation of Stochastic Inflation in General Relativity

The core of stochastic inflationary NR lies in the split of all scalar (and tensor) degrees of freedom into long-wavelength (infrared, IR) and short-wavelength (ultraviolet, UV) parts using a dynamically evolving coarse-graining window $W_k(t)$. For any field $X_k$, the split is defined by

$$X_k^{\rangle} = W_k X_k \qquad X_k^{\langle} = (1-W_k) X_k$$

with $W_k = \Theta(\sigma aH - k)$ for a sharp filter (alternatively a smooth function can be employed). The fundamental stochastic variable is the comoving curvature perturbation $\mathcal{R}$. All other metric and scalar variables can be written in any gauge in terms of $\mathcal{R}_k$ and its derivatives, exploiting the gauge-invariant structure of linearized perturbation theory and extending nonlinearly through their embedding into the full ADM or BSSN systems (Launay et al., 16 Dec 2025, Launay et al., 2024).

Crucially, the continuous flow of UV (quantum) modes into the IR is captured as a stochastic noise term. For the IR sector, the evolution equation acquires a stochastic source: $$\ddot{\mathcal{R}}_k^{\rangle} + H(3 - \epsilon_2)\dot{\mathcal{R}}_k^{\rangle} + \frac{k^2}{a^2} \mathcal{R}_k^{\rangle} = S_{\mathcal{R}}(k,t)$$ where $S_{\mathcal{R}}(k,t)$ is constructed from time derivatives of the window $W_k$ and the slow-roll parameters. Upon substituting the gauge-invariant IR variables back into the full nonlinear ADM or BSSN equations, the system becomes a coupled set of stochastic partial differential equations for all relevant dynamical fields.

2. Stochastic ADM and BSSN Evolution Equations

In the 3+1 decomposition, the ADM line element

$$ds^2 = -\alpha^2 dt^2 + \gamma_{ij} (dx^i + \beta^i dt)(dx^j + \beta^j dt)$$

is evolved with dynamical variables lapse $\alpha$, shift $\beta^i$, 3-metric $\gamma_{ij}$, extrinsic curvature $K_{ij}$, scalar field $\phi$, and its conjugate momentum $\Pi$. Stochastic sources derived from $S_{\mathcal{R}}$ are injected only into the dynamical equations; the Hamiltonian and momentum constraints remain unsourced. The principal stochastic RHSs in Fourier space are: $$\begin{aligned} \mathbf{S}_K(k,t) &= -\epsilon_1 S_{\mathcal{R}}(k,t) \, \boldsymbol{\alpha}_{\vec k},\ \mathbf{S}_{\tilde K_{ij}}(k,t) &= \epsilon_1 \left(\frac{1}{3} \delta_{ij} - \frac{k_i k_j}{k^2}\right) S_{\mathcal{R}}(k,t) \, \boldsymbol{\alpha}_{\vec k},\ \mathbf{S}_\Pi(k,t) &= \sqrt{2\epsilon_1} M_{\rm Pl} S_{\mathcal{R}}(k,t) \, \boldsymbol{\alpha}_{\vec k} \end{aligned}$$ with $\boldsymbol{\alpha}_{\vec k}$ complex Gaussian random variables with unit correlators.

For numerical stability, these equations are recast in the BSSN (Baumgarte-Shapiro-Shibata-Nakamura) formulation. The BSSN evolution variables $$(X, \tilde{\gamma}_{ij}, K, \tilde{A}_{ij}, \tilde{\Gamma}^i, \phi, \Pi)$$ evolve according to the standard first-order equations, with stochastic driving on $K$, $\tilde{A}_{ij}$, and $\Pi$ only (Launay et al., 16 Dec 2025, Launay et al., 2024). All constraint equations remain zero-noise.

3. Numerical Implementation: Discretization, Noise, and Stability

Implementation for 3+1 numerical relativity adopts a uniform cubic grid with $N^3$ points, typically $N=256$, and periodic boundary conditions. Spatial derivatives are discretized with fourth-order finite-difference stencils (5-point for gradients, 7-point for Laplacians). The time integrator is a fourth-order Runge-Kutta scheme with Courant ratio $\Delta t/dx \lesssim 0.4$ ensuring CFL stability (Launay et al., 16 Dec 2025).

Stochastic sources are generated and injected at every Runge-Kutta sub-step via the following workflow:

1. Evolve the real-space mode equations for $\mathcal{R}(\vec{x})$ and $\Pi_{\mathcal{R}}(\vec{x})$.
2. Compute their Fourier transforms, evaluate $S_\mathcal{R}(k,t)$ for each $k$ using the window derivatives, and multiply by appropriate $\alpha_k$.
3. Inverse Fourier transform to obtain real-space noise fields $S_K(\vec{x})$, $S_{A_{ij}}(\vec{x})$, $S_{\Pi}(\vec{x})$ to be added to the PDE right-hand sides.

Noise in time is implemented as white noise via Gaussian draws with variance $$\langle (\Delta W)^2 \rangle = \Delta t \times {\rm Amplitude}^2$$; in space, the window function $W_k$ ensures correlation length $\sim (\sigma aH)^{-1}$. Constraint enforcement is monitored at each step by calculating relative violations for the Hamiltonian and momentum constraint densities; violations remain $\lesssim 10^{-10}$ in slow-roll and $10^{-2}$ in ultra-slow-roll crossings, scaling as $O(S_\mathcal{R}^2)$ (Launay et al., 16 Dec 2025).

4. Inflationary Potentials, Gauge Choices, and Initial Conditions

Simulations have been performed for a range of inflationary potentials:

• Quadratic ($V(\phi) = \frac{1}{2}m^2\phi^2$),
• Inflection-point type,
• Strong-resonance (e.g., monodromy-induced features).

Initial conditions are set by classical realization of Bunch-Davies vacuum perturbations: for each mode, the quantum ladder operators are replaced by complex Gaussian variables, matching the two-point spectrum. The metric and inflaton field initial data are constructed so that the ADM Hamiltonian and momentum constraints are satisfied exactly to linear order, using algebraic relations between $\mathcal{R}$ and metric/matter perturbations in the specified gauge (usually geodesic or BSSN-compatible gauges). This procedure ensures realistic stochastic inhomogeneities at the simulation start, capturing the correct initial power spectrum and correlations (Launay et al., 10 Feb 2025).

5. Regimes of Evolution and Comparison to Linear Theory

Stochastic inflationary NR has been tested in both slow-roll (SR) and ultra-slow-roll (USR) scenarios, revealing the following behaviors (Launay et al., 16 Dec 2025):

• In SR, background quantities (e.g., Hubble parameter, slow-roll parameters) track the homogeneous Friedmann solution to $<0.1\%$, matter contrasts freeze at $O(10^{-5})$, and constraint violations remain insignificantly small.
• In USR, matter and curvature perturbations can grow to $O(1-10)$, enabling exploration of highly nonperturbative regimes inaccessible to standard cosmological perturbation theory.
• Power spectra evolution for $\mathcal{R}$ matches the linear Mukhanov–Sasaki result in SR and grows dramatically in USR, consistent with analytic expectations.

The formalism retains full gradient information and non-Gaussian statistics, enabling the study of rare, large fluctuations (e.g., peaks in $\Delta_\mathcal{R} \gtrsim 1$), nonperturbative density inhomogeneities, and initial data for subsequent radiation or preheating phases (Launay et al., 16 Dec 2025, Launay et al., 10 Feb 2025).

6. Extension Beyond Separate-Universe and Applications

Standard stochastic inflation in cosmology often employs a "separate-universe" approximation (ignoring spatial gradients) and a gauge-fixed $\Delta N$ (number of e-folds) prescription. Stochastic inflationary NR entirely removes these constraints: all gradient terms are retained, and all fields, including metric inhomogeneities, are evolved nonlinearly in real space as soon as modes cross into the IR. This framework enables:

• Direct simulation of nonperturbative gravitational collapse, such as primordial black hole formation,
• Computation of the full, non-Gaussian probability density function of $\mathcal{R}$,
• Real-time evolution of scalar-induced gravitational waves (via the stochastic driving of the anisotropic $\tilde{A}_{ij}$ sector),
• Reliable initial data for post-inflationary cosmic evolution, including preheating and Boltzmann code pipelines.

This approach is extensible to multifield models and can be systematically generalized to include graviton (tensor) modes, which are stochastically sourced and sourced into the evolution of traceless curvature components (Launay et al., 2024).

7. Outlook and Future Directions

Stochastic inflationary numerical relativity achieves, for the first time, the full unification of stochastic quantum diffusion, nonlinear gravitational backreaction, and spatial gradients in the evolution of inflationary universes. Real-space simulations with realistic power spectra, constraint-preserving initial conditions, and robust BSSN evolution lay the groundwork for precision studies of rare cosmological events, generation of primordial black holes, and calculation of non-Gaussian features emergent from genuinely nonperturbative physics.

A key ongoing development is the extension to continuous stochastic source injection (as opposed to all-at-once initialization), which will more accurately capture late-time quantum diffusion. Incorporation of tensor (graviton) stochasticity and implicit schemes for stiffer potentials are natural generalizations. These advances will enable the first truly general-relativistic, nonperturbative predictions for observables such as primordial black hole abundances, induced gravitational waves, and higher moments of curvature perturbation statistics (Launay et al., 16 Dec 2025, Launay et al., 2024, Launay et al., 10 Feb 2025).