Lagrange points redux

About a year ago, inspired by this GIF from the DSCOVR satellite, I wrote a post about the first Lagrange point of the Sun-Earth system, which is where DSCOVR is located. I wanted someday to return to the topic and get the locations of all the Lagrange points. This is the day.

Moon transit of Earth



Lagrange points are points in the orbital plane of a planet1 that orbit the sun with the same period as the planet. You might think you could put a satellite at any point along a planet’s orbital path and Kepler’s laws would ensure that it has the same period as the planet. But Kepler’s laws apply only to a two-body system. This is a three-body problem, in which the satellite’s motion is influenced by the gravitational pulls of both the sun and the planet. While there is no solution to the general three-body problem, the Lagrange points—so named because they were worked out by the 18th century natural philosopher Joseph-Louis Lagrange—represent special cases where the solution is possible.

In last year’s post, I showed how to find the first Lagrange point, L1, by balancing the two gravitational forces acting on it to create a centripetal acceleration that keeps a satellite at L1 in place. This approach works, but it’s a very non-Lagrangian way of solving the problem.

Lagrange was all about energy. He took Newtonian mechanics and recast it to eliminate the need to balance forces and inertias. In Lagrangian mechanics, you get solutions by taking derivatives of the kinetic and potential energy functions. It’s an elegant technique, well suited to the explosion of analysis on the Continent back at that time.

Two-body results

Let’s start by assuming we’ve already solved the two-body problem of a sun and its planet in a circular orbit. We’ll take their masses to be msm_s and mpm_p, respectively, and the distance between their centers to be RR. We’ll then introduce a nondimensional quantity, μ\mu, to represent the planet’s fraction of the total mass, MM. Thus,

M=ms+mpM = m_s + m_p mp=μMm_p = \mu M ms=(1μ)Mm_s = (1 - \mu)M

The center of mass of the two-body system—which astronomers call the barycenter because it sounds more scientific—is on a line between the two bodies a distance μR\mu R from the sun and (1μ)R(1 - \mu)R from the planet. Both the sun and the planet revolve about the barycenter with an angular speed ω\omega, where

ω2=GMR3\omega^2 = \frac{GM}{R^3}

The period is related to the angular speed through the relation

T=2πωT = \frac{2\pi}{\omega}

which leads to the well-known expression for Kepler’s Third Law, which states that the square of the period is proportional to the cube of the distance:

T2=4π2R3GMT^2 = \frac{4 \pi^2 R^3}{G M}

With these preliminaries out of the way, let’s move on to finding the Lagrange points. I want to start by pointing you to an excellent online resource, Richard Fitzpatrick’s Newtonian Dyanamics, which is available in both PDF and HTML format. Fitzpatrick, who teaches at the University of Texas at Austin (hook ’em), does a very nice job of explaining both the two-body problem and the restricted three-body problem. There’s one trick in particular that I stole directly from him to simplify a potential energy expression.


Reference frame

Here is our system of sun (yellow), planet (blue), and satellite (black) laid out on an x-yx\text{-}y coordinate system. We put the origin at the barycenter and the x-axisx\text{-axis} on the line between the sun and the planet. Furthermore, we’re going to have our coordinate system rotate at a constant angular speed of ω\omega, precisely matching the movement of the sun and the planet about the barycenter. This will be our reference frame for the analysis. The advantage of using a rotating reference frame is that the sun and planet are, by definition, motionless in this frame, and our search for Lagrange points is reduced to finding points where the satellite will be motionless, too.

Planet coordinate system

You may object to using a rotating reference frame.

A rotating reference frame isn’t inertial. That’s true.

You can’t do an analysis in a non-inertial reference frame. That’s not true.

Non-inertial reference frames are perfectly fine as long as you account for the acceleration terms correctly. This is the deeper truth behind d’Alembert’s Principle. Most of us learn d’Alembert’s Principle as simply moving the acceleration term in Newton’s Second Law over to the other side of the equation and treating it as an additional force.

𝐅=m𝐚𝐅m𝐚=0\mathbf{F} = m\: \mathbf{a} \quad \Longleftrightarrow \quad \mathbf{F} - m\: \mathbf{a} = 0

But d’Alembert works in an energy context, too.

Potential energy

In our rotating frame of reference, the potential energy of the satellite has three terms.

