In the last post, we talked about the stability of motion that started in the vicinity of the Lagrange points. We found that motion around any of the colinear points ($L_1$, $L_2$, or $L_3$) was unconditionally unstable and that motion around the triangular points ($L_4$ and $L_5$) was conditionally stable. In this post, we are going to talk about the kind of motion available. We are only going to touch on the basic concepts because the subject is quite involved with entire works have been devoted to the rich motions that result.

For our purposes, we are going to start with the linearized equations of motion

\[ \delta {\ddot x} – 2 \delta {\dot y} = (1+2A) \delta x \;, \]
\[ \delta {\ddot y} + 2 \delta {\dot x} = (1-A) \delta y \;, \]
and
\[\delta {\ddot z} = -A \delta z \;,\]
with
\[A = \frac{1-\mu}{d_1^3} + \frac{\mu}{d_2^3} \; .\]

Note that we are evaluating $A$ at whichever colinear point we are interested in, as if the motion were really occurring at that point. This assumption allows us to set $A$ to one specific value (which we’ll calculate for the Sun-Earth/Moon $L_2$ point below) and subsequently treat as a constant. In more careful considerations, we would evaluate $A$ at each point along the trajectory as the motion evolves. This ‘extra step’ is analytically very difficult, typically requiring Floquet or Magnus type solutions, and the preferred method is to employ numerical solutions programmed on a computer.

With that simplifying assumption understood, we start by tackling the motion in $z$ since it is not coupled to motion in the other two dimensions. The immediate solution is

\[ \delta z(t) = {\mathcal A}_{z} \cos ( \omega_z t + \phi_z ) \; \] with $\omega_z = \sqrt{A}$ . The phase $\phi_z$ and the amplitude ${\mathcal A}_z$ are arbitrary. Note that despite their formal similarity, the amplitude ${\mathcal A}_z$ is always distinguishable from the dynamic parameter $A$ by virtue of the subscript (even if the font-face difference is not always clear at a glance).

