Skip to content

Original ResearchVolume 3, Issue 1, Fall 2026

Five Places to Park a Spacecraft, Two of Which Are Stable for a Strange Reason

We located all five Lagrange points of the circular restricted three-body problem numerically, checked them against classical series and exact geometry, and found by eigenvalue and Monte Carlo integration exactly why two of them are stable: not because they sit in a well, but because the Coriolis force turns a hilltop into a whirlpool.

Written jointly by the Science Journaling Club

Study
Computational study
Review
Peer-edited by the club review board
Length
12 sectionsabout 31 min

Abstract

Two bodies orbit their common centre of mass in a circle. In the frame that rotates with them there are exactly five points where a third, massless body can sit still: three on the line through the two masses (L1, L2, L3), and two forming equilateral triangles with them (L4, L5). We located all five numerically for five mass ratios μ spanning five orders of magnitude, from Sun-Earth (μ = 3.00×10−6) to an illustrative equal-mass pair (μ = 0.5), by two independent numerical methods that agree with each other to machine precision (2.2×10−16) and with the classical small-μ series for the collinear points to four or more significant figures, and by exact algebra for the triangular ones (residual <4×10−16). Eigenvalues of the linearised equations of motion show L1, L2 and L3 unstable at every mass ratio tried, and show L4/L5 stable only below a critical mass ratio that our own numerical bisection places at μcrit = 0.03852085, against the classical Gascheau (1843) / Routh (1875) value of 0.03852090, a difference of 4.4×10−8 consistent with the square-root ill-conditioning of an eigenvalue problem at a genuine matrix degeneracy, not a bug. Direct integration of test particles (fixed fourth-order Runge-Kutta, Jacobi constant conserved to <4×10−13 relative up to first escape in every run) shows the same particle, same perturbation, same code, staying within 1000 times its starting displacement for the whole run at Earth-Moon's L4 and leaving it by t=15.0 at Pluto-Charon's L4, whose mass ratio (0.1087, 2.82× the critical value) exceeds threshold. A Monte Carlo experiment of 1,500 randomly directed perturbations per mass ratio, seeded at 20260914, put the escape probability at L4 at 0% ± 0% at 0.95×μcrit, 18.8% ± 1.0% at 0.98×μcrit, and 100% ± 0% by 1.20×μcrit, with the running estimate visibly converged by a few hundred trials. The stability of L4/L5 has nothing to do with sitting in a potential well: they sit on a hill of the effective potential at every mass ratio, stable or not, and what changes at μcrit is purely how the Coriolis force, which does no work and appears in no potential, organises the motion around that hill.

Five Points, Two Hills, One Surprise

Put a lighter body in orbit around two heavier ones that are themselves circling each other, and hunt for the places it could sit motionless relative to both. There are five. Three sit on the straight line running through the two heavy bodies. Two sit off to the side, each at the apex of an equilateral triangle whose base is the line joining the two masses. Lagrange found all five in 1772 solving a general three-body problem for prize money [1], apparently regarding the triangular ones as a mathematical curiosity with no physical system he expected to exhibit them. One hundred and thirty-four years later Max Wolf found an asteroid sitting at Jupiter's leading triangular point, and by now more than a million are known to sit at Jupiter's L4 and L5 [11]. Earth has two of its own, discovered in 2010 and 2021 [17]. The Kordylewski dust clouds, real cushions of interplanetary dust at the Moon's L5 point, were confirmed by imaging polarimetry in 2018 [12]. None of that was obvious from the algebra, and the reason it works is the part of this study we found funniest to compute.

Here is the joke, stated plainly, and the rest of the article is us taking it completely seriously. L4 and L5 sit on top of a hill of the gravitational-plus-centrifugal potential. Not a saddle, not a shoulder, an honest local maximum in the potential's own geometry. Drop a marble on a real hill and it rolls off. Drop a test particle at L4, under the right conditions, and it goes around and around forever, tracing a lopsided loop that never gets far. The trick has nothing to do with the shape of the hill. It is the Coriolis force, the same one that curves hurricanes and puts a spin on a dropped bolt aboard a rotating space station, and it appears nowhere in the potential we just used to build the hill, because it does no work and a force that does no work cannot show up in a potential energy at all. It shows up only in the equations of motion, as a term that depends on velocity, and it is entirely capable of steering a particle around the top of a hill it could never have climbed on its own.

We computed all five points for five mass ratios, checked every position against an independent numerical method and, where one exists, a closed-form or classical answer, classified the stability of each by the eigenvalues of the linearised dynamics, found the exact mass ratio at which the triangular points stop being stable, and then watched it happen directly by integrating actual trajectories and by throwing many random small nudges at the point and counting how many came back.

