our results indicate that phase mixing is more effective in Newtonian gravity than in MOND[^3]. It is then interesting to estimate in physical units the phase-mixing timescales of MOND systems. From Table 1 it follows that ${t_{\rm *n}}\simeq4.7\,({r_*}/\kpc)^{3/2}({{M_*}}/10^{10}
\Msun)^{-1/2}\Myr=29.8\kappa^{-3/4}({{M_*}}/10^{10} \Msun)^{1/4}\Myr$ for ${a_0}=1.2\times 10^{-10} {\rm m}\,{\rm s}^{-2}$. For example, in the case of model M1, adopting ${{M_*}}=10^{12} \Msun$ (and ${r_*}=\sqrt{G{{M_*}}/{a_0}}\simeq34\kpc$), shells in phase space are still apparent after $\sim8.3\Gyr$ ($\simeq44{t_{\rm dyn}}$). Simulation M1 might also be interpreted as representing a dwarf elliptical galaxy of, say, ${{M_*}}= 10^9 \Msun$ (and ${r_*}=\sqrt{G{{M_*}}/{a_0}}\simeq
1.1\kpc$). In this case $44{t_{\rm dyn}}\sim 1.5 \Gyr$. We conclude that in some MOND systems substructures in phase space can survive for significantly long times.

In addition to the ($r$,$\vr$) diagram, another useful diagnostic to investigate phase-space properties of gravitational systems is the energy distribution $N(E)$ (i.e. the number of particles with energy per unit mass between $E$ and $E+dE$; e.g., Binney & Tremaine 1987; Trenti & Bertin 2005). Independently of the force law, the energy per unit mass of a particle orbiting at ${{\bf x}}$ with speed $v$ in a gravitational potential $\phi({{\bf x}})$ is $E=v^2/2+\phi({{\bf x}})$, and $E$ is constant if $\phi$ is time-independent. In Newtonian gravity $\phi$ is 