COmoving Computer Acceleration (COCA): $N$-body simulations in an emulated frame of reference
Deaglan J. Bartlett Marco Chiarenza Ludvig Doeser Florent Leclercq




















deaglan.bartlett@iap.fr
ORCID: 0000-0001-9426-7723 https://orcid.org/0000-0001-9426-77230000-0001-9426-7723; Corresponding author.









florent.leclercq@iap.fr
https://www.florent-leclercq.eu/
ORCID: 0000-0002-9339-1404 https://orcid.org/0000-0002-9339-14040000-0002-9339-1404; Corresponding author.











§ ABSTRACT 

  
   N-body simulations are computationally expensive, so machine-learning-based emulation techniques have emerged as a way to increase their speed.
   Although fast, surrogate models have limited trustworthiness due to potentially substantial emulation errors that current approaches cannot correct for.
  
   To alleviate this problem, we introduce COmoving Computer Acceleration (COCA), a hybrid framework interfacing machine learning with an N-body simulator.
   The correct physical equations of motion are solved in an emulated frame of reference, so that any emulation error is corrected by design.
   This approach corresponds to solving for the perturbation of particle trajectories around the machine-learnt solution, which is computationally cheaper than obtaining the full solution, yet is guaranteed to converge to the truth as one increases the number of force evaluations.
  
   Although applicable to any machine learning algorithm and N-body simulator, this approach is assessed in the particular case of particle-mesh cosmological simulations in a frame of reference predicted by a convolutional neural network, where the time dependence is encoded as an additional input parameter to the network.
  
   We find that COCA efficiently reduces emulation errors in particle trajectories, requiring far fewer force evaluations than running the corresponding simulation without machine learning.
   As a consequence, we obtain accurate final density and velocity fields for a reduced computational budget.
   We demonstrate that this method shows robustness when applied to examples outside the range of the training data.
   When compared to the direct emulation of the Lagrangian displacement field using the same training resources, COCA's ability to correct emulation errors results in more accurate predictions.
  
   Therefore, COCA makes N-body simulations cheaper by skipping unnecessary force evaluations, while still solving the correct equations of motion and correcting for emulation errors made by machine learning.



COmoving Computer Acceleration (COCA): N-body simulations in an emulated frame of reference
    Florent Leclercq
    September 9, 2024
===========================================================================================









§ INTRODUCTION


N-body simulations represent the state-of-the-art numerical method for studying the dynamics of complex systems, including non-linear gravitational structure formation in the Universe <cit.>.
Such simulations can be incredibly computationally expensive to run <cit.>; hence, various machine learning (ML)-based approaches have been proposed to either remove the requirement to run physical simulations or to reduce the complexity of the simulator used.

The most straightforward application of ML methods is as surrogate models, which take initial conditions as inputs and emulate various features of the corresponding full N-body simulation as outputs.
For example, <cit.> were able to predict halo properties from the initial conditions given certain environmental properties.
<cit.> used Generative Adversarial Networks (GANs) to generate visually plausible three-dimensional cosmological density fields but encountered difficulties in reproducing the correct statistical distribution of physical density fields.
Going further, <cit.> built an emulator of cosmological density fields based on a combination of dimensionality reduction via principal component analysis and supervised ML.
<cit.> demonstrated that one can replicate the full result of particle-mesh (PM) simulations (i.e., a Lagrangian displacement field) using a deep neural network which takes the Zel'dovich approximation <cit.> displacement field as input.
This work was extended to tree-based N-body simulations by <cit.>.
Such emulators can replicate the power spectrum to the percent level up to k ≈ 1   h   Mpc^-1.
<cit.> further extended these works by predicting both the displacement and velocity fields through two separate neural networks and by incorporating the cosmological matter density information through the addition of a “style” parameter <cit.>.
The resulting emulator can reproduce power spectra and bispectra to within a few percent and achieves a similar level of cross-correlation with the true simulation run with the same initial conditions.
By adding a time variable as an additional style parameter, <cit.> were able to eliminate the need for two separate networks and produce an emulator capable of predicting N-body outputs as a function of redshift.
The speed and differentiable nature of particle-based emulators enable them to be integrated within field-based inference of initial conditions <cit.>, with posterior re-simulations indicating faithful reconstruction of the initial conditions.
Simulations including the effects of massive neutrinos <cit.> or modified gravity <cit.> can also be emulated using similar neural network techniques.

Instead of completely bypassing the N-body simulation, one can include ML corrections that capture unresolved physics in low-resolution, cheaper simulations.
For example, <cit.> introduced an additional effective force to PM simulations to capture unresolved forces between particles.
Their machine-learnt isotropic Fourier filter was extended by <cit.> to depend not only on time and wavenumber but also on cosmological parameters.
Super-resolution techniques based on GANs <cit.> and U-nets <cit.> have also been proposed, achieving power spectra correct within a few percent, as well as reasonable bispectra, void size functions, and halo mass functions correct to within 10% for halos down to ≈ 10^11  M_⊙ in mass.
Given the computationally demanding nature of hydrodynamical simulations, <cit.> introduced a light model (with only 𝒪(10) learnable parameters) to transform the output of a dark-matter-only simulation to one that resembles the hydrodynamical simulation run with the same initial conditions.
<cit.> also presented a light and interpretable neural network to produce halo catalogues from dark matter density fields.

Accuracy and interpretability are pivotal challenges in the application of machine learning to N-body simulations.
Despite the high reported accuracy of the methods reviewed above on various tests (mainly using summary statistics), none of these models can be expected to perfectly recover the truth.
Are ML-accelerated simulation algorithms sufficiently accurate to be used in real-world applications?
Without a ground-truth model to compare against during actual use (since such algorithms are designed to eliminate the need for it), current approaches have limited means of identifying the emulation error and cannot correct for it.
Since typical simulations usually also involve simplifying assumptions and approximations, perfectly accurate ML-based models may not be required for many purposes.
The question that arises is then that of the interpretability of ML, in order to control the approximation made with respect to a physical simulator.
Unfortunately, many ML algorithms, including (deep) neural networks, lack interpretability.
If machines predict something humans do not understand, how can we check (and trust) the results?

In this paper, we contend that addressing the lack of interpretability of ML is not always necessary to use an emulator of an expensive model while maintaining control over the degree of accuracy.
We elucidate this argument by constructing a framework in which emulation of N-body simulations is made an ML-safe task by physically rectifying emulation inaccuracies. By “ML-safe,” we mean systems that are reliable, robust, and trustworthy by construction.
The key idea is to find a mathematically equivalent form of the system's equations of motion, where we solve for the (not necessarily small) perturbation around the approximate solution provided by ML.
From a physical point of view, in N-body simulations obeying Newtonian dynamics, this is equivalent to solving the equation of motion in an emulated frame of reference.
Since the ML solution is designed to be approximately correct, computing corrections is numerically easier than evolving the full system, thus requiring fewer evaluations of the forces. 
Through the number and the temporal positions of force evaluations, the user controls the trade-off between speed and accuracy, ranging from fully trusting the ML solution by never correcting particle trajectories to correcting for ML emulation errors at any time step of the simulation.
The system has the theoretical guarantee of asymptotically converging to the physical solution as the number of force evaluations increases.

For gravitational N-body simulations of dark matter particles, we introduce the COmoving Computer Acceleration (COCA) approach to running cosmological simulations within an emulated frame of reference.
While traditional emulators aim to translate initial conditions into final particle positions, directly representing the non-linear dark matter distribution, COCA aims to emulate a frame of reference in which to run a physical simulation with lower computational cost.
The approach can be seen as a generalisation and improvement of the idea behind COmoving Lagrangian Acceleration (COLA) <cit.>.
As an illustration, we compare the results of COLA and COCA simulations when forces are evaluated through a particle-mesh (PM) scheme.
We find that using our ML-enhanced approach requires very few force evaluations (approximately 8, compared to 20 for COLA) to correct for emulation errors, yielding percent-level accurate power spectra, bispectra, and cross-correlation to a reference simulation.

This paper is organised as follows.
In <ref>, we review the COLA approach to N-body simulations, extend it to yield COCA, and describe the benefits of COCA in terms of computational efficiency. A more thorough description is provided in <ref>.
We introduce our emulator for the frame of reference in <ref> and describe the training procedure and validation metrics for the COCA simulations. 
In <ref>, we present our results: the performance of the emulator, the accuracy of COCA simulations as a function of the number of force evaluations, the generalisation to an example known to be outside the range of the training set, the comparison to a Lagrangian displacement field emulator, and a discussion of the computational performance.
We conclude in <ref>, discussing potential future extensions and applications of this study.



