COLA Method in Cosmology Simulations

• COLA method is a hybrid simulation technique that combines analytic 2LPT for large-scale evolution with numerical integration of residuals to efficiently model nonlinear cosmological structures.
• It dramatically reduces computational costs—by 100–1,000 times—while achieving percent-level accuracy in power spectra and halo statistics for galaxy surveys.
• The framework is extendable to modified gravity models, incorporating screening mechanisms to enable rapid, precise simulations for next-generation cosmological analyses.

The term "COLA method" ("COmoving Lagrangian Acceleration") denotes a suite of simulation techniques that enable accurate modeling of nonlinear large-scale structure formation in cosmology at a fraction of the computational cost required by full $N$-body simulations. The approach is particularly pivotal for next-generation galaxy surveys which depend on large ensembles of mock catalogues to extract cosmological constraints from the nonlinear regime, including extensions beyond general relativity such as modified gravity (MG) models ($f(R)$ gravity, nDGP, Horndeski). COLA strikes a balance between the speed of Lagrangian Perturbation Theory (LPT) and the accuracy of $N$-body codes, efficiently capturing both large-scale (linear/mildly non-linear) and small-scale (fully non-linear) structure within a hybrid framework (Fiorini, 2023, Winther et al., 2017, Ding et al., 2023, Izard et al., 2015).

1. Motivating Principles and Methodological Foundation

The core motivation of the COLA method is to address the prohibitive computational expense of full $N$-body simulations for cosmological analyses targeting the nonlinear regime, as required by Stage IV galaxy surveys. Traditional $N$-body codes evolve cold dark matter by integrating the full set of Newtonian equations with thousands of timesteps, achieving high accuracy on all scales but at dramatic computational cost. LPT, especially at second order (2LPT), provides an analytic description that is highly accurate on large scales (where bulk flows dominate), but fails on small scales due to nonlinear mode-coupling and shell-crossing.

The COLA framework exploits this scale separation by decomposing the particle displacement into an analytically tractable component (2LPT) and a residual nonlinear component requiring numerical integration:

$$x(q, t) = q + \Psi_\mathrm{LPT}(q, t) + \delta \Psi(q, t)$$

Here, $q$ is the initial Lagrangian position, $\Psi_\mathrm{LPT}$ the 2LPT displacement (analytically evolved), and $\delta \Psi$ the residual solved numerically. The evolution equations for $\delta \Psi$ are integrated only over the small nonlinear residuals, in a frame comoving with the LPT solution, using a particle-mesh (PM) N-body solver. This allows the use of far fewer timesteps (typically 10–40 rather than thousands), as the dominant large-scale evolution is handled analytically (Fiorini, 2023, Tassev et al., 2013, Izard et al., 2015).

2. Formalism and Computational Implementation

In the COLA scheme, particle positions are advanced according to:

$$\begin{align*} x(q, t) &= q + \Psi_\mathrm{LPT}(q, t) + \delta \Psi(q, t) \ \frac{d^2 \delta \Psi}{dt^2} + 2 H \frac{d\delta\Psi}{dt} &= -\nabla_x \Phi_\mathrm{residual} \end{align*}$$

where $H$ is the Hubble parameter, and $\Phi_\mathrm{residual}$ is the potential subtracting off the LPT contribution. The main numerical routine alternates between "kick" and "drift" steps in the residual frame, adding analytically computed LPT terms as source terms or anti-forces. Time-stepping is uniform in scale factor $a$ or a super-comoving time.

Key choices:

• Time-stepping scheme: Typically 10–40 steps, uniform in scale factor, $z_\mathrm{init}=20-50$ (Izard et al., 2015, Ding et al., 2023).
• LPT order: 2LPT is standard; 1LPT (Zel’dovich) degrades performance at large scales (Fiorini, 2023, Tassev et al., 2013).
• Force grid: PM-mesh resolution of $N_\mathrm{mesh}=2$–$3\times N_\mathrm{part}$ ensures subpercent accuracy in power spectrum to $k\sim1\,h/\mathrm{Mpc}$ (Izard et al., 2015).
• Initial conditions: 2LPT at high $z$ for displacement suppression of transients (Ding et al., 2023).

The computational cost of COLA is typically 100–1,000 times lower than full N-body, while retaining percent-level power spectrum accuracy for $k\lesssim 1\,h/\mathrm{Mpc}$ and recovering halo mass functions to a few percent for halos with $\gtrsim300$ particles (Izard et al., 2015, Ding et al., 2023).

3. Extensions to Modified Gravity and Screening Mechanisms

COLA has been extended to models with scale- and environment-dependent modifications to gravity, notably $f(R)$ and nDGP gravity. Linear theory is modified by scale-dependent growth rates, implemented in 2LPT by making the growth factors $k$-dependent, and nonlinearities are approximated via screening factors (Winther et al., 2017, Fiorini, 2023):