The Model We Actually Ran

Everything below is the circular restricted three-body problem (CR3BP): two bodies of mass 1−μ1-\mu and μ\mu, with μ=m2/(m1+m2)\mu = m_2/(m_1+m_2), on a fixed circular orbit about their barycentre, and a third body of no mass at all riding along in their gravity. Units are chosen so the separation of the two bodies is 1, their orbital angular rate is 1 (one orbit is 2π2\pi time units), and G(m1+m2)=1G(m_1+m_2)=1. In the frame that rotates with the two bodies, body 1 sits fixed at (−μ,0)(-\mu,0) and body 2 at (1−μ,0)(1-\mu,0), and the third body obeys

x′′−2y′=Ωxy′′+2x′=Ωyx'' - 2y' = \Omega_x \qquad y'' + 2x' = \Omega_y Ω(x,y)=12(x2+y2)+1−μr1+μr2,r1=(x+μ)2+y2,  r2=(x−1+μ)2+y2\Omega(x,y) = \tfrac12(x^2+y^2) + \frac{1-\mu}{r_1} + \frac{\mu}{r_2}, \quad r_1=\sqrt{(x+\mu)^2+y^2},\ \ r_2=\sqrt{(x-1+\mu)^2+y^2}

Ω\Omega is the effective potential: gravity from both masses, plus the centrifugal term from working in a rotating frame. The five Lagrange points are exactly its critical points, Ωx=Ωy=0\Omega_x=\Omega_y=0. The left-hand −2y′,+2x′-2y', +2x' terms are the Coriolis force. They do not come from Ω\Omega and they never will, because a force that always acts perpendicular to the velocity does zero work along any path, and a potential is by definition something whose gradient does the work. That single fact is the entire content of the section on why the corners survive, below, so we are going to keep repeating it until it stops sounding like a technicality.

The one conserved quantity of this system is the Jacobi constant, C=2Ω(x,y)−(x˙2+y˙2)C = 2\Omega(x,y) - (\dot x^2+\dot y^2), which we use throughout as a check on the numerical integration rather than as a claim about any Lagrange point's stability. The two are different things, and conflating them is a real and easy mistake (see the objection near the end).

The club has no telescope, no orbiting hardware, and no way to observe any of this. Every number below came out of a Python file we wrote, seeded once (seed 20260914, printed at the top of the run) and run start to finish on a laptop in 132 seconds. Real spacecraft near L1/L2 fly halo orbits and burn fuel to stay there; real Trojan-asteroid stability plays out over the age of the Solar System with resonances this two-body rotating frame does not contain. We say what our model leaves out as we go, and again in full near the end.

Locating the Points: Two Ways, Checked Against a Third

The three collinear points are roots of Ωx(x,0)=0\Omega_x(x,0)=0, a function with three singularities (harmless, they mark the two masses and one artefact of the algebra) and three zero crossings, at L1 between the masses, L2 beyond the smaller one, and L3 beyond the larger. We found them two ways. First, Newton's method from an initial guess built from the Hill radius rH=(μ/3)1/3r_H = (\mu/3)^{1/3}. Second, a coarse grid scan for sign changes of Ωx\Omega_x inside three intervals individually free of the two singularities, each refined by 80 rounds of bisection, a method that assumes nothing about where the root is beyond "not on top of a mass," and would raise loudly if it ever found the wrong number of roots in an interval. The two methods never disagreed by more than 2.2×10−16, machine epsilon for a number near 1, across every mass ratio we tried.

We then checked both against the classical small-μ\mu series (Szebehely 1967 [7]; Murray & Dermott 1999, §3.4–3.5 [8]):

xL1,L2=(1−μ)∓rH(1∓13rH−19rH2),xL3=−(1+512μ)x_{L1,L2} = (1-\mu) \mp r_H\Big(1 \mp \tfrac13 r_H - \tfrac19 r_H^2\Big), \qquad x_{L3} = -\Big(1+\tfrac{5}{12}\mu\Big)

Table 1 shows all three numbers side by side for every mass ratio. Agreement is excellent where the series is supposed to work (four to eight significant figures for Sun-Earth, Sun-Jupiter and Earth-Moon) and degrades smoothly as μ\mu grows, worst at the illustrative equal-mass case where the series is not even claiming to apply. That smooth, monotonic growth in disagreement (0.0002% at μ=10−5\mu=10^{-5} up to 2.5% at μ=0.1\mu=0.1, for L1) is itself a validation: an asymptotic series that failed randomly rather than growing smoothly with its own expansion parameter would mean we had the wrong formula, not merely one running out of range.

