Global Time Echoes: 25-Year Analysis of CODE Precise Clock Products

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, this work proposes 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 an expected consequence under the datum-projection transfer hypothesis, 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, temporal correlations and global light-time relations, while the locally measured speed of light remains invariant.

This investigation is grounded in empirical phenomenology: it systematically documents patterns in GNSS timing correlations over multi-decade timescales and interprets them 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 documented empirically in this work.

Formally, the temporal field is covariant and introduces no fundamental preferred frame; a non-zero environmental or cosmological field configuration may nevertheless define an observational covariance direction.

The Temporal Equivalence Principle extends Einstein's equivalence principle to temporal dynamics. Key predictions—velocity-dependent anisotropy and event responses—provide testable signatures that are investigated 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 relevant cosmological/background covariance 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: a conformal-sector covariance with characteristic decay length λ_T, whose velocity/orientation dependence is a channel response whose microscopic transfer law is not yet closed (§1.2)

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: The TEP interpretation of the GNSS covariance is conformal-dominant: the observed data confirm a finite, distance-structured covariance scale within the pre-specified TEP range (λ_T), and the absence of GM/r² scaling is consistent with the conformal-sector Temporal Topology. The exponential kernel strongly outperforms simple power-law models and is retained for cross-paper comparability, although Gaussian-family kernels provide the preferred empirical fit in the present dataset (ΔAIC = 12.8). The velocity-dependent anisotropy (r = −0.888) and geometric alignment (EW/NS = 2.16) are qualitatively consistent with a conformal shear projected through the GNSS processing transfer function, but their magnitudes are not yet derived from first principles—the microscopic transfer law mapping Σ_μ^obs to the observed anisotropy ratio and orbital-correlation coefficient remains open (§5). The observed 19% modulation is a network-level covariance-shape observable and cannot presently be compared directly with a microscopic conformal amplitude. The datum-projection mechanism provides a qualitative framework: the data are consistent with a processing architecture that suppresses common-mode amplitude information while retaining differential covariance structure. The magnitude and sign of this transfer must be calibrated using raw/RINEX data and synthetic end-to-end processing (§4.3.1.1). The disformal sector B(φ) is not tested by this open-path covariance channel and is reserved for dedicated closed-loop holonomy and one-way asymmetry experiments (Paper 8).

1.3 Previous Work

The multi-centre study established the phenomenon across independent analysis centres (CODE, IGS, ESA), reducing the likelihood of centre-specific processing artifacts over 2.5 years. That study demonstrated velocity-dependent anisotropy and anomalous planetary event responses, using rigorous cross-centre 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, core 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 Paper 1 (multi-centre study)'s comprehensive validation framework, which established cross-centre 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 long-baseline consistency 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-Centre 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 long-baseline consistency. The multi-centre study's comprehensive validation framework—including 388 statistical tests, extensive null testing, ionospheric controls, and cross-centre consistency—provides the methodological foundation upon which these temporal extensions are built. The convergence of independent findings across multiple analysis centres 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 consistent with the expected behavior of the datum-projection mechanism. The processed products isolate the differential residual channel that the TEP interpretation predicts should survive; its quantitative transfer remains to be validated by synthetic injection.

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 centres (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. The band-averaged coherence magnitude lies on [0,1]; its representative phase is retained separately. 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 absolute coherence amplitudes) 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 centres 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
• λ_T = Temporal Topology correlation 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. The effective degrees of freedom are approximated 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. With r = −0.888 and N_eff ≈ 11, the autocorrelation-corrected Pearson p-value is p ≈ 2.6×10⁻⁴. Separately, a Monte-Carlo surrogate test with 5,000,000 phase-randomised surrogates yields 0 exceedances, giving an empirical upper-bound p < 2×10⁻⁷ under the specified phase-randomisation null. These two p-values test different nulls (parametric t-distribution vs. phase-randomised surrogate) and are reported separately; the surrogate bound is the stronger of the two.

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), ±120 days is pre-declared 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 datum-projection 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. Two competing physical hypotheses are tested:

• Galactic Hypothesis: Modulation driven by Solar System's peculiar motion relative to the local standard of rest (~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, four grid resolutions were tested 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 dipole direction 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, two competing geometric models were explicitly evaluated 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-Centre Validation

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

3. Results: Multi-Signature Convergence Across a 25.3-Year Baseline

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 long-baseline consistency 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-pair observations 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, EW > NS correlation lengths are again observed. 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, Temporal Topology 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, the exponential decay kernel (C(d) = A·exp(-d/λ) + B) is retained for directional analyses. Alternative kernels are also evaluated 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, seven candidate spatial-correlation kernels were evaluated 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—is retained for cross-paper comparability rather than because it is selected by the present data (Δ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, it is retained 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 a 25.3-year baseline, testing whether this coupling is consistently recovered 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 a primary detection is found: 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])
• Parametric p-value (autocorr-corrected t-test, N_eff ≈ 11): p ≈ 2.6×10⁻⁴
• Surrogate test: 0/5,000,000 phase-randomised surrogates exceeded observed value (empirical p < 2×10⁻⁷ under the specified null)
• 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; t-test p ≈ 2.6×10⁻⁴ with N_eff ≈ 11; surrogate p < 2×10⁻⁷, 0/5M 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 is observed: 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, a baseline model was fit 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, surrogate p < 2 × 10⁻⁷ (0/5M)
• Seasonal model: an unconstrained annual harmonic can reproduce generic annual structure. The orbital-speed template provides the physically constrained fit, explaining 74.6% of the variance with phase and waveform fixed by celestial mechanics.

Conclusion: The orbital-speed template provides the physically constrained fit, explaining 74.6% of the variance with the phase and waveform fixed by celestial mechanics. Combined with directional anisotropy, same-phase hemispheric response and CMB-region alignment, this supports the TEP interpretation as velocity-dependent Temporal Shear rather than an unconstrained seasonal coincidence. A formal residual comparison reporting R², adjusted R², AIC/BIC, and held-out prediction error for both models on the same observations is a natural next step.

3.2.4 Control Analyses: Hemisphere & Space Weather Stratification

To rigorously test for environmental confounding, the 25-year dataset was stratified 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. Two competing physical hypotheses are tested:

• Galactic Hypothesis: The Solar System's peculiar motion relative to the local standard of rest (~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: A comprehensive sky search was performed 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). Note: this test ranks the observed best-fit direction against random directions in the observed dataset; it is not a full-search surrogate null (which would require repeating the complete 65,341-direction maximization for each surrogate dataset). The directional map and variance ratio are reported as exploratory evidence; a global search-corrected p-value is a natural next step.
• Bootstrap 95% CI: [0.568, 0.878] from 1,000 resamples, confirming robustness.
• Variance Ratio: The best-fit directional template explains ~10,300× more variance than Solar Apex (R²_best-fit/R²_Apex = 55.7% / 0.0054% ≈ 10,300). The exact CMB-dipole template explains ~6,600× more variance (35.7% / 0.0054%). Both strongly discriminate against the Solar Apex hypothesis.
• 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-region direction (18-22° separation) over the Solar-Apex direction (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 supports the CMB-directional template as the preferred covariance direction for temporal modulation.

Alternative Hypothesis Testing: The following alternatives were systematically ruled out:

• Random Coincidence: Permutation test shows p=0.0002 against random sky directions (not a full-search surrogate null).
• 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 the CMB dipole direction (r=0.744 at RA=190°, Dec=-5°; 22° from CMB dipole; directional-rank p=0.0002), not the Solar Apex direction (r=0.007; 89° from best fit). The best-fit directional template explains R²=55.7% of variance (1° resolution, RA=186°, Dec=-4°, 18.2° from CMB), while the exact CMB-dipole template explains R²=35.7%. Both strongly discriminate against the Solar Apex (R²≈0.005%), with variance ratios of ~10,300× (best-fit) and ~6,600× (exact CMB) respectively. The 5° grid resolution (2,701 directions tested) provides finer localization. This finding represents a seventh independent signature consistent with velocity-dependent temporal coupling. The grid scan fixes background speed and uses a direction-only predictor, so it identifies a preferred direction but does not establish that the coupling amplitude is controlled by the 369 km/s speed magnitude (§2.3.5). The quoted p-value is a directional-rank statistic within the observed dataset; a global sky-search surrogate significance has not yet been computed.

3.3 Planetary Event Responses

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

Note: A statistically robust phenomenological anomaly is reported 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; that set is reproduced here in the present pipeline and then scaled 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 correction across the 156 primary-window event tests).

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σ). Family-wide multiple-testing corrections were applied across all 156 tested events: 19 survive Bonferroni correction (α = 0.05/156 = 0.000321), 19 survive Holm step-down correction, 38 survive BH-FDR, and 28 survive dependence-aware BY-FDR. Under the complete null, only 0.05 Bonferroni survivors are expected. This high survival rate demonstrates robust, repeatable event responses:

3.3.2 Modulation Depth Distribution

Across all 156 tested events, the normalized modulation depths show a wide but meaningful distribution:

3.3.3 Gravitational Scaling Analysis

Objective: Test whether observed amplitudes follow planetary acceleration (M/r²) or tidal-gradient (M/r³) scaling proxies.

Key Finding: The observed GPS coherence modulations do not correlate with either gravitational (M/r²) or tidal (M/r³) 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 Centre Corrections: This analysis uses processed GNSS clock products from analysis centres (CODE, IGS, ESA) that apply systematic error corrections during data processing.

Potential Impact: If analysis centres 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.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.

Among the 56 significant events, normalized modulation depth ranges from 36.7–85.1%, with a typical value of approximately 55%. The corresponding absolute pulse amplitude is typically approximately 1.1% of the full coherence scale.

Interpretation: The typical absolute amplitude 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.4 Geophysical Couplings

Objective: Test TEP-motivated expectations for coupling to Earth's rotational dynamics, specifically the 18.6-year lunar nutation (expected R² ~ 0.6-0.7) and semiannual nutation (expected R² > 0.85) that were inaccessible in Paper 1's 2.5-year baseline. This represents the first test of genuinely novel expectations 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 the semiannual nutation cycle (R² = 0.904, p = 2.7×10⁻⁵), and a preliminary 18.6-year nutation-aligned component observed (R² = 0.641, 1.4 cycles; p = 0.010, pending red-noise surrogate validation), supporting Claim E. The semiannual component explains 90.4% of variance in the 12 phase-bin means. Notably, the annual period shows no significant coupling (R² = 0.0178, p = 0.92), 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: At the prescribed 433-day Chandler period, the network exhibits a weak but phase-coherent response across approximately 21 cycles (R² = 0.096, phase stability = 0.72). Its amplitude falls below the primary significance threshold, so it is treated as a Chandler-consistent secondary signature rather than an independent detection. The extended phase coherence at the known geophysical period distinguishes the response from random noise, which would not maintain phase coherence at the correct 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 consistently recovered or degrades over multi-decade 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 long-baseline recovery 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 a weak but phase-coherent response at the prescribed 433-day period (R² = 0.096, coherence = 0.72 across 21.2 complete cycles, p = 0.002), treated as a Chandler-consistent secondary signature rather than an independent detection. 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 convergent 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 is recovered consistently across the complete 25.3-year aggregate.

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) 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 (surrogate p < 2×10⁻⁷; t-test p ≈ 2.6×10⁻⁴ with N_eff ≈ 11). 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: alignment with the Solar Apex (~20 km/s, galactic motion) implies one physical mechanism, while alignment with the CMB dipole (~369 km/s, cosmological background 4-velocity) implies a different physical regime: CMB-frame-aligned covariance—a correlation between the TEP field configuration and the observable cosmological background state, without introducing a fundamental preferred frame.

