Calculating the power spectrum in stochastic inflation by Monte Carlo simulation and least squares curve fitting

This paper proposes a computationally efficient Monte Carlo method for calculating the curvature perturbation power spectrum in stochastic inflation by utilizing least squares curve fitting to avoid nested path generation and obtain a continuous spectrum estimate across a range of scales.

The Gist

The Big Picture: Predicting the Universe's "Fingerprint"

Imagine the early universe as a giant, rolling ball of dough. As it expands (a process called inflation), tiny quantum jitters—like tiny bubbles forming in the dough—get stretched out to become the seeds of galaxies, stars, and everything we see today.

Scientists want to know exactly how "bumpy" this dough was. They use a mathematical tool called the Power Spectrum to measure the size of these bumps at different scales. If the bumps are too big in certain areas, they might collapse into Primordial Black Holes.

The problem? When the universe is in a "diffusion-dominated" state (where the dough is so flat that the ball just jitters randomly instead of rolling smoothly), the math gets incredibly messy. We can't solve it with a simple formula; we have to simulate it.

The Old Way: The "Branching Tree" Nightmare

To simulate this randomness, scientists use a method called Monte Carlo simulation. Think of it like this:

1. The Trunk: You grow one long tree branch representing the history of the universe from the start to the end.
2. The Problem: To measure the "bumpiness" at a specific moment in time, you need to know what happens if the universe had taken a slightly different path at that exact moment.
3. The Old Solution (Nested Simulation):
   • You grow the main trunk.
   • At a specific point, you stop and say, "Okay, let's see what happens if we branch off here."
   • You grow 10,000 tiny branches from that single point to see the average outcome.
   • Then you go back to the main trunk, move to the next point, and grow another 10,000 branches.
   • You do this for every single point you care about.

The Analogy: Imagine you are trying to predict the weather for next Tuesday. The old method says: "Grow a tree for every single second of the day. At every second, cut the tree and grow 10,000 new trees to see if it rains."
The Result: This takes forever. The computer gets overwhelmed because the number of trees grows exponentially. It's like trying to count every grain of sand on a beach by digging a hole for every single grain.

