© The Authors 2025

1 Introduction

On the largest scales, the ΛCDM (Λ with cold dark matter) model provides a robust description of our Universe, where dark energy is represented as a cosmological constant (Λ) and dark matter as non-relativistic cold dark matter (Planck Collaboration VI 2020; Alam et al. 2021; Riess et al. 2022). The ΛCDM model also assumes General Relativity (GR), which offers an accurate representation of the Universe by incorporating small perturbations on a Friedman–Lemaître–Robertson–Walker (FLRW) background (Green & Wald 2014). The formation of the large-scale structure is driven by these small perturbations, which grow over time through gravitational interactions, giving rise to the cosmic web observed today (Peebles 1980).

While linear theory can analytically describe the growth of perturbations on large scales above ~100 Mpc (Yoo et al. 2009; Bonvin & Durrer 2011; Challinor & Lewis 2011), the non-linear nature of structure formation on smaller scales necessitates more advanced approaches. Higher-order perturbation theory, both in its Eulerian (Peebles 1980; Bernardeau et al. 2002) and Lagrangian (Zel’dovich 1970; Buchert & Ehlers 1993) formulations, allows for a more accurate understanding down to scales of ~20 Mpc. Below this scale, structure formation becomes highly non-linear, and the evolution of gravitationally interacting particles can only be accurately modeled through N-body simulations (Efstathiou et al. 1985).

Over time, numerous cosmological N -body codes have been developed (Couchman 1991; Kravtsov et al. 1997; Knebe et al. 2001; Teyssier 2002; Ishiyama et al. 2009; Potter et al. 2017; Springel et al. 2021; Garrison et al. 2021), with recent simulations including trillions of particles (Ishiyama et al. 2021; Euclid Collaboration: Castander et al. 2025). While significant effort has been invested in developing efficient Newtonian simulations for cosmology (see Angulo & Hahn 2022 for a review), many N -body codes for alternative gravity theories are derived from Newtonian codes. For instance, ECOSMOG (Li et al. 2012), ISIS (Llinares et al. 2014), RayMOND (Candlish et al. 2015), and Phantom of RAMSES (Lüghausen et al. 2015) are all based on the RAMSES code (Teyssier 2002). The former two implement the f (R) model (Hu & Sawicki 2007) and the nDGP model (Dvali et al. 2000), while the latter two implement MOND (modified newtonian dynamics, Milgrom 1983). Although it is possible to integrate modified gravity theories into TreePM codes like Gadget (Springel et al. 2021), as demonstrated by MG-Gadget (Puchwein et al. 2013), RAMSES has emerged as the preferred code for such implementations. This is due to the fact that alternative theories of gravity generally introduce additional fields governed by non-linear partial differential equations, which cannot be efficiently solved using standard tree-based methods. As a result, RAMSES, a particle-mesh (PM) code with adaptive-mesh refinement (AMR) and a multigrid solver, becomes an ideal choice for implementing these features. However, PM-AMR, TreePM, or Fast Multipole Method (FMM) codes typically require significant computational resources, even when highly optimised. This challenge is exacerbated in simulations of alternative gravity theories; for instance, Euclid Collaboration: Adamek et al. (2025) showed that f (R) or nDGP simulations can run up to ten times slower than their Newtonian counterparts.

To address these computational demands, researchers have explored PM codes (Knebe et al. 2001; Merz et al. 2005; Feng et al. 2016; Adamek et al. 2016; Klypin & Prada 2018) to reduce the cost at the expense of small-scale accuracy. These faster, albeit less precise, simulations are suitable for specific applications, such as producing large numbers of realisations for covariance matrix estimation or calculating the ratio (or boost) of specific statistical quantities relative to a reference case exploring large parameter spaces. Furthermore, exact PM codes can accurately reproduce structure formation at small scales given a higher resolution of the uniform mesh. Consequently, such codes have been ideal for developing solvers for modified gravity theories (see Llinares 2018 for a review), and have been used to study their impact relative to Newtonian simulations (Valogiannis & Bean 2017; Winther et al. 2017; Hassani & Lombriser 2020; Ruan et al. 2022; Hernández-Aguayo et al. 2022).

Given the need for speed, most N-body simulations have traditionally been written in compiled languages such as Fortran, C, or C++. In contrast, Python has become the most popular programming language in data science, owing to its straightforward syntax, rapid development speed, and extensive community- driven libraries, especially in astronomy. This popularity has created a gap between simulators and the broader scientific field. Despite its advantages, Python is often viewed as a slow language when used natively. To address this, significant effort has gone into developing efficient Python libraries, either as wrappers for C-based codes (such as NumPy, Harris et al. 2020) or through compiling Python code to machine language, as seen in Numba (Lam et al. 2015) and Cython (Behnel et al. 2011). Recently, the latter approach was used to develop the P³M code CONCEPT (Dakin et al. 2022), demonstrating the viability of Python for high-performance applications.

In this paper, we present PySCo¹ (Python Simulations for Cosmology), a cosmological PM N -body code written in Python and utilising the Numba library for increased performance and multithreading. The paper is organised as follows: Section 2 introduces the different models implemented in PySCo, including the modified gravity theories f (R) from Hu & Sawicki (2007), MOND (Milgrom 1983), parameterised gravity (Amendola et al. 2008), and dynamic dark energy (Chevallier & Polarski 2001; Linder 2003). Section 3 details the structure and algorithms implemented in the code, covering initial condition generation and various N-body solvers. In Section 4, we validate PySCo against other codes and analyse the impact of different numerical methods on the matter power spectrum. Finally, we conclude in Section 5.

2 Theory

2.1 Newtonian gravity

Let us consider only scalar perturbations on FLRW metric, in the Newtonian gauge (Ma & Bertschinger 1995) $$\text{d}{s}^2=a^2(\eta )\left[-(1+2\psi )\text{d}{\eta }^2+(1-2\varphi )\text{d}{\mathit{x}}^2\right],$$(1)

with a the scale factor, η the conformal time and ψ and ф the Bardeen potentials (Bardeen 1980). In GR, we have ψ = ф, and the Einstein equation gives $${\nabla }^2\varphi -3\mathscr{H}\left({\varphi }^{\prime }+\mathscr{H}\varphi \right)=4\pi G{a}^2\delta \rho ,$$(2)

with ℋ = aH the conformal Hubble parameter, ф′ the time derivative of the potential and ρ the Universe’s components density. It is common to apply the Newtonian and quasi-static (neglecting time derivatives) approximations, which involve neglecting the second term on the left-hand side of Eq. (2), as it is only significant at horizon scales. In the context of Newtonian cosmology, the Einstein equation takes the same form as the classical Poisson equation, with an additional dependence on the scale factor $${\nabla }^2\varphi =4\pi G{a}^2\delta {\rho }_m,$$(3)

with ρ_m the matter density. There are, however, cosmological General Relativity (GR) simulations (Adamek et al. 2016; Barrera-Hinojosa & Li 2020) that solve the full Einstein equations in the weak-field limit, taking into account gauge issues (Fidler et al. 2015, 2016). Additionally, there are methods to interpret Newtonian simulations within a relativistic framework (Chisari & Zaldarriaga 2011; Adamek & Fidler 2019). In this work, however, we focus on smaller scales where relativistic effects are negligible, and the Newtonian approximation remains well justified.

2.2 Dynamical dark energy

The w₀w_aCDM model provides a useful phenomenological extension of the standard ΛCDM framework by offering a dynamic, time-dependent description of dark energy. In this model, the cosmological constant is replaced by a variable dark energy component, which affects the formation of cosmic structures by modifying the universe’s expansion history. For a flat geometry (Ω_k = 0), the Hubble parameter H(z) is given by $$H(ɀ)=H_0\sqrt{{\text{Ω}}_{m,0}{(1+ɀ)}^3+{\text{Ω}}_{r,0}{(1+ɀ)}^4+{\Omega }_{\Lambda }(ɀ)}$$(4)