§ THEORY



For simplicity, some of the equations in this section are abridged. We reintroduce the omitted constants, temporal prefactors, and Hubble expansion in <ref>.



 §.§ Review of COLA


In a cosmological dark matter-only N-body code, one wishes to compute the final Eulerian positions of particles x, as a function of scale factor a, as they interact under gravity. If the initial comoving particle positions are q, then the Lagrangian displacement field is given by <cit.>

    Ψ(q,a) ≡x (a) - q.

One then must solve the equation of motion which reads schematically

    ∂_a^2 Ψ(q,a) = -∇Φ(x,a),

where the gravitational potential Φ satisfies the Poisson equation sourced by the density contrast field δ (x, a),

    ΔΦ(x,a) = δ(x,a).


In the perturbative regime, analytic solutions for Ψ(q,a) can be derived, which are known as Lagrangian Perturbation Theory <cit.>. These solutions are valid on large scales but become inaccurate once shell crossing occurs, making the approximation more reliable at early times.
This behaviour is illustrated in <ref>, where initially the LPT trajectories and the true trajectories are indistinguishable, but the discrepancy increases over time.

The temporal COmoving Lagrangian Acceleration <cit.> algorithm aims to separate the temporal evolution of large and small scales by evolving large scales using analytic LPT results and small scales using a numerical solver. This is accomplished by decomposing the Lagrangian displacement field into two components <cit.>:

    Ψ(q,a) ≡Ψ_LPT(q,a) + Ψ_res^COLA(q,a),

where Ψ_LPT(q,a) represents the LPT displacement field, and Ψ_res^COLA(q,a) denotes the residual displacement of each particle as observed from a frame comoving with an “LPT observer,” whose trajectory is defined by Ψ_LPT(q,a). Knowing Ψ_LPT(q,a), one does not need to solve for the full trajectory of the particle, but just the residual between the approximation and the truth (the blue arrows in <ref>). 

Using <ref>, it is possible to rewrite <ref> as 

    ∂_a^2 Ψ_res^COLA(q,a) = -∇Φ(x,a) - ∂_a^2 Ψ_LPT(q,a).

Therefore, one can view LPT as providing a new frame of reference within which we solve the equations of motion. 
The term ∂_a^2 Ψ_LPT(q,a) can be thought of as a fictitious force acting on particles, caused by our use of a non-inertial frame of reference.

Since particles experience lower typical accelerations in the LPT frame compared to the natural cosmological frame, solving the equation of motion numerically becomes a comparatively simpler task, requiring fewer time steps to achieve equivalent accuracy <cit.>. In particular, COLA has been demonstrated to always yield correct results at large scales, even with a small number of time steps (≤ 10), unlike a basic particle-mesh (PM) code. Given that <ref> is mathematically equivalent to <ref>, COLA is asymptotically equivalent to the corresponding standard N-body code (e.g. a PM code if forces -∇(Δ^-1δ) are evaluated via a standard PM technique), in the limit of an infinite number of time steps.



 §.§ COCA formalism





While the COLA formalism has proven effective in solving the equation of motion within the frame of reference defined by LPT, there is no requirement to adhere to LPT or any other analytic approximation in the decomposition given by <ref>.
According to the principle of Galilean invariance, the equation of motion can be solved in any frame of reference, provided appropriate fictitious forces are introduced for non-inertial frames.
Considering that the simplest scenario is the one where no motion occurs, we aim at finding a frame of reference where particles are nearly stationary. In such a frame of reference, solving the equation of motion numerically to reach a given level of accuracy becomes an easier problem than in COLA. This is the key insight that underpins the formalism proposed in this paper. We dub this approach COmoving Computer Acceleration (COCA).

We propose utilising a ML algorithm, such as a neural network, as an emulator to learn and predict the trajectories of particles in N-body simulations. Since the LPT frame already provides a good approximation of the trajectories, particularly on large scales, we opt to learn displacements relative to the LPT frame. Therefore, the emulator outputs a displacement field Ψ_ML(q,a) that approximates Ψ(q,a) - Ψ_LPT(q,a).

Rather than directly employing the emulator as a surrogate for simulation results, we use the frame of reference corresponding to the emulated trajectories in order to run a simulation. Following the same spirit as COLA, we split the Lagrangian displacement field into three contributions,

    Ψ(q,a) ≡Ψ_LPT(q,a) + Ψ_ML(q,a) + Ψ_res^COCA(q,a),

where Ψ_ML(q,a) is the ML contribution to the Lagrangian displacement field, and the residual displacement Ψ_res^COCA(q,a) represents the emulation error. Different contributions are shown schematically in <ref>.

Reframing <ref> using <ref>, the equation of motion for COCA contains an extra fictitious force with respect to COLA:

    ∂_a^2 Ψ_res^COCA(q,a) =    -∇Φ(x,a) - ∂_a^2 Ψ_LPT(q,a) 
       - ∂_a^2 Ψ_ML(q,a).


In COCA, the predicted displacement Ψ_LPT(q,a) + Ψ_ML(q,a) approximates the optimal frame of reference in which to solve the simulation (the one where all particles are at rest). Ideally, in case of perfect emulation, solving the equation of motion would result in no trajectory adjustment (Ψ_res^COCA(q,a) =0 for any a). Otherwise, numerically solving <ref> corrects the trajectories of particles to produce a more accurate solution.

We describe in more detail the COCA formalism in <ref>. Notably, while above we described the framework in terms of an emulated displacement field Ψ_ML(q,a), we show that we can equivalently define the new frame of reference by the momentum p≡dx/da of particles, so that

    p(a) ≡p_LPT(a) + p_ML(a) + p_res(a),

where p_LPT(a) and p_ML(a) denote momenta predicted by LPT and ML, respectively, and the residual momentum p_res(a) is determined by solving the equations of motion.



 §.§ Reducing the number of force evaluations


To integrate the equations of motion in the new frame of reference, we utilise a symplectic “kick-drift-kick” (leapfrog) algorithm <cit.>. 
With this method, the positions x, and momenta p, of the particles are updated at different times, typically with one momentum update between every two position updates.
A schematic illustration of the technique is given in <ref>, with the full details provided in <ref>.

At each momentum update (“kick”) we face two choices. One can assume that the emulated frame of reference is sufficiently accurate and thus update the particle momenta by simply evaluating the emulator, corresponding to following the “emulated” (purple) trajectory in <ref> (equivalent to assuming that g_δ(t^D) = 0 in equation (<ref>), using the notations of <ref>).
Alternatively, one may deem the emulation error significant and opt to correct the trajectory, aiming to bring the particles back to the “N-body” (black) trajectory in <ref>.
This correction involves evaluating gravitational forces between particles[In a PM scheme, which we employ in this paper, forces are computed by deriving the density field from particle positions through cloud-in-cell binning, solving the Poisson equation in Fourier space to obtain the gravitational potential, and then finite differencing the potential in configuration space to get the forces.] (g_δ(t^D) in <ref>) and using the complete form of the kick operator.
Correcting trajectories is more computationally expensive than simply following the emulated trajectories, so the number of force evaluations n_f, should be as small as possible, but large enough to correct for emulation errors.
During time steps without force evaluations, particles move according to trajectories defined by their respective frames of reference (Ψ_LPT for COLA and Ψ_LPT+Ψ_ML for COCA). Hence, n_f=0 corresponds to the LPT solution in COLA simulations and a purely emulated one in COCA simulations.

In COCA, the ability to reduce the number of force evaluations introduces an additional degree of freedom compared to PM/COLA simulations, where forces are evaluated at every time step. Force evaluations can in principle be placed at any of the time steps, however we find that concentrating all evaluations towards the end of the simulation, when structure formation is non-linear, typically yields the most accurate results.
Up until the first force evaluation, the COCA framework consists in predicting particle positions x and momenta p at specific times, functioning in a similar way as more traditional emulators <cit.>.
In <ref> we show an example of a kick-drift-kick scheme with ten time steps, with three force evaluations at time steps 8, 9 and 10. At all other time steps, momentum updates (kicks) rely solely on the chosen frame of reference.






§ EMULATION



We remind the reader that the field to be emulated, p_ML(q,a), is the residual momentum field with respect to the LPT momentum field, namely

    p(q,a) - p_LPT(q,a),

at any value of a corresponding to a kick time step (see <ref>).



 §.§ Training data



