Climate Formula

Earth’s climate-relevant orbital forcing arises from the gravitational interplay of all eight planets. Their orbital and rotational cycles synchronise over a common Solar System Resonance Cycle of 8H = 2,682,536 years, and every climate-relevant cycle on Earth therefore lands at an integer divisor of 8H. The canonical climate formula is a three-layer decomposition: L1 orbital forcing on the 8H integer-divisor lattice (32 components) + L2 climate-internal periodic response (the 405-kyr carbon-cycle thermostat and its harmonics) + L3 Heaviside step transitions at major Cenozoic boundary-condition changes (PETM, EOT, Mi-1, MMCT, iNHG, MPT). Sequential ridge regression fits the formula independently per climate regime, validated on four independent proxy records: LR04 benthic δ¹⁸O (Lisiecki & Raymo 2005), CENOGRID δ¹⁸O + δ¹³C (Westerhold 2020, 67 Myr Cenozoic), EPICA Dome C CO₂ (Bereiter 2015, 800 kyr ice-core composite), and CenCO2PIP atmospheric CO₂ (Consortium 2023, 66 Myr Bayesian multi-proxy synthesis).

This page is the architecture reference for the climate formula. Companion pages in the Climate section:

• Climate Summary — the higher-level synthesis statement this formula supports
• L1 Attribution — per-integer dual attribution (0 of 32 fully agree with Berger)
• Insolation Null Test — empirical demonstration that adding Berger insolation features yields ΔR² ≈ 0
• Related Work — position against peer-reviewed literature

Computational reproducibility scripts are in the 3d repository .

Canonical 3-layer Climate Formula vs LR04 δ¹⁸O on the post-MPT window (0–700 kyr BP). Red curve = formula (L1 lattice + L2 carbon thermostat + L3 step transitions, fit by sequential ridge regression); black = observed LR04 stack. R² = 0.87 within this regime.

Headline findings

1. R² up to 0.87 post-MPT LR04 — canonical 3-layer formula reaches R² = 0.87 (post-MPT LR04), 0.85 (EPICA CO₂), 0.76 (CenCO2PIP 0–66 Ma). Full per-regime table in §The canonical three-layer climate formula; cross-record synthesis in Climate Summary.
2. Every significant LR04 peak sits at an integer divisor of 8H. L1 has 32 integers — 25 canonical (Berger / Mars–Jupiter beats) + 6 precession-band sidebands + 1 Berger-quintet completion (n=141). Every one ALSO has an Earth-planet PLANET_CYCLES beat; 0 of 32 fully agree with Berger on both planet AND mechanism (26 name a different planet, 1 same planet/different mechanism, 5 have no Berger prediction). Full per-integer mapping in L1 Attribution.
3. Mars and Jupiter co-dominate the per-planet fingerprint. Jupiter has 4 direct L1 matches (n=16, 21, 39, 65); Mars has 3 (n=16, 21, 68). Mars’s exclusive match is n=68 (Mars.ICRF); Jupiter’s is n=39 (Jupiter.Peri_ecl). The Mars-Jupiter resonance lock at 8H/36 (Mars.Peri_ecl = Jupiter.AscNode = 74.5 kyr) is the cleanest planet-pair resonance in the framework. Mars also participates in n=25 (Mercury-Mars s₁−s₄ nodal, 100-kyr centroid), n=28 (Mars-Jupiter eccentricity, Berger 95k), and n=53 (Mars-Uranus s-beat). Detail in §Per-Planet Contributions.
4. The 100-kyr glacial cycle is an inclination-side eigenmode beat, not direct eccentricity forcing. Empirical centroid at the Mercury-Mars s₁−s₄ nodal beat at 107.3 kyr. Detail in §The 100-kyr Glacial Cycle.
5. Pre-MPT and post-MPT differ in climate sensitivity, not orbital forcing. Three per-regime fits: pre-iNHG R² = 0.43, iNHG-MPT 0.72, post-MPT 0.87. Detail in §Pre-MPT vs Post-MPT.
6. Forward projection: next natural glaciation peak ~60,500 AD (58,500 years from now); strongest in the next quarter-million years at ~198,500 AD; warmest interglacial at ~166,000 AD. Detail in §Forward Projection.

The canonical three-layer climate formula

Architecture

$C(t) \;=\; c_0 \;+\; \underbrace{\sum_{n \in N_{L1}} \left[ a_n \cos\!\left(\tfrac{2\pi n t}{8H}\right) + b_n \sin\!\left(\tfrac{2\pi n t}{8H}\right) \right]}_{\text{L1: orbital forcing}} \;+\; \underbrace{\sum_{p \in P_{L2}} \left[ \alpha_p \cos\!\left(\tfrac{2\pi t}{p}\right) + \beta_p \sin\!\left(\tfrac{2\pi t}{p}\right) \right]}_{\text{L2: carbon thermostat}} \;+\; \underbrace{\sum_{i \in T_{L3}} \gamma_i \cdot H(t - t_i)}_{\text{L3: boundary-condition steps}}$

where:

