N, M1, M2 and M4. We note that ${t_{\rm dyn}}\sim{G {{M_*}}^{5/2}(2|K+W|)^{-3/2}}$ in simulation N.

Following NLC06, the particles are spatially distributed according to equation (\[eqplum\]) and then randomly shifted in position (up to ${r_*}/5$ in modulus). This artificial, small-scale “noise” is introduced to enhance the phase mixing at the beginning of the collapse, because the numerical noise is small, and the velocity dispersion is zero (see also Section 4.2). As such, these fluctuations are not intended to reproduce any physical clumpiness.

All the simulations (realized with $N=10^6$ particles, and a grid with $\Nr=64$, $\Nth=16$ and $\Nph=32$) are evolved up $t=150{t_{\rm dyn}}$. In all cases the modulus of the center of mass position oscillates around zero with r.m.s $\lsim 0.1{r_*}$; similarly, the modulus of the total angular momentum oscillates around zero[^2] with r.m.s. $\lsim0.02$, in units of ${r_*}{{M_*}}{v_{\rm *n}}$ (simulations M$\kappa$ and N) and of ${r_*}{{M_*}}{v_{\rm *d}}$ (simulation D). $K+W$ in the Newtonian simulation and $W$ in the dMOND simulation are conserved to within $2\%$ and $0.6\%$, respectively. The volume-limited energy balance equation (\[eqbaltext\]) is conserved with an accuracy of $1\%$ in MOND simulations, independently of the adopted $V_0$. To estimate possible numerical effects, we reran one of the MOND collapse 