Here is a quick supplemental answer.
From the Neptune-Pluto Resonance subsection of the Wikipedia article Stability of the Solar System:
The Neptune–Pluto system lies in a 3:2 orbital resonance. C.J. Cohen and E.C. Hubbard at the Naval Surface Warfare Center Dahlgren Division discovered this in 1965. Although the resonance itself will remain stable in the short term, it becomes impossible to predict the position of Pluto with any degree of accuracy, as the uncertainty in the position grows by a factor e with each Lyapunov time, which for Pluto is 10–20 million years into the future. Thus, on the time scale of hundreds of millions of years Pluto's orbital phase becomes impossible to determine, even if Pluto's orbit appears to be perfectly stable on 10 MYR time scales (Ito and Tanikawa 2002, MNRAS).
The resonance they discovered means that for a moderate period of time (at least tens of millions of years) the relative motion between Neptune and Pluto will be repetitive, with $3\times T_{Neptune} \approx 2\times T_{Pluto}$. The paper cited in the other excellent answer is the report of this discovery by Cohen and Hubbard in 1965 mentioned here.
Figure 5 of that paper shows the motion of Neptune and Pluto in a rotating frame. The frame rotates with the average orbital rotation of Neptune, so over very long times you can see Neptune slowly rocks back and forth a bit:
If you like Python you can reproduce these fairly easy. The package Skyfield uses the same NASA JPL Ephemerides as the JPL Horizons site but is easier to use. You can see that 95% of this script is just making it look nicer and getting the position of the planets is just a few lines.
Here are the results for only 6,000 years, a much smaller period than shown in Cohen and Hubbard 1965, so Neptune only makes a small segment of its 25,000 cycle, near one end where it's "moving" slowly. The first set of plots are in inertial (non-rotating) J2000 ecliptic coordinates, and the second is rotating with the average orbital motion of Neptune, so that Neptune appears nearly fixed.
The first plot shows Neptune-Pluto separation versus calendar year. The minimum in this period seems to be around year -77, with a distance of 2.65 billion km or 17.73 AU.
def Rpos(pos, angle):
x, y, z = pos
ca, sa = np.cos(angle), np.sin(angle)
xr = x*ca - y*sa
yr = y*ca + x*sa
return np.vstack((xr, yr, z))
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from skyfield.api import Loader
loader = Loader('~/Documents/Skydata', verbose = True)
ts = loader.timescale()
de405 = loader('de405.bsp') # about 65 MB
de421 = loader('de421.bsp') # about 17 MB
de422 = loader('de422.bsp') # about 650 MB
de = de422
neptune = de['neptune barycenter']
pluto = de['pluto barycenter']
# years = np.arange(1600, 2201)
years = np.arange(-2999, 3001)
time = ts.utc(years, 1, 1) # January 1st of each calendar year
npos = neptune.at(time).ecliptic_position().km
ppos = pluto.at(time).ecliptic_position().km
npsep = np.sqrt(((npos-ppos)**2).sum(axis=0))
if True:
plt.figure()
plt.plot(years, npsep)
plt.show()
aukm = 149597870.700
print "minimum separation (km): ", npsep.min()
print "minimum separation (AU): ", npsep.min()/aukm
if True:
fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1)
for pos in [npos, ppos]:
x, y, z = pos
ax1.plot(x, y)
ax1.plot([0], [0], 'or')
ax1.set_xlim(-5E+09, 7E+09)
ax1.set_ylim(-5E+09, 7E+09)
ax2 = fig.add_subplot(1, 2, 2, projection='3d')
for pos in [npos, ppos]:
x, y, z = pos
ax2.plot(x, y, z, linewidth=1.0)
ax2.plot([0], [0], [0], 'or')
ax2.set_xlim(-5E+09, 7E+09)
ax2.set_ylim(-5E+09, 7E+09)
ax2.set_zlim(-6E+09, 6E+09)
ax2.view_init(elev=20., azim=-110)
plt.show()
nargperi = 276.3 * (np.pi/180.) # Neptune argument of perihelion (radians)
nT = 164.79 # Nepture orbital period (years)
nangle = np.arctan2(npos[1], npos[0])# - nargperi
nangle_mean = 2.*np.pi*np.mod(years/nT, 1.0) # this is a bit sloppy
nposr = Rpos(npos, -nangle_mean)
pposr = Rpos(ppos, -nangle_mean)
if True:
fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1)
for pos in [nposr, pposr]:
x, y, z = pos
ax1.plot(x, y)
ax1.plot([0], [0], 'or')
ax1.set_xlim(-7E+09, 7E+09)
ax1.set_ylim(-7E+09, 7E+09)
ax2 = fig.add_subplot(1, 2, 2, projection='3d')
for pos in [nposr, pposr]:
x, y, z = pos
ax2.plot(x, y, z, linewidth=1.0)
ax2.plot([0], [0], [0], 'or')
ax2.set_xlim(-7E+09, 7E+09)
ax2.set_ylim(-7E+09, 7E+09)
ax2.set_zlim(-7E+09, 7E+09)
ax2.view_init(elev=20., azim=-116)
plt.show()
