In this author's opinion, the most exciting recent prospect is the application of rescaling of time discussed above, particularly by Coomes et al. . The fact that they could shadow the Lorenz equations times longer using a rescaling of time is remarkable. Miller  has a diagram showing the distance between two nearby solutions with 8 particles, in which 2 of the particles form a temporary binary star with a highly elliptical orbit. There are very large spikes in the separation of the two solutions, representing the times of closest approach of the two stars, because at closest approach the stars speed up immensely. Thus, the phase space separation of the two nearby solutions temporarily increases dramatically. Without a re-scaling of time, it is unlikely that a shadowing algorithm would be able to shadow such an orbit, when in fact it is not unreasonable to say that a true solution that temporarily moves away from a numerical solution in this fashion can still be labelled a shadow. It is already known that glitches in N-body simulations occur most frequently near close approaches [62, 63, 33]. It will be interesting to see how often the glitches can be eliminated by a rescaling of time. It may be related to Chow & Palmer's -- perhaps p measures the number of timesteps in an encounter.
Another interesting point in regard to non-hyperbolicity is the observation that N-body systems have some purely imaginary eigenvalues (cf. footnote here ). This is not surprising, since there are components of regular motion in most N-body systems. It may be that these ``rotating'' directions can not be ``pinned to zero'' at the endpoints of the shadow. Perhaps it may be possible to handle this non-hyperbolic structure by allowing these ``rotating'' directions to have non-zero boundary conditions at the endpoints, unlike the exponential directions.
Quinlan and Tremaine [62, 63] ask the question of whether shadows, even if they exist, are typical of true orbits chosen at random. It would seem odd if they were not. Perhaps bi-shadowing may be applicable. Bi-shadowing studies not only the existence of a true solution close to a numerical solution, but the existence of possible numerical solutions close to arbitrary true solutions.
Many authors have demonstrated that the rate of exponential divergence decreases as a system becomes less collisional [47, 48, 49, 27, 63]. Hayes  demonstrates that a simple scaling of the N-body problem gives a shadow length inversely proportional to N, which seems rather disappointing. However, Quinlan and Tremaine, and Hayes, discuss the possibility that, even if shadowing becomes more difficult in higher dimensions, perhaps the fact that N-body systems become smoother and less collisional with higher N may compensate for this, allowing smooth systems with large numbers of particles to be shadowable.
There are many, many more interesting areas of further research. The author's Master's Thesis  discusses many of them. It also discusses other measures of error that may be used if shadowing turns out to be too stringent a measure of error.
Finally, there may be the possibility of proving rigorously that a shadow does not exist. The literature on this subject seems rather sparse; it is mentioned briefly in Adams et al.  and Iserles et al. .