============================================================================== SECTION 0. Mass ratios used in this study ============================================================================== master seed for every random draw in this file: 20260914 mu[Sun-Earth ] = 3.0027261364e-06 mu[Sun-Jupiter ] = 9.5344219113e-04 mu[Earth-Moon ] = 1.2150584133e-02 mu[Pluto-Charon ] = 1.0873684139e-01 mu[equal-mass pair (illustrative, not real) ] = 5.0000000000e-01 ============================================================================== SECTION 2. Collinear points: Newton vs grid-bisection vs classical series ============================================================================== system point Newton bisection series Newton-series Sun-Earth L1 0.9900274267 0.9900274267 0.9900274339 -7.164e-09 Sun-Earth L2 1.0100332739 1.0100332739 1.0100332677 6.198e-09 Sun-Earth L3 -1.0000012511 -1.0000012511 -1.0000012511 0.000e+00 Newton-bisection agreement (max |diff| over L1,L2,L3): 2.220e-16 Sun-Jupiter L1 0.9323758768 0.9323758768 0.9323913005 -1.542e-05 Sun-Jupiter L2 1.0688201588 1.0688201588 1.0688065474 1.361e-05 Sun-Jupiter L3 -1.0003972675 -1.0003972675 -1.0003972676 4.713e-11 Newton-bisection agreement (max |diff| over L1,L2,L3): 1.110e-16 Earth-Moon L1 0.8369151330 0.8369151330 0.8373676934 -4.526e-04 Earth-Moon L2 1.1556821598 1.1556821598 1.1552703297 4.118e-04 Earth-Moon L3 -1.0050626452 -1.0050626452 -1.0050627434 9.819e-08 Newton-bisection agreement (max |diff| over L1,L2,L3): 2.220e-16 Pluto-Charon L1 0.5926243971 0.5926243971 0.6008560250 -8.232e-03 Pluto-Charon L2 1.2625839159 1.2625839159 1.2546853411 7.899e-03 Pluto-Charon L3 -1.0452323446 -1.0452323446 -1.0453070172 7.467e-05 Newton-bisection agreement (max |diff| over L1,L2,L3): 1.110e-16 equal-mass pair (illustrative, not real) L1 -0.0000000000 -0.0000000000 0.0691484544 -6.915e-02 equal-mass pair (illustrative, not real) L2 1.1984061446 1.1984061446 1.1327538337 6.565e-02 equal-mass pair (illustrative, not real) L3 -1.1984061446 -1.1984061446 -1.2083333333 9.927e-03 Newton-bisection agreement (max |diff| over L1,L2,L3): 2.220e-16 The Earth-Moon and Sun-Earth rows agree with the classical series to four or more significant figures. The disagreement grows with mu (worst at the illustrative mu=0.5 case, where the series is not even meant to apply: it is a small-mu expansion), which is the expected signature of a genuine asymptotic series rather than a coincidence or a bug -- Section 4b below quantifies that growth. ============================================================================== SECTION 2b. Triangular points: exact geometry check ============================================================================== system point x y r1 r2 |Omega_x| |Omega_y| Sun-Earth L4 0.49999700 0.86602540 1.0000000000000 1.0000000000000 1.35e-16 2.71e-16 Sun-Earth L5 0.49999700 -0.86602540 1.0000000000000 1.0000000000000 1.35e-16 2.71e-16 Sun-Jupiter L4 0.49904656 0.86602540 1.0000000000000 1.0000000000000 1.57e-16 3.66e-16 Sun-Jupiter L5 0.49904656 -0.86602540 1.0000000000000 1.0000000000000 1.57e-16 3.66e-16 Earth-Moon L4 0.48784942 0.86602540 1.0000000000000 1.0000000000000 1.98e-16 3.68e-16 Earth-Moon L5 0.48784942 -0.86602540 1.0000000000000 1.0000000000000 1.98e-16 3.68e-16 Pluto-Charon L4 0.39126316 0.86602540 1.0000000000000 1.0000000000000 1.87e-16 3.33e-16 Pluto-Charon L5 0.39126316 -0.86602540 1.0000000000000 1.0000000000000 1.87e-16 3.33e-16 equal-mass pair (illustrative, not real) L4 0.00000000 0.86602540 1.0000000000000 1.0000000000000 0.00e+00 3.33e-16 equal-mass pair (illustrative, not real) L5 0.00000000 -0.86602540 1.0000000000000 1.0000000000000 0.00e+00 3.33e-16 r1 = r2 = 1 to machine precision in every case confirms the equilateral triangle, independent of mu, exactly as the algebra predicts: set r1=r2=1 in the two equilibrium conditions and mu cancels out completely, which is why this holds for every mass ratio at once rather than needing to be checked mu by mu. ============================================================================== SECTION 3. Eigenvalues of the linearised equations of motion ============================================================================== Sun-Earth (mu = 3.002726e-06) L1: eigenvalues [+2.5326+0.0000i, +0.0000+2.0864i, +0.0000-2.0864i, -2.5326+0.0000i] -> UNSTABLE (max Re = 2.533e+00) L2: eigenvalues [+2.4844+0.0000i, +0.0000+2.0571i, +0.0000-2.0571i, -2.4844+0.0000i] -> UNSTABLE (max Re = 2.484e+00) L3: eigenvalues [+0.0028+0.0000i, -0.0000+1.0000i, -0.0000-1.0000i, -0.0028+0.0000i] -> UNSTABLE (max Re = 2.808e-03) L4: eigenvalues [+0.0000+1.0000i, +0.0000-1.0000i, -0.0000+0.0045i, -0.0000-0.0045i] -> stable (linear) (max Re = 7.216e-16) L5: eigenvalues [+0.0000-0.0045i, +0.0000+0.0045i, -0.0000-1.0000i, -0.0000+1.0000i] -> stable (linear) (max Re = 7.481e-16) Sun-Jupiter (mu = 9.534422e-04) L1: eigenvalues [+2.6811+0.0000i, -0.0000+2.1777i, -0.0000-2.1777i, -2.6811+0.0000i] -> UNSTABLE (max Re = 2.681e+00) L2: eigenvalues [+2.3521+0.0000i, +0.0000+1.9772i, +0.0000-1.9772i, -2.3521+0.0000i] -> UNSTABLE (max Re = 2.352e+00) L3: eigenvalues [+0.0500+0.0000i, +0.0000+1.0008i, +0.0000-1.0008i, -0.0500+0.0000i] -> UNSTABLE (max Re = 5.001e-02) L4: eigenvalues [+0.0000+0.9968i, +0.0000-0.9968i, -0.0000+0.0804i, -0.0000-0.0804i] -> stable (linear) (max Re = 9.159e-16) L5: eigenvalues [+0.0000-0.0804i, +0.0000+0.0804i, -0.0000-0.9968i, -0.0000+0.9968i] -> stable (linear) (max Re = 1.003e-15) Earth-Moon (mu = 1.215058e-02) L1: eigenvalues [+2.9321+0.0000i, -0.0000+2.3344i, -0.0000-2.3344i, -2.9321+0.0000i] -> UNSTABLE (max Re = 2.932e+00) L2: eigenvalues [+2.1587+0.0000i, -0.0000+1.8626i, -0.0000-1.8626i, -2.1587+0.0000i] -> UNSTABLE (max Re = 2.159e+00) L3: eigenvalues [+0.1779+0.0000i, +0.0000+1.0104i, +0.0000-1.0104i, -0.1779+0.0000i] -> UNSTABLE (max Re = 1.779e-01) L4: eigenvalues [+0.0000-0.2982i, +0.0000+0.2982i, -0.0000-0.9545i, -0.0000+0.9545i] -> stable (linear) (max Re = 3.391e-16) L5: eigenvalues [+0.0000+0.9545i, +0.0000-0.9545i, -0.0000+0.2982i, -0.0000-0.2982i] -> stable (linear) (max Re = 4.857e-17) Pluto-Charon (mu = 1.087368e-01) L1: eigenvalues [+3.4119+0.0000i, +0.0000+2.6411i, +0.0000-2.6411i, -3.4119+0.0000i] -> UNSTABLE (max Re = 3.412e+00) L2: eigenvalues [+1.7884+0.0000i, +0.0000+1.6519i, +0.0000-1.6519i, -1.7884+0.0000i] -> UNSTABLE (max Re = 1.788e+00) L3: eigenvalues [+0.5226+0.0000i, +0.0000+1.0831i, +0.0000-1.0831i, -0.5226+0.0000i] -> UNSTABLE (max Re = 5.226e-01) L4: eigenvalues [+0.3929-0.8090i, +0.3929+0.8090i, -0.3929-0.8090i, -0.3929+0.8090i] -> UNSTABLE (max Re = 3.929e-01) L5: eigenvalues [+0.3929+0.8090i, +0.3929-0.8090i, -0.3929+0.8090i, -0.3929-0.8090i] -> UNSTABLE (max Re = 3.929e-01) equal-mass pair (illustrative, not real) (mu = 5.000000e-01) L1: eigenvalues [+3.7833+0.0000i, -0.0000+2.8834i, -0.0000-2.8834i, -3.7833+0.0000i] -> UNSTABLE (max Re = 3.783e+00) L2: eigenvalues [+1.1557+0.0000i, +0.0000+1.3289i, +0.0000-1.3289i, -1.1557+0.0000i] -> UNSTABLE (max Re = 1.156e+00) L3: eigenvalues [+1.1557+0.0000i, +0.0000+1.3289i, +0.0000-1.3289i, -1.1557+0.0000i] -> UNSTABLE (max Re = 1.156e+00) L4: eigenvalues [+0.6321-0.9484i, +0.6321+0.9484i, -0.6321-0.9484i, -0.6321+0.9484i] -> UNSTABLE (max Re = 6.321e-01) L5: eigenvalues [+0.6321-0.9484i, +0.6321+0.9484i, -0.6321-0.9484i, -0.6321+0.9484i] -> UNSTABLE (max Re = 6.321e-01) L1, L2 and L3 are unstable saddles at every mass ratio tried: they always carry one real positive and one real negative eigenvalue (plus a purely imaginary pair), because Oxx>0 and Oyy<0 there, so Oxx*Oyy<0 and the quartic below always has a positive root for lambda^2. That has nothing to do with mu -- collinear points are unstable no matter how the mass is split between the primaries. ============================================================================== SECTION 4. Critical mass ratio for L4/L5 stability ============================================================================== Classical Gascheau/Routh value: mu_crit = (1 - sqrt(23/27)) / 2 = 0.0385208965 equivalently a primary:secondary mass ratio of 24.959936 : 1 (commonly quoted as ~24.96 : 1) bisection on the discriminant 1-27mu(1-mu): mu_crit = 0.0385208965 bisection on the numerical eigenvalues (no quartic assumed): mu_crit = 0.0385208525 club (eigenvalue bisection) minus classical value: -4.396e-08 (1.14e-04% relative) This is a deterministic root-find, not a statistical estimate, so the right comparison is bisection resolution, not a z-score. The discriminant version (exact float64 arithmetic, no eigendecomposition) matches the classical value to 0 at printed precision. The eigenvalue version disagrees by about 4.4e-08, close to sqrt(machine epsilon) (~1.5e-8): exactly at mu_crit the matrix has a repeated eigenvalue, and the eigenvalue problem is provably ill-conditioned there (a perturbation in an eigenvalue near a degeneracy scales as the SQUARE ROOT of the perturbation in the matrix, not linearly). That is the expected size of this disagreement, not a bug, and the method still confirms the classical result to five significant figures using no L4-specific algebra at all, only the general linearisation matrix and numpy's eigenvalue solver. ============================================================================== SECTION 4b. Where the real Solar System pairs sit relative to mu_crit ============================================================================== Sun-Earth mu/mu_crit = 0.00008 BELOW mu_crit: L4/L5 linearly stable Sun-Jupiter mu/mu_crit = 0.02475 BELOW mu_crit: L4/L5 linearly stable Earth-Moon mu/mu_crit = 0.31543 BELOW mu_crit: L4/L5 linearly stable Pluto-Charon mu/mu_crit = 2.82280 ABOVE mu_crit: L4/L5 linearly UNSTABLE equal-mass pair (illustrative, not real) mu/mu_crit = 12.97997 ABOVE mu_crit: L4/L5 linearly UNSTABLE ============================================================================== SECTION 4c. Collinear-point series error grows with mu (as expected of a small-mu expansion) ============================================================================== mu L1 rel. err % L2 rel. err % L3 rel. err % 1.000000e-06 0.0000 0.0000 0.0000 1.000000e-05 0.0002 0.0002 0.0000 1.000000e-04 0.0024 0.0021 0.0000 1.000000e-03 0.0243 0.0205 0.0000 1.000000e-02 0.2464 0.2023 0.0013 3.000000e-02 0.7467 0.5999 0.0120 3.852090e-02 0.9615 0.7676 0.0198 1.000000e-01 2.5308 1.9597 0.1388 Relative error here is measured against each point's own offset from the nearer reference (mass or barycentre), which is the natural small parameter the series expands in. It grows monotonically with mu across four decades, confirming this is a genuine asymptotic series with a finite radius of usefulness, not an accidental match at one mass ratio. ============================================================================== SECTION 5. Test-particle integration near each point (dt=0.002, T=40.0, eps=1.0e-05, threshold=1.0e-02) ============================================================================== Same generic perturbation direction (normalised (1,1,0.5,-0.5)) and the same magnitude used at every point, so the comparison across points is not an artefact of choosing a favourable direction for one point and an unfavourable one for another. Earth-Moon (mu = 1.215058e-02, mu/mu_crit = 0.3154) L1: escaped past 1e-02 at t=2.69 (max dist reached 2.70e-01) Jacobi drift |dC/C0| max = 3.90e-15 L2: escaped past 1e-02 at t=3.60 (max dist reached 6.07e+00) Jacobi drift |dC/C0| max = 2.80e-16 L3: escaped past 1e-02 at t=27.34 (max dist reached 9.53e-02) Jacobi drift |dC/C0| max = 2.95e-16 L4: contained for the full run (max dist 1.99e-04) Jacobi drift |dC/C0| max = 2.97e-16 L5: contained for the full run (max dist 5.76e-05) Jacobi drift |dC/C0| max = 2.97e-16 Pluto-Charon (mu = 1.087368e-01, mu/mu_crit = 2.8228) L1: escaped past 1e-02 at t=2.33 (max dist reached 5.27e-01) Jacobi drift |dC/C0| max = 3.43e-15 L2: escaped past 1e-02 at t=4.26 (max dist reached 1.87e+01) Jacobi drift |dC/C0| max = 2.55e-16 L3: escaped past 1e-02 at t=11.74 (max dist reached 5.48e+02) Jacobi drift |dC/C0| max = 2.86e-16 L4: escaped past 1e-02 at t=15.01 (max dist reached 2.94e+00) Jacobi drift |dC/C0| max = 3.06e-16 L5: escaped past 1e-02 at t=20.76 (max dist reached 1.46e+00) Jacobi drift |dC/C0| max = 3.06e-16 Every collinear point escapes in both systems: mu decides nothing about L1/L2/L3, exactly as the eigenvalues in Section 3 said. L4 and L5 stay contained for the whole run in Earth-Moon (mu below mu_crit) and escape in Pluto-Charon (mu above mu_crit) -- the same perturbation, the same distance, the same code, a different outcome purely because of which side of 0.038521 the mass ratio falls on. The Jacobi-constant drift column is the integration's own honesty check: every value is many orders of magnitude smaller than the distances being measured, so the escapes above are real dynamics, not integration error masquerading as one. ============================================================================== SECTION 6. Monte Carlo escape probability at L4 (N=1500 trials per mass ratio, T=70, eps=1e-04, threshold=1e-02) ============================================================================== mass ratio mu mu/mu_crit escaped/N p_hat SE max |dC/C0| Earth-Moon 1.215058e-02 0.315 0/1500 0.0000 0.0000 5.94e-16 0.50 x mu_crit 1.926045e-02 0.500 0/1500 0.0000 0.0000 5.96e-16 0.85 x mu_crit 3.274276e-02 0.850 0/1500 0.0000 0.0000 5.98e-16 0.95 x mu_crit 3.659485e-02 0.950 0/1500 0.0000 0.0000 5.99e-16 0.98 x mu_crit 3.775048e-02 0.980 282/1500 0.1880 0.0101 4.50e-16 1.00 x mu_crit 3.852090e-02 1.000 1200/1500 0.8000 0.0103 6.00e-16 1.02 x mu_crit 3.929131e-02 1.020 1477/1500 0.9847 0.0032 6.00e-16 1.05 x mu_crit 4.044694e-02 1.050 1499/1500 0.9993 0.0007 6.00e-16 1.20 x mu_crit 4.622508e-02 1.200 1500/1500 1.0000 0.0000 6.01e-16 Pluto-Charon 1.087368e-01 2.823 1500/1500 1.0000 0.0000 4.59e-16 Escape probability rises from 0 deep inside the stable region to 1 well past mu_crit, with the transition centred close to mu_crit itself, purely from re-running the identical experiment with a different mass ratio: nothing about the perturbations, the threshold or the integration changed. The Jacobi-constant column confirms the escapes are dynamics rather than integration error, tracked only up to first threshold crossing per trial, exactly as in Section 5. ============================================================================== SECTION 6b. Convergence of the escape-probability estimate as trials accumulate ============================================================================== reusing the 1500-trial run already made for 0.98 x mu_crit in Section 6 above, in the order the trials were actually drawn, rather than drawing a fresh batch. case shown: 0.98 x mu_crit (mu = 0.037750) final estimate after 1500 trials: p_hat = 0.1880 +/- 0.0101 (1 SE) running estimate at selected trial counts: n=10 p_hat=0.3000 SE=0.1449 n=30 p_hat=0.2667 SE=0.0807 n=100 p_hat=0.2300 SE=0.0421 n=300 p_hat=0.2167 SE=0.0238 n=1000 p_hat=0.1860 SE=0.0123 The running estimate visibly stops moving by a few hundred trials and the SE band it carries shrinks as 1/sqrt(n), exactly as a binomial proportion should; Figure 5 in the article plots this trace directly, which is the convergence check this study's Monte Carlo component is required to show, not merely assert. ============================================================================== Done ============================================================================== total wall-clock time: 132.3 s