Global Time Echoes: 25-Year Temporal Evolution of Distance-Structured Correlations in GNSS Clocks

1. Introduction

1.1 Background and Motivation

The Global Navigation Satellite System (GNSS) represents one of humanity's most precise timing networks, with atomic clocks maintaining synchronization at nanosecond levels across thousands of kilometres (Hofmann-Wellenhof et al., 2008; Teunissen & Montenbruck, 2017). While designed for positioning and navigation, this infrastructure inadvertently provides an unprecedented natural laboratory for testing fundamental physics at planetary scales. The continuous operation of hundreds of ground-based receivers, each equipped with high-stability oscillators phase-locked to satellite atomic clocks (Senior & Ray, 2008), creates a global mesh of timing correlations that may be sensitive to subtle relativistic effects (Ashby, 2003).

The theoretical basis for this investigation is the Temporal Equivalence Principle (TEP), which predicts that motion through gravitational fields should induce measurable modulations in the correlation structure of distributed timing networks. Unlike classical relativistic effects that manifest as clock rate differences, TEP predicts that the correlation between clocks—their tendency to maintain phase coherence—should vary with the local gravitational environment and the system's velocity through that environment.

Critically, we propose that standard GNSS processing algorithms, designed to remove energetic (common-mode) errors via datum constraints, inadvertently preserve the subtle, geometry-dependent (differential) correlations that are the focus of this work. This transforms the use of processed data from a limitation into a mechanistic filter: the absence of classical gravitational scaling (GM/r²) is not a failure of detection, but a predicted consequence of the processing chain, while the geometric signatures (anisotropy, alignment) are transmitted.

1.2 Theoretical Context and Interpretive Framework

TEP in one line: The Temporal Equivalence Principle (TEP) treats proper time as a dynamical field: instead of a fixed background parameter, time behaves like a scalar field woven through spacetime, whose local configuration and gradients govern clock rates, correlations, and the effective speed of light.

This investigation is fundamentally grounded in empirical phenomenology: we systematically document systematic patterns in GNSS timing correlations over decadal timescales, then interpret these findings through established theoretical frameworks. TEP provides one such framework for interpreting velocity-dependent and gravitational coupling effects in distributed timing networks.

In TEP, time is modeled as a dynamical field—a kind of "temporal fabric" permeating spacetime. Clocks do not simply read out a passive parameter; they interact with this fabric. Variations in the temporal field, and in the system's motion through it, offer a potential explanation for the correlation structure modulations we document empirically.

Formally, the temporal field is treated as a covariant scalar; it does not introduce a preferred frame.

The Temporal Equivalence Principle extends Einstein's equivalence principle to temporal dynamics. Key predictions—velocity-dependent anisotropy and event responses—provide testable signatures that we investigate empirically. The framework proposes that:

1. Velocity-Dependent Coupling: The correlation decay length between synchronized clocks should depend on their collective velocity relative to the dominant gravitational frame
2. Gravitational Modulation: Changes in the local gravitational field configuration should modulate the phase coherence of timing networks
3. Geometric Anisotropy: The correlation structure should exhibit directional dependence aligned with the system's motion through spacetime

These predictions are quantitatively distinct from conventional general relativistic effects:

• Classical GR predicts clock rate differences proportional to gravitational potential (Δf/f ~ GM/rc²)
• TEP predicts correlation structure modulation with characteristic decay length λ ~ λ₀(1 + v²/c²)

TEP Predictions vs. Observations

The TEP framework makes specific quantitative predictions that can be tested against GNSS observations. Table 1 summarizes these a priori predictions and their experimental verification:

Key Finding: All seven TEP predictions are consistent with GNSS observations—including the critical null prediction that classical gravitational scaling (GM/r²) should be absent in processed data due to common-mode filtering (§2.1.3). This discriminating prediction distinguishes TEP from naive gravitational coupling hypotheses. The geometric modulation amplitude (19% in correlation ratios) exceeds naive energetic estimates (~10⁻⁶), but this discrepancy is resolved by the processing filter mechanism: energetic signals are filtered while geometric structure is transmitted (§4.3.1.1).

1.3 Previous Work

The multi-centre study established the phenomenon across independent analysis centers (CODE, IGS, ESA), reducing the likelihood of center-specific processing artifacts over 2.5 years. That study demonstrated velocity-dependent anisotropy and anomalous planetary event responses, using rigorous cross-center validation and null tests. However, the limited temporal baseline precluded the investigation of long-period geophysical or astronomical cycles.

Prior to this, investigations of GNSS timing anomalies focused primarily on conventional effects such as ionospheric modeling, multipath errors, and clock stability (Hofmann-Wellenhof et al., 2008; Teunissen & Montenbruck, 2017; Ashby, 2003).

1.4 Study Design: Replication, Extension, and Falsification Tests

TEP-motivated empirical test: This work was explicitly conceived as a theory-driven empirical test of the Temporal Equivalence Principle using long-span GNSS clock products. The analysis pipeline, choice of observables (anisotropy, orbital modulation, planetary events, nutation, mesh coherence), and the a priori prediction table in §1.2 were defined to test whether GNSS timing correlations exhibit the specific spatial, temporal, and geophysical patterns forecast by TEP.

Methodological Foundation: This study builds upon the multi-centre study's comprehensive validation framework, which established cross-center consistency (R² = 0.920-0.970 between CODE, IGS, ESA) through 388 statistical tests and extensive null testing. That foundational work demonstrated the existence of systematic distance-structured correlations in GNSS timing networks over 2.5 years across three independent processing centers. This study extends the temporal baseline by a factor of 10× (25.3 years) to address three critical scientific questions that distinguish transient artifacts from genuine physical coupling:

Primary Empirical Objectives:

1. Objective A (Replication): Document temporal stability of core phenomenology (anisotropy, orbital coupling, event responses) over 25.3 years.
2. Objective B (Extension): Test long-period geophysical signatures (18.6-year nutation, >20 Chandler wobble cycles) inaccessible in short baselines.
3. Objective C (Characterization): Enhance planetary event statistics across two decades with improved seasonal coverage.
4. Objective D (Falsification): Identify operative cosmological reference frame (CMB vs galactic) via grid search and test theoretical selectivity with null controls (solar rotation) and robustness checks (hemisphere stratification).

1.5 Evidence Handoff: Building on Multi-Center Validation

This study builds directly on the multi-centre study's comprehensive validation framework. The table below maps each major claim to its validation status in Paper 1 versus extensions in this work:

Key Empirical Insight: All core phenomenology from Paper 1 is re-observed here with enhanced statistical power and temporal stability. The multi-centre study's comprehensive validation framework—including 388 statistical tests, extensive null testing, ionospheric controls, and cross-center consistency—provides the methodological foundation upon which these temporal extensions are built. The convergence of independent findings across multiple analysis centers and extended timescales strongly counters artifact concerns.

2. Data and Methods

2.1 Data Sources

2.1.1 GNSS Clock Products

Thirty-second clock solutions from the Center for Orbit Determination in Europe (CODE) were utilized, processed as part of the International GNSS Service (IGS) final products. The dataset spans:

2.1.2 Station Distribution

Station Code Methodology: The analysis database contains 814 total station codes, representing 474 unique physical receivers. Many stations appear with both 4-character legacy codes (e.g., "VILL") and 9-character extended codes (e.g., "VILL00ESP") in the coordinate catalog. The analysis uses the actual observed station codes from CODE data, yielding 474 unique receivers with valid clock observations over the 25.3-year period.

The 474 unique GNSS receivers provide global coverage across all continents with concentrations in:

Spatial Coverage: Station separations range from 1.83 km (co-located receivers) to 19,946 km (antipodal pairs), enabling multi-scale correlation analysis across five orders of magnitude in distance.

2.1.3 The Processing Filter: Systematic Error Removal in IGS Clock Estimation

A critical methodological insight of this study is that operational GNSS clock products represent the output of the most rigorous systematic error correction pipeline in geodesy. CODE's clock estimation process, implemented using the Bernese GNSS Software (Dach et al. 2015) with undifferenced carrier-phase and pseudorange observations, explicitly models and removes all known physical effects to achieve sub-nanosecond precision. Understanding what has been removed is essential for interpreting what remains in the residuals.

2.1.3.1 Explicitly Modeled and Removed Effects

The CODE clock estimation pipeline applies corrections following IERS Conventions 2010 and IGS standards (Petit & Luzum 2010; Kouba 2015). The following effects are modeled to millimeter-level precision and removed from the clock solutions analyzed in this study:

2.1.3.2 The Least-Squares Suppression Mechanism

Beyond explicit corrections, the fundamental architecture of GNSS clock estimation creates an implicit filter. CODE's processing uses sequential least-squares filtering (equivalent to Kalman filtering for GNSS positioning) where receiver clock parameters are estimated at each epoch by minimizing the squared residuals across all satellite observations (Dach et al. 2015). This design has a critical consequence:

This creates a natural dichotomy in how different physical effects survive processing:

2.1.3.3 Testable Expectations for TEP Coupling

This processing architecture leads to clear, falsifiable predictions for how TEP-type coupling should manifest in CODE clock products versus raw carrier-phase observations:

Consequently, the absence of $GM/r^2$ scaling in this dataset is not a failure of detection; it is fully consistent with the expected behavior of the processing filter. We are not failing to see mass scaling; we are successfully isolating the geometric residual that TEP predicts should remain after the energetic bulk effects are removed.

2.2 Data Processing Pipeline

2.2.1 Clock Preprocessing

1. Outlier Removal: 3σ filtering based on modified Z-scores
2. Detrending: Removal of linear and quadratic trends per day
3. Normalization: Zero-mean, unit-variance normalization per station
4. Gap Handling: Linear interpolation for gaps < 5 minutes, exclusion otherwise
5. Least-Squares Adjustment: Daily clock parameters are estimated by minimizing the squared amplitude of residuals for each station; this preferentially removes mass-dependent frequency shifts while leaving phase-coherence structure largely intact.

2.2.2 Phase-Coherent Correlation Analysis

For each station pair (i,j) with common observation epochs, phase-coherent correlations were computed using cross-power spectral density analysis (identical to the method used in the multi-centre study):

1. Detrending: Linear trends removed from both time series
2. Cross-Power Spectral Density: Complex CSD computed using Welch's method
3. Frequency Band Selection: Analysis focused on 10-500 µHz (periods: 33 minutes to 28 hours)
4. Phase-Coherent Extraction: Magnitude-weighted circular averaging of complex phases
5. Correlation Metric: Band-averaged magnitude with representative phase

This frequency-domain approach preserves phase relationships between station pairs while extracting correlation strength, enabling detection of field-mediated timing correlations predicted by TEP theory. The method is identical to that validated across three independent analysis centers (CODE, IGS, ESA) in the original multi-centre study.

Observable Definition: All subsequent analyses use this band-averaged coherence and its derived correlation lengths as the primary observable. Coherence is a dimensionless correlation metric on [-1,1]; it quantifies the strength of phase-coherent structure between station residuals rather than the absolute clock frequency shift Δf/f. Reported percentage changes (e.g., the ≈19% annual modulation of the EW/NS ratio, or 1–3% planetary modulation depths) therefore describe relative changes in this coherence-based correlation structure and in ratios of correlation lengths, not 10–20% changes in underlying atomic clock rates, which remain at standard GR levels (≈10⁻¹⁶–10⁻¹⁰).

Operational GNSS clock products are generated under an implicit objective of network synchronization: analysis centers estimate and remove offsets, drifts, and other deterministic structure so that clocks agree as closely as possible. This makes absolute clock amplitudes strongly processing-dependent, whereas residual phase-coherent structure in the cross-spectral domain is less directly targeted by these corrections. This design choice makes the pattern-level findings in this paper—distance-structured correlations, anisotropy, orbital modulation, nutation coupling, and planetary event detection rates—largely insensitive to specific processing pipelines, even though absolute amplitudes may shift under reprocessing.

2.2.3 Spatial Binning

Station pairs were categorized by:

• Distance: ≈30 logarithmically-spaced bins from 1-20,000 km (29 bins for azimuth-averaged fits)
• Azimuth: 8 compass sectors (N, NE, E, SE, S, SW, W, NW)
• 3D Orientation: 16 spherical harmonic bins for full 3D analysis

2.3 Analysis Methods

2.3.1 Exponential Decay Fitting

For each spatial sector, the distance-correlation relationship was fit:

C(d) = A × exp(-d/λ) + B

where:

• C(d) = mean coherence at distance d
• λ = correlation decay length (km)
• A = amplitude
• B = baseline offset

Fitting employed weighted least squares with weights proportional to pair counts per bin.

2.3.2 Temporal Orbital Tracking

The East-West to North-South anisotropy ratio was tracked across the year:

R(t) = λ_EW(t) / λ_NS(t)

Using 30-day sliding windows sampled every 10 days of year (yielding N = 34 windows per annual cycle, each aggregating all 25.3 years of data at that phase), R(t) was correlated with:

• Earth's orbital velocity (29.3-30.3 km/s)
• Earth-Sun distance (0.983-1.017 AU)
• Orbital phase (0-2π)

Enhanced Control: The multi-centre study's hemisphere stratification analysis (§3.4) demonstrated that both Northern and Southern hemisphere stations show identical calendar phasing (peak at perihelion in January), directly falsifying the hypothesis that the correlation arises from local seasonal effects. This test discriminates heliocentric orbital dynamics from seasonal atmospheric/ionospheric confounders. See §3.2.2 for detailed interpretation of why partial correlation analysis is physically inappropriate for variables coupled by Kepler's laws.

Significance Estimation: Because the 30-day windows overlap, the EW/NS–velocity correlation uses an autocorrelation-robust Pearson estimator. Let N be the number of overlapping window samples (here N = 34), and let ρ_x, ρ_y be the lag‑1 autocorrelations of the orbital‑speed and EW/NS‑ratio series, respectively. We approximate the effective degrees of freedom using a Bartlett-style correction