The New Way: The "Smart Guess" (This Paper's Solution)

The authors, Koichi Miyamoto and Yuichiro Tada, came up with a clever shortcut that saves a massive amount of time. They combined two ideas: Smarter Branching and Curve Fitting.

1. The "Two-Path" Trick (No More 10,000 Branches)

Instead of growing 10,000 branches from a single point to find the average, they realized they only need two.

• The Analogy: Imagine you want to know how much a coin flip varies. You don't need to flip it 10,000 times to get a rough idea of the variance. You can flip it twice, see the difference, and use that as a "sample."
• How it works: At any point on the main trunk, they grow just two short branches. They measure the difference between them. Mathematically, this difference gives them a very good estimate of the "bumpiness" (variance) without needing a forest of branches.
• Benefit: This cuts the computational cost by a factor of 10,000 instantly.

2. The "Curve Fitting" Trick (No More Point-by-Point)

The old method calculated the answer for specific points (e.g., Tuesday at 9:00 AM, 9:01 AM, 9:02 AM...). If you wanted the answer for 1,000 different times, you had to run the simulation 1,000 times.

The new method is like drawing a smooth line through scattered dots.

• The Analogy: Instead of measuring the temperature at 1,000 specific spots in a room, you measure it at 10 random spots. Then, you use a computer to draw a smooth curve that fits those 10 points perfectly. Now, you have a formula that tells you the temperature at any spot in the room, not just the 10 you measured.
• How it works: They run their "two-path" simulation at random points in time. They collect the data and use a mathematical technique called Least Squares Curve Fitting to find a smooth function that describes the whole pattern.
• Benefit: They get the answer for the entire range of time at once, rather than calculating it point-by-point.

Why Does This Matter?

1. Speed: The old method was like trying to paint a mural by dipping a single brush into paint for every single pixel. The new method is like using a spray gun that covers the whole wall at once. It makes calculations thousands of times faster.
2. Complexity: This is a game-changer for Multi-Field Inflation. In the early universe, there might have been multiple "inflaton" fields (like multiple balls rolling in the dough) interacting with each other. The old method was too slow to handle this complexity. The new method is fast enough to simulate these complex, multi-field scenarios.
3. Accuracy: By using curve fitting, they avoid the "pixelation" errors that happen when you only look at specific points. They get a smooth, continuous picture of the universe's history.

Summary

The paper introduces a new way to simulate the early universe. Instead of growing a massive, computationally expensive forest of branching paths to measure cosmic fluctuations, they:

1. Grow just two branches at a time to get a quick estimate.
2. Use mathematical curve fitting to turn those scattered estimates into a smooth, complete map of the universe's "bumpiness."

This allows scientists to study complex, chaotic models of the early universe that were previously too difficult to calculate, potentially helping us understand how black holes formed and why the universe looks the way it does today.

Technical Summary

1. Problem Statement

The paper addresses the computational challenges associated with calculating the primordial curvature perturbation power spectrum, $P_\zeta(k)$, within the stochastic-$\delta N$ formalism. This formalism is essential for studying inflationary models where quantum diffusion dominates over classical slow-roll dynamics (e.g., in scenarios producing Primordial Black Holes or in multifield models like hybrid inflation).

Key Challenges Identified:

• Nested Monte Carlo Simulation: The standard numerical approach (Algorithm 1 in the paper) requires a "nested" simulation structure. To calculate the variance of the e-fold number $\langle \delta N^2 \rangle$ at a specific scale $k$, one must:
  1. Generate "trunk" paths from the initial state to the end of inflation (EOI).
  2. Identify the point on each trunk path corresponding to the scale $k$.
  3. Generate a large number of "branch" paths from that specific point to the EOI to estimate the local variance.
• Computational Cost: The total computational cost scales as $O(n_{\text{path},1} \cdot n_{\text{path},2} \cdot n_{PS})$, where $n_{PS}$ is the number of $k$-modes (scales) to be calculated. To achieve statistical accuracy, both $n_{\text{path},1}$ and $n_{\text{path},2}$ must be large (typically $10^4$), leading to a prohibitive cost of $O(1/\epsilon^4)$ for a target error $\epsilon$.
• Discrete Output: Traditional methods only yield $P_\zeta(k)$ at discrete, pre-selected values of $k$, requiring separate simulations for different scales.

2. Methodology

The authors propose a novel algorithm, MCLSFit (Monte Carlo Least Squares Fitting), which eliminates the nested structure and provides a continuous approximation of the power spectrum. The method consists of two main innovations:

A. Unnested Estimator for Variance

Instead of calculating the conditional variance $\langle \delta N^2 \rangle_X$ by generating many branches for every trunk path, the authors utilize an identity derived from the properties of independent random variables:
$$\langle \delta N^2 \rangle_X = \left\langle \frac{1}{2} \left( N^{(1)}_X - N^{(2)}_X \right)^2 \right\rangle$$
where $N^{(1)}_X$ and $N^{(2)}_X$ are the e-fold numbers of two independent paths branching from the same point $X$.

• Implementation: For each trunk path, the algorithm generates only two branches at a randomly chosen backward e-fold $N_{bk}$.
• Benefit: This reduces the computational complexity from $O(n_{\text{path},1} \cdot n_{\text{path},2})$ to $O(n_{\text{path},1})$, as the number of branches is fixed at 2 regardless of the desired statistical precision.

B. Least Squares Curve Fitting

Rather than calculating $P_\zeta(k)$ at discrete points $N_{bk,1}, \dots, N_{bk,n_{PS}}$, the method treats the problem as a regression task:

1. Data Generation: For each of $n_{\text{path}}$ trunk paths, a backward e-fold $N_{bk}$ is sampled uniformly from a target range $[N_{bk,l}, N_{bk,u}]$. Two branches are generated from this point to compute the estimator $Y = \frac{1}{2}(N^{(1)} - N^{(2)})^2$.
2. Regression: The dataset $\{(N_{bk,n}, Y_n)\}$ is fitted to a parametric function $f(N_{bk}, \theta)$ using least squares minimization:
   $$\theta_{\min} = \arg\min_\theta \sum_{n=1}^{n_{\text{path}}} \left( f(N_{bk,n}, \theta) - Y_n \right)^2$$
3. Derivation: The fitted function $\tilde{F}_{\langle \delta N^2 \rangle}(N_{bk}) = f(N_{bk}, \theta_{\min})$ approximates the cumulative variance. The power spectrum is then obtained analytically by differentiation:
   $$\tilde{P}_\zeta(N_{bk}) = \frac{d}{dN_{bk}} \tilde{F}_{\langle \delta N^2 \rangle}(N_{bk})$$
   This avoids finite-difference errors and provides a smooth function over the entire range of scales.

Supporting Method: MCBinAve

The authors also introduce MCBinAve, a binning and averaging technique. It groups the Monte Carlo samples into bins of $N_{bk}$ to produce a rough, non-parametric estimate. This serves as a diagnostic tool to guide the selection of the parametric function $f(N_{bk}, \theta)$ for the main MCLSFit algorithm.

3. Key Contributions

1. Algorithmic Efficiency: The proposed method reduces the computational cost from $O(n_{\text{path}}^2 \cdot n_{PS})$ to $O(n_{\text{path}})$. It eliminates the need for nested simulations and removes the linear dependence on the number of scales ($n_{PS}$).
2. Continuous Spectrum Reconstruction: Instead of discrete data points, the method outputs an analytical approximation function for the power spectrum, allowing for the calculation of $P_\zeta(k)$ at any scale within the fitted range without re-running simulations.
3. Error Estimation: By leveraging the covariance matrix of the fitted parameters (derived from the Jacobian of the fitting function), the method provides rigorous pointwise error estimates for both the variance and the power spectrum.
4. Applicability to Multifield Models: The approach is specifically designed to handle high-dimensional phase spaces (multifield inflation) where traditional PDE-based methods (like the adjoint Fokker-Planck equation) suffer from the "curse of dimensionality."

4. Numerical Results

The authors validated the method across four distinct inflationary models:

• Chaotic Inflation (Single-field, Slow-roll):

  • Result: The MCLSFit result matched the standard Mukhanov-Sasaki (MS) equation-based linear perturbation theory perfectly.
  • Significance: Demonstrated that the method works for standard slow-roll scenarios and correctly reproduces power-law spectra.

• Starobinsky's Linear-Potential Model (Single-field, Ultra-Slow-Roll):

  • Context: Features a sharp transition causing a large peak in $P_\zeta$ followed by oscillations.
  • Result: The method successfully captured the large peak. While the high-frequency oscillations on the baseline were smoothed out by the fitting function, the overall shape and peak magnitude were accurate.
  • Significance: Proved the method's ability to handle non-trivial features like USR-induced amplification.

• Flat Quantum Well (Single-field, Diffusion-dominated):

  • Context: A model where dynamics are purely driven by quantum diffusion. Known for a "leakage" effect where small-scale perturbations affect large scales.
  • Result: The method closely reproduced the analytical solutions from AV20, including the characteristic "leakage" tail for $N_{bk} > \langle N \rangle$.
  • Significance: Validated the method in a regime where quantum diffusion completely dominates, a regime where linear perturbation theory fails.

• Hybrid Inflation (Multifield):

  • Context: A two-field model with a waterfall transition, known for generating enhanced perturbations.
  • Result: The method produced a power spectrum with a single large peak, consistent with results from the adjoint Fokker-Planck approach (Tada-Yamada method).
  • Significance: Demonstrated the method's efficacy in multifield scenarios, where it is most needed due to the intractability of other numerical methods.

5. Significance and Conclusion

This paper presents a paradigm shift in the numerical calculation of stochastic inflation power spectra. By decoupling the estimation of variance from the specific scale $k$ and utilizing least squares fitting, the authors have created a tool that is:

• Computationally Scalable: It makes high-precision calculations for multifield models feasible on standard hardware (the hybrid inflation test took ~10 minutes vs. hours for nested methods).
• Robust: It provides continuous, differentiable outputs with quantified uncertainties.
• Versatile: It works across slow-roll, ultra-slow-roll, and diffusion-dominated regimes.

The authors conclude that while the choice of the fitting function requires some prior knowledge (or guidance from MCBinAve), the MCLSFit method is a powerful, general-purpose tool for exploring complex inflationary scenarios, particularly those involving primordial black hole production and multifield dynamics. They also suggest future extensions to calculate higher-order statistics (bispectra) using similar branching estimators.