The triangular points needed no root-finder at all. Set r1=r2=1r_1=r_2=1 in the two equilibrium conditions and μ\mu cancels out of both completely: the only points satisfying that are the two apexes of equilateral triangles on the mass segment, (12−μ,±32)(\tfrac12-\mu, \pm\tfrac{\sqrt3}{2}), independent of μ\mu by construction. We still checked it numerically: for every mass ratio tested, r1r_1 and r2r_2 at the computed point equal 1 to better than 4×10−16, and ∣Ωx∣,∣Ωy∣|\Omega_x|,|\Omega_y| are of the same order. This is not really a test of our code so much as a test of whether we can type an equation correctly, and we are including it anyway because a study that only reports the tests it might plausibly fail is quietly editing itself.

systemμpointNewtonbisectionseries/exactNewton−seriesμ/μcrit (L4/L5)
Sun–Earth, μ = 3.0027×10−6
Sun-Earth3.003e-6L10.99002742670.99002742670.9900274339−7.16e-90.00008 stable
Sun-EarthL21.01003327391.01003327391.0100332677+6.20e-9
Sun-EarthL3−1.0000012511−1.0000012511−1.00000125110.00e+0
Sun–Jupiter, μ = 9.5344×10−4
Sun-Jupiter9.534e-4L10.93237587680.93237587680.9323913005−1.54e-50.02475 stable
Sun-JupiterL21.06882015881.06882015881.0688065474+1.36e-5
Sun-JupiterL3−1.0003972675−1.0003972675−1.0003972676+4.71e-11
Earth–Moon, μ = 1.2151×10−2
Earth-Moon1.215e-2L10.83691513300.83691513300.8373676934−4.53e-40.31543 stable
Earth-MoonL21.15568215981.15568215981.1552703297+4.12e-4
Earth-MoonL3−1.0050626452−1.0050626452−1.0050627434+9.82e-8
Pluto–Charon, μ = 1.0874×10−1
Pluto-Charon1.087e-1L10.59262439710.59262439710.6008560250−8.23e-32.82280 UNSTABLE
Pluto-CharonL21.26258391591.26258391591.2546853411+7.90e-3
Pluto-CharonL3−1.0452323446−1.0452323446−1.0453070172+7.47e-5
equal-mass pair (illustrative only), μ = 0.5
equal-mass0.5L10.00000000000.00000000000.0691484544−6.91e-212.98 UNSTABLE
equal-massL21.19840614461.19840614461.1327538337+6.57e-2
equal-massL3−1.1984061446−1.1984061446−1.2083333333+9.93e-3

The equal-mass row is worth pausing on: Newton lands L1 exactly on the barycentre, x=0x=0, which is correct by symmetry (two identical masses must balance a test particle exactly between them) and not a bug, even though the small-μ\mu series, asked to describe a system about as far from small-μ\mu as this problem allows, misses it by 0.069. We kept the row in the table for exactly that reason.

Figure 1 draws what the numbers in that table actually look like: the effective potential itself, sliced along the Earth-Moon rotating axis, with L1, L2 and L3 marked exactly where the curve turns over. Every one of them is a local minimum along this one slice, which is part of why they look deceptively well-like from this angle alone. The next section shows why the missing second direction is the one that turns each of them into a saddle instead.

L1L2L3EarthMoon-1-0.500.51x (co-rotating frame, Earth-Moon separation = 1)effective potential Ω(x,0)
Figure 1. The effective potential along the Earth-Moon rotating axis. L1, L2 and L3 sit exactly where the curve turns over between the (capped, for display) columns that rise without bound at Earth and the Moon. Each is a local minimum along this one-dimensional slice; the section on why the corners are stable shows why the missing second (off-axis) direction is the one that makes them saddles instead.

One more picture puts the whole geometry side by side: the two mass points, the two equilateral triangles they define, and where L1 through L5 actually land, for Earth-Moon (filled, solid) and Pluto-Charon (open, dashed). Pluto-Charon's triangle apex sits noticeably closer to the barycentre than Earth-Moon's, because a bigger secondary drags the whole equilateral construction toward itself, exactly as (12−μ,±32)(\tfrac12-\mu,\pm\tfrac{\sqrt3}{2}) says it must.

