BullFrog Integrator: 2nd-Order PM Scheme

• BullFrog Integrator is a 2nd-order accurate particle-mesh scheme that enforces 2LPT consistency via a novel drift-kick-drift mechanism.
• It updates positions and velocities on-the-fly without storing explicit 2LPT displacement fields, ensuring rapid convergence to the exact Vlasov–Poisson solution.
• Comparative studies show BullFrog achieves target accuracy 2–3× faster than traditional FastPM and leapfrog integrators with no extra computational cost.

The BullFrog integrator is a second-order accurate, multi-step time integration scheme for particle-mesh (PM) N-body simulations in cosmological large-scale structure modeling. It is formulated as a drop-in replacement for standard leapfrog and FastPM integrators, extending their Lagrangian Perturbation Theory (LPT) matching from first (Zel’dovich Approximation, 1LPT) to second order (2LPT). The key innovation is the fully automatic and exact enforcement of 2LPT consistency in both positions and velocities at each time step, without requiring explicit computation or storage of 2LPT displacement fields. BullFrog achieves this by a novel choice of temporal weights in the drift–kick–drift (DKD) scheme, resulting in improved convergence to the exact cosmological Vlasov–Poisson solution at no extra computational cost (Rampf et al., 2024).

1. Theoretical Foundation and Motivation

Conventional leapfrog integrators—using the DKD or KDK splitting—are only second-order accurate in time and require many small steps to reproduce linear motion governed by the Zel’dovich approximation. FastPM improved upon this using D-time stepping and weights $(\alpha, \beta)$ that guarantee first-order LPT (1LPT/ZA) consistency: after each step, trajectories match linear theory up to higher-order corrections. BullFrog generalizes this paradigm to enforce 2LPT matching. Its temporal weights $(\alpha, \beta)$ are chosen such that, for each completed step, both positions and velocities exactly reproduce the solution from truncated 2LPT (before shell crossing) to errors of order $\mathcal{O}(\Delta D^3)$.

BullFrog can be interpreted as a multi-step realization of 2LPT: rather than evolving the gravitational field via a one-shot 2LPT forecast, BullFrog updates the field at each step, re-linearizing and building in nonlinear corrections "on-the-fly." This approach results in rapid convergence even as shell-crossing occurs and in regimes where 2LPT as an expansion breaks down, due to the UV completeness of BullFrog in time. Its global accuracy matches the theoretical optimum for second-order integrators in the growth factor $D$ (Rampf et al., 2024).

2. Second-Order Accurate Time-Stepping Formulas

The integrator employs the linear growth factor $D$ as the time variable. Let $x_n$ and $v_n$ be particle positions and velocities at step $n$. The DKD updates are:

• Half-drift: $x_{n+1/2} = x_n + \frac{\Delta D}{2} v_n$
• Kick: $$v_{n+1} = \alpha_n v_n + \beta_n [D_{n+1/2}^{-1} A(x_{n+1/2})]$$
• Half-drift: $(\alpha, \beta)$0

Here, $(\alpha, \beta)$1 with $(\alpha, \beta)$2, and the force is evaluated along a displaced Zel’dovich trajectory. The DKD weights are uniquely determined by requiring (i) $(\alpha, \beta)$3 for ZA consistency, and (ii) exact matching of the $(\alpha, \beta)$4 coefficients in the velocity update to those from 2LPT at $(\alpha, \beta)$5. These weights can be computed in closed form using the numerically precomputed growth functions $(\alpha, \beta)$6 and $(\alpha, \beta)$7, avoiding the Einstein–de Sitter approximation, which would otherwise spoil second-order convergence in a $(\alpha, \beta)$8CDM universe (Rampf et al., 2024).

3. Algorithmic Structure

The BullFrog algorithm operates as follows:

$x_n$3

Compared to standard leapfrog, modifications are limited to the computation of $(\alpha, \beta)$9 and the use of the D-time variable. No additional force calls, data movement, or storage of LPT grids are introduced. For codes already including FastPM (ZA-consistent) stepping, one replaces those weights with the BullFrog ones for immediate 2LPT consistency (Rampf et al., 2024).

4. Convergence Behavior and Error Analysis

Under BX-specific expansions, the BullFrog weights match those of FastPM to leading and next-to-leading order in $\mathcal{O}(\Delta D^3)$0: $\mathcal{O}(\Delta D^3)$1 Thus, global error in phase space after $\mathcal{O}(\Delta D^3)$2 steps scales as $\mathcal{O}(\Delta D^3)$3. After shell crossing, regularity of $\mathcal{O}(\Delta D^3)$4 is reduced due to caustics, but the error remains controlled ($\mathcal{O}(\Delta D^3)$5 away from caustics) and does not induce secular drift.

Numerical results demonstrate that residual power spectrum errors $\mathcal{O}(\Delta D^3)$6 scale as $\mathcal{O}(\Delta D^3)$7 for large $\mathcal{O}(\Delta D^3)$8 on large scales. For fixed target accuracy, BullFrog converges approximately $\mathcal{O}(\Delta D^3)$9–$D$0 faster than FastPM (Rampf et al., 2024).

5. Comparative Performance and Accuracy

BullFrog achieves superior step efficiency and accuracy compared with both LPT-based approximations and other D-time integrators:

Scenario	BullFrog Steps for <0.25% Error	FastPM Steps for Same Accuracy	COLA Dependence
$D$1, $D$2 and filtered ICs	4BF	50FPM	Requires careful tuning of $D$3
$D$4	10BF	50FPM	10COLA matches 10BF only with best $D$5
$D$6CDM (no cutoff), $D$7/Mpc	10BF: $D$8 at $D$9	64FPM required for 10BF accuracy	COLA with 100 steps matched only when $D$0 and $D$1 are tuned

Standard leapfrog with 2LPT initial conditions converges very slowly, requiring $D$2 steps for $D$3 accuracy. High-order LPT underestimates nonlinear power on small and mid scales. COLA, enforcing explicit 2LPT plus a leapfrog residual, exhibits non-monotonic convergence with step count and requires hyperparameter tuning for each case. In contrast, BullFrog has a monotonic convergence profile and no free hyperparameters save the number and spacing of steps (Rampf et al., 2024).

6. Practical Implementation and Operational Considerations

BullFrog starts from slaved initial conditions at $D$4 (Big Bang), with positions $D$5 and velocities $D$6, eliminating the need for external 2LPT initial condition generators. The growth functions $D$7 and $D$8 are numerically solved once during initialization; using the Einstein–de Sitter approximation for $D$9 is discouraged due to the importance of $x_n$0 corrections for second-order accuracy.

Application of a UV cutoff (filter) to the initial conditions is optional but delays shell crossing and provides a controlled environment for comparisons with LPT. BullFrog’s accuracy persists after shell crossing, indicating robustness in nonlinear regimes.

Computational and memory requirements are identical to standard second-order PM codes: two primary arrays for $x_n$1 and $x_n$2, plus mesh/FFT plans, with costs independent of step number. There are no hyperparameters aside from step count and spacing; accuracy improves monotonically by increasing the number of steps.

BullFrog is positioned as a general-purpose, highly accurate second-order integrator for cosmological N-body forward modeling, with applications in mock catalog generation, cosmological parameter inference, and effective field theory hybrid schemes (Rampf et al., 2024).