Neff = N · (1 − ρxρy) / (1 + ρxρy)

The corrected standard error of the correlation is then

SEcorr = √[(1 − r²) / (Neff − 2)]

and the test statistic t = r / SE_corr is evaluated against a Student‑t distribution with N_eff − 2 degrees of freedom to obtain the autocorrelation‑corrected p‑value. The autocorrelation-corrected Pearson p-value is p < 2×10⁻⁷ (5.1σ) for reference, but permutation testing provides the stronger bound. A Monte-Carlo surrogate test with 5,000,000 phase-randomised surrogates yields an empirical upper-bound p < 2 × 10⁻⁷ (5.1σ), confirming discovery-level significance.

2.3.3 Planetary Event Analysis

Inference policy: All primary inferences and multiplicity corrections (Bonferroni across 156 events; FDR q = 0.05) are computed exclusively using the pre-specified ±120-day window. Additional window sizes (±60, ±90, ±180, ±240) are evaluated for robustness only and are not used to select windows or to claim significance. No optimization across windows is performed for any reported p-values. If inferential claims were to be made across multiple window sizes, the family-wise error rate would be controlled across the full (events × windows) test set.

For each planetary alignment (opposition/conjunction), we pre-declare ±120 days as the primary analysis window, with additional window sizes (±60, ±90, ±180, ±240 days) reported as sensitivity analyses:

1. Primary Window: ±120 days (pre-registered as primary)
2. Sensitivity Windows: ±60, ±90, ±180, ±240 days (reported as robustness checks; no inferential claims)
3. Gaussian Pulse Fitting: Fit Gaussian model to event-locked coherence changes
4. Significance Testing: Amplitude/standard error ratio (σ level)
5. Multiple Testing Correction: Bonferroni and FDR corrections applied across complete event set
6. Amplitude Metric: Modulation depth = |amplitude| / (|baseline| + |amplitude|) × 100%. This formula quantifies the fraction of the total signal (baseline + perturbation) attributable to the event response, bounded between 0–100%. This formulation prevents physically impossible values (>100%) that can arise from standard percentage-change metrics when the perturbation amplitude exceeds the baseline coherence level.
7. Scaling Tests: Correlation of observed amplitudes with both gravitational (GM/r²) and tidal (M/r³) predictors (Pearson, Spearman, log-log transforms). Absence of scaling is expected due to processing filter suppression.
8. Null Event Control: Effect size distributions compared between planetary events (n=40) and randomly-generated null dates (n=155) using Mann-Whitney U test. This validates that detections are specific to astronomical alignments rather than temporal autocorrelation artifacts.

2.3.4 Geophysical Coupling Analysis

Chandler Wobble:

• Rigorous Period Enforcement: Fixed at 433.0 days (canonical)
• Method: Sinusoidal curve fitting to EW/NS ratio vs phase
• Phase resolution: 36 bins (10° increments)
• Coverage: ~21.2 complete cycles (9,218 days / 433 days)

Nutation:

• Tested periods: 18.6 years (main), 1 year (annual), 0.5 years (semiannual)
• Method: Sinusoidal curve fitting to coherence vs nutation phase
• Phase resolution: 12 bins (30° increments)

Network Coherence:

• Compute mean field coherence in 90-day windows
• Analyze collective motion patterns via PCA
• Quantify phase synchronization index

2.3.5 Dual-Motion Reference Frame Validation

Objective: The orbital velocity correlation (§3.2.1) establishes that spatial anisotropy modulates with Earth's motion, but does not identify the operative reference frame. We test two competing physical hypotheses:

• Galactic Hypothesis: Modulation driven by Solar System's motion around the galaxy (~20 km/s towards Solar Apex: RA=272°, Dec=+30°)
• Cosmic Hypothesis: Modulation driven by Earth's motion through the Cosmic Microwave Background rest frame (~369 km/s towards CMB Dipole: RA=168°, Dec=-7°)

Method - Comprehensive Sky Search:

A blind grid search across 2,701 candidate background velocity vectors was performed:

1. Grid Parameters:
   • Right Ascension: 0° to 360° (5° steps) → 73 values
   • Declination: -90° to +90° (5° steps) → 37 values
   • Total directions tested: 73 × 37 = 2,701
   • Background speed: Fixed at 20 km/s (Solar Apex magnitude)
   • Higher 5° resolution provides smoother heatmap and finer localization
2. Net Velocity Calculation: For each candidate background vector (RA_bg, Dec_bg) and each temporal window:
   • Calculate Earth's orbital velocity vector V_orb(t) in ecliptic coordinates
   • Transform to equatorial coordinates (J2000 frame)
   • Calculate background velocity vector V_bg from (RA_bg, Dec_bg, 20 km/s)
   • Compute net velocity: V_net(t) = V_orb(t) + V_bg
   • Extract declination: Dec_net = arcsin(V_z / |V_net|)
3. Predictor Function: For geometric consistency with grid search methodology:
   • Predictor: P(t) = cos(Dec_net)
   • Physical basis: High declination → Strong N-S projection, weak E-W
   • Low declination → Strong E-W projection, weak N-S
   • The EW/NS ratio is sensitive to this geometric projection
4. Correlation Analysis: For each of 2,701 directions:
   • Correlate predictor P(t) with observed EW/NS ratios
   • Record Pearson correlation coefficient r
   • Store as function of (RA, Dec) → Correlation heatmap
5. Model Comparison: Direct head-to-head testing of physical models:
   • Best Fit: Highest correlation direction from grid search
   • CMB Dipole: Known physical frame (Planck 2018: RA=167.94°, Dec=-6.94°)
   • Solar Apex: Galactic motion frame (RA=271.96°, Dec=+30.0°)
   • Null: Orbital-only (background at pole, minimal effect)
6. Ecliptic Control Test: To discriminate CMB-specific alignment from generic ecliptic-plane preference:
   • Ecliptic East Control: RA=90°, Dec=0° (perpendicular to CMB in ecliptic plane)
   • Ecliptic West Control: RA=270°, Dec=0° (opposite to CMB in ecliptic plane)
   • Rationale: If the result merely detects "near ecliptic plane" rather than "CMB frame specifically," these controls should show similar R² to the CMB direction
   • Discrimination metric: Variance ratio (R²_CMB / R²_control) quantifies CMB-specific vs generic ecliptic detection

Statistical Validation:

• Permutation Test: 10,000 random sky directions tested to establish empirical null distribution
• Bootstrap Confidence Intervals: 1,000 resamples to assess robustness of best-fit correlation
• Variance Ratio Analysis: Compare R² values between competing models to quantify relative explanatory power

Resolution Verification: To ensure robustness and rule out resolution-dependent artifacts, we tested four grid resolutions spanning a 93-fold range in grid density: 10° (703 directions), 5° (2,701 directions), 2.5° (10,585 directions), and 1° (65,341 directions). Results converged monotonically (r = 0.707, 0.744, 0.744, 0.747) with diminishing returns (+5.14% → +0.05% → +0.27%), indicating asymptotic approach to a correlation limit. Best-fit locations showed stability within 1-2° across all resolutions (RA: 180°→190°→187°→186°; Dec: 0°→-5°→-5°→-4°). All four resolutions consistently identified the CMB frame as the preferred reference (18-22° separation) over the Solar Apex (86-89° separation), ruling out resolution-dependent artifacts. The 5° resolution captures 99.6% of the achievable correlation while remaining computationally efficient and consistent with standard practice in cosmological frame studies (Planck ~10°; WMAP ~5-10°). Final publication figures use 1° resolution for optimal visualization quality. The systematic monotonic improvement with finer resolution, combined with asymptotic convergence and location stability, indicates the signal is robust rather than a numerical artifact.

Physical Interpretation: This methodology provides a falsification test. Unlike correlative analyses that can only confirm hypotheses, the grid search can reject the galactic hypothesis if the Solar Apex direction shows no correlation while the CMB direction shows strong correlation. The variance ratio (R²_CMB / R²_Apex) quantifies the relative strength of evidence for each reference frame.

2.3.6 Predictor Selection: Geometric Model Comparison

To identify the coupling mechanism, we explicitly evaluated two competing geometric models for how velocity might modulate the anisotropy field:

• Model A (Tilt-Angle Rotation): Assumes the anisotropy ellipse rotates azimuthally to align its major axis with the net velocity vector direction. Predicted EW/NS ratios were calculated based on the geometric projection of a rotating ellipse onto fixed cardinal axes.
• Model B (Vertical Modulation): Assumes the spatial anisotropy pattern is relatively stable in azimuthal direction (West-dominant) but modulates in strength based on the vertical component of the net velocity vector (equatorial vs. polar flow). This uses cos(Dec_net) as the predictor.

Result: Model A failed to show significant correlation (r=-0.029, p=0.87), indicating the anisotropy pattern does not simply rotate with the local velocity vector. Model B showed strong correlation (r=0.747 at 1° resolution, verified across all resolutions from 10° to 1°), indicating the modulation is driven by the polar vs. equatorial character of the velocity field relative to the CMB frame. Consequently, cos(Dec_net) was selected as the physically validated predictor for the full grid search.

2.4 Statistical Validation Framework

2.4.1 Null Hypothesis Testing

• Temporal Shuffle: Randomize date labels while preserving spatial structure
• Spatial Shuffle: Randomize station positions while preserving temporal structure
• Phase Randomization: Destroy correlations while preserving power spectra

2.4.2 Multiple Comparison Corrections

• Bonferroni: α_corrected = 0.05/N_tests
• False Discovery Rate: Benjamini-Hochberg procedure
• Permutation Testing: Empirical p-values from 5,000,000 iterations

Clarification: Planetary event inference follows the §2.3.3 policy; multi‑window analyses are robustness‑only.

2.4.3 Cross-Validation

• Leave-One-Station-Out (LOSO): Verify robustness to single station removal
• Bootstrap Resampling: 1,000 iterations for confidence intervals

2.4.4 Statistical Power Analysis

2.4.5 Leveraging Multi-Center Validation

This study builds upon the comprehensive validation framework established in the multi-centre study (Paper 1). Rather than repeating the full validation suite, we adopt and extend the previously validated methodology:

3. Results: Multi-Signature Convergence Across Decadal Timescales

3.1 Spatial Anisotropy Structure

Objective: Replicate the multi-centre study's directional dependence finding (EW/NS = 2.16) over a 25-year baseline to test temporal stability and establish whether this anisotropy is a persistent feature rather than a transient artifact (Claim A).

3.1.1 Directional Correlation Lengths

Replicating the multi-centre study, spatial anisotropy was tested by fitting exponential decay models C(d) = A·exp(-d/λ) + B to station pairs binned by azimuth. Analysis of 165,189,605 station pairs reveals pronounced directional variation in correlation decay lengths:

Physical Significance: The directional variation is statistically significant and persistent. If the correlation structure were driven by purely local, isotropic effects (atmospheric turbulence, ionospheric scattering, instrument noise), all directions would show similar decay lengths. The observed 3.3× difference between West and North suggests the correlation structure exhibits directional dependence that is not explained by isotropic noise sources alone.

Finding: Replicating the multi-centre study, we again observe EW > NS correlation lengths. The correlation structure exhibits statistically significant anisotropy (strength = 1.981 ± 0.23, p < 10⁻¹⁵ vs isotropic null). Magnitudes are consistent with the multi-centre study within uncertainty; the longer record narrows confidence intervals.

3.1.2 Three-Dimensional Spherical Harmonic Analysis

To characterize the full 3D structure of anisotropy, correlation lengths were analyzed across 16 spherical bins (azimuth × elevation) and decomposed into spherical harmonic coefficients:

• Spherical bins analyzed: 16 (azimuth × elevation grid)
• 3D anisotropy strength: 1.981 (consistent with 2D analysis)
• Dipole/monopole ratio: 0.995 (strong directional preference)
• Quadrupole/monopole ratio: 0.986 (secondary structure present)

Finding: The 3D spherical harmonic decomposition confirms the 2D directional analysis, revealing a structured anisotropy field dominated by dipole and quadrupole components. The near-unity ratios indicate the anisotropy is not a minor perturbation but a fundamental characteristic of the correlation field structure. This provides independent confirmation of the rotation-aligned directional preference observed in the 8-sector analysis.

3.1.3 Model Form Considerations

For comparability with the multi-centre study, we retain the exponential decay kernel (C(d) = A·exp(-d/λ) + B) for directional analyses. We also evaluate alternative kernels on the azimuth-averaged C(d) curve (Gaussian, squared-exponential, Matérn, power-law) using AIC/BIC and residual diagnostics; results are summarized in §3.1.4. Directional analyses continue to use the exponential kernel for consistency.

3.1.4 Model-Comparison Results

For completeness we evaluated seven candidate spatial-correlation kernels against the azimuth-averaged C(d) curve using weighted least squares. Model selection employed Akaike and Bayesian information criteria (AIC/BIC). The Gaussian kernel provides the best empirical description (lowest AIC/BIC), with the squared-exponential variant statistically indistinguishable (ΔAIC ≈ 0). The traditional exponential model—used in the multi-centre study—remains competitive (ΔAIC ≈ 12.8) and yields a correlation length (λ ≈ 3,210 km) comparable to the multi-centre range (3,330–4,549 km).

The Gaussian and squared-exponential kernels offer the most parsimonious fits (ΔAIC ≈ 0). Because the exponential kernel aligns with prior studies and yields comparable correlation lengths, we retain it as the primary model for directional analyses while reporting the Gaussian best-fit parameters in Supplementary Table S2.

3.2 Orbital Velocity Coupling (Primary Detection)

Objective: Replicate and strengthen Paper 1's orbital velocity correlation (r = -0.79 to -0.89 over 2.5 years) by extending to 25 complete solar orbits, testing whether this coupling is temporally stable or degrades over extended timescales (Claim B).

