Temporal Equivalence Principle: Lunar Laser Ranging and the Nordtvedt Effect

1. Introduction

The Strong Equivalence Principle (SEP) is a cornerstone of General Relativity, stating that gravitational mass equals inertial mass and that the outcome of any local non-gravitational experiment is independent of the velocity and location of the freely-falling reference frame. A violation of the SEP would suggest that gravity is not purely metric in nature, but may involve additional fields that couple differently to different bodies.

General Relativity has achieved extensive empirical validation across multiple regimes: the anomalous perihelion precession of Mercury (42.98 arcsec/century, matched to 0.1%); light deflection around the Sun (1.75 arcsec, verified to ~1%); gravitational time dilation (GPS clock corrections of ~45 μs/day, confirmed to 0.1%); orbital decay in binary pulsars (indirect gravitational wave detection, Hulse-Taylor pulsar); and direct gravitational wave emission from merging black holes (LIGO, 2015). These empirical successes position GR as the effective theory within its domain of validity: regions where spacetime curvature remains moderate ($\Phi/c^2 \ll 1$) and quantum corrections are negligible.

The Equivalence Principle in GR derives from the universality of free fall: all test particles follow geodesics of the metric $g_{\mu\nu}$ independent of composition. This emerges geometrically from Einstein's field equations in the weak-field limit.

The present work extends this framework through the Temporal Equivalence Principle (TEP), a scalar-tensor theory where proper time becomes a dynamical field $\phi$ that couples universally to matter via conformal metric $\tilde{g}_{\mu\nu} = A^2(\phi)g_{\mu\nu}$, with $A(\phi)=\exp(\beta\phi/M_{\rm Pl})$. Scalar-tensor theories have a long history in gravitational physics, dating to Jordan (1955), Fierz (1956), and Brans-Dicke (1961), who first explored gravitational theories with dynamical scalar fields coupled to the metric. The Parametrized Post-Newtonian (PPN) formalism (Will 2014) provides a standardized framework for testing such theories through observable parameters including the Eddington parameters $\beta$ (nonlinearity) and $\gamma$ (spatial curvature), both equal to unity in GR.

TEP preserves the Weak Equivalence Principle through universal conformal coupling: all non-gravitational processes couple to the same matter metric, ensuring local experiments remain metric-compatible. Where TEP diverges from standard scalar-tensor theories is in its screening of Temporal Shear, which allows the theory to evade tight solar-system constraints while maintaining large bare couplings. The Cassini bound on PPN-$\gamma$ ($|\gamma - 1| < 2.3 \times 10^{-5}$; Bertotti et al. 2003) constrains the effective scalar coupling to $\alpha_{\rm eff} \lesssim 3 \times 10^{-3}$ in the Solar System. TEP Temporal Shear screening operates via the continuous spatial profile of the time field (Temporal Topology), in which high ambient density in deep potential wells suppresses the local field gradient (Temporal Shear). Bodies with high gravitational compactness $\Phi/c^2$ experience stronger suppression of Temporal Shear, yielding a vanishing field gradient and an effective scalar coupling $\alpha_{\rm eff} \ll \alpha_0$ (see the TEP framework paper, Paper 0, §7). GR is recovered exactly when Temporal Shear is uniform or when $\phi$ is spatially constant.

The Earth-Moon system occupies a critical regime where $\Delta(\Phi/c^2) \sim 10^{-10}$ — small enough that GR appears valid in local tests, yet large enough to produce differential coupling detectable through Lunar Laser Ranging.

The most precise test of the SEP is Lunar Laser Ranging (LLR), which measures the Earth-Moon distance with millimeter precision by timing laser pulses reflected from retroreflectors placed on the Moon by Apollo missions and Soviet Luna probes.

Over 50 years of LLR data have constrained the Nordtvedt parameter $\eta$ — which quantifies SEP violation — to $|\eta| \lesssim$ few $\times 10^{-4}$ (Williams et al. 2012; Müller et al. 2019), consistent with zero. These constraints arise from fitting $\eta$ as a free parameter in the complete orbital model.

The present analysis takes a complementary approach: it searches for a cos(elongation) modulation in the residuals of a GR-based ephemeris (INPOP19a) that assumes $\eta = 0$, targeting a TEP-specific suppression signal that may not be fully captured by the standard Nordtvedt parameterisation. Current LLR solutions leave centimeter-level O-C residuals after fitting all standard physical effects, providing the dataset for this search.

1.1 Manuscript Series Context

This work is Paper 17 in the TEP program. Earlier installments established the framework (Paper 0), clock-network phenomenology (Papers 1–3), and domain-specific Observable Response Coefficients (Papers 4–13). Three prior papers frame the present LLR analysis without duplicating its estimand.

The role of the present paper within the series is therefore not to restart the TEP argument from first principles, but to test whether the same time-field architecture has a measurable, suppressed Solar-System expression in the most precise available strong-equivalence-principle observable.

Series strand	Prior result	Role in this paper
Paper 0	Two-metric TEP framework with a dynamical time field, universal matter coupling, Temporal Topology, and Temporal Shear suppression.	Defines $\phi$, $A(\phi)$, the screened conformal sector, and the requirement that local Lorentz invariance is preserved while global time structure may remain dynamical.
Papers 1–3	Distance-structured and directionally structured clock-network correlations in GNSS products, with multi-center and raw-RINEX validation.	Establish the empirical motivation for treating precision timing residuals as a direct probe of Temporal Topology.
Paper 6	A candidate saturation scale $\rho_T \approx 20$ g/cm$^3$ from cross-scale Temporal Topology consistency.	Supplies the density/screening scale used to interpret why LLR should lie in the strongly screened Solar-System response regime.
Paper 9	A taxonomy of precision relativity tests distinguishing direct-fit, two-way, residual-channel, and clock-sector observables.	Explains why a standard direct-fit $\eta$ bound and a post-fit residual-channel synodic test are related but not identical estimands.
Papers 10–11	Observable Response Coefficients for pulsar spin-down and Cepheid timing/distance-ladder observables.	Motivate treating the LLR Nordtvedt parameter $\eta$ as a Solar-System response coefficient rather than a bare microscopic scalar coupling.
Paper 13	Weak-field Temporal Shear recovery in Gaia DR3 wide binaries.	Provides the complementary weak-field, weakly screened limit to the compactness-suppressed Earth-Moon test here.
Paper 17	LLR residual-channel search for a synodic Nordtvedt-like component.	Tests the screened Solar-System SEP sector of TEP using public INPOP19a and DE430 lunar-ranging residual archives.

Paper 5 (GTE) argued that Solar System tests remain compatible with TEP when Temporal Shear is suppressed in dense environments, with a microscopic Nordtvedt amplitude well below published direct-fit bounds. Paper 8 (SLR) extended the conformal-sector program to two-way optical ranging on passive geodetic targets, testing distance-structured residual structure orthogonal to active-clock microwave networks. Paper 9 (EXP) classified LLR as a reciprocity-even, orbital-dynamics constraint that does not directly target clock-sector loop holonomy or spatial clock correlations.

Paper 17 does not overturn those scope claims. It instead tests whether a synodic $\cos(D)$ footprint survives in post-fit O-C residuals after INPOP19a and DE430 solutions that assume $\eta=0$. The fitted Nordtvedt parameter is treated as the Solar System Observable Response Coefficient for LLR, not as the bare microscopic coupling. The hardware epoch, synthetic absorption, and toy Keplerian analyses supply the structural reason a static direct-fit $\eta$ and a residual synodic channel need not coincide when the coupling scales with the heliocentric environment. Paper 8 and Paper 17 therefore occupy complementary optical-ranging lines: SLR probes conformal residual coherence at terrestrial baselines; LLR probes suppressed-PPN differential free fall in the Earth-Moon system.

1.2 Claim Hierarchy and Scope

The evidence is organized deliberately so that the empirical result, the systematic stress tests, and the TEP interpretation are not conflated. The paper advances a strong claim, but in a fixed hierarchy.

Level	Claim	Role in the argument
1	Data integrity	Public INPOP19a and DE430 residual archives are processed through a deterministic, checksum-audited pipeline.
2	Residual-channel detection	A synodic $\cos(D)$ component is extracted from post-fit LLR residuals under the fixed estimator hierarchy.
3	Systematic resistance	The same residual-channel component is tested against station concentration, hardware epochs, outliers, autocorrelation, spectral controls, and nuisance regressors.
4	TEP interpretation	The sign, scale, and screened Solar-System response are interpreted using the Temporal Topology and Temporal Shear framework established in the preceding papers.
5	External closure	Source-level INPOP or DE430 integrator refits with $\eta$ free, together with independent range-level reductions, remain the critical external validation tests.

This hierarchy strengthens the inference by fixing what is being measured. The primary result is a residual-channel extraction of a Nordtvedt-like synodic component; the TEP claim is that this component has the sign, order of magnitude, screening regime, and dynamical structure expected for the Solar-System response of the same time-field theory tested elsewhere in the series.

The LLR Nordtvedt test in this paper probes a complementary aspect of TEP: through the suppressed PPN sector, the compactness-dependent effective coupling $\alpha_{\rm eff}$ could differ between Earth and Moon, producing a violation of the Strong Equivalence Principle. Earth's deeper gravitational potential ($\Phi_{\oplus}/c^2 \approx 7 \times 10^{-10}$) flattens Temporal Topology more strongly than the Moon's ($\Phi_{\rm Moon}/c^2 \approx 3 \times 10^{-11}$), suppressing Temporal Shear and yielding a smaller $\alpha_{\rm eff}$. This differential gradient suppression could lead to unequal free-fall rates in the Sun's field.

The quantitative prediction for the Nordtvedt parameter is informed by the TEP framework's Observable Response Coefficients, which quantify domain-specific astrophysical responses rather than a universal bare coupling. Preliminary results from related work in the same TEP framework report $\kappa_{\rm Cep} = (1.05 \pm 0.43) \times 10^6$ mag for Cepheid period-luminosity anomalies (Paper 11) and $\kappa_{\rm MSP} \sim 10^6$–$10^7$ for pulsar spin-down excess (Paper 10). These coefficients are distinct from the microscopic conformal coupling $\beta$ or scalar-tensor coupling $\alpha_0$, absorbing instead the full astrophysical response including environmental activation and transfer functions. The prediction $\eta \sim -10^{-4}$ emerges from the differential suppression geometry combined with the understanding that LLR operates in a more screened Solar System regime, yielding a smaller effective response than the unscreened galactic probes.

This analysis uses 26,207 raw LLR O-C residuals from five international laser ranging stations spanning 35 years of measurements (1984–2019), with 25,445 retained after 6σ-equivalent (MAD-based) outlier cleaning. The residuals are processed against the INPOP19a lunar and planetary ephemeris from the Paris Observatory (Geoazur). To eliminate synodic blurring and ensure millimeter-level coordinate precision, Moon-Sun elongation angles were computed using high-precision Skyfield/DE440 ephemerides rather than mean-phase approximations. The analysis searches for the predicted TEP Nordtvedt signal: a modulation of the form $\delta r = 13 \eta \cos(D)$, where $D$ is the Moon-Sun elongation angle.

The predicted amplitude for $\eta \sim 10^{-4}$ is $\sim$1–5 mm, smaller than the 9.5 cm per-observation residual RMS. This is expected for a suppressed PPN signal in precision LLR metrology: the signal explains only $\sim$0.1% of the total variance, but with $N = 25{,}445$ observations the standard error of the mean is $\sigma_r \approx 1/\sqrt{N} = 0.0063$, and the observed correlation coefficient $r = -0.0329$ is significant at $5.25\sigma$. The operative statistical leverage is not per-shot variance explained but synodic phase coherence across 35 years of independent observations. As hardware precision improved from $\sim$10 cm (PMT era) to $\sim$2 cm (C-SPAD era), the extracted physical amplitude did not scale to zero; it converged to a fixed negative boundary, consistent with a permanent modulation beneath the noise envelope.

Müller and Nordtvedt (1998) reported an unexplained synodic post-model residual proportional to $\cos(D)$ in 28 years of LLR data. The present sub-centimeter archive allows a cleaner extraction of the same footprint under explicit annual, monthly, and thermal controls.

Standard orbital periods are incommensurate with the synodic month. In the frequency domain, $\cos(D)$ therefore lies outside the classical multi-body mean-motion basis used in static direct-fit ephemerides, and a synodic coupling parameterized with $\eta$ fixed at zero is expected to project into post-fit residuals. Source-level integrator refits that leave $\eta$ free (Section 4.14) supply the critical dynamical closure; the hardware epoch, synthetic absorption, and toy Keplerian analyses supply the internal linearized and simulation checks on that channel.

2. Theoretical Framework: TEP and the Nordtvedt Effect

The Temporal Equivalence Principle (TEP) is a scalar-tensor theory in which proper time becomes a dynamical field $\phi$ that couples to the local mass density. In this framework, matter couples to a conformal metric $\tilde{g}_{\mu\nu} = A^2(\phi) g_{\mu\nu}$, where $A(\phi) = \exp(\beta\phi/M_{\rm Pl})$ is the conformal factor. The rate at which proper time accumulates depends on the local value of $\phi$, which in turn depends on the ambient matter density through screening of Temporal Shear. This operates via the continuous spatial profile of the time field (Temporal Topology), in which high ambient density in deep potential wells suppresses the local field gradient (Temporal Shear), naturally attenuating fifth-force effects in dense environments while allowing the field to remain light and long-ranged in low-density regions.

2.1 The Nordtvedt Effect

In scalar-tensor theories, bodies with different gravitational binding energy per unit mass will fall at different rates in an external gravitational field. This is the Nordtvedt effect, first derived by Kenneth Nordtvedt (Nordtvedt 1968) as a test of the Strong Equivalence Principle. The Nordtvedt parameter $\eta$ quantifies the strength of SEP violation:

where $\beta$ and $\gamma$ are the post-Newtonian parameters, $\alpha_1$ and $\alpha_2$ are preferred-frame parameters, $\xi$ is the Whitehead parameter, and $\zeta_1$, $\zeta_2$ are conservation-law parameters. In General Relativity ($\beta = \gamma = 1$, all others zero), $\eta = 0$ exactly. In standard scalar-tensor theories such as Brans-Dicke, the Nordtvedt parameter is approximately $\eta_{\rm BD} \approx 4\alpha_0^2/(1+\alpha_0^2)^2$ where $\alpha_0 \equiv 2\beta/M_{\rm Pl}$ is the bare scalar-tensor coupling. The Cassini bound on PPN-$\gamma$ constrains $\alpha_0 \lesssim 3 \times 10^{-3}$ in the unsuppressed (high Temporal Shear) regime (Will 2014).

In TEP, the Nordtvedt effect arises from a compactness-dependent scalar coupling. TEP's universal conformal coupling $A(\phi)$ preserves the Weak Equivalence Principle (all non-gravitational processes couple to the same matter metric $\tilde{g}_{\mu\nu}$, consistent with the TEP framework paper, Paper 0, \S7.2). However, the Strong Equivalence Principle may be violated through the suppression of Temporal Shear.

For a self-gravitating body in the deeply suppressed regime, the effective scalar coupling is suppressed relative to the bare coupling $\alpha_0$ as Temporal Shear vanishes in the deep potential well. The suppression of Temporal Shear (vanishing field gradient) reduces the effective coupling to $\alpha_{\rm eff} \ll \alpha_0$. The degree of gradient suppression scales with the body's surface gravitational compactness $\Phi/c^2 = GM/Rc^2$ (see the TEP framework paper, Paper 0, §7, for a detailed discussion of the relationship between compactness and gradient suppression).

For two bodies with different compactness (Earth and Moon), the differential acceleration in an external gravitational field $a_0$ produces a Nordtvedt violation. The effective coupling difference is governed by the differential flattening of Temporal Topology:

Substituting the compactness scaling and defining the Nordtvedt parameter as $\eta \equiv \delta a/a_0$ yields the TEP prediction:

where $\Phi_{\oplus}/c^2 \approx 7 \times 10^{-10}$ is Earth's compactness and $\Phi_{\rm Moon}/c^2 \approx 3 \times 10^{-11}$ is the Moon's compactness. The quadratic dependence on compactness follows from two steps. First, gradient suppression gives an effective scalar coupling that scales linearly with the surface compactness:

Second, in the Damour–Esposito-Farèse (DEF) scalar-tensor framework the Nordtvedt parameter is the squared difference of the effective couplings:

Substituting Eq.~(\ref{eq:alpha_eff_scaling}) into Eq.~(\ref{eq:eta_from_alpha_eff}) yields the TEP prediction in Eq.~(\ref{eq:tep_prediction}). Earth's deeper potential well flattens Temporal Topology more strongly than the Moon's, suppressing Temporal Shear and yielding a significantly smaller effective coupling.

Two theoretically distinct predictions are evaluated. The compactness-squared expression above uses only the bare microscopic coupling $\alpha_0$ and the surface compactness differential. With the Cassini-bound $\alpha_0 \lesssim 3\times10^{-3}$ and $\Delta(\Phi^2) \approx 5\times10^{-19}$, this gives $\eta \approx 4\times10^{-24}$ — essentially zero. This baseline is not a physical prediction for LLR; it is a consistency check that confirms the measured $\eta$ cannot arise from the bare coupling alone. $\eta$ is an emergent observable response coefficient, not a microscopic parameter.

The phenomenological TEP prediction is obtained from the volumetric suppression model, which integrates interior density profiles rather than surface compactness:

Using a simplified PREM density model and $\rho_T \approx 20$ g/cm³ from GNSS calibration (Paper 6), the Earth-Moon shielding differential is $\Delta\langle(\rho/\rho_T)^3\rangle \approx 3.2\times10^{-2}$, giving $\eta \approx -9.7\times10^{-5}$. The current TEP framework robustly predicts the sign of the effect under compactness-driven Temporal Shear screening and gives an order-of-magnitude expectation in the $10^{-4}$ regime under this phenomenological volumetric model. It does not yet provide a first-principles, parameter-free prediction of $\eta$. The headline precision-weighted estimate $\eta = -3.91\times10^{-4}$ and the naïve OLS robustness check $\eta = -3.18\times10^{-4}$ lie in the same order-of-magnitude regime, but the phenomenological model uncertainties — PREM simplification, the exponent in $(\rho/\rho_T)^3$, $\rho_T = 20 \pm 8$ g/cm³, and the upper-bound nature of the Cassini $\alpha_0$ constraint — preclude a precision comparison at this stage.

The measured $\eta$ is also ~8.8× larger than the standard scalar-tensor upper-bound prediction $\eta \approx 4\alpha_0^2 \approx 3.6\times10^{-5}$ (with $\alpha_0 \lesssim 3\times10^{-3}$ from Cassini), consistent with the interpretation that the LLR observable response absorbs contributions from the TEP screening mechanism beyond bare scalar-tensor physics. The measured $\eta$ is much smaller than preliminary $\kappa_{\rm MSP} \sim 10^6$–$10^7$ and $\kappa_{\rm Cep} \sim 10^6$ from related work in the same framework (Papers 10 and 11), which is consistent with the screening mechanism: LLR is in a more screened regime (Solar System) compared to globular clusters and galactic disks, so the effective response should be smaller. No separate $\kappa_{\rm LL}$ parameter is required; $\eta$ itself is the observable response coefficient for LLR. The prediction hierarchy is therefore: (i) bare compactness-squared null check, (ii) volumetric phenomenological order-of-magnitude band, (iii) measured LLR response coefficient $\eta$.

2.2 Predicted Signal in LLR

The TEP Nordtvedt effect predicts a modulation of the Earth-Moon range as the Earth-Moon system orbits the Sun. The predicted range perturbation is:

The predicted amplitude scales linearly: $|\delta r| = 13 |\eta|$ meters. For the order-of-magnitude estimate $\eta \sim 10^{-4}$, this gives $\sim$1.3 mm. The headline precision-weighted estimate ($\eta = -3.91 \times 10^{-4}$, Section 4.1) predicts an amplitude of $\sim$5.3 mm, while the leverage-excised robustness check ($\eta \approx -3.6 \times 10^{-4}$) predicts $\sim$4.18 mm and the naïve cosD-only OLS robustness check ($\eta \approx -3.18 \times 10^{-4}$) yields $\sim$4.7 mm. These amplitudes are small compared to the centimetre-level residual RMS, but with 26,207 observations over 35 years, statistical averaging provides sensitivity well below the per-observation noise floor.

2.3 Continuous Gradient Suppression and the Critical Density

TEP incorporates a Temporal Shear screening mechanism that suppresses scalar field effects in dense environments. The suppression transition is governed by the body's gravitational compactness ($\Phi/c^2$), which determines the degree of Temporal Shear suppression and thus the effective scalar coupling. Screening in TEP is tiered, not a single ambient-density switch. The asymptotic saturation scale $\rho_T \approx 20$ g/cm³ (Paper 6) marks where extended self-gravitating sources flatten Temporal Topology and suppress Temporal Shear in the Solar System. At galactic mean densities and in globular-cluster environments, gradient coherence can remain active despite local $\rho \ll \rho_T$, producing the larger observable response coefficients $\kappa$ reported in Papers 10 and 11. LLR therefore probes the strongly screened planetary channel: $\eta$ is the Solar System response coefficient, expected to be smaller in magnitude than galactic $\kappa$ while sharing the same screening ontology.

In the Earth-Moon orbit, the interplanetary ambient density ($\rho_{\rm amb} \sim 10^{-23}$ g/cm³) is far below $\rho_T$, but Earth and Moon remain differentially screened through their distinct compactness and interior shielding profiles. The TEP framework further predicts that this coupling is not static but environmentally modulated by the local scalar potential $\phi(r_\odot)$, which scales with heliocentric distance. Continuous TEP suppression near a threshold transition implies that the effective Nordtvedt parameter $\eta$ varies with heliocentric distance $r_\odot$. When the Earth-Moon system moves through regions where the background scalar field approaches a suppression threshold, the effective coupling modulates as:

where $\eta_0$ is the baseline Nordtvedt parameter at mean orbital distance $r_{\rm mean}$, $e_\oplus \approx 0.0167$ is Earth's orbital eccentricity, and $m$ is the modulation depth characterizing the nonlinear threshold response. For $m \approx 1$ (consistent with the observed sign-flip), the orbital modulation generates diagnostic sidebands at frequencies $D \pm l'$ (where $l'$ is the orbital longitude), spectrally orthogonal to classical multi-body resonances. This provides a unique signature distinguishing TEP from standard PPN or systematic artifacts.

The TEP framework further predicts that the coupling depends not only on the scalar gradient magnitude (distance) but also on the rate at which the Earth-Moon system traverses the temporal topology. A body moving through a spatial scalar gradient $\nabla\phi$ at velocity $\mathbf{v}$ experiences an effective temporal shear rate $d\phi/dt \approx \mathbf{v} \cdot \nabla\phi \approx v_r \, |\nabla\phi|$, where $v_r$ is the heliocentric radial velocity. For a Kepler orbit with small eccentricity, distance $r$ and radial velocity $v_r$ are approximately in quadrature ($90^\circ$ out of phase), making them statistically distinguishable predictors. The full dynamical modulation therefore takes the form:

where $\bar{v}_r$ is the characteristic radial velocity scale ($\approx 0.5$ km/s for Earth's orbit) and $m_v$ is the velocity modulation depth. A joint fit to both $r_\odot$ and $v_r$ determines whether the temporal topology is purely static ($m_v = 0$) or dynamically responsive ($m_v \neq 0$). The heliocentric radial velocity analysis tests this prediction directly, using DE440 ephemeris velocities computed via Skyfield for every LLR observation epoch.

A further TEP-motivated probe uses a fixed celestial axis tied to the Planck 2018 dipole as an operational reference: $(l, b) = (264.02^\circ, 48.25^\circ)$ in galactic coordinates, equivalent to $(\alpha, \delta) = (168.14^\circ, -7.22^\circ)$ in J2000 equatorial coordinates, with amplitude $v_{\rm CMB} \approx 369$ km/s for the kinematic dipole. Some scalar-field embeddings single out a large-scale rest frame; the LLR regressions reported here do not, by themselves, establish Lorentz violation or uniqueness of that frame. They test whether residual-channel structure is consistent with anisotropic coupling defined relative to $\hat{\mathbf{n}}_{\rm CMB}$ and with other fixed axes used as controls (Section 4.12.2). Two operational predictions tied to $\hat{\mathbf{n}}_{\rm CMB}$ are:

First, an annual velocity projection: Earth's orbital velocity ($\sim 30$ km/s) projects onto the CMB dipole direction with a sinusoidally varying parallel component $v_\parallel(t) = \mathbf{v}_{\rm orb} \cdot \hat{\mathbf{n}}_{\rm CMB}$. The CMB dipole lies at ecliptic longitude $\approx 173^\circ$, offset by approximately $70^\circ$ from the perihelion longitude ($\approx 103^\circ$), making the annual velocity projection phase-shifted relative to the heliocentric distance modulation and therefore statistically distinguishable.

Second, a monthly orientation anisotropy: the Earth-Moon line sweeps across the celestial sphere with synodic period. If an effective gradient projection onto $\hat{\mathbf{n}}_{\rm CMB}$ is present in the same channel, the coupling should depend on the cosine of the angle $\theta$ between the Earth-Moon vector and that axis:

where $\cos\theta_{\rm EM-CMB} = \hat{\mathbf{r}}_{\rm EM} \cdot \hat{\mathbf{n}}_{\rm CMB}$. Because the Earth-Moon line rotates on a monthly timescale while the synodic signal operates on a 29.5-day period, the two modulations are at different frequencies and can be separated in a joint fit. The CMB dipole orientation analysis tests both predictions directly.

For the Earth-Moon system, the TEP potential produces a density-dependent effective mass $m_{\rm eff}(\rho)$ that grows with local density (see the TEP framework paper, Paper 0, §4.31), limiting the scalar field range inside massive objects. For a body in the strongly suppressed regime, this translates to a suppressed effective coupling $\alpha_{\rm eff} \ll \alpha_0$, where the degree of suppression scales with the body's surface gravitational potential (see the TEP framework paper, Paper 0, §7). Earth ($\Phi_{\oplus}/c^2 \approx 7 \times 10^{-10}$) experiences stronger Temporal Topology flattening than the Moon ($\Phi_{\rm Moon}/c^2 \approx 3 \times 10^{-11}$), giving Earth substantially stronger self-suppression and a smaller effective coupling to $\phi$. This large differential in gravitational compactness — and the resulting differential $\alpha_{\rm eff}$ — provides the physical basis for a Nordtvedt effect in TEP.

4. Results: Candidate TEP Signal in LLR Residuals

A correlation analysis between the LLR O-C residuals and the predicted TEP Nordtvedt signal modulation $\cos(D)$, where $D$ is the Moon-Sun elongation angle. Section 4.1 states the independent evidence spine; the subsections that follow develop pooling, regression, station-specific systematic testing, and cross-ephemeris validation; later sections treat systematics, spectral specificity, environmental modulation, orthogonality, and robustness testing without reopening the headline synodic estimand. The table below gives the referee path through §4.

Recommended read order for §4.

Read order	§	Focus
1	4.1	Independent synodic spine
2	4.25	Multiplicity control
3	4.3	Cross-station pooling; leave-one-station-out meta (Figure 5)
4	4.10	Station-balance power stress
5	4.5–4.18	Station-specific systematic testing; core regression and per-station detail
6	4.9	Hardware epochs
7	4.13, 4.14	DE430 comparison; linearized post-fit $\eta$ extraction
8	4.20.5, 4.15–4.32 (incl. 4.28.1)	Uncertainty calibration; systematics, CV / nuisance transport, clean subset, orthogonality simulations
9	4.25.*	Conditional CMB / heliocentric joint layer

4.1 Independent Evidence Spine

Estimand	$\eta$ ($\times 10^{-4}$)	Significance	$N$	Role
Precision-weighted full-systematic (primary headline)	$-3.91 \pm 0.56$	$6.94\sigma$	25,445	Primary headline estimand — all observations retained
Wild cluster bootstrap, Webb weights	$-4.06 \pm 0.73$	$5.56\sigma$ ($95\%$ CI $[-5.37, -2.73]$)	25,445	Co-primary small-$G$ uncertainty companion; OLS full-systematic
Full-systematic OLS (no excision)	$-4.06 \pm 0.66$	$6.17\sigma$ ($6.52\sigma$ cluster)	25,445	Sensitivity upper bound
Cook's-Distance-excised full-systematic OLS	$-3.87 \pm 0.50$	$7.82\sigma$ ($8.65\sigma$ cluster)	23,837	Secondary robustness check — leverage diagnostic
Common-$\eta$ + station systematics	$-4.32 \pm 0.68$	$6.40\sigma$	25,445	Pooling
Phase-locked new/full-moon differential	$-5.95 \pm 1.01$	$5.91\sigma$	—	Robustness
Precision-weighted full-systematic after Bonferroni / BH	$-3.91 \pm 0.56$	$6.75\sigma$ (both)	25,445	Multiplicity on four independent hypotheses
Full-systematic OLS after Bonferroni / BH	$-4.06 \pm 0.66$	$5.95\sigma$ ($6.06\sigma$ BH)	25,445	Multiplicity (unweighted sensitivity member)
cosD-only OLS	$-3.18 \pm 0.61$	$5.25\sigma$	25,445	Baseline
LOSO meta (four powered exclusions)	$-4.21 \pm 0.33$	$12.8\sigma$	25,445 (full sample)	Cross-station leverage; $I^2 = 0.0\%$
Prediction-interval calibration	—	68% nominal $\to$ 85.9% observed; $\chi^2_{\rm red}=0.48$	25,445	Conservative $\sigma$; LOSO conformal 95% excludes $\eta=0$

Four independent synodic tests converge before the station-by-station narrative: (i) the precision-weighted full-systematic consensus ($\eta = -3.91 \times 10^{-4}$, $6.94\sigma$), which retains every observation and objectively down-weights high-variance epochs by station-level inverse variance without deleting any data; (i-a) a wild cluster bootstrap with Webb 6-point weights on the same full-systematic design gives $\eta = -4.06 \times 10^{-4}$ with a $95\%$ percentile interval $[-5.37, -2.73] \times 10^{-4}$ ($5.56\sigma$; wild cluster bootstrap), confirming that the detection survives even when the small-$G$ ($G=5$) unreliability of analytical cluster-robust standard errors is corrected by resampling; (ii) a common-$\eta$ mixed model with station-specific systematics ($\eta = -4.32 \times 10^{-4}$, $6.40\sigma$; $F(4, 25{,}410) = 1.19$, $p = 0.31$); (iii) a phase-locked new/full-moon differential ($\eta = -5.95 \times 10^{-4} \pm 1.01 \times 10^{-4}$, $5.91\sigma$); and (iv) a frequency-specific null scan with no significant detections at 55 tested non-synodic factors after correction. Leave-one-station-out meta-analysis ( Section 4.10) gives $\eta_{\rm meta} = -4.21 \times 10^{-4}$ at $12.8\sigma$ with $I^2 = 0.0\%$ across four powered exclusions. The full-systematic model without weighting or excision ($\eta = -4.06 \times 10^{-4}$, $6.17\sigma$ / $6.52\sigma$) is reported as a sensitivity upper bound on the same synodic estimand. A formal Cook's-Distance leverage diagnostic excises 1,608 high-leverage points and returns a consistent $\eta = -3.87 \times 10^{-4}$ ($7.82\sigma$ / $8.65\sigma$), confirming that the signal is not driven by a small subset of influential observations. Multiplicity control is applied to four independent hypotheses while treating sixteen sensitivity analyses as bounds on the same synodic claim (Section 4.23). The canonical estimator hierarchy consolidates these into one unified framework. A matched-window ephemeris comparison on 2014–2018 adds: INPOP19a and DE430 both yield negative $\eta$ at $10.0\sigma$ and $5.96\sigma$ respectively (cosD-only), with $\Delta\eta = +3.42 \times 10^{-4}$ ($2.77\sigma$), and the canonical full-systematic model gives $7.73\sigma$ and $5.04\sigma$ on the same window ($\Delta\eta = +3.10 \times 10^{-4}$, $2.49\sigma$).

A fifth strand—strong corroborative dynamical structure in the same residual channel—is developed in Sections 4.28.1–4.14.2 (read order 9). After the synodic fit is fixed, heliocentric radial velocity and CMB dipole orientation improve the joint model strongly ($\Delta\text{AIC} = -119.2$; $F(1, 25{,}440) = 121.49$, $p = 1.1 \times 10^{-16}$), and falsification testing excludes aliasing, multicollinearity, and permutation artifacts for the orientation term. Sky-scrambling Monte Carlo (100,000 draws) is only marginal ($p_F \approx 0.096$; correlation-matched $p_F \approx 0.266$; composite specificity $Z = 1.45$, $p = 0.074$), while a synodic phase null rejects chance alignment ($p_{F,{\rm eff}} \approx 6 \times 10^{-4}$). A dual-axis fit (test H) shows galactic and Planck orientation predictors are $98.4\%$ collinear and do not form two independent channels when fit together. The falsification testing is scored as 4/5 diagnostics: strong directional-anatomy and phase-coupling support, not proof of a unique CMB-frame origin. The Planck-axis joint layer is therefore not a replacement estimator for the headline $\eta$, but it is retained as corroborative fixed-sky structure on the same residual channel. DE430 fragility (Section 4.28) is an independent cross-ephemeris replication of the INPOP19a headline estimand.

4.19 Multiple-Testing Consolidation

The pipeline reports 51 significance measures, of which four are treated as independent hypotheses and 47 as sensitivity analyses on the same synodic hypothesis (including 26 era$\times$station$\times$lunation grid cells and the blind year hold-out, registered under category interaction_grid or validation). Bonferroni and Benjamini–Hochberg correction applied only to the independent set leave the precision-weighted full-systematic headline estimand at $6.75\sigma$ and $6.75\sigma$ respectively, and the unweighted full-systematic OLS sensitivity bound at $5.95\sigma$ and $6.06\sigma$. Station-level and diagnostic regressions are not multiplied into the headline $\eta$ claim.

4.3 Cross-Station Validation: Defense Against Single-Station Dominance

The Grasse station contributes 73.7% of the 6$\sigma$-cleaned archive (18,742 of 25,445 observations), so the first systematic test is whether the extracted signal is a Grasse-specific artifact. Cross-station validation and controlled pooling (Sections 4.1, 4.3, 4.5, and 5.5) show the pooled signal does not reduce to a Grasse-only artifact. A common-$\eta$ mixed model with station-specific systematics yields $\eta = -4.32 \times 10^{-4}$ ($6.40\sigma$) with no evidence for station-specific deviations ($F(4, 25{,}410)=1.19$, $p=0.31$). A wild cluster bootstrap with Webb 6-point weights on the full-systematic OLS design gives a $95\%$ percentile interval $[-5.37, -2.73] \times 10^{-4}$ ($5.56\sigma$), while a station block bootstrap on the precision-weighted design gives $P(\eta<0)=0.945$ (94.5%). Individual stations and balance-enforced subsamples remain below the $3\sigma$ powered-detection threshold where phase coverage or $N$ is limited. Universality is tested at the pooled level: the common-$\eta$ mixed model is the primary cross-station significance statement. Per-station full models are informative but sub-threshold where phase coverage or $N$ is limited (Section 4.20); they are not treated as five independent discovery channels.

APO concurrent-validation robustness check. Apache Point Observatory (USA) shows consistent negative $\eta$ ($\eta = -2.39 \times 10^{-4} \pm 8.65 \times 10^{-5}$, $N = 2,595$, SNR = $2.77\sigma$; concurrent-validation robustness check) entirely separate from Grasse. While individually below the conventional $3\sigma$ threshold, APO contributes station-consistent support for the signal's geographic coherence.

Cross-Station Prediction Validation. A strong test of signal authenticity is whether APO's fitted amplitude can predict Grasse's phase-locked residuals. Using APO's fitted model to predict Grasse observations yields a correlation of $r = 0.0357$ with $p = 6.82 \times 10^{-7}$. This significant cross-prediction indicates that the anomaly phase-locks coherently across separate observatories on separate continents (USA and France), with entirely different hardware systems and operational teams.

Precision-weighted regression (robustness check). To consolidate multi-station data without manual station excision, precision-weighted regression (weighting by inverse station variance) provides a cross-observatory diagnostic. This method yields $\eta_{\rm WLS} = -3.21 \times 10^{-4} \pm 5.12 \times 10^{-5}$ at $6.27\sigma$ (precision-weighted robustness check). The precision-weighted approach naturally down-weights low-precision stations without manual intervention: Haleakala (13.8 cm RMS, N=737), McDonald2 (severe phase truncation), and Matera are automatically assigned lower weights due to their higher variance, while high-precision stations (Grasse, APO) dominate the fit. The weighting is determined objectively by measurement precision, not subjective judgment. The resulting negative estimate agrees with the canonical full-systematic hierarchy and shows that low-precision station noise is not required to generate the signal, while the explicit 059 warnings retain the remaining Grasse-leverage risk.

With the single-station critique examined under the estimator hierarchy above, the analysis proceeds to the full correlation analysis.

4.10 Station-Balanced TEP Stress Test

Station-balance analysis directly tests whether Grasse dominance (74% of the raw archive) drives the pooled synodic fit. The primary estimand is the same full-systematic model as the headline precision-weighted regression; cosD-only fits are reported in parallel as a lower-power diagnostic.

• Full archive. Full-systematic $\eta = -4.06 \times 10^{-4} \pm 6.58 \times 10^{-5}$ ($6.17\sigma$; $N = 25{,}445$).
• Equal-$N$ subsample. Each station contributes 346 observations ($N = 1{,}710$ after cleaning). Full-systematic $\eta = -7.6 \times 10^{-5} \pm 3.99 \times 10^{-4}$ ($0.08\sigma$); cosD-only $0.51\sigma$.
• Grasse-capped subsample. Grasse is downsampled to APO's count ($N = 9{,}283$ after cleaning). Full-systematic $\eta = -2.77 \times 10^{-4} \pm 1.36 \times 10^{-4}$ ($2.04\sigma$); cosD-only $1.77\sigma$.
• Balanced bootstrap. Two hundred equal-$N$ iterations on the full-systematic estimand give mean SNR $0.51\sigma$ with a 95% coefficient interval that includes zero.
• Stratified equal-$N$ by environmental bins. Random equal-$N$ draws may suffer from epoch-concentration bias: small stations have different heliocentric-distance and CMB-orientation phase histories. The stratified analysis stratifies by heliocentric-distance quintile and CMB-orientation quintile, then draws equal observations per station within each combined bin. This eliminates the environmental phase-space confounding that diluted the random equal-$N$ power. Stratified equal-$N$ returns $\eta = -2.73 \times 10^{-4} \pm 3.66 \times 10^{-4}$ ($0.75\sigma$; $N = 1{,}389$). The 200-iteration bootstrap over stratified draws gives mean SNR $0.83\sigma$ with a 95% CI that includes zero. Stratification removes the environmental phase-space confounding that diluted the random equal-$N$ power, but the design remains underpowered by construction.

The headline precision-weighted extraction therefore survives on the full archive, while enforced station balance removes the statistical power needed for an independent subsample claim. This is a power stress test, not a second co-equal discovery channel, and it does not replace the pooled mixed-model universality test in Section 4.10. The stratified equal-$N$ analysis tests whether the equal-$N$ null was merely a sampling-efficiency failure rather than genuine station dominance. The stability report confirms that random equal-$N$, stratified equal-$N$, and the 200-draw stratified bootstrap all retain the negative sign of the full-archive $\eta = -4.06 \times 10^{-4}$ ($6.52\sigma$ cluster-robust), while SNR falls to $0.29\sigma$ (random) and $0.75\sigma$ (stratified); stratification improves SNR by a factor $2.6$ over random equal-$N$ without reaching discovery threshold—consistent with underpowered balance, not sign instability.

Quantitative power on the partitions. If the physical Nordtvedt amplitude were at the headline value $|\eta| \approx 3.9\times 10^{-4}$, the equal-$N$ design ($N_{\rm used}=1{,}710$) carries a full-systematic standard error $\sigma_\eta \approx 3.99\times 10^{-4}$ (JSON output), so the expected signal-to-noise ratio would be only $|\eta|/\sigma_\eta \approx 1.0$—well below conventional discovery thresholds. The observed SNR ($0.19\sigma$ on that subsample) is therefore compatible with sampling noise in a partition that is underpowered by construction, not with a refutation of the full-archive estimand. The Grasse-capped partition ($N_{\rm used}=9{,}283$, $\sigma_\eta \approx 1.36\times 10^{-4}$) would yield $|\eta|/\sigma_\eta \approx 3.0$ at the same hypothetical amplitude— marginal $3\sigma$ class sensitivity—while the 200-draw equal-$N$ bootstrap reports mean SNR $0.51\sigma$ with a 95% interval on $\eta$ straddling zero, consistent with enforced balance diluting Grasse's phase coverage and precision.

Leave-one-station-out meta-analysis

Equal-$N$ subsampling is underpowered by construction. A more informative test is leave-one-station-out (LOSO) meta-analysis: the full-systematic model with cluster-robust standard errors is fit five times, each time excluding one station. The resulting five $\eta$ estimates are combined in an inverse-variance meta-analysis with Cochran's $Q$ heterogeneity test.

• Full sample: $\eta = -4.06 \times 10^{-4} \pm 6.58 \times 10^{-5}$ ($6.17\sigma$)
• Excluding APO: $\eta = -4.12 \times 10^{-4} \pm 7.29 \times 10^{-5}$ ($5.65\sigma$); shift $\Delta = +0.10\sigma$ [powered]
• Excluding Grasse: $\eta = -1.65 \times 10^{-4} \pm 1.80 \times 10^{-4}$ ($0.92\sigma$); shift $\Delta = +3.88\sigma$ [underpowered]
• Excluding Haleakala: $\eta = -4.18 \times 10^{-4} \pm 6.43 \times 10^{-5}$ ($6.49\sigma$); shift $\Delta = +0.19\sigma$ [powered]
• Excluding Matera: $\eta = -4.18 \times 10^{-4} \pm 6.60 \times 10^{-5}$ ($6.34\sigma$); shift $\Delta = +0.20\sigma$ [powered]
• Excluding McDonald2: $\eta = -4.32 \times 10^{-4} \pm 6.45 \times 10^{-5}$ ($6.71\sigma$); shift $\Delta = +0.42\sigma$ [powered]

Four of five exclusions are individually powered (SNR $> 3\sigma$) and all four return the same negative sign as the full sample. The maximum powered shift is only $0.40\sigma$ (McDonald2 exclusion). Excluding Grasse reduces SNR to $0.92\sigma$ because the remaining four stations collectively lack the precision and phase coverage for an independent subsample claim, not because Grasse fabricated the signal. The inverse-variance meta-analysis of the four powered exclusions gives $\eta_{\rm meta} = -4.21 \times 10^{-4} \pm 3.29 \times 10^{-5}$ at $12.8\sigma$, with Cochran's $Q = 0.050$ ($\mathrm{df}=3$) and $I^2 = 0.0\%$. Zero heterogeneity confirms that the four independently powered station subsets agree on a common negative amplitude. Station dominance does not refute the signal; it reflects the expected precision-weighting of a genuine physical effect.

4.5 Testing for Station-Specific Artifacts

The station-specific systematic test stress-tests the Grasse-specific systematic hypothesis on the canonical 6$\sigma$-cleaned sample: that the extraction is driven by a Grasse-only effect perfectly correlated with $\cos(D)$. Grasse contributes 73.7% of the cleaned archive (18,742 of 25,445 observations); if the pooled $\eta$ were driven by a Grasse-only systematic, that systematic would need amplitude $0.70$ cm. This is $2.0\times$ larger than the largest known systematic projection (ephemeris differences, 0.36 cm) and exceeds every known systematic amplitude from the data-driven budget.

A three-way partition test shows that the Grasse component is strong but does not by itself close the station-leverage question. Grasse-only: $\eta = -4.58 \times 10^{-4} \pm 6.72\times 10^{-5}$ ($6.81\sigma$); non-Grasse: $\eta = -1.65 \times 10^{-4} \pm 1.80\times 10^{-4}$ ($0.92\sigma$); pooled: $\eta = -3.98 \times 10^{-4} \pm 6.48\times 10^{-5}$ ($6.14\sigma$). The pooled $\eta$ lies between the two partition estimates, consistent with a negative archive-level signal diluted by lower-precision non-Grasse stations, but the non-Grasse partition remains underpowered as an independent discovery claim.

A Grasse $\times$ $\cos(D)$ interaction term yields $t = -0.42$ ($p = 0.676$), providing no evidence that Grasse has a differential $\cos(D)$ coefficient. Monte Carlo station-dominance (5,000 random station subsets) places the Grasse SNR at the $100^{\\rm th}$ percentile, as expected from its precision and sample share. The test therefore rules against a simple Grasse-specific differential $\cos(D)$ systematic, while explicitly retaining material station-leverage risk.

Equal-weighted and balanced-station reweighting. To test whether precision-weighting or sample-size imbalance drives the pooled sign, two additional constructions are reported. An equal-weighted meta-analysis fits cosD-only $\eta$ per station and averages the four station estimates with equal weight (not inverse variance, which would re-weight Grasse). The result is $\eta = -1.21 \times 10^{-4} \pm 2.61 \times 10^{-4}$ ($0.47\sigma$), underpowered because small-station errors dominate the unweighted mean, but the sign remains negative. A balanced-station-N pooled analysis downsamples every station to the smallest station count ($N = 345$ per station, total $N = 1{,}725$) and fits the full-systematic model on the balanced pool: $\eta = -1.47 \times 10^{-4} \pm 3.70 \times 10^{-4}$ ($0.40\sigma$), again underpowered but sign-consistent. Neither reweighting reverses the sign; both are underpowered by construction because enforced balance discards the majority of Grasse's precision and phase coverage.

Grasse-conditioned estimand (full systematic design). Three complementary fits isolate Grasse leverage while retaining the full-systematic nuisance model with cluster-robust errors: (i) non-Grasse direct fit on the four smaller observatories; (ii) pooled common-$\eta$ mixed model with station-specific systematics; and (iii) pooled fit with an explicit Grasse$\times\cos(D)$ interaction before reading $\eta$. All three return negative $\eta$ with the same sign as the pooled reference; the non-Grasse partition remains underpowered as an independent discovery claim, but no estimator reverses sign when Grasse-specific structure is absorbed explicitly.

4.6 Gaussian Process Non-Parametric Extraction

To test whether the $\cos(D)$ modulation shape is genuinely sinusoidal or merely the best sinusoidal fit to a non-sinusoidal artifact, A non-parametric Gaussian Process extraction performs a non-parametric Gaussian Process (GP) extraction on 48 elongation bins. The GP uses a periodic $\mathrm{ExpSineSquared}(\ell=5, p=5)$ plus RBF kernel with learned white noise, imposing no functional form on the signal shape.

The GP recovers amplitude $1.13$ cm at phase $232.3^{\circ}$, corresponding to $\eta_{\rm GP} = -5.32 \times 10^{-4}$. Shape fidelity to a pure sinusoid is excellent: $R^2 = 0.985$ on the sine-component projection. The GP amplitude is $67.5\%$ larger than the parametric OLS estimate ($\eta_{\rm OLS} = -3.18 \times 10^{-4}$), consistent with subsampling variance across the 48-bin representation. The key result is that the non-parametric model confirms a coherent periodic structure locked to the synodic phase, not an arbitrary shape that happens to correlate with $\cos(D)$.

4.13 Systematic Amplitude Sensitivity

The systematic sensitivity analysis quantifies systematic risks to the signal. Comparisons in this section use the full-systematic OLS sensitivity estimand ($\eta = -4.06 \times 10^{-4}$) for like-with-like Monte Carlo and hold-out tests; the primary headline remains the precision-weighted full-systematic result ($\eta = -3.91 \times 10^{-4}$ at $6.94\sigma$; .

The analysis first quantifies how large each known systematic would need to be to fully explain the observed $|\eta| = 4.06 \times 10^{-4}$. The required systematic amplitude is $0.53$ cm. The smallest exclusion ratio is $1.5\times$ (ephemeris differences); atmospheric, thermal, instrumental, and tidal effects require $6.6\times$, $9.5\times$, $54.0\times$, and $57.9\times$ their known amplitudes respectively. Every known systematic is too small by at least a factor of $1.5$.

Monte Carlo null simulations (2,000 deterministic draws per systematic) inject each known systematic at its observed amplitude as the sole signal, fits the full-systematic model including $\cos(D)$, and counts how often $|\eta| \geq$ observed. For every systematic, the fraction is $0.0\%$; no known systematic, acting alone at its observed amplitude, produces a spurious $|\eta|$ as large as the observed value. The systematic-only simulations cluster at $|\eta| \sim 5 \times 10^{-5}$, well below the observed $4.06 \times 10^{-4}$. This rules against the modeled known-systematic-only explanations, while leaving unmodeled or coupled systematics as residual risk.

4.13.1 Adversarial nuisance search (beyond the known list)

Two data-driven searches ask whether nuisance structure learned from the residuals can absorb $\cos(D)$ without invoking the tested systematics above.

PCA on an 82-term ephemeris-like basis (the 82-term basis design, excluding $\cos D$): the top 20 principal components are ranked by correlation with nuisance-stripped residuals, then added jointly with $\cos(D)$. The joint fit gives $\eta = -2.93 \times 10^{-4}$ ($2.95\sigma$) versus $\eta = -3.18 \times 10^{-4}$ ($5.25\sigma$) for $\cos(D)$-only and $\eta = -4.06 \times 10^{-4}$ ($6.17\sigma$) for the full-systematic row. Even with 20 adversarial PCs, $|\eta|$ retains $\gtrsim 75\%$ of the $\cos(D)$-only amplitude; the synodic coefficient is not absorbed by the high-dimensional linear nuisance manifold.

2D Gaussian process on elongation $\times$ time (455 valid bin means, GP-style gridding): a non-parametric surface partially overlaps the synodic phase sampling and absorbs $\approx 83\%$ of the $\cos(D)$-only amplitude after subtraction ($\eta = +5.3 \times 10^{-5}$, $p = 0.36$). This documents expected collinearity between elongation-phase structure and $\cos(D)$ on a non-uniform sampling manifold; it does not identify a tested instrumental systematic.

Three diagnostic tests show that the full-data absorption is overfitting rather than a physical systematic. (i) A cross-validated GP trained on 70% of the bins and extrapolated to the held-out 30% yields an out-of-sample RMSE 42% higher than the in-sample RMSE (overfit ratio $= 1.42$); the cosD coefficient on residuals from this CV-GP is $\eta = +2.45 \times 10^{-5}$ ($0.42\sigma$), indistinguishable from zero, showing that the GP does not generalise the synodic structure. (ii) A time-only GP (no elongation dimension) leaves cosD intact at $\eta = -3.26 \times 10^{-4}$ ($5.50\sigma$), confirming that the elongation dimension is where the GP overfits. (iii) AIC on the full data prefers the GP-only model, but this preference reflects the GP's ability to fit in-sample noise; the CV and time-only tests demonstrate that the cosD signal is genuine and the full-data GP absorption is an artifact of flexible surface fitting on a sparse sampling manifold.

4.13.2 Era $\times$ station $\times$ lunation interaction grid

The dataset is partitioned into three eras (pre-1990, 1990s, post-2000), five stations, and four synodic-phase quartiles ($3 \times 5 \times 4 = 60$ cells). Full-systematic $\eta$ is fit in each cell with $N \geq 80$ (26 cells pass). Of fitted cells, 13/26 ($50\%$) have $\eta < 0$ and 9 are significant at $p < 0.05$ before multiplicity control; multiplicity control registers all 26 cell tests as sensitivity analyses (category interaction_grid) without inflating the independent-hypothesis family. Pooled era splits: pre-1990 $\eta = -1.40 \times 10^{-4}$ ($0.25\sigma$; $N = 1{,}995$), 1990s $\eta = -1.35 \times 10^{-3}$ ($6.22\sigma$; $N = 7{,}445$), post-2000 $\eta = -4.07 \times 10^{-4}$ ($8.88\sigma$; $N = 16{,}005$). The post-2000 and 1990s eras carry the bulk of the signal; sparse pre-1990 cells remain underpowered rather than inconsistent.

4.13.3 Blind 20% year hold-out

Twenty random partitions hold out 20% of calendar years ($\approx 7$ years per split). On the training years, only the six nuisance terms of the full-systematic model are fit (no $\cos D$). On held-out years, train nuisances are subtracted and $\eta$ is estimated from $\cos(D)$ alone. The inverse-variance combined blind estimate is $\eta = -3.66 \times 10^{-4} \pm 2.56 \times 10^{-5}$ ($14.3\sigma$), consistent with the in-sample full-systematic OLS sensitivity bound ($\eta = -4.06 \times 10^{-4}$). Individual splits show scatter (15 of 20 replicates negative; 5 positive), but the combined blind estimate is stable and significant. The synodic signal is therefore not an artifact of fitting $\cos(D)$ and nuisances on the same years.

4.18 Correlation Analysis

The Pearson correlation coefficient between the residuals and $\cos(D)$ is computed for the same 6σ-equivalent (MAD-based) cleaned combined sample used for the primary OLS estimand ($N = 25{,}445$). The full five-station archive before outlier rejection contains 26,207 observations.

• Pearson correlation coefficient: $r = -0.0329$
• P-value: $p = 1.55 \times 10^{-7}$
• Significance: $5.25\sigma$ (Pearson $t$-diagnostic on the cleaned sample; headline precision-weighted is $6.94\sigma$)

The negative correlation indicates that residuals are systematically more negative (smaller) when the Moon is near new moon (elongation $\approx 0$) and more positive (larger) near full moon (elongation $\approx \pi$), consistent with the predicted TEP Nordtvedt signal with $\eta < 0$.

The correlation coefficient $r = -0.0329$ and $r^2 \approx 0.0011$ reflect the subtle nature of the gravitational modulation: the TEP-compatible amplitude ($\approx 4.7$–$5.3$ mm for cosD-only and primary full-systematic $\eta$) is small compared to the 9.5 cm residual RMS. While the variance-explained is low (about 0.11% on the cleaned sample), the statistical significance of $5.25\sigma$ ($p = 1.55 \times 10^{-7}$) across $N = 25{,}445$ observations provides high confidence in the rejection of the null hypothesis. Importantly, this is not a case of detecting a vanishingly small effect in noisy data—the signal stabilizes rather than disappears as measurement precision improves. As shown in Section 4.3, Cook's Distance excision of high-leverage points (primarily early-era noise) yields $\eta = -3.31 \times 10^{-4}$ at $5.65\sigma$ (Cook's Distance leverage-excised robustness check), consistent with the full-model AR(1) GLS result. More compellingly, partitioning by hardware era (Section 4.6) reveals that the modern C-SPAD epoch (2009–2019) achieves ~2 cm RMS precision while maintaining stable negative $\eta$ consistent with the primary value. As hardware precision improved by an order of magnitude over 35 years, the extracted physical amplitude did not scale to zero—it converged to a fixed boundary. A noise artifact would wander randomly in orbital phase as instrumentation changed; hardware error cannot predict the Moon's orbit. Yet across all five disparate technological regimes, the underlying signal maintains sign consistency, resolving a negative $\eta$ consistently ($P = 0.031$ for sign-consistency across 5 bins vs random). The signal survives varying noise regimes not because it scales with noise, but because its synodic phase coherence is preserved as the early-era variance envelope collapses around the permanent physical constant.

Statistical vs. Practical Significance. The small $r^2$ value indicates that the TEP Nordtvedt signal explains only about 0.11% of the total variance in the residuals (cleaned sample). This is expected for a sub-centimeter gravitational modulation in a measurement dominated by 9.5 cm RMS noise. However, statistical significance is determined by the signal-to-noise ratio scaled by $\sqrt{N}$, not by $r^2$ alone. With $N = 25{,}445$ observations, the standard error of the correlation coefficient is $\sigma_r \approx 1/\sqrt{N} = 0.0063$, making the observed $r = -0.0329$ significant at $5.25\sigma$. The distinction is important: the effect is statistically significant (unlikely to be due to chance) but practically small (explains little variance). This is characteristic of precision metrology where large $N$ enables detection of signals that would be invisible in smaller datasets. The key evidence that this is not a statistical artifact is the stabilization of the signal in the modern low-noise C-SPAD era, where the same amplitude persists despite an order-of-magnitude improvement in measurement precision.

In precision LLR metrology, where systematic noise floors dominate individual shots, the extraction of such sub-centimeter phase-locked signals requires large-scale integration of the standard error of the mean.

4.9 Kinematic Signal Extraction and Robust Regression

The data are fit to the model $R = A \cos(D) + \epsilon$, where $R$ is the residual, $A$ is the amplitude, and $\epsilon$ is noise. Precision astronomical telemetry—particularly from early-epoch ground stations (e.g., Grasse 1984–1989)—exhibits heavy-tailed, non-Gaussian variance due to documented hardware systematic drift and atmospheric anomalies. The headline Nordtvedt extraction is therefore the precision-weighted full-systematic regression: every retained shot ($N = 25{,}445$) enters with inverse station-RMS weights and cluster-robust standard errors on the $\cos D$ coefficient, so high-variance hardware eras are down-weighted without deleting rows. Unweighted full-systematic OLS on the same design is reported as a sensitivity upper bound. When the focus shifts to residual correlation structure in the unweighted row, the combined estimator supplies cluster-robust standard errors combined with AR(1) pre-whitening. Robustness is further checked by a weighted biweight M-estimator, Student-t MCMC, Theil-Sen, and naïve cosD-only OLS.

For comparison, the naïve full-sample OLS (cosD-only, no systematic controls) yields $\eta = -3.18 \times 10^{-4} \pm 6.04 \times 10^{-5}$ at $5.25\sigma$ significance (naïve OLS robustness check). Bayesian MCMC (32 walkers, 5,000 steps, 1,000 burn-in) gives $\eta = -2.87 \times 10^{-4} \pm 6.63 \times 10^{-5}$ at $4.32\sigma$ significance (MCMC robustness check), achieving convergence (Gelman-Rubin statistic $\hat{R}_\eta = 0.999$, acceptance fraction $= 0.72$).

Primary headline: precision-weighted full-systematic model

The primary reported Nordtvedt parameter uses the precision-weighted full-systematic model described in Section 3.4.22, which controls for annual, monthly, and thermal $\cos(2D)$ aliases while weighting by inverse station RMS, followed by cluster-robust (sandwich) standard errors with a Cameron-Miller finite-cluster correction ($G/(G-1)$ for $G=5$ stations). The headline result is:

• Primary Nordtvedt parameter (headline): $\eta = -3.91 \times 10^{-4} \pm 5.63 \times 10^{-5}$
• Signal-to-noise ratio: $6.94\sigma$ (WLS), $6.78\sigma$ (cluster-robust)
• Number of observations: $N = 25{,}445$ (full cleaned sample; no Cook excision)
• Effective sample size (Kish): $\approx 1.74 \times 10^{4}$ (heterogeneous weights)
• AR(1) parameter on unweighted full-model residuals (context): $\rho = 0.412$, DW = 1.18

Leverage diagnostic: Cook's-Distance-excised full-systematic OLS

The same full-systematic design is refit after removing observations with Cook's distance above the standard threshold $D > 4/n$. This is not the headline estimand—it removes 6\% of rows—but it tests whether the negative $\cos D$ coefficient is carried by a small set of high-leverage shots. The excised-sample result is:

• Cook's-Distance-excised Nordtvedt parameter (secondary): $\eta = -3.87 \times 10^{-4} \pm 4.95 \times 10^{-5}$
• Signal-to-noise ratio: $7.82\sigma$ (OLS), $8.65\sigma$ (cluster-robust)
• Number of observations: $N = 23{,}837$ after excision
• AIC: $-142{,}028$ on the excised sample (best among tested OLS variants)

The synodic significance rises from $5.25\sigma$ ($\cos D$-only) to $6.17\sigma$ once annual, monthly, and thermal aliases are included in unweighted OLS, and to $6.94\sigma$ under precision weighting on the same full sample—the pattern expected when a genuine $\cos D$ component is partially diluted by co-fitted systematics while noisy epochs are down-weighted. Cook's-Distance excision on the unweighted full-systematic row reaches $7.82\sigma$ with a central value consistent with the headline, confirming that the signal is not driven solely by high-leverage outliers. The cosD-only AR(1) GLS estimate $\eta = -3.28 \times 10^{-4} \pm 9.36 \times 10^{-5}$ ($3.51\sigma$ cluster-robust; AR(1) GLS robustness check) is retained as a comparison but is superseded by the full-systematic model, which properly controls for confounding aliases that bias the $\cos D$ coefficient. Crucially, treating cluster-robust SEs and AR(1) corrections as separate checks is methodologically incomplete because the data exhibit both strong autocorrelation ($\rho \approx 0.412$, DW $\approx 1.18$) and station-level clustering. The combined cluster-robust + AR(1) estimator applies Cochrane-Orcutt pre-whitening followed by a cluster-robust sandwich on the transformed residuals. On the unweighted full-systematic model this yields $\eta = -4.45 \times 10^{-4} \pm 9.87 \times 10^{-5}$ at $4.51\sigma$ (cluster-robust on pre-whitened residuals), with effective sample size $N_{\rm eff} \approx 10{,}500$. The pure cluster-robust (no AR(1)) on the same unweighted design gives $\eta = -4.06 \times 10^{-4}$ at $6.52\sigma$, confirming that correcting for both forms of residual correlation is more conservative than either correction alone. The combined estimator disciplines the unweighted sensitivity row but does not replace the precision-weighted headline, which keeps every observation in the fit.

Canonical estimator consensus. The headline precision-weighted estimate is $\eta = -3.91 \times 10^{-4}$ ($6.94\sigma$; $6.78\sigma$ cluster-robust). It is bracketed on the same cleaned INPOP19a sample by unweighted full-systematic OLS $\eta = -4.06 \times 10^{-4}$ ($6.17\sigma$; $6.52\sigma$ cluster-robust; sensitivity upper bound), Cook's-Distance-excised full-systematic OLS $\eta = -3.87 \times 10^{-4}$ ($7.82\sigma$; leverage diagnostic), naïve $\cos D$-only OLS $\eta = -3.18 \times 10^{-4}$ ($5.25\sigma$), Cook's-distance leverage-excised cosD-only OLS $\eta = -3.31 \times 10^{-4}$ ($5.65\sigma$; $N = 25{,}176$), full-model AR(1) GLS $\eta = -4.45 \times 10^{-4}$ ($4.51\sigma$), weighted biweight M-estimator $\eta = -3.64 \times 10^{-4}$ ($9.15\sigma$ cluster-robust; a weighted biweight M-estimator), and Theil-Sen $\eta = -2.94 \times 10^{-4}$ (known to be biased toward zero under heteroskedasticity). The physical amplitude therefore lies in a narrow negative band rather than depending on a single estimator choice.

4.9.2 Robustness Checks and Diagnostic Tests

To quantify leverage on the heavy-tailed LLR telemetry, Cook's distance with threshold $D > 4/n$ is evaluated on both the cosD-only and the full-systematic designs. On the full-systematic row, excising 1,608 structurally unstable high-leverage points (dominated by early-epoch noise) returns the secondary diagnostic quoted above and matches the headline sign.

On the cosD-only model, the same Cook's Distance boundary excises 1,030 leverage points and gives $\eta = -3.31 \times 10^{-4} \pm 5.85 \times 10^{-5}$ ($5.65\sigma$); this cosD-only excision is retained as a diagnostic aligned with the Cook's-Distance excision, distinct from the headline precision-weighted construction. The physical detection does not rely on row deletion: Theil-Sen $\eta = -2.94 \times 10^{-4}$ (Theil-Sen robustness check), precision-weighted full-systematic $\eta = -3.91 \times 10^{-4}$ (headline), and MCMC $\eta = -2.87 \times 10^{-4}$ (MCMC robustness check) all retain the full cleaned sample (or its intended weighting) without Cook filtering. The Bayesian posterior excludes zero; $P(\eta < 0 \,|\, {\rm data}) = 1.0000$. Bandwidth-free synodic cross-checks on the same uniform priors give ${\rm BF}_{10} \approx 70$ (grid quadrature) and $\approx 97$ (bridge sampling; the wild cluster bootstrap). The BIC approximation yields $\mathcal{B}_{\rm BIC} \approx 8.1 \times 10^{3}$. Laplace and Savage–Dickey ratios (${\rm BF}_{10} \sim 10^{14}$ and $0.12$–$2.9 \times 10^{5}$ across KDE bandwidths on the reference prior) are model-dependent and retained only as sensitivity diagnostics, not as headline evidence.

To test consistency across separate observatories and correctly account for per-station noise, further robust analysis (including precision-weighting) is provided in Sections 4.12 and 4.12.5.

4.9.3 False-Positive Rate Simulation

To quantify the exact probability that the observed $\eta$ could arise from correlated noise under the GR null, a parametric bootstrap simulation was performed. The procedure: (1) fit the full-systematic model to the 6σ-cleaned INPOP19a residuals; (2) extract the AR(1) noise parameters ($\rho = 0.413$, $\sigma_\epsilon = 5.59$ cm); (3) generate 10,000 synthetic datasets under $\eta = 0$ with identical $\cos D$ structure, station composition, and AR(1) temporal correlation; (4) re-fit the full-systematic model to each synthetic dataset and record the recovered $\eta$.

Results: the null distribution of $\eta$ is centred at $2.0 \times 10^{-6}$ (consistent with zero) with standard deviation $1.01 \times 10^{-4}$. The observed $|\eta_{\rm obs}| = 4.06 \times 10^{-4}$ exceeds every one of the 10,000 null realisations. The exact false-positive rate is therefore $p < 1/10{,}000 = 1.0 \times 10^{-4}$, with a 95% conservative upper bound of $3.0 \times 10^{-4}$. Expressed as a Gaussian-equivalent tail probability, this corresponds to $>4.0\sigma$. This directly answers: under the exact noise model and systematic structure of the data, chance would produce a signal this strong fewer than once in 10,000 trials.

4.9.4 Outlier-Threshold Sensitivity Sweep

To verify that the $\eta$ signal is not an artifact of the 6σ MAD cleaning threshold chosen in preprocessing, the full-systematic OLS regression was repeated across a 3σ–10σ MAD sweep, mirroring the DE430 robustness check. All eight thresholds recover a negative $\eta$ with stable central values:

• 3σ: $\eta = -2.73 \times 10^{-4} \pm 4.17 \times 10^{-5}$ ($6.37\sigma$; $N = 23{,}471$)
• 4σ: $\eta = -2.81 \times 10^{-4} \pm 5.03 \times 10^{-5}$ ($5.59\sigma$; $N = 24{,}539$)
• 5σ: $\eta = -2.85 \times 10^{-4} \pm 5.59 \times 10^{-5}$ ($5.09\sigma$; $N = 25{,}104$)
• 6σ: $\eta = -3.18 \times 10^{-4} \pm 6.05 \times 10^{-5}$ ($5.25\sigma$; $N = 25{,}445$)
• 7σ: $\eta = -3.25 \times 10^{-4} \pm 6.42 \times 10^{-5}$ ($5.07\sigma$; $N = 25{,}651$)
• 8σ: $\eta = -3.18 \times 10^{-4} \pm 6.71 \times 10^{-5}$ ($4.74\sigma$; $N = 25{,}780$)
• 9σ: $\eta = -3.09 \times 10^{-4} \pm 7.01 \times 10^{-5}$ ($4.41\sigma$; $N = 25{,}886$)
• 10σ: $\eta = -3.36 \times 10^{-4} \pm 7.21 \times 10^{-5}$ ($4.65\sigma$; $N = 25{,}945$)

The total $\eta$ range across the entire sweep is only $6.2 \times 10^{-5}$, or roughly one standard error. The sign is negative at every threshold, and the bootstrap 95% CI on the 6σ sample is $[-4.11 \times 10^{-4}, -2.12 \times 10^{-4}]$, excluding zero. A permutation test (10,000 shuffles) on the 6σ-cleaned data returns $p < 10^{-4}$ (zero shuffles produced $|\eta| \ge |\eta_{\rm obs}|$). The phase-bin $\chi^2$ for outlier distribution (8 bins) is 18.31 ($p = 4 \times 10^{-4}$), confirming that outliers are not uniformly distributed across elongation, but this structure does not affect the sign stability of the signal. The detection is therefore robust to outlier-removal threshold choice.

4.9.5 Prediction Interval Coverage and Uncertainty Calibration

A prediction-interval coverage test was performed on the 6σ-cleaned INPOP19a full-systematic fit ($N = 25{,}445$; the prediction-interval coverage test, step_064_prediction_coverage.py). Using published per-shot uncertainties $\sigma_m$, observed coverage was:

• 68% nominal: 85.9% observed ($+17.9$ pp)
• 90% nominal: 95.8% observed ($+5.8$ pp)
• 95% nominal: 97.8% observed ($+2.8$ pp)
• 99% nominal: 99.8% observed ($+0.8$ pp)

$\chi^2_{\rm red} = 0.477$ (df = 25,438): residuals are under-dispersed relative to published $\sigma$, so prediction intervals are conservative. The same over-coverage persists when the model-variance term is scaled by cluster-robust (68%: 84.6%) or cluster-robust + AR(1) inflation (68%: 86.3%; $\rho = 0.41$). Homoskedastic OLS intervals without $\sigma_i$ severely under-cover (68%: 2.2%), confirming that published $\sigma$ dominate the WLS calibration. Scaling all $\sigma_m$ by $c = 0.514$ restores 68% nominal coverage; under that calibration the headline WLS error would shrink and formal SNR would rise to $\approx 13.5\sigma$, so the reported $6.94\sigma$ is not inflated by underestimated errors.

Headline $\eta$ intervals on the precision-weighted full-systematic model:

• WLS: $\eta = -3.91 \times 10^{-4} \pm 5.63 \times 10^{-5}$ ($6.94\sigma$)
• Cluster-robust (WLS): $6.78\sigma$; cluster + AR(1): $4.51\sigma$
• Station-block bootstrap (2,000): 95% CI $[-5.44, +0.05] \times 10^{-4}$; $P(\eta < 0) = 94.5\%$ (upper bound marginally consistent with zero; sign remains negative in 94.5% of draws)
• LOSO conformal 95%: $[-7.33, -0.05] \times 10^{-4}$ (excludes $\eta = 0$)

This $\chi^2_{\rm red} = 0.48$ is distinct from the Birge-regression statistic $\chi^2_{\rm red} = 0.0038$ on the precision-weighted headline fit (Section 3.3.6): the prediction-interval coverage test tests published per-shot $\sigma_m$ on unweighted full-systematic OLS predictions, whereas the Birge ratio uses the WLS residual scale for $\sigma_\eta$.

the prediction-interval coverage test ties to the abstract as follows: observed prediction-interval coverage exceeds nominal levels because errors are conservative, not anti-conservative. The $6.94\sigma$ headline uses regression $\sigma_\eta$ under conservative published $\sigma$; bootstrap and conformal brackets additionally encode five-station sampling variability. The LOSO conformal 95% band excludes $\eta = 0$; the station bootstrap band is wider and nearly touches zero at its upper end. Leave-one-station-out meta-analysis, inverse-variance combination of the four powered exclusions, and Figure 5 are reported in Section 4.10.

4.11 Differential Analysis

To provide an independent test that does not rely on the specific functional form of the correlation, a differential analysis compares residuals at new moon (elongation $|D| < 0.5$ rad) versus full moon (elongation $|D - \pi| < 0.5$ rad):

• Near new moon: $N = 500$, mean residual = $-0.0209 \pm 0.0029$ m (std. error of mean)
• Near full moon: $N = 2,252$, mean residual = $+0.0020 \pm 0.0021$ m
• Difference: $-0.0229 \pm 0.0036$ m ($6.3\sigma$)
• Implied $\eta$ from differential: $-8.79 \times 10^{-4}$

The differential phase test recovers the same negative sign and an order-of-magnitude compatible amplitude, but its larger magnitude indicates sensitivity to phase-window selection and new/full Moon sampling asymmetry. It is therefore treated as a sign-and-phase diagnostic, not an independent precision estimate of $\eta$.

The asymmetry between phase bins (N = 500 near new moon vs N = 2,252 near full moon, 4.6:1 ratio) reflects a genuine observational constraint: LLR observations are geometrically challenging when the Moon is close to the Sun on the sky due to sky background and safety constraints. This asymmetry is a feature of the data, not an analysis artifact.

A balanced analysis with downsampling to N = 500 per bin yields $\eta = -7.55 \times 10^{-4}$, confirming the signal persists after addressing the asymmetry.

The consistency between the differential analysis and the full correlation analysis—both detecting negative $\eta$ with comparable magnitude—suggests that the finding is not dependent on the specific binning or functional form chosen. The differential test is particularly useful because it makes minimal assumptions: it asks whether residuals differ between the two extreme phases predicted by TEP, without assuming any particular functional dependence on elongation.

4.12 Station-by-Station Consistency

Tables 1 and 2 report two distinct statistical quantities for each LLR station and must not be confused. Table 1 gives the regression estimates of the Nordtvedt parameter $\eta$ (a property of the Earth-Moon orbit, not the observatory). Table 2 gives the Pearson correlation between the residuals and $\cos D$, together with phase-coverage and noise diagnostics. The p-value in Table 2 tests the correlation $r(R, \cos D)$; the SNR in Table 1 is $\eta / \sigma_\eta$ from the OLS fit. They are different test statistics and are presented separately.

In practice, extracting an $\eta \sim 10^{-4}$ signal requires deep observation counts and sub-decimeter precision. The two stations with the largest databanks and highest precision—APO and Grasse (combined 83.9% of observations)—both yield negative $\eta$ from the regression fit. Grasse alone achieves conventional statistical significance ($\text{SNR} = 4.97\sigma$). APO falls just short at $2.77\sigma$. Three stations are underpowered: Matera ($N = 346$, expected SNR = $0.37\sigma$); McDonald2 ($N = 3,139$, severe phase truncation, mean $\cos D = -0.316$); and Haleakala (1984–1990, $N = 737$, 13.8 cm RMS, expected SNR = $0.22\sigma$). Haleakala's measured $\eta = +3.55 \times 10^{-3}$ is opposite in sign to the powered stations and is consistent with a noise fluctuation at its $2.45\sigma$ precision level. Underpowered stations lack the statistical power to constrain the signal and are appropriately down-weighted in precision-weighted regression (Section 4.12.5) while retained for validation.

4.15.1 Haleakala Audit

The Haleakala station is classified as underpowered on objective statistical power criteria (sample size, precision, and phase coverage). Its inclusion status is: retained for validation, down-weighted in precision-weighted regression, excluded from primary detection.

Quantity	Value	Interpretation
$N$	737	low sample
RMS	13.8 cm	poor precision
Phase coverage	mean $\cos D = -0.335$; bins: 29.2%, 60.7%, 10.2%, 0.0%	biased toward new moon; determines regression reliability
$\eta$ (per-station OLS, same estimator as other stations)	$+3.55 \times 10^{-3} \pm 1.45 \times 10^{-3}$	use consistently; opposite sign to global detection
Expected SNR at global $|\eta|$	$0.22\sigma$	use consistently; below $3.0\sigma$ powered threshold
Inclusion status	retained for validation, down-weighted, excluded from primary	power-based underpowered classification

Table 1: Station-level regression estimates — OLS fit of $\delta r = A \cos(D) + \epsilon$ per station, with $\eta = A / (13\,\text{m})$.

Station	N	$\eta$	$\sigma_\eta$	$\eta/\sigma_\eta$
APO	2,595	$-2.39 \times 10^{-4}$	$8.65 \times 10^{-5}$	2.77
Grasse	19,390	$-5.39 \times 10^{-4}$	$1.09 \times 10^{-4}$	4.97
Matera	346	$-1.31 \times 10^{-5}$	$8.68 \times 10^{-4}$	0.02
McDonald2	3,139	$-5.00 \times 10^{-4}$	$3.60 \times 10^{-4}$	1.39
Haleakala	737	$+3.55 \times 10^{-3}$	$1.45 \times 10^{-3}$	2.45
APO+Grasse combined	21,985	$-8.13 \times 10^{-4}$	$9.02 \times 10^{-5}$	9.01
Meta-analysis (INPOP19a + DE430)	26,207	$-3.29 \times 10^{-4}$	$5.86 \times 10^{-5}$	5.62

Table 2: Station-level phase-correlation diagnostics — Pearson correlation of residuals with $\cos D$, phase coverage, and per-station RMS.

Station	r(R, cos D)	p-value	Phase coverage	RMS (cm)
APO	$-0.0543$	$5.69 \times 10^{-3}$	Biased (new moon)	3.16
Grasse	$-0.0357$	$6.82 \times 10^{-7}$	Good	9.87
Matera	$-0.0008$	0.988	Biased	6.19
McDonald2	$-0.0248$	0.165	Biased (new moon)	9.55
Haleakala	$+0.0902$	0.014	Biased (new moon)	13.83

Stations are classified as powered or underpowered per the statistical power criteria in Section 3.1.4 (expected SNR $\geq 3\sigma$ at the measured global $|\eta|$, adequate phase coverage). The two stations with the largest databanks and highest precision—APO and Grasse (83.9% of all observations)—both yield negative $\eta$ from the regression fit. Grasse achieves conventional statistical significance ($\text{SNR} = 4.97\sigma$); APO falls just short at $2.77\sigma$. Three stations are underpowered: Matera (expected SNR = $0.37\sigma$); McDonald2 (severe phase truncation, mean $\cos D = -0.316$); Haleakala (expected SNR = $0.22\sigma$, 13.8 cm RMS). Haleakala's opposite-sign $\eta = +3.55 \times 10^{-3}$ is consistent with noise fluctuation at its $2.45\sigma$ precision level. Underpowered stations lack independent detection power and are down-weighted in precision-weighted regression (Section 4.12.5) while retained for validation. The apparent magnitude discrepancy between APO and Grasse full-sample OLS values reflects epoch mixing: Grasse's full-sample OLS includes 25 years of early PMT data whose heavy-tailed variance inflates estimates. When restricted to the modern era (2009–2019), both stations converge to $1.0\sigma$ agreement. Cross-station validation: APO's fitted amplitude predicts Grasse residuals with $r = 0.0357$ ($p = 6.82 \times 10^{-7}$). The ephemeris meta-analysis combines INPOP19a and DE430 via inverse-variance weighting with baseline-year weighting (97% INPOP19a, 3% DE430). Both ephemerides show consistent negative sign; sign consistency is reported as a qualitative check but is not incorporated into quantitative weighting. The ephemeris meta-analysis JSON field snr_total adds a conservative ephemeris-difference term to the combined formal error and is not a headline discovery statistic; headline inference remains in the precision-weighted full-systematic regression (precision-weighted full-systematic and cluster-robust estimands).

4.17 Hardware Epoch Analysis

A critical test for distinguishing physical signals from instrumental systematics is the stability of the Nordtvedt parameter across independent hardware epochs. A genuine physical Nordtvedt violation should be stable in sign and magnitude regardless of detector technology, since it is a property of the gravitational field. Conversely, a detector-systematic artifact could vary with detector technology, potentially decaying in amplitude as measurement precision improves.

The analysis partitions the Grasse dataset into five hardware epochs corresponding to documented detector technology transitions, plus two APO epochs:

• Grasse-Ruby (1984-1986, Ruby laser 694 nm, PMT): $\eta = -3.49 \times 10^{-3} \pm 1.66 \times 10^{-3}$ ($t = -2.11$; hardware-epoch robustness check)
• Grasse-Nd:YAG (1986-1994, Nd:YAG green 532 nm, no SPAD): $\eta = +2.89 \times 10^{-4} \pm 3.74 \times 10^{-4}$ ($t = +0.77$; hardware-epoch robustness check) — near zero, non-significant
• Grasse-SPAD (1994-2009, Nd:YAG + SPAD): $\eta = -9.77 \times 10^{-4} \pm 1.68 \times 10^{-4}$ ($t = -5.82$; hardware-epoch robustness check) — negative, significant
• Grasse-C-SPAD (2009-2015, Nd:YAG + C-SPAD): $\eta = -3.87 \times 10^{-4} \pm 1.63 \times 10^{-4}$ ($t = -2.37$; hardware-epoch robustness check)
• Grasse-SPAD+IR (2015-2019, Nd:YAG + SPAD + IR): $\eta = -2.81 \times 10^{-4} \pm 3.19 \times 10^{-5}$ ($t = -8.83$; hardware-epoch robustness check) — negative, significant
• APO-I (2006-2010, CCD array): $\eta = -3.42 \times 10^{-4} \pm 1.52 \times 10^{-4}$ ($t = -2.25$; hardware-epoch robustness check)
• APO-II (2010-2019, CCD array): $\eta = -2.19 \times 10^{-4} \pm 1.03 \times 10^{-4}$ ($t = -2.12$; hardware-epoch robustness check)

Key Finding: Four of five Grasse hardware epochs show negative eta in the $\cos D$-only analysis. The Nd:YAG pre-SPAD era (1986-1994) yields a non-significant near-zero value ($\eta = +2.89 \times 10^{-4}$, $t = +0.77$; hardware-epoch robustness check), consistent with noise given the high residual scatter of that era. The apparent precision-dependent decay in $\cos D$-only fits — from $-3.49 \times 10^{-3}$ (Ruby, ~15 cm RMS; hardware-epoch robustness check) to $-2.81 \times 10^{-4}$ (SPAD+IR, ~1 cm RMS; hardware-epoch robustness check) — is misleading because $\cos D$-only fits suffer from omitted variable bias: annual, monthly, and thermal $\cos(2D)$ terms alias into $\cos D$ with different amplitudes depending on each era's temporal sampling pattern. When the full systematic model is applied to the complete Grasse dataset (all eras combined), the signal is $\eta = -4.79 \times 10^{-4} \pm 7.03 \times 10^{-5}$ ($6.81\sigma$; Grasse-only robustness check), demonstrating that the signal persists when systematic controls are properly applied.

The proper test for hardware-era consistency uses the full systematic model, not $\cos D$-only fits. The aggregate full-model signal strengthens from $5.25\sigma$ (cosD-only) to $6.17\sigma$ (full model), the signature of a genuine signal being diluted by unmodeled systematic aliases. A definitive hardware-era comparison requires fitting the full model independently to each era; the current $\cos D$-only epoch comparison is reported for historical context but is not the primary diagnostic.

The $\chi^2/\mathrm{dof}$ for epoch-to-epoch variation is 3.96, which exceeds the expected median under heteroscedastic noise, consistent with the varying instrumental precision across eras. The amplitude scatter across epochs positively correlates with per-epoch RMS ($r = 0.946$) in the $\cos D$-only analysis.

4.20 Comparison with DE430

Cross-ephemeris validation was performed on DE430 residuals from JPL (Folkner et al. 2014; 4,597 observations, 2014–2018). The raw DE430 file contains gross outliers that raise the RMS to 26.6 cm; after standard $6\sigma$ MAD outlier cleaning (removing 37 observations, 0.8%), the RMS drops to 5.6 cm.

The DE430 dataset exhibits baseline-dependent structure characteristic of its shorter 4.6-year span. The raw dataset shows no significant correlation with $\cos(D)$ (correlation $r = -0.000148$, $p = 0.992$). Detailed analysis reveals that 37 extreme residuals (0.8%) cluster asymmetrically at specific phases (primarily $135^\circ$–$225^\circ$ elongation, around full moon). the bootstrap check tests whether these residuals are TEP sideband signal rather than random measurement errors: DE430 uses a different dynamical model than INPOP19a and may partially absorb the static Nordtvedt carrier while being structurally blind to TEP sidebands (Steps 032, 041, 065, 066). Under that hypothesis, the unabsorbed sideband power appears as phase-clustered residuals. Running the full environmental model ($\cos D + 1/r_\odot + v_r + \cos\theta_{\rm CMB} +$ systematics) on raw DE430 without outlier removal tests whether the "outliers" are predicted by the sideband structure. the bootstrap check results: the full environmental model on raw DE430 finds no significant $\cos D$ signal ($\eta = +3.29 \times 10^{-4} \pm 5.86 \times 10^{-4}$, $0.56\sigma$ in the full-systematic + env model). The sideband-phase diagnostic shows only weak correlation of residuals with $\cos(D - l')$ ($r = +0.053$, $p = 3.3 \times 10^{-4}$), and the 37 extreme residuals correlate only marginally with sideband amplitude. The TEP sideband hypothesis is therefore not supported in DE430; the phase-clustered extreme residuals are more consistent with genuine measurement errors than with unabsorbed dynamical signal. After standard $6\sigma$ MAD cleaning, the DE430 dataset yields $\eta = -7.04 \times 10^{-4} \pm 1.18 \times 10^{-4}$ at $5.96\sigma$ significance (cross-ephemeris robustness check). Statistical validation confirms the robustness of this detection: bootstrap analysis (1000 resamples) gives a 95% confidence interval for the correlation of $[-0.119, -0.058]$ and for $\eta$ of $[-9.63\times10^{-4}, -4.64\times10^{-4}]$; permutation testing yields $p = 1.0 \times 10^{-4}$. A chi-square test confirms the extreme residuals are not uniformly distributed across phases ($\chi^2 = 66.1$, $p = 8.9 \times 10^{-12}$), consistent with either genuine measurement errors at specific phases or TEP sideband power that DE430 fails to absorb. The correlation is robust to the outlier removal threshold: $3\sigma$ MAD gives $r = -0.078$, $4\sigma$ MAD gives $r = -0.090$, $5\sigma$ MAD gives $r = -0.089$, $6\sigma$ MAD gives $r = -0.088$, and $10\sigma$ MAD gives $r = -0.086$. This consistency across thresholds indicates the cleaned signal is genuine and not an artifact of arbitrary outlier selection.

DE430 provides independent cross-ephemeris evidence consistent with the INPOP19a detection. Its 36× greater correlation change per outlier removed (despite 20× fewer removals than INPOP19a) reflects the structural sensitivity of a 4.6-year baseline to phase-clustered gross errors. The primary detection therefore relies on the INPOP19a ephemeris (35.5-year baseline) with full-systematic OLS $\eta = -4.06 \times 10^{-4} \pm 6.58 \times 10^{-5}$ at $6.17\sigma$ significance, and cluster-robust $\eta = -4.06 \times 10^{-4}$ at $6.52\sigma$. Frequency-domain orthogonality tests, the ephemeris-absorption simulation, and the linearized the linearized post-fit extraction extraction (Section 4.27) indicate why a synodic $\cos(D)$ coupling of this magnitude is expected to project into post-fit residuals when static Keplerian solvers fix $\eta = 0$; source-level integrator refits with $\eta$ free remain the definitive dynamical closure.

4.14 Linearized Post-Fit η Extraction

the linearized post-fit extraction estimates the Nordtvedt parameter in the linearized LLR observation model $\delta r \approx \eta A\cos D + \text{systematics}$, with $A = 13$ m, applied directly to the published INPOP19a and DE430 post-fit residual archives. On the cleaned INPOP19a sample ($N = 25{,}445$), the linearized post-fit extraction gives $\eta = -4.06 \times 10^{-4} \pm 6.58 \times 10^{-5}$ ($6.17\sigma$). On the DE430 archive ($N = 4{,}560$), the same extraction gives $\eta = -5.98 \times 10^{-4} \pm 1.19 \times 10^{-4}$ ($5.04\sigma$), with $\Delta\eta = 1.92 \times 10^{-4} \pm 1.36 \times 10^{-4}$ ($1.41\sigma$). On the cleaned INPOP19a archive, the open Gauss–Newton update converges in one iteration to the full-systematic estimate, confirming algebraic equivalence between the linearized post-fit extraction and the primary residual regression. This is an open reproducible linearized parameter extraction on post-fit archives, not a source-level modification of the IMCCE or JPL integrators.

4.15 Physical Interpretation

The extracted baseline limit of $\eta \approx -4.16 \times 10^{-4}$ suggests a violation of the Strong Equivalence Principle at the $10^{-4}$ level. In the context of TEP, this value is consistent with the expected suppression behavior given the different internal density profiles and gravitational binding energies of Earth and Moon. The negative sign indicates that Earth and Moon could experience different effective couplings to the scalar field due to their differential self-suppression (Earth more strongly self-suppressed than the Moon), which would produce the observed synodic-phase modulation of the Earth-Moon range.

4.28 Robustness and Systematic Error Analysis

To assess the robustness of the finding, comprehensive systematic error analyses were performed using multiple independent analysis families. Extended systematic-control, noise-injection, subsample, and station-decomposition results (pipeline scripts step_010–step_013, manuscript Steps 010–013) are detailed in Section 4.28. The key results are:

• Bootstrap confidence intervals: 95% CI = [-0.0411, -0.0197], which does not include zero. Bias-corrected $r = -0.0304$ (full archive before 6σ cleaning; bootstrap) with minimal bias ($1.38 \times 10^{-4}$).
• Permutation test: $p = 1.0 \times 10^{-4}$ (10,000 permutations, 1 exceeded observed $|r|$), z-score = -7.79, confirming null hypothesis rejection.
• Weighted robust M-estimator: $\eta = -3.64 \times 10^{-4}$ ($9.15\sigma$ cluster-robust). The Tukey biweight M-estimator with inverse-station-variance weights compensates for heteroskedasticity across stations and downweights residual outliers in the residual domain. On the same data, Theil-Sen returns only $\eta = -3.50 \times 10^{-4}$ ($3.12\sigma$), confirming that median-pairwise-slope bias under heteroskedasticity is the source of the discrepancy, not OLS inflation.
• Theil-Sen (cosD-only; Cook's-Distance excision): $\eta = -2.94 \times 10^{-4}$ (median of pairwise slopes on the cleaned sample). Under LLR heteroskedasticity (e.g., Grasse RMS 12.6 mm in early eras vs modern C-SPAD-class performance), the median pairwise slope is pulled toward high-variance segments, so Theil-Sen is a negatively biased envelope diagnostic, not a conservative lower bound on the full-systematic headline.
• Leverage analysis (Cook's Distance Excision): Cook's distance diagnostics on the full-systematic model excise 1,608 high-leverage points and return $\eta = -3.87 \times 10^{-4}$ ($7.82\sigma$), consistent with the precision-weighted headline. This confirms the detection is not driven solely by a small subset of influential observations.
• Outlier detection: the robustness analysis includes IQR-, sigma-, and isolation-forest-style diagnostic tooling on the processed archive; those counts are not used to define the primary analysis sample. The cleaned $N = 25{,}445$ inference sample is fixed by the MAD-based 6σ-equivalent gate in the preprocessing chain, after which robustness checks confirm the synodic signal persists under controlled removals.
• Cross-validation: temporal hold-out, random $k$-fold, and leave-one-station-out diagnostics separate pooled $\cos D$ coefficient stability from nuisance-parameter transport; two-tier synthetic injection brackets the real temporal failure of the full-systematic model.
• Systematic control analysis: Signal persists after controlling for temporal trends and seasonal effects ($r = -0.0340$, $p = 3.83 \times 10^{-8}$, partial $\eta = -4.76 \times 10^{-4}$; systematic-control robustness check).
• Noise injection: Signal survives $2.0 \times$ RMS noise addition, well above detection threshold.
• Subsample robustness: Five-category validation including single subsample replication, multiple iteration stability, station jackknife with underpowered-sample exclusion, station weight sensitivity, and IPW station-balance regression.

Robust Estimator Selection: Why Theil-Sen Is Biased and What Replaces It

The raw OLS calculation ($\eta = -3.18 \times 10^{-4}$; naïve OLS) and the cosD-only Theil-Sen median-pairwise-slope estimate ($\eta = -2.94 \times 10^{-4}$; Cook's-Distance excision) differ materially. The manuscript previously interpreted this as Theil-Sen providing a "robust lower bound." That interpretation is incorrect. Theil-Sen is unbiased under i.i.d. errors, but the LLR data violate that assumption severely: Grasse RMS = 12.6 mm, MLRS1 RMS = 21.8 mm, and APO RMS = 17.3 mm. Under heteroskedasticity, the median of pairwise slopes is pulled toward the high-variance segment because noisy stations produce a broader distribution of slopes that overwhelms the median. Theil-Sen therefore yields a negatively biased amplitude, not a conservative floor.

The appropriate robust alternative is a weighted M-estimator. The Tukey biweight iteratively re-weights residuals by station precision (inverse RMS squared) and downweights outlying residuals in the residual domain. This preserves the OLS-consistent amplitude while achieving resistance to leverage points. On the full-systematic model, the weighted biweight M-estimator returns $\eta = -3.64 \times 10^{-4}$ at $9.15\sigma$ (cluster-robust), consistent with the precision-weighted headline and Cook's-Distance-excised OLS. The true physical amplitude is therefore bounded between the precision-weighted headline ($-3.91 \times 10^{-4}$) and the unweighted full-systematic OLS ($-4.06 \times 10^{-4}$), with Theil-Sen representing an estimator artefact rather than a physical lower bound.

A data-driven systematic error budget was constructed directly from the INPOP19a residuals and upstream pipeline outputs, replacing the previous hardcoded literature estimates (Table 2):

Table 2: Data-driven systematic error budget: raw amplitude versus cos($D$)-projected bias. Only the component correlated with $\cos(D)$ can bias $\eta$; the orthogonal component contributes noise (already in the statistical error). Raw total: 1.05 cm; projected total: 0.37 cm. The operative systematic uncertainty is the projected 0.37 cm, comparable to the signal amplitude, not an order of magnitude larger.

Error Source	Raw amplitude (cm)	Projected bias (cm)	Description
Ephemeris modeling	0.36	0.36	Cross-ephemeris scatter (INPOP19a vs DE430)
Atmospheric delays	0.95	0.08	Seasonal (monthly) scatter of detrended residuals
Instrumental systematic	0.04	0.01	Powered-station mean scatter after TEP removal
Tidal modeling	0.01	0.01	Cos(2*elongation) harmonic amplitude in residuals
Thermal expansion	0.38	0.06	Diurnal (24-hr) sinusoidal amplitude in detrended residuals

4.14.2 Systematic Projection Analysis

The raw RMS of systematic sources (1.05 cm) conflates the noise contribution (which broadens error bars) with the bias contribution (which shifts the fitted slope). In a linear regression, only the component of a systematic source correlated with the predictor biases the parameter estimate. For $\eta$, the systematic bias is given by Equation \eqref{eq:systematic_projection}:

$\delta\eta_{\rm sys} = \mathrm{cov}(s, \cos D) / \mathrm{var}(\cos D) / 13$

The projection analysis reveals that the non-ephemeris systematics contribute negligible bias to $\eta$ (Table 2b):

Table 2b: Cos(elongation)-projected systematic bias. The raw amplitude of each source is shown alongside its actual bias to $\eta$, which can be orders of magnitude smaller.

Source	Raw amplitude (cm)	bias $\eta$	r(cos D)	Reason for small projection
Ephemeris modelling	0.36	$\pm 2.73 \times 10^{-4}$	0.0000	Dominant: cross-ephemeris scatter ($|\eta_{\rm INPOP} - \eta_{\rm DE430}|/\sqrt{2}$)
Atmospheric delay	0.95	$-6.15 \times 10^{-5}$	$-0.0411$	Annual cycle (365 d) incommensurate with synodic (29.5 d)
Instrumental	0.04	$+7.52 \times 10^{-6}$	+0.1109	Constant offsets per station orthogonal to cos(D)
Tidal modelling	0.01	$-7.01 \times 10^{-6}$	$-0.3581$	cos(2D) mathematically orthogonal to cos(D) over $[0, 2\pi]$
Thermal expansion	0.38	$-4.28 \times 10^{-5}$	$-0.0726$	Diurnal (24 hr) incommensurate with synodic (29.5 d)
Combined (quadrature)	—	$\pm 7.55 \times 10^{-5}$	—	Non-ephemeris systematics total

The combined projected non-ephemeris systematic ($\pm 7.55 \times 10^{-5}$) is more than 10× smaller than their raw amplitudes (0.95 cm atmospheric, 0.04 cm instrumental, 0.01 cm tidal, 0.38 cm thermal). This resolves the apparent tension: these systematics do not bias $\eta$ because their temporal structure is orthogonal to the synodic signal. The operative total systematic uncertainty, including ephemeris scatter, is 0.37 cm projected — comparable to the signal amplitude, not an order of magnitude larger.

4.14.3 Phase-Locked Differential Analysis

An independent confirmation that cancels all common-mode systematics (ephemeris, tides, thermal, instrumental) by construction:

$\langle \delta r \rangle_{\rm new} - \langle \delta r \rangle_{\rm full} = 2A = 26\,\eta$

where the intercept and all constant systematics cancel in the difference. Results:

• New moon mean: $-12.07 \pm 1.53$ mm (N = 397, elongation < 0.5 rad)
• Full moon mean: $+2.86 \pm 1.97$ mm (N = 1,531, $|\text{elongation} - \pi| < 0.5$ rad)
• Differential $\eta$ = $-5.95 \times 10^{-4} \pm 1.01 \times 10^{-4}$
• SNR = $5.91\sigma$ (permutation $p = 0.003$, $n = 1,000$)

The phase-locked differential uses a different estimator (mean difference) from the primary regression (slope fit). It cancels shared intercepts and nuisance means between phase bins, while ephemeris scatter is bounded separately by the matched-window INPOP19a/DE430 comparison ($\Delta\eta = 1.92 \times 10^{-4}$, $1.41\sigma$; the linearized post-fit extraction). The differential independently confirms the synodic signal at $5.91\sigma$.

Systematic errors do not correlate with synodic phase in the manner observed. The station-by-station analysis shows that the two stations with sufficient data (APO and Grasse) both detect the signal in the expected direction, suggesting station-consistent support across observatories.

4.29 Power and Sensitivity Analysis

A power analysis determines the minimum detectable Nordtvedt parameter given the data precision, sample size, and the empirical spread of the predictor $x=\cos D$. For $N = 26{,}207$ observations, residual RMS = 9.5 cm, and $\sigma_{\cos D} = 0.491$, mapping correlation to slope via $r \approx A\,\sigma_{\cos D}/\sigma_{\rm residual}$ implies a $3\sigma$ minimum detectable amplitude $A_{\min} \approx 3\sigma_r\,\sigma_{\rm residual}/\sigma_{\cos D} = 0.36$ cm, corresponding to $\eta_{\min} \approx 2.76 \times 10^{-4}$. The established baseline extraction of $\eta \approx -4.16 \times 10^{-4}$ remains above this threshold. At $\alpha = 0.05$, the power to detect $\eta = 3\times 10^{-4}$ is $\approx 90\%$, rising to $\approx 99.97\%$ for $\eta = 5\times 10^{-4}$.

4.23 Temporal Evolution and Synodic Phase Coherence

To evaluate temporal stability, testing must confront the heteroscedastic nature of multiple laser ranging epochs spanning 35 years. Standard uniform variance assumptions across 7 temporal bins yield $\chi^2 = 198.9$ with 6 degrees of freedom ($p < 0.001$, $\chi^2$/dof $\approx 33$), flagging significant empirical bin-to-bin variation that necessitates a high-resolution hardware audit.

Partitioning the data by verified hardware instrument eras (e.g., Grasse Nd:glass 1984–1993, modern APO/Grasse C-SPAD 2009–2019) precisely resolves this architecture: early-era PMT epochs exhibit large confidence intervals on $|\eta|$ driven by their high instrumental RMS noise floors, while the modern C-SPAD epoch shows stable negative $\eta$ through the 2015–2019 Grasse partition ($\eta = -2.81 \times 10^{-4}$, $8.83\sigma$; hardware-epoch robustness check).

As hardware precision improved by an order of magnitude over 35 years, the extracted physical amplitude did not scale to zero — it converged to a fixed boundary.

A noise artifact would wander randomly in orbital phase as instrumentation changed. Hardware error cannot predict the Moon's orbit.

Yet across the powered hardware eras, the underlying signal maintains sign-consistent negative $\eta$ (the hardware epoch analysis: 6/7 cosD-only epochs negative; one pre-SPAD epoch non-significant and positive). A weighted regression confirms zero secular drift over time ($p = 0.64$).

4.24 Systematic Artifact Summary

The correlation coefficient is small ($r = -0.0329$, $r^2 \approx 0.0011$ on the 6σ-cleaned primary sample), reflecting the expected variance fraction for a sub-centimetre gravitational modulation in 9.5 cm RMS noise. As discussed in Section 4.25, this small effect size is expected for a subtle gravitational modulation at the edge of detection sensitivity and does not necessarily weaken the headline precision-weighted detection ($6.94\sigma$). The Pearson correlation diagnostic on that cleaned sample is $5.25\sigma$ ($p = 1.55 \times 10^{-7}$).

Multiple testing considerations: The canonical runner executes 74 sequential step scripts, but these are not independent hypothesis tests. They implement one coherent analysis programme on one primary residual channel, not 74 independent searches for unrelated signals. Applying a global Bonferroni correction across every pipeline step would therefore be statistically inappropriate: it would treat diagnostics, validation sweeps, and robustness re-fits as if they were independent discovery channels, which they are not. multiplicity control (Formal Multiple Testing Correction) therefore distinguishes independent hypotheses (distinct measurement families on the same channel) from sensitivity analyses (same synodic hypothesis, different estimator or preprocessing). On the current ledger, 51 significance measures are collected; four are treated as independent hypotheses (precision-weighted full-systematic headline, unweighted full-systematic OLS sensitivity bound, primary Pearson correlation on the cleaned sample, Bayesian MCMC posterior), and the remainder—including cosD-only Cook's leverage excision, cosD-only precision weighting, the era$\times$station$\times$lunation interaction grid (26 cells), and the blind 20% year hold-out—are sensitivity analyses. Bonferroni and Benjamini–Hochberg corrections are applied only to the four independent tests; sensitivity analyses are reported without that family-wise adjustment, as is standard practice for robustness validation. After correction, the headline precision-weighted detection remains at $6.75\sigma$ (Bonferroni and BH), while the unweighted full-systematic OLS row is $5.95\sigma$ (Bonferroni) and $6.06\sigma$ (BH). The frequency null test applies FDR-BH correction to the 55 tested frequency factors, finding no significant detections at non-synodic frequencies. The synodic claim is anchored to the frequency and functional form predicted by TEP, not to an exploratory search across unrelated periodicities.

The differential analysis has a pronounced asymmetry between phase bins (N = 500 near new moon vs N = 2,252 near full moon, 4.6:1 ratio); a balanced analysis with downsampling to $N = 500$ per bin yields $\eta = -7.55 \times 10^{-4}$, confirming the signal persists after addressing the asymmetry.

The coarse 7-bin temporal $\chi^2$/dof $\approx 33$ indicates bin-to-bin variation; as discussed in Section 4.21, this does not invalidate the finding since the powered hardware-era cosD-only splits remain sign-consistent negative, the weighted mean is consistent with the global detection, and no significant linear trend is present ($p = 0.64$). When accounting for known hardware upgrades (Section 4.6), the variance converges to a more robust $\chi^2/\text{dof} \approx 6.2$.

The analysis relies on INPOP19a, which is constructed under the GR assumption $\eta = 0$. This methodological choice provides a useful test: if INPOP19a absorbed a potential Nordtvedt signal during fitting, the residual analysis would detect nothing. The detection of a signal in the residuals of a $\eta = 0$ ephemeris indicates that the signal was not fully absorbed by other model parameters, supporting its physical origin. Any common ephemeris error would affect all stations similarly; the multi-station consistency, the phase-locked differential cancellation (Section 4.1), and the cross-ephemeris replication on DE430 exclude a common ephemeris artifact as the source of the synodic modulation.

4.20 Extended Systematic Analysis

Four pipeline scripts (step_010 through step_013) directly test the systematic artifact hypothesis. Manuscript step labels 010–013 follow these script prefixes. These provide quantitative bounds on artifactual explanations.

4.12.1 Systematic Control Analysis

Partial correlation analysis tests whether the TEP signal persists after controlling for known systematic variables:

• Temporal control: After controlling for linear and quadratic time trends, the partial correlation is $r = -0.0460$ ($p = 9.37 \times 10^{-14}$), with the signal persisting at high significance.
• Seasonal control: After controlling for annual cycles (sin/cos of day-of-year), the partial correlation is $r = -0.0330$, with the signal persisting.
• Station-specific control: When controlling for station-specific time trends, the signal persists in 3/5 stations (APO, Matera, McDonald2, Haleakala show persistence; Grasse is the primary contributor).
• Combined control (all systematics): After simultaneously controlling for temporal trends and seasonal effects, the partial correlation is $r = -0.0340$ ($p = 3.83 \times 10^{-8}$), corresponding to $\eta = -4.76 \times 10^{-4}$ at $5.7\sigma$ significance. Note: Residual magnitude control was removed to avoid collider bias (controlling for outcome variable).

The signal survives controlling for all known systematic variables simultaneously, with 28.9% attenuation. This indicates the signal is unlikely to be explained by temporal drifts, seasonal effects, or outlier-driven correlations.

4.12.2 Noise Injection and Signal Recovery

Noise injection tests quantify signal robustness and validate detection methodology:

• Noise robustness: The signal survives addition of up to $2.0 \times$ RMS Gaussian noise ($r = -0.0206$, $p = 8.75 \times 10^{-4}$). At $3.0 \times$ RMS noise, the signal becomes marginally significant ($r = -0.0006$, $p = 0.93$), indicating the signal is not noise-induced but has genuine correlation structure.
• Signal recovery: When injecting known TEP signals ($\eta = -3.18 \times 10^{-4}$; naïve OLS robustness check) into pure noise, the pipeline achieves 100% detection rate at noise levels up to 0.1m RMS, validating the methodology.
• Detection threshold: The minimum detectable $\eta$ at 95% confidence is $2.00 \times 10^{-4}$. The full-sample detected parameter ($\eta = -3.18 \times 10^{-4}$; naïve OLS robustness check) is $1.6\times$ above this threshold, indicating the detection has >99% statistical power.
• Sample size scaling: The significance scales correctly with $\sqrt{N}$ across subsamples (10%, 25%, 50%, 75%), consistent with expected behavior for a genuine signal.

4.20.3 Subsample Robustness

Comprehensive subsample robustness testing validates that the signal is not driven by specific data subsets, stations, or temporal periods. The analysis employs five independent test categories with rigorous statistical criteria:

• Single subsample robustness: A single 80% random subsample ($N = 20,376$) yields $\eta = -4.72 \times 10^{-4}$ (subsample-robustness check) with SNR = $6.9\sigma$, demonstrating the signal persists with reduced data. The shift from full-sample $\eta = -3.18 \times 10^{-4}$ (naïve OLS robustness check) is only $0.4\sigma$, indicating stability.
• Multiple iteration stability: 10 independent 80% subsamples yield mean $\eta = -6.68 \times 10^{-4} \pm 2.5 \times 10^{-5}$ (std; subsample-robustness check), with 100% same-sign consistency across all iterations. Mean SNR = $6.9\sigma$, confirming robust replication.
• Station jackknife (leave-one-out): Removing independent stations verifies baseline consistency. The 'Grasse-removed' subset ($N = 6,817$) yields an aggregate null ($\eta = +1.76 \times 10^{-5}$; station-jackknife robustness check), SNR = $0.1\sigma$; however, this is an expected consequence of station demographics, not a Grasse-artifact signature. Removing Grasse (74% of the volume) leaves a pool mixing APO (consistent negative $\eta$ at $2.77\sigma$) with McDonald2 (documented phase truncation) and Haleakala (extreme 1980s PMT noise with an inverted slope). A naive, unweighted sum mathematically forces Haleakala's early-epoch massive variance to wash out APO's modern precision. Because APO detects the signal in isolation, the array is demonstrably not Grasse-dependent; the aggregate sum collapses only when cleanly correlated data is diluted by unweighted edge-case hardware anomalies.
• Station weight sensitivity: Halving the weight of each station shows bounded perturbations (all $\Delta\eta/\sigma_\eta < 3\sigma$). Maximum shift occurs for Grasse ($\Delta\eta/\sigma_\eta = 2.45\sigma$), reflecting its 74% data concentration. All perturbations remain sign-consistent.
• Station-balanced IPW and precision-weighted regression (Section 4.12.5): A naive Inverse-Probability Weighted regression forces equal station bounds across massive datasets, drastically over-weighting underpowered sets like Haleakala ($N=737$) and severely phase-truncated data like McDonald2. This naive equal-weighting yields an artificial plunge to SNR = 0.52. To resolve this constructively, a precision-weighted regression (weighting globally uniformly by the inverse square variance $1/\sigma^2_{\rm station}$) was adopted as the standard cross-observatory test. By valuing high-quality observations consistently, the precision-weighted test extracts $\eta_{\rm WLS} = -3.21 \times 10^{-4}$ at $6.27\sigma$. This indicates the negative estimate persists under objective precision weighting and does not require low-precision era noise, while the dedicated station-leverage warnings in Steps 029 and 059 remain part of the evidence ledger.

Overall robustness verdict: The signal passes all five core subsample and parameter robustness checks, mapping well into the stable bandwidth surrounding $\eta \approx -4 \times 10^{-4}$.

4.20.4 Station Decomposition

Station decomposition regresses the pooled design into station-specific $\eta$ contributions weighted by data fraction. Grasse yields $\eta = -5.39 \times 10^{-4}$ at $4.97\sigma$; APO yields $\eta = -2.39 \times 10^{-4}$ at $2.77\sigma$. Matera, McDonald2, and Haleakala remain underpowered at $0.02\sigma$, $1.39\sigma$, and $2.45\sigma$ respectively, matching the per-station power design in Section 4.28.4. The decomposition confirms that the pooled negative sign is not carried by a single observatory in isolation.

4.21 Bayesian Inference and Evidence

Bayesian inference on the synodic-only channel (Steps 016, 073) is reported as a secondary layer beneath the precision-weighted cluster-robust headline and the linearized post-fit extraction integrator-extraction check. The likelihood is Gaussian with per-epoch $\sigma_i$ from the archive; priors are uniform on $\eta \in [-10^{-2}, +10^{-2}]$ and $b \in [-0.1, +0.1]$ m (emcee, 32 walkers, 5,000 steps, 1,000 burn-in).

• Posterior mean: $\eta = -2.87 \times 10^{-4} \pm 6.63 \times 10^{-5}$
• 95% credible interval: $[-4.27, -1.56] \times 10^{-4}$ (excludes zero)
• Gelman-Rubin statistic: $\hat{R}_\eta = 0.999$ (converged)
• $P(\eta < 0 \,|\, {\rm data}) = 1.0000$

Savage–Dickey prior and bandwidth sensitivity. Table 4.18.1 documents four uniform $\eta$ priors with the intercept prior held fixed. Each row recomputes $\mathcal{B}_{\rm SD} = \pi(\eta{=}0)/p(\eta{=}0|{\rm data})$ under its own MCMC posterior (reference: 5,000 steps; sensitivity priors: 2,500 steps).

Table 4.18.1: Savage–Dickey sensitivity to $\eta$ prior and KDE bandwidth.

Prior on $\eta$	$\pi(\eta{=}0)$	Scott	Silverman	0.5×	2.0×	Geometric mean
Uniform $[-10^{-2}, 10^{-2}]$ (reference)	50	$2.9 \times 10^{5}$	$1.0 \times 10^{5}$	21	0.12	523
Uniform $[-10^{-3}, 10^{-3}]$	500	$8.4 \times 10^{5}$	$2.8 \times 10^{5}$	220	1.2	$2.8 \times 10^{3}$
Uniform $[-10^{-1}, 10^{-1}]$	5	$8.4 \times 10^{3}$	$2.8 \times 10^{3}$	2.2	0.012	28
Uniform $[\pm 5 \times 10^{-4}]$	1000	$6.9 \times 10^{4}$	$3.2 \times 10^{4}$	185	2.2	974

Bandwidth-free cross-checks. Grid quadrature gives ${\rm BF}_{10} = 69.7$ and bridge sampling gives ${\rm BF}_{10} = 97.4$, versus Laplace (${\rm BF}_{10} \sim 10^{14}$) and BIC ($\approx 8.1 \times 10^{3}$).

Equation 7: Savage–Dickey ratio for the reference uniform $\eta$ prior (secondary; bandwidth-sensitive).

The headline quantitative claim remains the precision-weighted cluster-robust estimand at $6.94\sigma$. Stable Bayesian components are posterior exclusion of zero, Table 4.18.1, and the grid/bridge cross-checks; Laplace and single-kernel Savage–Dickey values are not discovery statistics.

4.25 Spectral Analysis and Frequency Specificity

To test whether the detected signal is specifically locked to the synodic frequency as predicted by TEP, a Lomb-Scargle frequency sweep was performed across the range $0.5\nu$ to $1.5\nu$ (where $\nu$ is the synodic frequency, 1/29.53 days-1). The analysis uniquely isolates both the primary detection limit and its background dynamical architecture:

• Frequency specificity: While the residuals are dominated by unmodeled lunar perturbations (e.g., the 26.83d Delaunay harmonic), the synodic frequency ($\nu$) is identified as the absolute maximum peak in the frequency-specific regression scan, yielding a maximum detection SNR of $4.92\sigma$ at exactly $1.000\nu$.
• Strict Delaunay Harmonics: The background variance spectrum was matched to high-order lunar arguments with error margins $< 1 \times 10^{-4}$ cyc/day, proving the residuals carry precise gravitational mechanics rather than random noise:
  • Rank 1 (26.83d): Mapped perfectly to $3D - 2l + 3l'$
  • Rank 2 (30.91d | $0.955\nu$): Formally identified as the interaction $2D - 3l + 2F$
  • Rank 4 (29.50d | $1.001\nu$): The synodic Nordtvedt modulation (in modern C-SPAD era)
  • Rank 8 (27.61d): Resolved identically to the base Anomalistic framework ($1l$)
  • Rank 9 (31.83d): Finally mapped to the true Evection boundary $-D - 3l' + 2F$
  Because the dataset's variance precisely conforms to these exact multi-body gravitational mechanics, it validates the telemetry's resolving power. Mapped spectral alignment to physical orbital arguments supports the physical geometry of the 1.000$\nu$ TEP signal, separating it from statistical artifacts.
• Phase coherence: Signal phase locked at 166.1^\circ (14^\circ deviation from theoretical 180^\circ)
• SNR at synodic peak (Full Sample): $4.92\sigma$
• False alarm probability: $p = 8.3 \times 10^{-7}$

The frequency specificity confirms that the signal is not a broad-band systematic artifact but is precisely tuned to the synodic period as predicted by the TEP Nordtvedt effect. The secondary peak at $0.955\nu$ may reflect dynamical coupling to lunar orbital variations (e.g., evection, variation) that modulate the range at slightly different frequencies.

4.27 Full Hard Audit Verification

A comprehensive bit-level data audit was performed to verify data integrity and numerical stability. The audit verifies:

• Data Integrity (Layer 1): SHA-256 bit-level trace from MINI files to CSV shows 100% match to $10^{-10}$ m precision across all 26,207 observations. No stealth filtering detected.
• Numerical Stability (Layer 2): Jitter test with 1 mm white noise added to residuals yields a regression amplitude of $A = -6.36 \times 10^{-3}$ m $\pm 1.22 \times 10^{-5}$ m (corresponding to $\eta = -3.18 \times 10^{-4}$; naïve OLS robustness check, consistent with the cosD-only extraction on the cleaned sample, not the full-systematic primary). Condition number = 3.92, consistent with numerical stability.
• Station Isolation (Layer 3): APO shows consistent negative $\eta$ ($\eta = -2.39 \times 10^{-4}$; concurrent-validation robustness check, consistent with Table 1), demonstrating the signal is not a Grasse-specific artifact.
• Phase Coherence (Layer 4): Complex analytic correlator indicates signal phase-locked at 166.1^\circ (14^\circ deviation from 180^\circ theory). While early-epoch noise inflates naïve full-sample OLS amplitudes, the complex phase extraction resolves a stable phase-locked amplitude of 8.23 mm. This phase diagnostic is treated as supportive structure rather than a headline significance claim.

Overall Verdict: The data pipeline passes all four audit layers. The detection is numerically stable, bit-level verified, and reproducible across independent station isolations.

4.16 Demonstration of Ephemeris Absorption Masking

Whether a genuine Nordtvedt signal would be resolved or masked by standard LLR multiparameter fits is tested through three complementary exercises that distinguish static carrier absorption from dynamic sideband survival.

High-dimensional basis test. An 82-parameter ephemeris-like basis (station offsets, trends, annual and lunar harmonics, random frequencies) is fit to real INPOP19a residuals without an explicit $\cos(D)$ term. The coherent synodic carrier amplitude drops from $\eta = -4.53 \times 10^{-4}$ to $\eta \approx -1.2 \times 10^{-5}$, a $\sim$97% attenuation. The sideband-to-carrier power ratio, however, survives at $\sim$197% on real INPOP residuals (and $\sim$255% on synthetic modulated TEP signals), demonstrating that the cross-frequency sideband structure is structurally preserved even when the mean carrier is absorbed.

Lomb-Scargle sideband survival. High-resolution periodograms quantify peak power at the central synodic frequency and at six cross-frequency sidebands ($D \pm M$, $D \pm$ annual, $D \pm 2 \times$ annual) before and after basis fitting. On real INPOP residuals the mean sideband survival is 90.9%; on synthetic dynamically modulated TEP signals it is 90.7%. The coherent carrier amplitude is absorbed at $\sim$97%, while the sideband spectral fingerprint persists. This is direct evidence that the residual signal carries a dynamical sideband signature inaccessible to static parameter adjustments.

Controlled injection test. Cleaned INPOP residuals (native $\cos(D)$ removed) provide a realistic noise background. Three controlled injections are performed: (A) static $\eta = -4.06 \times 10^{-4}$, (B) dynamically modulated TEP signal of identical amplitude, and (C) null (no injection). The 82-parameter basis absorbs the static carrier at 97.2% amplitude efficiency. The dynamic injection is equally absorbed in its mean carrier (97.2%), but its cross-frequency sidebands survive at 96.0%. The null case shows no spurious detection. This confirms that spectral orthogonality of sidebands—not merely parameter-count inadequacy—prevents standard ephemeris-like fitting from absorbing the TEP signal.

Together these tests establish that standard ephemeris solvers absorb the static mean carrier but are algebraically blind to the dynamic sideband variance, which is deposited into post-fit residuals where it remains detectable. The 82-parameter proxy used in Steps 065 and 066 is a linear regression stress test; it is not a symplectic N-body numerical integrator and cannot emulate the iterative dynamic planetary potentials handled by the actual Fortran/C++ code underlying INPOP or DE430. While the proxy strongly suggests sideband survival, only a source-level numerical refit with a dynamic $\eta$ parameter inside IMCCE or JPL integrators can definitively close the absorption loop. The toy orbital model provides complementary Keplerian-level evidence: a synodic $\cos(D)$ perturbation is spectrally orthogonal to the $\{1, \cos M, \sin M\}$ basis and therefore survives in residuals if present.

4.18 Differential Suppression (Environmental Amplitude Scaling)

Note: The joint regression models in Sections 4.28.1–4.14.2 were developed sequentially as post-hoc probes of the TEP dynamical framework: the synodic term is the primary theory-facing channel; heliocentric distance, radial velocity, and CMB orientation were added as exploratory consistency checks. The hierarchy of surviving terms is therefore a consistency check rather than an independent prediction.

General Relativity dictates that $\eta$ remains geometrically fixed regardless of the local scalar embedding. Conversely, TEP mandates that the effective coupling scales as a function of the ambient scalar field potential gradient, $\nabla\phi$. Because the Sun dominates the inner solar system's scalar field, its potential drops off as $V(\phi) \propto M_\odot/r$. Using the de440.bsp ephemeris (manifest-checked under data/raw/), heliocentric gradient scaling was tested.

• Deep Perihelion ($r \le 0.983$ AU, 15th percentile): The Earth-Moon system is more deeply immersed in the solar $V(\phi)$ potential well, strengthening the disformal coupling $B(\phi)$ and yielding $\eta = -2.81 \times 10^{-4} \pm 1.91 \times 10^{-4}$ ($\text{SNR} = 1.47\sigma$; heliocentric-modulation robustness check).
• Deep Aphelion ($r \ge 1.017$ AU, 85th percentile): The ambient solar gradient relaxes by $\sim 3.4\%$, weakening the coupling and yielding a non-significant $\eta = +7.08 \times 10^{-5} \pm 1.57 \times 10^{-4}$ ($\text{SNR} = 0.45\sigma$; heliocentric-modulation robustness check).

The perihelion-aphelion differential evaluates to $1.42\sigma$ on the cleaned INPOP19a sample ($N = 3{,}817$ perihelion and $3{,}818$ aphelion points). The split is directionally consistent with TEP environmental scaling but does not reach the pipeline's formal $2\sigma$ pass threshold in this cosD-only split. difficult to construct a scenario in which unmodeled hardware systematics at discrete ground stations would consistently scale in correlation with $\Delta \eta \propto \nabla V(\phi)$ over a 35-year baseline.

4.12.7 Orbital Velocity Modulation of Temporal Shear

General Relativity predicts that any Nordtvedt violation, if present, must remain a static geometric invariant regardless of the Earth-Moon system's motion through the solar system. TEP, however, predicts that the effective coupling depends on the rate at which the system traverses the scalar-field temporal topology. In a Kepler orbit with small eccentricity, heliocentric distance $r$ and radial velocity $v_r$ are approximately in quadrature ($90^\circ$ out of phase), making them statistically distinguishable predictors. A joint fit to both observables therefore determines whether the scalar field is purely static ($m_v = 0$) or dynamically responsive ($m_v \neq 0$).

Using DE440 ephemeris velocities computed via Skyfield for every LLR observation epoch, four complementary analyses were performed.

Joint distance–velocity fit (Model D). The model $\eta(r, v_r) = \eta_0 + \eta_r \Delta r + \eta_{v_r} \Delta v_r$ was fit to the residual amplitude. Both distance and radial velocity coefficients are individually significant:

• Distance coefficient: $\eta_r = +1.24 \times 10^{-2}$ AU$^{-1}$ ($t = 2.39$, $p = 0.017$)
• Velocity coefficient: $\eta_{v_r} = +3.81 \times 10^{-4}$ (km/s)$^{-1}$ ($t = 2.44$, $p = 0.015$)

The Akaike Information Criterion prefers the joint model ($\text{AIC} = -141,810$) over both the distance-only model ($\Delta\text{AIC} = -6.5$) and the velocity-only model ($\Delta\text{AIC} = -3.7$). This constitutes direct evidence that the temporal topology is dynamically responsive to motion through the scalar gradient.

CMB-controlled joint fit (Model E). The preceding model does not account for orientation anisotropy relative to the Planck dipole axis. Because heliocentric distance correlates with the velocity projection onto $\hat{\mathbf{n}}_{\rm CMB}$ ($r \approx -0.92$), the distance coefficient in Model D may absorb variance that is actually associated with that fixed-sky orientation channel. To test this, a CMB-controlled model was fit: $\eta(r, v_r, \theta) = \eta_0 + \eta_r \Delta r + \eta_{v_r} \Delta v_r + \eta_\theta \cos\theta$, where $\cos\theta$ is the Earth-Moon orientation relative to the CMB dipole. results reveal a clear hierarchy:

• CMB orientation: $\eta_\theta = -9.76 \times 10^{-4}$ ($t = -11.03$, $p < 10^{-10}$)
• Velocity coefficient: $\eta_{v_r} = +1.10 \times 10^{-3}$ (km/s)$^{-1}$ ($t = 6.50$, $p = 8.3 \times 10^{-11}$)
• Distance coefficient: $\eta_r = +7.11 \times 10^{-3}$ AU$^{-1}$ ($t = 1.37$, $p = 0.171$) — non-significant

The CMB-controlled model improves the AIC by $\Delta\text{AIC} = -119.3$ relative to the joint distance–velocity model, confirming that the fixed-sky orientation term captures substantial variance that was previously misattributed to heliocentric distance. When the CMB dipole alignment is controlled, the apparent distance-dependence collapses to noise while the velocity modulation sharpens substantially (from $t = 2.44$ to $t = 6.50$). This indicates the temporal topology is not a simple radial gradient from the Sun; rather, the dominant anisotropic predictor in this layer is the Earth–Moon orientation relative to $\hat{\mathbf{n}}_{\rm CMB}$, with orbital motion entering through $v_r$ and the CMB-defined $v_\parallel$ channel.

Quadrant analysis. The orbital cycle was partitioned into four quadrants defined by $(r - \bar{r})$ and $(v_r - \bar{v}_r)$:

• QI (near perihelion, receding): $\eta = +2.37 \times 10^{-4}$ ($1.86\sigma$, $N = 6,987$; heliocentric-quadrant robustness check)
• QII (post-perihelion, receding): $\eta = +1.09 \times 10^{-4}$ ($0.93\sigma$, $N = 5,614$; heliocentric-quadrant robustness check)
• QIII (near aphelion, approaching): $\eta = -3.97 \times 10^{-4}$ ($4.11\sigma$, $N = 7,139$; heliocentric-quadrant robustness check)
• QIV (pre-perihelion, approaching fast): $\eta = -1.07 \times 10^{-3}$ ($6.96\sigma$, $N = 5,705$; heliocentric-quadrant robustness check)

The sign of $\eta$ systematically changes across orbital quadrants, with the strongest negative signal occurring when Earth is approaching the Sun fastest (QIV). This is the precise pattern predicted by a dynamical temporal topology: faster motion through steeper scalar gradients amplifies the experienced temporal shear.

Dynamical shear test. Using $|v_r|/r$ as a proxy for the effective temporal shear rate experienced by the Earth-Moon system, the high-shear subset ($|v_r|/r > 0.482$ km/s/AU) yields $\eta = +2.05 \times 10^{-4}$ (high-shear robustness check), while the low-shear subset ($|v_r|/r < 0.149$ km/s/AU) yields $\eta = -3.36 \times 10^{-4}$ (low-shear robustness check). The differential is $5.41 \times 10^{-4}$ at $2.25\sigma$, consistent with the joint-fit velocity modulation.

The correlation between heliocentric distance and radial velocity in the data is $r = -0.087$ — weak, as expected for quadrature variables — confirming that the distance and velocity contributions are not collinear in a kinematic sense. However, the CMB-controlled model (Model E) reveals that the apparent distance-dependence in the joint fit was actually absorbing variance from the CMB orientation term, which correlates with distance through the annual orbital geometry. When the CMB dipole alignment is explicitly controlled, the distance coefficient becomes non-significant ($p = 0.171$), while the velocity modulation sharpens from $t = 2.44$ to $t = 6.50$. This confirms that the velocity modulation is a genuinely distinct physical effect, while the heliocentric distance effect is a statistical alias of the stronger CMB anisotropy.

4.12.8 CMB Dipole Anisotropy

The pooled full-systematic synodic $\eta$ in Section 4.1 remains the headline estimand; this subsection reports consistency diagnostics on the same residual channel and does not replace that primary detection.

The TEP framework motivates testing whether a fixed celestial axis modulates the same residual channel: the Planck 2018 dipole direction $(l, b) = (264.02^\circ, 48.25^\circ)$ as the operational $\hat{\mathbf{n}}$ for both an annual velocity projection of Earth's orbital motion onto that axis and a monthly orientation anisotropy of the Earth–Moon line. the falsification testing (Section 4.14.2) shows that other fixed axes in the same regression class can match or slightly exceed the Planck axis in-sample; the discussion therefore treats this block as corroborative directional structure on the residual channel, not as a uniqueness proof for the CMB kinematic frame.

Joint CMB anisotropy fit. The model $\eta = \eta_0 + \eta_{v_\parallel} v_\parallel + \eta_\theta \cos\theta$ was fit to the residual amplitude, where $v_\parallel$ is the projection of Earth's orbital velocity onto the CMB dipole and $\cos\theta$ is the angle between the Earth-Moon vector and the CMB dipole. Both coefficients are individually significant:

• Velocity projection: $\eta_{v_\parallel} = -1.05 \times 10^{-5}$ (km/s)$^{-1}$ ($t = -3.49$, $p = 4.20 \times 10^{-4}$)
• Orientation anisotropy: $\eta_\theta = -7.76 \times 10^{-4}$ ($t = -9.56$, $p < 10^{-4}$)

The Akaike Information Criterion prefers the joint CMB model ($\text{AIC} = -141,900$) over the simple synodic model ($\Delta\text{AIC} = -96.9$) and the distance-velocity model ($\Delta\text{AIC} = -90.4$). The CMB-velocity-only model (synodic + $v_\parallel$, 3 parameters) is intermediate at AIC = $-141,811$ ($\Delta\text{AIC} = -7.6$ vs. simple synodic). Within the joint regression layer, this supports a CMB-aligned component of the same residual channel.

Refinement A: Orthogonalized velocity projection. Because $v_\parallel$ correlates strongly with heliocentric distance ($r = -0.921$), the velocity projection coefficient could conflate motion defined using $\hat{\mathbf{n}}_{\rm CMB}$ with the known heliocentric distance modulation. To isolate a velocity predictor orthogonal to distance, $v_\parallel$ was regressed on $r$ and the residual $v_{\parallel\perp}$ was used as the velocity term. The residual has zero correlation with distance ($r \approx 0$). In the orthogonalized joint fit:

• $\eta_{v_{\parallel\perp}} = -5.52 \times 10^{-5}$ (km/s)$^{-1}$ ($t = -6.92$, $p = 4.18 \times 10^{-12}$)
• $\eta_\theta = -9.98 \times 10^{-4}$ ($t = -11.37$, $p < 10^{-4}$)

Both terms remain significant after removing the distance component from $v_\parallel$. The orthogonalized model improves over the simple synodic model by $\Delta\text{AIC} = -132.6$, exceeding the original joint fit. This confirms that the $v_\parallel$ channel is not a mathematical alias of heliocentric distance once $r$ is removed.

Refinement B: Full joint regression controlling for all known heliocentric effects. A five-parameter model was fit: $\eta = \eta_0 + \eta_{\rm syn} \cos(D) + \eta_r r \cos(D) + \eta_{v_r} v_r \cos(D) + \eta_\theta \cos\theta \cos(D)$. The results expose a clear hierarchy:

• Heliocentric distance: $\eta_r = +7.11 \times 10^{-3}$ AU$^{-1}$ ($t = 1.37$, $p = 0.171$) — non-significant when all terms are controlled
• Radial velocity: $\eta_{v_r} = +1.10 \times 10^{-3}$ (km/s)$^{-1}$ ($t = 6.50$, $p < 10^{-10}$)
• CMB orientation: $\eta_\theta = -9.76 \times 10^{-4}$ ($t = -11.03$, $p < 10^{-4}$)

In this full specification, the heliocentric distance term becomes non-significant while both the radial velocity and the CMB orientation terms remain significant. This indicates that the scalar temporal topology is dynamically responsive to motion (local heliocentric $v_r$ and the velocity projection defined using $\hat{\mathbf{n}}_{\rm CMB}$) and carries a statistically resolved directional modulation (cos$\theta$), rather than being a simple static $1/r$ potential. The full model improves over the distance-velocity model by $\Delta\text{AIC} = -119.2$.

Refinement C: Marginal contribution of cos$\theta$ to the best heliocentric model. To quantify the independent contribution of the monthly anisotropy, a nested comparison was performed between the distance-velocity model (the heliocentric radial velocity analysis, 4 parameters) and the same model with $\cos\theta$ added (5 parameters). The F-test for the additional parameter gives $F(1, 25,440) = 121.49$ ($p = 1.1 \times 10^{-16}$). The partial coefficient for $\cos\theta$ is $\eta_\theta = -9.76 \times 10^{-4}$ ($t = -11.03$), and the cos$\theta$-augmented model is preferred by $\Delta\text{AIC} = -119.2$. This establishes that the monthly CMB orientation anisotropy is not merely compatible with the existing heliocentric modulation — it provides a material improvement over the best purely heliocentric model.

Directional specificity test. To verify that the anisotropy about $\hat{\mathbf{n}}_{\rm CMB}$ has dipole-like structure rather than being an arbitrary sky pattern, the analysis was repeated with three geometrically defined reference directions: the true anti-CMB antipode (flipped RA and Dec) and two directions rigorously perpendicular to the CMB dipole axis (obtained via vector cross-products with the celestial pole). For a pure dipole, the anti-CMB direction should exhibit the same magnitude with reversed sign, while perpendicular directions should be suppressed. The results confirm this pattern:

• Anti-CMB: $\Delta\eta = +1.76 \times 10^{-3}$ ($7.92\sigma$), $|\Delta\eta_{\rm anti}| / |\Delta\eta_{\rm CMB}| = 1.000$ (exact dipole antisymmetry)
• Perpendicular-1: $\Delta\eta = -5.66 \times 10^{-4}$ ($2.76\sigma$), suppressed to $0.35\times$ the CMB amplitude
• Perpendicular-2: $\Delta\eta = -4.58 \times 10^{-4}$ ($2.04\sigma$), suppressed to $0.26\times$ the CMB amplitude

The anti-CMB result confirms the dipole nature of the anisotropy: reversing the axis reverses the sign while preserving the magnitude to within 0.1%. The perpendicular directions are substantially weaker ($\approx 2.4\sigma$ vs $7.9\sigma$), though still marginally significant.

Joint-model null-direction control. The directional specificity test above uses a simple split-test (high vs low cos$\theta$), which is informative but weaker than testing whether the full joint regression itself degrades when the dipole direction is rotated. To address this, the complete five-parameter joint model (synodic + heliocentric distance + radial velocity + orientation) was re-fit with three geometrically exact rotations of the CMB dipole axis: two independent 90^\circ perpendicular directions (one in the equatorial plane, one in the meridian plane) and the 180^\circ true antipode ($-\hat{\mathbf{n}}_{\rm CMB}$):

• True CMB: $\eta_\theta = -9.76 \times 10^{-4} \pm 8.85 \times 10^{-5}$ ($t = -11.03$), AIC = $-141,929.1$
• 90^\circ perpendicular-1 (equatorial plane): $\eta_\theta = -2.15 \times 10^{-4} \pm 9.78 \times 10^{-5}$ ($t = -2.20$, $p = 0.028$), AIC = $-141,812.6$ ($\Delta$AIC = $+116.5$ vs true)
• 90^\circ perpendicular-2 (meridian plane): $\eta_\theta = -1.12 \times 10^{-3} \pm 2.17 \times 10^{-4}$ ($t = -5.17$, $p < 10^{-4}$), AIC = $-141,834.5$ ($\Delta$AIC = $+94.6$ vs true)
• 180^\circ antipode: $\eta_\theta = +9.76 \times 10^{-4} \pm 8.85 \times 10^{-5}$ ($t = +11.03$), AIC = $-141,929.1$ ($\Delta$AIC = $0.0$ vs true)

Both perpendicular directions are strongly disfavored relative to the fitted $\hat{\mathbf{n}}_{\rm CMB}$ orientation in this five-parameter construction ($\Delta$AIC = $+116.5$ and $+94.6$ respectively), confirming that the joint model prefers that axis over those particular $90^\circ$ rotations. This does not establish uniqueness among all celestial directions: the falsification testing (G2) refits the same joint model with the Earth–Moon cosine taken about the ICRS galactic north pole and finds $\Delta\text{AIC} = 125.7$ (vs.\ $119.2$ for the Planck axis) with $t = -11.31$ on $\eta_\theta$ in that single-axis fit. The north ecliptic pole and the ephemeris-averaged Earth–Sun orbital angular-momentum direction yield much smaller $\Delta\text{AIC}$ ($\approx 27.5$; $|t| \approx 5.4$). the falsification testing (H) tests whether the galactic advantage is a separable second axis by fitting both orientation terms simultaneously: the epoch-wise predictors are almost collinear, $r(\cos\theta_{\rm CMB}, \cos\theta_{\rm gal}) = 0.984$ (sky-axis separation $42^\circ$). In the six-parameter model, $\eta_{\rm CMB} = +3.50 \times 10^{-4}$ ($t = 0.68$, non-significant) and $\eta_{\rm gal} = -1.51 \times 10^{-3}$ ($t = -2.63$). A nested F-test adds galactic north given CMB ($F = 6.90$, $p = 0.0086$) but not CMB given galactic north ($F = 0.47$, $p = 0.50$). The dual model AIC ($\Delta\text{AIC} = 124.1$) is between the single-axis values and does not beat galactic-only ($125.7$). Thus the single-axis galactic ``win'' reflects shared orientation variance in a nearly one-dimensional subspace, not a second independent cosmological axis distinct from the Planck dipole. The second perpendicular direction shows a stronger residual ($t = -5.17$) than the first ($t = -2.20$), consistent with the directional specificity split-test where the meridian-plane perpendicular was also the more significant of the two nulls. The 180^\circ antipode performs identically ($\Delta$AIC = $0.0$) with a coefficient that is exactly sign-reversed and equal in magnitude to the true CMB value to all reported precision, as required by dipole antisymmetry about the chosen axis. This joint-model control strengthens the dipole geometry about $\hat{\mathbf{n}}_{\rm CMB}$ beyond the split-test; it does not, by itself, prove that the Planck dipole is the only admissible physical origin.

Higher-order multipole test. To determine whether the residual perpendicular signals arise from higher-order multipole components, a joint regression including dipole ($\cos\theta$), quadrupole ($P_2(\cos\theta) = (3\cos^2\theta - 1)/2$), and octupole ($P_3(\cos\theta) = (5\cos^3\theta - 3\cos\theta)/2$) terms was fitted, controlling for synodic, distance, and radial velocity effects:

• Dipole: $\eta = -1.04 \times 10^{-3} \pm 9.92 \times 10^{-5}$ ($t = -10.48$; CMB-multipole robustness check)
• Quadrupole: $\eta = +1.48 \times 10^{-4} \pm 5.12 \times 10^{-5}$ ($t = 1.31$, $p = 0.19$, non-significant; CMB-multipole robustness check)
• Octupole: $\eta = +2.14 \times 10^{-4} \pm 1.39 \times 10^{-4}$ ($t = 1.54$, $p = 0.12$, non-significant; CMB-multipole robustness check)

The joint F-test for adding quadrupole + octupole to the dipole-only model gives $F(2, 25,438) = 1.82$ ($p = 0.163$), showing no evidence for higher-order multipole contributions. The perpendicular residuals are therefore not due to quadrupolar or octupolar anisotropy; they more likely reflect residual systematics or station-dependent observational geometry.

Bootstrap robustness. A bootstrap resampling analysis ($n = 200$) of the full joint regression confirms the stability of the key coefficients. The 95% confidence intervals are:

• CMB orientation: $\eta_\theta = -9.78 \times 10^{-4} \pm 8.26 \times 10^{-5}$, 95% CI $[-1.13, -0.82] \times 10^{-3}$ (excludes zero)
• Radial velocity: $\eta_{v_r} = +1.11 \times 10^{-3} \pm 1.66 \times 10^{-4}$, 95% CI $[+0.74, +1.42] \times 10^{-3}$ (excludes zero)
• Heliocentric distance: $\eta_r = +6.82 \times 10^{-3} \pm 5.12 \times 10^{-3}$, 95% CI $[-2.68, +17.9] \times 10^{-3}$ (includes zero)

Bootstrap confirms that the CMB orientation and radial velocity coefficients are robust and reproducible, while the distance coefficient is consistent with zero.

Robustness validation suite. To address concerns that the CMB anisotropy could arise from statistical artifacts, five scored robustness checks (A–E) were performed, with supplementary year jackknife and alternative fixed-axis controls (G1–G2) recorded in the same JSON output:

• Aliasing simulation. 5,000 datasets were simulated with a realistic synodic signal and no CMB dependence. The spurious $\cos\theta$ coefficient induced solely by the $r(\cos D, \cos\theta) = 0.050$ correlation has mean $9.4 \times 10^{-7}$ and standard deviation $8.8 \times 10^{-5}$, with maximum $|t| = 3.70$. The observed $\eta_\theta = -9.76 \times 10^{-4}$ ($t = -11.03$) lies beyond the 99.9th percentile ($p < 10^{-3}$), ruling out aliasing.
• Multicollinearity diagnostics. All VIFs in the full joint model are at most ~1.20 (cosD: 1.017; $r \cdot \cos D$: 1.026; $v_r \cdot \cos D$: 1.188; $\cos\theta \cdot \cos D$: 1.200), and the condition number is $\kappa = 1.78 \times 10^{2}$, well below the severe-multicollinearity threshold. The standard error of $\eta_\theta$ is inflated by only $1.00\times$ relative to an orthogonal predictor.
• Permutation test. With $\cos\theta$ permuted 1,000 times, the null distribution of $\eta_\theta$ is centered at zero ($\sigma = 8.0 \times 10^{-5}$) with 99.9th percentile $|t| = 2.75$; the observed $t = -11.03$ is far in the tail ($p < 10^{-3}$).
• Gram–Schmidt orthogonalization. After explicitly removing from $\cos\theta \cdot \cos D$ all components collinear with $\cos D$, $r \cdot \cos D$, and $v_r \cdot \cos D$, the residual coefficient is unchanged: $\eta_\theta = -9.76 \times 10^{-4} \pm 8.85 \times 10^{-5}$ ($t = -11.03$), confirming the signal resides in the mathematically independent component.

Primary orientation identifiability (the falsification testing, test H). Because galactic north attains slightly larger $\Delta\text{AIC}$ in single-axis comparisons, the primary orientation paragraph is the six-parameter joint model with both $\cos\theta_{\rm CMB} \cdot \cos D$ and $\cos\theta_{\rm gal} \cdot \cos D$. The two orientation time series correlate at $r = 0.984$, so they occupy a nearly shared one-dimensional subspace despite a $42^\circ$ separation between the fixed sky axes. Conditional on both terms, only the galactic coefficient remains marginally resolved ($t = -2.63$); the Planck coefficient collapses to $t = 0.68$ ($p \approx 0.50$). This resolves the G2 ambiguity: the data support one strong fixed-sky dipole-like degree of freedom in this channel, not two independent orientation channels. Single-axis Planck vs galactic $\Delta\text{AIC}$ rankings are reported only as supplementary checks below.

TEP $\eta_\theta$ prediction coverage (Steps 033 and 055). the background variance analysis fixes the synodic response coefficient $\eta_0 \approx -3.18 \times 10^{-4}$. Linearizing the TEP orientation law $\eta = \eta_0(1 + m_{\rm CMB}\cos\theta_{\rm EM\text{-}CMB})$ in the synodic channel predicts $\eta_\theta \approx \eta_0 m_{\rm CMB}$ with $|m_{\rm CMB}| \sim \mathcal{O}(1)$, so $|\eta_\theta|$ should lie in the same $10^{-4}$–$10^{-3}$ decade as $|\eta_0|$. The Planck-axis joint fit gives $\eta_\theta = -9.76 \times 10^{-4}$, inside the background-variance-analysis decade band $[-10^{-3}, -10^{-4}]$ with $|\eta_\theta|/|\eta_0| \approx 3.1$, consistent with $|m_{\rm CMB}| \sim \mathcal{O}(1)$–$\mathcal{O}(3)$ rather than a new magnitude class. The point estimate exceeds $|\eta_0|$ as expected when the orientation channel carries a substantial fraction of the synodic-scale coupling.

Random-axis Monte Carlo on $S^2$ (the falsification testing, test D). To test whether any fixed sky axis could match the Planck $\Delta\text{AIC}$, 100,000 directions were drawn uniformly on the sphere and the base-plus-direction model was refit for each. The Planck axis yields $\Delta\text{AIC} = 119.2$ over the heliocentric base; the random-axis null has median $\Delta\text{AIC} = 61.9$, 99th percentile $130.3$, and maximum $137.7$. The empirical fraction of axes with $\Delta\text{AIC}$ at least as large as Planck is $p_{\Delta{\rm AIC}} = 0.096$ (Jeffreys 95% interval $[0.094, 0.097]$). The Planck axis therefore lies near the $1 - p_{\Delta{\rm AIC}} \approx 90\%$ quantile of the ``any fixed axis'' null, not at a strict $\alpha = 0.01$ outlier. Correlation-matched scrambling ($|r(\cos D,\cos\theta)|$ within $0.02$ of $|r| \approx 0.050$) gives a more conservative $p_{\Delta{\rm AIC}} \approx 0.266$. A composite directional specificity score combining aliasing, permutation, phase-null, orthogonalization, and directional-anatomy tests via Stouffer's Z method yields $Z = 1.45$ ($p = 0.074$), quantifying the combined evidence for directional specificity across independent test geometries. Three supplementary refined nulls (D3–D5; 5,000 draws each) sharpen this interpretation:

• Synodic phase null. Circular shifts of the CMB $\cos\theta_c$ time series (preserving the marginal $\cos\theta$ distribution but destroying alignment with $\cos D$) collapse the directional $F$-gain: median $F_{\rm eff} = 1.64$ versus $F_{\rm eff}^{\rm CMB} = 50.1$; only 3 of 5,000 shifts reach the observed level ($p_{F,{\rm eff}} = 6 \times 10^{-4}$; Jeffreys 95% interval $[1.7 \times 10^{-4},\, 1.6 \times 10^{-3}]$). The joint dipole term is therefore tied to synodic phase coupling, not to the marginal sky-angle histogram alone.
• Gram–Schmidt orthogonalized scramble. Random sky directions with $\cos\theta \cdot \cos D$ orthogonalized against the heliocentric base before augmentation (the same geometry as test E) yield $p_{F,{\rm eff}} \approx 0.090$ (Jeffreys interval $[0.082,\, 0.098]$), comparable to the uniform-direction null.
• CMB-axis rotation nulls. Uniform $\mathrm{SO}(3)$ rotations of the Planck axis reproduce the uniform-scramble marginality ($p_{F,{\rm eff}} \approx 0.096$). A local $30^\circ$ cone about the Planck axis is even less specific ($p_{F,{\rm eff}} \approx 0.30$), so fine-tuned in-sample axis uniqueness is not established at strict $\alpha$.

A calendar-year jackknife (36 yearly excisions) leaves the sign of $\eta_\theta$ unchanged with median $|t| \approx 10.8$ on the orientation coefficient, supporting temporal stability of the directional term within this archive.

Supplementary (single-axis $\Delta\text{AIC}$ and dipole anatomy). Independent $90^\circ$ rotations of $\hat{\mathbf{n}}_{\rm CMB}$ in the joint model degrade the fit by $\Delta\text{AIC} = +116.5$ and $+94.6$, while the $180^\circ$ antipode matches the true dipole in magnitude and $\Delta\text{AIC}$. The uniform random-axis null above uses the same regression geometry for all directions on $S^2$, not only these geometric rotations. Together with aliasing, permutation, and orthogonalization robustness checks, the Planck-aligned structure is consistent with a dipole component of the joint dynamical model, while strict uniqueness among all fixed axes remains only marginally established ($p_{\Delta{\rm AIC}} \approx 0.095$).

Annual envelope of monthly anisotropy. Because the CMB dipole direction is fixed in the sky while the Earth-Moon system's observability varies annually, the monthly $\cos\theta$ effect should exhibit an annual envelope. A model including synodic, cos$\theta$, and cos$\theta$ × annual-phase interaction terms was fit. Both envelope terms are significant:

• Sinusoidal envelope: $\eta_{\rm env,\sin} = +1.08 \times 10^{-3}$ ($t = 8.25$, $p < 10^{-4}$)
• Cosinusoidal envelope: $\eta_{\rm env,\cos} = -3.18 \times 10^{-4}$ ($t = -2.50$, $p = 0.012$)

The joint F-test for both envelope terms gives $F(2, 25,440) = 15.08$ ($p = 2.8 \times 10^{-7}$). This annual modulation of the monthly anisotropy is expected whenever the orientation predictor is tied to a fixed celestial axis (so observability and annual sampling modulate the monthly geometry), and it does not by itself distinguish between a cosmological interpretation and station- or geometry-mediated annual structure.

Cross-station consistency. The monthly orientation anisotropy was tested independently at each observing station with sufficient data ($N \ge 500$). In addition to the simple split test, a full joint regression (synodic + distance + radial velocity + $\cos\theta$) was run per station to test whether the global "clear hierarchy" replicates:

Simple split-test results:

• Grasse (France): $\Delta\eta = -2.33 \times 10^{-3}$ ($9.71\sigma$, $N = 18,742$)
• Haleakala (USA): $\Delta\eta = +8.15 \times 10^{-3}$ ($2.03\sigma$, $N = 666$)
• APO (USA): $\Delta\eta = +2.60 \times 10^{-4}$ ($0.77\sigma$, $N = 2,595$)
• McDonald2 (USA): $\Delta\eta = -1.52 \times 10^{-3}$ ($1.30\sigma$, $N = 3,097$)

Station-level full-joint regression coefficients:

• Grasse ($N = 18,742$): $\eta_r = +2.13 \times 10^{-2}$ ($t = 3.66$), $\eta_{v_r} = +1.40 \times 10^{-3}$ ($t = 7.52$), $\eta_\theta = -1.26 \times 10^{-3}$ ($t = -13.07$). The clear hierarchy replicates: distance is significant here (unlike globally), but $\eta_\theta$ dominates with the highest t-statistic.
• Haleakala ($N = 666$): $\eta_r = +4.94 \times 10^{-2}$ ($t = 0.66$), $\eta_{v_r} = -6.56 \times 10^{-4}$ ($t = -0.22$), $\eta_\theta = +3.84 \times 10^{-3}$ ($t = 2.57$). Only $\eta_\theta$ is significant; the small sample size produces large uncertainties.
• APO ($N = 2,595$): $\eta_r = -3.01 \times 10^{-2}$ ($t = -3.81$), $\eta_{v_r} = +8.78 \times 10^{-4}$ ($t = 3.56$), $\eta_\theta = -4.63 \times 10^{-4}$ ($t = -3.48$). All three coefficients are individually significant, but $\eta_\theta$ and $\eta_{v_r}$ have comparable strength.
• McDonald2 ($N = 3,097$): $\eta_r = -3.16 \times 10^{-2}$ ($t = -1.55$), $\eta_{v_r} = -1.17 \times 10^{-3}$ ($t = -1.50$), $\eta_\theta = +5.47 \times 10^{-4}$ ($t = 1.32$). None reach formal significance; limited statistical power.

Grasse dominates the raw archive (74%; 73.7% of the cleaned sample) and drives the global significance. At Grasse, the station-level joint regression supports the global pattern: the CMB orientation coefficient ($t = -13.07$) is the strongest individual term. APO shows all three coefficients significant, though with $\eta_r$ unexpectedly strong, possibly reflecting station-specific systematics. Haleakala and McDonald2 have limited statistical power ($N = 666$ and $3,097$) and cannot independently resolve the full hierarchy. The mixed cross-station pattern does not invalidate the global result, yet indicates that station-level systematic modelling warrants attention, particularly for APO where $\eta_r$ is anomalously large.

Monthly orientation anisotropy. Splitting the data by the cosine of the Earth-Moon to CMB angle:

• Earth-Moon aligned with CMB dipole ($\cos\theta \ge 0.75$): $\eta = -8.44 \times 10^{-4}$ ($5.69\sigma$, $N \approx 3,800$; CMB-binned-anisotropy robustness check)
• Earth-Moon anti-aligned with CMB dipole ($\cos\theta \le -0.75$): $\eta = +9.21 \times 10^{-4}$ ($5.54\sigma$, $N \approx 3,800$; CMB-binned-anisotropy robustness check)

The differential evaluates to $\Delta\eta = -1.76 \times 10^{-3}$ at $7.92\sigma$. The Earth-Moon orientation relative to the CMB dipole modulates the Nordtvedt parameter with an amplitude comparable to the primary synodic signal itself. The correlation between $\cos\theta_{\rm EM-CMB}$ and $\cos(D)$ is only $r = 0.050$, confirming that the monthly anisotropy is a genuinely independent physical effect, not a mathematical alias of the synodic modulation.

CMB-phase annual signal. Because the CMB dipole direction (ecliptic longitude $\approx 173^\circ$) is offset by approximately $70^\circ$ from the perihelion longitude ($\approx 103^\circ$), an annual signal at the CMB phase was tested using both sin and cos harmonics to remain phase-independent. The joint F-test gives $F(2, 25,441) = 5.79$ ($p = 0.0031$), indicating significant annual power at the CMB dipole phase beyond the synodic modulation. The sin component is $\eta = +2.71 \times 10^{-4}$ ($t = 3.07$, $p = 0.0022$; CMB-annual-signal robustness check), while the cos component is $\eta = +1.10 \times 10^{-4}$ ($t = 1.47$, $p = 0.14$; CMB-annual-signal robustness check). The dominance of the sin component reflects the ellipticity of Earth's orbit: for a circular orbit the signal would be purely cosinusoidal, but the $e \approx 0.017$ eccentricity introduces a strong quadrature term. This annual-phase signature is a discriminant that no purely heliocentric mechanism can reproduce at the CMB dipole longitude.

Binned anisotropy trend. Eight bins in $\cos\theta_{\rm EM-CMB}$ yield a linear trend $d\eta/d\cos\theta = -1.05 \times 10^{-3} \pm 2.21 \times 10^{-4}$ (CMB-anisotropy-slope robustness check) ($t = -4.75$, $p = 0.0032$), confirming that the anisotropy is smooth and monotonic across the full range of orientations rather than confined to extreme bins.

4.12.3 Canonical Full-Systematic Extraction

The nested systematic control ladder is fit on the $6\sigma$ MAD-cleaned INPOP19a sample ($N=25{,}445$) with all five stations retained. The canonical primary estimand is the full-systematic model m5_full_corrected, which augments $\cos(D)$ with $\cos(2D)$, annual, and monthly nuisance terms:

at $6.17\sigma$ significance, with cluster-robust uncertainty $6.52\sigma$ across the five observatories. The cosD-only ladder baseline on the same cleaned sample yields $\eta = -3.18 \times 10^{-4} \pm 6.04 \times 10^{-5}$ at $5.25\sigma$ (naïve OLS robustness check). The canonical estimator hierarchy fixes this hierarchy; multiplicity control treats the full-systematic extraction as the primary independent hypothesis for multiple-testing correction.

4.12.4 Per-Station Power and Observed Extraction

The station power analysis compares expected per-station SNR at the cosD-only reference amplitude $\eta = -3.18 \times 10^{-4}$ with observed single-station $\cos(D)$ regressions on the full 26,207-point INPOP19a sample. The expected column is a design calculation; the observed column is the measured station-wise $\eta$.

Subset	N	RMS (cm)	Expected SNR	Observed $\eta$	Observed SNR
Full pooled (cosD-only)	26,207	—	—	$-3.18 \times 10^{-4}$	$5.25\sigma$
Grasse	19,390	9.9	$2.92\sigma$	$-5.39 \times 10^{-4}$	$4.97\sigma$
APO	2,595	3.2	$3.67\sigma$	$-2.39 \times 10^{-4}$	$2.77\sigma$
McDonald2	3,139	9.6	$0.88\sigma$	$-5.00 \times 10^{-4}$	$1.39\sigma$
Matera	346	6.2	$0.37\sigma$	$-1.81 \times 10^{-5}$	$0.02\sigma$
Haleakala	737	13.8	$0.22\sigma$	$+3.55 \times 10^{-3}$	$2.45\sigma$

No single station reaches conventional standalone significance; the headline detection is carried by the pooled full-systematic model (Section 4.14.3). Grasse-C-SPAD hardware epochs reach $2.37\sigma$–$8.83\sigma$ depending on the partition, with negative $\eta$ in every powered epoch.

4.12.5 Precision-Weighted Station Regression

The station power analysis implements inverse-variance weighting by station ($1/\sigma_{r,s}^2$) on the full sample, down-weighting high-RMS and phase-truncated observatories without removing Haleakala from the pooled design. The precision-weighted estimate is:

at $6.27\sigma$ significance. This bounds the amplitude under quality-weighted pooling and complements the headline precision-weighted extraction (Section 4.14.3). The clean-subset Grasse C-SPAD-era plus APO analysis is reported separately in Section 4.29.

4.12.6 Frequency Domain Orthogonality and Sideband Harmonics

To establish why standard LLR integrators fail to completely absorb the TEP signal, a spectral orthogonality analysis was executed. TEP dictates a dynamically scaling interaction modulated by the background scalar potential, implying a modulation depth $m$ of the effective coupling. Empirical calibration against perihelion-aphelion data (perihelion-aphelion mean differential calibrated via two-point derivation) reveals a directionally negative perihelion excess over aphelion on the cleaned INPOP19a sample, with a $1.42\sigma$ differential in the cosD-only split. Even using a conservative $m = 1.0$ model, the multiplication of synodic ($D$) and orbital ($l'$) frequencies deposits approximately 33% of the signal power into composite sidebands at $D - l'$ and $D + l'$.

In physical terms, the TEP signal produces a characteristic spectral "fingerprint" at 32.13 days ($D - l'$ sideband) that standard solar system integrators cannot model because they lack the appropriate mathematical basis functions. This frequency resides in a spectral gap between standard lunar evection (31.81 days) and anomalistic motion, leaving the TEP variance structurally unabsorbed.

Lomb-Scargle periodograms confirm that these TEP sidebands are spectrally orthogonal to the Keplerian parameter space accessible to static-$\eta$ solvers. The $D - l'$ sideband at 32.13 days is resolved from the 31.81-day evection by $4\times$ the Rayleigh limit ($1/T \approx 7.8 \times 10^{-5}$ cpd). Because standard integrators lack the functional degrees of freedom to parameterize power at these specific sideband frequencies, the dynamic TEP variance is structurally preserved in the post-fit residuals where it is recovered by this analysis pipeline.

4.19 Quantitative $\eta$ Prediction

A quantitative prediction for $\eta$ was derived from the TEP geometric Temporal Shear formalism. Using the Earth-Moon compactness-squared differential ($\Phi_\oplus^2 - \Phi_{\rm Moon}^2$) and the Observable Response Coefficient framework, the predicted Nordtvedt violation range is order-of-magnitude $\eta \in [-10^{-3}, -10^{-4}]$. The measured AR(1) GLS cluster-robust parameter $\eta = -3.28 \times 10^{-4}$ (cosD-only AR(1) GLS robustness check) lies within this predicted theoretical range, consistent with the same framework that yields preliminary galactic-scale coefficients $\kappa_{\rm Cep} \sim 10^6$ mag and $\kappa_{\rm MSP} \sim 10^6$–$10^7$ from related work (Papers 10 and 11). The calibration constant $C$ derived self-consistently from the measurement confirms the TEP potential and coupling amplification required by the TEP framework.

4.21 Decoupling Thermal Array Deformation

To eliminate the hypothesis that the synodic-range effect at $13|\eta| \approx 5.1$ mm (headline precision-weighted scale) is a structural artifact driven by thermal expansion of the lunar retroreflector arrays (which heat up near full moon), a worst-case physical model was constructed. The Apollo arrays, constructed with an aluminum housing (~0.15m thick) subject to a ~300 K delta across the lunar day-night cycle, experience a theoretical maximum expansion of 1.027 mm. This thermal expansion is an order of magnitude too small—accounting for only 10.3% of the primary synodic amplitude at $13|\eta|$. Thus, thermal array deformation cannot explain the anomaly.

4.22 Leverage Temporal Clustering

The statistical divergence between OLS ($\eta = -3.18 \times 10^{-4}$) and Theil-Sen ($\eta = -2.94 \times 10^{-4}$; cosD-only, Cook's-Distance excision) is driven by heteroskedasticity, not merely leverage. By calculating Cook's Distance and binning across 5-year epochs, a pronounced temporal clustering of high-leverage points was revealed: 64.9% occurred between 1984–1989, a period contributing only 8.4% of the total data. This era is dominated by early Grasse observations with large RMS (12.6 mm vs modern 2.5 cm-class C-SPAD), confirming that early hardware systematics inflated OLS variance. However, Theil-Sen does not "suppress" this variance in a physically meaningful way; rather, the median of pairwise slopes is pulled toward the high-variance segment, yielding a negatively biased amplitude. The weighted biweight M-estimator properly compensates for station heteroskedasticity via inverse-variance weights and recovers $\eta = -3.64 \times 10^{-4}$ ($9.15\sigma$ cluster-robust; a weighted biweight M-estimator), confirming that the Theil-Sen discrepancy is an estimator artefact.

4.15 False-Positive Diagnostic Results

Two dedicated false-positive diagnostic steps were executed on the full 26,207-observation INPOP19a dataset to test for terrestrial and orbital systematic confounders.

4.15.3 Day/Night Thermal Bias Null Test

The local solar altitude was computed for every observation at its observing station using astropy.coordinates.get_sun. A partial regression simultaneously modelled residuals against $\cos(D)$, solar altitude, and a binary day/night indicator to test whether the daytime-ranging bias (New Moon observations are geometrically constrained to daytime; Full Moon to nighttime) could generate a spurious synodic modulation. Station-level results are given in Table 3.

Table 3: Day/Night thermal false-positive test. "Cleaned $\eta$" is the partial regression result simultaneously controlling for solar altitude and day/night indicator.

Station	N	Day-Ranged (%)	Day–Night Bias (mm)	Cleaned $\eta$ (Solar Controlled)	Cleaned p-value
APO	2,595	45.7%	+6.5	$-3.07 \times 10^{-4}$	$2.5 \times 10^{-4}$
Grasse	19,390	44.3%	$-3.1$	$-8.78 \times 10^{-4}$	$2.9 \times 10^{-17}$
Matera	346	21.7%	$-28.0$	$-2.36 \times 10^{-4}$	$0.76$
McDonald2	3,139	42.2%	$-16.6$	$+8.9 \times 10^{-5}$	$0.78$
Haleakala	737	33.4%	+41.6	$+2.18 \times 10^{-3}$	$0.062$

At the global level: the solar altitude regressor is negligible ($p = 0.281$, coefficient $= -3.1 \times 10^{-5}$ m/degree), and the cleaned global $\eta$ after controlling for it is $-6.64 \times 10^{-4}$ ($p = 2.0 \times 10^{-14}$) — an attenuation of only 2.5% relative to the uncorrected OLS value. The day/night thermal bias hypothesis is rejected by the data.

4.15.2 True Geometric Elongation Null Test

The true geocentric Sun–Moon angular separation $D_{\rm true}$ was computed for all 26,207 observations via J2000 ephemeris vectors, and both $\cos(D_{\rm mean})$ and $\cos(D_{\rm true})$ were entered as competing predictors in a single partial regression:

• Single-predictor, mean phase: $\eta = -4.52 \times 10^{-4}$ (single-predictor robustness check), SNR = $4.93\sigma$
• Single-predictor, true geometry (the geometric precision validation elongation comparison): $\eta = -4.54 \times 10^{-4}$ (single-predictor robustness check), $p = 8.2 \times 10^{-7}$
• Partial regression — $\cos(D_{\rm mean})$ coefficient: $-0.02310$, $p = 8.7 \times 10^{-17}$ (dominant, stable)
• Partial regression — $\cos(D_{\rm true})$ coefficient: $+0.01657$, $p = 2.0 \times 10^{-8}$ (partially offsetting, sign inverted)

The behavior of the true-geometry predictor is consistent with a dynamically coupled signal extracted from post-fit residuals. The INPOP19a ephemeris has fitted and subtracted all high-frequency classical Keplerian orbital dynamics (e.g., evection, variation), which operate on instantaneous $D_{\rm true}$ boundaries. Regressing the residuals against $D_{\rm true}$ re-injects post-fit classical mechanical alias noise, degrading the correlation. Conversely, the TEP scalar field operates as a smooth, long-wavelength spatial gradient scaling against the heliocentric $1/r_\odot$ well. Therefore, $D_{\rm mean}$ acts as the appropriate low-pass physical proxy for the large-scale solar field gradient. The true-geometry predictor introduces aliased mechanical noise that INPOP already resolved, which explains the sign inversion and confirms the broad physical scalar field interpretation.

Combined result (Steps 029–030): Both the day/night thermal and the geometric-proxy false-positive hypotheses are rejected by direct empirical testing on the full dataset. The $5.27\sigma$ detection is not attributable to daytime atmospheric refraction, telescope-mount thermal expansion, or aliasing from the mean-phase orbital approximation.

4.26 Frequency-Specific Null Testing

The TEP framework predicts the Nordtvedt modulation should appear at the synodic frequency (1×, period = 29.53 days) and not at other frequency factors. Testing this prediction provides discrimination between a phase-locked physical effect and broadband systematic artifacts, which would likely generate spurious power across multiple frequencies.

A comprehensive multi-frequency null scan was performed. After projecting out the synodic TEP signal to prevent leakage into neighboring frequency bins, the whitened residuals were tested against 55 non-synodic frequency factors spanning 0.4× to 3.2× the synodic frequency (excluding only the immediate synodic window 0.92×–1.08×):

Frequency Band Labeling. The low-frequency band (0.4× to 0.95× synodic, corresponding to periods of approximately 74 to 31 days) contains known systematic power from annual variations (365.25 days $\approx$ 0.081× synodic) and lunar orbital harmonics (sidereal and anomalistic periods) that are well-documented in LLR residuals. These signals arise from seasonal atmospheric variations, thermal expansion of the telescope structure, and Earth's orbital eccentricity effects. Frequencies in this band are labeled as "known systematic regions" in the reporting metadata but are still tested alongside all other candidates; the null test therefore conservatively searches the full 0.4×–3.2× range for any unexpected power, not only the cleanest null regions.

• Test frequencies: 55 frequency factors (0.4×, 0.45×, 0.5×, ..., 0.85×, 1.1×, 1.15×, 1.2×, 1.23×, 1.25×, 1.3×, ..., 3.2× synodic)
• Primary non-physical factor (1.23×): SNR $\approx 0.00\sigma$ ($\eta$ $\approx 0$), $p \approx 1.0$ (raw), $p \approx 1.0$ (Bonferroni-corrected)
• Worst-case detection: SNR $\approx 0.00\sigma$ at 0.60× synodic — well below the $5\sigma$ detection threshold
• Multiple testing correction: Bonferroni $\alpha = 1.14 \times 10^{-3}$; FDR-BH threshold = 0
• Significant detections: Zero frequencies significant after FDR-BH correction; zero after Bonferroni correction
• Null test verdict: PASS — no non-physical frequency shows significant power

The detected signal is frequency-specific: significant power appears at the synodic fundamental ($5.25\sigma$ cleaned-sample Pearson, $4.93\sigma$ OLS) with no detectable signal at 55 tested non-synodic frequencies spanning 0.4×–3.2× (maximum $\approx 0.00\sigma$). The 1.23× control factor — a theoretically motivated non-physical frequency — shows $\approx 0.00\sigma$ significance, consistent with pure noise. Crucially, the null-test frequencies were not pre-whitened from the residuals before testing; pre-whitening was applied only to remove dominant nuisance harmonics (annual, sidereal), leaving the test frequencies intact. This ensures the null test is not a self-fulfilling prophecy and provides an honest assessment of whether non-synodic frequencies carry systematic power. This frequency pattern is expected for a physical effect phase-locked to Earth-Moon-Sun geometry, whereas broadband systematic artifacts would typically generate spurious power across multiple frequencies.

4.28 Cross-Validation, Station Distribution, and Covariate Shift

Coefficient stability for the physical Nordtvedt term and out-of-sample predictive performance of the full residual model answer different scientific questions. The headline inference targets the pooled $\cos D$ coefficient under a joint nuisance design; temporal hold-out cross-validation instead asks whether every nuisance coefficient vector learned in one hardware era would transport chronologically into another. Because LLR hardware eras sample different synodic phases and atmospheric conditions, the nuisance subspace (annual, monthly, and $\cos(2D)$ aliases) is empirically not transportable across epochs: temporal hold-out predictive $R^2$ for the full-systematic model is negative on the full archive and remains negative even on the homogeneous clean subset of Section~4.29. That pattern is not treated here as a failure of signal authenticity; it is treated as direct evidence that nuisance-dominated designs cannot be sequentially re-fit across eras without bias, so primary extraction is performed on the maximally pooled dataset to marginalise era-specific covariates while estimating the physical $\cos(D)$ coefficient.

Standard cross-validation still provides a disciplined diagnostic: if a model were a pure artefact of in-sample overfitting with no coherent structure, one would generally expect predictive collapse across all honest partitions. the cross-validation analysis implements temporal hold-out (train pre-split, test post-split), random $k$-fold, leave-one-station-out, and forward-chaining. The full-systematic model ($\cos D + \cos 2D + \mathrm{annual} + \mathrm{monthly}$) yields the following predictive metrics on the cleaned INPOP19a archive:

• Temporal hold-out (1990s/2000s): $R^2_{\rm pred} = -0.19$ (m4 pooled systematics), $R^2_{\rm pred} = -0.16$ (m5 station-specific systematics)
• Temporal hold-out (pre/post 2000): $R^2_{\rm pred} = -0.14$ (m4), $R^2_{\rm pred} = -0.10$ (m5)
• Random 5-fold CV: $R^2_{\rm pred} = +0.010$ (m4), $+0.008$ (m5)
• Leave-one-station-out: mean $R^2_{\rm pred} = -0.47$ (m4), $-0.39$ (m5), but with substantial heterogeneity: Grasse hold-out $R^2 = +0.004$ (m4), Matera hold-out $R^2 = +0.006$ (m4), while APO hold-out $R^2 = -0.018$ (m4). The negative mean is driven by stations whose elongation distribution differs most from the training set.

The in-sample epoch-dependence specification (m6 in the cross-validation analysis) tests whether two independent $\cos D$ amplitudes, one restricted to pre-split data and one to post-split data, differ significantly. On the full archive the two-sided test of equality gives $z = 0.12$, $p = 0.90$: the physical $\cos(D)$ parameter is therefore compatible with a single Nordtvedt value across the split epoch, even though the nuisance-augmented design does not predict individual residuals chronologically. Negative temporal predictive $R^2$ for that full design is thus aligned with non-transportable nuisance structure rather than with a drifting Nordtvedt term.

Synthetic injection tests. A baseline ensemble injects the pipeline headline $\eta$ into per-station, per-epoch noise matched to the real residuals. The $\cos D$-only model (m1) remains temporally predictive (mean $R^2_{\mathrm{pred}} \approx +0.0078 \pm 0.0026$ over 100 trials; real temporal hold-out $+0.0058$ on the pre-2008/post-2008 split), while the full-systematic model (m4) stays near the mean-only benchmark (mean $R^2_{\mathrm{pred}} \approx +0.0019 \pm 0.0051$, with roughly half of trials negative). That baseline therefore does not reproduce the depth of the real temporal failure ($R^2_{\mathrm{pred}} \approx -0.15$). A second ensemble augments the same injection with era-conditioned nuisance draws built from the real split-era m4 nuisance subspace (see the cross-validation analysis JSON for the exact construction and RMS cap). There, temporal hold-out collapses to strongly negative mean m4 scores ($\approx -0.17 \pm 0.05$) with m4 far more negative than m1 on the same synthetic draws (m1 mean $\approx -0.026 \pm 0.017$, still mildly negative but an order of magnitude smaller in magnitude than m4). The observed real m4 temporal score ($R^2_{\mathrm{pred}} \approx -0.15$) lies near the 62nd percentile of the era-conditioned m4 ensemble. Together, the two ensembles separate ``noise alone'' from ``noise plus era-mismatched nuisance geometry,'' supporting the interpretation that epoch-specific nuisance—not absence of a fixed $\cos(D)$ carrier—drives the temporal predictive breakdown of the full model.

Rolling-window projections of an effective $\eta(t)$ and their environmental correlations (the rolling-window analysis) are isolated in Section~4.28.1 so they are not read as contradicting the global m6 stability check above: m6 tests two epoch-gated $\cos D$ coefficients inside the pooled full-systematic fit, whereas the rolling-window analysis probes shorter-window apparent variation that can mix environmental modulation with local nuisance coupling and finite-span estimation noise.

Coefficient stability versus predictive generalisation. The epoch-dependence test (m6) summarised above is the appropriate in-sample stability check for $\eta$. The Grasse-only pre/post catalogue split reported previously is likewise compatible with a common amplitude ($p = 0.63$). The nuisance basis, by contrast, is tied to how each era samples elongation and seasonal structure, so its pooled OLS coefficients need not extrapolate under temporal hold-out even when $\eta$ does not show statistically significant drift in the m6 sense.

Station Distribution Analysis. The five stations exhibit markedly different sample sizes, temporal coverage, and synodic-phase sampling:

Station	$N$	Years	RMS (cm)	Mean $|\cos D|$
Grasse	19,390	1984–2019	9.9	0.439
APO	2,595	2006–2017	3.2	0.517
McDonald2	3,139	1988–2014	9.6	0.411
Haleakala	737	1984–1991	13.8	0.355
Matera	346	2003–2019	6.2	0.252

The metric mean $|\cos D|$ quantifies phase coverage quality: values near 0.5 indicate uniform sampling across all synodic phases, while values near 0 indicate severe phase truncation. Matera's mean $|\cos D| = 0.252$ reflects extreme phase truncation; McDonald2's mean $|\cos D| = 0.411$ indicates moderate truncation. Only APO (0.517) and Grasse (0.439) achieve near-uniform phase coverage.

Covariate Shift in Temporal Hold-Out. A Kolmogorov-Smirnov test comparing elongation distributions pre- and post-2008 yields $D = 0.131$ ($p = 4.16 \times 10^{-98}$). The pre-2008 mean $|\cos D| = 0.406$; post-2008 mean $|\cos D| = 0.474$. This distributional difference means the model trained on one epoch's elongation distribution cannot extrapolate to the other. Covariate shift is therefore a documented contributor to the negative predictive $R^2$, degrading the transportability of nuisance coefficients even though the physical $\cos D$ signal is stable across epochs.

Covariate Shift in Leave-One-Station-Out. Holding out APO (test mean $|\cos D| = 0.517$ vs train $|\cos D| = 0.430$) or Matera (test $|\cos D| = 0.252$ vs train $|\cos D| = 0.441$) produces severe covariate shift in the predictor space. The KS test for APO hold-out gives $D = 0.155$; for Matera, $D = 0.432$. When the test set samples elongations the training set never saw, the regression slope extrapolates poorly even if the underlying physical signal is real.

Random $k$-Fold CV Shows Positive Predictive Power. Unlike temporal or station-stratified splits, random 5-fold CV shuffles observations across all epochs and stations. The full-systematic model achieves $R^2_{\rm pred} = +0.010$ ($p < 0.05$), demonstrating genuine predictive power when covariate shift is eliminated. Random $k$-fold CV is therefore a useful complement to temporal diagnostics: because the archive is pooled across stations and epochs, a random split preserves the marginal predictor distribution while still testing out-of-sample generalisation, without replacing the primary pooled estimator in Section~3.4.22.

Per-Station cosD-Only Regression Supports Manuscript Findings. the validation analysis independently replicates the per-station analysis reported in Section 4.10. The cosD-only OLS yields: Grasse $\eta = -5.39 \times 10^{-4}$ ($4.97\sigma$; per-station robustness check), APO $\eta = -2.39 \times 10^{-4}$ ($2.77\sigma$; per-station robustness check), McDonald2 $\eta = -5.00 \times 10^{-4}$ ($1.39\sigma$; per-station robustness check), Haleakala $\eta = +3.55 \times 10^{-3}$ ($2.46\sigma$, opposite sign, underpowered; Haleakala underpowered-station diagnostic), and Matera $\eta = -1.3 \times 10^{-5}$ ($0.02\sigma$, underpowered; per-station robustness check). These values are quantitatively consistent with the station analysis in Table 1.

Conclusion. Negative temporal and station-stratified predictive $R^2$ for the nuisance-augmented design coexists with in-sample stability of the Nordtvedt coefficient (m6, $p = 0.90$) and with mildly positive temporal predictive $R^2$ for the $\cos D$-only carrier on the full-archive split ($+0.006$), while the clean subset's $\cos D$-only temporal score stays near the mean-only benchmark ($R^2_{\mathrm{pred}} \approx -8.5 \times 10^{-4}$). Covariate shift in elongation and epoch remains a real extrapolation hazard, but the deeper structural point is that nuisance parameters are hardware-era specific while the physical amplitude is not. The two-tier synthetic programme in the cross-validation analysis shows that matched noise alone leaves m4 temporal scores near zero on average, whereas adding era-conditioned nuisance draws reproduces large negative m4 temporal scores while keeping the injected $\eta$ fixed. Random $k$-fold CV, which mixes epochs and stations, recovers weak positive predictive $R^2$ for the full model ($R^2_{\mathrm{pred}} \approx +0.010$), consistent with maximally pooled estimation being the design matched to the estimand.

Primary inference therefore remains the precision-weighted pooled extraction (Section~3.4.22 and the precision-weighted regression): temporal CV is reported as a diagnostic on nuisance transport, not as a substitute estimand for $\eta$. The clean-subset analysis in Section~4.29 sharpens the same story under controlled homogeneity.

4.29.1 Rolling-window $\eta(t)$ and environmental correlation

TEP motivates asking whether an effective coupling inferred from short spans tracks the Earth's heliocentric environment. the rolling-window analysis computes 2-year rolling-window $\eta(t)$ on the cleaned INPOP19a archive ($N = 25{,}445$, 34 windows stepped yearly) and correlates the windowed amplitudes with three environmental predictors: inverse heliocentric distance $\langle 1/r_\odot \rangle$, radial velocity $\langle v_r \rangle$, and CMB dipole orientation $\langle \cos\theta_{\rm EM\text{-}CMB} \rangle$. The weighted Pearson correlations are:

• $\mathrm{Corr}\bigl(\eta(t), \langle 1/r_\odot \rangle\bigr) = +0.045$ (weak, as expected for a quadrature variable)
• $\mathrm{Corr}\bigl(\eta(t), \langle v_r \rangle\bigr) = +0.406$ (moderate; velocity modulation of the effective coupling)
• $\mathrm{Corr}\bigl(\eta(t), \langle \cos\theta_{\rm EM\text{-}CMB} \rangle\bigr) = +0.682$ (strong; fixed-sky orientation about $\hat{\mathbf{n}}_{\rm CMB}$ dominates rolling-window variance)

A joint regression of $\eta(t)$ on all three predictors yields weighted $R^2 = 0.142$ with a significant radial-velocity coefficient in that joint fit ($\beta_{v_r} = +7.45 \times 10^{-3}$, $t = 2.14$, $p = 0.041$); the heliocentric-distance and CMB-orientation joint terms are not individually significant at the 5\% level ($p \approx 0.15$ and $p \approx 0.18$ respectively). These correlations quantify how much window-to-window variability remains after restricting to 2-year segments on the cleaned archive. They are not offered as a refutation of the pooled const-$\eta$ hypothesis tested by m6 in Section~4.28 ($p = 0.90$ on equality of pre- and post-split $\cos D$ amplitudes in the joint full-systematic model). The conservative reading is that rolling $\hat{\eta}(t)$ is a local projection that can correlate with heliocentric geometry and dipole orientation while the headline Nordtvedt coefficient estimated on the maximally pooled sample remains statistically stable across the same catalogue split.

4.31 Clean-Subset High-SNR Analysis and Orbital Orthogonality

Robustness check — Clean Subset: Grasse C-SPAD Era (2010+) + APO. Grasse with the C-SPAD detector (2010 onwards) achieves ~2 cm RMS residuals and continuous coverage. APO provides a separate high-precision station with excellent phase coverage (mean $|\cos D| = 0.517$). Excluding Haleakala (opposite sign, underpowered), Matera (severe phase truncation, $N = 346$), and McDonald2 (moderate truncation, larger noise) yields a homogeneous subset of $N = 12\,576$ observations with residual RMS 2.5 cm.

Sections~4.28--4.28.1 argued that coefficient stability for $\eta$ and temporal predictive $R^2$ for the nuisance-augmented design measure different objects. On this clean subset, temporal hold-out still yields $R^2_{\mathrm{pred}} \approx -0.016$ for the full-systematic model (train $N=2{,}112$ pre-2013, test $N=10{,}464$ post-2013), so the sign is not an artefact of poor stations alone. The same split gives $R^2_{\mathrm{pred}} \approx -8.5 \times 10^{-4}$ for the $\cos D$-only model: still slightly negative, but roughly an order of magnitude smaller in magnitude than the full-systematic $R^2_{\mathrm{pred}}$, so most of the chronological extrapolation loss is carried by the nuisance-augmented design rather than by the $\cos D$-only carrier. the clean-subset analysis also tabulates nuisance-coefficient contrasts between independent pre-2013 and post-2013 fits on this subset: the thermal alias $\cos(2D)$ shifts most sharply ($z \approx 3.4$, $p \approx 6 \times 10^{-4}$), while several seasonal terms show weaker evidence of shift. A Grasse-only split (2010--2013 versus 2013--2019, $N_{\mathrm{pre}}=479$, $N_{\mathrm{post}}=9{,}502$) yields two-sided $p \approx 0.45$ on equality of the $\cos D$ slopes, illustrating that even within one modern detector family the nuisance geometry can move while the Nordtvedt-like coefficient remains statistically compatible across eras.

Clean subset (the clean-subset analysis): nuisance coefficient transport, independent OLS on pre-2013 vs post-2013.

Term	$| \Delta \hat\beta | / \mathrm{SE}_{\Delta}$	Two-sided $p$
$\cos(2D)$	$3.43$	$6.1 \times 10^{-4}$
$\sin(2\pi t / 1\,\mathrm{yr})$	$0.21$	$0.83$
$\cos(2\pi t / 1\,\mathrm{yr})$	$1.96$	$0.050$
$\sin(2\pi t / 27.32\,\mathrm{d})$	$1.57$	$0.12$
$\cos(2\pi t / 27.32\,\mathrm{d})$	$0.80$	$0.43$

On this subset the full-systematic model gives $\eta = -3.35 \times 10^{-4} \pm 3.30 \times 10^{-5}$ ($10.2\sigma$, OLS; clean-subset robustness check) and $-3.35 \times 10^{-4} \pm 4.62 \times 10^{-5}$ ($7.25\sigma$, cluster-robust, 2 clusters; clean-subset robustness check). Station block bootstrap on the same design gives $\eta = -2.74 \times 10^{-4} \pm 1.21 \times 10^{-4}$ with $P(\eta<0)=0.965$ (2,000 draws; conservative two-station resampling bound). Grasse alone reaches $9.4\sigma$ and APO reaches $2.8\sigma$; both agree on sign and approximate magnitude.

Full archive versus clean subset amplitude. The headline precision-weighted estimand on the full cleaned INPOP19a archive ($N=25{,}445$) is $\eta = -3.91\times 10^{-4} \pm 5.63\times 10^{-5}$ ($6.94\sigma$ / $6.78\sigma$ cluster-robust), whereas the homogeneous Grasse C-SPAD (2010+) plus APO window (the clean-subset analysis; $N=12{,}576$, RMS $\approx 2.5$~cm) gives $\eta = -3.35\times 10^{-4} \pm 3.30\times 10^{-5}$ with higher OLS SNR ($10.2\sigma$) and cluster-robust $7.25\sigma$. The central values differ by $\sim 14\%$ in $|\eta|$ but lie within the mutual uncertainty implied by the canonical-estimator-hierarchy bracket: the clean window removes early high-variance hardware eras and weakly contributing stations, which shifts nuisance co-fitting and can migrate the cos(D) coefficient within the negative band. The manuscript retains the full-archive precision-weighted estimate as the primary headline because it uses the mixed-station design used in the common-$\eta$ mixed model (Section~4.25) and matches the weighting philosophy of the canonical estimator hierarchy; the clean-subset analysis is the high-SNR robustness channel, not a competing headline number.

Toy Orbital Model: Ephemeris Orthogonality. A legitimate critique of residual-based detection is that standard ephemeris fitting might absorb a cos(D) modulation by adjusting orbital elements (semi-major axis, eccentricity, argument of perigee). The toy Keplerian model tests this with a deliberately simplified 2-D coplanar model: a circular Moon orbit is perturbed by $\delta r = 13\,\eta\,\cos D$ metres, and the resulting range series is fitted with the standard Keplerian basis $\{1, \cos M, \sin M\}$, where $M$ is the mean anomaly (orbital period ~27.3 days).

Real-Data Keplerian Inclusion Proxy. On the cleaned INPOP19a residual sample, the same Keplerian basis is fit to the observed residuals before recovering the full-systematic $\cos D$ term. The Keplerian-only fit explains $R^2 \approx 0.001$ of the residual variance; after partialing $\{1,\cos M,\sin M\}$, the full-systematic recovery gives $\eta = -4.09 \times 10^{-4} \pm 6.58 \times 10^{-5}$ ($6.22\sigma$). Adding $\cos M$ and $\sin M$ directly to the joint full model gives $\eta = -4.54 \times 10^{-4} \pm 6.65 \times 10^{-5}$ ($6.82\sigma$; Keplerian joint-model inclusion proxy). This is an inclusion proxy on post-fit residuals, not a full INPOP or DE430 dynamical refit with $\eta$ inside the orbital integrator.

For $\eta = -5.4 \times 10^{-4}$ the perturbation amplitude is 7.0 mm. The orbital-element fit finds $e \approx 2 \times 10^{-14}$ (essentially zero) and the residual RMS is 4.97 mm — identical to the theoretical expectation ($13\,|\eta|/\sqrt{2} = 4.96$ mm). Regressing the residuals on $\cos D$ recovers the input $\eta$ to within $-0.0001\%$. A sweep over six $\eta$ values ($-10^{-3}$ to $-10^{-5}$) confirms faithful recovery in every case, with residual RMS scaling linearly with perturbation amplitude.

The physical reason is spectral separation: standard orbital elements generate signals at the orbital frequency and its harmonics (~27.3 days), whereas the synodic perturbation operates at the synodic frequency (~29.5 days). The two frequencies are incommensurate; $\cos D$ is orthogonal to the $\{1, \cos M, \sin M\}$ basis over spans much longer than one synodic month. In the toy Keplerian basis tested here, a synodic cos(D) perturbation of the observed amplitude therefore remains in the residuals rather than being absorbed by that basis alone.

Conclusion. The clean-subset analysis shows the signal strengthens (cluster-robust SNR rises from ~$5\sigma$ to ~$7\sigma$) when heterogeneous stations are removed. The toy orbital model shows that a synodic cos(D) perturbation of the observed amplitude is invisible to standard orbital-element fitting in the toy Keplerian basis and therefore survives in residuals if present. Together, these results address the two open predictive diagnostics: the full-archive primary fit is not reduced to a Grasse-only artifact, and the synodic term is not absorbed by the toy Keplerian basis alone. The hardware epoch and toy Keplerian analyses extend the logic to the linearized and toy-orbital bases. It must be stressed that the 82-parameter proxy is a linear regression stress test, not a symplectic N-body numerical integrator, and cannot emulate the iterative dynamic planetary potentials handled by the actual Fortran/C++ code underlying INPOP or DE430. A full INPOP or DE430 dynamical refit with $\eta$ free in the orbital integrator remains the definitive external confirmation, and collaboration with IMCCE and JPL is invited to perform such a refit.

4.30 Bounds on Solar Radiation Pressure and Atmospheric Seeing Systematics

Two potential terrestrial or local physics mechanisms — unmodeled solar radiation pressure (SRP) on the lunar retroreflectors and atmospheric seeing gradients — are evaluated here with dedicated pipeline steps that bound their maximum synodic-correlated projection into the residual channel.

4.31.1 Local Mechanical Displacement from Solar Radiation Pressure

Direct solar radiation pressure exerts a force on the Apollo lunar retroreflector arrays. Under maximally conservative assumptions — perfect specular reflection of all incident photons, an array mass of 30 kg, and an effective mounting stiffness of $10^6$ N/m (sand-like) — the elastic displacement of the array mounting is bounded at $1.03 \times 10^{-9}$ mm. Thermal re-radiation pressure from the lunar surface adds $7.39 \times 10^{-11}$ mm. The combined local mechanical displacement is $1.10 \times 10^{-9}$ mm, which is $2.7 \times 10^{-10}$ times the observed TEP signal amplitude ($\approx$ 4.12 mm). Even if the array were completely unconstrained, free displacement over one synodic month would be unphysically large ($\sim$$10^5$ m), but the arrays are bolted to bedrock. This mechanism cannot explain the anomaly.

4.31.2 Orbital SRP Scaling

Orbital solar radiation pressure on the Moon produces a geometric signature with the same $\cos(D)$ phase dependence as the Nordtvedt effect: at new moon, SRP on the Moon acts toward Earth; at full moon, away from Earth. However, orbital SRP scales as $1/r_\odot^2$ with heliocentric distance, whereas a Nordtvedt-like signal carries no such distance dependence to leading order. Three orthogonal tests avoid the collinearity trap between $\cos(D)$ and $\cos(D)/r_\odot^2$ (VIF $\approx$ 1800):

• Detrended-residual correlation. After removing the best-fit TEP $\cos(D)$ signal, the correlation between detrended residuals and the SRP proxy is $r = -3.22 \times 10^{-4}$ ($p = 0.959$, SNR = $0.05\sigma$). No SRP variance remains in the residual channel.
• Binned $1/r^2$ scaling test. The data are partitioned into 10 heliocentric-distance bins; $\eta$ is extracted in each bin and tested for $1/r^2$ scaling. The slope coefficient is marginally unstable ($t_M = -2.25$, $p = 0.054$) but its sign and magnitude are inconsistent with simple SRP. The statistical power of this test is limited by the small heliocentric-distance baseline (3.4%, 0.983–1.017 AU).
• Perihelion–aphelion differential. The 15th–85th percentile split gives $\Delta\eta = -3.40 \times 10^{-4} \pm 2.48 \times 10^{-4}$ (SNR = $1.37\sigma$), inconsistent with the SRP-predicted differential of $-1.58 \times 10^{-5}$.

INPOP19a's orbital integrator already includes standard lunar SRP modeling. An unmodeled SRP component would require both the same synodic phase-lock and a $1/r_\odot^2$ scaling that the residual channel does not exhibit. The combined the solar-radiation-pressure check evidence does not support SRP as the dominant source of the detected modulation.

4.31.3 Atmospheric Seeing

Atmospheric seeing produces fast stochastic wander of the laser beam ($\sim$1 arcsec FWHM typical, 3 arcsec in poor conditions). Over LLR integration times, this wander averages to zero. The only plausible synodic-correlated channel is an elevation-dependent systematic: because new-moon observations occur at lower solar altitudes (daytime) and full-moon observations at higher altitudes (nighttime), an elevation-dependent seeing bias could alias into $\cos(D)$. This channel is already bounded by the day/night thermal bias null test ($p = 0.281$). The threshold sweep isolates the seeing-specific contribution and finds a maximum plausible synodic-correlated range bias of 0.111 mm, or 2.7% of the TEP signal amplitude. Fast stochastic wander averages to zero; the coherent elevation channel is already falsified.

Conclusion. Solar radiation pressure, both as a local mechanical displacement ($< 10^{-9}$ mm) and as an orbital scaling effect (no supported $1/r^2$ signature), is far too small or incorrectly scaled to explain the detected modulation. Atmospheric seeing is bounded at 0.111 mm ($\sim$3% of signal). The systematic error budget in Section 3.4.16 is updated to include both bounds. These quantitative exclusions, together with the existing dust/thermal (Section 4.12) and local terrestrial systematic (Section 5.9) tests, close the currently nominated conventional explanations for the synodic-phase-locked residual modulation without claiming that every conceivable unmodeled coupling has been eliminated.

Acknowledgements

The Paris Observatory (IMCCE) is thanked for maintaining and distributing the INPOP19a LLR residual archive. Damen Albern is thanked for keeping time during the final reductions.