where a subscript zero indicates a present-day evaluation, and Ω_m, Ω_r, and Ω_λ represent the density fractions of matter, radiation, and dark energy, respectively. In this model, the dark-energy density evolves with redshift, following the relation $${\Omega }_{\Lambda }(ɀ)={\text{Ω}}_{\text{Λ},0}\exp \left\{{\displaystyle {\int }_0^{ɀ}\frac{3\left[1+w\left({ɀ}^{\prime }\right)\right]\text{d}{ɀ}^{\prime }}{1+{ɀ}^{\prime }}}\right\}.$$(5)

Using the widely adopted CPL parametrisation (Chevallier & Polarski 2001; Linder 2003), the dark-energy equation of state is expressed as $$w(ɀ)=w_0+w_a\frac{ɀ}{1+ɀ}.$$(6)

This simple modification does not alter the Einstein field equations nor the equations of motion, and it recovers the standard ΛCDM model when w₀ = −1 and w_a = 0.

2.3 MOND gravity

The Modified Newtonian Dynamics (MOND) theory, introduced by Milgrom (1983), was proposed as a potential solution to the dark matter problem by suggesting a deviation from Newtonian gravity in a Universe where the matter content is entirely in the form of baryons (for a detailed review, see Famaey & McGaugh 2012). In MOND, Newton’s second law is modified as follows $$\mu \left(\frac{\left|\mathit{𝑔}\right|}{{\mathit{𝑔}}_0}\right)\mathit{𝑔}={\mathit{𝑔}}_N,$$(7)

