Evolving only the observed spacetime region: perfectly parallel lightcone simulations with sCOLA

Introduction

Stage-IV galaxy surveys such as Euclid ([27] (Blanchard et al. 2020)), DESI ([22]), the Rubin Observatory Legacy Survey of Space and Time (Ivezić et al., 2019), Roman (Spergel et al., 2015) and Xuntian (Gong et al., 2026) map the large-scale structure of the Universe over tens of cubic gigaparsecs and over most of cosmic history. Their cosmological constraints increasingly come from the full shape of the clustering statistics rather than from a few compressed summaries (Kobayashi et al., 2022; Adame et al., 2025), which makes them sensitive to the non-linear regime of structure formation, where no exact analytical description exists.

Interpreting these data therefore relies on numerical $N$-body simulations of gravitational collapse ([5]). The mock catalogues built from them calibrate an analysis pipeline end to end, provide the covariance matrices of the summary statistics, whose noise otherwise biases the inferred parameters ([46, 23, 94] (Forero-Sánchez et al. 2025)), and serve as the forward model of simulation-based ([21, 24] (Makinen et al. 2022; Hahn et al. 2024; Elsner et al. 2020; Doeser et al. 2024; Nguyen et al. 2024)) inference. What a survey observes, however, is not the Universe at one epoch but the observer’s past lightcone: each object is seen at the time when the light we receive from it was emitted. A mock catalogue is accordingly assembled in stages — evolving the dark matter, identifying halos (e.g. [11]), assigning galaxies to them by a prescription calibrated on observations (e.g. [100]), and placing those galaxies on the dark matter lightcone (e.g. Monaco et al., 2025) — of which only the first is expensive: the others are comparatively cheap, and it is the cost of the dark matter evolution that decides whether an ensemble of mocks can be produced at all. The present paper is concerned with that step alone: producing the dark matter lightcone.

The scale required for these mock catalogues is set by the surveys themselves. For example, the Galactic north cap covered by the Euclid spectroscopic sample spans the redshift range $0.84 \leq z \leq1.88$ and a comoving volume of tens of $\mathrm{Gpc}^3$ ([28]). A lightcone of that region must cover the same range of epochs, and resolve the halos that host the observed tracers. One simulation of that size is already a major undertaking: state-of-the-art projects consume allocations of order $10^7$ CPU hours ([29] (Monaco et al. 2025; Bayer et al. 2025)). Covariance estimation and simulation-based inference, however, require large ensembles of such simulations. quire not one such realisation but thousands, and at that unit cost an ensemble is out of reach. What is needed is a method that delivers one survey-scale lightcone for a cost of order $10^3$ CPU hours.

The lightcones used today are produced with full $N$-body codes such as PKDGRAV3 [79], Gadget4 (Springel et al., 2021), RAMSES [95], GreeM [54], Abacus (Garrison et al., 2021) and HACC (Heitmann et al., 2019), which have produced the landmark runs of the field: the Horizon 4$\pi$ simulation (Teyssier et al., 2009), the DEUS Full Universe Run (Alimi et al., 2012), MICE-GC (Fosalba et al., 2015a,b; Crocce et al., 2015), the Euclid Flagship simulation [29], Uchuu (Ishiyama et al., 2021), Abacus-Summit (Maksimova et al., 2021) and Farpoint (Frontiere et al., 2022). They reach the accuracy the surveys demand, but at a unit cost that rules out large ensembles: most of these projects deliver a single realisation of the Universe. Such a realisation can be compared with the data, but it cannot on its own furnish the ensemble statistics an analysis needs (a covariance matrix, or the distribution of a summary statistic over realisations) which require many independent realisations of the same volume.

Faster alternatives trade accuracy for speed. Pinocchio [89] (Rizzo et al. 2017), used for the Gepetto mocks (Monaco et al., 2025), predicts halo formation from ellipsoidal collapse instead of solving the equations of motion at all. Where the equations of motion are kept, two distinct strategies reduce the number of force evaluations. The first is to build a time integrator informed by perturbation theory, so that a handful of steps suffices: this is the approach of FastPM (Feng et al., 2016), of the integrators of List & Hahn (2024), and of BullFrog [82]. The second is to change the frame of reference: COLA (COmoving Lagrangian Acceleration, [92]) and COCA (COmoving Computer Acceleration, Bartlett et al., 2025) solve for the residual displacement with respect to an approximation of the real trajectories of particles. In particular, COLA has already been used to produce mocks: cola_halo (Koda et al., 2016) fills the survey volume with remapped and replicated snapshots in redshift bins, while L-PICOLA [52] and ICE-COLA [58, 57] build the lightcone at run time. Machine-learning emulators go further still, replacing the dynamics by a map learnt from the initial conditions (He et al., 2019; Jamieson et al., 2023, 2025; Doser et al., 2024); most of them return snapshots at fixed redshift. Two recent exceptions are a generative model of halo lightcones (Zacharegkas et al., 2026) and a Lagrangian neural cellular automaton that integrates particle trajectories in time, from which a matter lightcone can be recorded on the fly [59]. Whichever approximation is adopted, the small-scale clustering of the resulting lightcone mocks departs from a full $N$-body reference and has to be corrected (Smith et al., 2022).

