Long-term planetary integration with individual time steps

We describe an algorithm for long-term planetary orbit integrations, including the dominant post-Newtonian effects, that employs individual timesteps for each planet. The algorithm is symplectic and exhibits short-term errors that are $O(\epsilon\Omega^2\tau^2)$ where $\tau$ is the timestep, $\Omega$ is a typical orbital frequency, and $\epsilon\ll1$ is a typical planetary mass in solar units. By a special starting procedure long-term errors over an integration interval $T$ can be reduced to $O(\epsilon^2\Omega^3\tau^2T)$. A sample 0.8 Myr integration of the nine planets illustrates that Pluto can have a timestep more than 100 times Mercury's, without dominating the positional error. Our algorithm is applicable to other $N$-body systems.