Competing Hypotheses:

Directional 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 a preferred directional template under the fixed-amplitude scan:

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, surrogate p<2×10⁻⁷). The CMB-dipole direction is strongly consistent with the observed anisotropy geometry.

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=0.006) 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 rules out a purely isotropic scalar-only model. It remains compatible both with a kinematic vector projection and with a directional scalar-field gradient such as Σ_μ = ∇_μ ln A, which is the central object in the conformal TEP sector. Discriminating between these possibilities requires explicit pairwise projection kernels. The hemisphere asymmetry, previously puzzling, is consistent with CMB-frame geometry under either interpretation.

Broader Significance:

This finding elevates the result from "orbital modulation is detected" to "the detection of CMB-frame-aligned covariance." This has profound implications:

1. Directional Alignment: The best-fit direction aligns with the CMB dipole (cosmological background 4-velocity), not the Solar Apex (galactic motion). The grid scan identifies a preferred direction; because background speed is fixed at 20 km/s and the predictor is direction-only (cos(Dec_net)), the scan does not fit the speed magnitude and therefore does not, by itself, establish that the coupling amplitude is controlled by 369 km/s rather than 20 km/s. The result is a directional identification with the CMB dipole, not a speed-magnitude measurement. This suggests CMB-frame-aligned covariance without introducing a fundamental preferred frame.
2. Reference Frame Question: If confirmed, the CMB dipole-direction 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-region directional result (r=0.747, directional-rank p<0.001, ~10,300× variance ratio over the Solar-Apex template) provides strong directional discrimination; global sky-search significance remains to be evaluated. This transforms the phenomenology from "velocity-dependent" to "CMB-directional," 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 dipole direction, 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-Directional Alignment (§4.1.2.1): The modulation aligns with the CMB dipole direction (RA=168°, Dec=-7°) with ~10,300× variance ratio over the Solar-Apex directional template. 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 dipole-direction 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: The analysis observes a massive correlation of r = -0.888; surrogate p<2×10⁻⁷ (0/5 million exceeded), with autocorrelation-adjusted t-test p≈2.6×10⁻⁴. This is physically inconsistent with a drifting draconitic signal. The high correlation demonstrates that 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

Event-locked windows around planetary alignments are analyzed 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 (family-wide correction across all 156 tested events)
• Expected false positives: 156 × 0.000321 = 0.05 events (not 8 under uncorrected α=0.05)
• Observed: 19 events survive family-wide Bonferroni correction (380× more than expected by chance)
• FDR control: 38 of 156 events survive Benjamini-Hochberg FDR at q = 0.05 (assuming independence); 28 of 156 survive BY-FDR (accounting for dependency)
• Bonferroni survival rate: 19/156 = 12% of all tested events survive family-wide correction (vs 0.05 expected under null)

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; 240-day window spans 2.73 orbital periods (273%)
• Venus: 225-day period; 240-day window spans 1.07 orbital periods (107%)
• Mars: 687-day period; 240-day window spans 0.35 orbital periods (35%)
• Jupiter: 4,333-day period; 240-day window spans 0.055 orbital periods (5.5%)
• Saturn: 10,759-day period; 240-day window spans 0.022 orbital periods (2.2%)

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.

Absolute Coherence Amplitudes: 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 consistent with the expected datum projection of common-mode components in IGS clock estimation. 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. This is resolved by explicitly defining the signal processing transfer function:

This hypothesised dual mechanism would create a natural dichotomy: energetic effects (proportional to gravitational potential) would be removed, while geometric effects (correlation structure, anisotropy patterns) would be preserved. The datum-projection mechanism is consistent with cross-centre consistency across three analysis centres using fundamentally different software implementations, though a definitive test requires synthetic signal injection through the actual processing chain:

The consistency across independent processing chains (R² = 0.920-0.970, Paper 1) supports that datum projection and common-mode suppression are reproducible features of GNSS processing methodology, not a post-hoc rationalization. If the "datum-projection suppresses mass signals" were merely an excuse, different software implementations would show different results. Cross-centre reproducibility shows that this suppression is not specific to one software implementation; it does not by itself establish the physical transfer mechanism.