EarthMoonL1L2L3L4L5L4′L1′
Figure 2. Both mass points and all five Lagrange points for Earth-Moon (filled, solid triangle) and Pluto-Charon (open, dashed triangle), same rotating frame, same scale (mass separation normalised to 1 in each case). Primed labels mark the Pluto-Charon points; unprimed labels mark Earth-Moon's, which sit further from the barycentre because Earth-Moon's mass ratio is smaller.

Why the Corners Are Stable and the Line Never Is

Linearise the equations of motion around any equilibrium (x0,y0)(x_0,y_0): write x=x0+ξx=x_0+\xi, y=y0+ηy=y_0+\eta, and keep only the linear terms. Writing Ωxx,Ωxy,Ωyy\Omega_{xx},\Omega_{xy},\Omega_{yy} for the second partials of the effective potential at the point,

ξ′′−2η′=Ωxxξ+Ωxyηη′′+2ξ′=Ωxyξ+Ωyyη\xi'' - 2\eta' = \Omega_{xx}\xi + \Omega_{xy}\eta \qquad \eta'' + 2\xi' = \Omega_{xy}\xi + \Omega_{yy}\eta

as a 4-by-4 first-order system in (ξ,η,ξ′,η′)(\xi,\eta,\xi',\eta'), whose matrix is

A=(00100001ΩxxΩxy02ΩxyΩyy−20)A=\begin{pmatrix}0&0&1&0\\0&0&0&1\\ \Omega_{xx}&\Omega_{xy}&0&2\\ \Omega_{xy}&\Omega_{yy}&-2&0\end{pmatrix}

Look at where the Coriolis force lives in that matrix: the +2+2 and −2-2 in the bottom two rows, off the diagonal of the velocity block. It contributes nothing to Ωxx,Ωxy,Ωyy\Omega_{xx},\Omega_{xy},\Omega_{yy}, which describe only the curvature of the hill. It only ever couples ξ′\xi' into the η\eta equation and η′\eta' into the ξ\xi equation. At the collinear points, Ωxx>0\Omega_{xx}>0 and Ωyy<0\Omega_{yy}<0 always (we verified this numerically at all five mass ratios), so ΩxxΩyy<0\Omega_{xx}\Omega_{yy}<0, and the characteristic quartic λ4+(4−Ωxx−Ωyy)λ2+(ΩxxΩyy−Ωxy2)=0\lambda^4+(4-\Omega_{xx}-\Omega_{yy})\lambda^2+(\Omega_{xx}\Omega_{yy}-\Omega_{xy}^2)=0 is guaranteed a positive root for λ2\lambda^2, hence a real, positive, growing eigenvalue. No amount of Coriolis coupling saves a genuine saddle; a saddle stays a saddle regardless of the frame's rotation rate, and that is the whole story for L1, L2 and L3 at every mass ratio we tried, including the ones far outside the small-μ\mu regime the series expansion trusts.

At L4 (and L5, its mirror image) something different happens. Ωxx=34\Omega_{xx}=\tfrac34, Ωyy=94\Omega_{yy}=\tfrac94, both positive, and Ωxy=±334(1−2μ)\Omega_{xy}=\pm\tfrac{3\sqrt3}{4}(1-2\mu), for every μ\mu at once, which we obtained by evaluating the general second-derivative formula at the exact point (12−μ,±32)(\tfrac12-\mu,\pm\tfrac{\sqrt3}{2}) rather than by asserting it. Both diagonal curvatures positive means the hill really is a hill in every direction of the potential alone: were the Coriolis terms switched off, L4 would collapse outward in all directions and there would be no argument left to make. With them switched back on, the characteristic equation collapses to

λ4+λ2+274μ(1−μ)=0\lambda^4 + \lambda^2 + \tfrac{27}{4}\mu(1-\mu) = 0

which has two purely imaginary conjugate eigenvalue pairs, centres rather than saddles, meaning small displacements orbit the equilibrium rather than fleeing it, exactly when both roots of the quadratic in λ2\lambda^2 are real and negative, which happens exactly when μ(1−μ)≤127\mu(1-\mu)\le\frac{1}{27}. That is the entire mechanism. The hill never flattens. What changes is whether the rotating frame's own bookkeeping term is strong enough, relative to how steeply the hill falls away, to bend an outward slide into a closed loop before it gets anywhere. Below the critical mass ratio it is. Above it, the hill wins.

The same argument shows up directly in the eigenvalues themselves. Plot every one, for every point, for both example systems, in the complex plane: L1, L2 and L3 straddle the real axis always, a real positive eigenvalue paired with its negative plus a purely imaginary pair, the unmistakable signature of a saddle. L4 sits entirely on the imaginary axis for Earth-Moon, a pure centre, and lifts off it into all four quadrants for Pluto-Charon, the unmistakable signature of a spiral.