• L1 = 32 integers N_L1 = {9, 12, 14, 16, 18, 20, 21, 22, 25, 28, 30, 31, 35, 38, 39, 48, 50, 53, 65, 66, 68, 73, 76, 96, 107, 110, 113, 120, 134, 141, 152, 185} (25 canonical Berger / Mars–Jupiter integers plus 6 precession-band sidebands surfaced by the all-integer MTM F-test scan plus n=141 Berger-quintet completion added 2026-05-28; sidebands in bold)
• L2 = 3 periods P_L2 = {404.5 kyr (silicate-weathering thermostat fundamental), 202.25 kyr (2nd harmonic), 134.83 kyr (3rd harmonic)} — confirmed carbon-amplified across LR04, CENOGRID, and EPICA
• L3 = up to 6 transitions T_L3 = {PETM (56 Ma), EOT (34 Ma), Mi-1 (23 Ma), MMCT (14 Ma), iNHG (2.7 Ma), MPT (1.0 Ma)}, applied only when inside the fit window; H(·) is the Heaviside step function
• t = age in kyr BP (positive = past, negative = future, t = 0 ≈ 2000 AD)
• 8H = 2,682,536 years (Solar System Resonance Cycle)
• C(t) = normalised δ¹⁸O / δ¹³C / CO₂ proxy depending on the fitted record

Sequential ridge regression, fit per regime

Coefficients are fit by sequential ridge regression: L1 first (with ridge regularization λ = 1 — see §9.5 of doc 92 for derivation ), then L2 on the residual, then L3 on the L1+L2 residual. Sequential fitting resolves the L1↔L2 collinearity at short windows (e.g. L2 135-kyr 3rd harmonic sits within one Rayleigh element of L1 lattice integer 8H/20 = 134 kyr at post-MPT 1000-kyr window); ridge resolves the L1↔L1 collinearity within the lattice itself (post-MPT condition number ~632, max VIF ~7×10⁴ under plain OLS).

The formula is fit independently per climate regime — the same lattice positions hold across regimes, but per-line amplitudes/phases vary. Per-regime fit results:

Interpretation: at Quaternary scales (post-MPT LR04, EPICA CO₂) L1 carries the heavy lifting (~83–87% R² from the lattice alone). At deep-time scales (CENOGRID, CenCO2PIP) L3 step transitions dominate (~60% of the variance over 67 Myr) — the lattice oscillations average toward zero across many cycles, but the Cenozoic story is the discrete step-changes at major transitions.

The same 8H lattice fitted to atmospheric CO₂ (EPICA Dome C; Bereiter 2015) — an independent proxy with independent chronology and entirely different recording physics from benthic δ¹⁸O. R² = 0.83 from L1 alone, 0.85 with L2 added: the orbital structure that drives ice volume also drives atmospheric carbon-cycle dynamics. The per-line carbon-amplification ratios (EPICA amp / LR04 amp) identify which lattice members manifest primarily through carbon-cycle response — see §10.3 of the doc 92 climate-formula architecture .

What survived investigation and what didn’t

The formula commits to 32 + 3 + 6 = 41 structural components. Other candidates were investigated and excluded:

• 13H = 4.36 Myr Boulila libration — investigated as a candidate long-period L2 line; failed cross-window stability (amplitude CV 42–50%, δ¹⁸O phase circular std 97.9° ≈ uniform random across CENOGRID windows). Not deployed.
• 9-Myr long-period carbon resonance — investigated; ΔR² CENOGRID δ¹³C = +0.078 with strong carbon-amplification ratio 4.2. Promoted to “investigated L2 candidate” only — cross-window stability test pending. Not yet deployed.
• 405-kyr 4th and 5th harmonics (101 and 81 kyr) — collinear with L1 lattice integers 8H/27 and 8H/33 respectively; absorbed by L1.
• 23 off-lattice Laskar eigenmode beats — all fail strict promotion criteria (ΔR² > 0.005 + carbon-amplification ratio > 0.5).
• Piecewise / polynomial / LOESS detrending — superseded by L3 step covariates (5–10× better R² on CENOGRID).
• Linear-response ODE models for L2 — failed (all L2 ↔ L1 amplitude correlations |r| < 0.4); L2 dynamics are nonlinear (threshold / hysteresis / saturating).
• L2-CO₂-feedback as separate components — captured implicitly: EPICA reveals strong obliquity-band CO₂ amplification (8H/66 ratio 15.79), but no separate L2 components are added; the obliquity-band L1 integers carry the signal and the per-line carbon-amplification ratio in the fit output is the diagnostic.

Full inclusion / exclusion inventory with sources: doc 92 §9 .

The 32 L1 lattice integers

Each n corresponds to a specific celestial-mechanics quantity, of two distinct kinds. (i) Eigenmode beats — Laskar 2004 secular eigenfrequency combinations (g_j apsidal, s_j nodal; k = Earth’s general precession in longitude, 50.29″/yr); these are composite, multi-planet modes, and any planet labels on g_j / s_j follow Berger’s convention (the model does not adopt the single-planet attribution — see Eigenfrequencies §“The eigenmodes are real…”). (ii) Direct planet cycles — rows that also read ”= 8H/N (model direct)” name the model’s own planet-specific period from the Solar System Resonance Cycle period table (e.g. n=39 = Jupiter’s ecliptic perihelion, n=65 = Jupiter ICRF = Saturn ecliptic); at those integers the model’s direct cycle and the eigenmode beat coincide. Amplitudes shown are illustrative LR04 OLS values for the 25 original integers; per-regime canonical fits give different amplitudes (see §“Sequential ridge regression, fit per regime” above).

