And now it’s all this

I just said what I said and it was wrong
Or was taken wrong

Next post Previous post

Lagrange points redux

August 17, 2016 at 6:23 PM by Dr. Drang

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.

Background

Lagrange

Lagrange points are points in the orbital plane of a planet¹ 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 $m_s$ and $m_p$, respectively, and the distance between their centers to be $R$. We’ll then introduce a nondimensional quantity, $\mu$, to represent the planet’s fraction of the total mass, $M$. Thus,

$$M=m_s+m_p$$ $$m_p=\mu M$$ $$m_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 $\mu R$ from the sun and $(1-\mu )R$ from the planet. Both the sun and the planet revolve about the barycenter with an angular speed $\omega$, where

$${\omega }^2=\frac{GM}{R^3}$$

The period is related to the angular speed through the relation

$$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:

$$T^2=\frac{4{\pi }^2{R}^3}{GM}$$

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.

Energy

Reference frame

Here is our system of sun (yellow), planet (blue), and satellite (black) laid out on an $x\text{-}y$ coordinate system. We put the origin at the barycenter and the $x\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.

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$$

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=-\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.²

In the expression for $U$,

• $m$ is the mass of the satellite,
• $r$ is its distance from the center of rotation (the barycenter),
• $r_s$ is its distance from the sun, and
• $r_p$ is its distance from the planet.

See the figure above for details.

Non-dimensional form

The first thing to do is substitute our previous expressions for $m_s$, $m_p$, and ${\omega }^2$ into the expression for $U$.

$$U=-\frac{GMm(1-\mu )}{r_s}-\frac{GMm\mu }{r_p}-\frac{GMm}{2R^3}{r}^2$$

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

$$r=\rho R,r_s={\rho }_{s}R,r_p={\rho }_{p}R$$

which allows us to write $U$ this way:

$$U=\frac{GMm}{R}[-\frac{1-\mu }{{\rho }_s}-\frac{\mu }{{\rho }_p}-\frac{1}{2}{\rho }^2]$$

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

$$u=-\frac{1-\mu }{{\rho }_s}-\frac{\mu }{{\rho }_p}-\frac{1}{2}{\rho }^2$$

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

