Accuracy of power spectrum measurements in dissipationless cosmological simulations

We exploit a suite of large N-body simulations (up to N=\(4096^3\)) performed with Abacus, of scale-free models with a range of spectral indices \(n\), to better understand and quantify convergence of the matter power spectrum in dark matter only cosmological N-body simulations. Using self-similarity to identify converged regions, we show that the maximal wavenumber resolved at a given level of accuracy increases monotonically as a function of time. At the \(1\%\) level it starts at early times from a fraction of \(k_\Lambda\), the Nyquist wavenumber of the initial grid, and reaches at most, if the force softening is sufficiently small, \(\sim 2 k_\Lambda\) at the very latest times we evolve to. At the \(5\%\) level accuracy extends up to slightly larger wavenumbers, of order \(5k_\Lambda\) at late times. Expressed as a suitable function of the scale-factor, accuracy shows a very simple \(n\)-dependence, allowing a straightforward extrapolation to place conservative bounds on the accuracy of N-body simulations of non-scale free models like LCDM. Quantitatively our findings are broadly in line with the conservative assumptions about resolution adopted by recent studies using large cosmological simulations (e.g. Euclid Flagship) aiming to constrain the mildly non-linear regime. On the other hand, we note that studies of the matter power spectrum in the literature have often used data at larger wavenumbers, where convergence to the physical result is poor. Even qualitative conclusions about clustering at small scales, e.g concerning the validity of the stable clustering approximation, may need revision in light of our results.