U=GmsmrsGmprp12m(rω)2U = -\frac{G m_s m}{r_s} - \frac{G m_p}{r_p} - \frac{1}{2} m (r \omega)^2

The first two terms are the gravitational potential energy due to the sun and the planet, respectively, and the third term is the centrifugal potential energy due to the rotating frame. The third term wouldn’t appear in a potential energy expression written for an intertial frame.2

In the expression for UU,

See the figure above for details.

Non-dimensional form

The first thing to do is substitute our previous expressions for msm_s, mpm_p, and ω2\omega^2 into the expression for UU.

U=GMm(1μ)rsGMmμrpGMm2R3r2U = -\frac{G M m (1 - \mu)}{r_s} - \frac{G M m \mu}{r_p} - \frac{G M m}{2 R^3} r^2

We’re starting to see some common terms we can factor out. We can do even better if we rewrite the rr terms using nondimensional variables,

r=ρR,rs=ρsR,rp=ρpRr = \rho R, \quad r_s = \rho_s R, \quad r_p = \rho_p R

which allows us to write UU this way:

U=GMmR[1μρsμρp12ρ2]U = \frac{GMm}{R} \left[ -\frac{1-\mu}{\rho_s} - \frac{\mu}{\rho_p} - \frac{1}{2}\rho^2 \right]

All of the terms with units have been factored out of the brackets into a constant scaling term. Finding the stationary points of UU now reduces to finding the stationary points of the nondimensional expression within the brackets, which we’ll call uu.

u=1μρsμρp12ρ2u = -\frac{1-\mu}{\rho_s} - \frac{\mu}{\rho_p} - \frac{1}{2}\rho^2

In effect, we’ve switched from the x-yx\text{-}y coordinate system of the figure above to the ξ-η\xi\text{-}\eta system shown below.

Nondimensional coordinate system

Using ρ\rho, ρs\rho_s, and ρp\rho_p makes for a compact expression, but it isn’t convenient for plotting, which is what I want to do to help find the stationary points3 of uu. We need to express uu in terms of ξ\xi and η\eta, which we get from the Pythagorean formulas

ρ2=ξ2+η2\rho^2 = \xi^2 + \eta^2 ρs2=(ξ+μ)2+η2\rho_s^2 = (\xi + \mu)^2 + \eta^2 ρp2=[ξ(1μ)]2+η2\rho_p^2 = [\xi - (1 - \mu)]^2 + \eta^2

So we end up with this,

u=1μ(ξ+μ)2+η2μ[ξ(1μ)]2+η212(ξ2+η2)u = -\frac{1-\mu}{\sqrt{(\xi + \mu)^2 + \eta^2}} - \frac{\mu}{\sqrt{[\xi - (1 - \mu)]^2 + \eta^2}} - \frac{1}{2}(\xi^2 + \eta^2)

which is a nasty mess, but we have computers to keep track of everything, so there’s no need to worry about losing terms.

Plotting the potential energy

Here’s the contour plot of uu as a function of ξ\xi (abscissa) and η\eta (ordinate). I’m plotting it for μ=0.1\mu = 0.1, because that’s a value that allows us to see all the Lagrange points. (For the Earth-Sun system, μ=0.000003\mu = 0.000003, which would put L1 and L2 so close to the Earth itself we wouldn’t be able to distinguish them at this scale.)

Contour lines for Lagrange points

The dirty yellow dot is the sun, the blue dot is the planet, the × is the barycenter, and the various crosses are the stationary points of uu. You can click on the plot to see a bigger version.

The contour lines represent equal spacing in the value of uu. They range from dark blue for the lowest points to dark red for the highest. We see that L1, L2, and L3 are colinear with the sun and planet and are at saddle points. L4 and L5 are at local maxima. The coordinates of the points, which I calculated using techniques we’ll get into later, are as follows:

Point ξ\xi    η\eta   
L1 0.609 0.000
L2 1.260 0.000
L3 -1.042 0.000
L4 0.400 0.866
L5 0.400 -0.866

The ξ\xi coordinates of L1, L2, and L3 pretty much have to be calculated numerically. There’s no nice closed-form solution to get those values. But there is a simple, non-computational way to get the positions of L4 and L5, and the clue is in the values you see in the table.