Before turning to the motion in the $x$-$y$ plane, let’s look at how to solve for the value of $A$. Since we are linearizing about a Lagrange point, we need to evaluate $A$ at that point. For the Sun-Earth/Moon $L_2$ point, we use the corresponding quintic equation
\[f(\gamma) = \gamma^5 + (3-\mu)\gamma^4 + (3-2\mu)\gamma^3 – \mu \gamma^2 – 2 \mu \gamma – \mu = 0 \; ,\]
with
\[\mu = \frac{m_{Earth}+m_{Moon}}{m_{Sun} + m_{Earth}+m_{Moon}} \approx 3.0404\times 10^{-6}\]
To find the root, we use the Newton-Raphson method, which requires us to use the derivative of the quintic, given by:
\[f'(\gamma) = 5 \gamma^4 + 4(3-\mu)\gamma^3 + 3(3-2\mu)\gamma^2 – 2\mu \gamma – 2 \mu \;\]
and then to update our current guess iteratively as

\[ x_{new} = x_{old} – \frac{f(\gamma)}{f'(\gamma)} ;. \]

Following this algorithm gives a numerical value of $A=3.8825$ with the corresponding frequency $\omega_z = 1.9704$ (in the appropriate non-dimensional units of the CR3BP).

We now turn to the in-plane motion. Start by assuming $\delta x = \alpha e^{\lambda t}$ and $\delta y = \beta e^{\lambda t}$ , where we have used the fact that the motion is coupled to have the time dependance in each of the in-plane degrees of freedom to match (i.e., $\lambda$ for both).

The coupled equations of motion take on the form of
\[ \alpha \lambda^2 – 2\ \beta \lambda = (1+2A) \alpha \]
and
\[\beta \lambda^2 + 2 \alpha \lambda = (1-A) \beta \; .\]
Next, recast these equations as matrix versions
\[ \left[ \begin{array}{cc} \lambda^2 – (1+2A) & -2 \lambda \\ 2 \lambda & \lambda^2 – (1-A) \end{array} \right] \left[ \begin{array}{c} \alpha \\ \beta \end{array} \right] = \left[ \begin{array}{c} 0 \\ 0 \end{array}\right] \; \]
and, since we want non-trivial solutions for $\alpha$ and $\beta$, we demand that the determinant of the matrix is zero. Carrying out the algebra leads to the characteristic equation with details similar to what was presented in the last post but expanded on here so that the presentation is self-containted
\[\lambda^4 + (2-A) \lambda^2 +(1+2A)(1-A) = 0 \; .\]
The solutions for $\lambda^2$ can then be written as:
\[\lambda^2 = \frac{A-2 \pm \sqrt{A(9A-8)}}{2} \; . \]
We know from our numeric work that $A>1$ so that the constant term, $c = (1+2A)(1-A)$ is negative. From this we know that the discriminant, $\sqrt{A(9A-8)}$, is positive and from Vieta’s formula’s for quadratic equations (namely $q_+ q_- = c/a$,), we conclude that the two roots for $\lambda^2$ have opposite signs and, following the notation of last post, we label these two roots as $q_+$and $q_-$. The subsequent sets of values for $\lambda$ we can categorize as
\[\lambda_{+1} = +\sqrt{q_+} \; ,\]
\[\lambda_{+2} = – \sqrt{q_+} \; ,\]

\[ \lambda_{-1} = i \sqrt{|q_-|} \; ,\]
and

\[\lambda_{-2} = – i \sqrt{|q_-|} \; .\]
The real pair of values, $\lambda_{+1}$ and $\lambda_{+2}$ lead to either exponentially growing or decreasing motion but, in either case, no oscillations and so can be ignored if we are looking for motion like JWST executes around a Lagrange point. The oscillatory motion arises from the imaginary pair with the corresponding frequency being given by

\[\omega_{xy} \equiv \sqrt{|q_-|} = \sqrt{\left|\frac{A – 2-\sqrt{A(9A-8)}}{2} \right|} \;.\]

Substituting in the value of $A$ we derived above gives
\[\omega_{xy} = 2.0427 \;. \]

Since the in-plane and out-of-plane frequencies are not equal (although they are relatively close), periodic motion is not possible and the result is a Lissajous pattern in which the motion slowly evolves toward or away from the line on which the Lagrange points lie (the Sun-Earth/Moon line in the case of JWST), depending on initial phase.

The final ingredient to model this Lissajous motion is to determine the scaling between motion along $x$ (denoted by $\alpha$ earlier) to that along $y$ (denoted by $\beta$). Either of the two dynamical equations will serve but the second one for $\delta {\ddot y}$ is slightly easier to work with. Plugging in the proposed solutions we get

\[\beta \lambda^2 + 2 \alpha \lambda = (1-A) \beta \;,\]

from which we get

\[\frac{\beta}{\alpha} = -\frac{2 \lambda}{\lambda^2 + A -1}\;.\]

Since we are focusing on oscillatory motion we will use $\lambda = \omega_{xy}$ (i.e. we throw away considerations of $\lambda_{-1}$ and $\lambda_{-2}$ being imaginary). Plugging in that value gives

\[ \frac{\beta}{\alpha} = – 3.187 \; .\]

The magnitude means that the amplitude along the $y$-axis is roughly 3 times larger than the corresponding $x$-axis and the minus sign serves only to alter the phase relationship between their motion.

A typical plot of this system is shown in the following plot, made several decades ago when I was the flight dynamics lead for the WMAP mission.

The RLP abbreviation stands for Rotating Libration Point and is the name NASA gave to the rotating coordinate system (at least at that time – names and fads change) for the rotating coordinate system in which the equations of motion are defined. This plot shows how motion around the Sun-Earth/Moon $L_2$ produces a Lissajous pattern. Each revolution around the Lagrange point (the origin) takes approximately 180 days in keeping with the corresponding frequencies being approximately 2.

We should close with two points in brief. First, since the equations of motion are non-linear, it is possible (and useful) to design motion that is periodic (rather than the quasi-periodic motion of the Lissajous) by increasing the amplitude of the $x$-$y$ and $z$ motions. At some point, roughly 10 times the values shown in the plot above, the in-plane and out-of-plane frequencies equal and it is possible to have periodic motion, which is often called a halo orbit. Second, one might wonder how close the actual solar system model can come to what is clearly an idealization. The answer is quite close as the following side-by-side comparison, created during my time on JWST shows.

In this plot, the Earth is the white pixel at the center and the green trail is the trajectory from launch onwards. The ‘squiggly’ segment between the Earth and the periodic halo is the ascent trajectory.

In this post we are going to look at motion about the equilibrium points.  From here forward, I will refer to these points as the Lagrange points both in honor of his achievement and to conform with common usage.  It is this aspect of the circular restricted three body problem – that an object can ‘orbit’ about a fictitious point in space – that fascinates many people and enables missions like ACE, WIND, WMAP, and, especially, JWST.  I personally designed Lagrange point orbits for ACE, WMAP, and JWST and had the pleasure of having each of those missions to ‘fly’ these trajectories and the transfers to them.

In some of the earlier posts in this series, we already established that there are 5 Lagrange points where an object placed at rest at their precise location will remain there since all the forces on the object (gravitational attraction by the massives as well as the centrifugal force in the rotating frame) negate each other.  The question is whether motion in the vicinity can be established.  In other words, can deviations between two trajectories, defined as

\[ \delta {\vec r} (t) = {\vec r}’(t) – {\vec r}(t) \; , \]

which we envision as starting small, remain so as time progresses or will they continue to deviate further and further in an unbounded way.  For our purposes, the former case will be dubbed stable while the later will be unstable.  This usage is consistent with, although vaguer than, the more precise modern definitions which require much more machinery.

To answer this stability question, we start with the equation of motion for the deviation $\delta {\vec r}(t)$, given in component form as

\[ \delta {\ddot x} – 2 \delta {\dot y} = \partial_x U({\vec r}’) – \partial_x U({\vec r}) \; , \]

\[ \delta {\ddot y} + 2 \delta {\dot x} = \partial_y U({\vec r}’) – \partial_y U({\vec r}) \; , \]

and

\[ \delta {\ddot z} = \partial_z U({\vec r}’) – \partial_z U({\vec r}) \; , \]

where

\[ U = \frac{1}{2} (x^2 + y^2) + \frac{\mu}{\sqrt{y^{2} + z^{2} + \left(\mu + x – 1\right)^{2}}} + \frac{1 – \mu}{\sqrt{y^{2} + z^{2} + \left(\mu + x\right)^{2}}} \; .\]

At this point, the equations are exact but not very helpful since the right-hand sides involve differences in the derivatives of the pseudopotential evaluated at different points in space.  We eliminate that awkwardness by linearizing the right-hand side about ${\vec r}$ so that we can express derivatives of $U(\vec r’)$ with respective to derivatives of $U(\vec r)$ using the standard Taylor series expansion.  To keep this expansion readable, we’ll tweak notation to use the comma abbreviation for partial differentiation (e.g., $U_{,x} \equiv \partial_x U$) to get

\[ \left[ \begin{array}{c} \delta {\ddot x} – 2\delta {\dot y} \\ \delta {\ddot y} + 2 \delta {\dot x} \\ \delta {\ddot z} \end{array} \right] = \left[ \begin{array}{cc} U_{,xx} & U_{,xy} & U_{,xz} \\ U_{,yx} & U_{,yy} & U_{,yz} \\ U_{,zx} & U_{,zy} & U_{,zz} \end{array} \right] \left[ \begin{array}{c} \delta x \\ \delta y \\ \delta z \end{array} \right] \; . \]

The evaluation of the derivatives are a bit painful to do by hand but straightforward using Sympy.  In symbolic form, the six independent second derivatives are:

\[ U_{,x x} = \frac{\mu \left(- 3 \mu – 3 x + 3\right) \left(- \mu – x + 1\right)}{\left(y^{2} + z^{2} + \left(\mu + x – 1\right)^{2}\right)^{\frac{5}{2}}} \\ – \frac{\mu}{\left(y^{2} + z^{2} + \left(\mu + x – 1\right)^{2}\right)^{\frac{3}{2}}} \\ + \frac{\left(1 – \mu\right) \left(- 3 \mu – 3 x\right) \left(- \mu – x\right)}{\left(y^{2} + z^{2} + \left(\mu + x\right)^{2}\right)^{\frac{5}{2}}} \\ – \frac{1 – \mu}{\left(y^{2} + z^{2} + \left(\mu + x\right)^{2}\right)^{\frac{3}{2}}} + 1 \; , \]

\[ U_{,x y} = \, – \frac{3 \mu y \left(- \mu – x + 1\right)}{\left(y^{2} + z^{2} + \left(\mu + x – 1\right)^{2}\right)^{\frac{5}{2}}} – \frac{3 y \left(1 – \mu\right) \left(- \mu – x\right)}{\left(y^{2} + z^{2} + \left(\mu + x\right)^{2}\right)^{\frac{5}{2}}} \; , \]

\[ U_{,x z} = \,- \frac{3 \mu z \left(- \mu – x + 1\right)}{\left(y^{2} + z^{2} + \left(\mu + x – 1\right)^{2}\right)^{\frac{5}{2}}} – \frac{3 z \left(1 – \mu\right) \left(- \mu – x\right)}{\left(y^{2} + z^{2} + \left(\mu + x\right)^{2}\right)^{\frac{5}{2}}} \; , \]

\[ U_{,y y} = \frac{3 \mu y^{2}}{\left(y^{2} + z^{2} + \left(\mu + x – 1\right)^{2}\right)^{\frac{5}{2}}} \\ – \frac{\mu}{\left(y^{2} + z^{2} \\ + \left(\mu + x – 1\right)^{2}\right)^{\frac{3}{2}}} \\+ \frac{3 y^{2} \left(1 – \mu\right)}{\left(y^{2} + z^{2} + \left(\mu + x\right)^{2}\right)^{\frac{5}{2}}} \\- \frac{1 – \mu}{\left(y^{2} + z^{2} + \left(\mu + x\right)^{2}\right)^{\frac{3}{2}}} + 1 \; , \]

\[ U_{,y z} = \frac{3 \mu y z}{\left(y^{2} + z^{2} + \left(\mu + x – 1\right)^{2}\right)^{\frac{5}{2}}} + \frac{3 y z \left(1 – \mu\right)}{\left(y^{2} + z^{2} + \left(\mu + x\right)^{2}\right)^{\frac{5}{2}}} \; , \]

and

\[ U_{,z z} = \frac{3 \mu z^{2}}{\left(y^{2} + z^{2} \\ + \left(\mu + x – 1\right)^{2}\right)^{\frac{5}{2}}} \\- \frac{\mu}{\left(y^{2} + z^{2} + \left(\mu + x – 1\right)^{2}\right)^{\frac{3}{2}}} \\+ \frac{3 z^{2} \left(1 – \mu\right)}{\left(y^{2} + z^{2} + \left(\mu + x\right)^{2}\right)^{\frac{5}{2}}} \\- \frac{1 – \mu}{\left(y^{2} + z^{2} \\+ \left(\mu + x\right)^{2}\right)^{\frac{3}{2}}} \; . \]

These forms would be necessary if we were to linearize about a reference trajectory as is often done in cutting edge studies but here we will content ourselves with evaluating these first at the colinear points and then at the equilateral points as a first step to understanding the stability of motion around the Lagrange points.

Colinear points

The trickiest part of evaluating the colinear points is that we don’t know the value of $x^*$ of any of the three points analytically even though we know that $y^*=0$ and $z^*=0$.  As was discussed in a previous post (Part 5: Finding the Equilibrium Points) the value of $x^*$ must come from solving a non-trivial quintic equation numerically.  Nonetheless, we can instruct Sympy to do much of the heavy lifting for us as follows.

With the data structure that contains the second partial derivatives (called ddU), we can use a substitution dictionary

to enter what we do know.  An important step before doing this is to confirm that the same substitution dictionary gives that the first partials evaluate to zero either directly (in the case of $U_{,y}$ and $U_{,z}$ or that we get the very quintic equation whose root we are using back in the case of $U_{,x}$.  Once a successful confirmation is obtained, we can then use it in the following snippet

to get the only non-zero derivatives as

\[ U_{,xx} = \frac{2 \mu}{r_{2}^{3}} + 1 + \frac{2 – 2 \mu}{r_{1}^{3}} \; , \]

\[ U_{,yy} = \frac{\mu}{r_{2}^{3}} + 1 + \frac{\mu – 1}{r_{1}^{3}} \; , \]

and

\[ U_{,zz} = – \frac{\mu}{r_{2}^{3}} + \frac{\mu – 1}{r_{1}^{3}} \; , \]

where

\[ r_1 \equiv \left| \mu + x^* \right| \; ,\]

\[ r_2 \equiv \left| \mu + x^* – 1\right| \; , \]

and $x^*$ is the value we get for the location of any of the colinear points (again see Part 5: Finding the Equilibrium Points for details).

Following Roy, we recognize that we can write these in an even more compact form and substitute them into the equations of motion for the deviation vector to get

\[ \delta {\ddot x} – 2 \delta{\dot y} = (1 + 2A) \delta x \; ,\]

\[ \delta {\ddot y} + 2 \delta{\dot x} = (1 – A) \delta y \; ,\]

and

\[ \delta {\ddot z} = – A \delta z \; ,\]

where

\[ A = \frac{1-\mu}{r_1^3} + \frac{\mu}{r_2^3} \; \]

is a real positive quantity with three choices depending on which of the colinear points we are at. 

The $z$ equation of motion is easy to solve, giving oscillatory behavior in that direction and since oscillatory behavior is bounded, we can conclude the motion to be stable.

Motion in the $x$ and $y$ directions are coupled and we have to examine the characteristic polynomial associated with the corresponding differential equations.  Defining a variational state as

\[ \delta {\bar s} = \left[ \begin{array}{c} \delta x \\ \delta y \\ \delta {\dot x} \\ \delta {\dot y} \end{array} \right] \; ,\]

the coupled $x$-$y$ motion takes on the form

\[ \frac{d}{dt} \delta {\bar s} = \left[ \begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\(1+2A) &  0 &  0 & 2 \\ 0 & (1-A) & -2 & 0 \end{array} \right] \delta {\bar s} \equiv M \delta {\bar s}\]

The characteristic polynomial, which comes from the determinant of $M – \lambda I_{4 \times 4}$, is

\[ \lambda^4 +(2-A) \lambda^2 + (1+ A -2 A^{2}) = 0 \; .\]

Roy states that

Using that fact and calling $c = 1 + A – 2A^2$, we know that the discriminant, $b^2 – 4ac$ (with $b=(2-A)$ and $a=1$) is positive and that the roots for $\lambda^2$ must be real numbers that we can call $q_{+}$ and $q_{-}$.  Also, since the product of these roots must equal $c/2a$, one of them must be positive and the other negative.  From this analysis we can conclude that of the four roots of $\lambda$, two of them must be real numbers with one positive and the other negative (corresponding to the $\pm$ choice in taking the square root) while the other two are purely imaginary.  The imaginary roots give oscillatory motion and the negative root gives exponentially decaying motion (in forward time) and so are stable.  The positive root give exponentially growing motion for the deviation and we are forced to conclude that the motion about the colinear Lagrange points are unstable.  We’ll defer the practical implications of that result for missions like JWST until the next post.

The equilateral Lagrange points are much easier to analyze.  The following Sympy snippet

The non-zero results are:

\[ U_{,xx} = 3/4 \; , \]

\[ U_{,xy} = \frac{3 \sqrt{3}}{4}(1 – 2 \mu) \; ,\]

\[ U_{,yx} = \frac{3 \sqrt{3}}{4}(1 – 2 \mu) \; ,\]

\[ U_{,yy} = 9/4 \; ,\]

and

\[ U_{,zz} = -1 \; ,\]

Following exactly the same procedure as above gives a characteristic polynomial of

\[ \lambda^4 + \lambda^2 + \frac{27}{4} \mu (1-\mu) = 0 \; . \]

If we want stable motion, we must have that the discriminant is zero or negative which occurs when

\[ 1 – 27\mu(1-\mu) \ge 0 \; .\]

Roy rewrites this inequality in terms of a small positive quantity $\epsilon$ and solves it for $\mu$ to get

\[ \mu = \frac{1}{2} \pm \sqrt{ \frac{23 + 4 \epsilon}{108} } \; .\]

He then argues that since $\mu \le 1/2$, we must choose the negative sign and that in the limit as $\epsilon$ approaches zero from above that $\mu < 0.0385$.  This is the condition for stability and it is meet in almost all cases in the solar system of interest (e.g., Sun-Jupiter, Sun-Earth/Moon, Earth-Moon, etc.).  Indeed, the existence of the Trojan asteroids that ‘hang’ at the Sun-Jupiter L4 (and L5) points are a strong confirmation of Lagrange’s original analysis.  Finally, the same analysis holds for L5 since the only thing that changes are the signs of $U_{,xy}$ and $U_{,yx}$ (which, of course, equal each other).  That sign change is ‘squared’ out when deriving the characteristic polynomial and so, except for the lagging versus leading location, L5 behaves like L4.

Up to this point we’ve looked at a variety of methods for understanding how motion can occur in the circular restricted three body problem.  At one extreme, the existence of the Jacobi constant and its significance, in combination with the pseudopotential, allowed us to make quasi-global statements about the spatial extents of allowed motion without saying anything, per se, about the trajectories themselves.  At the other extreme, the determination of the equilibrium points gives us some idea about very local structures and the trivial trajectories of the restricted bodies placed there.  But in both cases, we haven’t said anything about specific, non-trivial trajectories that can be realized with the CR3BP.  In this post and several ones hereafter we will address some of this topic. 

Unlike the Kepler problem, where we know from analysis that the allowed solutions must be one of the four conic sections, the CR3BP admits quite rich behaviors and it is impossible to catalog all the possibilities.  That said, one of the time-honored approaches is to look for periodic structures, sometimes called orbits or trajectories, interchangeably, within the CR3BP.  As Roy describes it:

One note of clarification is in order.  By ‘Poincaré’s conjecture’, Roy does not mean the problem from topology that is famous now (see Poincaré conjecture – Wikipedia) but rather the recurrence theorem of Poincaré from which he conjectures that phase space is, in some way, dense in periodic orbits.

We are assisted significantly in the hunt for periodic orbits by a dynamical symmetry exhibited by the equations of motion.  It is that symmetry that will concern us for the rest of this post.

To begin, let’s remind ourselves of the structure of the equations of motion which is:

\[ {\ddot x} – 2 {\dot y} = U_{,x} \; , \]

\[ {\ddot y} + 2 {\dot x} = U_{,y} \; , \]

and

\[ {\ddot z}  = U_{,z} \; . \]

The pseudopotential $U(x,y,z)$ is defined as

\[ U = \frac{1}{2} \left( x^2  + y^2 \right) \, – \frac{1-\mu}{d_1} \, – \frac{\mu}{d_2} \; , \]

where

\[ d_1 = \sqrt{ (x+\mu)^2 + y^2 + z^2 } \; \]

and

\[ d_2 = \sqrt{ (x -1 +\mu)^2 + y^2 + z^2 } \; . \]

By their very nature, periodic orbits look the same whether time is running forward or running backward.  Note that there is a philosophical point in the previous sentence that we will return to later.  So, we want to snoop out a symmetry operation that we can perform that leaves the equations of motion unchanged when $t \rightarrow -t$.

Under time-reversal, the position components remain the same, as one might expect or even insist upon since time concerns motion not relative positioning.  As a consequence, the pseudopotential, which only depends on the relative positions of the massive and the restricted body also remains unchanged.  The velocities are all reversed in sign, since $dt \rightarrow -dt$ while, by the same argument, the accelerations remain unchanged as $dt^2 \rightarrow (-dt)^2 = dt^2$.  The equations of motion are now subtly different only in the second terms on the left-hand side:

\[ {\ddot x} + 2{\dot y} = U_{,x} \; , \]

\[ {\ddot y} – 2{\dot x} = U_{,y} \; , \]

and

\[ {\ddot z}  = U_{,z} \; . \]

How does one fix the equations of motion to restore their form?  One might try making the spatial transformation of $x \rightarrow -x$ since:

\[ U_{,x} \rightarrow U_{,-x} = -U_{,x} \; , \]

\[ {\dot x} \rightarrow -{\dot x} \; , \]

and

\[ {\ddot x} \rightarrow -{\ddot x}  \; .\]

This transformation would result in the first equation having an overall minus sign, which cancels out, and the second equation having the Coriolis term flipping its sign back to the original signature.  However, both $d_1$ and $d_2$, would be different since:

\[ d_1 = \sqrt{ (x+\mu)^2 + y^2 + z^2 } \rightarrow \sqrt{ (-x+\mu)^2 + y^2 +z^2 } \; \]

and

\[ d_2 = \sqrt{ (x-1+\mu)^2 + y^2 + z^2 } \rightarrow \sqrt{ (-x-1+\mu)^2 + y^2 +z^2 } \; .\]

But, a moments reflection should convince us that the idea would work if the coordinate that was ‘reflected’ was $y$ instead of $x$ since $d_1$ and $d_2$ involve only $y^2$.

This dynamical symmetry, $t \rightarrow -t$ and $y \rightarrow -y$, is a special case of the much more general mirror theorem from $N$-body dynamics (see Roy’s Orbital Motion 3^rd Edition, Sec. 5.6). 

One implication of the mirror theorem in the CR3BP is that periodic orbits must come into the $x$-$z$ plane with a velocity perpendicular to the normal to the plane; in other words, at the plane crossing the motion must be only along the $y$ axis.

This observation allows for an entire apparatus to be built for generating periodic orbits by targeting some set of the initial conditions such that, at each $x$-$z$ plane crossing the trajectory is perpendicular.

While this hunt for periodic orbits is a fruitful pursuit for academic research, it often obscures the proper way to hunt for motion about the equilibrium points in a real problem (often referred to as the ephemeris problem).  There are two reasons for this that are interrelated but different.  First, and least important, is the fact that the ephemeris problem is only a broken shadow of the CR3BP.  There are more than two massive bodies in the solar system and none are in circular orbit around a common barycenter.  Second, and far more important philosophically, is that periodic orbits, by construction, cannot be entered or left.  A periodic orbit existed unaltered from the distant past and will exist unaltered into the distant future and, so, even if they are ‘close’ in some sense to all the other trajectories there is a limit to their use.  The concept of a periodic orbit serves the hunt for actual motion around the Lagrange point much in the way that broken symmetry serves in field theories – a guide for organizing the analysis but one must always remember it to be a guide.

This is a particularly important frame of mind when using orbital maneuvers to fashion a desired trajectory.  It is possible to jump onto a periodic orbit by the use of maneuvers, which serve to change the dynamical problem from one without man-made controls, where the desired trajectory is periodic, to one with controls where it is not periodic.  Using orbital maneuvers that target on a perpendicular plane crossing is a waste of time and is a fragile process since there are no real periodic orbits with or without maneuvers.  Undue focus on the mirror theorem makes one lose sight of it as a guide.

An example from my trajectory design work on JWST illustrates this point.  Consider the following figure.

Panels (a) and (b) show two separate periodic orbits in the CR3BP.  Being periodic means that despite their closeness, as shown in panel (c), they do not cross in phase space and so a transfer of the restricted body from red to black or vice versa is forbidden.  It is possible to ‘force’ the transfer by way of a maneuver somewhere in the vicinity of where the red trajectory, after it ‘loops back’, begins to cross the dotted line (which is the line joining the two massives in the rotating frame and so contains the $x$ axis).  However, in the ephemeris model, the elliptic motion of the Earth/Moon about the solar system barycenter allows to a maneuver free transition between red and black.  For programmatic reasons that had to do with the hardware JWST flew, the observatory did not take advantage of this dynamical ‘freebie’ but it was there for the taking nonetheless.  In my design work to make the JWST transfer work, I used the mirror theorem as a guide but I didn’t use the targeting technique suggested by it.  Instead, I used one of my own devising which allowed for incredibly rapid generation of trajectories – a feature that was particularly important given the number of times JWST’s launch date slipped.

In the last two posts, we established the general structure of the pseudopotential and the fact that the allowed motions are restricted to certain spatial regions depending on the value of the Jacobi constant $C_J$.  In addition, we qualitatively argued for the existence of five equilibrium points without giving a hard-and-fast proof that they exist or a constructive method for finding them.  While having a heuristic understanding of why the points arise is important, ultimately not knowing how to find them is a serious gap.  This post remedies that deficiency by showing definitively how to find these equilibrium points.

The condition of equilibrium means that if a test mass is placed at such a point with zero velocity (i.e., ${\dot x} = 0$, ${\dot y} = 0$, and ${\dot z}=0$) then it must remain there.  This condition, in turn, implies that there are no forces on the test mass and that ${\ddot x} = 0$, ${\ddot y} = 0$, and ${\ddot z} = 0$. 

Substituting these requirements into the equations of motion yields:

\[ \frac{\partial U}{\partial x} \equiv U_{,x} = 0 \; , \]

\[ \frac{\partial U}{\partial y} \equiv U_{,y} = 0 \; , \]

and

\[ \frac{\partial U}{\partial z} \equiv U_{,z} = 0 \; . \]

Thus, the process of finding the equilibrium points within a dynamical system transforms into a geometry problem of finding the critical points of the function $U(x,y,z)$ of three variables.

The first equation to start with is $U_{,z} = 0$, which, when written in full, becomes

\[ z \left[ \frac{1-\mu}{d_1^3} + \frac{\mu}{d_2^3} \right] = 0 \; .\]

Since $d_1$, $d_2$, and $\mu \leq \frac{1}{2}$, are all positive quantities.  Note that we ignore the important but almost impossible to achieve condition that either $d_1$ or $d_2$ equal zero because when such a ‘collision’ occurs, the equations of motion are no longer well-defined.  Of course, this is only a mathematically possibility since a real particle will encounter a planetary surface or stellar atmosphere long before getting to the ‘singularity’ of $d_1=0$ or $d_2 = 0$.  Thus, we can conclude that $z = 0$, a result consistent with our qualitative argument for placing the equilibrium points in the $x$-$y$ plane.

We now turn to the in-plane equations starting with the $y$ equation first to get:

\[ y \left[ 1 – \frac{1-\mu}{d_1^3} – \frac{\mu}{d_2^3} \right] = 0 \; . \]

There are two choices: the quantity within the braces $[\cdots] = 0$ or $y=0$.  Let’s examine them in order.

Next, let’s factor the $x$ equation to get

\[  x \left[ 1 – \frac{1-\mu}{d_1^3} – \frac{\mu}{d_2^3} \right] – \mu \frac{1-\mu}{d_1^3} + \mu \frac{1-\mu}{d_2^3} = 0 \; . \]

By our assumptions, the first set of terms are equal to zero leaving us to conclude

\[ \frac{1}{d_1^3} = \frac{1}{d_2^3} \equiv \frac{1}{d^3} ; .\]

Returning to the $y$-equation and using this result in the bracketed-term, we find

\[ \left( 1  – \frac{1-\mu}{d^3} – \frac{\mu}{d^3} \right) = \left( 1 – \frac{1-\mu+\mu}{d^3} \right) = 1 – \frac{1}{d^3} = 0 \; . \]

The only way to satisfy this condition is to set $d = 1$, which, when used in the definition of $d_1$ gives

\[ (x + \mu)^2 + y^2 = 1 \; \]

Likewise, the definition of $d_2$ gives

\[ (x+\mu-1)^2 + y^2 = 1 \; .\]

Combining gives

\[ (x+\mu)^2 = (x + \mu- 1)^2 \; , \]

which only has a consistent solution if, when we take the square root of both sides, we pick the positive root on one and the negative root on the other.  Doing so immediately gives

\[ x = \frac{1}{2} – \mu \; \]

and

\[ y = \pm \frac{\sqrt{3}}{2} \; .\]

These are the equilateral points L4 and L5 whose locations are off of the line joining the massives forming a equilateral triangle.

The remaining three points come from going back to the choice we had to make after requiring $z=0$ and this time requiring $y=0$ as well.  All we now have to work with is the master $x$-equation that reads

\[ x – \frac{(1-\mu)(x+\mu)}{|x+\mu|^3} + \frac{\mu(1-\mu-x)}{|x+\mu-1|^3} \; ,\]

where $d_1$ simplifies to $|x+\mu|$ and $d_2$ to $|x+\mu-1|$.   The master equation is actually three-fold based on the value of $x$; we get three different, specific equations for the three different regions shown in the figure below

and defined by:

• Region I: $ -\mu \leq x \leq 1 – \mu$
• Region II: $x \gt 1-\mu$
• Region III:  $x < -\mu$

In each region, the resulting equation is a fifth-order or quintic equation with no exact solution and only one real root.  There is no universal agreement about the numbering/labeling of these regions but the convention I am following is the single most common one being followed today by most (if not all) space agencies.  Since JWST flies in Region II, we will focus on the solution in that region.  The other two regions follow in similar manner.

In Region II, $x > 1 – \mu$ and so $|x+\mu| > 0$ and $|x-1+\mu| >0$.  It is convenient to define $\gamma \equiv x -1+\mu$, which is the distance from the secondary to the Region II equilibrium point (hereafter labeled L₂).  In terms of $\gamma$, the master equation in Region II becomes:

\[ \gamma + 1 – \mu \frac{1-\mu}{(\gamma+1)^2} – \frac{\mu}{\gamma^2} = 0 \; .\]

The next step, which is purely for notational convenience, is to define $\beta = 1-\mu$.  Next multiply the equation by $\gamma^2$ and $(\gamma+1)^2$ to arrive at

\[ \gamma (\gamma + 1)^2 \gamma^2 + \beta (\gamma+1)^2 \gamma^2 – \beta \gamma^2 – \mu (\gamma+1)^2 \; .\]

The algebra is manageable if we then add and subtract $2\beta \gamma + \beta$ allowing us to complete a square and then factor out $(\gamma + 1)^2$ from most of the terms.  These steps leave us at

\[ (\gamma^3 + \beta \gamma^2 – \beta – \mu)(\gamma + 1)^2 + \beta(2\gamma + 1) = 0 \; . \]

In the final step, expand and collect terms and then back-substitute to eliminate $\beta$ to get

\[ \gamma^5 + (3-\mu) \gamma^4 + (3-2\mu)\gamma^3 – \mu \gamma^2 -2\mu \gamma  -\mu = 0 \; .\]

Descartes’ rule of signs dictates that there is only one real root because all the coefficients above $\gamma^2$ are positive (recall $\mu \leq 1/2$) and all the remaining coefficients are negative.  For the Sun-Earth/Moon barycenter, $\mu \approx 3.04\times 10^{-6}$.  A plot of the polynomial above is

and the approximate value of the root is $\gamma_{L_2} = 0.01$.  Multiplying by 1AU yields that the L2 point is approximately 1.5 million kilometers from Earth, which is the well-known value.

A final note:  There are no instances of circular motion to be found in the solar system, so the actual motion, while close to the CRTBP, is really the elliptic restricted problem.  Because of that, it is common in most pieces of software to solve the quintic equation at each time step based on the instantaneous distance between the massives.  This leads to a ‘moving equilibrium’ point which may seem to be a contradiction but which, in practice, is fine as the motion of massives about their barycenter is small in its eccentricity.  The number of studies of the general elliptic, restricted three body problem is far more limited due to the fact that is has little practical application.  Nonetheless, the idea is dynamically interesting and may be explored in a future blog if space (and interest) permit.

Last month we made a preliminary exploration of the pseudopotential, mapping out the various large-scale features and discovering those characteristics that suggest the basic structure of the 5 equilibrium points.  The argument for their existence and their location was made heuristically by appealing to a simple physical problem (two fixed, co-equal masses) that was then ‘relaxed’ into a slightly more complex problem (two co-orbiting, co-equal masses).

This argument is compelling (and correct) and it provides the best guidance for understanding and exploring practical problem.  However, it lacks the quantitative aspects required to constructively find the equilibrium points and understand their dynamic stability.  Ordinarily this might be viewed as a deficiency that should be immediately remedied (and, in fact, it will be attended to in a sequel) but there is actually a lot more information that can be wrung from this heuristic viewpoint before moving detailed work.

The reason for lingering over these qualitative aspects is that the CRTBP cannot be solved analytically in any meaningful way.  Therefore, having qualitative limits on the possible motions is far more important than the usual ones available in the Kepler problem.

For example, energy conservation within the Kepler problem allows us to know that the semi-major axis and the period of resulting elliptical motion is fixed from orbit to orbit.  But, this result is not nearly as useful as knowing explicitly that the semi-major axis is given by

\[ a = \frac{ – \mu }{ {\mathcal E} } \;  , \]

where ${\mathcal E}$ is the specific energy of motion.  In addition, knowing that energy is conserved is not nearly as the vis viva formula for the specific energy

\[ {\mathcal E} = \frac{1}{2} v^2 – \frac{G m}{r} \; , \]

which is the sum of specific kinetic and potential energies.  By setting the specific energy to be zero we immediately arrive at the formula for the escape speed from any object, possessing  mass $m$, as

\[ v = \sqrt{ \frac{ 2 G m}{r} } \; .\]

These simple relations tell us quite a great deal about the allowed motions in the Kepler problem. 

Unfortunately, we can’t make a similar argument in the CRTBP since that dynamical system doesn’t possess a conservation of energy.  But it does possess a conserved Jacobi constant and an analogous argument can be made, although it lacks a direct one-to-one interpretation in an important way.

We start with the analog of vis viva in the CRTBP given by the (hopefully) now familiar statement that

\[ v^2 = 2 U – C_J \; . \]

Since the velocity ${\vec v}$ is expressed in the rotating frame, it isn’t directly related to the kinetic energy.  Nonetheless, wherever its magnitude becomes zero must be a limit on the spatial excursion of the solutions.  By plotting the contours of $U(x,y,z)$ for various values of $C_J$, we can get a visual identification of those regions around the massives where the restricted particle is free to move and of those regions forbidden to it.  In the following discussion it is important to remember that, because of the traditional sign convention used in its definition, larger positive values for the Jacobi constant mean that the motion is more tightly bound to one of the two massives. 

Obviously in making the requisite contour plots we can use any value of for the mass fraction $\mu$, but, to be concrete, we will be using $\mu = 0.05$ since this value is small enough so that we see behavior similar to what we might see in the solar system but large enough not to pose a challenge in making the plots.  For comparison, the Sun-Jupiter mass fraction is around 0.001 while the Earth-Moon mass fraction is approximately 0.01. 

The following figure shows what happens as we lower the value of $C_J$ over the range of values from 6.0 down to 2.9.  These values, which are irregularly spaced, were chosen in order to show some prominent features.  Regions shaded green are those regions in the $x$-$y$ plane where motion is allowed or, in other words, regions where $U \ge \frac{1}{2} C_J$; gray-shaded regions are forbidden.

For the initial value of $C_J = 6.0$, we see that particle motion is confined to the vicinity of either the primary (larger, circular region on the left) or to the secondary (smaller circular, region on the right).  As the value of $C_J$ shrinks (or, alternatively, as the particle’s motion gets more ‘energetic’), the regions deform and significant bulges appear on the line joining the massives.  In addition, regions outside of the massive become accessible as seen by the green corner regions ‘peaking’ into the individual subplots.  These regions grow as $C_J$ continues to decrease.

At roughly $C_J = 3.42$ the bulges along the line joining the massives just touch and begin to coalesce.  At this point, the restricted body can move from the vicinity of the primary to the secondary and vice versa. 

By the time the value of $C_J = 3.3$, a ‘throat’ or ‘gate’ has opened joining the previously inner region around the massives with the outer region.  With this configuration, even more exotic motion is possible with the restricted particle being free to move from the vicinity of either the primary or secondary to regions far away and with the ability to return.  Of course, nothing in these plots indicate the time required to make the journey – simply that it is possible.

Finally, as the value of $C_J$ gets even smaller, only two regions, roughly near the triangular equilibrium points, remain as forbidden.  They persist as forbidden regions until $C_J$ become small enough that they, too, vanish.

While the analogy between $C_J$ and energy is valid, the topological changes that come about indicate that there is something more than strict level surfaces of an energy-like parameter at work.  As subsequent posts will show, there is also a component of angular momentum that features in the CRTBP problem in a significant way, partly due to the analysis being done in the rotating frame and partly as a mirror of the importance of angular momentum in the two-body problem.  But that is a story for another post.

Last month we derived Jacobi’s pseudopotential $U$ for the circular restricted three body problem.  In the co-rotating frame, the pseudopotential takes the form:

\[ U = \frac{1}{2} \left( x^2 + y^2 \right) – \frac{1-\mu}{d_1} – \frac{\mu}{d_2}  \; ,  \]

where $d_1$ and $d_2$ are the distances between the restricted body and the location of the primary and secondary bodies, respectively, which are given by

\[ d_1 = \sqrt{ (x+\mu)^2 + y^2 + z^2 }  \; \]

and

\[ d_2 = \sqrt{ (x-1+\mu)^2 + y^2 + z^2 } \; . \]

Despite the fact that the underlying equations of motion are non-autonomous, there is a conserved quantity, called the Jacobi constant $C_J$, that exists in the co-rotating frame such that

\[v^2 = 2U  – C_J\; .\]

This very innocent looking expression holds a number of important results if we know how to mine it. 

We’ll start our interrogation with a visual inspection of the magnitude of $U(x,y,0)$, which corresponds to the value of the potential in the orbital plane of the massives.  We will also set the mass fraction $\mu=0.5$ meaning that the primary and secondary are of equal mass and that the barycenter or center-of-mass around which they orbit is halfway along on the line joining them from either one.

There are a number of things to note about this surface plot.  Where the distance from the origin is large the edges of the plot are curling upwards.  This is because at large values of $x$ or $y$ the quadratic term in the pseudopotential dominates and far from the origin, the surface will look, if one doesn’t squint too hard, like a paraboloid of revolution.  The behavior is a manifestation of centrifugal force which arises in the rotating frame.  Near each of the massives, the pseudopotential spikes sharply.  For the purposes of the plot, the maximum value has been capped at a value of 4.5 but the function itself is unbounded the closer the restricted particle gets to the center of the either massive. 

Due to the definition of the Jacobi constant $C_J = 2U – v^2$, and somewhat contrary to a conventional surface in the presence of uniform gravity, upward curvature in the surface plot of $U$ represents an instability.  This observation is particularly relevant along the line joining the massives.  Careful inspection will show that point half-way between them is a saddle point with upward curvature along the $x$-direction and downward curvature along $y$.  After placing the restricted particle at the origin, we can ask what will happen if that particle is displaced a bit in one direction or the other.  Displacement by a small amount in $x$ causes the restricted particle to fall towards whichever of the two massives is closer.  Displacement by a small amount in $y$ causes the restricted particle to oscillate back and forth through the origin. 

A better understanding of this dynamic behavior comes from first stepping back to a much simpler problem.  Imagine that we start with two equal masses that are fixed in space rather than being free to co-orbit their center-of-mass.  Next, let’s place our restricted body into this scenario in the region between the two bodies and examine the forces, ${\vec F}_1$ and ${\vec F}_2$.

When the restricted body reaches the origin, which is the only point of symmetry in this problem, the forces exerted by each of the massives are equal and opposite.  If we place the restricted body there at rest it will stay there at rest (i.e., in equilibrium) since no net force exists to move it anywhere else.  Should we displace the restricted body in a plane perpendicular to the line joining the massives (which we will call hereafter the center line) the component of the forces along the center line will cancel but the perpendicular components will add.  The action of the net force is to produce a restoring force back to the origin.  A particle so displaced will oscillate in that plane, starting from rest at the furthest distance from the origin it can reach (i.e., the turning points) and accelerating to the origin.  It’s kinetic energy will cause it to shoot past the origin and reach a turning point exactly the same distance but oppositely directed across the center line.  So far so good.  It is only when the particle is displaced along the center line that we see trouble in paradise.  In this scenario, the closer of the two massives now exerts a force with a component along the center line that is larger than the other.  The net force continues to accelerate the particle away from the origin and towards the closer of the two massives.  In this way, we can come to understand that the origin, which is the point halfway between the two massives (remember $\mu=0.5$ in this case) is an unstable equilibrium point and that the plane perpendicular to the center line and which contains the origin is a stable subspace in which the particle can oscillate in bounded motion.  But, as soon as the particle is displaced even minutely from this plane it will rush towards the closest of the two massives.

To transition from this idealized fixed mass scenario to the true circular restricted three body problem comes by regarding the frame as the rotating frame and taking into account the centrifugal force.  Since that force is always directed radially outward from the origin, 4 new equilibrium points arise as places where the inward pull of gravity due to both bodies balance the outward pointing centrifugal force.  Note that at the origin the centrifugal force is zero so the above argument is unaffected.  The following contour plot shows the shape of different values of $U(x,y,0)$ and the yellow dots denote schematically the location of the 5 equilibria.

As an example, the equilibrium point on the right that falls on the center line exists because both massives are tugging the particle to the left while the centrifugal force pulls it to the right.  The same hold true for the other equilibria with the obvious changes for direction.  The saddle curvature of $U(x,y,0)$ at the origin and the open ‘figure eight’ shape of the contour line surrounding the massives reflect the mixed equilibrium at the origin.  The contour plot indicates that the other two equilibrium points on the center are also mixed (also called unstable) equilibria.  The interpretation for the two points is not as obvious and we will deal with them (and the colinear points as the other three are called) in later posts. 

For now, we have a heuristic understanding of how equilibria arise at special locations in space based on the geometry of the pseudopotential.  Next month, we will quantitatively explore the interaction between the $U$ and the velocity of the particle in the rotating frame as governed by $2U = C_J + v^2$ where we will see how the geometry of $U$ and the value of Jacobi’s constant allows us to make global predictions as to the allowed regions in which the restricted particle can move.

Last month we formulated a model of and derived the equations of motion for the circular restricted three-body problem (CR3BP). In the CR3BP, two massive bodies, the primary and secondary (collectively called the massives) move in circular motion about their common center of mass, ignorant of the presence of a third body. Since this third body doesn’t influence or back-react onto the massives, it is called the restricted body and, for all practical purposes, it is massless.

The lack of back reaction means that the forces exerted by the massives on the restricted body result from a `schedule’ and do not arise naturally from the motion of all three. The restricted body is, therefore, driven by a external, time varying force making the inertial equations of motion non-autonomous.

Nonetheless, because of the requirement that the massives are moving in circular orbits, their motion has a fixed angular frequency (a feature unique in Keplerian orbits) and we can move to a co-rotating frame in which the equations of motion are autonomous. The cost of that simplification is that there are now ‘Coriolis’ type forces ($-2{\vec \omega}\times{\vec v}$ in basic classical mechanics) in the problem as identified by the presences of first-order time derivatives. The resulting equations of motion are:

\[ {\ddot x} – 2 {\dot y} – x = – (1-\mu) \frac{x+\mu}{d_1^3} – \mu \frac{x-1+\mu}{d_2^3} \; , \]

\[ {\ddot y} + 2 {\dot x} – y = – (1-\mu) \frac{y}{d_1^3} -\mu \frac{y}{d_2^3} \; , \]

and

\[ {\ddot z} = -(1-\mu) \frac{z}{d_1^3} – \mu \frac{z}{d_2^3} \; ,\]

with $d_1 = \sqrt{ (x+\mu)^2 + y^2 + z^2}$ and $d_2 = \sqrt{ (x – 1 + \mu)^2 + y^2 + z^2 }$.

Note that we’ve rearrange the result from last month so that the right-hand sides contain only the purely gravitational terms and that the terms linear in $x$ and $y$, which are centrifugal ($-{\vec \omega}\times({\vec \omega} \times {\vec r})$ in basic classical mechanics) and also the result of time derivatives in the co-rotating frame, are back on the left.

As discussed last month, despite the mixing of the in-plane inertial coordinates that results when moving to the co-rotating frame, the gravitational terms appear as if we had simply constructed them in the co-rotating frame to begin with. This is not a coincedence but a result of the fact that (unlike velocities) accelerations and forces are vectors that `live’ in the configuration space of the system.

From elementary classical mechanics, we know that gravitational forces are conservative and result from the taking a gradient of an appropriate potential. Therefore, we can expect that a similar potential, of the form

\[ {\tilde U} = -\frac{m_1}{d_1} – \frac{m_2}{d_2} = – \frac{1-\mu}{d_1} – \frac{\mu}{d_2} \; \]

should give rise to the right-hands sides by taking the appropriate derivatives.

We can verfiy this relatively easily by differentiating $U$ (easier that integrating the right-hand sides) and by organizing the derivatives. Defining

\[ d = \sqrt{ (x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2 } \; \]

we can easily verify that

\[ \frac{\partial}{\partial x} d= \frac{x-x_0}{d} \; \]

and

\[ \frac{\partial}{\partial x} \frac{1}{d} = -\frac{1}{d^2} \frac{\partial}{\partial x} d= – \frac{x-x_0}{d^3} \; , \]

with analogous expressions for the derivatives with respect to $y$ and $z$.

Following this program, we see that

\[ \frac{\partial}{\partial x} {\tilde U} = -\frac{(1-\mu)(x+\mu)}{d_1^3} – \frac{\mu(x-1+\mu)}{d_2^3} \; , \]

\[ \frac{\partial}{\partial y} {\tilde U} = -\frac{(1-\mu) y}{d_1^3} – \frac{\mu y}{d_2^3} \; , \]

and

\[ \frac{\partial}{\partial z} {\tilde U} = -\frac{(1-\mu) z}{d_1^3} – \frac{\mu z}{d_2^3} \; . \]

The next step is to recognize is to bring the centrifugal terms back to the right-hand side and to define a related potential

\[ U = \frac{1}{2} ( x^2 + y^2 ) + {\tilde U} = \frac{1}{2} ( x^2 + y^2 ) – \frac{1-\mu}{d_1} – \frac{\mu}{d_2} \;. \]

$U$ is known as the pseudo-potential of Jacobi, who first published its existence in 1836.

The equations of motion now take on the compact form (with $\partial_x U \equiv (\partial/\partial x) U$, etc.) of:

\[ {\ddot x} – 2 {\dot y} = \partial_x U \; , \]

\[ {\ddot y} + 2 {\dot x} = \partial_y U \; , \]

and

\[ {\ddot z} = \partial_z U \; . \]

The existence of the pseudo-potential enables us to make some remarkable predictions on the types of motion that can result from the CR3BP.  This is particularly useful given that the these equations have no general, closed form solution.

Start by taking each of the equations above and multiplying by the corresponding first derivative to get

\[ {\dot x}{\ddot x} – 2 {\dot x} {\dot y} = \partial_x U \cdot {\dot x} \; , \]

\[ {\dot y}{\ddot y} + 2 {\dot x} {\dot y} = \partial_y U \cdot {\dot y} \; , \]

and

\[ {\dot z}{\ddot z}  = \partial_z U \cdot {\dot z} \; . \]

We can see that by adding these three equations together we can eliminate the Coriolis terms and we are left with

\[ {\dot x}{\ddot x} + {\dot y}{\ddot y} + {\dot z}{\ddot z} = \partial_x U \cdot {\dot x} + \partial_y U \cdot {\dot y} + \partial_z U \cdot {\dot z} \; . \]

Both sides are total time derivatives that can be rewritten as

\[ \frac{1}{2} \frac{d}{dt} \left( {\dot x}^2 + {\dot y}^2 + {\dot z}^2 \right) = \frac{d}{dt} U \; .\]

Collecting the terms on the left-hand side and carrying out the integral gives

\[ \frac{1}{2} \left( {\dot x}^2 + {\dot y}^2 + {\dot z}^2 \right)-U = constant \; . \]

The form of the constant differs somewhat from book to book.  We will follow A. E. Roy here (except for the subscript) by expressing that constant as $-(1/2) C_J$ so that

\[ C_J  = 2 U – v^2 \; , \]

where $v^2 = {\dot x}^2 +{\dot y}^2 +{\dot z}^2$.  $C_J$ is universally referred to as Jacobi’s constant (even if the form of it isn’t always the same).

Next month, we’ll see how the presence of this constant enables us to understand the limits on the permitted motion without having to solve the equations of motion directly.  This global assessment is made possible by examining the geometry of the defining equation for $C_J$ above.

One of the most fascinating problems of classical physics is the circular, restricted, three-body problem (CR3BP). Not only is this problem full of rich behaviors that are not exhibited in the Kepler problem or in the problem of a single body moving in an external unchanging force, but the applications of the CR3BP are numerous. For example, the trajectory that the James Webb Space Telescope (JWST) moves along around the Sun-Earth/Moon L2 libration point can only be understood within the context of the CR3BP. Many other spacecraft have been found at these Libration Points, including Wind, ACE, SOHO, WMAP, and the soon-to-be-launched Nancy Grace Roman Space Telescope. (Side note: I’ve contributed to all of these missions, often as the flight dynamics lead, and JWST is currently flying the trajectory I designed for it).

In terms of the complexity and hierarchy of dynamics, the CR3BP positions itself at an interesting junction. It is far more complex than the Kepler problem but far simpler than the generic N-body problem making it at least tractable in many areas.

This post will be the first in a series that explores the CR3BP. Here we will establish the equations of motions, explaining along the way, why the word restricted gets attached. We will then show how moving to a rotating frame transforms the problem from being time-dependent (non-autonomous) to time-independent (autonomous). This rotating frame also has the benefit of revealing the existence of a pseudopotential that turns the equations of motion in the rotating frame into a Hamiltonian system (but we will have to wait until a future post to see that). At the end, I’ll comment on what makes the circular case much simpler than the elliptic, restricted, three-body problem.

The argument presented below is highly influenced by Chapter 5 of A E Roy’s book Orbital Motion, 3rd edition.

Our study begins with assuming that there are two massive bodies of mass $m_1$ and $m_2$ that are interacting gravitationally with each other. These bodies, referred to as the primary and the secondary, follow motion that is exactly described by the Kepler problem. We layer the additional requirement that the primary and secondary are in circular motion about their center-of-mass with a fixed distance between them of $L$. Under these assumptions, the bodies move at a constant angular rate

\[ n = \sqrt{\frac{G (m_1 + m_2)} {L^3}} \; . \]

Note the use of $n$ (for mean motion) in place of the more traditional elementary symbol $\omega$ for uniform circular motion is to honor the celestial mechanics use. The geometry can be further specified by fixing an inertial coordinate system $\{ {\hat X}, {\hat Y}, {\hat Z} \}$ at the center of mass (with ${\hat X}$ and ${\hat Y}$ in the orbital plane) and by denoting the position vectors of the two bodies with respect to the center of mass as ${\vec R}_1$ and ${\vec R}_2$.

In order to make the treatment universally applicable (e.g., to the Earth-Moon system as well as Sun-Jupiter), units are chosen to make the equations of motion non-dimensional. It is customary in all treatments to pick $L$ as the unit of distance between the masses and the total mass $M = m_1 + m_2$ as the unit of mass. The secondary’s mass is then denoted by the mass fraction parameter $\mu$ (i.e., $m_1 = 1 – \mu$ and $m_2 = \mu$). By its definition $\mu \leq 1/2$. There is no consensus on the choice of the unit of time and there are two reasonable choices: 1) set $G = 1$ so that the angular frequency $n=1$ (i.e., the period of motion is $2 \pi$) or 2) set the unit of time so that the period of motion is unity. We pick option 1) in order to more faithfully follow Roy.

At this point we can introduce the restricted body at a position ${\vec R} = X {\hat X} + Y {\hat Y} + Z {\hat Z}$ and form the equations of motion in the inertial coordinate system (note the mass of the restricted body $m$ plays no role if the forces are solely gravitational)

\[ {\ddot {\vec R}} = – \frac{(1-\mu) ({\vec R}-{\vec R}_1)}{| {\vec R}-{\vec R}_1|^3} – \frac{\mu ({\vec R}-{\vec R}_2)}{| {\vec R}-{\vec R}_2 |^3} \; .\]

While the above equation looks innocent enough, it is worth noting that the force is time-dependent since the primary and secondary execute their motion ignorant of the existence of the restricted body. This point can be driven home if we assume that the line joining the primary to the secondary is aligned with ${\hat X}$ at $t=0$ with the vector pointing to the secondary having a positive value. Their position vectors can be written explicitly as

\[ {\vec R}_1 = -(1-\mu) \left[ \cos t {\hat X} + \sin t {\hat Y} \right] \; \]

and

\[ {\vec R}_2 = \mu \left[ \cos t {\hat X} + \sin t {\hat Y} \right] \; . \]

As a check, one can confirm that at all times $m_1 {\vec R}_1 + m_2 {\vec R}_2 = 0$. Substituting in these expressions transforms the equations of motion (now written in component form] to:

\[ {\ddot X} = -(1-\mu)\frac{X+(1-\mu) \cos t}{D_1^3} + \mu \frac{X-\mu \cos t}{D_2^3} \; , \]

\[ {\ddot Y} = -(1-\mu)\frac{Y+(1-\mu) \sin t}{D_1^3} + \mu \frac{Y-\mu \sin t}{D_2^3} \; , \]

and

\[ {\ddot Z} = -(1-\mu)\frac{Z}{D_1^3} + \mu \frac{Z}{D_2^3} \; . \]

The distances $D_1$ \& $D_2$ between the primary \& secondary and the restricted body, respectively, are given by:

\[ D_1 = \sqrt{ \left[ X + \mu \cos t \right]^2 + \left[ Y + \mu \sin t\right]^2 + Z^2 } \; \]

and

\[ D_2 = \sqrt{ \left[ X – (1-\mu) \cos t \right]^2 + \left[ Y – (1- \mu) \sin t\right]^2 + Z^2 } \; .\]

This form clearly shows that the underlying equations are explicitly time-dependent. Of course, this makes sense since the motion of the restricted particle does not influence the motion of the primary and secondary. Their orbital motion acts like an external force that drives the system but is not part of the system. This is sometimes referred to as the restricted particle reacting to but not ‘back-reacting’ on the other two.

Despite this non-autonomous character, the CR3BP can be cast into a autonomous form. The key insight is to note that since the massive bodies move in circular motion about their center-of-mass, their angular rate is constant.  If we construct a co-rotating frame so that it rotates around ${\hat Z}$ and has it’s x-axis aligned with ${\hat X}$ at $t=0$ then $x_1 = -\mu$ and $x_2 = (1-\mu)$, with $y_1=y_2=z_1=z_2=0$ (with the lower case letters indicating the corresponding components of position measured by an observer moving with this frame). and, so, in a co-rotating frame they will appear fixed. We can relate the position components in the inertial frame with those in this co-rotating frame using

\[ X_i = x_i \cos(t) – y_i \sin(t) \; , \]

\[ Y_i = x_i \sin(t) + y_i \cos(t) \; , \]

and

\[ Z_i = z_i \; . \]

The $i$ index specifies the body in question: restricted $i$ is blank, primary $i=1$, and secondary $i=2$.

We now need to re-express the inertial equations of motion in the co-rotating frame.  We’ll tackle the terms appearing on the left-hand sides first and then right-hand sides second.

The time derivatives of the position components on the left-hand sides are a bit involved since the transformation to the co-rotating frame is time dependent.  Fortunately, only the in-plane components are involved, with the left-hand sides becoming:

\[ LHS_x = {\ddot x} \cos t – 2 {\dot x} \sin t – x \cos t – {\ddot y} \sin t – 2 {\dot y} \cos t + y \sin t \; \]

and

\[ LHS_y = {\ddot y} \cos t – 2 {\dot y} \sin t – y \cos t + {\ddot x} \sin t – 2 {\dot x} \cos t – x \sin t \; .\]

Note that the terms involving a single time derivative (either ${\dot x}$ or ${\dot y}$) appear with opposite signs in the two left-hand sides. This leads to taking linear combinations of the two to isolate the motion for $x$ and $y$, more cleanly:

\[ LHS_x \cos t + LHS_y \sin t = {\ddot x}-2 {\dot y}-x \; \]

and

\[ -LHS_x \sin t + LHS_y \cos t = {\ddot y} + 2{\dot x} – y \; . \]

When we look at the right-hand sides, we see a similar pattern where the in-plane, co-rotating coordinates $x_i$ and $y_i$ mix together when converting each of the independent inertial coordinates.  We can grind through substituting in the transformations above and then combining the right-hand sides in the same way as was done for the left-hand sides above (click here for a Jupyter notebook showing all the steps).  However, there is an easier way that involves more insight but far less industry.  Recognize that since the right-hand side doesn’t depend on time, it can be evaluated completely in the co-rotating frame.  In that frame the position of the restricted body is given by ${\vec R} = x {\hat x} + y {\hat y} + z {\hat z}$ and the specific forces (force/restricted mass) arise as if we hadn’t mixed the equations of motion. 

The distances between the massives and the restricted are:

\[ d_1 = \sqrt{ (x+\mu)^2 + y^2 + z^2 } \; \]

and

\[d_2 = \sqrt{ (x-1+\mu)^2 + y^2 + z^2} \; ,\]

where the switch has been made from upper case ‘D’ to lower case ‘d’ simply to emphasize the difference in structural forms.

\[ {\ddot x} – 2 {\dot y} = x – (1-\mu) \frac{x+\mu}{d_1^3} – \mu \frac{x-1+\mu}{d_2^3} \; , \]

\[ {\ddot y} + 2 {\dot x} = y – (1-\mu) \frac{y}{d_1^3} -\mu \frac{y}{d_2^3} \; , \]

and

\[ {\ddot z} = -(1-\mu) \frac{z}{d_1^3} – \mu \frac{z}{d_2^3} \; ,\]

This final form has successfully eliminated the non-autonomous nature of the inertial equations. This simplification was possible precisely because the co-rotating frame moves with a constant angular rate (i.e., there are no $\dot n$ in the transformation of the left-hand sides). This special case also explains why the difficulty in allowing the primary and secondary to move on elliptical paths is so much more pronounced.

Over the past two installments, we’ve seen that the percolation model, which has at its core simple, static probabilistic rules, exhibits three distinct behaviors.  For the occupation probability $p$ well below the critical probability $p_c$, the lattice is mostly empty with a few isolated islands of filled sites but with extremely little likelihood to have a spanning cluster joining opposites of the lattice.  For the occupation probability well above the critical probability, the lattice is mostly filled and, as a sort of mirror image, the lattice is mostly filled with a few isolated islands of vacant sites and with nearly perfect certainty that a spanning cluster is present.  The third behavior exists in the ‘local’ neighborhood near the critical probability where the observation evidence from computer experiments and the conventional wisdom indicates that many different length sizes exist.

In order to go from that qualitative picture summarized above to a quantitative structure, we need to define the measurements that we can make to the computational lattice.  The most important one is the cluster size distribution $n_s$ defined as the number $n$ of clusters with a size $s$.

The best way to understand $n_s$ is with a few sample lattices.  For the first example, let’s look at a percolation lattice of size $16 \times 16$ with an occupation probability of $p_{occ} = 0.4$. 

**image**

With a bit of patience one might convince oneself that there are 30 distinct clusters.  With a bit more effort, it becomes clear that there 13 distinct clusters with only one site and so we would assign $n_s(1) = 13$.  Likewise, there are 4 distinct clusters of size 2 so that $n_s(2) = 4$, and so on.

The full listing of Cluster ID versus Cluster Size for this lattice is:

note that a cluster with an ID of zero means that the site was not occupied. 

The last post introduced the basic model of site percolation in which an $M \times N$ lattice would have each of its sites occupied with a probability $p$.  It was demonstrated that by sweeping that probability up through a critical value $p_c$ there was a qualitative change in system.  Below $p_c$, there was a small likelihood of finding a spanning cluster, defined as a connected set of sites in which an unbroken path could be made between either the left and right sides of the lattice or the top and bottom of the lattice or both.  Above $p_c$ there it is almost certain that a spanning cluster will exist.  For a cluster to be connected, nearest neighbor sites must be occupied in the up-down direction or the right-left direction.

The percolation model finds itself in the center of a lot of analysis because it is a simple example of a phase transition, in this case what is known as a geometric phase transition.  Phase transitions are familiar in physics and we see numerous examples.  Perhaps the most common place one is $H_2 O$ going from solid form (ice) to liquid (water) to gas (steam).  Percolation provides a simple playground in which the basic ideas and theoretical foundations can be explored without bringing to account the vast machinery of statistical mechanics.   Despite this ‘toy’ status, percolation models are physically interesting in their own way as well.

An enormous amount of work went into understanding how phase transitions occur during the 1960s and 1970s and one of the most celebrated ideas is the notion that the system exhibits a type of length invariance at the transition.  The correlation length, properly defined, becomes infinite, and all length scales matter.

In this entry, I’ll try to show what is meant by this notion before turning, in subsequent columns, to the quantitative theory of critical exponents and scaling laws.  To this end, the results presented here will be more in the line of a numerical experimentalist presenting the findings of a variety of simulations without the presentation of much in the way of an underlying theory to explain the results from a compact set of premises.

These experiments consisted of running 50 trials of the percolation model for each of 4 lattice sizes: $L = 16$, $32$, $64$, and $128$ cells per linear dimension and over five separate occupation probabilities $p = 0.4$, $0.5$, $0.5927 \equiv p_c$, $0.7$, and $0.8$. 

Using the Hoshen-Kopelman cluster identification algorithm, the number of clusters in the lattice were labeled with a unique, sequential cluster id number.  These clusters labels afford the ability to both create visual representations of the lattice and perform quantitative analysis. 

We’ll start with the visual displays first and only at the end we’ll touch on the quantitative analysis, setting the stage for more detailed discussions in subsequent blogs.

While there are many ways to slice through the various combinations of lattice size and probability, the easiest one to start with is to look how cluster morphology changes as the lattice size changes at a given occupation probability.

For $p = 0.4$, which is well below the critical value $p_c$

there is a speckling of clusters throughout the whole lattice no matter how big the lattice size.  And, as shown last time, the probability of getting a spanning cluster is effectively zero.

Stepping the occupation probability up to the value $p = 0.5$ doesn’t significantly alter anything

although there are noticeably more clusters albeit somewhat larger ones; the lattices look grainer.

Appearances change remarkably at the critical value $p = p_c = 0.5927$

with the clusters both being bigger (each of the examples shown above have a spanning cluster) but also with the clusters being ‘rougher’ and ‘less compact’.  Having a way to quantify those observations will be the major focus of the numerical analysis that follows in the months to come.

Raising the occupation probability to $p = 0.7$

begins eliminating the rougher qualities seen when $p = p_c$ as the spanning cluster effectively crowds out the small clusters as it consumes more and more of the lattice.

By the time $p = 0.8$

the spanning cluster is leaving very little space for isolated clusters and, effectively, the entire lattice is consumed by a single, connected cluster that reaches all four edges.

It’s clear that at as the occupation probability is driven from low values through $p = p_c$ into high values, that there are three distinct behaviors. When $p << p_c$, the lattice supports small, relatively compact clusters whose presence gives it a rough appearance. When $ p >> p_c$, the lattice is consumed by essentially one large cluster, giving it a smooth or even appearance. At or around $p = p_c$, the lattice seems to have a mix of both small compact clusters and large spanning ones without either behavior dominating. This observation is what is meant by length invariance.

In the months to come, we’ll explore how to quantitatively describe these qualitative observations.