where g and g_N represent the MOND and Newtonian accelerations, respectively, and g₀ is a characteristic acceleration scale, approximately g₀ ≈ cH₀/[2π] ≈ 10⁻¹⁰m.s⁻². The function µ(x) is an interpolating function that governs the transition between the Newtonian regime (where gravitational force scales as r⁻²) and the MOND regime (where the force scales as r⁻¹), with r being the separation between masses. The interpolating function has the following limits $$\mu (x)=\{\begin{array}{ll}1,\hfill & x\gg 1,\hfill \\ x,\hfill & x\ll 1.\hfill \end{array}$$(8)

Similarly, the inverse interpolating function ν(y) can be defined as $$\mathit{𝑔}=\nu \left(\frac{{\mathit{𝑔}}_N}{{\mathit{𝑔}}_0}\right){𝑔}_N,$$(9)

where ν(y) follows the limits $$\nu (y)=\{\begin{array}{ll}1,\hfill & y\gg 1,\hfill \\ y^{-1/2},\hfill & y\ll 1.\hfill \end{array}$$(10)

In the MOND framework, the classical Poisson equation is modified as follows (Bekenstein & Milgrom 1984) $$\nabla \left[\mu \left(\frac{|\nabla \varphi |}{g_0}\right)\nabla \varphi \right]=4\pi G\delta \rho ,$$(11)

a formulation known as AQUAL, derived from a quadratic Lagrangian. In this paper, however, we consider the QUMOND (quasi-linear MOND, Milgrom 2010) formulation, where the non-linearity in the Poisson equation is re-expressed in terms of an additional effective dark matter fluid in the source term. The modified Poisson equations are $${\nabla }^2{\varphi }^{\text{N}}=4\pi G\rho ,$$(12) $${\nabla }^2\varphi =\nabla \left[\nu \left(\left|\nabla {\varphi }^{\text{N}}\right|/{𝑔}_0\right)\nabla {\varphi }^{\text{N}}\right],$$(13)

where ф^N and ф are the Newtonian and MOND potentials, respectively. This system of equations is more convenient to solve, as it involves only two linear Poisson equations. Moreover, numerical simulations have demonstrated that the AQUAL and QUMOND formulations yield very similar results (Candlish et al. 2015).

The missing component in the MOND framework is the specific form of the ν(y) function. Several families of interpolating functions have been proposed, each with different characteristics

• Simple function: $$\nu (y)=\frac{1}{2}+\frac{\sqrt{1+4/y}}{2},$$(14)

  which corresponds to the simple function proposed by Famaey & Binney (2005), equivalent to µ(x) = x/(1 + x).

• The n-family: $$\nu (y)={\left[\frac{1}{2}+\frac{\sqrt{1+4/y^n}}{2}\right]}^{1/n},$$(15)

  a commonly used parametrisation for the interpolating function (Milgrom & Sanders 2008). The case n = 2 is particularly well studied and is known as the standard interpolating function (Begeman et al. 1991). Additionally, Milgrom & Sanders (2008) introduced other functional forms

• The β-family: $$\nu (y)={\left(1-e^{-y}\right)}^{-1/2}+\beta e^{-y},$$(16)

• The γ-family: $$\nu (y)={\left(1-e^{-y^{\gamma /2}}\right)}^{-1/\gamma }+\left(1-{\gamma }^{-1}\right)e^{-y^{\gamma /2}},$$(17)

• The δ-family: $$\nu (y)={\left(1-e^{-y^{\delta /2}}\right)}^{-1/\delta },$$(18)

  which is a subset of the γ-family. While we focus on non- relativistic formulations of MOND for simplicity, it is important to acknowledge that various relativistic frameworks have been developed (Bekenstein 2006; Milgrom 2009; Skordis & Złośnik 2021), along with a recent generalisation of QUMOND (Milgrom 2023). These topics, however, are beyond the scope of this paper.

2.4 Parametrised gravity

A straightforward and effective phenomenological approach to modifying the theory of gravity is through the µ − σ parametrisation of the Einstein equations (Amendola et al. 2008). This is particularly useful when considering unequal Bardeen potentials ф ≠ ψ. Under the Newtonian and quasi-static approximations, the Einstein equations can be expressed as $${\nabla }^2\psi =4\pi G\mu (a)a^2\delta {\rho }_m,$$(19) $${\nabla }^2(\psi +\varphi )=8\pi \text{G}\mathrm{\Sigma }(a)a^2\delta {\rho }_m,$$(20)

where µ(a) represents the time-dependent ‘effective gravitational coupling’, which can be interpreted as a modification of the gravitational constant, and Σ(a) is the ‘light deflection parameter’. Since our focus is on the evolution of dark-matter particles, we only need to implement Eq. (19), which involves µ(a). In practice, the gravitational coupling µ(a) could be a function of both time and scale in Fourier space, µ(a, k), as in the ‘effective-field theory of dark energy’ (Frusciante & Perenon 2020). However, for simplicity, we prefer methods that can be solved numerically in both Fourier and configuration space. The inclusion of scale-dependent corrections will be considered in future work.

For the functional form of µ(a), we use the parametrisation from Simpson et al. (2013); Planck Collaboration XIV (2016); Planck Collaboration VI (2020); Abbott et al. (2019), which allows for deviations from GR during a dark-energy dominated era $$\mu (a)=1+{\mu }_0\frac{{\Omega }_{\Lambda }(a)}{{\text{Ω}}_{\text{Λ},0}},$$(21)

where µ₀ is the only free parameter, representing the gravitational coupling today.

2.5 f(R) gravity

In f (R) gravity, the Lagrangian extends the Einstein-Hilbert action (GR) by including an arbitrary function of the Ricci scalar curvature R. The total action is given as (Buchdahl 1970; Sotiriou & Faraoni 2010) $$S={\displaystyle \int {\text{d}}^4}x\sqrt{-𝑔}\left[\frac{R+f(R)}{16\pi G}+{ℒ}_m\right],$$(22)

where ℒ_m represents the matter Lagrangian, G is the gravitational constant, g is the determinant of the metric, and f (R) is the additional function of the curvature, which reduces to −2Λ in the standard ΛCDM model. A commonly used parametrisation of f (R) gravity is provided by Hu & Sawicki (2007), with the following functional form $$f(R)=-m^2\frac{c_1{\left(R/m^2\right)}^n}{c_2{\left(R/m^2\right)}^n+1},$$(23) $$\approx -\frac{c_1}{c_2}{m}^2+\frac{c_1}{c_2^2}{m}^2{\left(\frac{m^2}{R}\right)}^n,R\gg m^2,$$(24)

where n, c₁ and c₂ are the model parameters, and m represents the curvature scale, given by $m^2={\text{Ω}}_m{H}_0^2/c^2=8\pi G{\overline{\rho }}_{m,0}/\left(3c^2\right)$, where ${\overline{\rho }}_{m,0}$ is the current mean matter density. This model incorporates a Chameleon screening mechanism (Khoury & Weltman 2004; Burrage & Sakstein 2018) to suppress the fifth force caused by the scalar field f_R (also known as the ‘scalaron’). The scalaron is given by $$f_R=\frac{\text{d}f(R)}{\text{d}R}\approx -n\frac{c_1}{c_2^2}{\left(\frac{m^2}{R}\right)}^{n+1},$$(25)

which allows the theory to recover GR in high-density environments, ensuring consistency with Solar System tests. Observational evidence for dark energy as a cosmological constant imposes the constraint $$\frac{c_1}{c_2}=6\frac{{\text{Ω}}_{\text{Λ},0}}{{\text{Ω}}_{m,0}},$$(26)

and we can also express $$\frac{c_1}{c_2^2}=-\frac{3}{n}{\left[1+4\frac{{\text{Ω}}_{\text{Λ},0}}{{\text{Ω}}_{m,0}}\right]}^{n+1}{f}_{R0},$$(27)

where f_R0 is the present-day value of the scalaron, with current observational constraints from galaxy clusters indicating log₁₀ f_R0 < −5.32 (Vogt et al. 2025). In this framework, the Poisson equation is modified compared to its Newtonian counterpart. There is an additional term that depends on the scalaron field $${\nabla }^2\varphi =\frac{16\pi G{a}^2}{3}\delta \rho -\frac{1}{6}\delta R,$$(28) $${\nabla }^2{f}_R=-\frac{8\pi G{a}^2}{3c^2}\delta \rho +\frac{1}{3}\delta R,$$(29)

where the difference in curvature is given by $$\delta R=R-\overline{R}=\overline{R}\left[{\left(\frac{{\overline{f}}_R}{f_R}\right)}^{1/(n+1)}-1\right].$$(30)

Here, $\overline{R}$ represents the background curvature, expressed as $$\overline{R}=3m^2\left(a^{-3}+4\frac{{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_m}\right),$$(31)

and ${\overline{f}}_R$ is the background value of f_R $${\overline{f}}_R={\left(\frac{{\overline{R}}_0}{\overline{R}}\right)}^{n+1}{f}_{R,0}.$$(32)

The primary observational distinction between f (R) gravity and GR lies in the enhanced clustering on small scales, with the amplitude and shape of these deviations being dependent on f_R0, despite the two models sharing the same overall expansion history.

3 Methods

This section reviews the numerical methods used in PySCo, from generating initial conditions to evolving dark-matter particles in N-body simulations across various theories of gravity.

PySCo is entirely written in Python and uses the open-source library Numba (Lam et al. 2015), which compiles Python code into machine code using the LLVM compiler. This setup combines Python’s high development speed and rich ecosystem with the performance of C/Fortran. To optimise performance, PySCo relies on writing native Python code with ‘for’ loops, similar to how it would be done in C or Fortran.

Numba integrates seamlessly with NumPy (Harris et al. 2020), a widely used package for numerical operations in Python. Parallelisation in PySCo is simplified: by replacing ‘range’ with ‘prange’ in loops, the code takes advantage of multi-core processing. Numba functions are typically compiled just-in-time (JIT), meaning they are compiled the first time the function is called. Numba infers input and output types dynamically, supporting function overloading for different types.

In PySCo, however, most functions are compiled ahead-of- time (AOT), meaning they are compiled as soon as the code is executed or imported. This is because the simulation uses 32-bit floating point precision for all fields, allowing for AOT compilation. Since the simulation operates on a uniform grid, unlike AMR simulations, there is no need for fine-grained grids. Therefore, using 32-bit precision is sufficient and does not result in any loss of accuracy. Additionally, 32-bit floats improve performance by enabling SIMD (Single Instruction, Multiple Data) instructions, which the compiler implicitly optimises for.

3.1 Units and conventions

We adopt the same strategy as RAMSES and use supercomoving units (Martel & Shapiro 1998), where the Poisson equation takes the same form as in classical Newtonian dynamics but includes a multiplicative scale factor $${\nabla }^2\tilde{\varphi }=\frac{3}{2}a{\mathrm{\Omega }}_m(\tilde{\rho }-1),$$(33)

where a tilde denotes a quantity in supercomoving units, a is the scale factor, ф the gravitational potential, and ρ the matter density. We also define conversion units from comoving coordinates and super-conformal time to physical SI units $$\tilde{x}=\frac{x}{x_{*}},\text{\hspace{0.17em}\hspace{0.17em}d}\tilde{t}=\frac{\text{d}t}{t_{*}},\tilde{v}=v\frac{t_{*}}{x_{*}},$$(34) $$\tilde{\varphi }=\varphi \frac{t_{*}^2}{x_{*}^2},\tilde{c}=c\frac{t_{*}}{a{x}_{*}},$$(35)

where $\tilde{x},\text{d}\tilde{t},\tilde{v},\tilde{\varphi }$ and $\tilde{c}$ represent the particle position, time and velocity, gravitational potential and speed of light in simulation units. The conversion factors are defined as follows $$x_{*}=100a{L}_{\text{box}}/H_0,$$(36) $$t_{*}=a^2/H_0,$$(37) $${\rho }_{*}={\mathrm{\Omega }}_m{\rho }_c/a^3,$$(38)

as the length, time and density conversion units to km², seconds and kg/m³ respectively, where L_box is the box length in comoving coordinates and H₀ is the Hubble parameter today (in seconds). The particle mass is given by $$m_{\text{part\hspace{0.17em}}}={\rho }_{*}{x}_{*}^3/N_{\text{part\hspace{0.17em}}},$$(39)

where N_part is the total number of particles in the simulation.

3.2 Data structure

In this section, we discuss how PySCo handles the storage of particles and meshes using C-contiguous NumPy arrays. This approach was chosen for its simplicity and readability, allowing functions in PySCo to be easily reused in different contexts.

For particles, the position and velocity arrays are stored with the shape [N_part, 3], where the elements for each particle are contiguous in memory. This format is more efficient than using a shape of [3, N_part], particularly for mass assignment, where operations are performed particle by particle. To further enhance performance, particles are ordered using Morton indices (also known as z-curve indices) rather than linear or random ordering. Morton ordering improves cache usage and, thus, increases performance. This ordering is applied every N_reorder steps, as defined by the user, to maintain good data locality and avoid performance losses (see also Appendix E). We did not consider space-filling curves with better data locality properties (such as the Hilbert curve), because the associated encoding and decoding algorithms are much more computationally expensive, and Morton curves already provide excellent data locality. While more complex data structures are available (such as linked lists, fully-threaded trees, octrees or kdtrees), preliminary tests showed that using Morton-ordered NumPy arrays strikes a good balance between simplicity and performance, without the overhead of creating complex structures.

For scalar fields on the grid, arrays are stored with a shape $\left[N_{\text{cells}}^{1/3},N_{\text{cells}}^{1/3},N_{\text{cells}}^{1/3}\right]$ using linear indexing. Although linear indexing does not offer optimal data locality, it is lightweight and does not require additional arrays to store indices. Moreover, it works well with predictable (optimizable) memory-access patterns, such as those used in stencil operators. For vector fields, such as acceleration, the arrays have a shape $\left[N_{\text{cells}}^{1/3},N_{\text{cells}}^{1/3},N_{\text{cells}}^{1/3},3\right]$, similar to the format used for particle arrays, to maintain consistency and performance.

3.3 Initial conditions

In this section, we describe how PySCo handles the generation of initial conditions for simulations, although it can also read from pre-existing snapshots in other formats. PySCo can load data directly from RAMSES/pFoF format used in the Ray- Gal simulations (Breton et al. 2019; Rasera et al. 2022) or from Gadget format using the Pylians library (Villaescusa-Navarro 2018>). PySCo computes the time evolution of the scale factor, growth factors, and Hubble parameters, with the Astropy library (Astropy Collaboration 2022) and internal routines (see Appendix A).

To generate the initial conditions, the code requires a linear power spectrum P(k, z = 0), which is rescaled by the growth factor at the initial redshift z_ini. Additionally, PySCo generates a realisation of Gaussian white noise W, which is used to apply initial displacements to particles.

Some methods generate Gaussian white noise in configuration space, which is particularly useful for zoom simulations (Pen 1997; Sirko 2005; Bertschinger 2001; Prunet et al. 2008; Hahn & Abel 2013). However, in PySCo, the white noise is computed as $$W(\mathit{k})=A(\mathit{k})e^{i\theta (\mathit{k})},$$(40)

with A(k) an amplitude drawn from a Rayleigh distribution given by ${ℛ}_{\text{dist}}=\sqrt{-\text{\hspace{0.17em}ln\hspace{0.17em}}\mathcal{U}]0,1]}$ with U]0,1] a uniform random sampling between 0 and 1, and θ(k) = U]0,1]. We also fix $W\left(k\right)=\overline{W(-k)}$ where a bar denotes a complex conjugate, to ensure that the configuration-space field is real valued. In our case, since we consider a regular grid with periodic boundary conditions, we can generate the white noise directly in Fourier space. An initial realisation of a density field is computed using $${\delta }_{\text{ini}}(\mathit{k})={(}^2\sqrt{P(k)}W(\mathit{k}),$$(41)

with k = |k|. This ensures that we recover the Gaussian properties $$\langle {\delta }_{\text{ini\hspace{0.17em}}}(\mathit{k})\rangle =0,$$(42) $$\langle {\delta }_{\text{ini\hspace{0.17em}}}(\mathit{k}){\delta }_{\text{ini\hspace{0.17em}}}\left({\mathit{k}}^{\prime }\right)\rangle =(2\pi {)}^{3}P(\mathit{k}){\delta }_D\left(\mathit{k}+k^{\prime }\right),$$(43)

where δ_D is a Dirac delta and δ_ini(0) = 0. We have also implemented the option to use ‘paired and fixed’ initial conditions (Angulo & Pontzen 2016). This method greatly reduces cosmic variance by running paired simulations with opposite phases, at the cost of introducing some non-Gaussian features. The concept here is that instead of averaging the product of modes to match the power spectrum, the individual modes are set directly to ${(2\pi )}^{3/2}\sqrt{P(k)}$. In practice, the density field is used to compute the initial particle displacement from a homogeneous distribution, rather than directly sampling δ_ini . We use Lagrangian perturbation theory (LPT), with options for first- order 1LPT (also called ‘Zel’dovich approximation’, Zel’dovich 1970), second-order 2LPT (Scoccimarro 1998; Crocce et al. 2006), or third-order 3LPT (Catelan 1995; Rampf & Buchert 2012). The displacement field at the initial redshift z_ini up to third order is expressed as $$\text{Ψ}\left({ɀ}_{\text{ini}}\right)=D_{+}^{(1)}{\text{Ψ}}^{(1)}+D_{+}^{(2)}{\text{Ψ}}^{(2)}+D_{+}^{(3a)}{\text{Ψ}}^{(3a)}+D_{+}^{(3b)}{\text{Ψ}}^{(3b)}+D_{+}^{(3c)}{\text{Ψ}}^{(3c)},$$(44)

with D₊ ≡ D₊(z_ini)⁽¹⁾ is the linear (first order) growth factor at the initial redshift and Ψ⁽ⁿ⁾ are the different orders of the displacement field at z = 0, which can be written as (Michaux et al. 2021) $${\text{Ψ}}^{(1)}=-\nabla {\varphi }^{(1)},$$(45) $${\text{Ψ}}^{(2)}=\nabla {\varphi }^{(2)},$$(46) $${\text{Ψ}}^{(3a)}=\nabla {\varphi }^{(3a)},$$(47) $${\text{Ψ}}^{(3b)}=\nabla {\varphi }^{(3b)},$$(48) $${\text{Ψ}}^{(3c)}=\nabla \times A^{(3c)},$$(49) $${\varphi }^{(1)}={\nabla }^{-2}{\delta }_{\text{ini\hspace{0.17em}}},$$(50) $${\varphi }^{(2)}=\frac{1}{2}{\nabla }^{-2}\left[{\varphi }_{,ii}^{(1)}{\varphi }_{,jj}^{(1)}-{\varphi }_{,ij}^{(1)}{\varphi }_{,ij}^{(1)}\right],$$(51) $${\varphi }^{(3a)}={\nabla }^{-2}\left[\det {\varphi }_{,ij}^{(1)}\right],$$(52) $${\varphi }^{(3b)}=\frac{1}{2}{\nabla }^{-2}\left[{\varphi }_{,ii}^{(2)}{\varphi }_{,jj}^{(1)}-{\varphi }_{,ij}^{(2)}{\varphi }_{ij}^{(1)}\right],$$(53) $$A^{(3c)}={\nabla }^{-2}\left[\nabla {\varphi }_{,i}^{(2)}\times \nabla {\varphi }_{,i}^{(1)}\right],$$(54)

and $$\begin{array}{cc}{\nabla }^{(2)}{\varphi }^{(2)}=& {\varphi }_{,xx}^{(1)}\left({\varphi }_{,yy}^{(1)}+{\varphi }_{,zz}^{(1)}\right)+{\varphi }_{,yy}^{(1)}{\varphi }_{,zz}^{(1)}\\ & -{\varphi }_{,xy}^{(1)}{\varphi }_{,xy}^{(1)}-{\varphi }_{,xz}^{(1)}{\varphi }_{,xz}^{(1)}-{\varphi }_{,yz}^{(1)}{\varphi }_{,yz}^{(1)},\end{array}$$(55) $$\begin{array}{cc}{\nabla }^2{\varphi }^{(3a)}=& {\varphi }_{,xx}^{(1)}{\varphi }_{,yy}^{(1)}{\varphi }_{,zz}^{(1)}+2{\varphi }_{,xy}^{(1)}{\varphi }_{,xz}^{(1)}{\varphi }_{,yz}^{(1)}\\ & -{\varphi }_{,yz}^{(1)}{\varphi }_{,yz}^{(1)}{\varphi }_{,xx}^{(1)}-{\varphi }_{,xz}^{(1)}{\varphi }_{,xz}^{(1)}{\varphi }_{,yy}^{(1)}-{\varphi }_{,xy}^{(1)}{\varphi }_{,xy}^{(1)}{\varphi }_{,zz}^{(1)},\end{array}$$(56) $$\begin{array}{l}{\nabla }^2{\varphi }^{(3b)}=\frac{1}{2}{\varphi }_{,xx}^{(1)}\left({\varphi }_{,yy}^{(2)}+{\varphi }_{zz}^{(2)}\right)+\frac{1}{2}{\varphi }_{,yy}^{(1)}\left({\varphi }_{,xx}^{(2)}+{\varphi }_{,zz}^{(2)}\right)\\ +\frac{1}{2}{\varphi }_{,zz}^{(1)}\left({\varphi }_{,xx}^{(2)}+{\varphi }_{,yy}^{(2)}\right)-{\varphi }_{,xy}^{(1)}{\varphi }_{,xy}^{(2)}-{\varphi }_{,xz}^{(1)}{\varphi }_{,xz}^{(2)}-{\varphi }_{,yz}^{(1)}{\varphi }_{,yz}^{(2)},\end{array}$$(57) $$\begin{array}{cc}{\nabla }^2{A}_x^{(3c)}=& {\varphi }_{,xy}^{(2)}{\varphi }_{xz}^{(1)}-{\varphi }_{,xz}^{(2)}{\varphi }_{,xy}^{(1)}\\ & +{\varphi }_{,yz}^{(1)}\left({\varphi }_{,yy}^{(2)}-{\varphi }_{,zz}^{(2)}\right)-{\varphi }_{,yz}^{(2)}\left({\varphi }_{,yy}^{(1)}-{\varphi }_{,zz}^{(1)}\right),\end{array}$$(58) $$\begin{array}{cc}{\nabla }^2{A}_y^{(3c)}=& {\varphi }_{,yz}^{(2)}{\varphi }_{,yx}^{(1)}-{\varphi }_{,yx}^{(2)}{\varphi }_{,yz}^{(1)}\\ & +{\varphi }_{,xz}^{(1)}\left({\varphi }_{,zz}^{(2)}-{\varphi }_{,xx}^{(2)}\right)-{\varphi }_{,xz}^{(2)}\left({\varphi }_{,zz}^{(1)}-{\varphi }_{,xx}^{(1)}\right),\end{array}$$(59) $$\begin{array}{cc}{\nabla }^2{A}_z^{(3c)}=& {\varphi }_{,xz}^{(2)}{\varphi }_{,yz}^{(1)}-{\varphi }_{,yz}^{(2)}{\varphi }_{,xz}^{(1)}\\ & +{\varphi }_{,xy}^{(1)}\left({\varphi }_{,xx}^{(2)}-{\varphi }_{,yy}^{(2)}\right)-{\varphi }_{,xy}^{(2)}\left({\varphi }_{,xx}^{(1)}-{\varphi }_{,yy}^{(1)}\right).\end{array}$$(60)

In practice, we compute all derivatives directly in Fourier space, since configuration-space finite-difference gradients would smooth the small scales and thus create inaccuracies in the initial power spectrum. The second- and third-order contributions in the initial conditions can be prone to aliasing effects due to the quadratic and cubic non-linearities involved. To mitigate this, we apply Orszag’s 3/2 rule (Orszag 1971), as suggested in Michaux et al. (2021). The impact of this correction is minimal, by around 0.1% on the power spectrum at small scales, when using 3LPT initial conditions and for a relatively late start with z_ini ≈ 10). The initial position and velocity are then (up to third order) $$x_{ini}=q+{\text{Ψ}}^{(1)}{D}_{+}^{(1)}+{\text{Ψ}}^{(2)}{D}_{+}^{(2)}$$(61) $$+{\text{Ψ}}^{(3a)}{D}_{+}^{(3a)}+{\text{Ψ}}^{(3b)}{D}_{+}^{(3b)}+{\text{Ψ}}^{(3c)}{D}_{+}^{(3c)},$$(62) $$v_{\text{ini}}={\text{Ψ}}^{(1)}H{f}^{(1)}{D}_{+}^{(1)}+{\text{Ψ}}^{(2)}H{f}^{(2)}{D}_{+}^{(2)}$$(63) $$+{\text{Ψ}}^{(3a)}H{f}^{(3a)}{D}_{+}^{(3a)}+{\text{Ψ}}^{(3b)}H{f}^{(3b)}{D}_{+}^{(3b)}+{\text{Ψ}}^{(3c)}H{f}^{(3c)}{D}_{+}^{(3c)},$$(64)