-4-2024Re(λ)Im(λ)
Figure 3. Eigenvalues of the linearised equations of motion in the complex plane, both mass ratios, four points each (L5 omitted, a mirror image of L4 in every case). Filled markers are Earth-Moon, open markers Pluto-Charon. L1/L2/L3 straddle the real axis in both systems, unstable regardless of μ. L4 sits on the imaginary axis for Earth-Moon (a centre) and lifts into all four quadrants for Pluto-Charon (a spiral, two growing directions and two decaying ones).

A Short Gazetteer of Who Is Actually Parked Where

The unstable points are the crowded ones, which sounds backwards until you remember that unstable just means you have to work to stay, and plenty of missions want exactly the vantage point that comes with that work. The James Webb Space Telescope sits near the unstable Sun-Earth L2, 1.5 million kilometres out, in a halo orbit around the mathematical point rather than on it, and needs periodic thruster burns to stay there [13]. Gaia keeps the same address for the same reason, mapping a billion stars from a perch it could not hold without fuel [14]. SOHO does the same job ninety degrees around the other way, at the equally unstable Sun-Earth L1, watching the Sun from a spot that lets it see it continuously, all year, forever, provided someone keeps paying the propellant bill. Every one of these missions is sitting, deliberately, on top of a saddle, because a saddle near Earth with a clear view of something valuable is worth the station-keeping.

The stable points get no traffic from us and plenty from nature. More than a million Jupiter Trojans are catalogued at Jupiter's L4 and L5, no fuel, no burns, just the arrangement this article spent its middle third explaining [11]. Earth has two known Trojans of its own, found in 2010 and 2021 [17], too faint and too far to matter to anyone's mission planning and entirely real regardless. The Kordylewski dust clouds sitting at the Moon's L5 point are fainter still, confirmed only in 2018 by imaging polarimetry sensitive enough to see interplanetary dust that has been quietly collecting there, held by nothing but the same Coriolis argument, for who knows how long [12]. Nobody launched any of it. Nobody station-keeps it. It is just where a hilltop, done right, keeps things.

The Arithmetic

Gascheau (1843) and Routh (1875) worked out the classical threshold from the quartic above [3][4]: μ(1−μ)=127\mu(1-\mu)=\frac{1}{27} is a quadratic in μ\mu, and the root inside (0,1)(0,1) is

μcrit=1−23/272\mu_{\text{crit}} = \frac{1-\sqrt{23/27}}{2}

Compute it. 23/27=0.851851851‾23/27 = 0.851851\overline{851}. Its square root is 0.9229582069...0.9229582069.... One minus that is 0.0770417931...0.0770417931.... Halve it: μcrit=0.0385208965...\mu_{\text{crit}}=0.0385208965.... Equivalently, as a primary-to-secondary mass ratio, (1−μcrit)/μcrit=24.959936(1-\mu_{\text{crit}})/\mu_{\text{crit}} = 24.959936, commonly rounded to the "24.96 to 1" figure quoted in the literature.

We found the same number three ways in code: bisecting the discriminant 1−27μ(1−μ)1-27\mu(1-\mu) directly (agrees with the hand arithmetic above to 0 at printed precision, since both sides are the same closed-form expression evaluated in float64), and bisecting on the actual numerically diagonalised 4-by-4 matrix with no L4-specific algebra assumed anywhere, landing at 0.0385208525, a difference of 4.4×10−8. That residual is not noise we shrugged off: exactly at μcrit\mu_{\text{crit}} the matrix has a repeated eigenvalue, and it is a standard fact of matrix perturbation theory that an eigenvalue near a degeneracy moves as the square root of a perturbation in the matrix rather than linearly, so an error of order the square root of machine epsilon (2−52≈1.5×10−8\sqrt{2^{-52}}\approx1.5\times10^{-8}) is exactly the size this method should produce, and 4.4×10−8 matches that prediction well enough that we stopped looking for a bug.

0.0385209classical μcrit, Gascheau/Routh
0.0385209our discriminant bisection
0.0385209our eigenvalue bisection (5 s.f.)
24.96 : 1primary : secondary mass ratio

Figure 4 shows the threshold as a curve rather than a single number: the largest real part among L4's four eigenvalues, scanned over 240 mass ratios from 0 to 2.9× μcrit\mu_{\text{crit}}, computed by the same general eigenvalue routine used everywhere else in this study, with no shortcut taken for this figure specifically.