Using $\rho$, ${\rho }_s$, and ${\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 points³ of $u$. We need to express $u$ in terms of $\xi$ and $\eta$, which we get from the Pythagorean formulas

$${\rho }^2={\xi }^2+{\eta }^2$$ $${\rho }_s^2=(\xi +\mu {)}^2+{\eta }^2$$ $${\rho }_p^2=[\xi -(1-\mu ){]}^2+{\eta }^2$$

So we end up with this,

$$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 $u$ as a function of $\xi$ (abscissa) and $\eta$ (ordinate). I’m plotting it for $\mu =0.1$, because that’s a value that allows us to see all the Lagrange points. (For the Earth-Sun system, $\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.)

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 $u$. You can click on the plot to see a bigger version.

The contour lines represent equal spacing in the value of $u$. 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.

This is not a coincidence that just happens to work out when $\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

$${\rho }_s^2=(\xi +\mu {)}^2+{\eta }^2$$ $${\rho }_p^2=[\xi -(1-\mu ){]}^2+{\eta }^2$$

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

$$(1-\mu ){\rho }_s^2+\mu {\rho }_p^2={\xi }^2+{\eta }^2+\mu (1-\mu )$$

Therefore

$${\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 $u$ to get

$$u=-\frac{1-\mu }{{\rho }_s}-\frac{\mu }{{\rho }_p}-\frac{1}{2}[(1-\mu ){\rho }_s^2+\mu {\rho }_p^2-\mu (1-\mu )]$$

or, after rearranging

$$u=-(1-\mu )(\frac{1}{{\rho }_s}+\frac{{\rho }_s^2}{2})-\mu (\frac{1}{{\rho }_p}+\frac{{\rho }_p^2}{2})+\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 $u$ with respect to $\xi$ and $\eta$ in order to find the stationary points. We’ll use the chain rule to do it:

$$\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$$ $$\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 ${\rho }_s$ and ${\rho }_p$ are simple:

$$\frac{\partial u}{\partial {\rho }_s}=(1-\mu )(\frac{1}{{\rho }_s^2}-{\rho }_s)$$ $$\frac{\partial u}{\partial {\rho }_p}=\mu (\frac{1}{{\rho }_p^2}-{\rho }_p)$$

The easy way to get the partials of ${\rho }_s$ and ${\rho }_p$ with respect to $\xi$ and $\eta$ is to take the total differentials of the expressions for ${\rho }_s^2$ and ${\rho }_p^2$:

$$2{\rho }_{s}d{\rho }_s=2(\xi +\mu )d\xi +2\eta d\eta$$ $$2{\rho }_{p}d{\rho }_p=2[\xi -(1-\mu )]d\xi +2\eta d\eta$$

Dividing the top equation by $2{\rho }_s$ and the bottom by $2{\rho }_p$ gives us

$$d{\rho }_s=\frac{\xi +\mu }{{\rho }_s}d\xi +\frac{\eta }{{\rho }_s}d\eta$$ $$d{\rho }_p=\frac{\xi -(1-\mu )}{{\rho }_p}d\xi +\frac{\eta }{{\rho }_p}d\eta$$

which means

$$\frac{\partial {\rho }_s}{\partial \xi }=\frac{\xi +\mu }{{\rho }_s},\frac{\partial {\rho }_s}{\partial \eta }=\frac{\eta }{{\rho }_s}$$ $$\frac{\partial {\rho }_p}{\partial \xi }=\frac{\xi -(1-\mu )}{{\rho }_p},\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:

$$\frac{\partial u}{\partial \xi }=(1-\mu )(\frac{1}{{\rho }_s^2}-{\rho }_s)\frac{\xi +\mu }{{\rho }_s}+\mu (\frac{1}{{\rho }_p^2}-{\rho }_p)\frac{\xi -(1-\mu )}{{\rho }_p}=0$$ $$\frac{\partial u}{\partial \eta }=(1-\mu )(\frac{1}{{\rho }_s^2}-{\rho }_s)\frac{\eta }{{\rho }_s}+\mu (\frac{1}{{\rho }_p^2}-{\rho }_p)\frac{\eta }{{\rho }_p}=0$$

Simplifying a bit we get

$$\frac{\partial u}{\partial \xi }=(1-\mu )(\frac{1}{{\rho }_s^3}-1)(\xi +\mu )+\mu (\frac{1}{{\rho }_p^3}-1)[\xi -(1-\mu )]=0$$ $$\frac{\partial u}{\partial \eta }=(1-\mu )(\frac{1}{{\rho }_s^3}-1)\eta +\mu (\frac{1}{{\rho }_p^3}-1)\eta =0$$

Solving the equations

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

$$\eta [(1-\mu )(\frac{1}{{\rho }_s^3}-1)+\mu (\frac{1}{{\rho }_p^3}-1)]=0$$

This means that either

$$\eta =0$$

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

$$(1-\mu )(\frac{1}{{\rho }_s^3}-1)+\mu (\frac{1}{{\rho }_p^3}-1)=0$$

Let’s explore this condition. We’ll move the terms that don’t involve ${\rho }_s$ or ${\rho }_p$ to the other side of the equation.

$$\frac{1-\mu }{{\rho }_s^3}+\frac{\mu }{{\rho }_p^3}=(1-\mu )+\mu =1$$

An obvious solution to this equation is ${\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 $\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 ${\rho }_s^3$. We can multiply through by ${\rho }_s^3{\rho }_p^3$ to get rid of the fractions:

$$(1-\mu ){\rho }_p^3+\mu {\rho }_s^2={\rho }_s^3{\rho }_p^3$$

And then solve for ${\rho }_s^3$:

$${\rho }_s^3=\frac{(1-\mu ){\rho }_p^3}{{\rho }_p^3-\mu }$$

We plug this into the first stationary equation to get

$$(1-\mu )(\frac{{\rho }_p^3-\mu }{(1-\mu ){\rho }_p^3}-1)(\xi +\mu )+\mu (\frac{1}{{\rho }_p^3}-1)[\xi -(1-\mu )]=0$$

which simplifies first to

$$(1-\mu )\left[\frac{\mu }{1-\mu }(1-\frac{1}{{\rho }_p^3})\right](\xi +\mu )+\mu (\frac{1}{{\rho }_p^3}-1)[\xi -(1-\mu )]=0$$

and then to

$$(1-\frac{1}{{\rho }_p^3})(\xi +\mu )-(1-\frac{1}{{\rho }_p^3})[\xi -(1-\mu )]=0$$

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

$$(1-\frac{1}{{\rho }_p^3})\{(\xi +\mu )-[\xi -(1-\mu )]\}=0$$

With this, we can say either

$$1-\frac{1}{{\rho }_p^3}=0$$

or

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

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

$$1-\frac{1}{{\rho }_p^3}=0$$

and therefore ${\rho }_p=1$, which means ${\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

$$\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$$ $$\frac{\partial u}{\partial \eta }=\eta [(1-\mu )\left(\frac{1-{\rho }_s^3}{{\rho }_s^3}\right)+\mu \left(\frac{1-{\rho }_p^3}{{\rho }_p^3}\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 $\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:

L1

Let’s start with L1, where

$${\rho }_s=\xi +\mu =1-{\rho }_p,{\rho }_p=(1-\mu )-\xi$$

For very small values of $\mu$, ${\rho }_p$ will also be small, so it’s convenient to put the whole equation in terms of ${\rho }_p$:

$$(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-\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$$

or

$$3(1-\mu ){\rho }_p^3(1-{\rho }_p+\frac{{\rho }_p^2}{3})-\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 ${\rho }_p$ is relatively small when $\mu$ is very small to get a closed form approximate solution:

$${\rho }_p^3\approx \frac{\mu }{3(1-\mu )}$$

or

$${\rho }_p\approx \sqrt[3]{\frac{\mu }{3(1-\mu )}}$$

Notice that $\mu$ and ${\rho }_p$ are at different levels of “small.” The cube/cube root relationship means that $\mu$ is much smaller than ${\rho }_p$.

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

L2

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

$${\rho }_s=\xi +\mu =1+{\rho }_p,{\rho }_p=\xi -(1-\mu )$$

so

$$(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-\mu ){\rho }_p^3(1+{\rho }_p+\frac{{\rho }_p^2}{3})+\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 ${\rho }_p$. It’s the same as the approximation for L1:

$${\rho }_p^3\approx \frac{\mu }{3(1-\mu )}$$

or

$${\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 $\mu =0.1$, a numerical solution of the exact expression gives ${\rho }_s=0.360$ which corresponds to $\xi =1.260$ as given in the table above. The approximate solution is ${\rho }_s=0.333$, which again is pretty far off.

L3

Finally, we have L3, where we have to be careful with the signs. Because they’re distances, ${\rho }_s$ and ${\rho }_p$ are positive, but the coordinate $\xi$ is negative.

$${\rho }_s=-(\xi +\mu ),{\rho }_p=-[\xi -(1-\mu )]=1+{\rho }_s$$

In this case, we’ll write the first stationary equation in terms of ${\rho }_s$.

$$-(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, ${\rho }_s$ is going to be close to 1, so we can introduce a small value, $\delta$, such that ${\rho }_s=1-\delta$. That turns the stationary equation into

$$-(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-\mu )\delta (2-\delta {)}^2(1-\delta +\frac{{\delta }^2}{3})+\mu (7-12\delta +6{\delta }^2-{\delta }^3)(1-\delta {)}^2=0$$

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

$$\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

$$\xi =-1-\mu +\delta \approx -(1+\frac{5}{12}\frac{\mu }{1-\mu })$$

For $\mu =0.1$, a numerical solution of the exact expression gives $\delta =0.0584$ which corresponds to $\xi =-1.042$ as given in the table above. The approximate solution is $\delta =0.0648$. The percent error in this approximation for $\delta$ is comparable to that of the earlier approximations for ${\rho }_p$.

The Sun-Earth Lagrange points

As mentioned earlier, $\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.

• For L1, a numerical solution of the exact equation gives ${\rho }_p=0.00997$, while the approximate solution is ${\rho }_p=0.01000$.
• For L2, a numerical solution of the exact equation gives ${\rho }_p=0.01003$, while the approximate solution is ${\rho }_p=0.01000$.
• For L3, both the numerical solution of the exact equation and the approximate solution are $\delta =1.75\times {10}^{-6}$. You have to go four more decimal places to see a difference.

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 ${\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.

Stability

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.

Next post Previous post

Site search

Meta

• drdrang at leancrew
• Blog archive
• RSS feed
• JSON feed
• Mastodon
• GitHub repositories

Recent posts

Credits

This work is licensed under a Creative Commons Attribution-Share Alike 3.0 Unported License.

© 2005–2023, Dr. Drang