with H the Hubble parameter at z_ini, f_n the growth rate contribution to the n-th order, and x_uniform the position of cell centres (or cell edges, see also Appendix B.2) when N_part = N_cells. We internally compute the growth factor and growth rate contributions as described in Appendix A.

3.4 Integrator

In our simulations, we employ the second-order symplectic Leapfrog scheme, often referred to as Kick-Drift-Kick, to integrate the equations of motion for the particles. The steps in the scheme are as follows $$\begin{array}{cc}{\mathit{v}}_{i+1/2}={\mathit{v}}_i+{\mathit{a}}_i\text{Δ}t/2,& \text{Kick},\\ {\mathit{x}}_{i+1}={\mathit{x}}_i+{\mathit{v}}_{i+1/2}\text{Δ}t,& \text{Drift,}\\ {\mathit{v}}_{i+1}={\mathit{v}}_{i+1/2}+{\mathit{a}}_{i+1}\text{Δ}t/2,& \text{Kick},\end{array}$$(65)

where the subscript i indicates the integration step, while x, v and a are the particle positions, velocities and accelerations respectively. There are several ways to set the time step ∆t. Some authors use linear or logarithmic spacing, depending on the user input. In our case, we followed a similar strategy as RAMSES (Teyssier 2002) and use several time stepping criteria. The first criterion is based on a cosmological time step that guarantees the scale factor does not change by more than a specified amount (by default 2%, see also Appendix B.1), which is particularly effective at high redshift. The second criterion is based on the minimum free-fall time given by $$\text{Δ}{t}_{\text{ff}}=\sqrt{\frac{h}{\max \left(\left|a_i\right|\right)}}.$$(66)

