Perturbation-theory informed integrators for cosmological simulations

This paper introduces a class of perturbation-theory-informed time-stepping integrators that significantly outperform standard methods like FastPM in cosmological simulations by accurately reproducing density statistics with fewer timesteps, while also demonstrating that convergence is fundamentally limited by shell-crossing effects and that symplecticity is less critical for fast, approximate simulations.

The Gist

Imagine you are trying to predict the future of the entire universe. You want to know how galaxies form, how they cluster together, and how the cosmic web evolves over billions of years. To do this, scientists run massive computer simulations. They take a digital box of "dark matter" particles and let gravity pull them together, step by step, from the beginning of time until today.

The problem? These simulations are incredibly expensive. To get accurate results, you usually need to take millions of tiny, cautious steps. It's like trying to walk across a room by taking steps the size of a grain of sand. It's precise, but it takes forever.

This paper introduces a new way to walk across that room: taking giant, confident strides without losing your balance.

Here is the breakdown of how they did it, using some everyday analogies.

1. The Problem: The "Slow and Steady" Trap

In cosmology, the standard way to simulate the universe is like a hiker checking their map every few inches. This method (called a "symplectic integrator") is very reliable and keeps the energy of the system correct. However, because the universe expands and gravity changes over time, the hiker has to take tiny steps to avoid getting lost. If you want to simulate the whole history of the universe, this takes a massive amount of computer power.

2. The Insight: Knowing the "Script"

The authors realized that for the first part of the universe's history (before things get too chaotic), we actually know the answer. It's like knowing the script of a play before the actors start improvising.

In physics, there is a theory called Lagrangian Perturbation Theory (LPT). Think of this as a "cheat sheet" that tells you exactly where a particle should be if it's just drifting under the influence of gravity and the expansion of space, before particles start crashing into each other.

The standard simulation methods ignore this cheat sheet. They treat every step as a mystery, calculating forces from scratch every time. The authors asked: What if we built the cheat sheet directly into our walking steps?

3. The Solution: "Perturbation-Informed" Integrators

The team created a new family of algorithms (mathematical recipes) that "listen" to the cheat sheet. They call these $\Pi$-integrators.

Here is the analogy:

• Standard Method: You are walking in the dark. You take a step, feel the ground, calculate where you are, take another step, feel the ground again. You are very careful, but very slow.
• The New Method (FastPM, LPTFrog, PowerFrog): You have a GPS that knows the terrain perfectly for the first half of the journey. Instead of feeling the ground, you take a long stride based on where the GPS says you will be. You only check your footing when the GPS gets fuzzy (when particles start crashing into each other).

4. The New Characters: The "Frog" Family

The paper introduces several new "integrators" (walking styles), named after frogs because they jump (drift) and kick (accelerate):