That 0.866 you see for the η\eta value is the sine of 60°, and the 0.400 is exactly 0.1 less than the cosine of 60°. Remember that the sun is 0.1 to the left of the origin and the planet is 0.9 to the right of the origin. Putting this all together, we see that L4 is at the intersection of a 60° line up and out from the sun and a 60° line up and back from the planet. Similarly for L5, except that the lines are 60° down instead of up. Which means that L4 and L5 form equilateral triangles with the sun and the planet.

Equilateral triangles

This is not a coincidence that just happens to work out when μ=0.1\mu = 0.1. It’s true regardless of the mass distribution between the sun and the planet. In the next section, we’ll prove that, but the math gets messy. If you want to just take it on faith, skip this next section.

Equilateral triangle diversion

For the fearless few, we’re going to use that trick I found in Richard Fitzpatrick’s book. There’s nothing especially hard in this; it’s just a lot of tedious algebra, and I’m going to show all the steps. Textbooks usually don’t for reasons of space, but there’s a lot of space on a web page.

Recall that

ρs2=(ξ+μ)2+η2\rho_s^2 = (\xi + \mu)^2 + \eta^2 ρp2=[ξ(1μ)]2+η2\rho_p^2 = [\xi - (1 - \mu)]^2 + \eta^2

If we multiply the first of these by 1μ1 - \mu and second by μ\mu and add them together, we get (after some cancellation)

(1μ)ρs2+μρp2=ξ2+η2+μ(1μ)(1 - \mu)\rho_s^2 + \mu \rho_p^2 = \xi^2 + \eta^2 + \mu(1 - \mu)


ξ2+η2=ρ2=(1μ)ρs2+μρp2μ(1μ)\xi^2 + \eta^2 = \rho^2 = (1 - \mu)\rho_s^2 + \mu \rho_p^2 - \mu(1 - \mu)

We can substitute this into the compact expression for uu to get

u=1μρsμρp12[(1μ)ρs2+μρp2μ(1μ)]u = -\frac{1-\mu}{\rho_s} - \frac{\mu}{\rho_p} - \frac{1}{2} \left[ (1 - \mu)\rho_s^2 + \mu \rho_p^2 - \mu(1 - \mu) \right]

or, after rearranging

u=(1μ)(1ρs+ρs22)μ(1ρp+ρp22)+μ(1μ)2u = -(1 - \mu) \left(\frac{1}{\rho_s} + \frac{\rho_s^2}{2} \right) - \mu \left( \frac{1}{\rho_p} + \frac{\rho_p^2}{2} \right) + \frac{\mu (1 - \mu)}{2}

Equations for the stationary points

What good is this? Well, although it may not seem like it, it actually makes it a little easier to take the partial derivatives of uu with respect to ξ\xi and η\eta in order to find the stationary points. We’ll use the chain rule to do it:

uξ=uρsρsξ+uρpρpξ=0\frac{\partial u}{\partial \xi} = \frac{\partial u}{\partial \rho_s}\frac{\partial \rho_s}{\partial \xi} + \frac{\partial u}{\partial \rho_p}\frac{\partial \rho_p}{\partial \xi} = 0 uη=uρsρsη+uρpρpη=0\frac{\partial u}{\partial \eta} = \frac{\partial u}{\partial \rho_s}\frac{\partial \rho_s}{\partial \eta} + \frac{\partial u}{\partial \rho_p}\frac{\partial \rho_p}{\partial \eta} = 0

The partial derviatives with respect to ρs\rho_s and ρp\rho_p are simple:

uρs=(1μ)(1ρs2ρs)\frac{\partial u}{\partial \rho_s} = (1 - \mu) \left( \frac{1}{\rho_s^2} - \rho_s \right) uρp=μ(1ρp2ρp)\frac{\partial u}{\partial \rho_p} = \mu \left( \frac{1}{\rho_p^2} - \rho_p \right)

The easy way to get the partials of ρs\rho_s and ρp\rho_p with respect to ξ\xi and η\eta is to take the total differentials of the expressions for ρs2\rho_s^2 and ρp2\rho_p^2:

2ρsdρs=2(ξ+μ)dξ+2ηdη2 \rho_s\; \mathrm{d}\rho_s = 2(\xi + \mu)\;\mathrm{d}\xi + 2\eta\; \mathrm{d}\eta 2ρpdρp=2[ξ(1μ)]dξ+2ηdη2 \rho_p\; \mathrm{d}\rho_p = 2[\xi - (1 - \mu)]\;\mathrm{d}\xi + 2\eta\; \mathrm{d}\eta

