  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.