μcritEarth-MoonPluto-Charon(0.393, off scale)00.511.522.500.10.20.30.4mass ratio μ / μcritmax Re(λ) at L4
Figure 4. The largest real part of the L4 eigenvalues, scanned over mass ratio and computed by the same general eigenvalue routine used throughout. It is exactly zero for every μ below the classical threshold (to eigenvalue-solver precision, <10−15) and rises smoothly above it. Earth-Moon sits deep in the flat part; Pluto-Charon's growth rate, 0.393 per time unit, is off the top of this axis entirely.

Working Notes From the Club Table

Meeting 1
First attempt at the grid-scan root finder used one big array from −3 to 4 with the two singular points simply masked out. It found four "roots" for every system instead of three. Spent forty minutes convinced the physics was wrong before noticing the masked-out gap made two genuinely unrelated points look adjacent, so the sign change across the gap got treated as a real root. Fixed by scanning three separate intervals, one per point, each provably free of a singularity. The bug produced a spurious root near x=−μx=-\mu every single time, which in hindsight was the tell.

Meeting 2
Single-trajectory integration near Pluto-Charon's L3 reported a Jacobi-constant drift of 500-something, which would mean the integrator was worthless. Traced it to exactly what you would guess once you say it out loud: an unstable trajectory that has already escaped keeps flying near a mass point, the fixed step size stops resolving the close approach, and by then we do not need that part of the trajectory anyway. Now we only score the drift up to a short margin past the first threshold crossing. Drift there is <4×10−13 relative, every case, full stop.

Meeting 2, later
First Monte Carlo design used a starting displacement of 10−5 and a threshold at 1000× that, matching the single-trajectory demo. At the mass ratios near the threshold, where the whole point of the experiment lives, essentially nothing escaped inside a runtime we could afford: the growth rate there is small enough that even a trajectory aligned as well as possible with the unstable direction could not multiply a displacement by 1000 in the time we had budgeted. Loosened the ratio to 100× instead. That one number is the entire reason Figure 5 below shows a transition instead of a flat line at zero.

Meeting 3
Ran the Monte Carlo at μ exactly equal to the classical critical value out of curiosity, expecting close to 0% escape since the linear theory says the real part of every eigenvalue is exactly zero there. Got 80.0%. Not a bug: exactly at the critical mass ratio the 4-by-4 matrix has a repeated eigenvalue with only one independent eigenvector, a defective matrix, and the general solution of a defective linear system contains a term proportional to tt times the oscillation, not just the oscillation. A term that grows linearly in time still eventually crosses any fixed threshold. So "exactly at the edge" is, in this particular and slightly delightful sense, not actually a safe place to sit. We had not planned to discover that; the code found it for us.

Meeting 4
Seed fixed at 20260914, printed at the top of the run, never touched again. Every number in this article comes from one execution of the file as it now stands, 132.3 seconds on a laptop.

Watching It Happen: Test Particles in Time

Eigenvalues describe what an infinitesimal displacement does. To see the actual claim in action we integrated real trajectories: a fourth-order Runge-Kutta scheme we wrote, fixed step Δt=0.002\Delta t=0.002, started at each of the five points plus a perturbation of 10−510^{-5} in a single fixed generic direction (normalised (1,1,0.5,−0.5)(1,1,0.5,-0.5) in (ξ,η,ξ′,η′)(\xi,\eta,\xi',\eta')), run for 40 time units (about 6.4 orbital periods), and scored by the first time the distance from the starting point exceeds 10−210^{-2}, one thousand times the initial nudge.

systempointoutcomeescape timemax distance reachedJacobi drift |ΔC/C0|
Earth-Moon (μ/μcrit = 0.315, below threshold)
Earth-MoonL1escaped2.692.70e-13.90e-15
Earth-MoonL2escaped3.606.07e+02.80e-16
Earth-MoonL3escaped27.349.53e-22.95e-16
Earth-MoonL4containedn/a1.99e-42.97e-16
Earth-MoonL5containedn/a5.76e-52.97e-16
Pluto-Charon (μ/μcrit = 2.823, above threshold)
Pluto-CharonL1escaped2.335.27e-13.43e-15
Pluto-CharonL2escaped4.261.87e+12.55e-16
Pluto-CharonL3escaped11.745.48e+22.86e-16
Pluto-CharonL4escaped15.012.94e+03.06e-16
Pluto-CharonL5escaped20.761.46e+03.06e-16

