They say "It's never too late to forget high school physics" and that's no doubt what had happened here. I've revised this post based on comments below, we should keep them there to maximize my embarrassment. ;-)
Since there's discussion below this answer which currently states
Contrary to what Uwe's answer seems to say, in neither a two- nor a three-body system is the ordinary energy or momentum of each body separately conserved. Even in the two-body problem, energy and momentum flow between two bodies via the gravitational force. So I disagree with that interpretation of the Wikipedia quote.
and I called the "momentum flow between two bodies" into question I thought I would look at this in a less sophisticated, less cerebral way, and instead turn to python.
Here's a two body orbit calculator. I used $m_1, m_2 = 0.2, 0.8$ and balanced the velocities to zero out center of mass motion.
With the parameter f set to 1.0 they are circular orbits and each body's angular momentum is constant. Since these are circles it also means that the magnitude of each of their linear momenta are constant.
With f set to 0.5 they are in elliptical orbits, and while each separate body's angular momentum rises and falls, we know that the sum $m_1 \mathbf{v_1} \times \mathbf{r_1} + m_2 \mathbf{v_2} \times \mathbf{r_2} = 0$ must hold.
Though I'm still not 100% comfortable with momentum flow between two bodies nor energy flow, it's certainly true that the linear momenta move in opposite manners in order to maintain conservation of momentum.
Likewise there is exchange between each body's kinetic energy and their shared potential energy, but I don't necessarily see energy "flowing" from one body to the other.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
def deriv(X, t):
x1, x2, v1, v2 = X.reshape((4, -1))
a1 = -(x1-x2) * m2 * (((x1-x2)**2).sum())**-1.5
a2 = -(x2-x1) * m1 * (((x2-x1)**2).sum())**-1.5
return np.hstack((v1, v2, a1, a2))
m1, m2 = 0.2, 0.8
f = 0.5
X0 = np.array([0.8, 0, -0.2, 0, 0, f*0.8, 0, -f*0.2])
times = np.arange(0, 20, 0.01)
answer, info = ODEint(deriv, X0, times, full_output=True)
x1, x2, v1, v2 = answer.T.reshape(4, 2, -1)
p1, p2 = m1 * v1, m2 * v2
L1, L2 = m1 * np.cross(x1, v1, axisa=0, axisb=0), m2 * np.cross(x2, v2, axisa=0, axisb=0)
KE1, KE2 = 0.5 * m1 * (v1**2).sum(axis=0), 0.5 * m2 * (v2**2).sum(axis=0)
PE = - m1 * m2 / np.sqrt(((x2-x1)**2).sum(axis=0))
Etot = KE1 + KE2 + PE
if True:
plt.figure()
plt.subplot(5, 1, 1)
plt.plot(x1[0], x1[1])
plt.plot(x2[0], x2[1])
plt.plot([0], [0], '.k')
plt.plot(x1[0][0], x1[1][0], 'ok')
plt.plot(x2[0][0], x2[1][0], 'ok')
xmin, xmax = plt.xlim()
plt.xlim(xmin-0.05, xmax+0.05)
ymin, ymax = plt.ylim()
plt.ylim(ymin-0.05, ymax+0.05)
plt.gca().set_aspect('equal')
plt.subplot(5, 1, 2)
for thing in (x1[0], x1[1], x2[0], x2[1]):
plt.plot(times, thing)
plt.title('x1, y1, x2, y2')
plt.subplot(5, 1, 3)
for thing in (p1[0], p1[1], p2[0], p2[1]):
plt.plot(times, thing)
plt.title('px1, py1, px2, py2')
plt.subplot(5, 1, 4)
plt.plot(times, L1)
plt.plot(times, L2)
plt.ylim(0, 0.14)
plt.title('L1, L2')
plt.subplot(5, 1, 5)
plt.plot(times, KE1)
plt.plot(times, KE2)
plt.plot(times, PE)
plt.plot(times, Etot)
plt.title('E1, E2, PE, Etot')
plt.show()
