In this author's opinion, there is currently no conclusive evidence that the MMV assumption is correct. That it is generally accepted seems mostly to stem from the heuristic observation that N-body simulation results often look quite reasonable, and in fact sometimes reproduce results predicted by theory. For example, theoretically predicted phenomena in globular clusters such as mass segregation, core collapse, binary star formation rates, and relaxation times all seem to be reproduced reasonably faithfully by good N-body simulations [38]. Astronomers have also been successful in producing realistic-looking spiral structure in collisionless disks, but it is not clear that the simulations produce this structure for the same reason that real spiral galaxies do [68].
One of the strongest reasons for suggesting caution in the interpretation of results of N-body simulations is that N-body systems display sensitive dependence on initial conditions. This means that, at the particle level, two solutions whose phase-space initial conditions differ by an arbitrarily small amount will diverge away from each other exponentially with time. Since a numerical simulation constantly introduces small errors (cf. section 1) to the solution, we are guaranteed that our numerical solution diverges exponentially away from the true solution with the same initial conditions. This is the case even if quantities such as energy are conserved to arbitrary accuracy, since there is an infinite number of qualitatively different phase-space trajectories occupying the same energy hyper-surface. Indeed, this exponential magnification of small errors led one astronomer to exclaim, ``This view calls into question the very need for accurate simulation.'' [36]
Miller [58] seems to have been the first to demonstrate
that N-body systems display exponential divergence of nearby
trajectories, even though all first integrals were well conserved and
the evolution of the systems appeared plausible. The e-folding time
is defined to be the interval of time it takes for a quantity to increase
by a factor of e, Euler's constant (2.7182...). The e-folding
time 
of the exponential divergence
seems to be of order the crossing time, 
, although it increases
(i.e., the rate of exponential divergence decreases) as the system becomes
less collisional. This has been
both derived analytically
[36, 37, 27] and measured
experimentally
[36, 46, 47, 27]. Dejonghe
and Hut [21] showed that this exponential divergence
can happen even in the smallest possible case of N=3, where they
demonstrated that nearby initial conditions can be pushed apart by as
much as a factor of 
during a complex 3-body encounter,
although 
was more the norm. (Their method was to keep the
displacement in the linear regime by letting the trajectories diverge,
then ``pulling them together'' along the shortest path, then letting
them diverge some more, etc. Thus they were able to measure
divergences larger than could be represented by the floating point system.)
Lecar [53]
co-ordinated a study between many researchers, each of whom
independently computed the solution to a 25-body problem with
identical initial conditions. They measured the radius of the
body from the centre of mass, the moment of inertia,
and the average squared radial & tangential velocities. They found
that different algorithms and computers gave results in which these
measurements differed by as much as 100%. Although none of these
statistics are the same as the ones mentioned above as being reliable,
it is disturbing to note that some statistical measures appear to be
reliable while others are not.
Kandrup [46] apparently confuses (according to [27]) exponential divergence of nearby trajectories with relaxation, and provides an unconvincing argument that an inadequately defined phenomenon he calls ``mixing'' occurs on a crossing timescale, concluding that this mixing is a hitherto unknown form of relaxation.
Kandrup, Smith, and Willmes
[47, 48, 49]
performed hundreds of numerical experiments to
study the details of how exponential divergence is distributed amongst
individual particles, and how it behaves with respect to various
parameters. They found that the distribution of e-folding times
amongst particles had a long tail, indicating that the global
divergence was dominated by a small number of particles that, on
average, had a divergence rate about 4 times higher than the average,
and this distribution changed little with changing initial conditions.
(Interestingly, they did not mention which particles had the
highest divergence rate; this author presumes that they would be the
particles that underwent more or harder collisions than average.) They
found that the distribution of perturbations amongst particles 
after 4 crossing times was approximately Gaussian in 
, and that the softening parameter had little effect on these
results until it was greater than the average inter-particle
separation, at which point divergence of nearby trajectories became
weaker. They found that when they had particles of unequal mass, the
e-folding time was independent of mass. This seems
surprising, since it would be plausible to expect the e-folding times
to increase with increasing mass, because heavier particles move more
slowly and thus suffer fewer encounters, and the encounters are mostly
with less massive particles. Perhaps this is offset by the encounters
lasting longer with low velocity, and the fact that massive particles
sink to the centre of the cluster where they can encounter other
massive particles more often. Finally, Kandrup et al. postulate that
the exponential divergence occurs not only due to the magnifying effect
of collisions, but also due to the ``pulling apart'' force exerted by
the 
tidal force.
In contrast, Goodman, Heggie, and Hut [27] argue
that the exponential divergence is solely a result of encounters,
and is greatest in particles that undergo the strongest encounters.
They postulate that a completely smooth gravitational potential would
display no exponential divergence of nearby trajectories, i.e.,
. They build a very elegant model that explains the
exponential divergence based on magnification of perturbations during
encounters, and derive that is of order 
, which
is nearly independent of N for large N if softening is maintained
at the average inter-particle separation, or if softening is not used
at all. However, if N is increased with constant softening,
then 
for large N. This is partly the
basis on which they base their conclusion that 
for a smooth potential. Based on individual trajectory divergence,
they derive a bound on the average energy error of a single star
per crossing time of
, and postulate that a collisionless
simulation is ``accurate enough'' if individual particle energies are
conserved to this tolerance.
Heggie [37]
studied many aspects of chaos in stellar N-body systems. He provides
yet another derivation of 
, similar to Goodman et al.
[27], and uses numerical experiments to estimate
for large N. Since the derivation is based on
the magnification of perturbations during encounters, he concludes that
smooth potentials are good for modelling orbits, but not their
divergence. Finally, he experiments with integraton accuracy and finds
that the only important statistic that seems highly sensitive to
integration accuracy is the escape energy and number of escapers: with
lax integration tolerances, stars tend to escape the cluster more
often, and with higher energy, whereas, with high integration accuracy,
fewer stars escape and with less energy. This is reasonable, since
escapes usually happen as a result of close approaches, which can tax
the integrator and lead to unreal 2-body encounters if too low an
accuracy is used. This seems to be a strong
argument for the use of regularization in simulations of collisional
systems. Furthermore, since escape rates and escape energies are
notoriously difficult to predict theoretically [37],
it seems important that
collisions be handled correctly in order to get a reliable measures of
this phenomenon.
One surprising observation is that, to this author's knowledge, nobody has tried doing numerical experiments of exponential divergence in a cold disk, which may help us understand exponential divergence in collisionless systems.