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 . 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 .
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.'' 
Miller  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  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  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  apparently confuses (according to ) 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  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  studied many aspects of chaos in stellar N-body systems. He provides yet another derivation of , similar to Goodman et al. , 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 , 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.