Several integers correspond directly to specific planet cycles from the Solar System Resonance Cycle period table: n=9 (Mercury Axial = AscNode, Cassini lock), n=16 (Mars Axial / Jupiter Obliquity), n=21 (Mars Obliquity / Jupiter Axial — “Mars–Jupiter Axial-Obliquity Swap”), n=39 (Jupiter ecliptic perihelion), n=65 (Saturn ecliptic = Jupiter ICRF perihelion — the Law 6 resonance lock), n=68 (Mars ICRF perihelion), n=110 (Venus ICRF = Venus Obliquity). The LR04 peaks at n=12, n=35 and n=53 are generated by multi-planet beats rather than direct period matches. The n=66 integer is the obliquity-band arithmetic-mean cycle length — a cycle-counting artifact at full LR04 (Jensen’s inequality), not a direct planet-cycle match. The 6 italicised sideband rows (n=96, 107, 110, 134, 152, 185) are the precession-band sidebands added to the canonical lattice from the all-integer MTM F-test scan (each significant at α=0.05 individually but collinear with the 25 canonical lines at full-record resolution — gain is small on full LR04 but accumulates in regime windows). The bold-italic seventh row (n=141, k+g₃ Earth climatic precession at ~19 kyr) completes the Berger precession quintet — subthreshold in LR04 but 3σ-significant in the Cheng full Asian-monsoon speleothem record (3.60× median); it closes the Wigley 1976 / Berger 1978 combination tone 1/95 ≈ 1/19 − 1/23.7 by beating with n=113 to produce the 95-kyr eccentricity peak at n=28. All 32 land on integer divisors of 8H because 8H is the natural synchronisation period of the solar system.

In-app Explorer

Closure: no orphan peaks off the 8H lattice

The strongest possible test of the 8H framework is its closure: do all significant spectral peaks in LR04 land on integer divisors of 8H, or are there orphan peaks at non-integer positions that would imply forcing from outside the planetary-eigenmode framework? A peak at, say, n = 43.5 with non-negligible amplitude would falsify the framework.

We performed this closure test by fitting all 200 integer divisors of 8H jointly to LR04 (joint OLS, R² = 0.443) and scanning the residual for non-integer power. The result:

• Residual amplitude at integer positions: 0.000 (machine zero, by construction)
• Noise floor at random non-integer positions: 0.029 median, 0.127 (95th percentile)
• 14 residual peaks above noise threshold sitting >0.3 from any integer — every single one between two adjacent integer divisors in or near the formula

The biggest “orphan” (n = 65.45, period 40.99 kyr, amp 0.544) sits between integers 65 (k+s₃ Earth obliquity, 41.27 kyr) and 66 — the cycle-counting arithmetic mean position of the non-stationary obliquity band (a Jensen’s-inequality phenomenon: when the obliquity band’s cycle length varies, the spectrum’s amplitude-weighted mean of 1/T does not equal 1/⟨T⟩). The other 13 orphans similarly sit between adjacent integer beats — they are the expected fingerprint of cycle-length dispersion (LR04 obliquity cycles vary 31–59 kyr) and spectral leakage between close integer signals.

Crucially: zero orphan peaks land in “empty” regions of the 8H lattice — no peaks at positions like n = 12.7, n = 42.3, n = 80, or n = 140 that would suggest genuinely off-lattice forcing. The 8H integer-divisor framework captures the full frequency structure of LR04’s significant spectral content.

This closure result is one of several independent empirical confirmations of the framework, alongside the 405-kyr-absence test, the bispectral coupling absence (see Eccentricity attribution has specific empirical headwinds below), and the broader test suite documented in Framework validation and the 405-kyr investigation below.

Framework validation and the 405-kyr investigation

Beyond the closure test above, the framework has been subjected to a battery of independent falsifiable tests. The combined picture: most predictions hold strongly; one specific signal — the 405-kyr long-eccentricity cycle — does not fit the orbital lattice and is best explained as a climate-system internal phenomenon rather than a planetary orbital beat.

Super-cycle hypothesis on deep-time geology: NULL

The Plio-Pleistocene epochs each align closely with 1×8H = 2.68 Myr, raising the question: does this generalize? A pre-registered test asked whether 20 major Phanerozoic geological/climatic events (mass extinctions, period boundaries, Cenozoic climate transitions) cluster preferentially at integer multiples of 8H (or of H, 8× tighter).

Result: NULL. The 8H super-cycle hypothesis is rejected at deep-time scales. The Plio-Pleistocene 1×8H + 1×8H pattern is most parsimoniously a climate-response amplification artifact combined with statistical coincidence, not a deep-time geological clock. The framework describes orbital geometry that paces possible climate transitions, not a structural clock that causes geological events.

Fourteen follow-up tests on framework predictions

A battery of fourteen independent falsifiable tests (A–N) covering chronology validation, statistical baselines, bispectral coupling, external-publication convergence, cross-validation, phase prediction, cross-proxy replication, 67-Myr Cenozoic generalization, per-line Thomson MTM F-test significance, time-frequency centroid stability, and head-to-head against standard Laskar 2004 placements. Combined result: 16 clean positives, 2 partials, 5 nulls — each null traces cleanly to one of four mechanistic causes (Rayleigh resolution at T < 8H, MPT amplitude non-stationarity, proxy-specific spectral concentration, non-linear climate lags); none falsifies the framework. Full per-test reproduction in doc 91 (3d repo) .

The 405-kyr investigation