Same code, same perturbation direction, same magnitude, same threshold. L1, L2 and L3 escape in both systems, usually within a handful of time units. L4 and L5 sit through the entire 40 time units in Earth-Moon and leave, respectively, at t=15.0t=15.0 and t=20.8t=20.8 in Pluto-Charon. The only thing that changed between the two blocks of that table is which side of 0.0385 the mass ratio falls on. The Jacobi-constant column is the honesty check on the integrator itself: every value is at the level of double-precision roundoff, which means the escapes above are the dynamics doing what the dynamics do, not an artefact of an integrator that cannot add.

How Sure Are We? A Monte Carlo Coin Flip at the Edge

One trajectory per point shows that two outcomes exist. It cannot show where the line between them is, because whether one particular fixed perturbation escapes by a fixed time depends on how much of that one direction happens to line up with the growing mode, and a single direction is a sample size of one. So we drew many random directions instead: for each of ten mass ratios from deep in the stable region to well past threshold, 1,500 independent unit vectors in the four-dimensional (ξ,η,ξ′,η′)(\xi,\eta,\xi',\eta') phase space, each scaled to a displacement of 10−410^{-4} and integrated for 70 time units, scored as escaped if the distance from L4 ever exceeded 10−210^{-2} (a threshold-to-perturbation ratio of 100, chosen, per the working notes, because the original 1000× ratio left the interesting mass ratios at 0% for lack of runtime). Every draw came from one seeded generator (seed 20260914), and it is the only randomness anywhere in this study.

mass ratioμμ/μcritescaped / Np̂SE
Earth-Moon1.2151e-20.3150 / 15000.00000.0000
0.50×μcrit1.9260e-20.5000 / 15000.00000.0000
0.85×μcrit3.2743e-20.8500 / 15000.00000.0000
0.95×μcrit3.6595e-20.9500 / 15000.00000.0000
0.98×μcrit3.7750e-20.980282 / 15000.18800.0101
1.00×μcrit3.8521e-21.0001200 / 15000.80000.0103
1.02×μcrit3.9291e-21.0201477 / 15000.98470.0032
1.05×μcrit4.0447e-21.0501499 / 15000.99930.0007
1.20×μcrit4.6225e-21.2001500 / 15001.00000.0000
Pluto-Charon1.0874e-12.8231500 / 15001.00000.0000

The escape probability climbs from a clean 0% to a clean 100% inside a factor of roughly 1.2 in mass ratio, and the climb is centred almost exactly on μcrit\mu_{\text{crit}}, not because we tuned it there but because that is where the linear growth rate genuinely turns on. The one entry that looks odd, 80.0% exactly at the classical threshold, is the defective-matrix effect from the working notes: linear-in-time growth from a repeated eigenvalue with a single eigenvector, not exponential growth, but growth all the same, and over 70 time units enough of the 1,500 random directions had a large enough coefficient on that term to cross the threshold. We did not expect that number and do not think we would have found it without actually running the integration; see the objection below.

final p̂ = 0.188110100100000.20.40.6trials accumulated n (log scale)running escape probability p̂
Figure 5. The running Monte Carlo estimate of L4's escape probability at 0.98×μcrit, plotted against the number of trials drawn so far (log scale), with a shaded band of ±1 standard error. The estimate visibly stops moving by a few hundred trials, and the band narrows the way 1/n1/\sqrt{n} predicts, converging to 18.8% ± 1.0% after all 1,500 trials.

The Strongest Objection We Could Make

The objection runs like this. Every number in this study is a solution of equations that were fully solved algebraically before we were born. The characteristic quartic at L4/L5 has a closed form; the Gascheau/Routh threshold has a closed form; the collinear-point series is a hundred and fifty years old. Running numerical code against a problem with a known answer checks whether we can type the equations correctly and whether numpy's linear algebra library works, and neither of those was ever in serious doubt. A club member who had memorised Szebehely's book could have predicted every number in Table 1 to five digits with a pocket calculator, before we wrote a line of code.

We think this objection is largely correct, and we think two parts of the study survive it anyway.

The triangular-point verification (the r1=r2=1 rows above, and the closed-form quartic itself) genuinely is only a code check, and we said so where it happens rather than dressing it up as a discovery, because it was never one. Its value is in what it would have caught: a sign error in Ωxy\Omega_{xy}, say, or a mismatched convention for which body sits at −μ-\mu versus 1−μ1-\mu, would have shown up immediately as r1 ≠ r2 or a nonzero Ωx\Omega_x at the "L4" the code reported. It caught nothing, which is itself informative and not nothing.