For the application of COCA described in this work, we chose to emulate the frame of reference in a cubic box of length 128   h^-1  Mpc with N^3 = 64^3 dark matter particles, resulting in a final density field at a=1 on a grid with a resolution of Δ x = 2   h^-1  Mpc. This resolution is approximately the same as that used by <cit.>.[
The number of force evaluations needed in COCA to achieve a given accuracy likely depends on the accuracy of the frame of reference emulator and, therefore, on the resolution. We leave the investigation of such effects to future work.
]
Since the focus of this paper is the time evolution of the fields, we adopt fixed cosmological parameters equal to the best-fit values (TT,TE,EE+lowE+lensing+BAO) from Planck 2018 <cit.>: Ω_ b = 0.0490, Ω_ m = 0.3111, h = 0.6766, τ = 0.0561, n_ s = 0.9665, and σ_8 = 0.8102. We assume a flat Universe and a non-evolving equation of state for dark energy.

Although the COCA formalism can be applied to any method of computing the forces between particles (PM, P^3M, tree-based, etc.), for this paper, we chose to work with a PM force solver, utilising a modified version of the publicly available  code[<https://simbelmyne.florent-leclercq.eu/>] <cit.>.
For our simulations, we generated initial conditions at a scale factor a = 0.05 using second-order LPT and solved the equations of motion using COLA with 20 time steps equally spaced in a and a PM grid of size 64^3 <cit.>.
Although we have verified that this initial scale factor and number of time steps are appropriate to give converged results for all k ≤ 1   h   Mpc^-1, the “reference” against which we compare in testing refers to a COLA simulation with the same setup, except with 100 time steps equally spaced in a.
At each time step of the simulations, we output the difference between the computed momentum of the particles p and the LPT momentum p_ LPT, which is the quantity we must emulate.

We produce 100 simulations for training, 50 for validation, and a further 50 for testing.
This is a sufficiently small number of training simulations that re-training with a different resolution or specifications does not require significant computational resources.
While one could potentially achieve higher accuracy for the emulator with more training simulations, the aim of this paper is primarily to demonstrate how to correct for emulation errors rather than to produce the optimal emulator. Therefore, we find 100 training simulations to be sufficient for our purposes.
For each simulation, we use all 20 output snapshots, resulting in a total of 2000 fields for training. In addition, we use 1000 fields for validation and 1000 for testing.



 §.§ Scaling of momenta




In <ref>, we plot a slice of the field p - p_ LPT as a function of scale factor, for one of our test simulations.
From visual inspection, we find that the large-scale spatial structure of the field to be emulated does not change significantly as a function of time, particularly at late times, but its magnitude does.
We therefore choose to rescale the momenta to be emulated by defining

    p_ ML( q, a )
        ≡
        D(a) ℋ(a) ϖ (a) p̃_ ML( q, a ),

where p̃_ ML is defined to have a standard deviation of unity, D(a) is the linear growth factor, ℋ(a) is the conformal Hubble parameter in units of h, and ϖ(a) is a time-dependent function which we wish to approximate.
Our emulator is designed to directly predict p̃_ ML, and thus this scaling has the benefit of standardising the output, since p̃_ ML has zero mean and standard deviation unity.



To find an approximation for ϖ(a), we compute the standard deviation of the 2000 training p - p_ LPT fields and fit these as a function of a using the ESR <cit.> symbolic regression code.
We use a mean squared error loss function and allow functions to be comprised of addition, multiplication, subtraction, division, the power operator, as well as free constants, θ, and the scale factor a.
Upon inspecting the fitted equations, we find that a power law provides a sufficiently simple yet accurate approximation for our purposes:

    ϖ (a) ≈( θ_0 a )^θ_1,

with parameters θ_0 = 1.1415174 and θ_1 = 2.3103984, which yields a root mean squared error of 1.5 × 10^-3.
We compare this fit to the training data in <ref>, from which we see that it accurately reproduces ϖ(a) at all scale factors. Note that any error in this fit can be compensated for by the emulator, so a perfect fit is not required.



 §.§ Neural network architecture



To emulate p̃_ ML, we utilise a U-net/V-net architecture <cit.>, with a similar implementation to <cit.>.
Our model consists of three resolution levels connected in a “V” shape, using two downsampling (by stride-2 2^3 convolutions) and two upsampling (by stride-1/2 2^3 convolutions) layers. At each level, we apply a 3^3 convolution and, as in a V-Net, we apply a 1^3 convolution as a residual connection <cit.> within each block. A batch normalisation is applied after each convolution, which is followed by a leaky ReLU activation function with a negative slope of 0.01. Each layer has 64 channels, except the input (1), output (3), and those after concatenations (128).

For every convolutional layer, we introduce a “style” parameter[
A style parameter in a neural network is an additional input at each layer that encodes dependence on an important feature. In our case, the scale factor a encodes the dependence on clustering at various cosmological times. For more details, we refer the interested reader to Eq. (1)–(3) of <cit.>.
] 
(borrowing the nomenclature from StyleGAN2, ), where each convolutional kernel is multiplied by an array of the same dimension as the layer's input and with values equal to the style parameter.
Since we are producing a time-dependent emulator, we use the scale factor a as our style parameter.
Our network is implemented and trained using a modified version of the [<https://github.com/eelregit/map2map/>] package and  <cit.>.

The input to the emulator is the redshift-zero linear density field (at a=1). This contrasts with <cit.>, <cit.>, <cit.> and <cit.>, who use the three displacement fields predicted by first-order LPT as input. Given that the latter is computed deterministically from the former, both fields contain the same amount of information. However, the linear density field requires three times less memory to store, making our approach more memory efficient.

We also note that we achieve good performance for our emulator with only N^3=64^3 voxels in the input field, whereas <cit.> uses N^3=128^3 for a similar resolution, resulting in a further factor of 8 reduction in memory requirements.
The smaller input size necessitates one fewer resolution layer in our neural network architecture, thus reducing the number of parameters in our model to 2.4×10^6, compared to 3.4×10^6 in <cit.>.
It also requires less padding of the input field: we use periodic padding of 24 voxels, compared to 48 in <cit.>.

Regarding the dependence of the emulation on cosmology, we expect the sensitivity of p - p_LPT to cosmological parameters to be relatively small, since long-range features should be captured in p_ LPT.
Moreover, we choose to use the linear density field as input instead of the white noise field from which it is produced. This way, our emulator of p - p_LPT only depends on Ω_ m, as the equations of motion depend solely on this parameter.
The dependence on all other cosmological parameters is contained in the linear power spectrum, which is used to transform the white noise field into the linear density field. 
Thus, adding only Ω_ m as a second style parameter to the network would be sufficient to capture the dependence of the framework on cosmological parameters. For simplicity, we fix Ω_ m and save this extension for future work.
Omitting Ω_ m as a second style parameter also enables us to test the robustness of the COCA framework in the case of cosmological parameter misspecification, and hence check for ML-safety.
We discuss this aspect in <ref>.

To summarise, our architecture is similar to that of <cit.>, with three main differences:
(i) we use a single channel (linear density) input rather than three channels (LPT displacements or velocities);
(ii) we have three resolution levels instead of four (since we work with N^3=64^3 grids as opposed to N^3=128^3); and
(iii) we include a as a style parameter <cit.> and fix Ω_ m.



 §.§ Training



As our loss function, we choose

    Loss≡log L_1 + log L_2,

where

    L_n ≡∑_q∑_i
            {[( p_ LPT + p_ ML)_i]^n - 
          [( p_ true)_i ]^n
          }^2,

and the sum runs over the Lagrangian coordinates of the particles, q, and the three Cartesian components, i ∈ x, y, z.

This functional form is partially inspired by <cit.>.
The L_1 term matches p_ML to the residual momenta p_true - p_LPT, whereas the L_2 term ensures that the full momentum field (including the LPT contribution) matches p_true.
<cit.> found that terms similar to L_2 are required to correctly predict redshift-space distortions. We leave the investigation of redshift-space distortions in COCA for future work.
Both terms of our loss function use the mean square error between the fields in Lagrangian coordinates. Unlike <cit.>, we do not include any term in Eulerian coordinates. Given the computational and memory requirements to use the displacement fields in Eulerian coordinates, and the good performance already achieved with our choice, we decided to omit such additional terms.

We use the Adam optimiser <cit.> with decoupled weight decay (AdamW) <cit.>, an initial learning rate of 1.5×10^-4, a weight decay coefficient of 8×10^-3, and parameters β_1=0.85, β_2 = 0.994, and ϵ = 3×10^-9.
The learning rate is reduced on a plateau by a factor of 0.35 when the loss does not improve by more than 10^-3 over 50 epochs. After a change in learning rate, we apply a cooldown of 30 epochs before the scheduler resumes normal operation. We use a batch size of 5 and train on a single V100 GPU, which has 32 GB of RAM.
The entire time for generating the training, validation, and test simulations (for which we use 40 Intel Xeon Gold 6230 cores) and training was 120 hours, corresponding to 277 epochs, by which time the training and validation losses had plateaued.



 §.§ Validation metrics



To quantitatively determine the accuracy of the COCA simulations, we compute the dark matter density field δ using a cloud-in-cell estimator <cit.> and the velocity field v using a simplex-in-cell estimator <cit.>. To work with a scalar field rather than a vector field, we compute the divergence of the velocity field in Fourier space, ∇·v.[
The velocity potential is usually of greater physical interest than the divergence of the velocity field. However, in Fourier space, they are related by a factor of 1/k^2, and since we only compare the ratio of auto and cross spectra at a given k, all quantities shown will be identical for both. Thus, we compute only the divergence.
]

For both fields φ∈{δ, ∇·v}, we compute the (auto) power spectrum P_φ (k) defined by

    ⟨φ( k) φ(  k^') ⟩≡(2π)^3 _ D(  k +  k^') P_φ (k),

where _ D is a Dirac delta distribution. 
For all simulations, we compute the ratio of power spectra between the simulation of interest and the reference.
We also compute the cross spectrum P_φ_a φ_b (k) between the test simulation and the reference simulation, defined by

    ⟨φ_a ( k) φ_b (  k^') ⟩≡(2π)^3 _ D(  k +  k^') P_φ_a φ_b (k).

Thus, we obtain the cross-correlation coefficient

    r_φ_a φ_b( k )
        = P_φ_a φ_b (k)/√(P_φ_a (k) P_φ_b (k)) .

One can interpret 1-r^2 as the fraction of the variance in the prediction that is not explained by the reference.
In schemes such as carpool <cit.>, where one combines exact and approximate simulations, 1-r^2 is proportional to the required number of simulations. Hence, improving r^2 can dramatically reduce the required computational resources. 
Just comparing r can hide the importance of improving the cross-correlation: for example, improving r from 0.9 to 0.99—a change of 0.09—corresponds to explaining an additional 17% of the variance at that scale. 
For these reasons, in all figures, we plot r^2 rather than r, since it is more meaningful.
All two-point statistics are computed using .

To assess higher-order statistics, we also compute the bispectrum B (k_1, k_2, k_3) of the density field, defined by

    ⟨δ( k_1 )δ(  k_2 ) δ(  k_3 ) ⟩≡(2π)^3 _ D( ∑_i=1^3 k_i  ) B (k_1, k_2, k_3),

and, to factor out dependence on scale and cosmological parameters, the reduced bispectrum,

    Q (k_1, k_2, k_3)
        ≡B (k_1, k_2, k_3)/P_1 P_2 + P_2 P_3 + P_3 P_1,

with P_i ≡ P_δ(k_i) for i ∈{ 1,2,3 }.
We consider two different configurations in this work, which are designed to be approximately the same as those used in <cit.>. First, we consider a “squeezed” bispectrum, consisting of an isosceles triangle configuration with one small wavenumber, k_ℓ = 9.8 × 10^-2  h  Mpc^-1, and two sides of equal but varying size, k_1=k_2=k_ s. For our second configuration, we fix two of the wavenumbers, k_1 = 0.1   h  Mpc^-1 and k_2 = 1.0   h  Mpc^-1, and vary the angle θ between them.
All bispectrum calculations are performed using  <cit.>.



§ RESULTS





 §.§ Emulator performance




In <ref>, we plot slices of the input δ_linear, output p_ML, target p-p_LPT, and emulation error p_res≡p-p_LPT-p_ML for one of our test simulations, where all fields are evaluated at a=1. 
The target fields are obtained by running COLA simulations with 20 time steps, using initial conditions matching those of the test simulations, and saving the residuals between the calculated and LPT momenta.
As described in <ref>, the input is the linear density field, comprising a single channel, whereas the output prediction is a three-component vector for each Lagrangian grid point. Since we are learning the residuals between the true momentum and the LPT prediction, correlations observed in p-p_LPT  are highly localised, reflecting the accurate capture of large-scale modes by LPT.

Visually, there is a notable correlation between p - p_LPT (second column of <ref>) and p_ ML (third column of <ref>). Leveraging the linear density field and scale factor information, the emulator accurately identifies the spatial structure of the p_ ML field. The small emulation errors indicate its capability to predict magnitudes as well. We observe that the emulation error is signal-dependent, resulting in larger values of p_ res in the regions where |p_ ML| is large. These regions are highly nonlinear and appear as the simulation progresses. It is noteworthy that the emulation errors become particularly visible when visualising the final quantities in <ref>, given their lesser prevalence at earlier times.



To quantify the magnitude and time-dependence of the emulation error, we plot the root mean squared error (RMSE) between the true p - p_ LPT and p_ ML, as predicted by the emulator, as a function of the scale factor in <ref>.
We present the mean and standard deviation of the RMSE across 50 test simulations.
At early times, the trajectories of the particles are well described by perturbation theory. Thus, even though LPT is not a perfect description of the dynamics, the emulator can easily correct for the error, maintaining a relatively constant RMSE of less than 20   km   s^-1 for a < 0.5.
We observe a slight decrease in RMSE between a=0.2 and a=0.4, which is understandable given <ref>: initially, the field exhibits a high degree of small-scale structure, which becomes less significant and approximately constant over time, making p - p_ LPT easier to predict during this period.
Beyond a ≈ 0.5, the small-scale dynamics become highly non-linear, making it more challenging for the emulator to predict the correct frame of reference.
Consequently, we observe the behaviour schematically illustrated in <ref>: the emulation error grows at late times, approximately doubling between a=0.5 and a=1. 
<cit.> and <cit.> found similar issues with predicting virialised motions within collapsed regions due to their chaotic and random nature. It is precisely these emulation errors that we aim to correct using the COCA framework.



 §.§ COCA performance





We now turn to testing the use of our frame of reference emulator within a cosmological simulation. To do this, for each realisation of initial conditions in our test set, we run a reference simulation (see <ref>) as well as COCA and COLA simulations with varying specifications. For these runs, we use 20 time steps between a=0.05 and a=1, spaced linearly in scale factor, but we vary the number of force evaluations n_ f.
After some experimentation, we found that the best strategy to maximise the statistics described in <ref> is to place all force evaluations at the end of the simulation.
This is expected, as the dynamics become more non-linear at later times, making it crucial to accurately resolve particle trajectories during these periods, especially since the emulation error is also largest at these times (see <ref>).

In <ref>, we plot a slice of the final density field for one of the reference simulations in our test set, as well as the corresponding COCA simulations with n_ f = 10 and n_ f = 20, and their respective residuals relative to the reference. Both COCA simulations accurately recover the overall structure of the cosmic web, with correctly positioned filaments and nodes. With the smaller number of force evaluations, there is a small residual in the final density, but this has almost completely disappeared when n_ f = 20.



To assess the relative performance of COCA and COLA and to determine the optimal number of force evaluations, in <ref> we plot the fractional error on the matter power spectra and the cross-correlation coefficient for the a=1 density field as a function of wavenumber for both simulation frameworks on the test set.
As a sanity check, we verify that both COLA and COCA achieve similar performance when performing a force evaluation at each of the 20 time steps.
Our first observation is that COCA performs dramatically better than COLA, even when using few force evaluations. It is unsurprising that with n_ f=0 COLA performs poorly, as this is merely the LPT prediction, which is known to be a poor description at this redshift and on these scales. In contrast, we find that COCA with n_ f=0 is already extremely accurate: purely following the trajectories of the emulated frame of reference (n_ f=0) produces a behaviour practically identical to running a COLA simulation with n_ f=12 force evaluations. The matter power spectrum of the emulated field is 99% accurate up to k ≈ 0.3   h   Mpc^-1, with r^2(k) > 0.99 up to k ≈ 0.6   h   Mpc^-1. 
One would expect that, if the training simulations and evolution used a higher-accuracy gravity solver (e.g. P^3M or a tree-based approach), COCA would outperform COLA. However, it is not possible to check this conjecture in this example, since both the frame of reference emulator and COCA solver are based on PM forces.
Despite the good predictions of the emulator, we see that the relative error on the power spectrum increases to more than 10% at k=1   h   Mpc^-1 when n_f=0.
However, the error is reduced to less than 1% up to k = 0.5   h   Mpc^-1 by adding just 8 force evaluations, and less than 1% up to k = 1   h   Mpc^-1 with 10 force evaluations, both in terms of the power spectrum and phase accuracy.
This feature highlights the benefit of the COCA framework: we use machine learning to provide good approximations to the true solution and can run a physical simulation to correct for any errors made, using far fewer force evaluations than is ordinarily required.



The same behaviour is observed when considering the three-point statistics. In <ref> we plot the bispectrum for the COCA simulations in the configurations outlined in <ref>. As with the power spectrum, reasonable agreement with the reference is achieved without any force corrections, with errors of the order of 5-10%. However, with just 8 force evaluations, one achieves close to perfect agreement with the reference for almost all configurations considered, with the only discrepancy occurring for k_ s > 1   h   Mpc^-1.



We now evaluate the accuracy of the simulated velocity fields by plotting the error on the power spectrum and cross-correlation coefficient for the final velocity potential in <ref>.
Velocity fields are very poorly predicted for all COLA simulations that skip force evaluations, with an under-prediction of power beyond k ≈ 0.1  h   Mpc^-1, and an over-prediction as one approaches the Nyquist frequency of our simulations.
The cross-correlation between the COLA velocities and the reference is also very low, with practically zero correlation at k = 1   h   Mpc^-1 when no force evaluations are used, and with r^2(k) ≈ 0.5 at this scale for n_ f = 12.
This is unsurprising, since this latter case is equivalent to initialising a COLA simulation with an LPT prediction at a redshift of z=1.5 and using 12 time steps; one would not expect the initial conditions of such a simulation to be reasonable, as this is well beyond the validity of LPT.
However, this problem is alleviated if one uses an emulated frame of reference. Using only the emulator (n_ f = 0) reduces the error on the velocity field power spectrum to approximately 5% at a a=1, which, although still reasonably large, is much smaller than what is found with COLA with up to n_ f=12.
The advantage of the COCA framework is particularly evident when varying n_ f, as the addition of just 6 force evaluations practically eliminates this error, reducing it to 1%.
Similarly, we find that the COCA fields are much more correlated with the reference, even when using far fewer force evaluations, with r^2(k) > 0.8 for all k ≲ 1   h   Mpc^-1 and for any number of force evaluations. The degree of correlation improves as one increases n_ f.

In summary, we find that our emulator can reasonably recover the density and velocity fields even without any correction. However, emulation errors of up to 𝒪(10%) remain, but these can be reduced to the sub-percent level with just 8 force evaluations. Thus, COCA is able to correct for mistakes made in the emulation of particle trajectories by running a simulation in the corresponding frame of reference.



 §.§ COCA with misspecified cosmological parameters





One of the key motivations behind the COCA framework is the concept of ML-safety. Although emulation techniques have previously been applied to predict the results of dark matter simulations <cit.>, there may be concerns that the emulated solutions might not match the truth if the initial conditions or cosmological parameters are “unusual,” i.e., unlike the training data.
The capacity of emulators to extrapolate was tested by <cit.> in the context of well-understood simple matter distributions that had not been seen during training. Furthermore, <cit.> found that their emulator performed well with initial conditions containing significantly less power than their training examples.
However, with regular emulators, it is not possible to test all possible configurations, and thus, in general, one can only hope that the model extrapolates well to the application of interest.
In contrast, COCA uses a frame of reference emulator but solves the fundamental equations of motion. Therefore, any extrapolation mistake made in the emulation should be automatically corrected, unlike with the use of an emulator alone.

Our frame of reference emulator was trained using simulations run at a single cosmology. To test its out-of-distribution behaviour and its use in the COCA formalism, we ran 50 additional test simulations with a different set of cosmological parameters: Ω_ b = 0.03, Ω_ m = 0.35, h = 0.7, n_ s = 0.99, and σ_8 = 0.9.
These parameters are chosen to be relatively extreme, yet still within the support of a moderately wide prior that could be used for a cosmological analysis.
We note that we use the correct cosmological parameters for producing the initial density field, obtaining the LPT displacement fields, and solving the equations of motion; the only place where cosmological parameters are misspecified is in the prediction of p_ ML.

In <ref>, we plot the fractional error on the power spectra and the cross-correlation coefficients for the density and velocity fields in COCA simulations. The reference is the COLA simulations run with the same initial conditions, which are not subject to model misspecification in this scenario.
Despite the relatively extreme cosmological parameters, the uncorrected fields (n_ f=0) yield reasonable power spectra and cross-correlation. The mean error on the density power spectrum is approximately 20% by k=1   h   Mpc^-1 with r^2(k) > 0.93 at these scales, while the velocity power spectrum has slightly smaller errors—around 10%—and r^2(k) ≈ 0.85 by k=1   h   Mpc^-1.
This moderate agreement with the truth is enabled by using the initial density field rather than the white noise field as the input to the emulator (see section <ref>).
Indeed, even if the initial density appears different from that of the training simulations, the emulator does not have to predict the relevant initial matter power spectrum, which contains the entire dependence on all cosmological parameters except Ω_ m.
Additionally, since one expects p_ ML to be sourced only by local contributions in Lagrangian coordinates, the sensitivity to cosmological parameters should be relatively small.
Similar moderately accurate extrapolation behaviour has also been observed in other cosmological simulation emulators <cit.>.

Despite the moderate performance of this emulator in the presence of cosmological parameter misspecification, without any force evaluations (i.e., with n_f = 0), the error on the matter power spectrum would be too large for current cosmological analyses <cit.>.
Therefore, relying solely on an emulator of particles' trajectories (i.e., a frame of reference emulator with n_f = 0) as a forward model would produce inappropriate results and would not be a safe use of machine learning.
However, trajectories can be rectified in the COCA framework by evaluating gravitational forces and solving for the residual displacements with respect to the emulated frame of reference.
In our test, using just 8 force evaluations is sufficient to achieve percent-level agreement in both P(k) and r^2(k) for the density field, for all k < 1   h   Mpc^-1 (see <ref>). The same conclusion is true for the velocity field up to k ≈ 0.6   h   Mpc^-1.
Thus, with only a small additional computational cost, we can convert an unsafe use of machine learning in cosmology into a well-behaved one, even when the emulator is applied outside the range of its training data. This is one of the main benefits of COCA compared to traditional emulators of N-body simulation results.



 §.§ COCA versus a Lagrangian displacement field emulator




In this work, we advocate for using an emulator for the frame of reference in N-body simulations, as it allows for correcting emulation errors by introducing force evaluations. This approach contrasts with previous emulators <cit.>, which directly predict the Lagrangian displacement field, i.e., the simulation output. This section compares the relative accuracy of these two approaches as a function of the number of force evaluations in COCA.

To investigate this question, we train a time-dependent emulator for the residual displacement field Ψ_ ML (q, a), defined as the difference between the true Lagrangian displacement field Ψ (q, a) and that predicted by LPT, Ψ_ LPT (q, a). We opted to train a new Ψ-emulator rather than directly compare COCA to existing literature results to minimise the impact of differences in gravity solvers, training set sizes, architecture choices, and training procedures. For a fair comparison, we trained our Ψ-emulator using the same simulations as for the frame of reference emulator, employing the same architecture and training procedure outlined in <ref> (with p replaced by Ψ in the loss function).

In a similar manner as before, we begin by normalising the target variable by defining the function ψ(a) such that

    Ψ_ ML (q, a) ≡ψ(a) Ψ̃_ ML (q, a),

where Ψ̃_ ML (q, a) has unit standard deviation. Applying symbolic regression to the function ψ(a), we find that it is well approximated by 

    ψ(a) ≈ a^ϕ_0 + ϕ_1,

with ϕ_0 = 1.2412539 and ϕ_1 = -0.05402543, yielding a root mean squared error of 5 × 10^-4.
We take the linear density field as input, but this time output Ψ̃_ ML (q, a).
We evaluate this emulator on the same test simulations as for the frame of reference emulator, and convert the returned Lagrangian displacements into an Eulerian density field using a cloud-in-cell estimator.
We compute the power spectrum, cross-correlation coefficient, and bispectra (see <ref>) and plot these in <ref>, where we compare against COCA without force evaluations (n_ f=0; using solely the frame of reference emulator) and both n_ f=4 and n_ f=8.
We perform this analysis for both the fiducial and misspecified cosmology (see <ref>).

When compared to COCA with n_ f = 0, the Lagrangian displacement field emulator more accurately recovers the reference density field. The power spectra of the two methods are relatively similar with fiducial cosmological parameters, but the difference becomes more pronounced when the cosmology is misspecified. For all other metrics, the Ψ-emulator produces summary statistics that are closer to the reference.
This behaviour is expected: the Ψ-emulator is designed to optimise the prediction of dark matter particle positions through its loss function, naturally resulting in an accurate density field. In contrast, the frame of reference emulator in COCA aims to match particle momenta p. Consequently, without force evaluations, emulation errors in p accumulate over time, reducing the quality of the final density field.

Although the Ψ-emulator performs better than COCA with no force evaluations, there is no way to correct its errors, meaning its performance cannot be improved. Conversely, in COCA, force evaluations can be added to correct the errors made by the frame of reference emulator.
<ref> shows that the addition of only four force evaluations results in performance nearly identical to that of the Ψ-emulator for the bispectra, and better results for two-point statistics with residual errors reduced by a factor of 1.4.
Residual errors almost entirely disappear when eight force evaluations are used in COCA: the final power spectrum P(k) has approximately four to five times smaller errors than the one derived from the Ψ-emulator at all scales for both cosmologies.
Thus, even with very limited additional computations beyond the emulation, the COCA framework outperforms a Lagrangian displacement field emulator.

We note that our displacement emulator is slightly less accurate than that of <cit.>, who also emulated a PM-like output. They achieved errors on P(k) of 0.8% and 4% at k=0.4   h   Mpc^-1 and k=0.7   h   Mpc^-1, respectively, whereas our emulator is accurate to 3% and 9% at these scales.
We attribute this discrepancy to our use of fewer training simulations (2,000 fields compared to 10,000 in ), the need for our emulator to learn time-dependence (only 100 of the 2,000 training fields are at a=1), and <cit.> employing a more optimised architecture and training schedule.
As mentioned in <ref>, since the aim of this paper is to demonstrate how to correct for emulation errors rather than produce the optimal emulator, we chose not to increase the number of training simulations or fine-tune the architecture, as our emulator is already of similar quality to those in the literature.
If a frame of reference emulator with performance similar to that of <cit.> were used in COCA, fewer time steps would be needed to correct for emulation errors, thereby achieving the same theoretical guarantees with reduced computational expense.



 §.§ Timing tests




To assess the computational performance of COCA, in <ref> we show the required amount of CPU/GPU time for each of the stages of the framework. To perform the timing tests, we use an Intel Xeon Gold 6230 processor with 40 CPU cores and an Nvidia V100 GPU. We compare the results of running COLA with 20 time steps (and 20 force evaluations) and COCA with the same number of time steps, but with varying numbers of force evaluations (0, 8, and 20). All other settings are identical to those in section <ref>.

The COLA simulation takes approximately 61 CPU-seconds to run, with almost half of this time spent on cloud-in-cell binning (converting particle positions to the density field).
In our test, running COCA with n_ f=0 is approximately four times faster. 
In the current implementation, the emulation and the simulation codes are disjoint, with the frame of reference being written to and then read from disk at each kick time step (a process responsible for 77% of the CPU/GPU time in this case).
Separating emulation and N-body evolution is not a fundamental requirement of the COCA framework: one could emulate on the fly, which would effectively reduce the input/output time to zero for a slightly higher memory cost (two p_ML fields need to be kept in memory for each kick operation, see <ref>.
Such an approach would make the n_ f=0 case 18 times faster than the COLA simulation, with this emulator.

To enable the safe use of an ML emulator, COCA relies on including a finite number of force evaluations.
Naturally, if one uses the same number of force evaluations, then COCA is more computationally expensive than COLA, since it must perform the same steps as COLA but with an additional emulation stage.
However, because the ML correction makes the frame of reference more accurate than LPT alone, the number of force evaluations can be reduced to approximately 8 (see <ref>).
With the current implementation of separate emulation and N-body evolution codes, the cost of COCA with n_f=8 is approximately two-thirds of the cost of the COLA simulation. If emulation were done on the fly, the time required for inputs (36% of the total time) would be eliminated, making COCA 2.3 times faster than COLA.

These timing improvements are expected to become more dramatic if COCA were extended to include a more accurate gravity solver, such as a P^3M or tree-based code.
The computational expense for computing the forces in these codes is significantly higher than in the PM-based model used in this work. Therefore, reducing the number of force evaluations would dramatically improve run time. We leave such an investigation to future work.



§ DISCUSSION AND CONCLUSION



In this paper, we have introduced COmoving Computer Acceleration (COCA), a hybrid formalism involving ML and N-body simulations. Unlike previous works that directly emulate the simulation output, COCA solves the dynamics of an N-body simulation using a machine-learnt frame of reference.
COCA can be seen as an improvement of COLA, which solves the dynamics of an N-body simulation in the LPT frame of reference.
By virtue of the principle of Galilean invariance, equations of motion can be solved in any frame of reference, making COCA a ML-safe framework.
COCA is the first framework to use physics to determine the (otherwise uncorrected) emulation error in N-body simulations using ML and correct for it.

The concept behind COCA is entirely independent of the N-body solver and of the ML emulation algorithm used. For this proof-of-concept, we employed a PM approach to solving the equations of motion and a V-net architecture for the frame of reference, with fixed cosmological parameters.
We have demonstrated that after the ML-prediction of the optimal frame of reference (the one in which all particles are at rest), running the N-body simulation corrects for potential emulation errors in the particle trajectories.
We have quantitatively shown that the number of force evaluations required to achieve a given accuracy is reduced compared to COLA. The frame of reference emulator achieves between 1% and 10% accuracy when used in isolation, but only eight force evaluations are needed to reduce emulation errors to the percent level, compared to a 100-time step COLA simulation. Therefore, COCA can be utilised as a cheap N-body simulator.
Furthermore, with eight force evaluations, COCA is four to five times more accurate than a Lagrangian displacement field emulator, when the frame of reference emulator and the Lagrangian displacement field emulator are trained using the same computational resources. This increased accuracy is due to COCA's ability to correct emulation errors and represents one of the main advantages of this framework compared to the direct emulation approach explored in earlier literature.

In <ref>, we demonstrated that our frame of reference emulator is moderately robust to changes in cosmological parameters (despite training at a fixed cosmology). However, the COCA framework can correct for extrapolation errors arising from applying the emulator outside the range of validity of the training simulations.
Even when the frame of reference is inaccurate (because the ML training/prediction and the N-body evolution use different cosmological parameters), we found that percent-level accuracy can be reached on final density and velocity fields up to k ≈ 0.6   h   Mpc^-1.
Thus, the COCA framework provides ML-safety even when models are required to extrapolate.
There is no fundamental reason why the emulator cannot depend on cosmological parameters, and future implementations of COCA can include these as additional style parameters of the neural network.

Our example focused on relatively small simulation volumes (with a side length of 128 h^-1  Mpc) compared to those required for modern-day surveys (typically several gigaparsecs in length). With the current memory limitations of GPUs, it is not possible to emulate the entire volume with a single emulator at the desired resolution.
As a workaround, <cit.> splits the volume into several padded sub-boxes and treats each one separately, relying on sequential predictions for particle displacements in each sub-box to cover the full volume.
Similarly, in COCA, one could predict the frame of reference for particles in each sub-box, and then solve the equations of motion in each sub-box independently.
This idea relates to the algorithm introduced by <cit.> for perfectly parallel N-body simulations using spatial comoving Lagrangian Acceleration (sCOLA).
There, a tiling of the simulation volume is used, and the evolution of tiles is spatially decoupled by splitting the Lagrangian displacement field into large and small-scale contributions.
In sCOLA, the frame of reference used in the evolution of tiles is given by LPT, but it could be easily replaced by a frame of reference including both LPT and an ML contribution, as introduced in this paper.
Such an approach would overcome the memory limitations of GPUs, which currently limit COCA to small simulation volumes.
An additional benefit of this approach would be the inexpensive generation of light-cones. Indeed, when using a tiling approach as with sCOLA, only one tile needs to be evolved to a redshift of zero; the tiles farthest from the observer only need to be evolved until they intersect the light-cone at higher redshifts.

It is important to emphasise that the specific implementation details used in this work are not requirements but just an example.
For instance, one could use a perturbation theory-informed integrator for the equations of motion <cit.>, an approach complementary to COLA for fast generation of approximate cosmological simulations.
Furthermore, instead of using training simulations run with a PM gravity solver, one could learn the frame of reference given simulations with higher force accuracy, for example using a P^3M or tree-based gravity solver.
Subsequently, solving the equations of motion in the emulated frame of reference with the same solver would result in simulations with similar accuracy to those of the training set, but with significantly reduced computational cost.
The guarantee that any emulation mistakes are removed asymptotically as the number of force evaluations increases—a central feature of COCA—will remain.
This ML-safety cannot be guaranteed through direct emulation of P^3M or tree-based simulations.
As with COLA, the COCA framework could be adapted to include more extended physical models, such as neutrinos, which induce a scale-dependent growth factor <cit.>.
Finally, although in this work we have focused on gravitational N-body simulations in a cosmological context, the approach of solving equations of motion in an emulated frame of reference could be applied to any kind of simulation involving interacting particles (e.g., electrodynamics, hydrodynamics, radiative transfer, magnetohydrodynamics).
We generally expect a reduction in computational demands while retaining physical guarantees of convergence to the truth.

Benefiting from its modest computational cost, COCA could be used in analyses of cosmological data using fully non-linear models. It could straightforwardly be used as a forward model in implicit likelihood inference algorithms such as delfi <cit.>, bolfi <cit.>, selfi <cit.>, or the LtU-ili pipeline <cit.>.
As COCA is an ML-safe framework, its use as a forward model cannot bias the inference result.
We also note that using a V-net emulator and a PM force solver, the entire COCA framework is differentiable. For the emulation of the frame of reference, differentiability is achieved via automatic differentiation. For the N-body evolution, differentiable PM simulators already exist <cit.>. Building upon these, future work could be dedicated to writing a differentiable COCA solver, which could be used in Bayesian large-scale structure inference using an explicit field-level likelihood <cit.>.

Machine learning offers great promise in the acceleration of forward modelling in the physical sciences.
The output of any ML model is usually an approximation with inevitable emulation errors. In this paper, we have shown that emulation errors are correctable in gravitational N-body simulations.
By solving the correct physical equations while using the ML solution as an approximation, one can exploit the speed of ML while retaining the safety of more traditional methods.  







§ THE ACTUAL EQUATIONS





 §.§ Model equations with COCA


Using the notations of <cit.> and <cit.>, we consider dark matter particles with positions x and momenta p in comoving coordinates. Denoting the scale factor as a and the over-density field as δ, the equations to be solved are:

    xa    =    D(a) p  with  D(a) ≡1/a^2 ℋ(a), 
    pa    =    K(a) ∇( Δ^-1δ)   with  K(a) ≡ -3/2Ω_m^(0)ℋ^(0)2/a ℋ(a),

where ℋ(a) ≡ a^' / a is the conformal Hubble factor (where a prime denotes a derivative with respect to conformal time) and Ω_m^(0) is the matter density parameter at the present time (a=1). For simplicity, we note ∇_x = ∇, Δ_x = Δ and δ(x,a) = δ.

With x(a) = x_LPT(a) + x_ML(a) + x_res(a) (denoting the Langrangian perturbation theory, machine-learnt, and residual contributions to the position, respectively), 
we note for each contribution y∈{LPT, ML, res}:

    x_ya≡D(a) p_y  and  p_ya = / a( 1/D(a)x_ya) ≡ - K(a) V[x_y](a),

where the differential operator V[·](a) is defined by

    V[·](a) ≡ - 1/K(a)/ a( 1/D(a) ·a) .

Analogously, one writes the momenta as p(a) = p_LPT(a) + p_ML(a) + p_res(a), and thus Eqs. (<ref>) and (<ref>) take the form

    xa    =    D(a) {p_res(a) + p_LPT(a) + p_ML(a) }, 
    p_resa    =    K(a) {[ ∇( Δ^-1δ) ](a) + V[x_LPT](a) +  V[x_ML](a) }.


The analytical properties of LPT are <cit.>:

    x_LPT(a)     =    q - D_1(a) Ψ_1 + D_2(a) Ψ_2 , 
    D(a) p_LPT    =     - D_1aΨ_1 + D_2aΨ_2 , 
    V[x_LPT](a)     =     - D_1(a) Ψ_1 + [D_2(a) - D_1^2(a)] Ψ_2 ,

where Ψ_1 and Ψ_2 are the time-independent first and second order displacements, with corresponding growth factors D_1 and D_2.
This gives

    xa    =    D(a) {p_res(a) + p_ML(a) } - D_1aΨ_1 + D_2aΨ_2, 
    p_resa    =    K(a) {[ ∇( Δ^-1δ) ](a) - D_1(a) Ψ_1 + [ D_2(a) - D_1^2(a)]  Ψ_2 + V[x_ML](a) }.


Furthermore, for any arbitrary positive function u of a, we can rewrite

    xa    =    D(a) u(a) {1/u(a)×p_res(a) + 1/u(a)×p_ML(a) } - D_1aΨ_1 + D_2aΨ_2, 
    p_resa    =    u(a)a{K(a)/ u(a)/ a×[ [ ∇( Δ^-1δ) ](a) - D_1(a) Ψ_1 + [ D_2(a) - D_1^2(a) ] Ψ_2 + V[x_ML](a) ] }.




 §.§ Time stepping with COCA


In this paper, we adopt the second order symplectic “kick-drift-kick” algorithm, also known as the leapfrog scheme <cit.>, to integrate the equations of motion, for a series of n+1 time steps t(a) between t_0=t(a_i) and t_n+1=t(a_f). This algorithm relies on integrating the model equations on small time steps and approximating the momenta and accelerations that appear in the integrands (the part between curly brackets in the model equations) by their value at some time within the interval.

The discrete versions of the COCA model equations (equations (<ref>)–(<ref>) or (<ref>)–(<ref>)) give the Drift and Kick operators for COCA:

    D(t_i^D,t_f^D,t^K) :         x(t_i^D) ↦x(t_f^D) = x(t_i^D) + α_p(t_i^D, t_f^D, t^K) p_res(t^K) - [ D_1 ]_t_i^D^t_f^DΨ_1 + [ D_2 ]_t_i^D^t_f^DΨ_2 
                         +   α_p(t_i^D, t_f^D, t^K) p_ML(t^K), 
    K(t_i^K,t_f^K,t^D) :         p_res(t_i^K) ↦p_res(t_f^K) = p_res(t_i^K) + β_δ(t_i^K, t_f^K, t^D) ×
                        {[ ∇(Δ^-1δ)](t^D) - D_1(t^D) Ψ_1 + [ D_2(t^D) - D_1^2(t^D) ] Ψ_2 + g_ML(t^D) } .

Using equations (<ref>)–(<ref>), the standard discretisation of operators <cit.> gives the time prefactors as <cit.>,

    α_p(t_i^D, t_f^D, t^K) ≡∫_t_i^D^t_f^DD(t̃)  t̃ = ∫_t_i^D^t_f^Dt̃/t̃^2ℋ(t̃),   β_δ(t_i^K, t_f^K, t^D) ≡∫_t_i^K^t_f^KK(t̃)  t̃ = -3/2Ω_m^(0)ℋ^(0)2∫_t_i^K^t_f^Kt̃/t̃ℋ(t̃) .

The arbitrary function u of a appearing in equations (<ref>)–(<ref>) can be used to improve upon the standard discretisation of operators <cit.>. Indeed, if during the time step, the terms between brackets in equations (<ref>)–(<ref>) are closer to constants than the terms between brackets in equations (<ref>)–(<ref>), the approximation will hold better.
Therefore, using equations (<ref>)–(<ref>), the modified discretisation of operators gives the time prefactors, for any positive function u of t, as <cit.>,

    α_p(t_i^D, t_f^D, t^K) ≡1/u(t^K)∫_t_i^D^t_f^DD(t̃) u(t̃)  t̃,   β_δ(t_i^K, t_f^K, t^D) ≡[ u(t_f^K) - u(t_i^K) ] ×K(t^D)/[ u(t)/ t ] (t^D) .

In this paper, consistently with earlier literature, we use u(t) ≡ a^n_LPT with n_LPT = -2.5 <cit.>.

The ML frame of reference gives particles an acceleration g_ML(t^D) which should satisfy

    ∫_t_i^K^t_f^KK(t) V[x_ML](t)  t = ∫_t_i^K^t_f^Ku(t)t×{K(t)/ u(t) /  tV[x_ML](t) } t ≈β_δ(t_i^K, t_f^K, t^D) g_ML(t^D) .


In the standard discretisation, the integral can be approximated by using the value of V[x_ML](t) at t^D (assuming it is constant during the time step), giving β_δ(t_i^K, t_f^K, t^D) V[x_ML](t^D) with the definition of β_δ(t_i^K, t_f^K, t^D) given in equation (<ref>). In the modified discretisation, the integral can be approximated using the value of K(t)/ u(t) /  tV[x_ML](t) at t^D, giving also β_δ(t_i^K, t_f^K, t^D) V[x_ML](t^D) but with the definition of β_δ(t_i^K, t_f^K, t^D) given in equation (<ref>). In both cases, we get

    g_ML(t^D) ≡V[x_ML](t^D) .

But the integral is also:

    ∫_t_i^K^t_f^KK(t) V[x_ML](t)  t = ∫_t_i^K^t_f^K - /t( 1/D(t)x_ML/t)  t =  ∫_t_i^K^t_f^K - d p_ML/t t = p_ML(t_i^K) - p_ML(t_f^K),

which gives the alternative form

    g_ML(t^D) ≡1/β_δ(t_i^K, t_f^K, t^D)[ p_ML(t_i^K) - p_ML(t_f^K) ].

As such, to use the COCA Kick and Drift operators (Eqs. (<ref>) and (<ref>)), one does not require to emulate both p_ML and g_ML, but one only needs a single emulator (for p_ML), which is evaluated at the kick time steps.

In the end, the time evolution between t_0 and t_n+1 is achieved by applying the following operator to the initial state {x(t_0),p(t_0) }:

    L_+(t_n+1) E(t_n+1,t_0) L_-(t_0),

where E(t_n+1,t_0) is the operator given by (see <ref>)

    K(t_n+1/2,t_n+1,t_n+1)D(t_n,t_n+1,t_n+1/2) [ ∏_i=0^n K(t_i+1/2,t_i+3/2,t_i+1) D(t_i,t_i+1,t_i+1/2) ] K(t_0,t_1/2,t_0),

and L_± will be defined in <ref>.



 §.§ Generic Drift and Kick operators for PM, COLA and COCA


The difference between the COCA Kick and Drift operators and the corresponding COLA operators <cit.> is the last term in each operator. Therefore, we can introduce generic operators, valid for both COLA and COCA: for any external momentum p_ext and acceleration g_ext,

    D(t_i^D,t_f^D,t^K) :         x(t_i^D) ↦x(t_f^D) = x(t_i^D) + α_p(t_i^D, t_f^D, t^K) p_res(t^K) + α_LPT1(t_i^D, t_f^D, t^K) Ψ_1 + α_LPT2(t_i^D, t_f^D, t^K) Ψ_2 
                         +   α_ext(t_i^D, t_f^D, t^K) p_ext(t^K), 
    K(t_i^K,t_f^K,t^D) :         p_res(t_i^K) ↦p_res(t_f^K) = p_res(t_i^K) + β_δ(t_i^K, t_f^K, t^D) g_δ(t^D) + β_LPT1(t_i^K, t_f^K, t^D) Ψ_1 + β_LPT2(t_i^K, t_f^K, t^D) Ψ_2 
                         +   β_ext(t_i^K, t_f^K, t^D) g_ext(t^D) ,

where

    α_LPT1(t_i^D, t_f^D, t^K)    ≡    - [ D_1 ]_t_i^D^t_f^D, 
    α_LPT2(t_i^D, t_f^D, t^K)    ≡   [ D_2 ]_t_i^D^t_f^D, 
    α_ext(t_i^D, t_f^D, t^K)    ≡   α_p(t_i^D, t_f^D, t^K)  for COCA or  0  for COLA,
    p_ext(t^K)     =    p_ML(t^K)  for COCA or 0 for COLA,
    g_δ(t^D)    ≡   [ ∇(Δ^-1δ)](t^D),
    β_LPT1(t_i^K, t_f^K, t^D)     =     - β_δ(t_i^K, t_f^K, t^D) D_1(t^D), 
    β_LPT2(t_i^K, t_f^K, t^D)     =    β_δ(t_i^K, t_f^K, t^D) [ D_2(t^D) - D_1^2(t^D) ], 
    β_ext(t_i^K, t_f^K, t^D)     =    β_δ(t_i^K, t_f^K, t^D)  for COCA or  0  for COLA,
    g_ext(t^D)     =    g_ML(t^D)  for COCA or 0 for COLA.


We note that these operators also remain valid for a standard PM algorithm, by setting α_LPT1(t_i^D, t_f^D, t^K), α_LPT2(t_i^D, t_f^D, t^K), β_LPT1(t_i^K, t_f^K, t^D), and β_LPT2(t_i^K, t_f^K, t^D) to zero.



 §.§ Machine learning prediction of the frame of reference in COCA


The goal in COCA is to find the frame of reference in which p_res is as small as possible. Therefore, the machine needs to predict:


  * At any “kick” time step t^K,

    p_ML(t^K) ≡1/α_p(t_i^D, t_f^D, t^K)[ x(t_f^D)- x(t_i^D) - α_LPT1(t_i^D, t_f^D, t^K) Ψ_1 - α_LPT2(t_i^D, t_f^D, t^K) Ψ_2  ] ≡p_res^COLA(t^K),

that is the momentum residual p_res^COLA(t^K) of COLA <cit.>.

  * At any “drift” time step t^D,

    g_ML(t^D)    ≡    - 1/β_δ(t_i^K, t_f^K, t^D)[ β_δ(t_i^K, t_f^K, t^D) g_δ(t^D) + β_LPT1(t_i^K, t_f^K, t^D) Ψ_1 + β_LPT2(t_i^K, t_f^K, t^D) Ψ_2 ] 
        =    1/β_δ(t_i^K, t_f^K, t^D)[ p_res^COLA(t_i^K) - p_res^COLA(t_f^K) ] ≡ - g_res^COLA(t^D),

that is the residual acceleration g_res^COLA(t^D) of COLA <cit.>, up to a minus sign.


From equations (<ref>) and (<ref>), we see that it is sufficient for the machine to predict the momentum residual p_res^COLA(t^K) of COLA at any “kick” time step t^K, as the accelerations g_res^COLA(t^D) can be derived from the momenta.



 §.§ Initial and final momenta of particles in COCA


In the initial conditions, we have p(t_0) = p_LPT(t_0) + p_ML(t_0), which means that the momentum residual in the COCA frame of reference, p_res(t_0) = p(t_0) - p_LPT(t_0) - p_ML(t_0), should be initialised to zero. Furthermore, if initial conditions are generated with LPT, the ML contribution is p_ML(t_0) = 0 initially and we recover p_res(t_0) = p(t_0) - p_LPT(t_0), as in COLA.

At the end, the momentum p_LPT(t_n+1) + p_ML(t_n+1) of the COCA frame of reference has to be added to p_res(t_n+1) to recover the full momentum of particles, p(t_n+1). These operations correspond respectively to the L_-(t_0) : p(t_0) ↦p_res(t_0) and L_+(t_n+1) : p_res(t_n+1) ↦p(t_n+1) operators, given by

    L_±(t) :   p(t) ↦p(t) ±p_LPT(t) ±p_ML(t) = p(t) ±1/D(t)( -D_1tΨ_1 + D_2tΨ_2 ) ±p_ML(t) .






§ STATEMENT OF CONTRIBUTION


DJB implemented the frame of reference emulator, ran the COCA and COLA simulations, produced the plots, and contributed to the analysis and interpretation of the results and the writing of the paper.
MC contributed to the model equations and the validation of COCA.
LD provided code for the bispectrum analyses, advised on the design of the emulator, and contributed to the interpretation of ML performance results.
FL conceived and designed the study, wrote the model equations, modified the  code to accept any arbitrary frame of reference as input, contributed to the analysis and interpretation of results, edited the manuscript, advised early-career authors, and secured funding.
All authors read and approved the final manuscript.



§ ACKNOWLEDGEMENTS


We thank
Guilhem Lavaux, Natalia Porqueres, Benjamin Wandelt, and Ewoud Wempe
for useful comments and suggestions.
DJB and LD are supported by the Simons Collaboration on “Learning the Universe.”
FL and MC acknowledge financial support from the Agence Nationale de la Recherche (ANR) through grant INFOCW, under reference ANR-23-CE46-0006-01. LD acknowledges support by the Swedish Research Council (VR) under the project 2020-05143 – “Deciphering the Dynamics of Cosmic Structure".
This work has received funding from the Centre National d’Etudes Spatiales (CNES).
This work was done within the Aquila Consortium (<https://www.aquila-consortium.org/>).

For the purposes of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.



§ REFERENCES