• $f(R)$ (chameleon mechanism): The additional scalar degree of freedom is included via a scale-dependent force equation with an environment-dependent screening efficiency: $$\epsilon_\mathrm{screen} = \min[1, |3 f_R(a)/2\Phi_N|]$$.
• nDGP (Vainshtein mechanism): Screening is handled by smoothing the density field and modulating the fifth-force with an analytic factor, imposing the correct limits at both large (linear) and small (nonlinear) scales.
• Practical implementation: Screening factors are calibrated to one full N-body run per model, then used for efficient grid searches (Fiorini, 2023, Winther et al., 2017).

The screening approximation enables MG models to be simulated within COLA with only modest ($2$–$3\times$) overhead relative to COLA–GR, delivering percent-level accuracy in clustering and halo statistics up to $k\sim1$–$3\,h/\mathrm{Mpc}$ (Winther et al., 2017).

4. Performance, Accuracy, and Validation Across Applications

COLA’s accuracy has been validated across a diverse array of cosmological observables and structure formation statistics:

• Matter and velocity power spectra: Errors of 1–2% up to $k=1$–$3\,h/\mathrm{Mpc}$ for both $\Lambda$CDM and MG models (Fiorini, 2023, Winther et al., 2017, Izard et al., 2015).
• Two- and three-point statistics: COLA recovers the monopole and multipoles of the two-point correlation function within $2$–$4\%$ of $N$-body for $k\lesssim0.4\,h/\mathrm{Mpc}$, and the three-point function (bispectrum) to $\sim$10% in $f(R)$; errors in the latter mainly stem from the screening approximation not capturing complex nonlinear interactions, but are typically washed out by galaxy bias (Fiorini, 2023).
• Void and void-galaxy cross-correlation functions: Real-space and redshift-space statistics of cosmic voids and their cross-correlation with galaxies are reproduced to within 1$\sigma$ error bars (Fiorini, 2023).
• Robustness in emulator and ML contexts: COLA mock catalogues, when combined with subhalo abundance matching (SHAM) and tuned on clustering data, produce percent-level agreement in a variety of clustering and higher-order statistics across cosmological models, consistently demonstrating resilience to moderate cosmology variations (Ding et al., 2023, Gordon et al., 2024).
• Computational speed: Wall-time for a standard resolution COLA run ($1024^3$ particles, $800\,h^{-1}$Mpc box, 30 steps) is $0.4$ h on $256$ cores, compared to $\sim150$ h for a full Gadget-2 run (Ding et al., 2023).

The accuracy of COLA at the relevant scales can be further enhanced by post-processing (halo abundance matching, velocity recalibration) to reach subpercent precision where demanded by survey systematics budgets (Izard et al., 2015).

5. Practical Recommendations and Optimal Settings

The literature converges on best-practice parameter choices for cosmological analyses at percent-level accuracy (Fiorini, 2023, Izard et al., 2015, Ding et al., 2023):

Parameter	Value/Recommendation	Purpose/Effect
Number of steps	30–40, uniform in $a$ from $z_\mathrm{init}\sim 19$–$50$	Ensures stability and accuracy to $k\sim1\,h/\mathrm{Mpc}$
PM grid	$N_\mathrm{mesh}=2$–$3\times N_\mathrm{part}$	Subpercent force accuracy for halo and power spectrum stats
Initial redshift	$z_\mathrm{init}\gtrsim 19$ (2LPT)	Transient suppression, faithful displacement initialization
Screening fudge	Calibrated per MG model	Ensures correct linear-to-nonlinear interpolation
Recalibration	Halo abundance, bias correction, velocity rescaling (if needed)	Post-processing for $<1\%$ halo bias and RSD multipoles

Further extension to massive neutrinos, generalized dark energy backgrounds (e.g., $w_0w_a$ models), and Horndeski/MG theories is possible via appropriate modification of the analytic LPT and force kernels (Fiorini, 2023).

6. Limitations, Extensions, and Frontier Applications

Notable limitations pertain to the scale of applicability: percent-level accuracy is generally maintained only up to $k\sim1\,h/\mathrm{Mpc}$; smaller scales require higher force and time resolution or hybrid schemes. The spherical screening approximation does not capture all nuances of highly nonlinear or chameleon-dominated regimes (e.g., deep voids in $f(R)$); for such cases, further development of nonlinear scalar field solvers or hybrid methods is warranted (Fiorini, 2023, Winther et al., 2017).

Recent literature demonstrates COLA’s value as the nonlinear core for cosmological emulators, permitting rapid parameter-space exploration and enabling Stage-IV inference in extended models with a fraction of the CPU cost of full N-body suites (Fiorini, 2023, Gordon et al., 2024). This positions COLA as a practical engine for cosmological analyses that demand massive mock generation and robust theoretical error control.

7. Summary and Impact

The COLA method is a landmark in computational cosmology, allowing scalable mock generation and high-fidelity nonlinear prediction for galaxy survey science. It fuses analytic 2LPT with efficient PM-residual integration, incorporates screening-modified gravity efficiently, and has been validated for diverse nonlinear probes, including the matter power spectrum, galaxy clustering, void statistics, and emulators for cosmological inference (Fiorini, 2023, Ding et al., 2023, Winther et al., 2017, Izard et al., 2015). COLA's balance of speed, accuracy, and extensibility makes it foundational for both traditional and next-generation analyses in precision cosmology.