with h the cell size. We select the smallest value between this criterion and the cosmological time step. Additionally, we implemented a third criterion based on particle velocities ∆t_vel = h/max(|v_i|), though we found this value often exceeds ∆t_ff in practice. This approach ensures that the time step dynamically adapts to the structuration of dark matter in the simulation. To further refine the time step, we multiply it by a user-defined Courant-like factor.

An interesting prospect for the future is the use of integrators based on Lagrangian perturbation theory (LPT) that could potentially reduce the number of time steps required while maintaining high accuracy (as suggested by Rampf et al. 2025). However, such methods couple the integration scheme with specific theories of gravity through growth factors. Since these factors may not always be accurately computed (for example, in MOND), we prefer to maintain the generality of the standard leapfrog integration scheme for now. We may explore the implementation of such LPT-based integrators in the future for specific theories of gravitation.

3.5 Iterative solvers

To displace the particles, we first need to compute the force (or acceleration). There are various algorithms available for this purpose, either computing the force directly from a particle distribution, or determining the gravitational potential from which the force can subsequently be derived. When using the gravitational potential approach, the force can be recovered by applying a finite-difference gradient operator g = −∇ϕ. Given that our results are sensitive to the order of the operator used, we have implemented several options for central difference methods, each characterised by specific coefficients. These coefficients are detailed in Table 1. We aim to solve the following problem $$ℒu=f$$(67)

Where u is unknown, f is known and L is an operator. For the classical Poisson equation, u ≡ ϕ, f ≡ 4πGρ and L ≡^∇2 with the seven-point Laplacian stencil $${\nabla }^2{u}_{i,j,k}=\frac{1}{h^2}\left(L_{i,j,k}(u)-6u_{i,j,k}\right),$$(68)

with the subscripts i, j, k the cell indices and L_i,j,k(u) = u_{i + 1,j,k} + u_{i,j +1, k} + u_{i,j, k +1} + u_{i−1,j, k} + u_{i,j −1, k} + u_{i,j,k − 1}.

Lastly, f is the density (source) term of Eq. (67), which is directly estimated from the position of dark-matter particles. In code units, the sum of the density over the full grid must be $${\displaystyle \sum _{i,j,k}}{\tilde{\rho }}_{i,j,k}=N_{\text{part}},$$(69)