The empirical 405-kyr line is the prominent long-eccentricity cycle in pre-Pleistocene climate records. Standard Milankovitch attributes it to the Laskar g₂−g₅ “Venus-Jupiter apsidal beat”. In this framework’s planet motions, Venus and Jupiter have different apsidal periods than Laskar’s eigenfrequencies (Venus: 8H/6 = 447,089 yr retrograde; Jupiter: 8H/39 = 68,783 yr; their actual beat is ~58–79 kyr, not 405 kyr). So the standard attribution does not carry over. The 405-kyr investigation asks: what is it, then?

Empirical measurement. CENOGRID full Cenozoic (67 Myr span) places the line at 404.5 kyr as a narrow peak (FWHM = Rayleigh resolution — a clean single-frequency line). Position is stable across Eocene, Oligocene, Paleocene intervals to within 0.7%.

Mathematical proof: 405 kyr is unreachable on the 8H lattice. A systematic search across all pair and triplet beats among the 46 cycles in the 8H period table finds zero combinations within ±3% of 405 kyr. The closest matches cluster at 8H/7 = 383 kyr (5.4% below) and 8H/6 = 447 kyr (10.4% above). Mathematically: every beat between period-table cycles is itself a 8H/integer fraction; 405 kyr sits in the gap between adjacent integers, and no combination can reach it.

Carbon-cycle amplification in CENOGRID δ¹³C. The δ¹³C/δ¹⁸O amplitude ratio at 405 kyr is 1.53× in the full Cenozoic — compared to 0.49 at obliquity and 0.66 at precession. The 405-kyr signal lives preferentially in the carbon record — the signature predicted by the silicate-weathering thermostat resonance at ~400 kyr (Walker, Hays & Kasting 1981; Pälike et al. 2006 “heartbeat of the Oligocene climate system”).

Phase analysis: entrained internal oscillator. δ¹³C and δ¹⁸O phase residuals from a linear orbital model correlate at only r = 0.21. About 96% of phase variance evolves independently between the two proxies. The 405-kyr cycle behaves as an entrained internal oscillator — the carbon cycle has its own ~400-kyr resonance, with long-period orbital eccentricity supplying energy that synchronises but does not rigidly drive it.

Sub-harmonic (“17 × precession”) alternatives are ruled out by three independent tests, and the g₄−g₃ 2.4-Myr beat shows no comparable carbon-amplification (ratio 0.20 vs 1.53 for 405-kyr) — refining the interpretation: the silicate-weathering thermostat has a narrow resonance peak near 400 kyr, not a broad long-period amplification. Full investigation detail in doc 91 (3d repo) .

Architectural conclusion: orbital geometry + climate physics

The combined findings produce the three-layer picture committed to in the canonical formula above:

Layer 1 — Orbital motions on the 8H integer-divisor lattice. Axial precessions, perihelion advances, ascending-node regressions, obliquity oscillations, eccentricity-cycle wobble periods. All of the form 8H/N for integer N. The 32 framework integers carry the bulk of variance at Quaternary scales (see the per-regime ridge fit table above for full R² breakdown). Beats between any period-table cycles also land on the 8H lattice (mathematical closure).

Layer 2 — Climate-internal periodic responses. Carbon-cycle silicate-weathering thermostat (τ ≈ 400 kyr) — resonantly produces the empirical 405-kyr line, plus its 2nd (202 kyr) and 3rd (135 kyr) harmonics, in δ¹³C records. These three lines are deployed in the canonical formula’s L2 layer. Ice-sheet hysteresis (post-MPT) is real climate physics but not modelled as a periodic L2 component; its signature shows in per-line carbon-amplification ratios of the L1 lattice integers (e.g. 8H/22 ratio 12.84 on CENOGRID, 8H/66 ratio 15.79 on EPICA CO₂).

Layer 3 — Discrete boundary-condition transitions. PETM (56 Ma), EOT (34 Ma), Mi-1 (23 Ma), MMCT (14 Ma), iNHG (2.7 Ma), MPT (1.0 Ma) — six Heaviside step terms. L3 carries the dominant variance at deep-time scales (~60% of CENOGRID δ¹⁸O over 67 Myr).

The 8H integer-divisor structure is a Layer-1 statement about orbital motions. The 405-kyr cycle is not a Layer-1 phenomenon — no combination of period-table cycles produces it. It is a Layer-2 climate-system phenomenon: the silicate-weathering thermostat resonates at ~400 kyr (a property of Earth’s carbon-cycle physics, not orbital geometry), entrained by long-period orbital eccentricity forcing — without any specific Venus-Jupiter or other planetary-pair attribution required.

Cycles observed in climate records that don’t coincide with 8H/n positions and aren’t already in the L2 thermostat family — the 13H and 9-Myr candidates being the documented cases — remain investigated but not yet deployed.

Per-Planet Contributions

For each climate peak, cross-referencing against the full period table (8 planets × 6 cycle types) reveals which planets contribute directly versus via eigenmode beats. The per-integer dual attribution (Berger vs Holistic top-1 beat) is fully tabulated in L1 Attribution; this section focuses on aggregate per-planet patterns.

Direct L1 matches per planet

LR04 full window. A “direct L1 match” = the planet cycle exactly equals an L1 lattice integer:

Mars and Jupiter co-dominate. Both have multiple direct L1 matches; Jupiter has 4, Mars has 3. Mercury and Venus each have 2 cycles converging on a single L1 integer (Mercury’s Cassini lock at n=9, Venus’s ICRF=Obliq coincidence at n=110). Earth’s six cycles all sit off the L1 lattice — a structural property of the framework.