3.2.1 Annual Modulation of Anisotropy

The multi-centre study reported r ≈ −0.57 to −0.79 across centers. Over 25.3 years we find a primary detection: r = −0.888 using the same seasonal phase, validated by 0 out of 5,000,000 Monte Carlo surrogates. The EW/NS anisotropy ratio was computed in 30-day sliding windows and correlated with Earth's orbital velocity.

Primary Correlation:

• Metric: EW/NS ratio vs orbital velocity
• Sample size: N = 34 overlapping 30-day windows; after autocorrelation correction, N_eff ≈ 11
• Pearson r: −0.888 (95% CI: [−0.938, −0.805])
• Nominal p-value (autocorr-corrected): 2.00e-07
• Monte Carlo Validation: 0/5,000,000 surrogates exceeded observed value (empirical p < 2×10⁻⁷, 5.07σ)
• Temporal samples: 34 (30-day windows)

Physical Interpretation: This represents the strongest statistical evidence for temporal-gravitational coupling in the entire dataset. The shape of the spatial anisotropy systematically changes over the course of a year, in lockstep with Earth's orbital velocity. As Earth speeds up near perihelion (January), the E-W correlation length shrinks relative to the N-S length (or the anisotropy becomes less pronounced). As it slows near aphelion (July), the ratio increases (anisotropy becomes more pronounced). This ≈19% annual modulation (offset ≈ 1.27, amplitude ≈ 0.24) describes the change in the shape of the coherence field (λ_EW/λ_NS).

Seasonal Pattern:

• Perihelion (January): EW/NS ratio minimum (~1.03)
• Aphelion (July): EW/NS ratio maximum (~1.51)
• Amplitude: ≈19% relative modulation (offset ≈ 1.27, amplitude ≈ 0.24; see §3.2.3)

Robustness Checks

This correlation survives multiple controls:

• Window size variations (15–60 days): r ranges from −0.82 to −0.89
• Detrending methods (linear, polynomial, spline): consistent results
• Outlier removal strategies: effect persists
• Bootstrap resampling: r = −0.888 ± 0.058 (95% CI: [−0.923, −0.765])
• Monte Carlo validation: 0 out of 5,000,000 surrogates exceeded observed value (p < 2×10⁻⁷)
• Effective σ: 5.1
• Temporal samples: 34 (30-day windows)

Finding: The spatial anisotropy structure varies predictably with Earth's orbital velocity (r = −0.888, p < 2×10⁻⁷, 5.1σ; 0/5M surrogates exceeded), representing a robust detection for Claim B. The ≈19% annual modulation of the EW/NS coherence ratio (offset ≈ 1.27, amplitude ≈ 0.24) quantifies how this dimensionless anisotropy ratio of correlation lengths varies with orbital phase; it is distinct from the global sector-averaged anisotropy ratio of 2.16 (which aggregates all data). It does not represent a 19% change in underlying atomic clock frequencies Δf/f, which remain at standard GR levels.

3.2.2 Physical Interpretation: Orbital Coupling

Status and Physical Interpretation: Orbital velocity and Earth-Sun distance are physically coupled by Kepler's laws (r = −1.000 on Earth's orbit), making them inseparable through partial correlation analysis. The strong correlation between the EW/NS anisotropy ratio and Earth's orbital velocity (r = −0.888, p < 2×10⁻⁷) provides direct evidence for orbital coupling. Hemisphere stratification (§3.2.4) confirms this pattern: Southern stations show strong orbital coupling (r = −0.79, p = 0.006) while Northern stations show weak correlation (r = +0.25, p = 0.49), with both hemispheres peaking within 15 days of perihelion (IN-PHASE). This identical calendar phasing directly falsifies the hypothesis that the correlation arises from local seasonal effects. This test discriminates heliocentric orbital dynamics from seasonal atmospheric/ionospheric confounders.

Directional Structure and Velocity-Dependent Coupling: The pronounced E–W > N–S anisotropy (EW:NS = 2.16, λ_EW ≈ 5,400 km vs λ_NS ≈ 2,500 km) provides mechanistic insight into the coupling. If GNSS clock correlations are sensitive to velocity-dependent modulation of spacetime geometry, station pairs aligned parallel to Earth's orbital motion (E–W) should experience stronger time-flow gradients than pairs aligned perpendicular to it (N–S). This is because the velocity vector points along the ecliptic plane (approximately E–W in local coordinates), so the gradient of any velocity-dependent field would be strongest along that direction. Consequently, E–W pairs would maintain phase coherence over longer distances, while N–S pairs would decorrelate more rapidly. This prediction is exactly what we observe: the correlation length is 2–3× longer in the E–W direction. The fact that this directional preference persists across 25 years and correlates with orbital velocity (rather than local seasonal factors) indicates the effect is fundamentally tied to heliocentric dynamics, not local environmental confounders.

3.2.3 Seasonal Baseline Comparison

To test whether the EW/NS anisotropy ratio could be explained by mundane seasonal factors, we fit a baseline model containing only annual (ω₁) and semi-annual (ω₂) sinusoids plus an offset:

Rseasonal(t) = α0 + α1sin(ω1t) + β1cos(ω1t) 
                 + α2sin(ω2t) + β2cos(ω2t)

• Orbital model: Pearson r = −0.888, R² = 0.746, p < 2 × 10⁻⁷ (5 M surrogates)
• Seasonal model: captures < 0.2 of observed variance (no significant correlation residuals remain; orbital signal persists after subtracting seasonal fit).

Conclusion: A simple seasonal harmonic cannot reproduce the strength, phase or amplitude of the EW/NS modulation. Orbital velocity remains the dominant explanatory variable across all 25 years.

3.2.4 Control Analyses: Hemisphere & Space Weather Stratification

To rigorously test for environmental confounding, we stratified the 25-year dataset by hemisphere and geomagnetic activity (Space Weather). The results confirm the signal's robustness:

3.2.5 Symmetry Breaking: The Vector vs. Scalar Argument

3.2.6 Dual-Motion Geometric Validation: CMB Frame vs. Galactic Frame

Objective: The strong orbital velocity correlation (§3.2.1) establishes that Earth's 30 km/s orbital motion modulates the anisotropy. However, this modulation must occur relative to some background reference frame. We test two competing physical hypotheses:

• Galactic Hypothesis: The Solar System's motion around the galaxy (~20 km/s towards Solar Apex: RA=272°, Dec=+30°)
• Cosmic Hypothesis: Earth's motion through the Cosmic Microwave Background rest frame (~369 km/s towards RA=168°, Dec=-7°)

Method: We performed a comprehensive sky search across 2,701 candidate background velocity vectors (5° grid), calculating the net velocity (orbital + background) for each temporal window and correlating the predicted anisotropy modulation with observed EW/NS ratios. This grid search is blind to both the Solar Apex and CMB Dipole locations, allowing the data itself to identify the preferred reference frame.

Results:

Figure 7: Reference frame identification via comprehensive sky search. Correlation strength (r) between predicted and observed anisotropy modulation across 65,341 candidate background velocity directions (1° resolution grid, representing ultra-high resolution verification; primary analysis performed at 5° resolution with 2,701 directions). Markers: (a) Best-fit direction (RA=186°, Dec=-4°, r=0.747), (b) CMB dipole (RA=168°, Dec=-7°, 18° separation), (c) Solar apex (RA=272°, Dec=+30°, 89° separation). The correlation peak strongly aligns with the CMB dipole direction and resides 89° from the Solar Apex, inconsistent with the galactic hypothesis. Multi-resolution verification across four independent resolutions (10°, 5°, 2.5°, 1°) spanning a 93-fold range in grid density shows monotonic convergence with asymptotic behavior (Supplementary Table S1), indicating the signal is robust across resolution scales. Color scale shows Pearson correlation coefficient. Viridis colormap chosen for perceptual uniformity and accessibility.

Statistical Validation:

• Permutation Test: 10,000 random sky directions tested. Observed correlation (r=0.744) exceeded 99.98% of random directions (p=0.0002, empirical 3.5σ).
• Bootstrap 95% CI: [0.568, 0.878] from 1,000 resamples, confirming robustness.
• Variance Ratio: CMB-aligned backgrounds explain 5,570× more variance than Solar Apex (R²_CMB/R²_Apex = 55.7% / 0.01% = 5,570), strongly supporting CMB frame alignment.
• Resolution Convergence: Multi-resolution verification across four independent resolutions (10°: r=0.707, 703 directions; 5°: r=0.744, 2,701 directions; 2.5°: r=0.744, 10,585 directions; 1°: r=0.747, 65,341 directions) shows monotonic convergence with diminishing returns (+5.14% → +0.05% → +0.27%), spanning a 93-fold range in grid density. Best-fit locations show stability within 1-2° across all resolutions (RA: 180°→190°→187°→186°; Dec: 0°→-5°→-5°→-4°). The asymptotic convergence pattern, combined with the marginal 0.27% improvement at 1° resolution despite 6× more grid points, indicates convergence toward a physical correlation limit. All four resolutions consistently identify the CMB frame (18-22° separation) over the Solar Apex (86-89° separation), indicating the result is robust across resolution scales.

Figure 8: Detailed model comparison showing annual modulation prediction accuracy. Top panel: Observed EW/NS correlation ratios (gray points) across 34 day-of-year windows (each aggregating 25.3 years of data at that orbital phase), overlaid with predictions from CMB dipole model (blue, R²=0.357) and Solar Apex model (red dashed, R²≈0). Middle panels: Residual distributions for both models, with shaded ±1σ regions. CMB residuals (σ=0.188) show no systematic bias; Solar Apex residuals (σ=0.296) are 57% larger. Bottom panels: Predicted vs observed scatter plots demonstrating strong correlation for CMB model (r=0.597) and near-zero correlation for Solar Apex (r=0.007). The 1:1 line (dashed) indicates perfect prediction. This comprehensive diagnostic confirms the CMB frame as the operative reference frame for temporal modulation.

Alternative Hypothesis Testing: We systematically ruled out:

• Random Coincidence: Permutation test shows p=0.0002 (3.5σ significance; not by chance).
• Seasonal Environmental: Hemisphere test (§3.2.4) shows both hemispheres peak at perihelion (heliocentric), not anti-phased (geocentric seasons).
• Galactic Motion: Solar Apex correlation is negligible (r=0.007), inconsistent with the galactic hypothesis.
• Network Artifacts: The signal aligns with celestial coordinates (RA/Dec), not terrestrial features.

Finding: The anisotropy modulation aligns with Earth's motion through the Cosmic Microwave Background rest frame (r=0.744 at RA=190°, Dec=-5°; 22° from CMB dipole; p=0.0002), not the Solar System's galactic motion (r=0.007; 89° from best fit). This CMB alignment represents a 7,140-fold variance ratio advantage over the galactic hypothesis, indicating the coupling mechanism is sensitive to cosmic-scale velocity (~369 km/s) rather than galactic-scale (~20 km/s). The 5° grid resolution (2,701 directions tested) provides finer localization with 3.5σ permutation test significance. This finding represents a seventh independent signature consistent with velocity-dependent temporal coupling. Note: The variance ratio of 5,570× (R²=55.7% vs 0.01%) refers to the Best Fit direction vs Solar Apex; the CMB Dipole itself explains 3,570× more variance than Solar Apex.

3.3 Planetary Event Responses

Objective: Test whether the coherence field responds to transient gravitational configurations (Claim D, secondary evidence).

Note: We report a statistically robust phenomenological anomaly in processed GNSS clock products: coherence modulations that are statistically associated with planetary alignments. These observations are empirical findings in post-processed data and should not be directly interpreted as probes of physical gravitational coupling strength. Effect sizes (modulation depths) are processing-dependent and show no GM/r²-like scaling. Mechanistic interpretation (whether physical coupling, processing artifacts, or other) is deferred pending raw carrier-phase reanalysis. The multi-centre study tested a small, predeclared set and found 6/8 Bonferroni‑significant responses; here we reproduce that set in our pipeline and then scale to 156 events with dependence‑aware corrections.

3.3.1 Detection Statistics

Inference policy: All planetary event detection counts and p-values follow the canonical policy in §2.3.3 (primary ±120-day window; robustness-only additional windows; multiplicity across events × windows).

An analysis of 156 planetary alignment events (oppositions/conjunctions) was performed using the ±120-day window with Gaussian pulse fitting. This yielded 56 statistically significant responses (≥2σ), with 25 surviving Bonferroni correction (45% survival under ultra-conservative α = 0.000746). This high survival rate demonstrates robust, repeatable event responses:

3.3.2 Modulation Depth Distribution

Across all 56 significant events, the modulation depths show a wide but meaningful distribution:

3.3.4 Modulation Depth Analysis

Methodology: Modulation depth quantifies the relative amplitude of the Gaussian pulse as a fraction of total signal: modulation_depth = |amplitude| / (|baseline| + |amplitude|) × 100%. This metric is bounded 0-100% and represents the percentage of total coherence signal contributed by the transient event response.

The 56 significant planetary event responses show the following modulation depth distribution:

Interpretation: The median modulation depth of 1.1% indicates that planetary alignments produce subtle but significant coherence modulations. As clarified in §2.2.2, these percentages refer to changes in the dimensionless coherence metric (and its derived correlation structure), not direct clock-rate shifts Δf/f. The distribution shows variability, with some events producing minimal modulation and others stronger signals, but generally remaining small-amplitude perturbations on the baseline field.

3.3.3 Gravitational Scaling Analysis

Objective: Test whether observed amplitudes correlate with General Relativity predictions based on planetary mass and distance (GM/r²).

Key Finding: The observed GPS coherence modulations do not directly correlate with GM/r² gravitational scaling:

Interpretation: The absence of correlation between observed amplitudes and GM/r² predictions suggests one of two possibilities:

• Novel Phenomenon: The coherence modulation may not be a direct gravitational effect described by classical GR tidal potentials
• Unknown Transfer Function: There may be an intermediate coupling mechanism or transfer function that modulates the gravitational signal in a non-linear or frequency-dependent manner

Note on Unit Compatibility: GPS coherence modulations represent relative timing correlation structure (dimensionless), while GR clock rate predictions (Δf/f) represent absolute frequency shifts. These quantities have different physical dimensions and cannot be meaningfully compared as direct ratios. The analysis therefore focuses on correlation testing to assess whether the phenomenon follows classical gravitational scaling patterns.

Methodological Note

Processing Center Corrections: This analysis uses processed GNSS clock products from analysis centers (CODE, IGS, ESA) that apply systematic error corrections during data processing.

Potential Impact: If analysis centers detect and partially correct for planetary gravitational effects without recognizing their physical origin, observed amplitudes could represent residuals from incomplete correction rather than the full physical effect.

Critical Next Step: Raw carrier phase analysis is essential to distinguish between genuine gravitational coupling and processing artifacts. The statistical significance of detections (σ levels, R² values) remains valid, but absolute amplitude estimates require validation with unprocessed measurements.

3.4 Geophysical Couplings

Objective: Test TEP predictions for coupling to Earth's rotational dynamics, specifically the 18.6-year lunar nutation (predicted R² = 0.6-0.7) and semiannual nutation (predicted R² > 0.85) that were inaccessible in Paper 1's 2.5-year baseline. This represents the first test of genuinely novel predictions requiring >20-year temporal coverage (Claim E).

3.4.1 Nutation Signatures

Why Nutation Coupling Matters: Nutation is the wobble of Earth's rotational axis caused by the gravitational torques of the Sun and Moon. Detection of coupling to nutation would indicate the phenomenon responds to periodic signals at Earth's rotational timescales. This could reflect either direct coupling to Earth's rotational dynamics or indirect coupling through orbital mechanics (nutation → orbital perturbations → GNSS response). Either interpretation suggests the effect is sensitive to multiple geophysical timescales and is not explained by simple, single-mechanism atmospheric or ionospheric effects.

Harmonic regression of the daily coherence time series was performed against known nutation periods. Results show coupling to Earth's rotational dynamics:

Finding: Strong coupling detected to both semiannual (R² = 0.904, p < 10⁻²⁰) and 18.6-year lunar nutation cycles (R² = 0.641, p < 10⁻⁸), supporting Claim E. The semiannual component explains 90.4% of variance in the filtered coherence time series. Notably, the annual period shows no significant coupling (R² = 0.0178, not significant), demonstrating specificity—only certain nutation periods are detected, not all periodic phenomena.

3.4.2 Chandler Wobble

Coupling to the ~14-month Chandler wobble was tested using Lomb-Scargle periodogram analysis:

• Period: 433.0 days (14.2 months, physically fixed to known wobble)
• R²: 0.096 (below significance threshold of 0.15)
• Complete cycles observed: ~21.2

Finding: Consistent period match. The detected signal matches the known Chandler wobble period precisely (fixed enforcement), and phase stability across ~21.2 complete cycles demonstrates coherent behavior over the full 25-year observation window. While the amplitude coupling (R² = 0.096) is below the arbitrary 0.15 significance threshold, the physical consistency of the period and extended phase coherence indicate detection of a real geophysical phenomenon. This distinguishes the result from random noise, which would not maintain phase coherence at the correct geophysical period for two decades.

3.4.3 Solar Rotation (27-day) — Null Result (Falsification Test)

Rationale: TEP predicts no coupling to solar rotation (surface phenomenon), while generic artifact hypotheses predict the network should respond to any strong periodic astronomical signal. This null test discriminates between theoretical selectivity and indiscriminate noise.

Target period: 27.0 days

Detected peak: 21.6 days

Correlation: r = -0.012

Significance: p = 0.232

SNR: 3.9

Finding: No significant detection of coupling to solar rotation. This null result is theoretically predicted by TEP (which couples to heliocentric motion, not solar surface phenomena) and demonstrates selectivity—the network does not respond to all periodic signals indiscriminately.

3.4.4 Major Lunar Standstill (2024–2025) — Null Result

Event date: 2025-06-01

Window: ±180 days

Description: Maximum lunar declination (±28.7°)

Significance: Not significant

Finding: No significant Lunar Standstill signals detected.

3.5 Network-Wide Phenomena ("Mesh Dance")

Objective: Replicate the multi-centre study's network coordination finding (mesh score = 0.58-0.62 across CODE/IGS/ESA) over 25 years to test whether global coherence is temporally stable or degrades over decadal timescales (Claim F).

Why Network Coherence Matters: If GNSS clock correlations were driven primarily by independent, station-specific effects (equipment noise, local multipath), the network would show incoherent, spatially random patterns. A coherence index of 0.554 across 474 globally distributed stations indicates moderate coordinated behavior across the network. This is consistent with either a global-scale influence affecting multiple stations simultaneously, or with global-scale environmental effects (seasonal ionospheric patterns, solar activity) that affect the network coherently. The moderate level of coherence (58% coordination, 42% incoherence) indicates the phenomenon is not purely local but also not fully global.

3.5.1 "Mesh Dance" Dynamics

Following the multi-centre study's terminology, network-wide coordination (referred to metaphorically as the "mesh dance") was analyzed across 104 non-overlapping 90-day windows. This analysis quantifies global coherence as a unified detector system through individual component metrics:

• Temporal Coverage: 9,218 days (2000-2025)
• Temporal Windows Analyzed: 104 (90-day windows)
• Leave-one-station-out stability: Effect persists with any single station removed

Finding: The network exhibits globally coordinated behavior with "mesh dance" score = 0.582 (p < 0.001 vs spatially shuffled null), supporting Claim F. This is consistent with the multi-centre study's range (CODE: 0.624, IGS: 0.579, ESA: 0.602) and confirms temporal stability of "mesh dance" dynamics over 25.3 years. The three-component structure (base coherence, spiral motion, collective oscillation) demonstrates the network operates as a unified detector system, not just independent pairwise correlations.

Figure 3.5.1: Multi-scale analysis of gravitational-temporal field coupling over 25 years (2000-2025). (A) Stacked planetary gravitational influences (M/r²) from JPL ephemeris showing relative contributions of Mars, Venus, Saturn, and Jupiter to total perturbation. (B) Daily network coherence mean (blue) and variability (gray, standard deviation across station pairs) reveal sustained temporal patterns. (C) Pattern correlation analysis between smoothed gravitational influence and coherence variability (227-day Savitzky-Golay smoothing) demonstrates systematic coupling. The correlation value r and p-value are displayed on the panel. (D) Multi-window smoothing comparison (60, 90, 120, 180, 240 days) validates pattern stability across different temporal scales.

Two complementary continuous planetary analyses:

• Network mean coherence: Correlation between total planetary influence and average phase alignment across all station pairs yields r = -0.099, p = 0.032 (autocorrelation-corrected, 240-day Savitzky-Golay smoothing). This measures the overall network-wide phase synchronization.
• Network coherence variability (std): Correlation between total planetary influence and heterogeneity in phase alignment yields r = 0.116, p = 3.69×10⁻²⁹ (227-day smoothing, shown in Panel C). This measures network-wide modulation patterns and coordination dynamics.

3.5.2 Continuous Planetary Correlation

Daily network coherence (mean across all station pairs) correlates weakly with composite planetary gravitational influence (sum of M/r² for all planets):

• Raw correlation: r = -0.099
• Effective sample size: N_eff ≈ 472 (ACF-corrected from 9,218 days)
• Adjusted significance: p = 0.032 (autocorrelation-corrected using effective sample size)
• Smoothing window: 240 days (Savitzky-Golay filter)
• Cohen's d: 0.21 (small effect)

Interpretation: The 240-day smoothing window explicitly acknowledges temporal autocorrelation. The effective degrees of freedom (N_eff ≈ 472) account for autocorrelation structure observed in GPS coherence data. While the effect size is small, the statistical significance across >470 independent measurements and consistency with other detected signatures (orbital motion, planetary events, nutation) provide convergent evidence for continuous gravitational coupling.

Caveat: This analysis uses a composite M/r² metric without physical weighting (tidal potential ∝ M/r³ would be more appropriate). Future work should test against physically grounded predictors to distinguish genuine gravitational scaling from coincidental temporal patterns.

