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
- Files
- LaTeX source22 kB
- Analysis code37 kB
- Raw output18 kB
- Interactive model
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 and , with , 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 time units), and . In the frame that rotates with the two bodies, body 1 sits fixed at and body 2 at , and the third body obeys
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, . The left-hand terms are the Coriolis force. They do not come from 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, , 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 , 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 . Second, a coarse grid scan for sign changes of 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- series (Szebehely 1967 [7]; Murray & Dermott 1999, §3.4–3.5 [8]):
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 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 up to 2.5% at , 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 in the two equilibrium conditions and cancels out of both completely: the only points satisfying that are the two apexes of equilateral triangles on the mass segment, , independent of by construction. We still checked it numerically: for every mass ratio tested, and at the computed point equal 1 to better than 4×10−16, and 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 | μ | point | Newton | bisection | series/exact | Newton−series | μ/μcrit (L4/L5) |
|---|---|---|---|---|---|---|---|
| Sun–Earth, μ = 3.0027×10−6 | |||||||
| Sun-Earth | 3.003e-6 | L1 | 0.9900274267 | 0.9900274267 | 0.9900274339 | −7.16e-9 | 0.00008 stable |
| Sun-Earth | L2 | 1.0100332739 | 1.0100332739 | 1.0100332677 | +6.20e-9 | ||
| Sun-Earth | L3 | −1.0000012511 | −1.0000012511 | −1.0000012511 | 0.00e+0 | ||
| Sun–Jupiter, μ = 9.5344×10−4 | |||||||
| Sun-Jupiter | 9.534e-4 | L1 | 0.9323758768 | 0.9323758768 | 0.9323913005 | −1.54e-5 | 0.02475 stable |
| Sun-Jupiter | L2 | 1.0688201588 | 1.0688201588 | 1.0688065474 | +1.36e-5 | ||
| Sun-Jupiter | L3 | −1.0003972675 | −1.0003972675 | −1.0003972676 | +4.71e-11 | ||
| Earth–Moon, μ = 1.2151×10−2 | |||||||
| Earth-Moon | 1.215e-2 | L1 | 0.8369151330 | 0.8369151330 | 0.8373676934 | −4.53e-4 | 0.31543 stable |
| Earth-Moon | L2 | 1.1556821598 | 1.1556821598 | 1.1552703297 | +4.12e-4 | ||
| Earth-Moon | L3 | −1.0050626452 | −1.0050626452 | −1.0050627434 | +9.82e-8 | ||
| Pluto–Charon, μ = 1.0874×10−1 | |||||||
| Pluto-Charon | 1.087e-1 | L1 | 0.5926243971 | 0.5926243971 | 0.6008560250 | −8.23e-3 | 2.82280 UNSTABLE |
| Pluto-Charon | L2 | 1.2625839159 | 1.2625839159 | 1.2546853411 | +7.90e-3 | ||
| Pluto-Charon | L3 | −1.0452323446 | −1.0452323446 | −1.0453070172 | +7.47e-5 | ||
| equal-mass pair (illustrative only), μ = 0.5 | |||||||
| equal-mass | 0.5 | L1 | 0.0000000000 | 0.0000000000 | 0.0691484544 | −6.91e-2 | 12.98 UNSTABLE |
| equal-mass | L2 | 1.1984061446 | 1.1984061446 | 1.1327538337 | +6.57e-2 | ||
| equal-mass | L3 | −1.1984061446 | −1.1984061446 | −1.2083333333 | +9.93e-3 | ||
The equal-mass row is worth pausing on: Newton lands L1 exactly on the barycentre, , which is correct by symmetry (two identical masses must balance a test particle exactly between them) and not a bug, even though the small- series, asked to describe a system about as far from small- 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.
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 says it must.
Why the Corners Are Stable and the Line Never Is
Linearise the equations of motion around any equilibrium : write , , and keep only the linear terms. Writing for the second partials of the effective potential at the point,
as a 4-by-4 first-order system in , whose matrix is
Look at where the Coriolis force lives in that matrix: the and in the bottom two rows, off the diagonal of the velocity block. It contributes nothing to , which describe only the curvature of the hill. It only ever couples into the equation and into the equation. At the collinear points, and always (we verified this numerically at all five mass ratios), so , and the characteristic quartic is guaranteed a positive root for , 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- regime the series expansion trusts.
At L4 (and L5, its mirror image) something different happens. , , both positive, and , for every at once, which we obtained by evaluating the general second-derivative formula at the exact point 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
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 are real and negative, which happens exactly when . 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.
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]: is a quadratic in , and the root inside is
Compute it. . Its square root is . One minus that is . Halve it: . Equivalently, as a primary-to-secondary mass ratio, , 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 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 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 () 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.
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× , computed by the same general eigenvalue routine used everywhere else in this study, with no shortcut taken for this figure specifically.
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 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 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 , started at each of the five points plus a perturbation of in a single fixed generic direction (normalised in ), run for 40 time units (about 6.4 orbital periods), and scored by the first time the distance from the starting point exceeds , one thousand times the initial nudge.
| system | point | outcome | escape time | max distance reached | Jacobi drift |ΔC/C0| |
|---|---|---|---|---|---|
| Earth-Moon (μ/μcrit = 0.315, below threshold) | |||||
| Earth-Moon | L1 | escaped | 2.69 | 2.70e-1 | 3.90e-15 |
| Earth-Moon | L2 | escaped | 3.60 | 6.07e+0 | 2.80e-16 |
| Earth-Moon | L3 | escaped | 27.34 | 9.53e-2 | 2.95e-16 |
| Earth-Moon | L4 | contained | n/a | 1.99e-4 | 2.97e-16 |
| Earth-Moon | L5 | contained | n/a | 5.76e-5 | 2.97e-16 |
| Pluto-Charon (μ/μcrit = 2.823, above threshold) | |||||
| Pluto-Charon | L1 | escaped | 2.33 | 5.27e-1 | 3.43e-15 |
| Pluto-Charon | L2 | escaped | 4.26 | 1.87e+1 | 2.55e-16 |
| Pluto-Charon | L3 | escaped | 11.74 | 5.48e+2 | 2.86e-16 |
| Pluto-Charon | L4 | escaped | 15.01 | 2.94e+0 | 3.06e-16 |
| Pluto-Charon | L5 | escaped | 20.76 | 1.46e+0 | 3.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 and 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 phase space, each scaled to a displacement of and integrated for 70 time units, scored as escaped if the distance from L4 ever exceeded (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 | μ | μ/μcrit | escaped / N | p̂ | SE |
|---|---|---|---|---|---|
| Earth-Moon | 1.2151e-2 | 0.315 | 0 / 1500 | 0.0000 | 0.0000 |
| 0.50×μcrit | 1.9260e-2 | 0.500 | 0 / 1500 | 0.0000 | 0.0000 |
| 0.85×μcrit | 3.2743e-2 | 0.850 | 0 / 1500 | 0.0000 | 0.0000 |
| 0.95×μcrit | 3.6595e-2 | 0.950 | 0 / 1500 | 0.0000 | 0.0000 |
| 0.98×μcrit | 3.7750e-2 | 0.980 | 282 / 1500 | 0.1880 | 0.0101 |
| 1.00×μcrit | 3.8521e-2 | 1.000 | 1200 / 1500 | 0.8000 | 0.0103 |
| 1.02×μcrit | 3.9291e-2 | 1.020 | 1477 / 1500 | 0.9847 | 0.0032 |
| 1.05×μcrit | 4.0447e-2 | 1.050 | 1499 / 1500 | 0.9993 | 0.0007 |
| 1.20×μcrit | 4.6225e-2 | 1.200 | 1500 / 1500 | 1.0000 | 0.0000 |
| Pluto-Charon | 1.0874e-1 | 2.823 | 1500 / 1500 | 1.0000 | 0.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 , 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.
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 , say, or a mismatched convention for which body sits at versus , would have shown up immediately as r1 ≠ r2 or a nonzero at the "L4" the code reported. It caught nothing, which is itself informative and not nothing.
The Monte Carlo result at exactly , 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 -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 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, 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 ) or the location of 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.txtExpect 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
- Raw output
- lagrange-points-output.txt18 kB
Lagrange, J.-L. (1772). Essai sur le problème des trois corps. Prize essay, Académie Royale des Sciences; reprinted in Œuvres de Lagrange, vol. 6.
Jacobi, C. G. J. (1836). Sur le mouvement d'un point et sur un cas particulier du problème des trois corps. Comptes Rendus de l'Académie des Sciences de Paris 3, 59–61.
Gascheau, G. (1843). Examen d'une classe d'équations différentielles et application à un cas particulier du problème des trois corps. Comptes Rendus de l'Académie des Sciences 16, 393–394.
Routh, E. J. (1875). On Laplace's three particles, with a supplement on the stability of steady motion. Proceedings of the London Mathematical Society s1-6(1), 86–97.
Poincaré, H. (1892). Les Méthodes Nouvelles de la Mécanique Céleste, vol. 1. Gauthier-Villars, Paris.
Barrow-Green, J. (2008). The Three-Body Problem. In T. Gowers, J. Barrow-Green & I. Leader (eds.), The Princeton Companion to Mathematics, pp. 726–728. Princeton University Press.
Szebehely, V. (1967). Theory of Orbits: The Restricted Problem of Three Bodies. Academic Press, New York.
Murray, C. D. & Dermott, S. F. (1999). Solar System Dynamics. Cambridge University Press. ISBN 978-0-521-57597-3.
Koon, W. S., Lo, M. W., Marsden, J. E. & Ross, S. D. (2011). Dynamical Systems, the Three-Body Problem and Space Mission Design. Marsden Books.
Gutzwiller, M. C. (1998). Moon-Earth-Sun: the oldest three-body problem. Reviews of Modern Physics 70(2), 589–639. doi:10.1103/RevModPhys.70.589
Krishnaswami, G. S. & Senapati, H. (2019). An introduction to the classical three-body problem. Resonance 24, 87–114. doi:10.1007/s12045-019-0760-1
Slíz-Balogh, J., Barta, A. & Horváth, G. (2018). Celestial mechanics and polarization optics of the Kordylewski dust cloud in the Earth-Moon Lagrange point L5, Part I. Monthly Notices of the Royal Astronomical Society 480(4), 5550–5559. doi:10.1093/mnras/sty2049
NASA (2024). Orbit: Webb Space Telescope. science.nasa.gov/mission/webb/orbit, accessed September 2026.
Prusti, T. et al. (Gaia Collaboration) (2016). The Gaia mission. Astronomy & Astrophysics 595, A1. doi:10.1051/0004-6361/201629272
Pitjeva, E. V. & Standish, E. M. (2009). Proposals for the masses of the three largest asteroids, the Moon-Earth mass ratio and the astronomical unit. Celestial Mechanics and Dynamical Astronomy 103(4), 365–372. doi:10.1007/s10569-009-9203-8
Brozović, M. & Jacobson, R. A. (2024). Post-New-Horizons orbits and masses for the satellites of Pluto. The Astronomical Journal 167(6), 256. doi:10.3847/1538-3881/ad39f0
Hui, M.-T., Wiegert, P. A., Tholen, D. J. & Föhring, D. (2021). The second Earth Trojan, 2020 XL5. The Astrophysical Journal Letters 922(2), L25. doi:10.3847/2041-8213/ac37bf