The Mars-Jupiter resonance lock

A second structural finding: Mars.Peri_ecl = Jupiter.AscNode = 8H/36 (= 74.5 kyr). Two distinct cycles on two different planets converge on the same period. The L1 lattice has n=35 nearby (76.6 kyr, generated by a 3-term Earth.Axial(104) − Mercury.ICRF(93) + Saturn.Obliq(24) beat), one Rayleigh element away from the 8H/36 resonance.

Mars also contributes to multiple climate peaks via eigenmode beats:

• n=25 (107.3 kyr, the 100-kyr-band centroid) — the s₁−s₄ Mercury-Mars nodal beat
• n=28 (95.8 kyr, the classical “Berger 95k”) — the g₄−g₅ Mars-Jupiter eccentricity beat
• n=53 (50.6 kyr, near Mars.Ecc = 8H/52) — the Mars.AscNode(64) − Uranus.AscNode(11) s-beat

Neptune and Uranus appear only pre-MPT

The pre-MPT window (1,200–3,000 kyr BP) reveals three additional peaks not in LR04 full that involve outer planets:

• n = 14 (191.6 kyr) — g₂−g₈ Venus-Neptune eccentricity beat
• n = 30 (89.4 kyr) — g₃−g₇ Earth-Uranus eccentricity beat
• n = 38 (70.6 kyr) — s₈−s₃ Neptune-Earth nodal beat

Neptune’s eigenmodes are the slowest (g₈ ≈ 0.67 ″/yr, |s₈| ≈ 0.69 ″/yr); beats with inner planets produce relatively long periods (~70–200 kyr) with small amplitudes. Pre-MPT, ice sheets were smaller and less hysteretic — slow outer-planet beats registered above the noise floor. Post-MPT large ice sheets impose a ~100-kyr hysteretic response that dominates and filters out the slower outer-planet signals.

Cross-planet obliquity validation

The model’s obliquity-period predictions for the inner solar system match three independent published references to within 2.5 % with zero free parameters:

Mercury’s 0.09 % match against the independent Cassini-state forced-obliquity calculation (Bills & Comstock used Cassini-state dynamical theory, corroborated by Yseboodt & Margot 2006, Peale 2005, Bois & Rambaux 2007) is the model’s tightest cross-validation against non-Holistic published references.

The 100-kyr Glacial Cycle

The 100-kyr glacial cycle is paleoclimatology’s most-debated feature. The data point clearly to an inclination-side / orbital-plane origin via planet-pair eigenmode beats rather than direct eccentricity forcing.

What the spectrum shows

LR04’s 100-kyr band is a broad single peak spanning ~80–125 kyr — not the split structure (95k + 125k) that direct eccentricity forcing would produce. The energy-weighted centroid sits at n = 25 (= 107.3 kyr in 8H units), with strong adjacent contributions at n = 28 (95.8 kyr) and weaker contributions at n = 22 (~122 kyr).

Empirical centroid: Mercury-Mars nodal beat

The dominant 100-kyr-band integer n = 25 corresponds to the s₁ − s₄ Mercury-Mars nodal beat in eigenmode terms (predicted 25.11, error 0.11). This is a nodal (orbital-plane) beat between two inner rocky planets, not an eccentricity beat. Consistent in spirit with Muller & MacDonald 1997’s “inclination, not eccentricity” framing, with a specific eigenmode identification not previously stated.

Eccentricity attribution has specific empirical headwinds

Direct-eccentricity attribution faces three discriminating empirical failure modes that the inclination-side framework does not:

Theoretical proposal: H/3 second-obliquity component

Within the inclination-side family, the model proposes Earth’s own inclination precession (H/3 = ~111,772 years = 111.77 kyr) produces a second obliquity component, contributing alongside the Mercury-Mars planet-pair nodal beat. The mechanism (no dust required):

Inclination precession (H/3) → second obliquity component at H/3 → standard Milankovitch insolation forcing → ice sheets

Every step after “second obliquity component” is standard Milankovitch physics; the mechanism needs no new forcing.

At T = 1.2 Myr the Rayleigh resolution at P = 110 kyr is ΔP ≈ 10 kyr, so 95k / 100k / 112k are spectrally collinear. H/3 lies within one Rayleigh element of the empirical 107-kyr centroid and remains a candidate within the inclination-side family — but cannot be singled out from the Mercury-Mars nodal beat at this data length. The empirical signal is consistent with this family; the specific H/3 second-obliquity contribution is theoretical.

Earth’s Fibonacci H/n positions are empty in the climate spectrum

The model’s Earth precession periods (Fibonacci divisors of H — H/3, H/5, H/8, H/13, H/16) all sit at near-zero amplitude in LR04 (0.014 to 0.078, none above 0.085). This is structurally consistent with eigenfrequencies §“The eigenmodes are real…”: the Fibonacci H/n periods are Earth’s geometric precession in its own orbital frame; climate observes Earth in the heliocentric frame where planetary-induced orbital-plane motion adds in, shifting the observed signal to integers adjacent to the Fibonacci anchors (Berger’s k+g_j and k+s_j sub-peaks structure).

Climate couples to ecliptic-frame periods, not ICRF-frame periods