None of these methods addresses the inefficiency specific to lightcones. A lightcone is usually assembled by stacking replicated copies of the box into concentric shells [63] (Blaizot et al. 2005; Fosalba et al. 2008; Fosalba et al. 2015b), or by recording lightcone crossings on the fly and interpolating between time steps [51] (Evrard et al. 2002; Merson et al. 2013; Smith et al. 2017; Ramakrishnan et al. 2025). Both methods must replicate the box whenever the simulated volume is smaller than the observed one (e.g. [52, 57], as is the case even for the largest runs to date: the Euclid Flagship simulation replicates its $3.6\ h^{-1}$ Gpc box to cover the sky out to $z=3$ [29], HalfDome tiles its $3.75\ h^{-1}$ Gpc box out to $z=4$ (Bayer et al., 2025), and the Gepetto mocks use a $1.2\ h^{-1}$ Gpc box out to $z=2$ (Monaco et al., 2025). The same structures then recur in different directions and at different redshifts, biasing higher-order statistics and small-area surveys notably; reorienting the lightcone inside the box only delays the first repetition [18] (Ramakrishnan et al. 2025). Both approaches also evolve the whole volume to redshift zero, long after distant particles crossed the lightcone: much of the computation is spent on particles already observed. Shrinking the simulation domain as the lightcone contracts recovers part of that waste [72, 39], at the price of reshaping it during the run.

In this paper we use sCOLA (Tassev et al., 2015; Leclercq et al., 2020), an extension of COLA that splits an arbitrarily large volume into small tiles evolved independently (here 216 000 tiles of side $64\ h^{-1}$ Mpc). That independence is what makes the tiling well suited to lightcone construction. A tile can be stopped as soon as the past lightcone has left it, each of its particles having been recorded at its own crossing time, without any reference to the rest of the volume — a restriction in time. And only the tiles intersecting the survey footprint need be evolved at all — a restriction in space. Together, they make the evolved spacetime region coincide with the observed one, up to the tile buffers: no computation is spent on particles that have already been observed, nor on volume that will never be observed. What dedicated lightcone schemes obtain by reshaping the simulated domain during the run is therefore obtained here at no cost, as a by-product of the decomposition. Measured against a full $N$-body reference run with Gadget4 from the same initial conditions, an sCOLA lightcone is about three orders of magnitude cheaper, at the price of approximate dynamics below the scales resolved by the particle mesh. We also introduce two improvements of the sCOLA formalism of Leclercq et al. (2020) at the boundaries of the tiles, where its approximations are least accurate. Together they allow substantially smaller buffers, and therefore a smaller fraction of the computation spent on volume that is afterwards discarded. Section II recalls the sCOLA paradigm, introduces our improvements to it, and details the lightcone-building procedure. Section III describes the simulations we ran. Section IV reports the accuracy and the cost of the method, which we discuss in Section V. We conclude in Section VI; implementation details and parameter studies are provided in the appendices.

Methods

We begin by reviewing the tCOLA and sCOLA methods in § 2 A and II B, respectively. In § 2 C, we discuss the improvements to sCOLA accuracy introduced by this paper, and in § 2 D we describe our lightcone-building procedure.

Review of tCOLA

The temporal COmoving Lagrangian Acceleration (tCOLA, Tassev, Zaldarriaga & Eisenstein [92]) method is a hybrid approach designed to combine the accuracy of $N$-body simulations with the efficiency of analytical approximations. In this framework, the equations of motion for a particle are solved in a reference frame that co-moves with the approximate trajectory of the particle, as predicted by Lagrangian Perturbation Theory (LPT).

In the following, we adopt the notations of Leclercq [66]; Leclercq et al. (2020); Bartlett et al. (2025), where complete derivations can be found. Let $\tau$ denote the conformal time and $a$ denote the scale factor, normalised to unity at the present time, i.e., $a_0 = 1$. The initial and final comoving positions of particles are denoted by $\boldsymbol{q}$ and $\boldsymbol{x}$, respectively. The canonical momenta with respect to the Hubble flow, scaled by $a^2$, are defined as

The Drift coefficient $\mathcal{D}(a)$ is defined as

with $\mathcal{H}(a) \equiv\mathrm{d}a/(a\mathrm{d}\tau)$ being the conformal expansion rate.

Second-order Lagrangian Perturbation Theory (2LPT, Bouchet et al. (1995)) is built upon two potentials. The first-order potential $\phi^{(1)}$ is sourced by the initial overdensity field $\delta_i$, while the second-order potential $\phi^{(2)}$ is sourced by the initial gravitational tidal field:

From these potentials, the time-independent LPT displacement fields $\boldsymbol{\Psi}_1$ and $\boldsymbol{\Psi}_2$ are defined as $\boldsymbol{\Psi}_i \equiv\nabla_q\phi^{(i)}(\boldsymbol{q})$ for $i = 1, 2$. The 2LPT positions and momenta are then given by

where $D_1$ and $D_2$ are the first- and second-order growth factors, respectively.

In the tCOLA framework, the full particle position of any particle is decomposed as (Tassev & Zaldarriaga [91])

where $\boldsymbol{x}_{\mathrm{res}}(a)$ is the residual displacement with respect to the 2LPT trajectory. We also define $\boldsymbol{p}_{\mathrm{res}}$ such that $\mathrm{d}\boldsymbol{x}_{\mathrm{res}}/\mathrm{d}a = \mathcal{D}(a)\boldsymbol{p}_{\mathrm{res}}$. The time evolution of the system is governed by the Drift and Kick equations (Leclercq [66]):

where the Kick coefficient $\mathcal{K}(a)$ is defined by

and a superscript $(0)$ denotes quantities at the present time. The natural variables in tCOLA are $\boldsymbol{x}$ and $\boldsymbol{p}_{\mathrm{res}}$. The $L_-$ and $L_+$ operators perform the change of frame of reference by converting the physical momenta of particles to the residual momenta, and vice versa. They are defined by subtracting or adding the 2LPT momenta to the physical momenta, i.e.

The continuous tCOLA equations (8) and (9) are discretised to allow integration over a finite number of time steps. We use $a$ as the time variable, and choose a symplectic Kick-Drift-Kick (KDK, Quinn et al. (1997)) integrator, which alternates between Kick and Drift steps. In the following, the initial and final scale factors for a Drift step are denoted by $a_i^{D}$ and $a_f^{D}$, respectively, and by $a_i^{K}$ and $a_f^{K}$ for a Kick step. The COLA method usually employs a modified discretisation scheme, introduced by Tassev, Zaldarriaga & Eisenstein [92], involving a free parameter $n_{\mathrm{LPT}}$ (see Leclercq et al. (2020), Appendix B). Following earlier experimentations, we choose $n_{\mathrm{LPT}}=-2.5$, which best reproduces the results of full gravity.

Gravitational forces (the term $\nabla_{x}\left(\Delta_{x}^{-1}\delta\right)$ in the Kick equation (9)) can, in principle, be evaluated using any numerical prescription. So far, implementations of tCOLA have relied on evaluating forces on a particle-mesh (PM) grid, which is the approach adopted here. To do so, we obtain the density field $\delta$ by mass assignment of particles onto a regular mesh using the Cloud-in-Cell (CiC, Hockney & Eastwood, [49]) method, from which the gravitational potential $\Phi= \Delta_{x}^{-1}\delta$ is computed spectrally (see § 2 B), and the force field $\nabla_{x}\Phi$ is derived using finite differences.

Since the residual $x_{\mathrm{res}}$ with respect to the LPT trajectory is typically much smaller than the total displacement, fewer time steps are required to integrate equations (8) and (9) than the standard equations of motion to achieve results at the same accuracy. With respect to tree codes, tCOLA leads to a factor of $\sim100$ reduction in the number of time steps required to reach 1% agreement down to scales of comoving wavenumber $k = 0.5\ h\ \mathrm{Mpc}^{-1}$ ({Tassev, Zaldarriaga & Eisenstein, 2013; Leclercq, Jasche & Wandelt, 2015; Izard, Crocce & Fosalba, 2016; Koda et al., 2016}).

Review of sCOLA

The spatial COmoving Lagrangian Acceleration framework (sCOLA, originally introduced by Tassev et al., 2015) extends COLA by using a comoving Lagrangian frame of reference to spatially decouple far and near fields in space, so that any single small volume embedded inside a larger simulation can reproduce the results of a standard $N$-body run. The idea was later extended and refined by Leclercq et al. (2020) to simulate arbitrarily large volumes by splitting the simulation domain into small, independently evolved tiles. The key challenge is to make each tile self-contained: particles near the border of a tile are gravitationally influenced by particles outside it, and particles should be able to physically move in or out of the tile during the simulation. sCOLA addresses these two issues by introducing two approximations: (i) surrounding each tile with a buffer region that pushes boundary artefacts away from the domain of interest, and (ii) imposing non-periodic boundary conditions on the gravitational potential $\Phi$ to encode the gravitational influence of the full matter distribution on each tile. The result is a simulation scheme that is perfectly parallel and asynchronous: each tile can be evolved independently (Leclercq et al., 2020).

We consider a full simulation volume of comoving size $L$, with $N_{\mathrm{p}}$ particles per axis and $N_{\mathrm{pm}}$ PM grid cells per axis, divided into $N_{\mathrm{t}}^{3}$ cubic tiles, each of size $L^{\mathrm{tile}}$ with $N_{\mathrm{p}}^{\mathrm{tile}} \equiv\frac{L^{\mathrm{tile}}}{L}N_{\mathrm{p}}$ particles per axis (left panel of Figure 1; see also Figure 2 and Table 1 in Leclercq et al. 2020). Each tile is surrounded on all sides by a Lagrangian buffer of initial width $L^{\mathrm{buffer}}$. The Eulerian volume occupied by the tile and its buffer forms the sCOLA box, which is evolved by the $N$-body solver. The generation of the initial conditions in sCOLA is slightly different from the traditional framework; indeed, the generation of the initial positions and momenta of particles in each sCOLA box can be included in the perfectly parallel section of the algorithm by directly computing them from a chunk of the 2LPT potentials $\phi^{(1)}$ and $\phi^{(2)}$ read locally (see Figure 1 of Leclercq et al. 2020). This implies that memory-expensive operations needed on the full box (e.g. fast Fourier transforms) are limited to calculating the LPT potentials.

In a standard periodic PM simulation, the Poisson equation $\Delta_{x}\Phi= \delta$ is solved globally via fast Fourier transforms, which assume periodic boundary conditions on the gravitational potential $\Phi$. In contrast, in sCOLA, each sCOLA box is a finite sub-volume of a larger universe. Therefore, assuming periodicity would incorrectly wrap the gravitational influence of the sCOLA box on itself and be unphysical (see Tassev et al., 2015, for a discussion of the resulting artefacts). Instead, following Leclercq et al. (2020), we impose Dirichlet boundary conditions on the faces of sCOLA boxes, prescribing the value $\Phi_{\mathrm{BCs}}$ of $\Phi$. Solving the Poisson equation $\Delta_{x}\Phi= \delta$ with Dirichlet boundary conditions require the use of discrete sine transforms (DST) with ghost cells immediately outside the active domain of the PM grid (Leclercq et al., 2020, Appendix C).

Boundary conditions $\Phi_{\mathrm{BCs}}$ should be as accurate as possible: errors in $\Phi_{\mathrm{BCs}}$ propagate inward over time, affecting the quality of the gravitational potential field throughout the box, particularly near the sCOLA box edges. Unfortunately, having access to the exact $\Phi_{\mathrm{BCs}}$ would require solving the simulation in the full box, which is what sCOLA seeks to avoid. Thus, to preserve the perfectly parallel nature of the algorithm, in Leclercq et al. (2020), the boundary conditions are approximated under the linearly-evolving potential (LEP) approximation (Brainerd, Scherrer & Villumsen, [15]; Bagla & Padmanabhan [6]):

In this approximation, the spatial structure of the potential is frozen to the initial linear field of $\phi^{(1)} = \Phi_{\mathrm{i}}$ (whose computation is required for the initial conditions), with only its amplitude rescaled by the linear growth factor $D_{1}(a)$. All non-linear effects, such as gravitational clustering, are neglected.

Improvements of sCOLA

In this work, we introduce two elements that improve the accuracy of the sCOLA algorithm, as described by Leclercq et al. (2020): (1) a LPT buffer to approximate the inflow of particles into sCOLA boxes, and (2) second-order boundary conditions on the gravitational potential.

LPT buffer. As the simulation advances in time, particles that were initially outside the sCOLA box should be able to cross its boundary and enter the corresponding Eulerian volume. The standard sCOLA approach considers particles only within the initial sCOLA box and therefore misses this inflow, producing an approximate interior density contrast $\delta_{\mathrm{int}}$ that sources the Poisson equation. This approximation produces artefacts near the boundary of sCOLA boxes and affects modes larger than the sCOLA box in the full simulation. Of the two approximations involved in sCOLA ($\delta\approx\delta_{\mathrm{int}}$ and $\Phi_{\mathrm{BCs}} \approx\Phi_{\mathrm{LEP}}$), Leclercq et al. (2020) Section IV.C and Figure 7 found this one to be the dominant one. To address this shortcoming, we introduce for each tile a LPT buffer, defined as a few additional Lagrangian layers of particles surrounding the sCOLA box (see Figure 1). Within the LPT buffer, particles evolve solely under their 2LPT trajectories, i.e. no gravitational forces are computed for them. They provide an inexpensive but physically motivated proxy for the inflow of matter from outside the sCOLA box and improve the approximation that defines $\delta_{\mathrm{int}}$ (see Appendix A for more details). The LPT displacements $\boldsymbol{\Psi}_1$ and $\boldsymbol{\Psi}_2$ of these particles are generated from the 2LPT potentials $\phi^{(1)}$ and $\phi^{(2)}$ when the sCOLA box is initialised, within the perfectly parallel section of the algorithm.

Second-order boundary conditions. After introducing the LPT buffer, our tests showed that the LEP approximation at the boundary of sCOLA boxes ($\Phi_{\mathrm{BCs}} \approx\Phi_{\mathrm{LEP}}$) became dominant over the approximation of the density contrast ($\delta\approx\delta_{\mathrm{int}}$). For this reason, in this work, we also replace the LEP approximation with an expression derived from second-order standard perturbation theory for the boundary conditions $\Phi_{\mathrm{BCs}}$. The new approximation, derived in Appendix B, reads:

where the second-order correction is:

This expression involves the second-order LPT potential $\phi^{(2)}$ and the squared norm of the first-order displacement $\boldsymbol{\Psi}_1$, both quadratic in $\phi^{(1)}$. Computation of $\Phi_{\mathrm{BCs}}(\boldsymbol{x},a)$ is inexpensive, since $\phi^{(2)}$ and $\boldsymbol{\Psi}_1$ are computed for each particle of the sCOLA box upon initialisation. Unlike the LEP approximation, the new $\Phi_{\mathrm{BCs}}(\boldsymbol{x},a)$ has a time-evolving spatial structure: the second-order term captures some of the effects of non-linear structure formation on the gravitational potential. Taken together, the two improvements described above increase the accuracy of the sCOLA algorithm (see Appendix C for a detailed study). Therefore, to reach a given accuracy, they allow for smaller buffer regions around sCOLA tiles than reported by Leclercq et al. (2020): typically $L_{\mathrm{buffer}} \gtrsim16h^{-1}\,\mathrm{Mpc}$ instead of $\gtrsim40h^{-1}\,\mathrm{Mpc}$, see Appendix D.

Lightcone building

In a simulation, the observer sees its past lightcone, which is defined by the ensemble of null geodesics intersecting both the source and the observer’s worldlines. An exact treatment of lightcone building using a relativistic ray-tracing framework is provided by MهاgrATHÉA-PATHFINDER (Breton et al., 2019; Breton & Reverdy, [16]). Following previous work (e.g. Howlett, Manera & Percival, [52]; Euclid Collaboration et al., [29]), we use a simplified geodesic equation allowing us to construct the lightcone on-the-fly during the simulation, rather than in post-processing. We define the radius of the lightcone from the observer at any past scale factor $a$ as the comoving distance $R_{\mathrm{L}}(a)$ (Hogg, [50]):

with $\mathcal{H}(a)$ the conformal expansion rate. Let us consider a particle at position $\boldsymbol{x}_{p}(a)$, which is at a comoving distance $R_p(a) \equiv\lVert\boldsymbol{x}_{p}(a)-\boldsymbol{x}_{\mathrm{obs}}\rVert$ of the observer located at $\boldsymbol{x}_{\mathrm{obs}}$. It will be seen by the observer when it crosses its past lightcone, i.e. at a time $a_L$ such that

This equation is the lightcone crossing condition. Defining $\Delta R \equiv R_L - R_p$, the crossing condition reads $\Delta R(a_L) = 0$. We assume that $\Delta R$ is a continuous function of $a$. Since particles move slower than light, the particle distance $R_p(a)$ cannot decrease faster than the horizon $R_L(a)$ as $a$ increases, so $\Delta R(a)$ is a decreasing function of $a$ and $\Delta R(1) \le0$ (because $R_L(1) = 0$). For observable particles (particles within the observable radius $R_L(0)$), we have $\Delta R(0) \ge0$. Therefore, for all observable particles, there exists a scale factor $a_L$ such that $\Delta R(a_L) = 0$.

Since particle trajectories are integrated in discrete Drift steps, positions $\boldsymbol{x}$ are only defined at discrete scale factors $a_i^{\mathrm{D}}, a_f^{\mathrm{D}}$. A lightcone crossing occurs between $a_i^{\mathrm{D}}$ and $a_f^{\mathrm{D}}$ when:

Thus, at each time step, for all particles that have not yet crossed the lightcone, we monitor the sign of $\Delta R$. When it becomes non-positive, the particle is flagged, and its position $\boldsymbol{x}_L$ and residual momentum $\boldsymbol{p}_{\mathrm{res},L}$ at $a_L$ are recovered by interpolation. For periodic simulations, the evaluation of $\Delta R(a_f^{\mathrm{D}})$ and interpolation of $\boldsymbol{x}_L$ are done before wrapping particle positions across the box boundaries to avoid unphysical discontinuities in the interpolation.

For a flagged particle, the crossing time $a_L$ and position $\boldsymbol{x}_L$ are obtained by linear interpolation between $a_i^{\mathrm{D}}$ and $a_f^{\mathrm{D}}$: using $\Delta R(a_L)=0$,

The residual momenta $\boldsymbol{p}_{\mathrm{res}}$ are defined at Kick times. Therefore, the residual momentum at lightcone crossing $\boldsymbol{p}_{\mathrm{res},L}$ is interpolated between the two Kick times $a_i^{\mathrm{K}}$ and $a_f^{\mathrm{K}}$ bracketing $a_L$, i.e. such that $a_i^{\mathrm{K}} < a_L \le a_f^{\mathrm{K}}$.

The LPT momentum $\boldsymbol{p}_{\mathrm{LPT}}$ at $a_L$ (from equation (6)) is then added back via the $L_+$ operator to recover the full momentum $\boldsymbol{p}$. Figure 2 summarises the above interpolation procedure.

For a monolithic tCOLA simulation with the observer inside the box, the full simulation is run from the initial time to $a_0=1$, and lightcone outputs $\boldsymbol{x}_L,\boldsymbol{p}_L$ are collected progressively as particles cross the lightcone. In the case of sCOLA, the simulation of each box can be stopped as soon as the lightcone has left the sCOLA box. Indeed, once the lightcone has passed the corner of the box closest to the observer, all particles in the tile have already crossed the lightcone and their outputs $\boldsymbol{x}_L,\boldsymbol{p}_L$ have been collected. Tiles farther from the observer are crossed by the lightcone at earlier times (smaller $a$), so they require fewer time steps, as illustrated in Figure 3. In Figure 4, we show the average fraction of time steps needed as a function of the size of the box. For example, sCOLA lightcone building for a $3.84\,h^{-1}\,\mathrm{Gpc}$ octant with 10 time steps linearly spaced in $D_1$ yields a 57% reduction in the average number of time steps with respect to evolving all of the sCOLA boxes to $a_0=1$.

Gadget4, monolithic tCOLA, and sCOLA handle the periodic wrapping of particles differently when building a lightcone, leading to differences near the edges. Gadget4 replicates the box, which means that there can be duplication of particles in the produced snapshot. Monolithic tCOLA wraps particles during their evolution, so particles can cross the full box boundary and be wrapped before being interpolated. Finally, sCOLA does not wrap particles during evolution but does so at the end when gathering the tiles into a single snapshot. Due to these differences, we remove a safe border of $30\,h^{-1}\,\mathrm{Mpc}$ around the box edges before any comparative analysis.

Simulations

In Section III A, we describe the two simulation volumes used throughout this work and the generation of initial conditions. In Section III B, we present the reference Gadget4 run; in Section III C, the monolithic tCOLA run, which isolates the error due to the sCOLA tiling from that due to the COLA time-stepping; and in Section III D, the sCOLA runs themselves.

General setup

We use two simulation setups, Main and Mini. The Main simulations have a box side length $L = 3.84\,h^{-1}\,\mathrm{Gpc}$, $N_{\mathrm{p}} = 3840$ particles per dimension, and $N_{\mathrm{pm}} = 2N_{\mathrm{p}} = 7680$ PM grid cells per dimension. The particle mass is $m_{\mathrm{p}} = 8.62 \times10^{10}\,h^{-1}\,\mathrm{M}_{\odot}$.

The box is slightly larger than that of the Euclid Flagship simulation (Euclid Collaboration et al. [29]) ($L_{\mathrm{Euclid\ Flagship}} = 3.60\,h^{-1}\,\mathrm{Gpc}$) but the resolution is coarser ($N_{\mathrm{p}}^{\mathrm{Euclid\ Flagship}} = 16\,000$). To compare sCOLA to monolithic tCOLA, we also ran Mini simulations, with the same mass and spatial resolution as the Main simulations but half their box side length: $L^{\mathrm{Mini}} = 1.92\,h^{-1}\,\mathrm{Gpc}$, $N_{\mathrm{p}}^{\mathrm{Mini}} = 1920$, $N_{\mathrm{pm}}^{\mathrm{Mini}} = 3840$.

The initial conditions were generated with the MUSIC2-monofonIC code (Hahn et al. (2020)), modified to write to disk the quantities required by Simbelmynë and Gadget4: the initial overdensity field $\delta_{\mathrm{i}}$, the 2LPT potentials $\phi^{(1)}$ and $\phi^{(2)}$, and the Lagrangian positions $\boldsymbol{q}$. monofonIC relies on CLASS (Lesgourgues [70]) for the initial power spectrum. The initial conditions were generated at scale factor $a_{\mathrm{i}} = 0.05$. The cosmological parameters are those of Planck Collaboration et al. [78]: $\Omega_{\mathrm{m}} = 0.311$, $\Omega_{\mathrm{b}} = 0.049$, $\Omega_{\Lambda} = 0.689$, $H_0 = 67.7\,\mathrm{km\,s^{-1}\,Mpc^{-1}}$, $n_{\mathrm{s}} = 0.967$, and $A_{\mathrm{s}} = 2.10 \times10^{-9}$.

For both tCOLA and sCOLA, we used the KDK leapfrog integrator, with 10 time steps evenly spaced in $D_1(a)$ (following Rampf, List & Hahn [82]). Table I lists the scale factors of the corresponding Kick and Drift operations.

Reference with Gadget4

We used Gadget4 (Springel et al. (2021)) as the reference $N$-body code against which we assess the accuracy and computational cost of sCOLA. It was started from the particles generated by MUSIC2-monofonIC at $a_{\mathrm{i}} = 0.05$. We used the TreePM force computation of Gadget4 with a $3840^3$ PM grid (PMGRID), a force-split scale of 1.25 mesh cells (ASMTH), and a comoving softening length of $50\,h^{-1}\,\mathrm{kpc}$ (1/20 of the mean inter-particle separation); the tree was opened with the relative criterion at ErrTolForceAcc = 0.005, the time-step accuracy parameter was ErrTolIntAccuracy = 0.025, and long- and short-range forces were integrated on a common time step (TREEPM_NOTIMESPLIT). We ran Gadget4 on the Irene cluster of the Très Grand Centre de calcul du CEA (TGCC). The Main simulation was run on 192 nodes (128 CPUs and 228 GB of RAM per node), in many short jobs of about one hour each. In total, it consumed 3 379 000 CPU hours (CPUh). We also ran a Mini reference simulation with Gadget4, using the same parameters except for the PM grid, set to $1920^3$; it consumed 128 200 CPUh.

Monolithic tCOLA

To separate the error due to the sCOLA approximations from that due to the reduced number of COLA

caption: Scale factors of the Kick and Drift operations of the time-stepping used in this article. Each operation runs between two consecutive scale factors (for instance the first Kick is between 0.0500 and 0.0868).

:::

time steps, we ran a monolithic tCOLA simulation with the Simbelmynë code (Leclercq, Jasche & Wandelt, 2015; Leclercq et al., 2020). We used an upgraded version of the code, which includes the sCOLA improvements and the lightcone-building capability described in § 2. Since the code does not support multi-node computation, the number of particles in a monolithic run is limited by the memory of a single node ($1.02 \times10^{10}$ particles for a 2 TB node). We therefore ran tCOLA only in the Mini setup. The tCOLA Mini simulation was run on the Infinity cluster of the Institut d’Astrophysique de Paris (IAP), on a single node with 128 CPUs and 2 TB of RAM, for 1 h 36 min of wall-clock time (205 CPUh). To disentangle the effects of the time integration from those of the force model, we also ran a tCOLA Mini simulation with 100 time steps (1 527 CPUh).

sCOLA

For the Main sCOLA simulations, we used $N_{\mathrm{tiles}} = 60$ tiles per dimension, for a total of 216 000 tiles; and for the Mini sCOLA simulations, $N_{\mathrm{tiles}} = 30$, for a total of 27 000 tiles. Each tile thus had a side length $L^{\mathrm{tile}} = 64 h^{-1}\,\mathrm{Mpc}$ and $N_{\mathrm{p}}^{\mathrm{tile}} = 64$ particles per dimension. We used buffers of width $L^{\mathrm{buffer}} = 16 h^{-1}\,\mathrm{Mpc}$ ($N_{\mathrm{p}}^{\mathrm{buffer}} = 16$). Each sCOLA box is evolved in full, but only the particles of its tile are kept; the number of particles evolved therefore exceeds the number kept by the oversimulation factor (Leclercq et al., 2020)

which equals 3.38 for the configuration used here. Around the buffer, we added an LPT buffer of width $L^{\mathrm{LPT\ buffer}} = 8 h^{-1}\,\mathrm{Mpc}$ ($N_{\mathrm{p}}^{\mathrm{LPT\ buffer}} = 8$), in which particles move only under their 2LPT displacements. Appendix D justifies the choice of these sCOLA parameters.

We ran the Main simulation with sCOLA twice: once up to $a_0 = 1$, and a second time with lightcone-only output, in which case each sCOLA box is stopped as soon as the lightcone has left it (§ 2 D). The sCOLA Main and Mini simulations were run asynchronously on the Infinity cluster of the IAP, on whichever CPUs were available at the time. Each sCOLA box used 4 CPUs and 150 MB of RAM. We also ran sCOLA simulations with 100 time steps for both Mini and Main, and both up to $a_0 = 1$ and lightcone-only.

Results

We now compare our sCOLA lightcone to the reference Gadget4 simulation, and assess the cost of producing it. Section IV A inspects the density and velocity fields directly, where artefacts of the tiling would be most conspicuous. Section IV B turns to summary statistics, which quantify the agreement scale by scale and allow the residual differences to be attributed to individual approximations. Section IV C reports the computational cost of the method, including the savings obtained by stopping each tile as soon as the lightcone has left it.

Density and velocity fields

Figure 5 compares a slice of the Main simulation lightcone density field from sCOLA and Gadget4. We show the quantity $\log_{10}(2 + \delta)$ to enhance the visual contrast between voids and high-density regions. The mean absolute difference in $\log_{10}(2 + \delta)$ between the two fields is 0.036, and no discontinuity is visible at tile boundaries: the sCOLA tiling preserves the continuity of the density field despite the independent evolution of the tiles.

Figure 6 compares a $Y$-$Z$ slice of the line-of-sight velocity field from sCOLA and Gadget4. The velocity field was reconstructed using the Simplex-in-Cell (SiC) projection (Abel, Hahn & Kaehler, [1]; Hahn, Abel & Kaehler, [43]; Hahn, Angulo & Abel, [44]; Leclercq et al., 2017), and it is shown as $v_{\parallel}$ in km s$^{-1}$. sCOLA recovers the large-scale velocity structure of Gadget4, with a mean absolute difference in $v_{\parallel}$ of 9.67 km s$^{-1}$. A few tile-related velocity defects remain, for instance around (800, 480) $h^{-1}\,\mathrm{Mpc}$ in Figure 6. We traced them to particle displacements exceeding the buffer width, and checked that they disappear when the buffer is enlarged to 32 $h^{-1}\,\mathrm{Mpc}$.

Power spectra and bispectra

Figure 7 compares two-point statistics in Fourier space for sCOLA, monolithic tCOLA, and Gadget4, allowing us to separate three sources of discrepancy: (i) the reduced number of time steps (COLA), (ii) the force model (PM versus TreePM), and (iii) the spatial tiling (sCOLA). With 10 time steps, sCOLA reaches 1% agreement with Gadget4 on the matter power spectrum up to $k \simeq0.4 h\,\mathrm{Mpc}^{-1}$, and cross-correlations above 0.999 up to $k \simeq0.6 h\,\mathrm{Mpc}^{-1}$.

The sCOLA power spectrum shows a mild bump with respect to tCOLA near the fundamental mode of the sCOLA box ($k_{\mathrm{sCOLA}}\simeq0.065\,h\,\mathrm{Mpc}^{-1}$, i.e. $2\pi/L_{\mathrm{sCOLA}}$ with $L_{\mathrm{sCOLA}}=L^{\mathrm{tile}}+2L^{\mathrm{buffer}}=96\,h^{-1}\,\mathrm{Mpc}$, the side of the sCOLA box and not of the tile), associated with the finite size of the sCOLA boxes. For the chosen configuration, its amplitude remains below 1%. The bump can be further reduced by changing the sCOLA configuration (see Appendix D): with $L^{\mathrm{tile}}=192\,h^{-1}\,\mathrm{Mpc}$ and $L^{\mathrm{buffer}}=16\,h^{-1}\,\mathrm{Mpc}$, for instance, its amplitude falls below 0.2%.

Increasing the number of time steps from 10 to 100 in sCOLA and tCOLA reduces the overall difference from Gadget4 and pushes the small-scale divergence to smaller scales. For sCOLA, the wavenumber up to which the cross-correlation exceeds 0.999 moves from $k\simeq0.6\,h\,\mathrm{Mpc}^{-1}$ to $k\simeq0.9\,h\,\mathrm{Mpc}^{-1}$. The divergence at small scales is caused by the limited resolution of the PM force model, as opposed to the TreePM method of Gadget4.

We also computed the equilateral bispectrum (see Appendix § 4.7). sCOLA agrees with Gadget4 at the $\sim2\%$ level up to $k\simeq0.4\,h\,\mathrm{Mpc}^{-1}$, confirming that it captures the non-Gaussian features of the density field. The velocity divergence power spectrum agrees with Gadget4 at the 1.5% level up to $k \simeq0.6\ h\ \mathrm{Mpc}^{-1}$(see Appendix F). The power spectra and bispectra were computed on the lightcone, so the Main and Mini simulations probe different redshift ranges, $0 \leq z \leq10.9$ and $0 \leq z \leq1.75$ respectively; this difference in redshift coverage explains the difference between the Main and Mini spectra, notably at small scales.

Computational cost

The Gadget4 Main and Mini simulations cost 3 379 000 CPUh and 128 200 CPUh respectively. The monolithic tCOLA Mini simulation cost 205 CPUh, which corresponds to an average of 9.5 CPU seconds per million particles per force step (hereafter CPU sec / M part / step). Overall, on the Mini volume, tCOLA reduces the computational cost by a factor of 625 with respect to Gadget4. This reduction is due to (i) the reduced number of time steps allowed by COLA and (ii) the use of a PM force model.

For the Main simulation up to $a_0 = 1$, an sCOLA box

Project	Gravity Model	$N_{\mathrm{p}}^3$	$L\ (h^{-1}\,\mathrm{Mpc})$	$m_{\mathrm{p}}\ (h^{-1}\,M_{\odot})$	CPU hours / million $N_{\mathrm{p}}$
Uchuu ((Ishiyama et al., 2021))	GreeM	$12\,800^3$	$2\,000$	$3.27 \times10^8$	$\sim9.54$
Euclid Flagship ([29])	PKDGRAV3	$16\,000^3$	$3\,600$	$1.00 \times10^9$	$\sim16.6^{a}$
HalfDome ((Bayer et al., 2025))	FastPM	$6\,144^3$	$3\,750$	$1.95 \times10^{10}$	$2.63$
Gepetto ((Monaco et al., 2025))	Pinocchio	$2\,160^3$	$1\,200$	$1.52 \times10^{10}$	$0.020$
This work	Gadget4	$3\,840^3$	$3\,840$	$8.62 \times10^{10}$	$59.68$
This work	sCOLA	$3\,840^3$	$3\,840$	$8.62 \times10^{10}$	$0.022$
This work	tCOLA	$1\,920^3$	$1\,920$	$8.62 \times10^{10}$	$0.029$

ran on average for 11.2 s (43.5 CPU seconds), corresponding to a total computational cost of 2 690 CPUh. This amounts to 15.6 CPU sec / M part / step when only the particles kept, i.e. those of the tile, are counted (CPU sec / M part kept / step), above the monolithic tCOLA value because of the oversimulation factor. When all particles evolved are counted, i.e. including the buffer, the cost is approximately 4.6 CPU sec / M part / step, below the tCOLA value because of the imperfect CPU efficiency of the monolithic run on large grids. The memory footprint of sCOLA was 170 MB / M part.

The cost per million particles evolved per time step is expected to be nearly independent of the sCOLA configuration; taking it to be the ∼4.6 CPU sec / M part / step measured for all of the configurations studied here, the total cost of an sCOLA simulation is proportional to the oversimulation factor (here 3.38). In this work, we chose the smallest tiles that achieve an accuracy of 1% on the power spectrum, $L^{\mathrm{tile}} = 64\,h^{-1}\,\mathrm{Mpc}$, to maximise the number of independent tiles, and hence the degree of parallelism. Larger tiles, for instance $L^{\mathrm{tile}} = 192\,h^{-1}\,\mathrm{Mpc}$ at fixed $L^{\mathrm{buffer}} = 16\,h^{-1}\,\mathrm{Mpc}$, would have given an oversimulation factor of 1.59 and divided the total budget by 2.1 (see Appendix D).

Because each sCOLA box is stopped as soon as the lightcone has left it, the boxes of the lightcone run ran for 5.4 s on average (20.8 CPU seconds), for a total of 1 300 CPUh: a 52% reduction with respect to the run to $a_0 = 1$. This reduction is consistent with the predicted 57% reduction in the number of force steps (Figure 3 and Figure 4). The difference is explained by the overhead due to I/O, which does not scale with the number of steps. Overall, the computational cost of the sCOLA lightcone run is lower than that of our Gadget4 run by a factor of 2 600. Four effects drive this reduction: (i) the reduced number of time steps allowed by COLA, (ii) the use of a PM force model, (iii) the absence of global communications between nodes and (iv) the restriction of the evolution to the observed spacetime region. It is partly offset by the oversimulation factor, which multiplies the number of particles evolved by 3.38. These contributions can be separated using the runs of Section III. On the Mini volume, effects (i), (ii) and (iii) together reduce the cost by a factor of 625 with respect to Gadget4. Restricting the evolution to the observed spacetime region adds a factor 2 690/1 300 = 2.1. The tiling itself costs a net factor of the oversimulation factor of 3.38 being partly compensated by the better CPU efficiency of the small sCOLA boxes (4.6 against 9.5 CPU sec / M part / step). The remaining factor of 3.3 is the higher per-particle cost of Gadget4 on the Main volume than on the Mini (59.7 against 18.1 CPUh per million particles), and is therefore a property of the reference code rather than of sCOLA. Together, $625 \times2.1 \times3.3/1.6 \simeq2\,600$.

tbl:II compares the computational cost of our simulations with that of other large cosmological runs. sCOLA reaches a cost per particle comparable to that of approximate methods such as Pinocchio (0.022 versus 0.020 CPUh per million particles) while retaining full $N$-body dynamics.

Discussion

We first ask what controlling the two sCOLA approximations has bought (§ 4 A), and how a tiling should be chosen in practice (§ 4 B). We then consider what restricting the evolved volume to the observed spacetime region makes possible (§ 4 C).

Accuracy improvements with respect to previous work

Compared to Lelercq et al. (2020), the present work improves the accuracy of sCOLA by better controlling its two approximations: the LPT buffer accounts for the inflow of particles, and the second-order gravitational potential improves the boundary conditions. In previous work, the sCOLA matter power spectrum differed from that of tCOLA at large scales (Figure 6 of Lelercq et al. (2020)); this difference disappears with the LPT buffer and the second-order gravitational potential (Figure 8). The gain also allows smaller tiles and buffers: whereas Lelercq et al. (2020) needed a sCOLA configuration of $L^{\mathrm{tile}} = 125\,h^{-1}\,\mathrm{Mpc}$ and $L^{\mathrm{buffer}} = 48.8\,h^{-1}\,\mathrm{Mpc}$ to reach $\sim1\%$ agreement with tCOLA on the matter power spectrum, a configuration with $L^{\mathrm{tile}} = 64\,h^{-1}\,\mathrm{Mpc}$ and $L^{\mathrm{buffer}} = 16\,h^{-1}\,\mathrm{Mpc}$ now reaches 1% agreement with Gadget4.

Accuracy improvements are directly visible in the density and velocity fields. In previous work, the density and velocity differences between sCOLA and tCOLA were of the order of $\pm0.4$ in $\log_{10}(2+\delta)$ and $\pm100$ km s$^{-1}$, respectively (Figure 3 and Figure 4 of Lelercq et al. (2020)); here, between sCOLA and Gadget4, they are of the order of $\pm0.15$ in $\log_{10}(2+\delta)$ and $\pm50$ km s$^{-1}$ (see the colour bar scales of Figure 5 and Figure 6). Although the tiling artefacts of Lelercq et al. (2020) are no longer visible by eye, these residual differences are still preferentially located close to the tile edges and in high-density regions. As reported in Section IV A, these tile-related defects of the velocity field disappear when the buffer is enlarged; Appendix D shows more generally how the agreement with tCOLA improves with buffer size.

The velocity field also shows low-amplitude, spatially extended differences, which we interpret as a misestimation of bulk flows (of the order of $\pm10$ km s$^{-1}$). Both the density and velocity differences decrease with distance from the observer, as structures at higher redshift are less non-linear.

The residual differences between sCOLA and Gadget4 can be attributed to three approximations:

1. The sCOLA tiling leads to a bump at the scale of the sCOLA box in the density power spectrum, the velocity divergence power spectrum, and the

density bispectrum. Its amplitude is controlled by the tile and buffer sizes.

2. The reduced number of time steps in COLA leads to differences at wavenumbers $0.1 \lesssim k \lesssim1\,h\,\mathrm{Mpc}^{-1}$ in the density power spectrum, the velocity divergence power spectrum, and the density bispectrum. These differences can be mitigated by increasing the number of time steps.

3. The PM force model is responsible for the divergence at wavenumbers $k \gtrsim1\,h\,\mathrm{Mpc}^{-1}$. Removing it requires a more accurate force model.

Computational strategy and scalability

A key quantity in sCOLA is the oversimulation factor of equation (22), which quantifies the extra volume simulated because of the buffers. The computational cost is roughly proportional to this ratio, which must therefore be balanced against accuracy requirements. Given the available RAM, the target mass resolution, and the accuracy requirements, one can use the Pareto fronts of Figure 11 and Figure 12 to find the $L^{\mathrm{tile}}$ and $L^{\mathrm{buffer}}$ that minimise the oversimulation factor. For instance, with 2 GB of RAM per sCOLA box, the maximum number of particles per box is $230^3$; at the mass resolution of this work ($8.62 \times10^{10}\,h^{-1}\,M_{\odot}$) and for a required accuracy of 1%, the configuration that minimises the oversimulation factor is then $L^{\mathrm{tile}} = 192\,h^{-1}\,\mathrm{Mpc}$ and $L^{\mathrm{buffer}} = 16\,h^{-1}\,\mathrm{Mpc}$, with an oversimulation factor of 1.59. With this configuration, the Main lightcone simulation would then cost an estimated 600 CPUh.

The tiling strategy enables asynchronous and perfectly parallel computations (see Leclercq et al. 2020 for a detailed discussion). Unlike traditional large $N$-body simulations, this approach does not require tightly synchronised high-performance computing resources. Instead, it allows a flexible use of distributed computing resources and is resilient to node failure. On the other hand, it does not exploit the fast inter-node communication capabilities typical of large supercomputers. The memory footprint of a single sCOLA box is also small, 150 MB for the configuration used here, so that sCOLA could in principle run on a wide range of hardware, and even be distributed over volunteer computing resources.

Lightcone simulations and applications

Cosmological surveys measure the properties of galaxies over a fraction of the sky and within a range of redshifts. Their footprints do not fill their bounding cube (the smallest cube containing the survey geometry). For instance, the Galactic north cap covered by Euclid spectroscopic data, defined as a truncated spherical cone sector covering the redshift range $0.84 \leq z \leq1.88$ and a Galactic latitude $b > 23^{\circ}$ (Euclid Collaboration et al. [28]), has a volume of $42.9\,(h^{-1}\,\mathrm{Gpc})^3$, only 29% of the volume of its bounding cube (Appendix § 4.9).

The observed spacetime region is therefore further restricted by the survey geometry: of the past lightcone contained in the simulation box, only the part covered by the footprint is of interest. The sCOLA tiling makes it possible to evolve only the tiles that overlap that region, so that the cost of a simulation is reduced accordingly. Evolving only the observed spacetime region thus combines a restriction in time, set by the lightcone, and a restriction in space, set by the survey geometry. This spatial restriction also removes the need to replicate the box: since only the tiles overlapping the footprint are evolved, the bounding cube can be made as large as the survey requires while the evolution is paid for over the observed volume alone. The initial conditions are the exception: the LPT potentials must still be computed over the entire cube (§ 2 B), so the memory footprint of that step is set by the bounding volume and not by the survey.

With an allocation of $\sim10^7$ CPUh, comparable to that of state-of-the-art projects (Euclid Collaboration et al. [29]; Monaco et al. 2025; Bayer et al. 2025), one could run $\sim17\,000$ Main lightcone simulations in the $L^{\mathrm{tile}} = 192\,h^{-1}\,\mathrm{Mpc}$, $L^{\mathrm{buffer}} = 16\,h^{-1}\,\mathrm{Mpc}$ configuration. Extrapolating to the Euclid north cap volume, the same budget would allow for $\sim25\,000$ north cap simulations (Appendix § 4.9) at the mass resolution of this work. On the other hand, improving the mass resolution to resolve the smallest halos containing tracers detectable by the survey would require many more particles (as in Euclid Collaboration et al. [29]) and therefore a much higher cost per simulation.

Combined with a prescription for populating the dark matter field with halos and galaxies, sCOLA therefore opens the possibility of producing large ensembles of mock catalogues. Such ensembles are required to estimate the covariance matrices of two-point and higher-order statistics (Hartlap, Simon & Schneider [46]; Kitaura et al. 2016; Forero-Sánchez et al. 2025), to perform simulation-based inference (SBI, also known as implicit-likelihood inference, ILI) (Hahn et al. 2024; Bairagi, Wandelt & Villescusas-Navarro [7]), to test modified gravity models (Fiorini et al. 2021; Ruan et al. 2022), and to probe relativistic effects such as gravitational redshift (Wojtak, Hansen & Hjorth [101]; Alam et al. 2017; Zhu et al. 2017). All these applications need mocks that are deep as well as numerous: a box shallower than the survey suffers from super-sample covariance (Bayer et al. 2023) and, once replicated to fill the volume, from spurious correlations between its copies (Chen & Yu [18]).

Conclusion

Lightcone simulations at survey scale have so far been affordable one at a time. We have shown that sCOLA lifts that constraint: because the volume is tiled into independent boxes, each box can be stopped as soon as the observer’s past lightcone has left it, and only the boxes intersecting the region of interest need be evolved at all, so that evolving only the observed spacetime region becomes a property of the tiling rather than a scheme added on top of it. We also improved the tiling itself, with an LPT buffer that approximates the inflow of particles into a box and a second-order approximation to the gravitational potential at its boundaries; together they remove the large-scale inaccuracies observed by Leclercq et al. (2020) and allow substantially smaller buffers and tiles, which is what makes the method cheap as well as perfectly parallel.

With respect to Gadget4, sCOLA reaches percent-level accuracy in the matter power spectrum up to $k \simeq0.4\,h\,\mathrm{Mpc}^{-1}$, 2% agreement in the equilateral bispectrum up to the same scale, and 1.5% agreement in the velocity divergence power spectrum up to $k \simeq0.6\,h\,\mathrm{Mpc}^{-1}$, at a computational cost reduced by a factor of 2600. Four effects drive this reduction:

1. Changing the frame of reference: sCOLA is an extension of COLA, in which we solve the equations of motion of particles in the reference frame of their 2LPT trajectories; fewer time steps are then needed to reach a given accuracy than in a standard $N$-body integration (Tassev et al. (2015)).

2. Using a particle-mesh force model: gravitational forces are evaluated on a particle-mesh grid rather than with a tree, which makes each time step cheaper, at the price of accuracy below the scales resolved by the mesh.

3. Avoiding global communications: thanks to appropriate approximations at the tile boundaries (a buffer of simulated particles, an LPT buffer, and Dirichlet boundary conditions on the gravitational potential), each tile can be simulated independently of all the others, with no communication at all: the algorithm is perfectly parallel (Leclercq et al. (2020)).

4. Evolving only the observed spacetime region: we simulate only the tiles that intersect the region of interest, and evolve each of them only until the lightcone has left it (the contribution of the present paper).

The perfect parallelism and low cost of sCOLA make it particularly well suited to generating the large ensembles of simulations required by modern cosmological analyses: for instance, we estimate that $\sim25\,000$ lightcone simulations of the Euclid Galactic north cap at a mass resolution of $8.62 \times10^{10}\,h^{-1}\,\mathrm{M}_{\odot}$ would require a budget of $\sim10^7$ CPUh. Such ensembles are essential for covariance estimation, simulation-based inference, and tests of fundamental physics.

Three directions follow naturally. The halo and galaxy statistics of sCOLA lightcones, in redshift space, remain to be validated against a full $N$-body reference, since it is halos, not the matter field, that mock catalogues are built from; there, ten time steps may not suffice, since even a 20-step COLA integrator, which describes the density field around the most massive clusters accurately, does not yield converged virial masses for them (Stopyra et al. (2024)). The accuracy limit set by the particle-mesh force model, visible at $k \gtrsim1\,h\,\mathrm{Mpc}^{-1}$, could be pushed to smaller scales by a more accurate short-range gravitational solver. Finally, the present framework could be combined with an emulated frame of reference to further reduce its computational cost without sacrificing asymptotic convergence to the physical solution (see Bartlett et al. (2025)). With these in hand, producing the ensembles of galaxy mock catalogues that Stage-IV surveys call for becomes a question of allocation rather than of method.

LPT buffer

Leclercq et al. (2020) quantified the error caused by the sCOLA approximations on the density field $\delta\simeq\delta_{\mathrm{int}}$ and on the gravitational boundary conditions $\Phi_{\mathrm{BCS}} \simeq\Phi_{\mathrm{LEP}}$ (their Fig. 7). The dominant source of error was the approximation on the density field $\delta\simeq\delta_{\mathrm{int}}$. The discrepancy between $\delta_{\mathrm{int}}$ and $\delta$ is mostly due to the lack of particle inflow while the Lagrangian sCOLA box deforms (see Figure 1). This flow imbalance leads to defects in the gravitational potential near the boundaries of the box, which propagate as the simulation progresses. To address this issue, we implemented an additional buffer of particles around the sCOLA box, the LPT buffer, to approximate the missing inflow.

Particles in this LPT buffer follow the 2LPT trajectories of equation (5); no forces are computed for them. We compute the positions of the LPT-buffer particles on the fly, and project their mass onto the PM density grid using the same density assignment scheme as for particles of the sCOLA box (CiC here). Formally, using the notation of Leclercq et al. (2020), where $\mathrm{B}^{\mathrm{sCOLA}}$ denotes the density assignment operator, $\{\boldsymbol{x}(a)\}_{\mathrm{sCOLA}}$ the set of particles of the sCOLA box, and $\{\boldsymbol{x}(a)\}_{\mathrm{LPT\ buffer}}$ that of the LPT buffer, we replaced the prescription (Leclercq et al. (2020), equation (19))

by

Figure 1 shows the overall structure of the new tiling scheme: each sCOLA box is composed of a tile and a buffer of particles (with PM trajectory), surrounded by an LPT buffer.

The LPT buffer must be wide enough to contain the particles that can enter the sCOLA box in the course of the run: its width should therefore exceed the distance a particle travels between the initial time and $a_0 = 1$. That distance can be estimated within perturbation theory: at leading order the displacement is $-D_1(a)\boldsymbol{\Psi}_1$, so that the variance of one of its Cartesian components at $a_0 = 1$ is

where $P_{\mathrm{L}}$ is the linear matter power spectrum at $a_0 = 1$; the second-order term contributes only a small correction to $\sigma_{\Psi}$ (Bouchet et al., 1995). For the cosmology of Section III A, equation (25) gives $\sigma_{\Psi} \simeq5.9\,h^{-1}\,\mathrm{Mpc}$ per Cartesian component. The LPT buffer adopted in this work, $L_{\mathrm{buffer}}^{\mathrm{LPT}} = 8\,h^{-1}\,\mathrm{Mpc}$, that is $N_{\mathrm{p}}^{\mathrm{LPT\,buffer}} = 8$ layers of particles (Section III D), therefore covers about $1.4\sigma_{\Psi}$, and so accounts for the inflow of the large majority of the particles that cross the boundary of the sCOLA box. The residual effect of this approximation on the power spectrum is quantified in Appendix C.

Second-order gravitational potential approximation term 2

We start the derivation from the Vlasov–Poisson system in the single-stream regime

where $\Phi$ absorbs the usual normalisation $4\pi G a^2(\tau)\bar{\rho}(\tau)$.

We consider the case of an Einstein–de Sitter universe. In second-order Eulerian (standard) perturbation theory, the density contrast is $\delta(\mathbf{k},\tau) = D_1(\tau)\delta_1(\mathbf{k}) + D_1^2(\tau)\delta_2(\mathbf{k})$, with

where

We aim at building Eulerian perturbation theory at second order in the gravitational potential, such that

with $\Phi_1$ rather than $\delta_1$ as the fundamental field. The first term corresponds to the linearly-evolving potential (LEP) approximation. The potential $\Phi_1$ and the linear density contrast $\delta_1$ are linked by the Poisson equation:

therefore $\Phi_1$ coincides with $\phi^{(1)}$, the gravitational potential that appears in first-order Lagrangian perturbation theory and obeys (in Lagrangian coordinates) $\Delta_q\phi^{(1)}(q) = \delta_1(q)$. Furthermore, equations (B1) and (B5) require

From equations (B2) and (B6), we obtain the formal expression of $\Phi_2(\mathbf{k})$ as a convolution of $\Phi_1(\mathbf{k})$:

where the kernel for the gravitational potential $K_2$ is related to the kernel for the density contrast $F_2$ by:

with $k^2 \equiv\lVert\mathbf{k}_1+\mathbf{k}_2\rVert^2$. In an Einstein–de Sitter universe, $K_2$ is a purely spatial kernel (as is $F_2$).

We now express $K_2$ in terms of the kernels that correspond to the operators $\phi^{(2)}$ and $\lVert\boldsymbol{\Psi}_1\rVert^2$, which are two quadratic invariants of $\phi^{(1)}$ that appear in LPT:

In Fourier space, in Lagrangian coordinates, we have

so the kernel for $\Delta_q\phi^{(2)}$ is

Since $\left[\Delta\phi^{(2)}\right](k) = -k^2\phi^{(2)}(k)$, the kernel for $\phi^{(2)}$ itself is

Furthermore, $(\boldsymbol{\Psi}_1)_\ell(k) = ik_\ell\phi^{(1)}(k)$, so the kernel for $\lVert\boldsymbol{\Psi}_1\rVert^2 = \boldsymbol{\Psi}_1\cdot\boldsymbol{\Psi}_1$ is

Let us seek constants $A$ and $B$ such that

This equation is successively equivalent to:

using $\frac{k^2}{k_1k_2}=\frac{k_1}{k_2}+\frac{k_2}{k_1}+2\mu$. Matching the independent terms in $1$, $\mu$ and $\mu^2$ gives

which are mutually consistent and yield

The final result is

which corresponds in configuration space to the following expression for the gravitational potential at second order:

This expression is valid in Lagrangian coordinates. The boundary conditions are, however, required at the Eulerian positions $x$ of the faces of the sCOLA box, whereas $\phi^{(2)}$ and $\Psi_1$ are known on the Lagrangian grid. Using the expression above with $x$ in place of $q$ amounts to neglecting the Jacobian of the transformation between the two sets of coordinates, which only affects the result at third order. Indeed, the displacement $x-q$ is of first order at leading order (equation (5)); therefore, any field $F$ that is already quadratic in $\phi^{(1)}$ satisfies

whose first correction is of third order. The same counting holds for the derivative operators, which satisfy $\partial/\partial x_i=\partial/\partial q_i+O(x-q)$ and thus modify the equations defining $\phi^{(1)}$ and $\phi^{(2)}$ one order beyond the order at which they are used. Since third-order terms are in any case discarded in the second-order expansion of $\Phi$, keeping track of the Jacobian would not improve the accuracy of the approximation.

Accuracy gains from the LPT buffer and the second-order boundary conditions

With the LPT buffer added, sCOLA has two approximations with respect to tCOLA: (i) the trajectories of inflow particles are 2LPT trajectories rather than tCOLA

trajectories, and (ii) the boundary conditions on the gravitational potential are those of the approximation $\Phi_{\mathrm{BCS}}(a)$ rather than the exact $\Phi(a)$. We tested these two approximations separately, by replacing the trajectories of inflow particles by tCOLA trajectories, and the boundary conditions by the tCOLA potential. We found that once the LPT buffer is enabled, the dominant source of error is the boundary conditions on the potential; this is why we sought an approximation at second order.

Figure 8 quantifies the gain brought by each improvement on its own and by the two combined, in a simulation of $N_{\mathrm{p}}^3=640^3$ particles in a box of $L=640\,h^{-1}\,\mathrm{Mpc}$. The sCOLA configuration is $L^{\mathrm{tile}}=64\,h^{-1}\,\mathrm{Mpc}$, $L^{\mathrm{buffer}}=16\,h^{-1}\,\mathrm{Mpc}$ and $L^{\mathrm{LPT\ buffer}}=8\,h^{-1}\,\mathrm{Mpc}$. Without either of them, sCOLA is within 3% of monolithic tCOLA on the power spectrum and the cross-correlations diverge at $k\simeq0.2\,h\,\mathrm{Mpc}^{-1}$. The LPT buffer alone brings the power spectrum within 2% but leaves the cross-correlations unchanged, which shows that the missing inflow is not the only source of error. The second-order gravitational potential alone also brings the power spectrum within 2%, and pushes the divergence of the cross-correlations to $k\simeq0.5\,h\,\mathrm{Mpc}^{-1}$. With both, the sCOLA power spectrum converges to that of tCOLA at large scales and stays within 1% of it, and the cross-correlations diverge at $k \simeq0.8\,h\,\mathrm{Mpc}^{-1}$.

sCOLA configurations: tile and buffer optimal choices

The comparisons presented in this appendix are carried out on a dedicated test volume, which we refer to as the Micro simulation: a box of side $L^{\mathrm{Micro}} = 960\,h^{-1}\,\mathrm{Mpc}$ with $N_{\mathrm{p}}^{\mathrm{Micro}} = 960$ particles and $N_{\mathrm{pm}}^{\mathrm{Micro}} = 1\,920$ PM grid cells per dimension, that is, the same mass and force resolutions as the Main and Mini simulations of Section III ($d_{\mathrm{p}} \equiv L/N_{\mathrm{p}} = 1\,h^{-1}\,\mathrm{Mpc}$, $m_{\mathrm{p}} = 8.62 \times10^{10}\,h^{-1}\,\mathrm{M}_{\odot}$). The cosmological parameters, the initial conditions, the starting scale factor $a_{\mathrm{i}} = 0.05$ and the 10 KDK time steps are also those of Section III; the simulation is evolved to $a_0 = 1$ and no lightcone is built. This volume is small enough to be run monolithically, which provides the tCOLA reference $P_{\mathrm{ref}}$ used below, and its side is divisible by all the tile sizes considered here.

To compare sCOLA configurations quantitatively, we define the distance of a configuration to the tCOLA reference as the sum of a distance on the power spectrum and one on the cross-correlations:

with $\{k_i\}$ the wavenumber bins, $k_s \equiv1.0\,h\,\mathrm{Mpc}^{-1}$ the maximum wavenumber considered, $N_k$ the number of modes below $k_s$, $P$ and $P_{\mathrm{ref}}$ the sCOLA and tCOLA power spectra, and $R$ their cross-correlation.

sCOLA has two more parameters than a monolithic tCOLA simulation: the size of the tile $L^{\mathrm{tile}}$ and the size of the buffer $L^{\mathrm{buffer}}$. We keep the particle resolution $d_{\mathrm{p}} \equiv\frac{L}{N_{\mathrm{p}}}$ (or equivalently the mass resolution) and the PM grid resolution $d_{\mathrm{pm}} \equiv\frac{L}{N_{\mathrm{pm}}}$ fixed. When exploring the $(L^{\mathrm{tile}}, L^{\mathrm{buffer}})$ space, four constraints must be satisfied:

1. The number of tiles per axis $N_{\mathrm{tiles}} = \frac{L}{L^{\mathrm{tile}}}$ must be an integer;

2. The number of particles per axis in the tile $N_{\mathrm{p}}^{\mathrm{tile}} = \frac{L^{\mathrm{tile}}}{d_{\mathrm{p}}}$ must be an integer;

3. The number of particles per axis in the buffer $N_{\mathrm{p}}^{\mathrm{buffer}} = \frac{L^{\mathrm{buffer}}}{d_{\mathrm{p}}}$ must be an integer;

4. The sCOLA box PM grid size $N_{\mathrm{pm}}^{\mathrm{sCOLA}} = \frac{L^{\mathrm{tile}}+2L^{\mathrm{buffer}}}{d_{\mathrm{pm}}}$ must have a simple prime factor decomposition for the DST operations to be fast. More precisely, the potential in a sCOLA box is obtained with type-I DSTs, whose Dirichlet boundary planes lie one cell outside the stored grid: the grid holds $N_{\mathrm{pm}}^{\mathrm{sCOLA}}-1$ cells and the logical size of the transform, which governs its cost, is $N_{\mathrm{pm}}^{\mathrm{sCOLA}}$. We therefore require this number to be 7-smooth, $N_{\mathrm{pm}}^{\mathrm{sCOLA}} = 2^{\alpha_2}3^{\alpha_3}5^{\alpha_5}7^{\alpha_7}$, and to be dominated by its factors of 2, in the sense that $\alpha_3 + 2\alpha_5 + 3\alpha_7 < \alpha_2$, each odd prime being weighted by the exponent of the power of 2 closest to it (3 $\simeq2^2$, 5 $\simeq2^2$, 7 $\simeq2^3$). When the requested grid does not satisfy this criterion, Simbelmÿne substitutes for it the smallest admissible size at or above it, at fixed PM cell size $d_{\mathrm{pm}}$, so that the force resolution is left unchanged. The substitution is never downwards: an admissible size below $N_{\mathrm{pm}}^{\mathrm{sCOLA}}$, even when it is the closest one, would leave the PM grid too small to cover the sCOLA particle box, and its outermost particles would then receive no force and deposit only part of their mass.

For $L^{\mathrm{tile}} \in[30.0, 200.0]\,h^{-1}\,\mathrm{Mpc}$, $L^{\mathrm{buffer}} \in[8.0,48.0]\,h^{-1}\,\mathrm{Mpc}$, and the Micro volume defined above ($d_{\mathrm{p}} = 1\,h^{-1}\,\mathrm{Mpc}$, $d_{\mathrm{pm}} = 0.5\,h^{-1}\,\mathrm{Mpc}$), 72 sets of parameters satisfy these constraints. Running all of them is unnecessary: the quantities we are after — the trend of the distance to tCOLA with $L^{\mathrm{tile}}$ and $L^{\mathrm{buffer}}$, and the Pareto fronts of these traces — are smooth in both parameters, so a coarser sampling resolves them just as well. We therefore ran 35 configurations, thinning the $(L^{\mathrm{tile}},L^{\mathrm{buffer}})$ grid while keeping the full extent of both ranges. Figure 9 and Figure 10 show the power spectra and cross-correlations of a subset of these configurations compared to tCOLA; Figure 11 and Figure 12 show, for all of them, the Pareto fronts of the distance to tCOLA against the oversimulation factor and the RAM usage.

We compare the distance $D_{\mathrm{ref}}$ to tCOLA, defined in equation (50), to two cost measures: the oversimulation factor of equation (22), which is equivalent to the CPU seconds per million particles kept per step (CPU sec / M part kept / step); and the RAM required, which is equivalent to the number of particles in the sCOLA box, since both the particle arrays and the PM grids scale with the number of particles. The first equivalence holds if the cost per million particles evolved per step, counting those of the buffer, is the same for every configuration. We did not time the Micro configurations individually, so we take this as an assumption; its scale is set by the only sCOLA timing at our disposal, the $\sim4.6$ CPU sec / M part / step measured on the Main boxes (Section IV C). Plotting each cost measure against the distance metric yields the Pareto fronts of Figure 11 and Figure 12. Minimising the oversimulation ratio (CPU sec / M part kept / step) and minimising the RAM usage at a given accuracy are conflicting objectives: the Pareto optimum for the former is the largest tiles, whereas for the latter it is the smallest tiles. Increasing the buffer size moves along the Pareto front. The effect on the Pareto fronts of the particle resolution $d_{\mathrm{p}}$, fixed here to $1.0\,h^{-1}\,\mathrm{Mpc}$ (or equivalently of the mass resolution), is beyond the scope of this article.

The Pareto fronts of Figure 11 and Figure 12 can be used to choose the most efficient configuration as follows. First, determine the RAM available per sCOLA job (with 4 CPUs per job, the RAM associated with these 4 CPUs, for instance 2 GB). Then, set a target accuracy in terms of the distance metric (for instance, 1% accuracy corresponds approximately to $D_{\mathrm{ref}} = 10^{-2}$). These two constraints delimit an allowed region in the bottom-left corner of the RAM Pareto plot (Figure 12): at the mass resolution of this work, 2 GB allow at most $230^3$ particles per sCOLA box. Finally, among the configurations in this region, choose the one with the smallest oversimulation factor in Figure 11 (the lowest point). Here this is $L^{\mathrm{tile}} = 192\,h^{-1}\,\mathrm{Mpc}$, $L^{\mathrm{buffer}} = 16\,h^{-1}\,\mathrm{Mpc}$, with an oversimulation factor of 1.59, the configuration discussed in § 4. Both Pareto fronts also show where the returns from enlarging the buffer stop: at every tile size, the distance to tCOLA improves sharply up to buffer $\simeq16\,h^{-1}\,\mathrm{Mpc}$ and only marginally beyond, while the oversimulation factor keeps growing. With the LPT buffer and the second-order boundary conditions of § 2 C enabled, buffers of $L^{\mathrm{buffer}} \gtrsim16\,h^{-1}\,\mathrm{Mpc}$ are therefore sufficient. In this work, we focused on the configuration $L^{\mathrm{tile}} = 64\,h^{-1}\,\mathrm{Mpc}$, $L^{\mathrm{buffer}} = 16\,h^{-1}\,\mathrm{Mpc}$ to demonstrate that such small tiles and buffers reach 1% accuracy on the power spectrum; larger tiles would however have been more efficient in terms of computing resources.

Equilateral bispectra

We used a fork of the DensityFieldTools code (Flöss, [33]) to compute the equilateral bispectra of the dark matter density fields of the Mini and Main simulations, using the following estimator:

with $\delta_{\mathrm{shell}}^{\mathrm{bin}}$ the overdensity field restricted to one $k$ bin, computed with a fast Fourier transform and its inverse; $N_{\triangle}^{\mathrm{bin}}$ the number of equilateral triangles in the considered bin; $N$ the number of grid cells per dimension used for the Fourier transforms; $V$ the volume of the simulation box; and $B^{\mathrm{bin}}$ the estimated equilateral bispectrum of this bin.

Figure 13 compares the equilateral bispectra of sCOLA, tCOLA and Gadget4. sCOLA agrees with Gadget4 at the $\sim2\%$ level up to $k \simeq0.4\,h\,\mathrm{Mpc}^{-1}$.

The approximations have effects comparable to those observed on the power spectrum: the sCOLA tiling introduces a slight bump around the fundamental mode of the sCOLA box, $k_{\mathrm{sCOLA}}$; the reduced number of time steps leads to an overall increase of the bispectrum with a divergence at small scales; and the PM force model results in a divergence around $k \sim1\,h\,\mathrm{Mpc}^{-1}$.

Velocity power spectra

We deposit the lightcone velocity field on a grid with the Simplex-in-Cell (SiC) projection (Hahn, Angulo & Abel, [44]; Leclercq et al. (2017)) already used in Section IV A, on the same grid as the density field, and compute its divergence $\theta$ in Fourier space with fast Fourier transforms. Figure 14 compares the velocity divergence power spectra of sCOLA, tCOLA and Gadget4, together with their cross-correlations with the overdensity field $\delta$. The density–divergence cross-correlations $P^{\theta\delta}$ and the divergence power spectrum $P^{\theta\theta}$ of sCOLA agree with Gadget4 at the $1.5\%$ level up to $k \simeq0.6\,h\,\mathrm{Mpc}^{-1}$. For both $P^{\theta\delta}$ and $P^{\theta\theta}$ the effects of the approximations are the same as on the matter power spectrum and bispectrum.

Volume and cost of a Euclid north cap simulation

This appendix derives the two numbers quoted in § 4 C: the volume of the Euclid Galactic north cap relative to its bounding cube, and the cost of simulating that region with sCOLA.

We model the north cap as a truncated spherical cone sector. With the cosmological parameters of Section III A, the redshift range $0.84 \leq z \leq1.88$ translates into the comoving range $R_{\min} \leq R \leq R_{\max}$, with $R_{\min}=2015\,h^{-1}\,\mathrm{Mpc}$ and $R_{\max}=3470\,h^{-1}\,\mathrm{Mpc}$, and the Galactic latitude cut $b>b_{\max}=23^\circ$ into the opening angle $\theta_{\max}=90^\circ-b_{\max}=67^\circ$. Integrating $2\pi r^{2}\sin\theta$ over the sector gives its volume,

Since the sector has an axis of symmetry, it fits most tightly in a cube whose main diagonal carries that axis. Let $\alpha\equiv\arccos(1/\sqrt{3})\simeq54.7^\circ$ be the angle between that diagonal and any edge of the cube. Since $\theta_{\max}>\alpha$, the direction of each edge falls inside the sector, which therefore reaches the distance $R_{\max}$ along it; in the opposite direction it extends only to $-R_{\max}\cos(\theta_{\max}+\alpha)$, on its rim. The side of the bounding cube is the sum of the two,

We find $L_{\mathrm{cube}} = 5295\ h^{-1}\ \mathrm{Mpc}$, so that the cap occupies a fraction $s_{\mathrm{survey}} = V/L_{\mathrm{cube}}^{3} = 29\%$ of its bounding cube. The cost of simulating the cap with sCOLA follows. Only the tiles overlapping the footprint need to be evolved, which saves the factor $s_{\mathrm{survey}}$, and over this redshift range the lightcone reduces the average number of time steps by a further factor $s_{\mathrm{lightcone}} \simeq45\%$, estimated from Figure 4. Writing $N_{\mathrm{cube}} = (L_{\mathrm{cube}}/d_{\mathrm{p}})^{3}$ for the number of particles in the bounding cube at the particle resolution $d_{\mathrm{p}} = 1\ h^{-1}\ \mathrm{Mpc}$, and $c$ for the cost per particle and per time step, the cost of one north cap simulation is

With $n_{\mathrm{steps}} = 10$, the oversimulation factor $r \simeq1.6$ of the $L^{\mathrm{tile}} = 192\ h^{-1}\ \mathrm{Mpc}$, $L^{\mathrm{buffer}} = 16\ h^{-1}\ \mathrm{Mpc}$ configuration (equation (22)), and $c = 1.28 \times10^{-9}$ CPUh, which is the 4.6 CPU seconds per million particles per time step measured in Section IV C, this gives $C \simeq400$ CPUh.

Statement of contribution

MA implemented the lightcone construction, the LPT buffer and the second-order boundary conditions in the sCOLA module of Simbelmynë, ran the simulations, performed the analysis, produced the figures, and wrote the first draft of the paper. FL conceived and designed the study, wrote the Simbelmynë code and the original sCOLA algorithm, did the analytic calculation of the second-order gravitational potential approximation term $\Phi_{2}$, contributed to the analysis and interpretation of the results, edited the manuscript, supervised MA, and secured funding. GL contributed to the design of the algorithm, in particular the LPT buffer, co-supervised MA, contributed to the interpretation of the results, contributed to writing the paper, and secured funding. All authors read and approved the final manuscript.

Declaration of generative AI use

Generative AI assistants were used to help develop and review parts of the Simbelmynë code, and to improve the clarity, style, and grammar of the manuscript. All AI-assisted code was reviewed by the authors and subjected to the same testing and validation procedures as manually written code. The authors take responsibility for the scientific content of the paper.

Acknowledgements

We thank Oliver Hahn, Andrew Jaffe, Cornelius Rampf and Benjamin Wandelt for useful discussions, and Ludvig Doeser for comments on the draft. This work was done within the Aquila Consortium (https://www.aquila-consortium.org/). This work has made use of the Infinity Cluster hosted by the Institut d’Astrophysique de Paris; we appreciate the efforts of Stéphane Rouberol for running the cluster smoothly. Infinity was supported by the Domaine d’Intérêt Majeur (DIM) ORIGINES of the Région Île-de-France. This work was granted access to the resources of the Très Grand Centre de calcul du CEA (TGCC) under the allocation A0190416910 made by GENCI. MA is funded by an AMX doctoral fellowship from the École polytechnique. MA, FL and GL acknowledge financial support from the Agence Nationale de la Recherche (ANR) through the grant INFOCW, under reference ANR-23-CE46-0006-01. FL and GL acknowledge financial support from the Centre National d’Études Spatiales (CNES), and from the Simons Collaboration on “Learning the Universe”.

This work made use of Gadget4 (wwwmpa.mpa-garching.mpg.de/gadget4, (Springel et al., 2021)), MUSIC2-monofonIC (github.com/cosmo-sims/monofonIC, (Hahn et al., 2020)), CLASS (class-code.net, [70]) and Simbelmynë (git.aquila-consortium.org/Aquila-Consortium/simbelmyne, (Leclercq, Jasche & Wandelt, 2015)), as well as NumPy ([98]), SciPy ((Virtanen et al., 2020)) and Matplotlib ([53]). The figures use the scientific colour maps of Crameri ((2021)) and the CMasher package ([97]).