Datum projection and common-mode suppression are reproducible features of GNSS processing methodology; the specific transmission of a TEP covariance signal remains a falsifiable transfer hypothesis to be tested by synthetic signal injection.

Predicted Behavior in Processed vs Raw Data:

Mechanism of Transmission (Why Signal Survives):

The survival of these correlation structures has a concrete candidate mechanistic explanation in 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 detected in this analysis occupy this "transmitted" differential channel.

Supporting Evidence from GNSS Literature: The GNSS community has independently documented the exact phenomena detected in this analysis—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 measured in this analysis (λ_T = 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 preliminary nutation-aligned component is observed over 1.4 cycles, providing candidate evidence for sensitivity to multi-decadal astronomical structure pending red-noise validation.

Chandler Wobble (14-month polar motion):

• R² = 0.096. At the prescribed 433-day Chandler period, the network exhibits a weak but phase-coherent response across approximately 21 cycles (phase stability = 0.72). Its amplitude falls below the primary significance threshold (R² > 0.15), so it is treated as a Chandler-consistent secondary signature rather than an independent detection. 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

The convergence of independent evidence streams on systematic patterns is summarized:

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-Directional Alignment identifies a preferred directional template near the CMB dipole (not the Solar Apex direction), resolving the hemisphere asymmetry. The present fixed-magnitude directional scan does not measure the relevant background-speed amplitude (§2.3.5)
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-directional alignment (§4.1.2.1) is particularly significant: it transforms the phenomenology from "velocity-dependent modulation" to "CMB-frame-aligned covariance candidate," providing strong directional discrimination between the CMB-dipole and Solar-Apex templates (Solar Apex rejected at ~10,300× 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 test the datum-projection mechanism 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 temporal field (Temporal Shear), not to absolute potential. Temporal Shear suppression in dense environments therefore decouples the observed geometric amplitudes from simple GM ⁄ r² expectations. Future raw-carrier analyses will be required to test mass-scaling directly.

A critical validation of long-baseline recovery 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; this study confirms and extends to 56/156 events ≥2σ, of which 19 remain family-wide 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.3-year baseline vs 2.5 strengthens velocity correlation inference

New (Long-Period Access):

• 18.6-Year Nutation: Preliminary — R² = 0.641, p = 0.010 (1.4 cycles; pending red-noise validation; multi-centre study: inaccessible)
• Long-baseline Stability: Confirms signatures across 2+ solar cycles

Replicated (Confirms Temporal Stability):

• Network Covariance Score (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 reproducible across the processed products examined so far over 2.5 years. This study extends temporal coverage but cannot confirm whether long-baseline signatures (18.6-year nutation) are reproducible across other processing centres until IGS and ESA archives extend backward. Cross-centre agreement establishes robustness across processed-product pipelines, but does not alone determine the physical origin because the centres share observations, conventions and some modeling assumptions.

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 convergent signatures:

4.3.1.1 Order-of-Magnitude Considerations

The 19% geometric modulation observed in EW/NS ratios is a network-level anisotropy ratio rather than a fractional clock-frequency amplitude. A naive v²/c² estimate yields ~10⁻⁸ for orbital speeds, but this quantity is not directly comparable to the observed anisotropy ratio because the two are different observables: one is a microscopic energetic amplitude, the other is a macroscopic covariance-shape parameter projected through the GNSS measurement and processing chain. The magnitude and sign of the GNSS transfer function remain to be calibrated with raw/RINEX data and end-to-end signal injection. For the detailed architecture, see §2.1.3.2.

Context: A naive kinematic benchmark at orbital velocity is v²/c² ≈ 10⁻⁸. The observed 19% anisotropy ratio is a geometric covariance-shape parameter, not a direct measurement of this microscopic amplitude. Because the 19% observable is a network-level anisotropy ratio rather than a fractional clock-frequency amplitude, it cannot presently be compared directly with a microscopic conformal coupling estimate. The measurement-transfer law mapping Σ_μ^obs to the observed anisotropy ratio remains open (§5).

Resolution: Information Resonance vs. Force Coupling

This distinction provides a possible qualitative interpretation of the apparent conflict between the "stiff" response (19% modulation) and the "weak" coupling (no mass scaling):

• Amplitude = Energy (Suppressed): Simple energetic couplings ($GM/r^2$) are consistent with absorption by datum estimation (see §2.1.3); their absence in differential covariance is expected under a priori relativistic modelling.
• Phase = Information (Transmitted): The network detects the information content of the gravitational configuration—specifically, the synchronization geometry.

This could account for why the observed response does not follow a simple planetary-mass hierarchy. The network may not be weighing the planets (force coupling); it may be 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-directional alignment, planetary events) are reproducible across the processed products examined so far; their raw-data origin remains to be established, and 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.087 (p = 0.522)
• Amplitude vs Gravity (M/r²): r = -0.042 (p = 0.759)

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-Transfer Hypothesis (§2.1.3): The absence of mass scaling is consistent with the hypothesised high-pass transfer nature of GNSS processing described in §2.1.3. If the energetic component (proportional to $GM/r^2$) is suppressed by least-squares adjustment, the surviving signal would be phase-coherence geometry, which depends on alignment rather than mass. This interpretation remains a hypothesis pending synthetic signal injection through the actual processing chain.

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, explicit null models were evaluated and their qualitative predictions compared 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 ~10,300× 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)

The analysis qualitatively assesses 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 Diversity and Cross-Centre Reproducibility: 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 software-diverse products from three processing centres, 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 (19 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% anisotropy ratio is a macroscopic covariance-shape parameter, not directly comparable to the microscopic v²/c² ~ 10⁻⁸ energetic amplitude. The transfer function between these observables remains to be derived.
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 datum-projection transfer, 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 Centre 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-centre 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-centre 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 centre 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 cross-centre reproducibility

Current Assessment

Robust findings (independent of processing):

• Statistical significance: 56/156 events ≥2σ, 19 survive family-wide Bonferroni correction (α = 0.05/156)
• Detection rates (primary ±120-day window): Mercury 42.5%, Venus 18.8%, Mars 33.3%, Jupiter 34.8%, Saturn 28.0% (§3.3.1)
• Absence of GM/r² scaling (see §3.3)
• Multi-centre 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 centre 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 covariance score (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 consistent with the datum-projection transfer hypothesis: GNSS least-squares processing removes network-wide frequency offsets (where gravitational mass scaling would manifest), while preserving phase-coherence structure (where the observed correlations reside). This provides a candidate explanation for why amplitude scaling is absent 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 Centre (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 long-baseline 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-Directional Alignment: The modulation aligns with the CMB dipole direction (RA=168°, Dec=-7°) with ~10,300× variance ratio over the Solar Apex directional template—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 dipole-direction geometry provides the primary geometric discriminant.
3. Planetary Event Analysis - Physical Predictor Modeling: The planetary event analysis documents 56 statistically significant event-associated modulations (19 surviving family-wide 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 recovered consistently across the 25.3-year aggregate and are 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 a strong preliminary 18.6-year nutation-aligned component (R² = 0.641, observed over 1.4 cycles; p = 0.010, pending red-noise surrogate validation) 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]; surrogate p < 2×10⁻⁷, 0/5M exceeded; t-test p ≈ 2.6×10⁻⁴ with N_eff ≈ 11), 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 diversity and cross-centre reproducibility and this study's 25.3-year aggregate consistency, 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 (19 surviving family-wide Bonferroni correction at α = 0.05/156, 19 surviving Holm step-down correction, 38 surviving BH-FDR, 28 surviving dependency-aware BY-FDR), demonstrating repeatability across 25 years. The 19 Bonferroni survivors—versus 0.05 expected under the complete null—indicate 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; preliminary R² = 0.641 for 18.6-year, pending red-noise validation), network covariance score (index = 0.582), and persistent anisotropic structure (EW/NS = 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 long-baseline recovery 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 cross-centre reproducibility 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 long-baseline recovery (25.3 years), multi-signature consistency (orbital, planetary, nutation, network), Paper 1's cross-centre validation (R²=0.92-0.97 across software-diverse products from three processing centres), and Monte Carlo validation (0/5M surrogates) is consistent with systematic gravitational coupling, pending independent confirmation. The orbital velocity correlation (r=-0.888, surrogate p < 2×10⁻⁷) persists across a 25.3-year baseline with statistically robust validation. The 18.6-year nutation coupling (R²=0.641, observed over 1.4 cycles; p=0.010 pending red-noise surrogate validation) provides preliminary evidence for 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-directional alignment). This selectivity profile characterizes the network as an inertial interferometer measuring velocity-dependent correlation geometry, not a gravimeter measuring force. The CMB-directional alignment identifies the preferred covariance direction at cosmic scales: alignment with the cosmological background 4-velocity (CMB dipole direction), with the speed-scale magnitude to be tested in future work.

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-Frame-Aligned Covariance: The identification of the CMB dipole direction as the preferred reference direction (rejecting the Solar Apex at ~10,300× variance ratio) suggests that local temporal correlations align with the cosmological background 4-velocity. This identifies the GNSS network as revealing a secondary covariance-frame signature. The grid scan identifies a preferred direction; it does not by itself establish that the coupling amplitude is controlled by the 369 km/s speed magnitude (§2.3.5).
  • 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 Evidence for CMB-Directional Covariance: This reports directional evidence for CMB-aligned covariance 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 (Paper 0, §10.2) rather than network correlation C(r) (Paper 0, §10.3)
   • Challenge: Separating TEP signal from ~10⁸× larger classical effects (GR, atmosphere)
   • Outcome: If null, favours a conformal-dominant Earth-scale branch and places a stronger bound on the non-exact/disformal (B) sector; if positive, provides evidence for one-way non-reciprocal transport
2. Processing Centre 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-Centre 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 recovered consistently across the 25.3-year aggregate and extend to long-period geophysical phenomena. The orbital velocity correlation (r = -0.888, p < 2×10⁻⁷) and the preliminary 18.6-year nutation-aligned component (R² = 0.641, pending red-noise validation) 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 diversity and cross-centre reproducibility (three different packages: Bernese, GIPSY-OASIS, ESA proprietary achieving R²>0.92 over 2.5 years) combined with this study's 25.3-year aggregate consistency and successful falsification tests (solar rotation null) strengthen the case for physical coupling. The datum-projection mechanism provides a concrete candidate mechanism: datum constraints and common-mode estimation can absorb network-wide offsets while retaining differential phase structure. Cross-centre reproducibility is consistent with this mechanism, but synthetic signal injection through the actual processing pipelines is required to calibrate and validate the transfer quantitatively. 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-centre study to investigate long-baseline consistency and long-period phenomena:

Supplementary Materials

Supplementary Table S1: Multi-Resolution Convergence Analysis

Verification of CMB-region directional alignment across four grid resolutions to demonstrate robustness and rule out resolution-dependent artifacts. The p-values are directional-rank statistics within the observed dataset; global sky-search surrogate significance has not yet been computed.

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-centre study while extending temporal coverage to enable investigation of long-baseline consistency and long-period geophysical phenomena.

References

Data Sources

@misc{Smawfield_TEP_GNSS_Longspan_2025,
author = {Matthew Lukin Smawfield},
title = {Global Time Echoes: 25-Year Analysis of CODE Precise Clock Products (v0.19 Cairo)},
year = {2025},
publisher = {Zenodo},
doi = {10.5281/zenodo.17517141},
url = {https://doi.org/10.5281/zenodo.17517141},
note = {Preprint}
}