A complementary structural observation explains which per-planet periods appear in the climate record. Earth’s climate responds to orbital forcing in the ecliptic frame — the plane of Earth’s own orbit — not the inertial ICRF (fixed to distant quasars). When another planet’s gravitational coupling registers in the paleoclimate spectrum, it appears at that planet’s ecliptic-frame perihelion period, not its ICRF-frame period:

• Mars’s ecliptic perihelion (8H/36 = 74.5 kyr) sits close to the LR04 n=35 peak at 76.6 kyr (permutation-null z ≈ +4.4 on the post-MPT window); the n=35 peak is generated by a 3-term beat (Earth.Axial − Mercury.ICRF + Saturn.Obliq), with Mars near the position. The 8H/36 period exactly matches Jupiter’s ascending-node period — a Mars-Jupiter resonance lock.
• Saturn’s ICRF perihelion (8H/169 = 15.9 kyr) — a pure inertial-frame period with no ecliptic-frame counterpart — shows no detectable LR04 signal (z ≈ −1, below the noise floor).
• The apparent exceptions where an ICRF period does appear (Earth ICRF = H/3 inclination precession; Jupiter ICRF = 8H/65, the period at which it drives the obliquity beat) are exactly the cases where the ICRF period coincides with an independently climate-active ecliptic-frame feature (Earth’s inclination cycle; the k+s₃ obliquity beat, which is also Saturn’s ecliptic perihelion).

This is physically expected: the ICRF has no dynamical coupling to the solar system’s internal secular motions, while Earth’s climate is set entirely by ecliptic-frame geometry — eccentricity, the direction of perihelion, and the tilt of the orbital plane. It is also why the dynamical (ecliptic-frame) attribution of each planet’s perihelion (Jupiter 8H/39, Saturn −8H/65; see Fibonacci Laws — Law 6) is the climate-relevant one, rather than the kinematic Fibonacci anchor.

Pre-MPT vs Post-MPT

Canonical formula on full LR04 with three independent per-regime fits stitched together (pre-iNHG / iNHG-MPT / post-MPT). Stitched R² = 0.93. The amplitude growth from left to right is the MPT story: orbital forcing is essentially stationary across 5.3 Myr, but ice-sheet hysteresis amplification grew dramatically after iNHG (~2.7 Ma) and again after MPT (~1 Ma). The L3 step transitions encode this directly.

Same forcing, different response

Three structural transitions are visible in the full LR04 record:

The orbital forcing — planetary eigenmodes, Earth’s precession, obliquity oscillation — is essentially constant over the past 5 Myr. What changed is the climate system’s sensitivity: long-term CO₂ decline, growing Northern Hemisphere ice sheets, and ice-sheet hysteresis crossing the Willeit threshold (~1 Ma).

The L1 lattice positions are stationary across the full 5 Myr (orbital eigenmodes don’t move). Within a single regime the per-line amplitudes are also stationary; across regimes they change to reflect the changing climate sensitivity — pre-MPT amplitudes differ from post-MPT amplitudes, which is why the canonical formula is fit per regime (see the per-regime table above). The conflated single fit on full LR04 forces one amplitude set across both regimes and reaches only R² = 0.25, versus the 3-regime per-regime fits which reach much higher per-window R².

MPT amplitude growth pattern

Windowed amplitude analysis across the MPT (pre-MPT 1,500–2,500 kyr vs post-MPT 0–1,000 kyr) gives concrete post/pre growth ratios:

The 41-kyr peak decreased — consistent with the standard “ice-sheet saturation silences the obliquity pacemaker” framing (Willeit 2019). The 100-kyr and 23-kyr bands both grew, with the 23-kyr growth concentrated in one specific sub-peak (k+g₅) rather than spread across the precession band.

Forward Projection

Forward projection within the post-MPT regime. Past (left of ‘today’) overlays the canonical formula on LR04 data; future (right of ‘today’) is the formula extrapolated forward 250 kyr. The next predicted natural glaciation peak sits at ~60,500 AD (58,500 years from now), with the warmest interglacial at ~166,000 AD and the strongest glaciation in the window at ~198,500 AD. Forward projection stays in scope only within the current regime (no scheduled boundary-condition shift in the next 250 kyr).

The next 250,000 years (canonical formula)

The canonical 3-layer formula is fit on the post-MPT regime alone (0–1 Ma) and extrapolated forward 250 kyr. Within this regime the fit reaches R² = 0.87 (LR04 δ¹⁸O) and R² = 0.85 (EPICA CO₂); forward extrapolation stays in-scope as long as the post-MPT regime continues. Ridge regression bounds the peak-to-peak forward range to ~4.7 ‰ normalized (plain OLS would extrapolate to ~24 ‰, which is unphysical).

Predicted glacial maxima:

Predicted interglacial peaks (between glacial maxima):

The 9-kyr LGM offset in past validation (predicted ~29 kyr BP vs. observed ~20 kyr BP) suggests the surface-temperature response can lag the orbital clock by several kyr; the actual surface peak-warmth and glacial-onset dates may differ from the orbital values above.

Comparison with other published forecasts

Several independent frameworks predict the natural length of the current interglacial. All converge on the qualitative answer that the next ~50–100 kyr will be an exceptionally long interglacial by Pleistocene standards, but differ in mechanism:

Where the frameworks agree. All five identify that we sit in an unusually weak eccentricity minimum (analogous to MIS 11 or MIS 19c) and predict a longer-than-typical interglacial. Without anthropogenic forcing, the consensus next-onset range is ~50–80 kyr.

