Accuracy of power spectra in dissipationless cosmological simulations

We exploit a suite of large \emph{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. 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-3 k_\Lambda$ at the very latest times we evolve to. At the $5\%$ level, accuracy extends up to 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 extrapolation to place conservative bounds on the accuracy of \emph{N}-body simulations of non-scale free models like LCDM. We note that deviations due to discretization in the converged range are not well modelled by shot noise, and subtracting it in fact degrades accuracy. 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 remark that conclusions about small scale clustering (e.g. concerning the validity of stable clustering) obtained using PS data at wavenumbers larger than a few $k_\Lambda$ may need revision in light of our convergence analysis.