Dividing the top equation by 2ρs2 \rho_s and the bottom by 2ρp2 \rho_p gives us

dρs=ξ+μρsdξ+ηρsdη\mathrm{d}\rho_s = \frac{\xi + \mu}{\rho_s} \mathrm{d}\xi + \frac{\eta}{\rho_s} \mathrm{d}\eta dρp=ξ(1μ)ρpdξ+ηρpdη\mathrm{d}\rho_p = \frac{\xi - (1- \mu)}{\rho_p} \mathrm{d}\xi + \frac{\eta}{\rho_p} \mathrm{d}\eta

which means

ρsξ=ξ+μρs,ρsη=ηρs\frac{\partial \rho_s}{\partial \xi} = \frac{\xi + \mu}{\rho_s}, \qquad \qquad \frac{\partial \rho_s}{\partial \eta} = \frac{\eta}{\rho_s} ρpξ=ξ(1μ)ρp,ρpη=ηρp\frac{\partial \rho_p}{\partial \xi} = \frac{\xi - (1 - \mu)}{\rho_p}, \qquad \quad \frac{\partial \rho_p}{\partial \eta} = \frac{\eta}{\rho_p}

Now we have all the pieces needed to build the equations for the stationary points:

uξ=(1μ)(1ρs2ρs)ξ+μρs+μ(1ρp2ρp)ξ(1μ)ρp=0\frac{\partial u}{\partial \xi} = (1 - \mu) \left( \frac{1}{\rho_s^2} - \rho_s \right) \frac{\xi + \mu}{\rho_s} + \mu \left( \frac{1}{\rho_p^2} - \rho_p \right) \frac{\xi - (1 - \mu)}{\rho_p} = 0 uη=(1μ)(1ρs2ρs)ηρs+μ(1ρp2ρp)ηρp=0\frac{\partial u}{\partial \eta} = (1 - \mu) \left( \frac{1}{\rho_s^2} - \rho_s \right) \frac{\eta}{\rho_s} + \mu \left( \frac{1}{\rho_p^2} - \rho_p \right) \frac{\eta}{\rho_p} = 0

Simplifying a bit we get

uξ=(1μ)(1ρs31)(ξ+μ)+μ(1ρp31)[ξ(1μ)]=0\frac{\partial u}{\partial \xi} = (1 - \mu) \left( \frac{1}{\rho_s^3} - 1 \right)(\xi + \mu) + \mu \left( \frac{1}{\rho_p^3} - 1 \right)[\xi - (1 - \mu)] = 0 uη=(1μ)(1ρs31)η+μ(1ρp31)η=0\frac{\partial u}{\partial \eta} = (1 - \mu) \left( \frac{1}{\rho_s^3} - 1 \right) \eta + \mu \left( \frac{1}{\rho_p^3} - 1 \right) \eta = 0

Solving the equations

The second equation is the key. First, we can factor out the η\eta:

η[(1μ)(1ρs31)+μ(1ρp31)]=0\eta \left[ (1 - \mu) \left( \frac{1}{\rho_s^3} - 1 \right) + \mu \left( \frac{1}{\rho_p^3} - 1 \right) \right] = 0

This means that either

η=0\eta = 0

which is what leads us to L1, L2, and L3 (we’ll get to that later), or

(1μ)(1ρs31)+μ(1ρp31)=0(1 - \mu) \left( \frac{1}{\rho_s^3} - 1 \right) + \mu \left( \frac{1}{\rho_p^3} - 1 \right) = 0

Let’s explore this condition. We’ll move the terms that don’t involve ρs\rho_s or ρp\rho_p to the other side of the equation.

1μρs3+μρp3=(1μ)+μ=1\frac{1 - \mu}{\rho_s^3} + \frac{\mu}{\rho_p^3} = (1 - \mu) + \mu = 1

An obvious solution to this equation is ρs=ρp=1\rho_s = \rho_p = 1, which will work for all values of μ\mu. What we don’t know, though, is whether that’s the only solution for η0\eta \ne 0. To see if it is, we have to combine this result with the first stationary equation.