Where the frameworks differ — the mechanism. Berger-Loutre, Loutre-Berger, and Tzedakis attribute the 100-kyr glacial cycle to direct eccentricity forcing modulated by precession; Ganopolski adds nonlinear CO₂-ice-sheet hysteresis as the dominant amplifier; the Holistic Climate Formula attributes the empirical 100-kyr-band centroid to the Mercury-Mars s₁−s₄ nodal beat at 107.3 kyr (a planet-pair orbital-plane coupling, see §100-kyr-cycle above), with eccentricity-attribution failing three discriminating empirical tests (405-kyr term absent, no bispectral 95k+125k coupling, wrong-family centroid).

Where the Holistic formula is distinctive. It commits to a fixed integer-divisor lattice (positions established from observed orbital motion, not fitted to climate), making it falsifiable in a way the standard frameworks are not — climate peaks at non-lattice frequencies would refute the L1 layer. The closure test confirms zero such peaks in LR04.

Validation against past glacials

The LGM offset of ~9 kyr is expected: pure orbital forcing predicts the insolation drive, not the ice-sheet response. Ice sheets carry memory; peak ice volume lags peak orbital cooling by several kyr (standard climate physics). The MIS 6 match within 2 kyr is excellent.

Pacing shift: the next 250 kyr departs from 100-kyr regime

The canonical forward projection shows intervals between predicted glacial peaks that are not ~100 kyr — the projection has two strong glaciations 138 kyr apart with weak intermediate wiggles, qualitatively different from the past 700 kyr of regular ~100-kyr pacing:

By contrast the past ~700 kyr in LR04 was dominated by ~100-kyr pacing (LGM → MIS 6 → MIS 8 → MIS 10 → MIS 11, mean interval ~100 kyr). The canonical projection shows two patterns simultaneously: (a) at the strong-event level, the two major glaciations sit 138 kyr apart — longer than the late-Pleistocene 100-kyr cycle; (b) at the local-maximum level, intermediate wiggles alternate between obliquity-band (43–48 kyr) and precession-band (22–25 kyr) intervals. Neither the regular ~100-kyr regime nor a clean ~41-kyr-paced “41-kyr world” describes the projection cleanly — it is a transitional/mixed pattern.

Why: The L1 layer is a sum of 32 sinusoidal components at orbital-eigenmode periods, each with stationary per-regime amplitude. The L1 layer is exactly 8H-periodic by construction, so the L1 orbital signal in the next 250 kyr is byte-identical to the L1 orbital signal in 2.43–2.68 Ma BC — the late-Pliocene “41-kyr world”. On top of that orbital backbone, the L2 405-kyr carbon thermostat and the post-MPT regime amplitudes modulate the response: in the next 250 kyr the obliquity-band integers (n=65, 66, 68 near 41 kyr) are partially active while the 100-kyr-band integers (n=22, 25, 28) are partially out of phase with each other; the L2 405-kyr family adds long-period modulation. The result is a mixed projection: the orbital component matches a 41-kyr-paced backbone but the carbon-cycle response adds longer-period structure. The same qualitative finding — a long current interglacial — was reached independently by Berger & Loutre 2002 .

Climate response is a separate question. Whether Earth’s surface climate actually tracks this mixed orbital pacing depends on ice-sheet hysteresis. Three scenarios: (A) ice sheets shrink below the Willeit 2019 threshold and the climate response re-couples to the orbital clock (a “return to 41-kyr world” at the L1 backbone level); (B) post-MPT hysteresis persists and the climate keeps producing ~100-kyr cycles even with weaker orbital triggers (“skipped” or subdued glacials, which matches the very-mild intermediate peaks in the projection); (C) anthropogenic CO₂ overrides orbital pacing entirely for the near-term (Ganopolski 2016). The 8H formula speaks to the orbital half only.

Empirical analogue: the late Pliocene (8H ago)

The canonical formula’s L1 layer is exactly periodic at 8H — every L1 component is an integer divisor of 8H, so C_L1(t) = C_L1(t + 8H) by construction (verified to floating-point precision). The L2 carbon-thermostat family (404.5 / 202.25 / 134.83 kyr) and L3 step terms are not 8H-periodic, but L1 carries the bulk of the orbital signal at sub-Myr scales. This means the orbital component of the forward projection in the next 250 kyr closely resembles the orbital signal in 2.43–2.68 Ma BC, making the late-Pliocene LR04 record a direct empirical analogue for how Earth’s climate responded to the same orbital forcing we’re about to enter.

Recomputed under the canonical 32-component multi-proxy ridge formula (per-regime fits: iNHG-MPT for the Pliocene window, post-MPT for the modern window):

(Reproduce with scripts/milankovitch_late_pliocene_analogue.py; output saved to data/milankovitch-late-pliocene-analogue.json.)

Two findings:

1. The orbital signal is empirically 41-kyr-paced in the late-Pliocene window — LR04 glacial-peak intervals cluster around the obliquity band (mean 39.5 kyr), and the dominant LR04 spectral period in the window sits at 41.3 kyr. By contrast the modern post-MPT window is 100-kyr-band-dominated (125 kyr). At the L1 orbital level, the next 250 kyr backbone is byte-identical to this 41-kyr-paced Pliocene signal — making the late Pliocene LR04 record a direct empirical anchor for what an L1-orbital-coupled response looks like.
2. The climate envelope reachable by orbital forcing alone is modest. Late-Pliocene LR04 raw detrended std is 0.257 ‰, vs 0.464 ‰ post-MPT — the post-MPT climate response is 1.8× more amplified than late Pliocene. Without large hysteretic ice sheets the late Pliocene climate tracked the orbital signal more directly and didn’t reach the extreme glacial-interglacial swings of the late Pleistocene (raw range −0.51 to +0.52 ‰ in late Pliocene vs −1.09 to +0.83 ‰ post-MPT). If ice sheets shrink below the Willeit threshold and climate re-couples to the L1 orbital clock, natural climate amplitude over the next 250 kyr would resemble the late Pliocene more than the late Pleistocene.

This empirical anchor strengthens Scenario A above: the late Pliocene happened, the LR04 record exists, and it shows what an orbital-coupled climate response to this exact L1 orbital backbone looks like — moderate-amplitude, 41-kyr-paced, without the extreme glacial swings of the post-MPT regime.

Limitations

Per the scope warning at the top of this page, the formula models orbital forcing (L1) plus the L2 silicate-weathering carbon-cycle thermostat plus L3 Cenozoic step transitions. It tells you when the orbital clock makes glaciation possible; the actual ice-volume response depends on the full climate system (ice-sheet hysteresis, CO₂ amplification beyond L2, internal variability, regional asymmetries, anthropogenic CO₂ — none modelled here).

Within-regime descriptive power is high; cross-regime predictive power is zero — formulas trained on pre-MPT data give catastrophically negative R² when applied to post-MPT, and vice versa. The post-MPT amplitude/phase structure is not derivable from pre-MPT; the lattice itself is regime-independent but the response amplitudes evolve with boundary conditions (CO₂, ice-sheet area). Forward projection therefore stays in scope only within the post-MPT regime (no scheduled boundary-condition shift in the next 250 kyr).

Anthropogenic CO₂ may delay the next natural glaciation by 50+ kyr in moderate-emission scenarios (Ganopolski et al. 2016 ) or up to ~100 kyr in high-emission scenarios. This is not modelled here. The Caillon et al. 2003 ice-core analysis argues CO₂ lags temperature as a feedback rather than driving it; if that view is correct, the orbital projection above largely holds regardless of CO₂ trajectory.

Methodology highlights

Reproducibility scripts are in the 3d repository’s scripts directory  (Python implementations of every spectral, beat-search, and bootstrap analysis described above).

Data sources

• LR04 stack — Lisiecki, L. E. & Raymo, M. E. (2005), Paleoceanography 20, PA1003. 57 globally distributed benthic δ¹⁸O records, 0–5,320 kyr BP. Primary spectral analysis; per-regime L1+L2+L3 fits at post-mpt / inhg-mpt / pre-inhg / lr04-full windows.
• CENOGRID composite — Westerhold, T. et al. (2020), Science 369, 1383. Cenozoic benthic δ¹⁸O + δ¹³C global composite spanning 0–67 Myr BP, ~5-kyr binning, multi-site radiometric anchor points. Deep-time generalization; demonstrates L3 step-component dominance at Cenozoic scale.
• EPICA Dome C CO₂ — Bereiter, B. et al. (2015), Geophys. Res. Lett. 42, 542. Antarctic ice-core atmospheric CO₂ composite spanning 0–805 kyr BP, native ~2-kyr resolution. Cross-proxy validation: same orbital lattice fits atmospheric CO₂ at R² = 0.83 with 32-line L1 alone (0.85 total with L2 carbon thermostat).
• CenCO2PIP CO₂ — Consortium 2023, Science 382, eadi5177. Bayesian multi-proxy CO₂ synthesis spanning 0–66 Ma at 100-kyr resolution. Deep-time CO₂ validation; L3 step amplitudes recover canonical Cenozoic CO₂ history (PETM, EOT, Mi-1, MMCT, iNHG, MPT) independently from the fit.
• Cheng 2016 Asian speleothem composite — Cheng, H. et al. (2016), Nature 534, 640. U-Th-dated Asian Monsoon δ¹⁸O, 0–640 kyr BP. Non-tuned chronology-bias control (places the dominant peak at the same FFT bin as orbitally-tuned LR04, refuting a systematic ~10 % dating offset).

See also

Climate cluster (this section’s siblings):

• Climate Summary — the synthesis this formula supports
• L1 Attribution — per-integer Berger vs Holistic mapping (the dual-attribution finding)
• Insolation Null Test — empirical anchor for the “lattice subsumes Berger” claim
• Related Work — position relative to recent peer-reviewed literature

Connected pages elsewhere on the site:

• Eigenfrequencies — the secular eigenmodes (g_j, s_j) used throughout this page, with the model-frame vs Berger-attribution comparison
• Obliquity & Inclination — Earth’s two-component obliquity construction (axial tilt + inclination tilt)
• Eccentricity — Earth’s eccentricity oscillation and the 100-kyr problem statement
• Why Earth Is Special — Earth’s reference frame duality
• Supporting Evidence — broader validation across multiple model domains
• Scientific Background — historical context for the 100-kyr problem
• Reproducibility scripts  — Python implementations of the LR04 spectral analyses, beat searches, and bootstrap routines used to produce the findings on this page
• Solar System Resonance Cycle period table — the 8 planets × 6 cycle types reference table used for direct-match cross-referencing