Monthly Archive: July 2026

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.