where the density in a given cell is computed using the ‘nearestgrid point’ (NGP), ‘cloud-in-cell’ (CIC) or ‘triangular-shaped cloud’ (TSC) mass-assignment schemes $$W_{\text{NGP}}\left(x_i\right)=\{\begin{array}{cc}1& \text{\hspace{0.17em}if\hspace{0.17em}}\\ 0& \text{\hspace{0.17em}otherwise,\hspace{0.17em}}\end{array}\left|x_i\right|<0.5,$$(70) $$W_{\text{CIC}}\left(x_i\right)=\{\begin{array}{cc}1-& \left|x_i\right|\\ 0& \end{array}\begin{array}{c}\text{if}\\ \text{otherwise}\end{array}\left|x_i\right|<1.0,$$(71) $$W_{\text{TSC}}\left(x_i\right)=\{\begin{array}{ccc}0.75-x_i^2\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}& \left|x_i\right|<0.5,\\ {\left(1.5-\left|x_i\right|\right)}^2/2\hfill & \text{\hspace{0.17em}else\hspace{0.17em}if\hspace{0.17em}}& 0.5<\left|x_i\right|<1.5,\\ 0\hfill & \text{\hspace{0.17em}otherwise}.& \end{array}$$(72)

Here, x_i is the normalised separation between a particle and cell positions, scaled by the cell size. In a three-dimensional space, this implies that a dark-matter particle contributes to the density of one, eight, or twenty-seven cells depending on whether the NGP, CIC, or TSC scheme is employed, respectively. For consistency, we use the same (inverse) scheme to interpolate the acceleration from cells to particles.

3.6 The Jacobi and Gauss-Seidel methods

Let us consider $ℒ=\left[\begin{array}{ccccc}{\ell }_{11}& {\ell }_{12}& {\ell }_{13}& \dots & {\ell }_{1n}\\ {\ell }_{21}& {\ell }_{22}& {\ell }_{23}& \dots & {\ell }_{2n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\ell }_{n1}& {\ell }_{n2}& {\ell }_{n3}& \dots & {\ell }_{nn}\end{array}\right],u=\left[\begin{array}{c}{u}_1\\ u_2\\ ⋮\\ u_n\end{array}\right]$ and $f=\left[\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_n\end{array}\right]$ The Jacobi method is a naive iterative solver which, for Eq. (67), takes the form $${\ell }_{11}{u}_1^{\text{new}}+{\ell }_{12}{u}_2^{\text{old}}+{\ell }_{13}{u}_3^{\text{old}}+\cdots +{\ell }_{1n}{u}_n^{\text{old}}=f_1,$$(73) $${\ell }_{21}{u}_1^{\text{old}}+{\ell }_{22}{u}_2^{\text{new}}+{\ell }_{23}{u}_3^{\text{old}}+\cdots +{\ell }_{2n}{u}_n^{\text{old}}=f_2,$$(74) $$\begin{array}{l}{\ell }_{31}{u}_1^{\text{old}}+{\ell }_{32}{u}_2^{\text{old}}+{\ell }_{33}{u}_3^{\text{new}}+\cdots +{\ell }_{3n}{u}_n^{\text{old}}=f_3,\\ ⋮\end{array}$$(75) $${\ell }_{n1}{u}_1^{\text{old\hspace{0.17em}}}+{\ell }_{n2}{u}_2^{\text{old\hspace{0.17em}}}+{\ell }_{n3}{u}_3^{\text{old\hspace{0.17em}}}+\cdots +{\ell }_{nn}{u}_n^{\text{new\hspace{0.17em}}}=f_n,$$(76)

where the superscripts ‘new’ and ‘old’ refer to the iteration. This means that for one Jacobi sweep, the new value is directly given by that of the previous iteration.

In practice, we use the Gauss-Seidel method instead, which has a better convergence rate and lower memory usage. This method is akin to Jacobi’s but enhances convergence by incorporating the most recently updated values in the computation of the remainder. It follows $${\ell }_{11}{u}_1^{\text{new\hspace{0.17em}}}+{\ell }_{12}{u}_2^{\text{old\hspace{0.17em}}}+{\ell }_{13}{u}_3^{\text{old\hspace{0.17em}}}+\cdots +{\ell }_{1n}{u}_n^{\text{old\hspace{0.17em}}}=f_1,$$(77) $${\ell }_{21}{u}_1^{\text{new\hspace{0.17em}}}+{\ell }_{22}{u}_2^{\text{new\hspace{0.17em}}}+{\ell }_{23}{u}_3^{\text{old\hspace{0.17em}}}+\cdots +{\ell }_{2n}{u}_n^{\text{old\hspace{0.17em}}}=f_2,$$(78) $${\ell }_{31}{u}_1^{\text{new\hspace{0.17em}}}+{\ell }_{32}{u}_2^{\text{new\hspace{0.17em}}}+{\ell }_{33}{u}_3^{\text{new\hspace{0.17em}}}+\cdots +{\ell }_{3n}{u}_n^{\text{old\hspace{0.17em}}}=f_3\text{,\hspace{0.17em}}$$(79) $$\begin{array}{l}⋮\\ {\ell }_{n1}{u}_1^{\text{new\hspace{0.17em}}}+{\ell }_{n2}{u}_2^{\text{new\hspace{0.17em}}}+{\ell }_{n3}{u}_3^{\text{new\hspace{0.17em}}}+\cdots +{\ell }_{nn}{u}_n^{\text{new\hspace{0.17em}}}=f_n,\end{array}$$(80)

which seems impossible to parallelise because each subsequent equation relies on the results of the previous one. However, we implement a strategy known as ‘red-black’ ordering in the Gauss-Seidel method. This technique involves colouring the cells like a chessboard, where each cell is designated as either red or black. When updating the red cells, we only utilise information from the black cells, and vice versa. This approach is equivalent to performing two successive Jacobi sweeps—one for the red cells and another for the black cells. For the Laplacian operator as expressed in Eq. (68), the Jacobi sweep for a given cell can be formulated as follows $$u_{i,j,k}^{\text{new\hspace{0.17em}}}=\frac{1}{6}\left(L_{i,j,k}\left(u^{\text{old\hspace{0.17em}}}\right)-h^2{f}_{i,j,k}\right).$$(81)

In scenarios involving a non-linear operator, finding an exact solution may not be feasible. In such cases, we linearise the operator using the Newton-Raphson method, expressed as $$u^{\text{new\hspace{0.17em}}}=u^{\text{old\hspace{0.17em}}}-\frac{ℒ\left(u^{\text{old\hspace{0.17em}}}\right)}{\partial ℒ/\partial u^{\text{old\hspace{0.17em}}}},$$(82)

This approach is commonly referred to as the ‘Newton Gauss- Seidel method’, allowing for more effective convergence in nonlinear contexts.

3.7 Multigrid

In practice, both the Jacobi and Gauss-Seidel methods are known for their slow convergence rates, often requiring hundreds to thousands of iterations to reach a solution and typically unable to achieve high accuracy. To overcome these limitations, a popular and efficient iterative method known as ‘multigrid’ is employed. This algorithm significantly accelerates convergence by solving the equation iteratively on coarser meshes, effectively addressing large-scale modes. The multigrid algorithm follows this procedure (Press et al. 1992): we first discretise the problem on a regular mesh with grid size h as $${ℒ}_h{u}_h=f_h,$$(83)

which can be solved using the Gauss–Seidel method, and where L_h is the numerical operator on the mesh with grid size h which approximates L. ũ_h denotes our approximate solution and v_h is the ‘error’ on the true solution $$v_h=u_h-{\tilde{u}}_h,$$(84)

and the ‘residual’ is given by $$d_h=f_h-ℒ{\tilde{u}}_h.$$(85)

Depending on whether the operator ℒ is linear or non-linear, different multigrid schemes will be employed to solve the equations effectively.

3.7.1 Linear multigrid

Considering the case where ℒ is linear, meaning that ℒ_h(u_h – ũ_h) ℒ_hu_{h –} ℒ_hũ_h, we have the relation $${ℒ}_h{v}_h=d_h.$$(86)

From there, the goal is to estimate v_h to find u_h. We use numerical methods to find the approximate solution ${\tilde{v}}_h$ by solving $${ℒ}_h{\tilde{v}}_h=d_h,$$(87)

using Gauss-Seidel. The updated approximation for the field is $${\tilde{u}}_h^{\text{new\hspace{0.17em}}}={\tilde{u}}_h+{\tilde{v}}_h.$$(88)

The issue is that the approximate operator L_h is usually local and finite-difference based, for which long-range perturbations are slowly propagating and therefore very inefficient computationally. The multigrid solution to this issue is to solve the error on coarser meshes to speed up the propagation of long-range modes. First, we use the ‘restriction’ operator R which interpolates from fine to coarse grid $$d_H=ℛd_h,$$(89)

where H = 2h is the grid size of the coarse mesh. We then solve Eq. (87) on the coarser grid to infer ${\tilde{v}}_H$ as $${ℒ}_H{\tilde{v}}_H=d_H.$$(90)

We then use the prolongation operator 𝒫, which interpolates from coarse to fine grid $${\tilde{v}}_h=𝒫{\tilde{v}}_H,$$(91)

and we finally update our approximation on the solution $${\tilde{u}}_h^{\text{new\hspace{0.17em}}}={\tilde{u}}_h+{\tilde{v}}_h.$$(92)

We have provided a brief overview of the multigrid algorithm using two grids, but in practice, to solve for ${\tilde{v}}_H$ , we can extend the scheme to even coarser meshes. This results in a recursive algorithm where the coarser level in PySCo contains 16³ cells. There are several strategies for navigating the different mesh levels, commonly referred to as V, F, or W-cycles. The V cycle is the simplest and quickest to execute, although it converges at a slower rate than the F and W cycles (see Appendix C).

For the restriction and prolongation operators, the lowest- level schemes employ averaging (where the field value on the parent cell is the mean of its children’s values) and straight injection (where child cells inherit the value of their parent), respectively. However, as noted by Guillet & Teyssier (2011), to minimise inaccuracies in the estimation of the final solution, we use a higher-order prolongation scheme defined as $$𝒫=\{\begin{array}{ll}27/64,\hfill & x_{0,1,2}<0.5,\hfill \\ 9/64,\hfill & x_{0,1}<0.5<x_2<1.5,\hfill \\ 3/64,\hfill & x_0<0.5<x_{1,2}<1.5,\hfill \\ 1/64,\hfill & 0.5<x_{0,1,2}<1.5,\hfill \\ 0,\hfill & \text{\hspace{0.17em}otherwise\hspace{0.17em}},\hfill \end{array}$$(93)

where x = (x₀, x₁, x₂) with x₀ ≤ x₁ ≤ x₂, is the separation (normalised by the fine grid size) between the centre of the fine and coarse cells.

We consider our multigrid scheme to have converged when the residual is significantly lower than the ‘truncation error’, defined as: $$\tau \equiv {ℒ}_{h}u-f_h,$$(94)

and which can be estimated as (Press et al. 1992; Li et al. 2012) $${\tilde{\tau }}_h\equiv {ℒ}_H\left(ℛ{\tilde{u}}_h\right)-ℛ{ℒ}_h{\tilde{u}}_h.$$(95)

We consider that we reached convergence when $$\left|d_h\right|<\alpha {\tilde{\tau }}_h,$$(96)

where |d_h | is the square root of the quadratic sum over the residual in each cell of the mesh, α is the stopping criterion. It is noteworthy that (Knebe et al. 2001) proposed using an alternative approach for estimating the truncation error $${\tau }_{h,K01}=𝒫\left({ℒ}_H\left(ℛ{\tilde{u}}_h\right)\right)-{ℒ}_h{\tilde{u}}_h,$$(97)

which we can approximate by $${\tau }_{h,K01}\approx 𝒫\left(ℛf_h\right)-f_h,$$(98)

to reduce computational time, when the source term is non-zero. The relation between these two approaches is τ_h ≈ 0.1τ_h,K03. In practice, we use the first estimation τ_h. Since the Jacobi and Gauss-Seidel methods are iterative, an initial guess is still required. While the multigrid algorithm is generally insensitive to the initial guess in most cases, some small optimisations can be made. For instance, if we initialise the full grid with zeros, one Jacobi step directly provides $$u_{i,j,k}^{\text{ini}}=-\frac{h^2}{6}{f}_{i,j,k},$$(99)

which is used to initialise u_h for the first simulation step, and v_H the error on coarser meshes.

In addition to the Gauss-Seidel method, RAMSES also incorporates ‘successive over-relaxation’ (SOR), which allows the updated field to include a contribution from the previous iteration, expressed as ϕⁿ⁺¹ = ωϕⁿ⁺¹ + (1 − ω)ϕⁿ, where n denotes the iteration step, and ω is the relaxation factor. Typically, we perform two Gauss-Seidel sweeps before the restriction (termed pre-smoothing) and one sweep after the prolongation (postsmoothing), controlled by the parameters Npre and Npost. We set ω = 1.25 by default (see also Appendix C), similarly to Kravtsov et al. (1997).

A pseudo code for the V cycle is depicted in Fig. 1, where the components that require modification based on the gravitational theory are the smoothing and residual functions. As noted by Press et al. (1992), omitting SOR (setting ω = 1) for the first and last iterations can enhance performance since this allows the prolongation operator (line 15) to operate on only half the mesh, since the other half will be directly updated during the first Gauss-Seidel sweep. Similarly, the restriction operator can also be applied to only half of the mesh because, following one Gauss-Seidel sweep, the residual will be zero on the other half by design. It is thus feasible to combine lines 6 and 7 (the residual and restriction operators) to compute the coarser residual from half of the finer mesh. In this example, we use a zero initial guess for the error at the coarser level, but in principle, we could initialise it with a Jacobi step (as shown in Eq. (99) for Newtonian gravity).

This optimisation, however, should be approached cautiously, as the residual is zero up to floating-point precision. In certain scenarios, this can significantly impact the final solution. For example, in PySCo, we primarily utilise 32-bit floats, rendering this optimisation less accurate for the largest modes. Consequently, we ultimately decided to retain only the optimisation concerning the prolongation on half the mesh when employing the non-linear multigrid algorithm, as we do not utilise SOR in this context (see Appendix C).

Fig. 1 Pseudo code in Python for the multigrid V-cycle algorithm. We highlight in cyan the lines where the theory of gravitation enters.

3.7.2 Non-linear multigrid

If ℒ is a non-linear operator instead, then we need to solve for $${ℒ}_h\left({\tilde{u}}_h+v_h\right)-{ℒ}_h{\tilde{u}}_h=d_h.$$(100)

Contrarily to the linear case in Eq. (90) where we only need to solve for the error at the coarser level, we now use $${ℒ}_H{\tilde{u}}_H={ℒ}_H\left(ℛ{\tilde{u}}_h\right)+ℛd_h,$$(101)

where we need to store the full approximation of the solution ũ at every level, hence the name ‘Full Approximation Storage’. Finally, we update the solution as $${\tilde{u}}_h^{\text{new\hspace{0.17em}}}={\tilde{u}}_h+𝒫\left({\tilde{u}}_H-ℛ{\tilde{u}}_h\right).$$(102)

3.8 Fast Fourier transforms

We implemented three different fast Fourier transforms (FFT) procedures to compute the force field, most of which were already implemented in FastPM (Feng et al. 2016).

– FFT: The Laplacian operator is computed through the Green’s function kernel $${\nabla }^{-2}=-k^{-2}{W}_{\text{MAS}}^{-2}(k),$$(103)

where W_MAS(k) is the mass-assignment scheme filter, given by $$W_{\text{MAS}}(k)={\left[{\displaystyle \prod _{d=x,y,z}\mathrm{sinc}}\left(\frac{{\omega }_d}{2}\right)\right]}^p,$$(104)

with ω_d = k_dh between [−π, π] and where p = 1,2,3 for NGP, CIC and TSC respectively (Hockney & Eastwood 1981; Jing 2005).

– FFT_7pt: instead of using the exact kernel for the Laplacian operator, we use the Fourier-space equivalent of the seven-point stencil Laplacian (Eq. (68)), which reads $${\nabla }^{-2}=-{\left[{\displaystyle \sum _{d=x,y,z}{\left(h{\omega }_d\mathrm{sinc}\frac{{\omega }_d}{2}\right)}^2}\right]}^{-1},$$(105)

where no mass-assignment compensation is used. For both FFT and FFT_7pt methods, the force is estimated through the finite- difference stencil as shown in Table 1.

– FULL_FFT: The force is directly estimated in Fourier space through the differentiation kernel $$\nabla {\nabla }^{-2}=-ik{k}^{-2}{W}_{\text{MAS}}^{-2}(k).$$(106)

As we will see in Section 4.2, this naive operator can become very inaccurate when the simulation has more cells than particles.

3.9 Newtonian and parametrised simulations

In Eq. (99), we showed how to provide a generic initial guess for the multigrid algorithm. However, leveraging our understanding of the underlying physics can enable us to formulate an even more accurate first guess, thus reducing the number of multigrid cycles required for convergence. In the context of N -body simulations, we anticipate that the potential field will closely resemble that of the preceding step, especially when using sufficiently small time steps. This allows us to adopt the potential from the previous step as our initial guess for the multigrid algorithm. While this approach facilitates faster convergence to the true solution, it necessitates storing one additional grid in memory. Moreover, it is important to note that, in the linear regime in Newtonian cosmology, the density contrast evolves as a function of redshift with a scale-independent growth factor. Consequently, we can optimise our first guess by rescaling the potential field from the previous step according to the following equation $$\tilde{\varphi }\left({ɀ}_1\right)=\frac{\left(1+{ɀ}_1\right)D_{+}\left({ɀ}_1\right)}{\left(1+{ɀ}_0\right)D_{+}\left({ɀ}_0\right)}\tilde{\varphi }\left({ɀ}_0\right),$$(107)

where z₀ and z₁ denote the initial and subsequent redshifts, respectively. This rescaling is performed for every time step of the simulation, except for the initial time step, ensuring that we maintain an efficient and accurate estimate for our initial guess in the multigrid process. More details can be found in Appendix C.

Fig. 2 MOND interpolating functions for the different families and parameters shown in Sect. 2.3.

3.10 MOND simulations

The classical formulations of MOND (such as AQUAL and QUMOND) have already been implemented in several codes (Nusser 2002; Llinares et al. 2008; Angus et al. 2012; Candlish et al. 2015; Lüghausen et al. 2015; Visser et al. 2024). In PySCo, we further allow, if specified by the user, for a time-dependent acceleration scale g₀ → a^Ng₀, which basically delays (or accelerates) the entry of perturbations in the MOND regime, and which is set to N = 0 by default. In Fig. 2, we show the behaviour of the interpolating functions ν(y) described in Sect. 2.3. We see that, as expected from Eq. (13), ν(y) → 1 when y ≫ 1, meaning that we recover Newtonian gravity in a regime of strong acceleration.

To implement the right-hand side of Eq. (13), we employ a method analogous to that described by Lüghausen et al. (2015), and using the same notation $$\begin{array}{l}{\nabla }^2{\varphi }^{\text{N}}=\frac{1}{h}[{\nu }_{B_x}\nabla {\left({\varphi }^{\text{N}}\right)}_{B_x,x}-{\nu }_{A_x}\nabla {\left({\varphi }^{\text{N}}\right)}_{A_x,x}\\ +{\nu }_{B_y}\nabla {\left({\varphi }^{\text{N}}\right)}_{B_y,y}-{\nu }_{A_y}\nabla {\left({\varphi }^{\text{N}}\right)}_{A_y,y}\\ +{\nu }_{B_{ɀ}}\nabla {\left({\varphi }^{\text{N}}\right)}_{B_{ɀ},ɀ}-{\nu }_{A_{ɀ}}\nabla {\left({\varphi }^{\text{N}}\right)}_{A_{ɀ},ɀ}]\end{array}$$(108)

where B_i and A_i are the points at +0.5h e_i and −0.5h e_i respectively, with e_i ∈ {e_x, e_y, e_z}, the unit vectors of the simulation. We have also defined $${\nu }_{B_x}=\nu \left(\frac{\sqrt{{\left[\nabla {\left({\varphi }^{\text{N}}\right)}_{B_x,x}\right]}^2+{\left[\nabla {\left({\varphi }^{\text{N}}\right)}_{B_x,y}\right]}^2+{\left[\nabla {\left({\varphi }^{\text{N}}\right)}_{B_x,ɀ}\right]}^2}}{g_0}\right),$$(109)

with ∇(ϕ^N) the i-th component of the force (with a minus sign) at the position B_x . The force components are estimated as $$\nabla {(\varphi )}_{B_x,x}=\frac{{\varphi }_{1,0,0}^{\text{N}}-{\varphi }_{0,0,0}^{\text{N}}}{h},$$(110) $$\nabla {(\varphi )}_{B_x,y}=\frac{\left({\varphi }_{1,1,0}^{\text{N}}-{\varphi }_{1,-1,0}^{\text{N}}\right)+\left({\varphi }_{0,1,0}^{\text{N}}-{\varphi }_{0,-1,0}^{\text{N}}\right)}{4h}$$(111) $$\nabla {(\varphi )}_{B_x,ɀ}=\frac{\left({\varphi }_{1,0,1}^{\text{N}}-{\varphi }_{1,0,-1}^{\text{N}}\right)+\left({\varphi }_{0,0,1}^{\text{N}}-{\varphi }_{0,0,-1}^{\text{N}}\right)}{4h},$$(112)

and similarly for other points. We note that for B_x , x (as well as for B_y, y and B_z, z), we perform a three-point central derivative with half the mesh size, differing from the approach taken by (Lüghausen et al. 2015), who implemented a non-uniform five-point stencil. This decision was made to ensure that in the Newtonian case (that is, ν = 1), we exactly recover the sevenpoint Laplacian, maintaining consistency with the Laplacian operator employed in our multigrid scheme. Consequently, we opted to retain three-point derivatives for the other components as well (although we can still use a different stencil order when computing the acceleration from the MOND potential).

Additionally, we exclusively solve Eq. (108) using either the multigrid or FFT_7pt solvers. Using the FFT solver presents challenges because, although we deconvolve ϕ^N by the massassignment scheme kernel, the uncorrected mesh discreteness in the force computation introduces inaccuracies that can significantly affect the matter power spectrum.

3.11 f(R) simulations

In supercomoving units, the f (R) field equations from Eqs. (28)– (29) are given by (Li et al. 2012) $${\nabla }^2\tilde{\varphi }=\frac{3}{2}{\mathrm{\Omega }}_{m}a(\tilde{\rho }-1)-\frac{{\tilde{c}}^2}{2}{\nabla }^2{\tilde{f}}_R,$$(113) $${\nabla }^2{\tilde{f}}_R=-\frac{1}{{\tilde{c}}^2}{\mathrm{\Omega }}_{m}a(\tilde{\rho }-1)+\frac{1}{3{\tilde{c}}^2}R{a}^4\left[{\left(\frac{{\overline{f}}_R}{{\tilde{f}}_R}\right)}^{1/(n+1)}-1\right]$$(114)

where ${\tilde{f}}_R=a^2{f}_R{\tilde{c}}^2/c^2$. While Eq. (113) is linear and can be solved by standard techniques, Eq. (114) is not and needs some special attention. To this end, Oyaizu (2008) used a nonlinear multigrid algorithm (see Section 3.7.2) with the Newton- Raphson method (as shown in Eq. (82)), making the change of variable $u\equiv \ln \left({\tilde{f}}_R/{\overline{f}}_R\right)$ to avoid unphysical zero-crossing of f_R. Bose et al. (2017) proposed instead that for this specific model, one could perform a more appropriate change of variable, that we generalise here to $u\equiv {\left({\tilde{f}}_R/{\overline{f}}_R\right)}^{1/(n+1)}$. We can then recast Eq. (114) as $$u^{n+1}+pu+q=0,$$(115)

where $$p=\frac{h^2}{6{\tilde{c}}^2{\overline{f}}_R}\left[{\mathrm{\Omega }}_{m}a(1-\tilde{\rho })-\frac{a^4\overline{R}}{3}\right]-\frac{1}{6}{L}_{i,j,k}\left(u^{n+1}\right),$$(116) $$q=\frac{a^4{h}^2\overline{R}}{18{\tilde{c}}^2{\overline{f}}_R}.$$(117)

We note that q is necessarily negative (because ${\overline{f}}_R<0$), which is useful to determine the branch of the solution for Eq. (115).

– Case n = 1 : As noticed in Bose et al. (2017), when making the change of variable $u=\sqrt{-f_R}$, the field equation could be recast as a depressed cubic equation, $$u^3+pu+q=0,$$(118)

which possesses the analytical solutions (Ruan et al. 2022) $$u=\{\begin{array}{ll}{(-q)}^{1/3},\hfill & p=0,\hfill \\ -\left[C+{\mathrm{\Delta }}_0/C\right]/3,\hfill & p>0,\hfill \\ -\left[C+{\mathrm{\Delta }}_0/C\right]/3,\hfill & p<0\text{\hspace{0.17em}and\hspace{0.17em}}{\mathrm{\Delta }}_1^2-4{\mathrm{\Delta }}_0^3>0\hfill \\ -2\sqrt{{\mathrm{\Delta }}_0}\cos (\theta /3+2\pi /3)/3,\hfill & \text{\hspace{0.17em}else\hspace{0.17em}}\hfill \end{array}$$(119)

with ${\mathrm{\Delta }}_0=-3p,{\mathrm{\Delta }}_1=27q,C={\left[\frac{1}{2}\left({\mathrm{\Delta }}_1+\sqrt{{\mathrm{\Delta }}_1^2-4{\mathrm{\Delta }}_0^3}\right)\right]}^{1/3}$ and $\cos \theta ={\mathrm{\Delta }}_1/\left(2{\mathrm{\Delta }}_0^{3/2}\right)$.