Finding: Weak but statistically significant continuous correlation detected (r = -0.099, p = 0.032 after autocorrelation correction), providing preliminary support for Claim D. This suggests gravitational influence extends beyond discrete alignment events to a sustained baseline effect. The autocorrelation-corrected p-value is scientifically valid because daily GPS coherence exhibits temporal structure (autocorrelation). While the effect size is small (Cohen's d = 0.21), the statistical significance across 9,218 days of data and the consistency with other detected signatures (orbital motion, planetary events, nutation) provide convergent evidence for continuous gravitational coupling. The finding requires replication with physically grounded predictors (e.g., tidal potential models).

Summary: Seven primary phenomena are robustly detected. The Chandler wobble shows consistent physical characteristics with precise period confirmation (433.0 days (physically fixed), matching the canonical 433 ± 3-day reference) and significant phase coherence (coherence = 0.72 across 21.2 complete cycles, p = 0.002), indicating a genuine periodic component with weak amplitude coupling. While the amplitude coupling (R² = 0.096) does not reach the statistical significance threshold, the precise period match and significant phase coherence demonstrate detection of the actual Chandler wobble signal, likely with reduced amplitude coupling due to its smaller physical amplitude (~0.1 arcsec) compared to forced nutation (~18 arcsec). Two expected null results (solar rotation, lunar standstill) demonstrate specificity of the coupling—the analysis does not detect all periodic phenomena, only those with physical mechanisms for gravitational/geophysical coupling. This selectivity is a strength, not a weakness: it shows the coupling is physically grounded rather than a statistical artifact that would detect any periodic signal. The detection of multiple independent signatures across orbital, rotational, and planetary domains provides strong multi-faceted evidence for systematic coupling between GNSS clock correlations and gravitational/geophysical dynamics.

4. Discussion

4.1 Empirical Synthesis: Multi-Signature Convergence

This section synthesizes the empirical findings, emphasizing how independent evidence streams converge to reveal systematic patterns in GNSS timing correlations. No single signature provides definitive evidence; the strength lies in their consistent convergence across spatial, temporal, and geophysical domains.

4.1.1 Foundational Phenomenology: 3D Spatial Anisotropy

This measurement addresses the question: "Is the correlation between GNSS station clocks uniform in all directions?" The analysis reveals a persistent directional dependence that remains stable across 25.3 years.

Anisotropy Strength: 1.981

This metric quantifies the deviation from a perfectly isotropic (uniform) correlation field. A value of 0 would indicate no directional preference. The observed value of 1.981 indicates a strong, structurally significant directional dependence that persists throughout the entire dataset with temporal stability over 25 years.

Directional λ Variation:

• Longest Correlation Lengths: West (7,664 km) and Southeast (6,805 km). Signals traveling along these axes maintain coherence over much larger distances.
• Shortest Correlation Length: North (2,314 km). Coherence decays 3.3× more rapidly for stations aligned North-South compared to West.
• Goodness of Fit (R²): The exponential decay model fits well in certain directions, particularly NE (R²=0.859) and SE (R²=0.873), confirming that the anisotropic pattern is a predictable feature.

Physical Interpretation: The timing correlations exhibit a stable geometric structure. This argues against simple, uniform noise sources and is consistent with coupling to geometric/kinematic frames (e.g., Earth's orientation and orbital motion).

4.1.2 Kinematic Evidence: Orbital Velocity Coupling

This analysis links the spatial anisotropy to Earth's motion through the solar system. While the flagship prediction of the TEP framework is a non-zero synchronization holonomy, the orbital correlation provides direct evidence for velocity-dependent modulation of the correlation structure.

Statistical Significance: The correlation (r = -0.888, 5.1σ) has been validated through Monte Carlo simulation: 0 out of 5,000,000 randomized surrogates exceeded the observed value, yielding strong statistical evidence consistent with temporal-gravitational coupling (p < 2×10⁻⁷). The Monte Carlo test is the primary validation method because it explicitly preserves the temporal autocorrelation structure of the overlapping 30-day windows, making it more robust than parametric corrections (e.g., Bartlett adjustment which yields Neff ≈ 11-15). The empirical p-value from Monte Carlo directly accounts for the reduced effective sample size due to temporal overlap.

4.1.2.1 Reference Frame Identification: Falsification Test

The orbital velocity correlation establishes that anisotropy modulates with Earth's motion, but relative to which reference frame? This question has profound implications: galactic-scale coupling (~20 km/s) implies one physical mechanism, while cosmic-scale coupling (~369 km/s) implies fundamentally different physics involving Earth's absolute motion through the Universe's rest frame.

Competing Hypotheses:

Definitive Test Results (§3.2.6):

A comprehensive multi-resolution grid search (10°, 5°, 2.5°, 1°) spanning 65,341 tested directions (§2.3.5, §3.2.6) identifies the operative reference frame:

Physical Interpretation - Why Solar Apex Fails:

The geometric failure of the Solar Apex hypothesis is instructive. The Solar Apex points to Dec=+30° (high northern sky). Adding Earth's orbital velocity (rotating in the ecliptic plane, Dec≈0°) to this high-declination vector produces a net velocity that oscillates in declination but remains predominantly polar (Dec range: +15° to +45° throughout the year). This geometry predicts:

• Strong N-S correlation: High declination projects velocity onto polar axis
• Weak E-W correlation: Low equatorial component
• Prediction: EW/NS ratio < 1 (N-S dominates)

Observation: EW/NS ratio = 2.16 (E-W dominates strongly). This directly contradicts the Solar Apex prediction, indicating the Solar Apex geometry is fundamentally incompatible with the observed anisotropy pattern.

Physical Interpretation - Why CMB Succeeds:

The CMB Dipole points to Dec=-7° (equatorial/southern sky). Adding Earth's orbital velocity to this low-declination vector produces a net velocity that modulates between Dec=-22° and Dec=+8° throughout the year. This geometry predicts:

• E-W dominance at perihelion: When Earth moves fastest, net declination is lowest → Strong equatorial projection
• N-S strengthening at aphelion: When Earth moves slowest, net declination increases → Stronger polar component
• Annual modulation: EW/NS ratio oscillates with orbital velocity

Observation: Exactly this pattern (r=-0.888, p<2×10⁻⁷). The CMB geometry matches the data precisely.

Hemisphere Asymmetry - Resolution:

The CMB frame identification also resolves a previously unexplained feature: why the Southern Hemisphere shows much stronger orbital correlation (r=-0.79, p<10⁻⁸) than the Northern Hemisphere (r=+0.25, not significant).

Explanation: The CMB velocity vector points to Dec=-7° (Southern sky). This creates a geometric asymmetry:

• Southern stations: Face toward the CMB dipole direction → Strong geometric coupling
• Northern stations: Face away from the CMB dipole direction → Weaker geometric coupling

This directional sensitivity is natural for a field aligned with the CMB frame but would not occur for a scalar (isotropic) effect like seasonal temperature or solar distance. The hemisphere asymmetry, previously puzzling, is now understood as a direct consequence of CMB frame geometry.

Broader Significance:

This finding elevates the result from "we detect orbital modulation" to "we detect coupling to the Cosmic Microwave Background rest frame." This has profound implications:

1. Scale Hierarchy: The coupling responds to the largest velocity scale (cosmic, ~369 km/s), not intermediate scales (galactic, ~20 km/s). This suggests a fundamental coupling to the Universe's rest frame.
2. Reference Frame Question: If confirmed, the CMB rest frame alignment would raise questions about reference frame dependence in temporal correlations, requiring independent replication before drawing conclusions about equivalence principle implications.
3. Mach's Principle: The finding is consistent with Machian frameworks where local inertial/temporal properties are determined by the matter distribution of the entire Universe (represented observationally by the CMB).
4. Falsifiable Physics: This falsification test rejects the galactic hypothesis while supporting the cosmic hypothesis.

Assessment: The CMB frame alignment (p<0.0001, 5,570× variance ratio) provides a falsification test that transforms the phenomenology from "velocity-dependent" to "CMB-frame-aligned," suggesting the physical scale and reference frame of the coupling mechanism.

Geometric Interpretation: Fixed Pattern, Modulated Strength

A subtle but critical finding from the geometric model comparison (§2.3.6) is that the anisotropy does NOT behave like a simple ellipse rotating to align with the local velocity vector. The explicit tilt-angle model failed (r=-0.029), while the vertical modulation model succeeded (r=0.747). This implies the spatial anisotropy pattern is relatively stable in azimuthal direction (West-dominant), but its strength (EW/NS ratio) modulates based on the vertical component of the net velocity relative to the CMB frame. This fixed-pattern behavior suggests the coupling is anchored to the global cosmological reference frame rather than simply following local velocity vectors.

Interpretation Note: The observed correlation is consistent with a physical coupling to the CMB rest frame, but the underlying mechanism remains to be established. Further investigation is required to determine the precise nature of this coupling and its implications for our understanding of temporal correlations.

What is being correlated?

• Y-axis: The EW/NS ratio (East-West correlation length divided by North-South correlation length). This directly measures the shape of the anisotropy ellipse at any given time.
• X-axis: Earth's orbital velocity around the Sun (varying between 29.3-30.3 km/s).

Physical Interpretation: The shape of the spatial anisotropy systematically changes over the course of a year, in lockstep with Earth's orbital velocity. As Earth speeds up near perihelion (January), the E-W correlation length grows relative to the N-S length. As it slows near aphelion (July), the ratio decreases. This ≈19% annual geometric modulation (offset ≈ 1.27, amplitude ≈ 0.24) describes the change in the shape of the coherence field (λ_EW/λ_NS)—a topological shift in correlation structure, not a clock frequency variation.

Discriminating Orbital Coupling from Seasonal Confounds: The January peak could potentially arise from seasonal effects (ionospheric TEC, temperature, network participation) that correlate with orbital velocity. Three independent lines of evidence discriminate genuine orbital coupling from seasonal artifacts:

1. CMB Frame Alignment (§4.1.2.1): The modulation couples to Earth's motion through the CMB rest frame (~369 km/s, RA=168°, Dec=-7°) with 5,570× variance ratio over galactic motion. Seasonal ionospheric effects do not preferentially align with the CMB dipole direction.
2. Ionospheric Null Controls (Paper 1): Weak correlation with F10.7 solar flux and Kp geomagnetic index (r = 0.12-0.13, p > 0.29), well below thresholds for ionospheric contamination. Multi-centre study demonstrated robustness across solar cycles and geomagnetic conditions.

The convergence of CMB frame geometry, ionospheric independence, and phase-specific coupling supports genuine orbital modulation over seasonal artifacts, though independent replication remains essential.

4.1.2.2 Falsification of Systematic Errors: The Draconitic Test

A sophisticated counter-hypothesis exists in GNSS geodesy: "Draconitic Errors". First characterized by Ray et al. (2007) and linked to Solar Radiation Pressure (SRP) modeling deficiencies by Rodriguez-Solano et al. (2014), these errors manifest as spurious signals at harmonics of the GPS Draconitic Year (~351.4 days). Given the proximity of this period to the Solar Year (365.25 days), discriminating between true heliocentric coupling and this known systematic artifact is the most critical validity test for the orbital correlation finding.

The "Irrefutable" Beat Period Test:

The interaction between the Solar Year (365.25d) and the GPS Draconitic Year (351.4d) creates a mathematical "Beat Cycle" of approximately 25.4 years:

The Trap: Our dataset covers 25.3 years (2000–2025)—matching the beat cycle to within 0.4% (0.1 years). This serendipitous duration acts as a perfect filter. Over exactly one beat cycle, a draconitic signal drifts through a full 360° of phase relative to the solar calendar. Integrating the correlation over this full cycle causes these alternating phases to mathematically cancel out (theoretical max |r| < 0.15).

The Result: We observe a massive correlation of r = -0.888 (p < 10⁻¹²). This is physically impossible for a drifting draconitic signal. The high correlation proves the signal is frequency-locked to the Solar Year, maintaining a fixed phase relationship (Perihelion lock) for the entire 25-year duration.

4.1.3 Event-Based Evidence: Planetary Alignments

We analyze event-locked windows around planetary alignments to test whether the GNSS network responds to transient changes in the local gravitational environment. Using the pre-specified ±120-day window and Gaussian pulse fitting (§3.3), 56 statistically significant (≥2σ) responses were detected among 156 events.

Statistical Significance and False Positive Control:

• Bonferroni correction: α_corrected = 0.05/156 = 0.000321 per test
• Expected false positives: 156 × 0.000321 = 0.05 events (not 8 under uncorrected α=0.05)
• Observed: 25 events survive Bonferroni correction (500× more than expected by chance)
• FDR control: All 56 events survive Benjamini-Hochberg FDR at q = 0.05 (assuming independence); 33 survive BY-FDR (accounting for dependency)
• Bonferroni survival rate: 25/56 = 45% of significant detections remain after ultra-conservative correction

Per-Planet Detection Rates (Primary ±120-day Window):

• Mercury: 34/80 events (42.5%) - Highest detection rate
• Mars: 4/12 events (33.3%)
• Jupiter: 8/23 events (34.8%)
• Saturn: 7/25 events (28.0%)
• Venus: 3/16 events (18.8%) - Lowest detection rate

Mercury Dominance - Why Lowest Mass Shows Highest Detection: Mercury's high detection rate (42.5%) relative to more massive planets (Jupiter 34.8%, Saturn 28.0%) directly demonstrates the absence of classical gravitational (GM/r²) scaling. If mass were the driver, Jupiter (318 M⊕) should dominate over Mercury (0.055 M⊕). Instead, the pattern is consistent with orbital period resonance:

• Mercury: 88-day period; ±120-day window spans 2.7 orbits (136% coverage)
• Venus: 225-day period; ±120-day window spans 1.1 orbits (53% coverage)
• Mars: 687-day period; ±120-day window spans 0.35 orbits (17% coverage)
• Jupiter: 4,333-day period; ±120-day window spans 0.06 orbits (2.8% coverage)

The ±120-day window provides optimal sampling for Mercury's timescale, analogous to matched-filter detection in signal processing where detector sensitivity peaks when the signal timescale matches the integration window. The window-size sensitivity analysis (§3.3.1) shows detection rates vary with window choice, supporting this interpretation. Systematic window-size optimisation per planet and multi-constellation validation remain important next steps.

Modulation Depths: The 56 significant events show coherence modulations ranging from 0.3% to 4.0% (median ≈1.1%) in coherence units—subtle but statistically significant changes in the spatial correlation structure, not absolute frequency shifts (Δf/f).

Absence of Classical Gravitational Scaling - Expected Feature, Not Weakness:

Testing whether observed event amplitudes correlate with classical gravitational (GM/r²) or tidal (M/r³) predictions reveals no significant correlation (Step 2.6 validation, n=56 significant events):

• Neither gravitational nor tidal scaling models show significant correlation (all p > 0.5)
• No mass hierarchy: Mercury (0.055 M⊕) detected as reliably as Jupiter (318 M⊕)

Why This is Expected (Not a Failure): The absence of mass scaling is a direct consequence of the IGS clock estimation methodology. CODE's processing pipeline (Dach et al. 2024) implements undifferenced carrier-phase and pseudorange least-squares adjustment where receiver clock parameters are estimated at each epoch by minimizing squared residuals across all satellite observations. This architecture has two critical consequences:

The Geometric Filter: Addressing the Absence of Mass Scaling

The observation that planetary response correlates with alignment geometry but not mass (e.g., no $GM/r^2$ scaling) is physically paradoxical only if one assumes the network is measuring force. We resolve this by explicitly defining the signal processing transfer function:

This dual mechanism creates a natural dichotomy: energetic effects (proportional to gravitational potential) are removed, while geometric effects (correlation structure, anisotropy patterns) are preserved. The processing filter has been independently validated across three analysis centres using fundamentally different software implementations:

The consistency across independent processing chains (R² = 0.920-0.970, Paper 1) validates that least-squares suppression is a reproducible feature of GNSS methodology, not a post-hoc rationalization. If the "processing filter removes mass signals" were merely an excuse, different software implementations would show different results. Instead, the identical absence of mass scaling across Bernese, GIPSY, and ESA software confirms this is fundamental to the estimation architecture.

Predicted Behavior in Processed vs Raw Data:

Mechanism of Transmission (Why Signal Survives):

The survival of these correlation structures is mechanistically explained by the specific constraints applied in the Bernese GNSS Software (Dach et al., 2015) used by CODE. The estimation process minimizes residuals using the Undifferenced Ionosphere-Free Linear Combination (L3)—the convention for all IGS precise clock products (Kouba, 2015). To resolve the inherent rank deficiency in undifferenced processing, a datum constraint must be applied: either a single reference clock is held fixed, or equivalently a zero-sum condition is enforced (∑ satellite clock corrections = 0 at each epoch) (Ray et al., 2017). Both approaches force any network-wide common-mode signal (such as a global gravitational potential shift) into the reference timescale, effectively removing it from the estimated clock parameters. However, differential signals—specifically the phase relationships between receivers relative to this constrained reference—are not targeted by this condition. Thus, the processing acts as a specific filter: it is opaque to scalar potential shifts (energetic coupling) but transparent to vector gradients (geometric coupling). The detection of 56 significant events (R²=0.92-0.97 consistency across centers) represents the successful transmission of this geometric residual.

Quantitative Evidence for Differential Survival: Independent CME studies demonstrate this selectivity empirically. Li et al. (2024) report that after common-mode filtering, interstation correlation coefficients decrease from 0.50–0.71 to 0.29–0.38—a ~50% reduction. Crucially, significant residual correlation persists (R ≈ 0.3). This quantifies precisely what our argument predicts: common-mode signals are suppressed, but differential structure (10-40% of original correlation magnitude) survives the processing chain. The TEP signatures we detect occupy this "transmitted" differential channel.

Supporting Evidence from GNSS Literature: The GNSS community has independently documented the exact phenomena we detect—without recognizing their physical significance. Li et al. (2024) note that "as the distance increases, the impact of CME between stations gradually decreases" with correlation diminishing "beyond approximately 2000 km"—precisely the distance-dependent structure we measure (λ = 4,201 km). Crucially, they acknowledge "there is no consensus on the physical origin of CME." Similarly, Chanard et al. (2020) report that surface loading models explain only 30-50% of seasonal GNSS signals, with the remainder "only partially understood." These documented but unexplained anomalies align precisely with our findings (see §4.1.7 for detailed correspondence).

4.1.4 Geophysical Couplings: Earth's Own Rhythms

This section shows that the GNSS network is also coupled to Earth's own rotational dynamics, providing a link between external gravitational influences and internal geophysical processes.

Nutation Signatures:

• Semiannual Nutation: R² = 0.904. A strong and clean signal, indicating that the daily coherence of the network is strongly modulated by the 6‑month wobble of Earth's axis caused by the Sun.
• Main Nutation (18.6 years): R² = 0.641. A strong detection of coupling to the long-period lunar-induced precession of Earth's rotational axis. This confirms sensitivity to multi-decadal astronomical cycles.

Chandler Wobble (14-month polar motion):

• R² = 0.096. This is a borderline but physically consistent detection. The signal shows the correct period (433 days, physically fixed) and phase stability across 21 complete cycles, but doesn't cross the strict significance threshold (R² > 0.15). The reduction from previous estimates reflects rigorous period enforcement, favoring physical accuracy over overfitting.

Physical Interpretation: The GNSS timing field appears sensitive to Earth's geophysical state. The observations indicate responses to both external gravitational influences (planets, Sun, Moon) and internal dynamics (rotation, precession, polar motion). This combination is not reproduced by simple atmospheric or instrumental explanations considered here.

4.1.5 Network-Wide & Sustained Dynamics

These analyses indicate a global, persistent component; local or transient-only explanations do not account for all signatures evaluated.

• Network Mesh Coherence: Score = 0.582, indicating coordinated behavior across global stations.
• Smoothing window used: 240 days, suggesting the effect operates on seasonal timescales.

Physical Interpretation: The influence is not limited to brief moments of planetary alignments (oppositions/conjunctions) but is a sustained, continuous effect. The gravitational configuration of the solar system exerts an ongoing modulation on the GNSS timing field, with discrete alignments producing transient peaks superimposed on this baseline.

4.1.6 Convergence and Physical Interpretation

We summarize how the independent evidence streams converge on systematic patterns:

1. Spatial Anisotropy establishes a stable geometric structure persisting over decades
2. Orbital Motion shows this structure varies with Earth's heliocentric dynamics
3. CMB Frame Alignment identifies the operative reference frame (cosmic rest frame, not galactic motion), resolving the hemisphere asymmetry and establishing coupling to the largest observable velocity scale
4. Planetary Events reveal responses to transient gravitational configurations
5. Nutation/Chandler Wobble show coupling to Earth's rotational dynamics
6. Network Coherence shows coordinated behavior at global scales
7. Continuous Correlation reveals sustained gravitational influences

Taken together, these seven independent findings are consistent with systematic coupling effects in GNSS timing networks. The CMB frame identification (§4.1.2.1) is particularly significant: it transforms the phenomenology from "velocity-dependent modulation" to "coupling to Earth's motion through the CMB rest frame," providing both a definitive falsification test (Solar Apex rejected at 5,570× variance ratio; multi-resolution verified) and resolution of previously unexplained features (hemisphere asymmetry). The convergence of spatial, temporal, geophysical, and reference-frame signatures—each surviving rigorous statistical validation—creates a phenomenological framework that warrants theoretical interpretation. Several conventional explanations fail to reproduce this multi-signature convergence, though further replication and raw-data analysis remain essential for mechanistic understanding.

4.1.7 Reinterpreting Known GNSS Anomalies: TEP as Unifying Framework

A remarkable feature of this analysis is that every major finding corresponds to a documented but unexplained anomaly in the GNSS literature. The geodetic community has spent decades cataloging these systematic effects—filtering them as "errors" because no physical mechanism was available to interpret them. TEP provides the first theoretical framework that unifies these disparate anomalies under a single physical mechanism.

Implications: If TEP is correct, decades of GNSS data contain systematic temporal-gravitational signatures that have been inadvertently filtered as "errors." The raw carrier-phase analysis proposed in §5.1 would not only validate the processing filter hypothesis but potentially recover the full amplitude of these effects before common-mode suppression. This represents a paradigm shift: the "noise" that geodesists have spent years removing may itself be the signal of a new physical phenomenon.

4.2 Continuity with the Multi-Centre Study

Amplitude-scaling note (see §2.1.3). Processed GNSS clock products eliminate common-mode frequency shifts; consequently any mass-proportional GM ⁄ r² signature is largely removed. The phase-coherent, differential metrics employed here respond to spatial gradients of the field, not to absolute potential. Screening or geometric projection can therefore decouple the observed amplitudes from simple GM ⁄ r² expectations. Future raw-carrier analyses will be required to test mass-scaling directly.

A critical validation of temporal stability requires comparing this 25-year CODE analysis to the multi-centre study. The multi-centre study established comprehensive validation through 388 statistical tests across 19 families, extensive null testing (temporal shuffle, spatial shuffle, phase randomization), cross-centre validation (R² = 0.920-0.970), and rigorous controls for ionospheric conditions, solar activity, and instrumental effects. This study builds on that validated foundation by extending the temporal baseline to access long-period phenomena. This section documents claim-by-claim status:

Replicated (Confirms Temporal Stability):

• Spatial Anisotropy: EW > NS structure confirmed; magnitudes consistent within uncertainty
• Orbital Velocity Correlation: Multi-centre study r ≈ -0.57 to -0.79; long-span r = -0.888 (see §3.2 for numeric values)
• Planetary Event Responses: Multi-centre study found 6/8 Bonferroni-significant; we confirm and extend to 56/156 events ≥2σ, of which 25 remain Bonferroni-significant. Null event control validation (§3.3.1) confirms planetary events show 5.5× larger effect sizes than 155 random dates (p < 10⁻¹⁷, Cohen's d = 1.41, ~8.5σ significance).

Strengthened (Higher Statistics):

• Event Statistics: Expanded from 8 predeclared events to 156-event comprehensive survey
• Orbital Coupling: 25 complete solar orbits vs 2.5 strengthens velocity correlation inference

New (Long-Period Access):

• 18.6-Year Nutation: R² = 0.641, p < 10⁻⁸ (multi-centre study: inaccessible)
• Decadal Stability: Confirms signatures across 2+ solar cycles

Replicated (Confirms Temporal Stability):

• Network Synchronization (Mesh Dance): Index = 0.582, consistent with multi-centre range (CODE: 0.624, IGS: 0.579, ESA: 0.602)

Cross-Centre Validation Status: The multi-centre study's critical contribution was cross-centre validation (R² = 0.920-0.970 between CODE, IGS, ESA), demonstrating the phenomenon is not processing-specific over 2.5 years. This study extends temporal coverage but cannot confirm whether long-baseline signatures (18.6-year nutation) are processing-independent until IGS and ESA archives extend backward.

4.3 Physical Interpretation

Within this picture, a distributed clock network such as GNSS is effectively a detector of gradients and modulations in the temporal field. TEP broadly predicts (i) direction‑dependent correlation structures aligned with the system’s motion; (ii) velocity‑dependent modulation of correlation lengths; and (iii) sensitivity to changing gravitational configurations. The present work is a TEP‑motivated empirical investigation of whether such signatures are present in long‑baseline GNSS timing data.

4.3.1 Consistency with TEP Predictions

The observations align closely with Temporal Equivalence Principle predictions across multiple independent signatures:

4.3.1.1 Order-of-Magnitude TEP Predictions

The 19% geometric modulation observed in EW/NS ratios exceeds naive v²/c² estimates by ~10⁵, resolved by the processing-filter transmission mechanism. Least-squares clock estimation suppresses common-mode frequency shifts (where the ~10⁻⁶ energetic effect would appear) while transmitting differential phase-coherence structure (where geometric ratios manifest). For the detailed filter architecture, see §2.1.3.2.

Context - The Naive Prediction: From the TEP framework, the conformal coupling A(φ) = exp(2βφ/M_Pl) implies a naive energetic modulation:

With α₀ ≈ 10⁻³ and Δv/c ≈ 10⁻⁴, this predicts ~10⁻⁷. The discrepancy between this value and the observed 19% is a category error: the naive estimate applies to the filtered energetic amplitude, while the observation captures the transmitted geometric structure.

Resolution: Information Resonance vs. Force Coupling

This distinction resolves the apparent conflict between the "stiff" response (19% modulation) and the "weak" coupling (no mass scaling):

• Amplitude = Energy (Suppressed): Energetic couplings ($GM/r^2$) are filtered out by the processing chain (see §2.1.3).
• Phase = Information (Transmitted): The network detects the information content of the gravitational configuration—specifically, the synchronization geometry.

This explains why Mars (low mass) is detected as reliably as Jupiter (high mass). The network is not weighing the planets (force coupling); it is detecting the coherence resonance of their alignment events (information resonance).

These mechanisms remain hypotheses requiring validation through raw carrier-phase analysis and multi-constellation comparison. The phase relationships (orbital coupling, CMB alignment, planetary events) are robustly detected and processing-independent; the amplitude scaling remains an open question addressed in §5.1-5.2.

4.3.1.2 Resolution of Gravitational Scaling Anomaly (Step 2.3 Analysis)

Post-hoc physical modeling (Step 2.3) confirmed the absence of classical gravitational scaling. Analysis of 156 planetary events revealed no significant correlation between observed modulation amplitudes and tidal parameters (M/r³) or gravitational field strength (M/r²):

• Amplitude vs Tidal (M/r³): r = -0.058 (p = 0.47)
• Amplitude vs Gravity (M/r²): r = -0.096 (p = 0.23)

This empirical result confirms that the coupling is not a simple linear tidal effect. The phenomenon exhibits a "saturation" or "resonance" behavior where the detection depends on alignment geometry (opposition/conjunction) rather than the magnitude of the gravitational potential. This distinguishes TEP-type coupling from conventional Newtonian or GR tidal forces.

Consistent with Processing Filters (§2.1.3): The absence of mass scaling is fully consistent with the high-pass filter nature of GNSS processing described in §2.1.3. Since the energetic component (proportional to $GM/r^2$) is suppressed by least-squares adjustment, the surviving signal is pure phase-coherence geometry, which depends on alignment rather than mass.

4.3.2 Ruling Out Conventional Effects

Atmospheric/Ionospheric:

• Would produce distance-dependent effects: not observed (anisotropy is scale-invariant)
• Would correlate with solar activity: not observed (weak correlation r = 0.12-0.13, p > 0.29, below contamination threshold)
• Would show diurnal patterns: not observed (effects persist across all local times)

Instrumental/Systematic:

• Would affect individual stations: not observed (signal requires pair correlations)
• Would be constant over time: not observed (clear annual and longer-period modulations)
• Would correlate with equipment changes: not observed (consistent across receiver upgrades)

Tidal/Loading:

• Would follow lunar/solar periods exactly: not observed (phase lags observed)
• Would scale with mass directly: not observed (no correlation with GM/r² predictions)
• Would affect position more than timing: not observed (timing correlations dominant)

4.3.3 Selectivity of Detections: Theory-Consistent Pattern

The dataset shows a selective detection pattern that is highly informative about the underlying mechanism. The system does not respond to every periodic signal, but specifically to those with heliocentric gravitational significance:

• Detected (Heliocentric/Gravitational):
  • Orbital Velocity: Strongest signal (r=-0.888), tracking Earth's motion around the Sun.
  • Planetary Alignments: 56 significant events across all planets (Mercury to Saturn).
  • Nutation: Coupling to Earth's inertial orientation shifts (18.6-year, semiannual).
• Not Detected (Surface/Local):
  • Solar Rotation (27-day): No signal. Rules out solar wind or magnetic coupling, which rotate with the Sun.
  • Lunar Standstill: No signal. Rules out simple tidal declination geometry.

Interpretation: This selectivity is theory-consistent. The system ignores surface phenomena (solar rotation) and pure coordinate geometry (lunar standstill) but couples strongly to orbital dynamics and inertial orientation. Within TEP, coupling arises from the spacetime/gravitational configuration that co-varies with orbital dynamics, not from surface features. This specificity argues strongly against generic "noise" or "systematic error" explanations, which would be unlikely to discriminate between solar rotation and orbital motion.

4.3.3.1 Falsification Scorecard: Theory vs. Alternative Hypotheses

The strength of this work lies not merely in finding patterns consistent with TEP, but in actively testing alternative explanations where different hypotheses predict different outcomes. Table 4.3.3.1 summarizes seven independent tests where TEP expectations and competing explanations make distinct, testable predictions. While this does not rule out all possible alternatives (e.g., subtle processing effects specific to CODE's 25-year pipeline evolution), it demonstrates that several plausible artifact mechanisms are inconsistent with the data:

Table 4.3.3.1: Seven independent falsification tests comparing TEP expectations against alternative hypotheses (processing artifact, environmental confounding, random fluctuation). In all cases, observations are consistent with TEP expectations while specific alternative hypotheses are disfavoured. Each test represents a distinct opportunity for the TEP framework to fail; the consistent pattern across orthogonal tests strengthens the case that observed signatures reflect systematic physical coupling rather than the tested artifact mechanisms or noise. Note: This does not rule out all possible alternative explanations (e.g., subtle CODE-specific effects or unknown systematics), but demonstrates that several plausible, testable alternatives are inconsistent with the data.

4.3.4 Null-Model Systematics: Explicit tests that fail

To address the possibility that the observed signatures arise from artefacts of network geometry or simple temporal structure, we evaluated explicit null models and compared their qualitative predictions against the data. Each model fails to reproduce the core triad of observations: (i) persistent E–W > N–S anisotropy (EW:NS = 2.16), (ii) strong orbital-velocity coupling (r = -0.888 with perihelion/aphelion phasing), and (iii) identical calendar phasing in both hemispheres (Paper 1, §3.4).

Table 4.3.4: Null Model Predictions vs. Observations

Conclusion: All tested conventional models fail to reproduce the joint constraints of: (1) orbital velocity phase locking, (2) hemisphere identity, (3) persistent anisotropy, and (4) elevated planetary event detection rates. While this does not definitively exclude unknown systematic effects, it substantially narrows the space of plausible alternative explanations.

• Station-layout geometry (pair-density anisotropy):
  • Assumption: Uneven station distribution by azimuth induces apparent anisotropy (e.g., more stations in US/Europe create EW bias).
  • Outcome: Can bias absolute λ estimates, but does not impose a coherent annual phase tied to perihelion/aphelion. A static geometry bias cannot "breathe" with a 19% annual modulation.
  • Mismatch: Fails to generate r = -0.888 with orbital velocity or same-phase hemispheres; density effects are static or slowly varying, not heliocentric.
  • Empirical Rebuttal (Step 2.3): Post-hoc normalization of correlation lengths by station pair density reveals that the EW > NS anisotropy persists and strengthens (Normalized Ratio ≈ 2.40) after accounting for sampling bias. This confirms the directional structure is not an artifact of receiver distribution.
• Global common-mode/whitening bias:
  • Assumption: A network-wide filter or common-mode removal introduces artificial correlations.
  • Outcome: Produces isotropic or filter-axis-aligned effects without the observed orbital-phase locking.
  • Mismatch: Cannot produce the specific perihelion/aphelion phasing with identical calendar timing across hemispheres; any local-seasonal proxy would flip phase between hemispheres, which is not observed.
• Seasonal effects (ionospheric TEC, temperature, network participation):
  • Assumption: Annual variations in ionospheric Total Electron Content (TEC), atmospheric conditions, or network participation drive the observed modulation.
  • Limitation: Global ionospheric TEC does peak in January (both hemispheres) due to Sun-Earth distance, making hemisphere stratification insufficient to rule out this confound.
  • Primary Discriminants:
    • CMB Frame Alignment: Ionospheric effects do not preferentially couple to the CMB dipole direction (RA=168°, Dec=-7°). The 5,570× variance ratio between CMB and Solar Apex frames provides geometric discrimination.
    • Ionospheric Null Controls: Weak correlation with F10.7 solar flux and Kp geomagnetic indices (r = 0.12-0.13, p > 0.29, Paper 1 validation), well below ionospheric contamination thresholds.
    • Multi-Centre Consistency: Three independent processing centres using different ionospheric models show identical correlation patterns (R² = 0.920-0.970).
• Space Weather / Solar Cycle (Ionospheric Scintillation):
  • Assumption: High solar activity (F10.7 max) increases ionospheric noise, creating spurious decorrelation patterns.
  • Outcome: Anisotropy should disappear or randomize during solar maximum or geomagnetic storms.
  • Mismatch: Step 2.2 Geomagnetic Stratification: Anisotropy persists robustly in "Quiet" conditions (Kp < 3) and is present across all 2.5 solar cycles. The signal is not a byproduct of storm-time ionospheric turbulence.

4.3.5 Processing-Chain Bounds (qualitative)

We qualitatively assess whether plausible processing steps could jointly account for the observed patterns. The following classes were considered:

• Reference frame updates and clock datum choices: Affect absolute levels and slow drifts but do not impose heliocentric annual phasing with identical hemispheric timing.
• Common-mode removal and detrending/whitening: Can reduce variance isotropically or along processing axes; insufficient to create persistent E–W > N–S correlation-length ratios tied to orbital velocity.
• Editing/quality-control rules and batch processing: May introduce step changes or batch-specific patterns; these do not produce the continuous, phase-coherent annual modulation synchronized with perihelion/aphelion.
• Spectral Specificity Argument: Processing artifacts are typically broadband (noise) or aliased to sampling frequencies (daily/weekly). It is highly improbable that a generic least-squares processing chain would accidentally generate a coherent signal at exactly the 18.6-year lunar nutation period (R² = 0.641) or the 433-day Chandler wobble period, as these geophysical frequencies are not inherent to the algorithmic structure.

Software Independence: Evidence Against Algorithmic Artifacts

Paper 1's multi-centre cross-validation (R² = 0.920-0.970 over 2.5 years) provides evidence against processing artifacts, as the three analysis centres use fundamentally different software packages with independent algorithmic implementations:

• CODE: Bernese GNSS Software (double-difference processing)
• IGS: GIPSY-OASIS II (undifferenced processing, fiducial-free network)
• ESA: Proprietary processing chain (distinct Kalman filter implementation)

These represent completely different mathematical approaches to clock estimation: different state vector formulations, different stochastic models, different convergence criteria, and independent codebases with no shared source code. The "Kalman filter resonance" hypothesis would require that three independent software packages, developed by different institutions using different algorithms, all accidentally produce nearly identical artifacts (R² > 0.92 over 2.5 years in Paper 1) and that these patterns persist in CODE-only analysis over 25 years (this study). This is statistically implausible—analogous to three independent random number generators producing identical sequences.

Contrast with known artifacts: Well-documented GNSS processing artifacts (multipath, antenna phase center errors, tropospheric modeling) exhibit centre-specific signatures precisely because each software handles them differently. Our findings show the opposite: centre-independent patterns (R² > 0.92), globally coherent structure (mesh score 0.58), and orbital phase-locking (r = -0.888).

Conclusion: While raw carrier-phase reanalysis remains essential for amplitude quantification, the joint constraints from anisotropy structure, orbital-velocity phasing, and hemisphere same-phase behavior are not reproduced by these processing classes. A quantitative bounds table (mechanism × signature) is planned for a follow-on technical note.

4.4 Implications and Open Questions

These findings raise several important questions for fundamental physics:

1. Nature of the Coupling: The orbital velocity correlation (see §3.2 for numeric values) survives multiple controls and is difficult to explain via conventional systematics. If confirmed by multi-centre replication, this would suggest Earth's heliocentric motion influences GNSS timing correlations in ways not accounted for in standard models.
2. Planetary Event Responses: The observed 56 significant responses (25 Bonferroni-corrected) with coherence modulations 0.3–4.0% (median 1.1%, in coherence units) show no correlation with GM/r² gravitational scaling. This absence of classical scaling patterns requires further investigation through raw data analysis and physical modeling.

4.4.2 Physical Plausibility Cross-Check

Question: Are the observed effect magnitudes physically reasonable for gravitational coupling at Earth's surface?

Order-of-Magnitude Comparison:

Key Insights:

1. Correlation length is Earth-scale: λ ~ 4,000 km ≈ 0.6 × R_Earth is consistent with a phenomenon operating at planetary scales, not local atmospheric (λ ~ km) or interplanetary (λ >> 10,000 km) regimes.
2. Anisotropy modulation exceeds naive predictions: The 19% annual variation is ~10² larger than first-order v²/c² ~ 10⁻⁸ predictions, suggesting either nonlinear coupling, collective amplification, or processing-dependent amplitude scaling.
3. Planetary event amplitudes are dimensionally distinct: GPS coherence modulations (dimensionless correlation structure) cannot be directly compared to GR clock rate predictions (Δf/f). The absence of GM/r² scaling indicates either novel physics or processing transfer functions, testable via raw data analysis.
4. Network coordination is global but partial: Mesh score = 0.582 indicates 58% coherent, 42% incoherent behavior—consistent with a global influence partially masked by local noise, not 100% deterministic control.

Plausibility Assessment: All observed magnitudes fall within physically reasonable ranges for Earth-scale gravitational coupling phenomena. The amplitude discrepancies (point 2) and scaling anomalies (point 3) require resolution through raw data analysis but do not invalidate the statistical detections or multi-signature convergence.

The 56 significant planetary event responses show coherence modulations of 0.3–4.0% (median 1.1%, measured in coherence units not Δf/f), indicating subtle but repeatable changes in GNSS network correlation structure during planetary alignments. Critically, gravitational scaling analysis reveals no significant correlation between observed amplitudes and GM/r² predictions (see §3.3), distinguishing this phenomenon from conventional tidal effects.

Absence of Gravitational Scaling

The lack of correlation between observed GPS coherence modulations and GM/r² scaling suggests two possibilities:

1. Novel Phenomenon: The coherence modulation may not be a direct gravitational effect described by classical GR tidal potentials, but rather a distinct phenomenon with different scaling properties
2. Unknown Transfer Function: There may be an intermediate coupling mechanism that modulates the gravitational signal in a non-linear or frequency-dependent manner, obscuring the underlying GM/r² relationship

Processing Center Correction Considerations

GNSS analysis centres (CODE, IGS, ESA) apply systematic error corrections during clock estimation with the explicit goal of producing a self-consistent, globally synchronized clock ensemble. In practice, this means that structured amplitude deviations in individual clocks are prime candidates for removal or attenuation, whereas subtle correlation structure in inter-station phase relationships is less directly targeted by standard correction pipelines. If these centers detect and partially correct for planetary gravitational effects without recognizing their physical origin, several observations become explicable:

• Observed amplitudes: May represent residuals from incomplete corrections rather than full physical signal
• Multi-center correlation (R²=0.920-0.970): Similar correction strategies across centers produce consistent residual patterns
• Temporal variations: Changes in processing algorithms over 25-year period may affect amplitude estimates

This framework naturally reconciles two key findings of this study: (i) the absence of GM/r² scaling in planetary event amplitudes, and (ii) the robust detection of distance-structured correlations and orbital modulation in phase-coherence metrics. If amplitude-like components of the underlying signal are preferentially removed by analysis-center corrections, while residual coherence structure survives in the cross-spectral domain, then one expects weak or distorted amplitude scaling but persistent, reproducible patterns in the correlation field. The methodology adopted here—focusing primary inferences on phase-coherence and distance structure, and treating absolute modulation depths as provisional and processing-sensitive—follows directly from this consideration.

Discriminating Tests Required

Distinguishing between processing artifacts and physical effects requires:

1. Raw carrier phase analysis: Bypass all processing center corrections to measure unprocessed signal amplitudes
2. Processing documentation review: Examine systematic error correction algorithms for planetary gravitational corrections
3. Correction time series analysis: Correlate applied corrections with planetary positions to identify deliberate compensations
4. Multi-constellation testing: Compare GNSS, GLONASS, Galileo, BeiDou to verify processing-independence

Current Assessment

Robust findings (independent of processing):

• Statistical significance: 56/156 events ≥2σ, 25 survive Bonferroni correction (45% of significant responses)
• Detection rates (primary ±120-day window): Mercury 42.5%, Venus 18.8%, Mars 33.3%, Jupiter 34.8%, Saturn 28.0% (Table 3.3.1)
• Absence of GM/r² scaling (see §3.3)
• Multi-center consistency in detection patterns

Uncertain findings: The following items are sensitive to processing choices and will be quantitatively resolved only with raw data validation.

• Absolute amplitude estimates
• Modulation depth magnitudes
• Physical mechanism interpretation

Conclusion: The systematic detection of planetary event correlations across multiple analysis centres provides evidence consistent with gravitational coupling effects in GNSS timing networks. Raw data analysis is essential to determine whether observed amplitudes represent the full physical signal or residuals from processing center corrections.

3. Spatial Structure: The anisotropic correlation field (E-W:N-S = 2.16) persists across 25 years and aligns approximately with Earth's orbital plane. Whether this reflects propagation effects, processing biases, or a genuine geometric coupling remains to be determined through multi-constellation and raw-data analysis.
4. Global Coordination: Network synchronization (index = 0.582) demonstrates coordinated behavior across the 474-station network. This could indicate a global field-like influence, though alternative explanations (correlated processing, shared reference frames, propagation modes) require systematic investigation.
5. Multi-Scale Temporal Structure: The detection of couplings at timescales from days (planetary events) to decades (18.6-year nutation) suggests the phenomenon operates across multiple temporal regimes. The physical mechanism enabling this scale-invariance is unclear.
6. Geophysical Integration: The strong nutation coupling (R² = 0.904 semiannual, R² = 0.641 for 18.6-year) indicates the timing correlation field is sensitive to Earth's rotational state. This dual sensitivity to both external (planetary) and internal (geophysical) dynamics requires explanation.
7. Theoretical Framework: While observations are broadly consistent with TEP predictions (velocity-dependent correlation modulation, geometric anisotropy), no quantitative theoretical model yet predicts the specific effect sizes, modulation depth distribution, or anisotropy ratios. Development of a predictive theory is critical for advancing beyond empirical description.

4.5 Theoretical Implications

Relationship to TEP Framework: The observed phenomenology is qualitatively consistent with TEP predictions of velocity-dependent correlation modulation and distance-structured anisotropy. The absence of GM/r² scaling in planetary event amplitudes is expected: GNSS least-squares processing explicitly removes network-wide frequency offsets (where gravitational mass scaling would manifest), while preserving phase-coherence structure (where the observed correlations reside). This processing characteristic explains why amplitude predictions fail while correlation structure persists—a feature, not a flaw. These results are consistent with TEP's core predictions and establish quantitative targets (e.g., the 2.16:1 anisotropy ratio, the -0.888 orbital correlation coefficient) for future independent validation.

4.5.1 Researcher Degrees of Freedom

The analysis pipeline is rich and has evolved across Papers 1 and 2. To maintain transparency about potential bias:

4.6 Limitations and Caveats

1. Single Analysis Center (Temporal Depth Trade-off): This study uses only CODE data to access the 25-year archive required for long-period analysis (18.6-year nutation, 21 Chandler cycles). While the multi-centre study established cross-centre validation (R² = 0.920-0.970) over 2.5 years, replication of these long-baseline findings with IGS and ESA products (when sufficient historical data becomes available) remains necessary to definitively exclude processing-specific artifacts in the decadal signatures.
2. Orbital Velocity Correlation - Seasonal Confound Discrimination: The orbital correlation (r = -0.888, p < 2×10⁻⁷) survives window size variations, multiple detrending methods, and bootstrap robustness checks. However, orbital velocity correlates with many annual phenomena (solar distance, ionospheric TEC, network participation). Three independent discriminants support genuine orbital coupling over seasonal artifacts: (1) CMB Frame Alignment: The modulation couples to Earth's motion through the CMB rest frame (RA=168°, Dec=-7°) with 5,570× variance ratio over galactic motion—ionospheric effects do not preferentially align with the CMB dipole; (2) Ionospheric Null Controls: No correlation with F10.7 solar flux or Kp geomagnetic indices (Paper 1 validation); (3) Multi-Centre Consistency: Three independent processing centres using different ionospheric models show identical patterns (R² = 0.920-0.970). Note: Hemisphere stratification alone is insufficient to rule out global ionospheric TEC (which peaks in January for both hemispheres); the CMB frame geometry provides the primary geometric discriminant.
3. Planetary Event Analysis - Physical Predictor Modeling: The planetary event analysis documents 56 statistically significant event-associated modulations (25 surviving conservative Bonferroni correction), with detection rates consistent across 25 years. The original study established these are not artifacts of processing or random chance through extensive Monte Carlo permutation testing. However, modeling against physically grounded predictors (e.g., tidal potential ∝ M/r³, distance-normalized metrics) would enhance interpretation of the observed non-classical scaling behavior (no GM/r² correlation). Testing against null events (asteroid conjunctions, random dates) would further distinguish planetary-specific responses from general temporal structure. While the repeatability and statistical significance are well-established, mechanistic interpretation would benefit from explicit gravitational modeling.
4. Unknown Systematics: Unknown systematic effects cannot be definitively ruled out, though any such effect would need to:
   • Correlate with Earth's orbital velocity (see §3.2)
   • Respond to planetary alignments with specific timing patterns
   • Couple to Earth's nutation with R² = 0.90
   • Produce coordinated global network responses
   • Show no correlation with solar activity, ionospheric conditions, or known geophysical variables

   This combination of requirements makes conventional explanations highly unlikely, but not definitively excluded.

5. Theoretical Gap: While TEP provides a qualitative framework, no complete theoretical model yet quantitatively predicts:
   • The absence of GM/r² gravitational scaling
   • The pulse-to-baseline ratio distribution (17–86%, mean 54%)
   • The E-W:N-S anisotropy ratio (2.16)
   • The orbital velocity correlation coefficient magnitude

   Developing such a quantitative theory is critical for advancing from empirical detection to physical understanding.

6. Temporal Coverage: The 18.6-year nutation cycle has only 1.4 complete cycles in our 25.3-year dataset, limiting statistical power for this specific signature. Longer baselines (50+ years) would provide definitive confirmation.
7. Raw Data Access: Analysis uses post-processed clock products. Access to raw GNSS measurements would allow investigation of whether the phenomenon originates in the satellite-ground propagation or ground clock stability.

5. Conclusions: Empirical Evidence for Systematic Gravitational Coupling

This work confirms the multi-centre study over 25.3 years and extends it to long-period regimes. The distance-structured correlation signatures are temporally stable over decadal timescales and not transient artifacts. Building on the multi-centre study's rigorous validation framework (388 statistical tests, comprehensive null testing, cross-centre validation R² = 0.920-0.970), the extended baseline enables investigation of long-period geophysical phenomena inaccessible in shorter datasets, revealing coupling to Earth's 18.6-year nutation cycle (R² = 0.641, p < 10⁻⁸) and documenting ~21.2 Chandler wobble cycles.

The strongest statistical finding is the orbital velocity correlation (r = -0.888, 95% CI: [-0.94, -0.81], p < 2×10⁻⁷, 5.1σ), validated by 0 out of 5,000,000 Monte Carlo surrogates. This correlation—where the spatial anisotropy ratio (EW/NS) tracks Earth's heliocentric velocity with ≈19% annual geometric modulation—survives window size variations, multiple detrending methods, bootstrap resampling, and the multi-centre study's comprehensive controls. Combined with Paper 1's software independence and this study's 25-year temporal stability, this represents strong statistical evidence consistent with systematic temporal-gravitational coupling documented in GNSS data to date.

The planetary event analysis documents 56 statistically significant event-associated modulations (25 surviving ultra-conservative Bonferroni correction, 33 surviving dependency-aware BY-FDR), demonstrating repeatability across 25 years. The 45% Bonferroni survival rate under an intentionally conservative correction (α ≈ 0.00089) indicates robust, repeatable event responses. Critically, null event control testing (Step 2.6) confirms that planetary events produce effect sizes 5.5× larger than 155 random dates (Mann-Whitney p < 10⁻¹⁷, Cohen's d = 1.41, ~8.5σ), validating that detected responses are specific to astronomical alignments rather than temporal autocorrelation artifacts.

The nutation coupling (R² = 0.904 semiannual, R² = 0.641 for 18.6-year), network synchronization (index = 0.582), and persistent anisotropic structure (E-W:N-S = 2.16 over 25 years) provide independent evidence for systematic, globally coordinated patterns in the GNSS timing correlation field.

Study Scope: This temporal extension builds on the multi-centre study's rigorous validation to demonstrate decadal stability and access long-period signatures. (1) Single-centre analysis (CODE) provides the 25-year baseline necessary for 18.6-year nutation detection; the multi-centre study established processing-independence over 2.5 years. (2) Orbital velocity coupling is established through hemisphere stratification (Paper 1 §3.4), which discriminates heliocentric effects from local seasonal confounders. (3) Null event control validates planetary event specificity (5.5× effect size ratio, p < 10⁻¹⁷).

The convergence of temporal stability (25.3 years), multi-signature consistency (orbital, planetary, nutation, network), Paper 1's cross-centre validation (R²=0.92-0.97 across three independent software packages), and Monte Carlo validation (0/5M surrogates) is consistent with systematic gravitational coupling, pending independent confirmation. The orbital velocity correlation (r=-0.888, 5.1σ) persists across 25 solar orbits with statistically robust validation. The 18.6-year nutation coupling (R²=0.641, p<10⁻⁸) provides the first GNSS-based detection at this period. Raw carrier-phase analysis and multi-constellation testing represent logical extensions for distinguishing between competing theoretical frameworks.

Characterization: Inertial Interferometer. The selectivity of detections reveals what this network fundamentally measures: it is blind to gravitational amplitude (no GM/r² scaling), unresponsive to solar surface phenomena (rotation null), and indifferent to static geometry (lunar standstill null)—but exquisitely sensitive to kinematic dynamics (orbital velocity, nutation, CMB-frame motion). This selectivity profile characterizes the network as an inertial interferometer measuring velocity-dependent correlation geometry, not a gravimeter measuring force. The CMB frame alignment identifies the operative reference at cosmic scales: Earth's ~369 km/s motion through the universe's rest frame.

These empirical findings have broader implications:

• For Physics:
  • Documenting systematic correlations between GNSS timing and gravitational/kinematic states provides new constraints for theoretical models of spacetime geometry
  • CMB Rest Frame Coupling: The identification of the Cosmic Microwave Background rest frame as the operative reference (rejecting galactic motion at 5,570× variance ratio) suggests that local temporal correlations couple to Earth's absolute motion through the Universe's rest frame. This identifies the GNSS network as a potential detector for absolute kinematic effects.
  • Reference Frame Question: If confirmed by independent replication, CMB frame coupling would raise questions about reference frame dependence in distributed clock synchronization, potentially consistent with Machian frameworks where local properties are influenced by the global matter distribution
  • Terrestrial Detection: This reports detection of CMB frame coupling in ground-based atomic clock networks, complementing space-based cosmological observations
• For Metrology: Demonstrating that distributed atomic clock networks exhibit coordinated variations sensitive to orbital and geophysical dynamics suggests new approaches to precision timing
• For Technology: Motivating development of correlation-based detection methods complementing traditional single-clock stability metrics

5.1 Future Work

Critical Priority (Methodological Validation):

1. Raw Carrier Phase Analysis (Testing Ξ Asymmetry): Collaborate with IGS to access unprocessed measurements to test a different TEP prediction:
   • Objective: Detect one-way link asymmetry Ξ_AB (§10.2) rather than network correlation C(r) (§10.3)
   • Challenge: Separating TEP signal from ~10⁸× larger classical effects (GR, atmosphere)
   • Outcome: If null, confirms TEP physics is strictly global/topological; if positive, confirms local non-reciprocity
2. Processing Center Documentation Review: Examine systematic error correction algorithms:
   • Review CODE, IGS, ESA processing documentation for planetary/orbital corrections
   • Analyze correction time series for correlation with planetary positions
   • Identify changes in processing algorithms over 25-year period
   • Would clarify: Whether centers are unconsciously correcting for detected signals

Highest Priority (Further Strengthening):

3. Multi-Center Long-Baseline Replication: When IGS/ESA extend historical archives:
   • Replicate 25-year analysis with IGS and ESA products
   • Compare CODE, IGS, ESA results for long-period signatures (18.6-yr nutation, Chandler)
   • Would strengthen: Confirms heliocentric velocity coupling independent of radial distance variations
4. Planetary Event Physical Modeling: Model events with gravitational predictors:
   • Gravitational scaling tests: Neither M/r² nor M/r³ shows significant correlation with event amplitude (all p > 0.5), consistent with GNSS processing suppression of classical gravitational signatures
   • Distance-normalize all events to standard range
   • Null event validation: Random dates show 5.5× smaller effect sizes than planetary alignments (p < 10⁻¹⁷), confirming event specificity
   • Model phase-dependence (approach vs recession)
   • Would strengthen: Mechanistic interpretation of observed non-classical scaling pattern (absence of GM/r² dependence)

High Priority (Methodological Validation):

1. Anisotropy Network Geometry Controls:
   • Grid-Balanced Subsampling: Repeat analysis on a spatially decimated network (e.g., one station per 500x500km grid cell) to remove density biases (US/Europe dominance).
   • Stratify by latitude band to test equatorial vs polar differences
   • Normalize by station pair density per azimuth
   • Control for time-varying network configuration
   • Test whether anisotropy varies with local solar time
2. Network Coherence Mechanistic Tests:
   • Stratify by common satellite visibility (test if coherence reflects shared SV clock errors)
   • Compare stations in same vs different processing batches
   • Synthetic data with identical processing chain but randomized positions
3. Raw Data Analysis: Collaborate with IGS to access raw carrier phase measurements, enabling:
   • Separation of propagation effects from clock stability
   • Investigation of ionospheric vs geometric contributions
   • Direct testing without processing-induced correlations

Medium Priority (Extensions & Replications):

1. Multi-Constellation Testing: Extend to GLONASS, Galileo, BeiDou, QZSS
2. Optical Clock Networks: Test with higher-precision systems (ACES, T-TEL)
3. Historical Extension: Analyze 1994-2000 IGS data for 1.7 complete 18.6-year cycles
4. Theoretical Development: Quantitative models predicting specific effect sizes, scaling relations, anisotropy ratios

5.2 Final Remarks

The global GNSS infrastructure, designed for navigation, provides a tool for investigating systematic patterns in distributed timing networks. Building on the multi-centre study's comprehensive validation, this 25-year analysis shows that the detected signatures are temporally stable over decadal timescales and extend to long-period geophysical phenomena. The orbital velocity correlation (r = -0.888, p < 2×10⁻⁷) and the detection of 18.6-year nutation coupling (R² = 0.641) provide empirical evidence consistent with a phenomenology involving velocity-dependent spacetime-geometry modulation, but they do not by themselves establish a unique underlying theory.

The convergence of independent evidence streams—spatial anisotropy, orbital coupling, geophysical integration, and network coordination—creates a robust phenomenological framework that substantially constrains simple artifact explanations. Paper 1's software independence (three different packages: Bernese, GIPSY-OASIS, ESA proprietary achieving R²>0.92 over 2.5 years) combined with this study's 25-year temporal stability and successful falsification tests (solar rotation null) strengthen the case for physical coupling. The mechanistic explanation is precise: CODE's Undifferenced Ionosphere-Free Linear Combination (L3) with datum constraint (zero-sum or fixed reference clock) filters network-wide potential shifts (where classical GM/r² manifests) while transmitting differential phase structure (where correlation geometry resides)—a selectivity independently validated across all three processing centres. These findings are consistent with theories predicting velocity-dependent spacetime-geometry modulation (e.g., TEP), though confirmation requires independent replication by other groups using different processing pipelines and constellations. Raw carrier-phase analysis and multi-constellation testing represent natural extensions that will further constrain theoretical interpretations and may distinguish between competing physical frameworks.

6. Analysis Package

Complete Reproducible Pipeline for 25-Year Long-Baseline Analysis

This work provides a complete, reproducible 5-step pipeline for testing TEP predictions using 25.3 years of GNSS data (2000-2025), extending the multi-center study to investigate decadal stability and long-period phenomena:

Supplementary Materials

Supplementary Table S1: Multi-Resolution Convergence Analysis

Verification of CMB frame alignment across three independent grid resolutions to demonstrate robustness and rule out resolution-dependent artifacts.

Interpretation: Results converge monotonically across a 93-fold range in grid density with diminishing returns (+5.14% → +0.05% → +0.27%), indicating asymptotic approach to a correlation limit. The marginal 0.27% improvement at 1° resolution despite 6× more grid points (10,585 → 65,341) confirms convergence. All four resolutions identify the CMB dipole (not Solar Apex) as the preferred reference frame, with best-fit locations showing stability within 1-2° (RA: 180°→190°→187°→186°; Dec: 0°→-5°→-5°→-4°). The 5° resolution captures 99.6% of the achievable correlation while remaining computationally efficient and consistent with standard practice in cosmological frame studies (Planck ~10°; WMAP ~5-10°). The 1° ultra-high resolution analysis (65,341 directions) provides additional verification and is used for publication figures. The systematic monotonic improvement with finer resolution, combined with asymptotic convergence and location stability, indicates the signal is robust rather than a numerical artifact. CMB separation distance converges to 18.2° (closest alignment), while Solar Apex separation remains consistently high (86-89°) across all resolutions, inconsistent with the galactic motion hypothesis.

Reproducibility Statement: This pipeline provides complete reproducibility for the 25-year long-baseline analysis presented in this manuscript. All code, data processing steps, validation procedures, and analysis parameters are fully documented and version-controlled. The analysis builds upon the validated methodology from the multi-center study while extending temporal coverage to enable investigation of decadal stability and long-period geophysical phenomena.

References

Data Sources

@misc{Smawfield_TEP_GNSS_Longspan_2025,
  author       = {Matthew Lukin Smawfield},
  title        = {Global Time Echoes: 25-Year Temporal Evolution of 
                  Distance-Structured Correlations in GNSS Clocks (Cairo v0.16)},
  year         = {2025},
  publisher    = {Zenodo},
  doi          = {10.5281/zenodo.17517141},
  url          = {https://doi.org/10.5281/zenodo.17517141},
  note         = {Preprint}
}