Let’s start by solving for ρs3\rho_s^3. We can multiply through by ρs3ρp3\rho_s^3 \rho_p^3 to get rid of the fractions:

(1μ)ρp3+μρs2=ρs3ρp3(1 - \mu) \rho_p^3 + \mu \rho_s^2 = \rho_s^3 \rho_p^3

And then solve for ρs3\rho_s^3:

ρs3=(1μ)ρp3ρp3μ\rho_s^3 = \frac{(1 - \mu) \rho_p^3}{\rho_p^3 - \mu}

We plug this into the first stationary equation to get

(1μ)(ρp3μ(1μ)ρp31)(ξ+μ)+μ(1ρp31)[ξ(1μ)]=0(1 - \mu) \left( \frac{\rho_p^3 - \mu}{(1 - \mu) \rho_p^3} - 1 \right)(\xi + \mu) + \mu \left( \frac{1}{\rho_p^3} - 1 \right)[\xi - (1 - \mu)] = 0

which simplifies first to

(1μ)[μ1μ(11ρp3)](ξ+μ)+μ(1ρp31)[ξ(1μ)]=0(1 - \mu) \left[ \frac{\mu}{1 - \mu} \left(1 - \frac{1}{\rho_p^3} \right) \right](\xi + \mu) + \mu \left( \frac{1}{\rho_p^3} - 1 \right)[\xi - (1 - \mu)] = 0

and then to

(11ρp3)(ξ+μ)(11ρp3)[ξ(1μ)]=0\left(1 - \frac{1}{\rho_p^3} \right) (\xi + \mu) - \left(1 - \frac{1}{\rho_p^3} \right) [\xi - (1 - \mu)] = 0

Once again, we can factor out a common term and simplify:

(11ρp3){(ξ+μ)[ξ(1μ)]}=0\left(1 - \frac{1}{\rho_p^3} \right) \left\{ (\xi + \mu) - [\xi - (1 - \mu)] \right\} = 0

With this, we can say either

11ρp3=01 - \frac{1}{\rho_p^3} = 0


(ξ+μ)[ξ(1μ)]=0(\xi + \mu) - [\xi - (1 - \mu)] = 0

But the second of these is impossible because the ξ\xi and μ\mu terms cancel, leaving 1=01 = 0. So the only solution for η0\eta \ne 0 is

11ρp3=01 - \frac{1}{\rho_p^3} = 0

and therefore ρp=1\rho_p = 1, which means ρs=1\rho_s = 1, confirming our guess about the equilateral triangle solution for L4 and L5.

Determining the colinear positions

OK, now that we’ve confirmed the equilateral triangle postions for L4 and L5, let’s explore the colinear positions, L1, L2, and L3.

The two equations that must be satisfied for every Lagrange point are

uξ=(1μ)(1ρs3ρs3)(ξ+μ)+μ(1ρp3ρp3)[ξ(1μ)]=0\frac{\partial u}{\partial \xi} = (1 - \mu) \left( \frac{1 - \rho_s^3}{\rho_s^3} \right)(\xi + \mu) + \mu \left( \frac{1 - \rho_p^3}{\rho_p^3} \right)[\xi - (1 - \mu)] = 0 uη=η[(1μ)(1ρs3ρs3)+μ(1ρp3ρp3)]=0\frac{\partial u}{\partial \eta} = \eta \left[ (1 - \mu) \left( \frac{1 - \rho_s^3}{\rho_s^3} \right) + \mu \left( \frac{1 - \rho_p^3}{\rho_p^3} \right) \right] = 0

(If you’re wondering where these equations came from, it’s because you skipped over the previous section. The path to enlightenment is not easy, grasshopper.)

An obvious condition that solves the second equation is η=0\eta = 0. That’s the value of η\eta for L1, L2, and L3. All we need to do then is pull three solutions for ξ\xi out of the first equation. We’ll refer to this layout of the points to specialize the equation for each of the points:

Colinear points


Let’s start with L1, where

ρs=ξ+μ=1ρp,ρp=(1μ)ξ\rho_s = \xi + \mu = 1 - \rho_p, \qquad \rho_p = (1 - \mu) - \xi

For very small values of μ\mu, ρp\rho_p will also be small, so it’s convenient to put the whole equation in terms of ρp\rho_p:

(1μ)(1(1ρp)3(1ρp)2)μ(1ρp3ρp2)=0(1 - \mu) \left( \frac{1 - (1 - \rho_p)^3}{(1 - \rho_p)^2} \right) - \mu \left( \frac{1 - \rho_p^3}{\rho_p^2} \right) = 0

Expanding and collecting terms gives

(1μ)(3ρp(1ρp+ρp2/3)(1ρp)2)μ((1ρp)(1+ρp+ρp2)ρp2)=0(1 - \mu) \left( \frac{3\rho_p (1 - \rho_p + \rho_p^2/3)}{(1 - \rho_p)^2} \right) - \mu \left( \frac{(1 - \rho_p)(1 + \rho_p + \rho_p^2)}{\rho_p^2} \right) = 0


3(1μ)ρp3(1ρp+ρp23)μ(1ρp)3(1+ρp+ρp2)=03 (1 - \mu) \rho_p^3 \left( 1 - \rho_p + \frac{\rho_p^2}{3} \right) - \mu (1 - \rho_p)^3 (1 + \rho_p + \rho_p^2) = 0

Most numerical equation solving routines will have no trouble with this equation, but as I said earlier, there is no simple closed-form solution for it. We can, however, take advantage of the fact that ρp\rho_p is relatively small when μ\mu is very small to get a closed form approximate solution:

ρp3μ3(1μ)\rho_p^3 \approx \frac{\mu}{3 (1 - \mu)}


ρpμ3(1μ)3\rho_p \approx \sqrt[3]{\frac{\mu}{3 (1- \mu)}}

Notice that μ\mu and ρp\rho_p are at different levels of “small.” The cube/cube root relationship means that μ\mu is much smaller than ρp\rho_p.

For μ=0.1\mu = 0.1, a numerical solution of the exact expression gives ρs=0.291\rho_s = 0.291 which corresponds to ξ=0.609\xi = 0.609 as given in the table above. The approximate solution is ρs=0.333\rho_s = 0.333, which is pretty far off, mainly because ρs\rho_s just isn’t small enough.


The determination of L2 follows the same pattern. For this position, with the point beyond the planet,

ρs=ξ+μ=1+ρp,ρp=ξ(1μ)\rho_s = \xi + \mu = 1 + \rho_p, \qquad \rho_p = \xi - (1 - \mu)


(1μ)(1(1+ρp)3(1+ρp)2)+μ(1ρp3ρp2)=0(1 - \mu) \left( \frac{1 - (1 + \rho_p)^3}{(1 + \rho_p)^2} \right) + \mu \left( \frac{1 - \rho_p^3}{\rho_p^2} \right) = 0

After expanding, collecting, and rearranging as we did above, we get

3(1μ)ρp3(1+ρp+ρp23)+μ(1ρp3)(1+ρp)2=0-3 (1 - \mu) \rho_p^3 \left( 1 + \rho_p + \frac{\rho_p^2}{3} \right) + \mu (1 - \rho_p^3) (1 + \rho_p)^2 = 0

As with L1, this can be solved numerically without much trouble, but there is a decent closed-form approximation for small μ\mu and ρp\rho_p. It’s the same as the approximation for L1:

ρp3μ3(1μ)\rho_p^3 \approx \frac{\mu}{3 (1 - \mu)}


ρpμ3(1μ)3\rho_p \approx \sqrt[3]{\frac{\mu}{3 (1- \mu)}}

This puts the L2 position about as far outside the planet’s orbit as L1 is inside the planet’s orbit.

For μ=0.1\mu = 0.1, a numerical solution of the exact expression gives ρs=0.360\rho_s = 0.360 which corresponds to ξ=1.260\xi = 1.260 as given in the table above. The approximate solution is ρs=0.333\rho_s = 0.333, which again is pretty far off.


Finally, we have L3, where we have to be careful with the signs. Because they’re distances, ρs\rho_s and ρp\rho_p are positive, but the coordinate ξ\xi is negative.

ρs=(ξ+μ),ρp=[ξ(1μ)]=1+ρs\rho_s = -(\xi + \mu), \qquad \rho_p = -[\xi - (1 - \mu)] = 1 + \rho_s

In this case, we’ll write the first stationary equation in terms of ρs\rho_s.