• FastPM: This was an existing method that used the "Level 1" cheat sheet (Zel'dovich approximation). It was good, but the authors wanted to do better.
• LPTFrog: This one uses the "Level 2" cheat sheet. It predicts not just where the particle is going, but how the curve of its path is bending. It's like knowing the road curves ahead, so you steer slightly before you get there.
• TsafPM: A clever twist on the standard method that forces it to obey the cheat sheet rules.
• PowerFrog: The star of the show. This is the "super-integrator." It doesn't just look at the immediate future; it adjusts its steps so that as you get closer to the beginning of time, it perfectly matches the most advanced physics predictions (2LPT). It's the most accurate "giant stride" walker.

5. The Catch: When the Road Gets Bumpy

There is a limit to how far you can stride. In the early universe, particles move smoothly like a flowing river. But eventually, they crash into each other (a moment called "shell-crossing"). The river becomes a chaotic waterfall.

The paper proves a fascinating fact: Once the water turns into a waterfall, no amount of fancy math helps.
If the particles crash into each other, the smoothness of the path breaks. The authors show that even the most advanced "PowerFrog" cannot take giant steps once the chaos begins. At that point, you must go back to taking tiny, careful steps. However, since the "giant strides" cover the smooth, early part of the universe so efficiently, you save a huge amount of time overall.

6. The Results: Speed vs. Accuracy

The authors tested these new methods on massive simulations (using real data from the Quijote and Camels projects).

• The Old Way: Needed 100+ steps to get a decent picture of the universe's structure.
• The New Way (PowerFrog): Needed only 2 to 8 steps to get almost the same accuracy!

It's the difference between taking 100 small steps to cross a room versus taking 8 giant leaps. The new methods reproduce the "power spectrum" (the statistical map of how galaxies are distributed) with incredible accuracy, even with very few steps.

Summary

This paper is about smart shortcuts. Instead of brute-forcing the simulation with millions of tiny steps, the authors built the known physics of the early universe directly into the stepping algorithm.

• Before: "I don't know where I'm going, so I'll take a tiny step and check."
• Now: "I know the script for the first act of the play. I'll take a giant leap, and I'll only slow down when the actors start improvising."

This allows cosmologists to run thousands of simulations in the time it used to take to run one, helping them understand the universe much faster.

Technical Summary

1. Problem Statement

Large-scale cosmological $N$-body simulations are essential for understanding cosmic structure formation but are computationally expensive. To explore cosmological parameter spaces and generate covariance matrices for upcoming surveys (e.g., Euclid, LSST), researchers require fast yet accurate simulations.
The primary bottleneck is the trade-off between computational cost and accuracy. While force calculation algorithms (like Tree-PM or FMM) have been optimized, the number of timesteps required to achieve a given accuracy remains high. Standard symplectic integrators (like the leapfrog/Verlet method) treat the cosmological expansion and gravitational acceleration as generic time-dependent forces, requiring many small steps to reproduce the relatively simple dynamics of the early universe (where gravitational acceleration and Hubble expansion nearly cancel).
The paper addresses the need for time-stepping schemes that incorporate knowledge of cosmological perturbation theory (PT) to reduce the number of timesteps required, particularly in the pre-shell-crossing regime (early times/large scales).

2. Methodology

Theoretical Framework

The authors analyze the Vlasov-Poisson system governing collisionless cold dark matter in an expanding universe. They focus on Drift-Kick-Drift (DKD) integrators, which are standard symplectic schemes.

• Symplecticity & Accuracy: They prove that all canonical DKD integrators are symplectic regardless of the time coordinate. They derive conditions for these integrators to achieve global second-order accuracy, noting that time-dependence in the Hamiltonian introduces specific error terms.
• Shell-Crossing Limit: A key theoretical contribution is Proposition 3, which demonstrates that after "shell-crossing" (when particle trajectories cross and the density field becomes multi-stream), the acceleration field loses regularity. This limits the global order of convergence for any integrator to 3/2 for planar collapse, rendering high-order integrators (e.g., 4th order) futile in the post-shell-crossing regime.

Perturbation-Theory Informed Integrators ($\Pi$-integrators)

The authors introduce a new class of integrators called $\Pi$-integrators.

• Variable Transformation: They re-parameterize the equations of motion using the growth factor $D_+$ as the time variable and define a modified momentum variable $\Pi = d_D \mathbf{X}$ (the rate of change of position with respect to the growth factor).
• Zel'dovich Consistency: They define a property called Zel'dovich consistency, where an integrator must exactly reproduce the analytic Zel'dovich solution (First-Order LPT) in 1D prior to shell-crossing, regardless of the timestep size.
• Construction: They derive a necessary and sufficient condition for the coefficient functions $p$ and $q$ (governing the momentum update) to satisfy this consistency.

New Integrator Schemes

Based on this framework, the authors construct and compare several integrators:

1. FastPM (Baseline): The existing scheme by Feng et al. (2016). The authors prove it is the unique integrator that is both Zel'dovich consistent and symplectic.
2. LPTFrog: A new scheme derived by assuming particles follow a quadratic trajectory (matching Second-Order LPT or 2LPT) within a timestep. It is Zel'dovich consistent but not symplectic (it relates to contact geometry).
3. TsafPM: A scheme constructed by taking the standard symplectic kick factor and modifying the drift factor to satisfy Zel'dovich consistency.
4. PowerFrog: The most advanced integrator. It is constructed to match the 2LPT asymptotic behavior as the scale factor $a \to 0$ (early times). It is designed to capture higher-order corrections better than FastPM or LPTFrog.

3. Key Results

One-Dimensional Tests

• Pre-Shell-Crossing: In 1D, Zel'dovich-consistent integrators (FastPM, LPTFrog, PowerFrog) reproduce the exact analytic solution with zero error (up to round-off) in a single timestep. Standard symplectic integrators require many steps to converge to this solution.
• Post-Shell-Crossing: Once shell-crossing occurs, the advantage of PT-informed integrators diminishes. All integrators converge at the theoretical limit of order 3/2, confirming that higher-order schemes offer no benefit in this regime due to the singularity in the acceleration field.
• Decaying Modes: When initial conditions include decaying modes (which $\Pi$-integrators neglect by design), the PT-informed integrators no longer yield the exact solution but still perform comparably to standard second-order integrators.
• Energy Balance: Non-symplectic integrators (like LPTFrog and PowerFrog) do not exhibit systematic energy drifts larger than symplectic methods in fast simulations, provided timesteps are sufficiently small. The Hubble drag term appears to stabilize the energy evolution.

Two-Dimensional Tests

• In 2D, the Zel'dovich approximation is no longer exact.
• FastPM minimizes error relative to 1LPT (Zel'dovich) but fails to capture 2LPT/3LPT corrections.
• PowerFrog significantly outperforms others, producing displacement fields that are closer to 3LPT than to 2LPT, effectively capturing higher-order physics in a single large timestep.

Three-Dimensional Realistic Simulations (Quijote & Camels Suites)

The authors tested the integrators on realistic $\Lambda$CDM cosmologies using the Quijote (large box, $1 \, h^{-1}\text{Gpc}$) and Camels (small box, $25 \, h^{-1}\text{Mpc}$) simulation suites.

• Performance: In the regime of 1–100 timesteps, the new $\Pi$-integrators (PowerFrog, TsafPM, LPTFrog) significantly outperform standard symplectic integrators and even the FastPM scheme.
• Statistics: They accurately reproduce the power spectrum and equilateral bispectrum of the density field. For example, with only 8 timesteps, PowerFrog achieves cross-power spectrum correlations $>97\%$ on scales up to the Nyquist frequency, whereas standard integrators fail to correlate with the reference solution.
• Convergence: In the small Camels box (highly non-linear), convergence is limited to $\sim 3/2$ due to early shell-crossing, consistent with theory. In the large Quijote box, convergence is closer to order 2, likely due to the smoothing effect of the Particle-Mesh (PM) force resolution masking the singularities.

4. Significance and Impact

1. Efficiency for Simulation-Based Inference: The new integrators allow for "fast approximate simulations" that are orders of magnitude faster than standard methods while maintaining high accuracy. This is critical for simulation-based inference (e.g., differentiable cosmology, FlowPM) where thousands of simulations are needed to estimate cosmological parameters.
2. Theoretical Insight: The paper provides a rigorous proof that symplecticity is not the primary driver of accuracy for fast, approximate cosmological simulations with large timesteps. Instead, matching the physical dynamics (via LPT) is more important.
3. Optimal Integration Strategy: It demonstrates that incorporating perturbation theory directly into the integrator (rather than as a separate correction step like COLA) yields superior performance without the need for hyperparameter tuning.
4. Practical Application: PowerFrog is presented as a drop-in replacement for standard integrators in cosmological codes, offering the best balance of accuracy and speed for large-scale structure formation studies, particularly in the linear and mildly non-linear regimes.

In conclusion, List and Hahn successfully bridge the gap between analytical perturbation theory and numerical integration, providing a new class of integrators that drastically reduce the computational cost of cosmological $N$-body simulations without sacrificing statistical accuracy.