PySCo: A fast Particle-Mesh $N$-body code for modified gravity simulations in Python

TL;DR

PySCo delivers a fast Python-based particle-mesh N-body code capable of simulating ΛCDM, dynamical dark energy, and a range of modified gravity theories (including $f(R)$, MOND, and parametrised gravity) using multigrid or FFT solvers. It combines LPT-based initial conditions, scalable data structures, and multiple solvers to enable rapid exploration of gravity and dark-energy parameter space, validated against RAMSES, e-MANTIS, and Phantom of RAMSES. The results demonstrate that with sufficient small-scale resolution, the late-time power spectrum is robust to initial redshift within ~0.1% for appropriate LPT order, while outlining tradeoffs between gradient order, solver type, and mass-assignment schemes. PySCo thus offers a practical, emulator-ready tool for exploring cosmological models, generating covariance matrices, and training data for ML applications, with a clear pathway to GPU acceleration and broader gravity theories in future work.

Abstract

We present PySCo, a fast and user-friendly Python library designed to run cosmological $N$-body simulations across various cosmological models, such as $Λ$CDM and $w_0w_a$CDM, and alternative theories of gravity, including $f(R)$, MOND and time-dependent gravitational constant parameterisations. PySCo employs Particle-Mesh solvers, using multigrid or Fast Fourier Transform (FFT) methods in their different variations. Additionally, PySCo can be easily integrated as an external library, providing utilities for particle and mesh computations. The library offers key features, including an initial condition generator based on up to third-order Lagrangian Perturbation Theory (LPT), power spectrum estimation, and computes the background and growth of density perturbations. In this paper, we detail PySCo's architecture and algorithms and conduct extensive comparisons with other codes and numerical methods. Our analysis shows that, with sufficient small-scale resolution, the power spectrum at redshift $z = 0$ remains independent of the initial redshift at the 0.1\% level for $z_{\rm ini} \geq$ 125, 30, and 10 when using first, second, and third-order LPT, respectively. Although the seven-point Laplacian method used in multigrid also leads to power suppression on small scales, this effect can largely be mitigated when computing ratios. In terms of performance, PySCo only requires approximately one CPU hour to complete a Newtonian simulation with $512^3$ particles (and an equal number of cells) on a laptop. Due to its speed and ease of use, PySCo is ideal for rapidly generating vast ensemble of simulations and exploring parameter spaces, allowing variations in gravity theories, dark energy models, and numerical approaches. This versatility makes PySCo a valuable tool for producing emulators, covariance matrices, or training datasets for machine learning.

Figures (18)

• Figure 1: Pseudo code in Python for the multigrid V-cycle algorithm. We highlight in cyan the lines where the theory of gravitation enters.
• Figure 2: MOND interpolating functions for the different families and parameters shown in Sec. \ref{['sec:theory_mond']}.
• Figure 3: Ratio of the matter power spectrum at $z = 0$ for different starting redshift $z_{\rm ini}$ (in coloured lines), with respect to a reference simulation with $z_{\rm ini} = 150$. We use 1LPT, 2LPT or 3LPT initial conditions from left to right and three-, five- and seven-point gradients from top to bottom panels. We use the FFT solver in any case (see Section \ref{['sec:fft_methods']}). In each subplot the top and bottom panels show simulations which use TSC and CIC mass-assignments schemes respectively. Grey shaded area denotes a 1% discrepancy w.r.t the reference power spectrum.
• Figure 4: Same as Fig. \ref{['fig:initial_conditions_gradients']}, but varying the $N$-body solver instead of gradient order, which is kept to seven-point (except for the FULL_FFT solver which does not make use of finite gradients).
• Figure 5: Top panel: Ratio of the power spectrum for three-, five- and seven-point gradient in blue, orange and green lines respectively using multigrid, with respect to RAMSES (with AMR disabled and also using multigrid). Because RAMSES does not stop exactly at $z = 0$ we rescaled its power spectrum the linear growth factor. The grey shaded area show the $\pm 1\%$ limits. Bottom panel: Ratio of the power spectrum using the PySCo FFT_7pt solver w.r.t PySCo multigrid. The dark grey shaded area show the $10^{-5}$ limits. In any case we use a TSC algorithm.