(1μ)[1ρs3ρs2]μ[1(1+ρs)3(1+ρs)2]=0-(1 - \mu) \left[ \frac{1 - \rho_s^3}{\rho_s^2} \right] - \mu \left[ \frac{1 - (1 + \rho_s)^3}{(1 + \rho_s)^2} \right] = 0

In this case, ρs\rho_s is going to be close to 1, so we can introduce a small value, δ\delta, such that ρs=1δ\rho_s = 1 - \delta. That turns the stationary equation into

(1μ)[1(1δ)3(1δ)2]μ[1(1+(1δ))3(1+(1δ))2]=0-(1 - \mu) \left[ \frac{1 - (1 - \delta)^3}{(1 - \delta)^2} \right] - \mu \left[ \frac{1 - (1 + (1 - \delta))^3}{(1 + (1 - \delta))^2} \right] = 0

which looks like a real mess, but as before we expand, collect, and rearrange to get

3(1μ)δ(2δ)2(1δ+δ23)+μ(712δ+6δ2δ3)(1δ)2=0-3 (1 - \mu) \delta (2 - \delta)^2 \left( 1 - \delta + \frac{\delta^2}{3} \right) + \mu (7 - 12\delta + 6\delta^2 - \delta^3)(1 - \delta)^2 = 0

Ignoring the higher-order terms in δ\delta, we get the approximation

δ712μ1μ\delta \approx \frac{7}{12} \frac{\mu}{1 - \mu}

In this case, μ\mu and δ\delta are at about the same order of “small.”

Using this approximation, the ξ\xi coordinate is

ξ=1μ+δ(1+512μ1μ)\xi = -1 - \mu + \delta \approx - \left( 1 + \frac{5}{12} \frac{\mu}{1 - \mu} \right)

For μ=0.1\mu = 0.1, a numerical solution of the exact expression gives δ=0.0584\delta = 0.0584 which corresponds to ξ=1.042\xi = -1.042 as given in the table above. The approximate solution is δ=0.0648\delta = 0.0648. The percent error in this approximation for δ\delta is comparable to that of the earlier approximations for ρp\rho_p.

The Sun-Earth Lagrange points

As mentioned earlier, μ=0.000003\mu = 0.000003 for the Sun-Earth system. With such a small value of μ\mu, the approximations developed above should be pretty accurate. Let’s see.

As expected, the approximations are quite good. Probably not good enough for NASA, but good enough for a blog post.

The real value of the approximate formulas is not for computation, it’s for insight. By seeing how ρp\rho_p and δ\delta scale with μ\mu, we get a sense of how the positions of the colinear Lagrange points change with changing mass distributions.


It’s often said that L4 and L5 are the stable Lagrange points. This seems wrong, because those points are at local maxima of the potential energy, not local minima, and stability is associated with minima. My understanding is that the stability comes from Coriolis forces, which tend to keep objects in orbit around L4 and L5. We didn’t include a Coriolis term in our potential energy expression because our analysis was designed to find places where the satellites would be stationary in our rotating frame of reference. Coriolis forces arise only when a body is moving relative to the rotating frame.

I may look into redoing the analysis with a Coriolis term. Check back in another year.

Update 08/18/2016 8:23 AM
The Trojan asteroids are clustered around the L4 and L5 positions of the Sun-Jupiter system. They got a mention from Jason Snell and Stephen Hackett on this week’s episode of their Liftoff podcast, which I just listened to this morning. The plan of the proposed Lucy space mission is to visit five of the Trojan satellites.

A tip from Jeff Youngstrom on Twitter led me to this remarkable page by Petr Scheirich, which has a wealth of graphics related to comets and asteroids, including this animation of the Trojan (green) and Hilda (red) groups as viewed in a reference frame that rotates with Jupiter.

Petr Scheirich animation

The animation covers, I believe, one Jovian year. The in-and-out movement of Jupiter represents its elliptical orbit from perihelion to aphelion, and you can track the orbits of at least some of the green dots around the L4 and L5 positions.

  1. Although we tend to be most interested in the Sun-Earth Lagrange points, there are similar points for every sun-planet combination and for every planet-moon combination, too. 

  2. And it’s not a coincidence that it looks like a kinetic energy term with the sign changed. D’Alembert strikes again! 

  3. Stationary points are where the function is at a local maximum, minimum, or saddle point. They’re the points where the slopes of the function’s surface are zero.