The Monte Carlo result at exactly μcrit\mu_{\text{crit}}, 80.0% escaped rather than the roughly-zero we expected from "the real parts are exactly zero there," was not predictable from the closed-form linear theory alone. Nobody handed us that number. It required actually running the nonlinear integration long enough, and then noticing that a defective matrix's solution includes a secular tt-times-oscillation term that a glance at "purely imaginary eigenvalues" will not tell you about. That is squarely the kind of thing a numerical experiment finds and a closed-form derivation, read too quickly, will not. We would defend that result as the actual contribution of this study, small as it is against a body of hundred-year-old mathematics.

A second objection resists us less comfortably: everything here is the idealised two-body-plus-massless-third-body planar circular problem, and every real system we plugged into it (Earth-Moon, Sun-Earth, Sun-Jupiter, Pluto-Charon) has nonzero eccentricity, other massive bodies nearby, and (for the ones with spacecraft) solar radiation pressure that this model has never heard of. We address that directly below rather than argue our way past it here.

What This Does Not Tell You, and Where We'd Get a Different Answer

The threshold μcrit=0.03852\mu_{\text{crit}}=0.03852 is exact for the model we actually solved: two bodies, exactly circular, exactly two-dimensional, a third body of exactly zero mass. Real orbits are not exactly circular. The Earth-Moon system has an orbital eccentricity around 0.055 and the true equations acquire time-dependent terms this frame does not carry at all; the elliptic restricted three-body problem is a different, harder problem, and we did not solve it. Nothing here should be read as a claim about the actual, eccentric Earth-Moon system's Lagrange-point stability beyond "the circular approximation is probably a reasonable starting point," which is a weaker claim than the numbers in this article might look like they are making.

Planar versus three-dimensional. We integrated in the orbital plane only. Real spacecraft near L1 and L2 fly three-dimensional halo orbits specifically because the planar problem's L1/L2 instability is already fast enough to matter; nothing about our stability classification for L4/L5 would change by adding a third dimension (the out-of-plane linearisation decouples from the in-plane one at these points and is separately, mildly stable), but we did not verify that decoupling ourselves and are reporting it as a textbook fact rather than a result of this study [8].

Fixed-step, non-symplectic integration. A proper long-term stability study of a formally stable point would use a symplectic integrator, which conserves a shadow Hamiltonian exactly rather than approximately, precisely because generic integrators like ours slowly leak or gain energy over very long runs even when short-run Jacobi-constant drift looks fine. Our runs are short (40 to 70 time units, a handful to a few dozen orbital periods) and the drift we measured is at machine precision throughout, so we do not think this changes any conclusion in this article, but a study of Trojan-asteroid stability over the age of the Solar System (billions of periods) absolutely could not use our integrator and expect the same answer.

Two bodies only. The real Solar System is not a two-body system anywhere. Jupiter's Trojans are perturbed by Saturn; the Moon is perturbed by the Sun to a degree that matters for real trajectory design. Had we included even one more massive body, μcrit\mu_{\text{crit}} itself would stop being a fixed number and become a statement about resonances and secular perturbations, which is most of what a professional mission-design textbook on this subject actually is [9]. We picked the two-body idealisation because it is the version with an exact answer to check ourselves against, which was the point of this particular study; a club project on Trojan-asteroid long-term stability would need the fuller problem and could not reuse this code as-is.

The threshold-to-perturbation ratio in the Monte Carlo. We chose 100× after the first choice (1000×) produced an uninformative flat zero near the interesting mass ratios, as the working notes admit. A different ratio would move every escape-probability number in the Monte Carlo table, though not the qualitative picture (a sigmoid centred near μcrit\mu_{\text{crit}}) or the location of μcrit\mu_{\text{crit}} itself, which came from the eigenvalues and not from the Monte Carlo at all.

Reproducing this

Everything in this article comes from one file and one command. You need Python 3.12 and numpy; no other package is used anywhere in the analysis.

python lagrange-points.py > lagrange-points-output.txt

Expect about two minutes on a modern laptop; ours took 132.3 s on numpy 2.4.2 and Python 3.12.3. The master seed is 20260914, hard-coded at the top of the file and used for nothing except the Monte Carlo section's random directions, so the whole output is deterministic and every number in this article can be reproduced exactly by re-running the file. Numbers differing from ours by more than the printed standard errors, or by more than floating-point roundoff in the deterministic sections, mean something is wrong, and we would like to hear about it. The interactive model runs the same effective-potential and stability calculations live in your browser, with sliders for mass ratio and for the size and direction of a perturbation at any of the five points.

References

Reproduce these numbers

From the analysis folder, run:

python lagrange-points.py > lagrange-points-output.txt
Analysis code
lagrange-points.py37 kB