Global Time Echoes: Distance-Structured Correlations in GNSS Clocks

Key Theoretical Predictions

1. Spatial correlation structure: Clock frequency residuals should exhibit exponential distance-decay correlations $C(r) = A \cdot \exp(-r/\lambda) + C_0$
2. Correlation length range: For screened scalar fields in modified gravity, $\lambda$ typically ranges from $\sim$1,000 km (strong screening, $m_\phi \sim 10^{-4}$ km$^{-1}$) to $\sim$10,000 km (weak screening, $m_\phi \sim 10^{-5}$ km$^{-1}$), corresponding to Compton wavelengths $\lambda_C = \hbar/(m_\phi c)$ of potential screening mechanisms
3. Universal coupling: The correlation structure should be independent of clock type and frequency band (within validity regime)
4. Multi-center consistency: Independent analysis centers should observe the same correlation length $\lambda$
5. Falsification criteria: $\lambda < 500$ km or $\lambda > 20,000$ km would rule out screened field models; a coefficient of variation across centers >20\% would indicate systematic artifacts

Four-Step Analysis Pipeline

Step 1: Data Acquisition

1. 1.0 Provenance Documentation: Establishes computational provenance with version tracking and execution logging
2. 1.1 TEP Data Acquisition: Acquires GNSS clock data from three independent analysis centers (CODE, IGS, ESA) covering 1 January 2023 to 30 June 2025
3. 1.2 Coordinate Validation: Validates station coordinates against ITRF2014 with ECEF validation and spatial analysis

Establishes computational provenance and acquires GNSS clock data from three independent analysis centers with rigorous quality controls.

Step 2: Core Analysis

1. 2.0 TEP Correlation Analysis: Implements phase-coherent cross-spectral density analysis in the 10-500 µHz band, extracting phase-alignment index correlations and fitting exponential decay models
2. 2.1 Aggregate Geospatial Data: Consolidates station pair correlations with geographic metadata and distance binning
3. 2.2 Geospatial Temporal Analysis: Analyzes correlations with astronomical events, orbital tracking, anisotropy patterns, and temporal field dynamics

Detects exponential correlation decay patterns using phase-coherent cross-spectral density analysis.

Step 3: Validation Suite

1. 3.0 Cross-Validation Suite: Block-wise validation (monthly/spatial), Leave-One-Station-Out (LOSO), Leave-One-Day-Out (LODO), and block bootstrap analyses
2. 3.1 Robust Block Bootstrap: Bootstrap confidence intervals with 1000 iterations
3. 3.2 TEP Null Tests: Null hypothesis testing through three independent scrambling approaches: distance randomization (100 iterations), phase randomization (20 iterations), and station identity randomization (20 iterations), with statistical significance assessment via permutation testing and z-score analysis
4. 3.3 Methodology Validation: Bias characterization and geometric artifact detection
5. 3.4 Geographic Bias Validation: Continental and hemispheric bias analysis
6. 3.5 Ionospheric Validation: Realistic ionospheric control analysis
7. 3.6 Multi-Band Frequency Analysis: Spectral characterization across 13 frequency bands (10-1500 µHz) with exponential model fitting, frequency specificity assessment, and systematic effects quantification through control bands
8. 3.7 Multiple Comparison Corrections: Comprehensive Bonferroni, FDR, and Family-wise Error Rate corrections across 231 statistical tests

A comprehensive validation suite including cross-validation, null tests, bias characterization, and comprehensive multiple comparison corrections across 231 statistical tests.

Step 4: Advanced Analysis and Visualization

1. 4.0 TEP Advanced Analysis: Model comparison (7 models tested), circular statistics, and advanced correlation metrics
2. 4.1 TEP Visualization: Publication-quality figures for correlation patterns and multi-center consistency
3. 4.2 Synthesis Figure: Comprehensive synthesis of observational evidence
4. 4.3 High-Resolution Astronomical Events: Eclipse and supermoon coherence modulation analysis
5. 4.4 Gravitational Temporal Field Analysis: Multi-window planetary gravitational correlation analysis with timescale strengthening
6. 4.5 Detailed Diurnal Analysis: Hourly temporal analysis with seasonal pattern characterization
7. 4.8 Multi-Band Visualization: Publication-quality figures for spectral analysis, frequency specificity assessment, and gravitational enhancement patterns across 13 frequency bands

Advanced analyses including multi-band spectral characterization, eclipse effects, gravitational correlations, and detailed visualizations.

Pipeline Validation Framework

• Multi-Center Consistency: All analyses performed across three independent GNSS analysis centers (CODE, IGS, ESA) with systematic cross-validation to assess coefficient of variation
• Statistical Rigor: Null testing through data scrambling, cross-validation methods, and multiple comparison corrections applied throughout the pipeline
• Reproducibility: Complete computational provenance tracking, version control, execution logging, and open-source analysis scripts ensure full reproducibility

Implementation Details: For complete usage instructions, configuration options, runtime estimates, and step-by-step execution guide for all 23 pipeline steps, see Section 6 (Analysis Package) at the end of this document. The complete source code and documentation are available at github.com/matthewsmawfield/TEP-GNSS.

Authoritative data sources

• Station coordinates: International Terrestrial Reference Frame 2014 (ITRF2014) via IGS JSON API and BKG services, with mandatory ECEF validation
• Clock products: Official .CLK files from CODE (AIUB FTP), ESA (navigation-office repositories), and IGS (BKG root FTP)
• Quality assurance: Hard-fail policy on missing sources; zero tolerance for synthetic, fallback, or interpolated data

Dataset characteristics

• Data type: Ground station atomic clock correlations
• Temporal coverage: 1 January 2023 to 30 June 2025 (912 days) with date filtering applied
  • IGS: 912 files processed (complete coverage)
  • CODE: 912 files processed (complete coverage)
  • ESA: 912 files processed (complete coverage)
• Spatial coverage: The analysis is based on data from 392 unique GNSS stations. These stations are sourced from the overlapping networks of three analysis centers: CODE (345 stations), IGS (316 stations), and ESA (289 stations). After combining the lists from these centers and filtering for stations with continuous data in the analysis period, the final set of 392 unique stations is obtained. Hemisphere distribution for the 392 analyzed stations is 252 Northern / 115 Southern (68.7%/31.3%, ratio = 2.19), with 367 stations successfully matched to coordinate data. For reference, the master station catalog contains 766 unique stations globally, with hemisphere distribution of 559 Northern / 207 Southern (73%/27%).
• Data volume: 62.7 million station pair cross-spectral measurements
• Analysis centers: CODE (39.0M pairs), ESA (10.8M pairs), IGS (12.9M pairs). Station pair counts vary across centers due to different station network sizes and center-specific data availability and quality criteria.

File counts reflect actual processed files within the 912-day analysis window (1 January 2023 to 30 June 2025) after date filtering.

Data Quality Assurance

A thorough quality validation across 62.73 million station pair measurements demonstrates strong data integrity essential for phase-coherent analysis:

• Filtering efficiency: 99.958% overall data retention across all centers (CODE: 99.944%, IGS: 99.964%, ESA: No removals)
• Phase quality: Healthy boundary clustering (1.62-1.67% of values within ±0.05 rad distance to ±π under wrap) indicates proper phase wrapping without accumulation artifacts. The expected probability for uniform phase distribution is 0.1/(2π) = 1.59%, confirming the observed values are consistent with proper phase processing.
• Temporal completeness: Files present all 912 days; per-day pairs vary by center
• Data integrity: Zero duplicate pairs, 25 stations excluded due to missing coordinates (367/392 retained), zero distance outliers (post-filtering)
• Phase distribution: Uniform coverage across full ±π range enables unbiased phase-alignment index analysis

The boundary clustering metric provides critical validation that phase wrapping is handled correctly—values should distribute uniformly across the ±π range with minimal accumulation at boundaries. The observed 1.62-1.67% boundary clustering matches the expected probability 0.1/(2π) = 1.59% for uniform distribution with ±0.05 rad distance to ±π under wrap, confirming proper phase processing across all analysis centers.

Distance Quality Filtering

Systematic distance outlier removal ensures geodesic calculation accuracy:

• CODE: 21,679 outliers removed (<1 km or >20,000 km) from 39,068,754 initial pairs (0.056%)
• IGS Combined: 4,633 outliers removed from 12,879,380 initial pairs (0.036%)
• ESA Final: 0 outliers removed from 10,809,085 initial pairs (0.000%)

Post-filtering, all analysis centers show zero distance outliers, ensuring clean exponential decay fitting without contamination from geodesic calculation errors or data artifacts.

Temporal Data Density Patterns

Daily pair count variation across the 912-day analysis window differs between centers, reflecting their operational characteristics:

• CODE: CV of daily pair counts = 6.2% (highly consistent: 33,372-48,776 pairs/day)
• IGS Combined: CV of daily pair counts = 58.4% (variable: 990-29,618 pairs/day, reflects multi-center combination strategy)
• ESA Final: CV of daily pair counts = 14.1% (consistent: 10,731-16,290 pairs/day)

The higher temporal variation in IGS Combined reflects its nature as a weighted combination of multiple analysis center solutions, with varying contributor availability over time. Despite this operational difference, IGS achieves R² = 0.966 (primary pooled fit on bin means) and λ consistent with independent centers, demonstrating robustness to temporal sampling patterns.

Core methodology

1. Cross-spectral density computation: For each station pair (i, j), compute complex CSD from clock residual time series
2. Magnitude-weighted phase correlation: Extract phase-alignment index as cos(weighted_phase) where phase is magnitude-weighted circular average
3. Frequency band selection: Analyze 10-500 µHz (periods: 33 minutes to 28 hours) where GNSS clock noise shows characteristic low-frequency behavior
4. Dynamic sampling: Compute actual sampling rate from timestamps (no hardcoded assumptions)

Why magnitude-weighted phase correlation works

The TEP signal manifests as correlated fluctuations with consistent phase relationships. Standard coherency normalizes by amplitude, but GNSS processing suppresses spatially correlated amplitudes precisely when phases are aligned, causing coherency to collapse toward zero even for genuine TEP signals.

Magnitude-weighting preserves signal quality

The approach uses both magnitude and phase information optimally:

• Magnitude weighting: Stronger frequency components (higher |CSD|) get more influence in determining the representative phase through weighted circular averaging
• Phase preservation: Circular statistics correctly handle phase wrapping while incorporating magnitude-based quality weighting
• Amplitude invariance: Final phase-alignment index cos(weighted_phase) is normalized, making it robust to amplitude suppression from GNSS processing

Key insight: Magnitudes are used as quality weights for phase calculation, then amplitude-invariant phase-alignment index is extracted from the weighted phase. This preserves both pieces of CSD information while being robust to processing artifacts.

Physical interpretation of the phase-based approach

For two zero-mean, wide-sense stationary clock residual processes $x_i(t), x_j(t)$, the cross-spectrum $S_{ij}(f)$ is the Fourier transform of the cross-correlation $R_{ij}(\tau)$ (Wiener–Khinchin):

Under TEP, each clock's fractional frequency $y_k(t)$ receives a common field contribution $y_k(t) \propto \phi(\mathbf{x}_k,t)$ plus local noise. In the 10-500 µHz band, any propagation delay across baselines ($\leq 15{,}000$ km) is negligible relative to the periods (33 minutes–28 hours):

Hence, the physically expected inter-station phase is $\approx 0$ in this band; the information lies in how tightly phases cluster, not in a systematic lag. Writing the unit phasor $U_{ij}(f)=S_{ij}(f)/|S_{ij}(f)|$, the metric uses $\mathrm{Re}\{U_{ij}(f)\}= \cos(\arg S_{ij}(f))$. When averaged over pairs within a distance bin, this estimates the circular mean of phases. If the within-bin phase distribution is von Mises $\mathrm{VM}(\mu\!\approx\!0, \kappa(r))$, the expected value is

If the underlying field has exponential spatial covariance, $\mathrm{Cov}[\phi(\mathbf{x}),\phi(\mathbf{x}+\mathbf{r})]\propto e^{-r/\lambda}$, then the concentration $\kappa(r)$ (and thus the circular mean above) inherits an exponential distance-decay, matching the fitted form.

This phase-only approach is robust to amplitude artifacts because it normalizes each $S_{ij}$ to unit magnitude before averaging (amplitude invariance). It distinguishes genuine spatial organization from mathematical artifacts through: (i) comprehensive randomization testing (distance, phase, and station scrambling), which destroys the spatial correlation structure in null tests while preserving it in genuine data; and (ii) replication across independent processing chains (CODE, IGS, ESA) with different systematic vulnerabilities. Standard magnitude-based metrics ($|\mathrm{CSD}|$ or band-averaged real coherency) discard this directional information and therefore cannot detect the distance-structured phase relationships central to TEP.

Model comparison and selection

To validate the theoretical exponential decay assumption, a model comparison is performed using information-theoretic criteria:

• Models tested: Seven correlation functions including Exponential, Squared Exponential (Gaussian RBF), Power Law, Power Law with Cutoff, and Matérn (ν=1.5, 2.5)
• Selection criteria: Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC)
• Methodology: Each model fitted using weighted nonlinear least squares with full uncertainty propagation
• Validation: Cross-center consistency analysis to ensure robust model selection

Exponential model fitting

• Model: $C(r) = A \cdot \exp(-r/\lambda) + C_0$
  • $C(r)$: Mean phase-alignment index at distance $r$
  • $A$: Correlation amplitude at zero distance
  • $\lambda$: Characteristic correlation length (km)
  • $C_0$: Asymptotic correlation offset
• Distance metric: Geodesic distance on WGS-84 (Karney), computed via GeographicLib
• Rationale: For ground-to-ground baselines, geodesic separation tracks propagation-relevant geometry; results are unchanged ($\leq 1$–$2\%$) versus ECEF-chord distances at continental scales
• Distance binning: Logarithmic bins spanning 50 to 13,000 km
• Fitting method: Weighted nonlinear least squares with physical bounds
• Weights: Number of station pairs per distance bin

Uncertainty quantification and independence

• Bin-level bootstrap resampling: 1000 iterations with replacement at distance-bin level, quantifying correlation function fitting uncertainty while preserving network topology
• Station-block bootstrap: 50 iterations systematically removing ~30% of stations per resample, testing robustness to network structure and high-connectivity station bias (yields shifted point estimates due to reduced long-distance coverage)
• Resampling unit: Distance bins for correlation uncertainty; station blocks for structural robustness
• Effective sample size: Accounting for spatial correlations between overlapping pairs
• Independence validation: Station pair non-independence addressed through LOSO cross-validation and block-wise validation (Step 5.5)
• R² interpretation: All reported R² values refer to fits on distance-bin means (N_eff ≈ 25–28), not the 62.7 million individual station pairs
• Random seeds: Sequential 0-999 for reproducibility

Null test validation

• Distance scrambling: Randomize distance labels while preserving correlation values (100 iterations)
• Phase scrambling: Randomize phase relationships while preserving magnitudes (20 iterations)
• Station scrambling: Randomize station assignments within each day (20 iterations)
• Statistical assessment: Permutation p-values computed from null distributions, z-scores for effect size quantification
• Significance threshold: α = 0.05 with correction for multiple testing where appropriate

Dynamic Event Analysis

Methodologically consistent analysis of astronomical events (eclipses, supermoons) using identical phase-alignment index algorithms to test dynamic field responses.

Earth Motion Coupling

Temporal tracking of correlation anisotropy patterns to detect coupling with Earth's orbital motion, rotation, and polar wandering.

Directional Anisotropy Analysis with Distance Distribution Guardrails

Analysis of directional anisotropy in correlation patterns using azimuth-sector decomposition to test for Earth-motion-aligned effects. The λ_EW/λ_NS ratio provides a key diagnostic for rotation-aligned anisotropy.

Distance Distribution Matching Methodology

Critical Guardrail: To prevent bias from differing distance distributions across azimuth sectors, the analysis implements distance distribution matching through two complementary approaches:

• Distance-weighted sector analysis: Each azimuth sector's correlation function is weighted by the inverse of its distance distribution density to ensure equal representation across distance ranges
• Matched-distance subsampling: For each sector, pairs are subsampled to match the global distance distribution, ensuring identical distance sampling patterns across all azimuth sectors

Implementation Details

• Azimuth sectors: Eight 45° sectors (N, NE, E, SE, S, SW, W, NW) with sector assignment based on great-circle azimuth between station pairs
• Distance distribution matching: Global distance histogram computed from all pairs; each sector's pairs weighted or subsampled to match this reference distribution
• Exponential fitting per sector: C(r) = A·exp(-r/λ) + C₀ fitted to each sector's distance-matched data using weighted least squares
• Anisotropy metrics: λ_EW/λ_NS ratio computed from E-W sector average (E, W) versus N-S sector average (N, S)
• Validation: Distance distribution uniformity across sectors verified through Kolmogorov-Smirnov tests (p > 0.05 required for valid analysis)

Bias Prevention

This approach prevents systematic bias in λ_EW/λ_NS ratios that could arise from:

• Geographic sampling bias: Different continents having different azimuth orientations relative to global station network
• Distance sampling bias: Certain azimuth sectors having systematically shorter or longer baseline distributions
• Network topology bias: Station network geometry creating non-uniform distance sampling across azimuth sectors

Result: The distance distribution matching ensures that any observed λ_EW/λ_NS anisotropy reflects genuine directional effects rather than sampling artifacts.

Gravitational Correlation Analysis

Multi-window smoothing analysis correlating planetary gravitational influences with temporal coherence patterns using NASA/JPL ephemeris data.

Planetary Enhancement Factor Definition

The planetary enhancement factor quantifies the observed TEP signal strength relative to gravitational field expectations, providing a mass-independent comparison across different celestial bodies.

Worked Examples: Mercury and Jupiter

Network Coherence Analysis

Investigation of collective network behavior to detect unified detector responses and systematic motion signatures.

Temporal Modulation Studies

High-resolution diurnal and seasonal analysis to detect systematic variations in correlation strength with Earth's orientation and orbital position.

Energy-vs-Velocity Discrimination Analysis

To distinguish between energy-based and velocity-based scaling mechanisms in TEP coupling, a discrimination metric is computed as the simple arithmetic difference of Pearson correlations between gravitational-temporal field patterns and Earth motion parameters:

where $r_E$ is the correlation coefficient for energy-based scaling (gravitational potential energy) and $r_V$ is the correlation coefficient for velocity-based scaling (orbital velocity). This metric quantifies the relative preference between energy and velocity coupling mechanisms, with positive values indicating energy-based scaling preference and negative values indicating velocity-based scaling preference.

Implementation: The discrimination is calculated independently across three analysis centers (CODE, ESA Final, IGS Combined) and then aggregated using bootstrap resampling (1000 iterations) to provide robust confidence intervals. With n=3 centers, the bootstrap CI is mostly illustrative; the analysis reveals discrimination = r_E − r_V = −0.057 (bootstrap 95% CI: [−0.143, +0.030]), indicating no statistically significant preference between energy-based and velocity-based scaling mechanisms, supporting complex multi-mechanism coupling in TEP physics.

TID Exclusion: Ionospheric Contamination Control

Goal: Bound ionospheric contamination by excluding time periods with elevated Traveling Ionospheric Disturbance (TID) activity and quantifying the change in fitted correlation quality.

• Proxy for TIDs: High-resolution temporal metrics (Step 4.3 wavelet and Hilbert instantaneous frequency summaries) aggregated to daily TID activity indices.
• Exclusion rule: Remove days whose TID index exceeds the 75th percentile threshold; retain the remaining days for recomputing mean coherence.
• Effect size: Report percentage change in mean coherence relative to the original (unfiltered) series: Δ% = 100 · (coherence_retained − coherence_original) / coherence_original.
• Operational band context: Assessment pertains to the primary analysis band (10–500 µHz), overlapping canonical TID periods (10–180 minutes), so exclusion uses external temporal structure rather than frequency separation.

Implementation details and file paths: See Appendix 7.2.

Principal Findings

Phase-coherent cross-spectral analysis of 62.7 million quality-filtered station pair measurements from 392 unique GNSS stations reveals systematic distance-dependent correlations in atomic clock residuals. The patterns exhibit substantial multi-center consistency and demonstrate sensitivity to Earth's motion, gravitational fields, and temporal structure. For a detailed breakdown of station counts, see Section 2.2.

These observations suggest systematic sensitivity of global GNSS networks to large-scale spatial structure, Earth's motion, gravitational fields, and temporal dynamics. The diurnal analysis shows patterns with characteristics that could be consistent with certain theoretical predictions. Comprehensive validation demonstrating statistical signal authenticity is presented in Section 4.2.

Summary: This group establishes the fundamental correlation parameters and model validation, demonstrating consistent exponential decay patterns across independent analysis centers and optimal model selection through comprehensive statistical comparison.

Correlation Parameters by Analysis Center

Note: λ values are derived from Leave-One-Station-Out (LOSO) cross-validation for maximum robustness. R² values in this table refer to sensitivity subset fits on distance-bin means (N_eff ≈ 25–28), not individual pairs from the 62.7 million station pair measurements. These sensitivity subset R² values (e.g., 0.78–0.88 for IGS) are distinct from the primary pooled fit R² values (e.g., 0.966 for IGS) reported elsewhere. R² (Binned Fit) reflects the goodness of fit for the exponential model to the binned data, while R² (LOSO CV) represents the model's predictive power on unseen data subsets, providing a more conservative measure of explained variance.

Multi-Center Consistency

• Correlation length ranges: See comprehensive λ + CI Summary Table above for complete values across all validation methods
• Station-block bootstrap robustness: All centers maintain stable λ values within respective bin-level bootstrap CIs despite removing ~30% of stations per resample
• Average λ (unweighted): 3,880 km (within theoretical predictions: 1,000–10,000 km)
• Fit quality (Binned): R² = 0.920–0.970 across all centers (fits to distance-bin means, N_eff ≈ 25–28)
• Processing independence: PPP vs network processing independence demonstrated through station-block bootstrap: λ values remain within respective bin-level bootstrap CIs despite removing ~30% of stations per resample, confirming neither high-connectivity stations nor geographic clustering bias affects the TEP signature.

Model Performance (Akaike Information Criterion)

Key findings: Exponential family models (exponential, Matérn) consistently outperform alternatives. Systematic preference for exponential over squared exponential models (ΔAIC = 5.5-16.6) supports screened scalar field interpretation. All best-fit models achieve R² > 0.92, confirming robust exponential decay characteristics.

Summary: This group examines how correlation patterns depend on environmental factors including elevation, geomagnetic latitude, and directional anisotropy, providing evidence for environmental screening effects and spatial field structure.

Quintile Elevation Stratification Results

Elevation dependence: Correlation length shows systematic variation with elevation quintile (Q1: 3,174 km, Q4: 7,688 km, 2.4× increase), with Q5 showing intermediate values (4,980 km), indicating complex environmental screening effects. The R² values remain consistently high (0.73-0.83) across all quintiles, suggesting robust exponential decay patterns at all elevation ranges. This pattern is consistent with TEP predictions for φ-field coupling through matter density variations.

Geomagnetic-elevation stratification: Combined analysis reveals enhanced correlation structure when both elevation and geomagnetic latitude effects are considered simultaneously, with equatorial high-elevation stations showing the strongest correlations (λ > 5,000 km) and polar sea-level stations showing the shortest (λ < 2,000 km).

Key Anisotropy Observations

• Distance-dependent decay: Clear exponential decay with distance across all centers
• Longitude dependence: Systematic variations with longitude difference (40-80° and 120-160° ranges)
• Multi-center consistency: Reproducible patterns across independent processing strategies
• Intercontinental coherence: Correlation preservation at distances >6,000 km
• Directional variation range: Eight-sector analysis yields typical λ_EW/λ_NS ≈ 0.55–2.03, with peaks up to ~4.75 (IGS Combined center), with robust sample sizes (N = 2.9M-8.3M pairs per sector)

Independence safeguards for directionality

Distance-distribution matching is enforced across azimuth sectors. Sector-level Kolmogorov–Smirnov tests require p > 0.05; all 8/8 sectors pass for each analysis center. Minimum KS p across sectors: CODE 0.563, ESA Final 0.766, IGS Combined 0.705. See Table S1.

Table S1 (inline): Sector-level KS tests on distance distributions confirm that anisotropy estimates are not driven by geometry; complete per-sector p-values are available in analysis outputs.

Summary: This group investigates the coupling between TEP correlations and Earth's various motions, including orbital velocity, rotation, and polar wandering, revealing systematic energy hierarchy relationships and interference patterns (combination and beat frequencies).

Orbital Velocity Correlations

Physical pattern: Negative correlation indicates that higher orbital velocities (perihelion, ~30.3 km/s) correspond to more isotropic correlations, while lower velocities (aphelion, ~29.3 km/s) show stronger directional anisotropy. Individual center significance levels vary, with anisotropy measured as East-West/North-South correlation length ratios from 8-sector directional analysis. Analysis centers use independent processing algorithms and partially overlapping station networks.

Temporal independence: Analysis utilizes 30-day correlation windows with 15-day smoothing spans and effective degrees of freedom N_eff = 28-25 bins across centers, ensuring temporal decorrelation intervals exceed autocorrelation timescales.

Seasonal periodicity: 365.25-day periodicity synchronized with Earth's orbital motion detected across all centers (amplitude variations: 36-55%).

Network Coherence Foundation

Gravitational Energy Hierarchy Validation

Systematic analysis investigates whether TEP detection strength correlates with gravitational energy scales versus kinematic velocities, revealing sophisticated coupling mechanisms involving Earth-Moon gravitational dynamics.

*The ~10²⁸ J gravitational modulation scale represents the effective energy scale of Earth-Moon gravitational field variations induced by Chandler wobble. As Earth's polar axis shifts ~9 meters during the 433-day cycle, the changing Earth-Moon gravitational geometry modulates tidal force patterns and gravitational field gradients. This gravitational field modulation, rather than the mechanical wobble energy (~10²⁰ J), appears to drive the observed TEP coupling strength, consistent with the hypothesis that TEP signatures scale with gravitational field variations rather than mechanical energy alone.

Key Discovery: Moon-Chandler Coupling Association

The unexpectedly strong Chandler wobble TEP signature (R² = 0.377–0.471, |r| = 0.61–0.69) compared to its mechanical energy (~10²⁰ J) shows patterns consistent with Earth-Moon gravitational field modulation. During the 433-day Chandler cycle:

• Polar axis shifts ~9 meters from geographic poles
• Changes Earth-Moon gravitational geometry and tidal force patterns
• Modulates gravitational field gradients at ~10²⁸ J scale
• Associates with stronger TEP coupling than mechanical wobble energy alone would predict

Quantitative sketch (order-of-magnitude): The ~9 m polar wander corresponds to a reorientation of Earth's figure by δθ ≈ 9 m / R_⊕ ≈ 1.4×10⁻⁶ rad. The lunar tidal differential acceleration at Earth's surface is ≈ 1.1×10⁻⁶ m s⁻²; a reorientation by δθ modulates its projection along a fixed station baseline by O(δθ) via a cosine projection, i.e., a fractional change of ~10⁻⁶ in the common-mode driver at 433 days. In the phase framework, the distance-binned phase-alignment index C(r) = 𝔼[cos(arg S_ij)] satisfies C ≈ ½κ in the small-signal regime, with κ ∝ A_common(r)²; therefore a fractional geometric modulation δ produces ΔC/C ≈ 2δ = O(10⁻⁶). This links the 433-day Earth–Moon geometry change to a parts-per-million level modulation of bin-mean phase-alignment index at fixed r—small per bin but statistically resolvable after aggregation—consistent with the detected Chandler-phase correlation and the Chandler±annual sidebands below.

Earth Motion Interference Patterns

All three centers consistently detect four Earth motion interference patterns, demonstrating network sensitivity to the complex interplay of terrestrial rotation, orbital motion, and polar axis wandering:

Energy vs Velocity Scaling Analysis

Methodology: Discrimination computed as the simple arithmetic difference of Pearson correlations: Discrimination = r(energy) - r(velocity), calculated independently across three analysis centers (CODE, ESA Final, IGS Combined) and then averaged.

Scientific Significance: The non-significant discrimination indicates TEP coupling involves gravitational field geometry changes and multi-body dynamics (Earth-Moon system) beyond simple mechanical energy or kinematic velocity scales. This pattern supports sophisticated TEP physics with multiple coupling pathways operating simultaneously, consistent with theoretical predictions of non-integrable time transport mechanisms.

Cross-Planetary Results

*Statistically significant after advanced sham controls (randomized orbital phase assignments) and comprehensive multiple comparison corrections (see Supplement: Multiple Comparison Corrections Summary, Table S1).

Key observation: Venus (4.05 complete orbits) shows strongest correlation, supporting hypothesis that orbital completeness enhances signal detectability in orbital phase analysis.

Summary: This group analyzes correlations between TEP signals and planetary gravitational influences, astronomical events, and mass scaling relationships, revealing complex coupling mechanisms and dynamic field responses.

Multi-Window Timescale Analysis

Gravitational-temporal correlations strengthen systematically with longer analysis windows:

• 31-day window: r = -0.362, p = 1.3×10⁻²⁹
• 59-day window: r = -0.400, p = 2.5×10⁻³⁶
• 91-day window: r = -0.428, p = 7.8×10⁻⁴²
• 119-day window: r = -0.454, p = 1.2×10⁻⁴⁷
• 179-day window: r = -0.477, p = 7.2×10⁻⁵³
• 227-day window: r = -0.503, p = 3.3×10⁻⁴ (autocorr-corrected; raw p = 1.5×10⁻⁵⁹)
• Raw correlation: r = -0.164, p = 6.0×10⁻⁷

Pattern: Correlation strength increases systematically with smoothing window, indicating gravitational-temporal coupling operates on seasonal to annual timescales. Note: The 227-day window p-value is autocorrelation-corrected per §2.3.

Individual Planetary Signatures

• Jupiter: r = -0.256, p = 3.6×10⁻¹⁵ (strongest individual effect)
• Sun: r = -0.230, p = 2.2×10⁻¹² (dominant mass influence)
• Mars: r = -0.111, p = 8.2×10⁻⁴ (opposition cycle effects)
• Venus: Complex distance-dependent coupling
• Saturn: Minimal correlation (incomplete orbital coverage)

Anti-phase seasonal anti-correlations: Higher gravitational influence correlates with lower temporal phase-alignment index, consistent with TEP predictions of non-integrable time transport.

Note: These long-term gravitational correlations (seasonal timescales with multi-month smoothing windows) reflect different physical mechanisms than the short-term opposition event responses examined in Section 3.4.3 (±30-60 day event windows). The anti-correlation pattern observed here operates on annual cycles, while opposition events show transient, center-specific responses to gravitational alignments. Both phenomena may coexist, representing coupling at different timescales.

Solar Eclipse Effects

Analysis of five solar eclipses (2023-2025) using 30-second resolution CLK data reveals significant phase-alignment index modulations across 294,572 station pair measurements. Eclipse effects range from 18-87% of baseline TEP signal strength, representing major network responses to ionospheric disruptions:

Eclipse response patterns: Partial eclipses produce the strongest network phase-alignment index modulation (mean: 4.54×10⁻², 87% of baseline), followed by annular eclipses (mean: 2.48×10⁻², 48% of baseline). A geometry-matched control yields the same ordering, suggesting this hierarchy may reflect TEP correlations responding primarily to ionospheric gradient steepness rather than absolute disruption magnitude. Partial eclipses create broader spatial gradients across the global network, while total eclipses with narrow totality paths produce more localized but less network-wide effects.

Supermoon Perigee Effects

Analysis of 11 supermoon events (2023-2025) reveals coherence modulations:

• CODE: -0.211 ± 0.220%
• IGS: +0.298 ± 0.534%
• ESA: +0.419 ± 0.300%
• Range: -3.17% to +2.04% (ΔR² = 0.85-0.93 relative to control-band baseline at 1000–1500 µHz)

Opposition Event Results Across Analysis Centers

Observed Patterns

• Jupiter oppositions: Strong positive effects (up to +31.86%), particularly in ESA and IGS centers
• Saturn oppositions: Predominantly negative effects (-13.85% to +1.84%)
• Mars opposition: Variable center-specific responses (-14.79% to +6.82%)
• Center variations: Different processing approaches yield different sensitivities to short-term gravitational events

Note: These short-term event responses differ from long-term gravitational correlations (Section 3.4.1), reflecting different physical timescales and mechanisms. The center-specific variations suggest different processing approaches yield different sensitivities to short-term gravitational events, with IGS and ESA showing stronger responses to Jupiter oppositions while CODE exhibits more consistent patterns across planetary events.

Summary: This group provides comprehensive signal characterization across frequency bands, assessment of interference sources, and detailed temporal analysis, establishing the authenticity and nature of the observed correlations.

Key Multi-Band Findings

TID Impact Assessment Results

Temporal Pattern Analysis

High-resolution Hilbert instantaneous frequency analysis reveals systematic temporal effects across multiple frequency bands:

• Solar Rotation (27-day): Moderate effects (0.51-0.60) with significant coupling in IGS Combined (p = 0.009)
• Lunar Month (29-day): Variable response with ESA Final showing anomalously low sensitivity (0.068 vs 0.60-0.93)
• Venus 2f harmonic (112-day): Strongest temporal effects across all centers (1.23-1.73), representing the Venus 2f harmonic (224.7 d / 2 = 112.35 d; observed 112.0 d; 0.3% error). This provides direct evidence for gravitational-temporal field coupling to planetary orbital motion
• Lunar Harmonic (20-day): Consistent coupling (0.12-0.39) corresponding to ~2/3 lunar month (19.86d ≈ 19.69d theoretical), representing lunar tidal harmonics rather than Jupiter-Saturn beat

Key Finding: The 112-day signal is identified as the Venus 2f harmonic (224.7 d / 2 = 112.35 d; observed 112.0 d; 0.3% error), providing quantitative evidence for TEP field coupling to planetary gravitational dynamics. TID exclusion analysis reveals substantial ionospheric contamination with 21-23% coherence improvement potential when excluding high TID activity periods (CODE: 22.36%, ESA Final: 21.28%, IGS Combined: 22.66%), demonstrating significant signal enhancement after ionospheric correction.

Venus 2f Harmonic Identification

The 112-day period represents the Venus 2f harmonic (224.7 d / 2 = 112.35 d; observed 112.0 d; 0.3% error), providing direct evidence for TEP field coupling to planetary gravitational dynamics:

Aliasing Control: Annual and Chandler wobble sideband aliases excluded through Monte-Carlo surrogate phase-scrambling with spectral prewhitening at 365.25-day and 433-day periods before 112-day period estimation.

Physical Significance: The Venus 2f harmonic identification transforms what initially appeared as a "mysterious" signal into quantitative evidence for gravitational-temporal field coupling. The 0.3% error between observed (112.0d) and predicted (112.35d) periods provides direct validation of TEP field response to planetary orbital motion. The strongest temporal effects (1.23-1.73 days standard deviation difference between event-locked and background periods) occurring at this frequency suggest Venus orbital dynamics modulate the temporal field structure with measurable harmonic resonance.

Cross-Validation: This finding aligns with Step 3.6 multiband analysis showing the 112-day period as dominant across all analysis centers, and supports the non-ionospheric origin through TID exclusion analysis. The systematic ranking of temporal effects (ESA Final < CODE < IGS Combined) correlates with TID impact scores, suggesting Venus coupling strength varies systematically across analysis centers.

Day-Night Coupling Strength

Systematic day-night variations in correlation strength indicate coupling to Earth's rotation and orbital position:

• CODE: 1.019× stronger during day (1.9% enhancement) - global-mean ratio
• IGS: 1.076× stronger during day (7.6% enhancement) - global-mean ratio
• ESA: 1.060× stronger during day (6.0% enhancement) - global-mean ratio

Diurnal Pattern Characteristics

Temporal analysis scope: Comprehensive hourly temporal analysis across three independent analysis centers reveals systematic diurnal variations consistent with φ-field coupling predictions:

Key Finding: ESA Final demonstrates enhanced φ-field sensitivity with 2× greater coupling strength compared to CODE, while IGS Combined exhibits the strongest diurnal modulation (7.6% day-night variation), indicating center-dependent sensitivity to temporal field dynamics.

Convergent Lines of Evidence

The comprehensive analysis reveals seven independent observational domains supporting systematic network sensitivity to spacetime structure:

1. Baseline spatial correlations: Exponential decay with a primary correlation length of 3,330–4,549 km, consistent across independent processing chains, and bootstrap validation ranges of 2,986–3,438 km (CODE), 2,346–2,757 km (ESA Final), and 2,616–2,851 km (IGS Combined)
2. Spectral characterization: Broadband correlation structure (R² > 0.85 from 10-200 µHz) with gravitational enhancement in tidal bands (λ = 4,677 km) but persistent post-tidal signals (30-40 µHz shows R² = 0.946), suggesting universal coupling rather than frequency-selective contamination
3. Environmental modulation: Systematic variations with elevation and geomagnetic latitude indicating screening effects
4. Earth motion coupling: Correlations with orbital velocity, Chandler wobble, and interference patterns (combination and beat frequencies)
5. Network coherence: Moderate unified detector behavior (0.493-0.560) with strong base mesh coherence (0.643-0.644) and comprehensive multi-frequency coupling patterns
6. Gravitational correlations: Planetary influences with significant timescale strengthening (31-day to 227-day)
7. Temporal and dynamic modulation: Diurnal, seasonal, eclipse, and supermoon responses

Theoretical context: These observations find natural interpretation within the Temporal Equivalence Principle, which predicts scalar field variations should couple to atomic transition frequencies through conformal metric structure g̃_μν = A(φ)g_μν + B(φ)∇_μφ∇_νφ.

The convergence of multiple independent observational domains, combined with consistent reproduction across analysis centers, indicates that global GNSS networks exhibit sensitivity to large-scale spacetime structure at previously unexplored scales.

Multi-Layered Validation Approach

1. Null Hypothesis Testing: Statistical Validation

Test Design: Three independent scrambling approaches applied systematically across all analysis centers to validate signal authenticity through controlled randomization:

• Distance scrambling: Randomize station pair distances while preserving phase relationships (100 iterations)
• Phase scrambling: Randomize phase values while preserving distance structure (20 iterations)
• Station scrambling: Randomize station assignments within temporal blocks (20 iterations)

Statistical Significance Assessment: Comprehensive multiple comparison correction analysis across all 231 statistical tests demonstrates robust statistical findings:

Comprehensive Analysis (231 statistical tests across 18 families):

• Major test families: Multiband analysis (40), Band diagnostics (40), Anisotropy orbital (24), Model comparisons (18), Advanced analysis (15), Eclipse analysis (15), Bootstrap cross-method (12), Hilbert-IF astronomical (12)
• Bonferroni correction: 137/231 tests remain significant (59.3%) with ultra-stringent α = 0.000216 (0.05/231 tests)
• Benjamini-Hochberg correction: 168/231 tests remain significant (72.7%) with FDR q = 0.05
• Family-wise correction: 154/231 tests remain significant (66.7%) across analysis families
• Primary TEP findings: All 3 primary exponential fits maintain significance after all corrections
• Complete pipeline validation: 231 tests span every aspect from data quality to astronomical events, including null hypothesis validation and bootstrap consistency checks

Interpretation: The comprehensive multiple comparison correction analysis across 231 statistical tests provides strong evidence for theory-predicted signal authenticity. The analysis tests the established Temporal Equivalence Principle theoretical framework, which predicted exponential correlation decay with characteristic length scales before data analysis. With comprehensive validation spanning data quality metrics, null hypothesis testing, bootstrap validation, band-specific diagnostics, astronomical period analysis, and complete pipeline verification, this represents rigorous statistical validation of theory-driven predictions. The 59-73% survival rate after ultra-conservative corrections (α = 0.000216) indicates that the core TEP theoretical predictions are statistically robust across 18 independent validation approaches, though alternative explanations remain to be fully excluded.

2. Cross-Center Consistency: Independent Processing Chains

Test: Compare results across CODE, IGS Combined, and ESA Final using different algorithms, software implementations, and station network configurations.

Result: λ = 3,330-4,550 km with CV of correlation lengths across centers = 18.2% consistency across fundamentally different processing approaches.

Interpretation: The convergence of independent processing methodologies on consistent physical parameters provides strong evidence against center-specific systematic artifacts. While all centers analyze the same underlying IGS network data, their different software implementations, quality control procedures, and solution strategies would be expected to produce different artifacts if the signal were methodological in origin. The observed consistency supports a genuine physical phenomenon.

3. Multi-Band Frequency Specificity: Universal Coupling Discrimination

Test: Analyze correlation patterns across 13 independent frequency bands (10-1500 µHz, 506M pairs for CODE) to distinguish universal broadband coupling from frequency-selective mechanisms including classical tidal contamination.

Result: Cross-center validation demonstrates remarkable consistency - ESA Final (R² = 0.970), CODE (R² = 0.920), and IGS Combined (R² = 0.966) achieve nearly identical optimal band performance, eliminating systematic biases. Broadband correlation structure with R² > 0.85 maintained from 10-200 µHz (CODE shows CV of R² across bands = 2.9% across this range). Post-tidal 30-40 µHz band exhibits R² = 0.946 (mean across centers), ranking among strongest bands rather than showing expected drop-off. Enhancement ratios of 1.26-1.90× (mean: 1.58×) across all centers fall well below the >3-5× threshold characteristic of frequency-selective tidal contamination. Control bands (1000-1500 µHz) show R² = 0.618 (mean), indicating reduced model fit quality (ΔR² ≈ 0.33) compared to TEP band.

Interpretation: The cross-center validation achievement, combined with frequency analysis, provides compelling evidence consistent with universal conformal coupling rather than frequency-selective contamination. The remarkable consistency across independent processing methodologies - ESA Final, CODE, and IGS Combined achieving nearly identical multiband patterns - eliminates center-specific systematic biases and represents a crucial scientific validation milestone. The absence of sharp spectral features characteristic of classical tidal contamination, combined with persistent post-tidal signals (R² = 0.946) and modest enhancement ratios (1.58× vs. expected >3-5× for tidal contamination), supports the TEP interpretation. Tidal band gravitational enhancement (λ = 4,627 km mean) indicates φ-field response to gravitational gradients, while broadband persistence suggests coupling operates universally across frequency scales. The quantitative separation between TEP (R² ≈ 0.95) and control (R² ≈ 0.618) bands establishes robust discrimination criteria for physical correlations vs. systematic effects.

4. Signal-Artifact Discrimination: Null Test Validation

Test: Apply comprehensive null tests using randomized distance, phase, and station permutations to quantify systematic artifact thresholds.

Result: Substantial separation between real TEP signal (R² = 0.920-0.970) and null test artifacts (R² = 0.03-0.06) in the 10-500 µHz frequency band, representing ΔR² = 0.86-0.94. Control band analysis shows ΔR² = 0.334 between TEP and control frequencies (TEP R² = 0.952 vs control R² = 0.618 at 1000-1500 µHz, mean across centers).

Interpretation: The substantial separation between genuine TEP correlations and systematic artifacts provides strong evidence against processing contamination. Null tests establish quantitative thresholds: artifacts consistently produce R² ≤ 0.06, while TEP signals achieve R² ≥ 0.920. The frequency-dependent separation (stronger at TEP frequencies, weaker at control frequencies) demonstrates signal authenticity rather than broadband systematic effects. This discrimination validates the phase-coherent methodology and confirms that observed correlations represent genuine field-structured coupling rather than common-mode processing artifacts.

5. Geometric Controls: Network Topology Cannot Explain the Pattern

Test: Apply identical analysis to synthetic data with same station geometry using four distinct scenarios: uniform noise, Gaussian noise, distance-independent signals, and weak geometric patterns.

Result: Maximum synthetic R² = 0.068 across all scenarios compared to observed R² = 0.92-0.97, representing ΔR² = 0.85-0.90. All synthetic scenarios produce R² ≤ 0.07, providing clear discrimination from real TEP signals.

Interpretation: Station distribution and network geometry are insufficient to explain the observed correlation patterns. The substantial difference in model performance (ΔR² = 0.85-0.90) provides strong evidence that geometric artifacts do not account for TEP correlations. This synthetic data analysis establishes quantitative baselines that must be exceeded for genuine physical effects, which the TEP signals consistently achieve.

6. Ionospheric Independence: Bounding Contributions with External Proxies

Test: To robustly bound ionospheric contributions, a two-stage approach is employed. First, the TEP signal is correlated with external space weather indices (Kp, F10.7) to test for direct coupling. Second, a high-resolution TID exclusion analysis is applied to quantify the remaining, non-linear contributions.

Result: Direct correlation with space weather indices is weak (r = 0.12-0.13, p > 0.29), indicating no primary linear relationship. The subsequent TID exclusion analysis, designed to capture more complex interactions, reveals substantial ionospheric contamination with 21-23% coherence improvement potential when excluding high TID activity periods. This analysis further identified the Venus 2f harmonic (see Section 3.3.2), providing direct evidence for gravitational coupling.

Interpretation: Ionospheric contamination represents a significant component (~22%) of the observed signal variability, with substantial improvement potential through TID exclusion. The weak proxy correlations with space weather indices, combined with the systematic improvement found via TID analysis, indicate complex ionospheric coupling that can be mitigated through temporal filtering. The clear detection of a planetary orbital harmonic further strengthens the interpretation of gravitational-temporal coupling.

7. Distribution Neutrality: No Statistical Bias

Test: Apply equal-count binning to eliminate right-skewed distance distribution effects.

Result: 94.8-97.9% signal preservation (average 96.7%) when switching from logarithmic to equal-count binning. R² differences remain small (0.021-0.052), with equal-count binning producing R² = 0.933 (CODE), 0.894 (IGS), 0.914 (ESA), demonstrating robust exponential correlation independent of distance distribution effects.

Interpretation: The correlations represent genuine physical signals, not statistical artifacts of the distance distribution.

8. Bias Characterization: Clear Discrimination Thresholds

Test: Generate realistic GNSS noise scenarios and quantify methodological bias through random data simulation (50 iterations).

Result: Clear signal-to-bias separation (TEP R² = 0.953 vs maximum artifact R² = 0.050, ΔR² = 0.903). Random data simulations confirm the bias envelope is well below the TEP signal strength.

Interpretation: Clear quantitative thresholds distinguish genuine signals from artifacts. The discrimination factor provides a robust statistical foundation for signal authenticity claims.

9. Binning & Weighting Sensitivity: Methodological Robustness

Test: To ensure the results are not artifacts of specific methodological choices, a sensitivity analysis was performed as part of the main analysis pipeline (Step 4.0). The analysis was repeated by sweeping three independent parameters: the number of distance bins attempted (25, 40, 80), the binning strategy (logarithmic, linear, quantile/equal-count), and the weighting scheme used in the exponential fit (weighted by pair count, by sqrt(pair count), or unweighted).

Result: The key parameters (λ and R²) demonstrate good stability. The correlation length λ remains within a tight cluster (~4350-4450 km) across different bin counts and strategies, with R² consistently exceeding 0.91. The analysis confirms that weighting the fit by the number of pairs per bin is appropriate, but the result is not highly sensitive to the specific weighting function (count vs. sqrt(count)).

Interpretation: The stability of the results across a wide range of analysis parameters confirms that the detected exponential correlation is a fundamental property of the data, not an artifact of the chosen methodology. This provides strong evidence against concerns of methodological bias.

10. Phase Distribution Quality: Proper Phase Processing Validation

Test: Validate phase boundary handling through boundary clustering analysis to detect phase wrapping artifacts or accumulation errors.

Methodology: For proper phase processing, plateau_phase values should distribute uniformly across the ±π range. Excessive boundary clustering (values accumulating near ±π limits) would indicate phase wrapping artifacts or numerical errors. Boundary clustering is quantified as the percentage of phase values within ±0.05 radians distance to ±π under wrap.

Result: Boundary clustering rates of 1.62-1.67% across all centers match theoretical expectations for uniform phase distribution (expected ~1.6% for ±0.05 rad distance to ±π under wrap on uniform distribution), with phase values distributing evenly across all octants of the phase circle.

Interpretation: The healthy boundary clustering rates (1.62-1.67%, matching theoretical expectations) combined with balanced ±π distribution validate proper phase processing across all centers. This confirms that the phase-alignment index metric operates on properly wrapped phase values without accumulation artifacts, providing an essential foundation for phase-coherent correlation analysis. The multi-center consistency (CV of boundary clustering rates = 1.5%) further validates methodological robustness.

11. Volume Independence: Robustness Across Data Densities

Test: Assess whether correlation parameters depend on data volume by comparing centers with vastly different pair counts.

Data volumes: CODE: 39.0M pairs, IGS: 12.9M pairs, ESA: 10.8M pairs (3.6× volume ratio between largest and smallest)

Result: Despite 3.6-fold difference in data volume, all centers demonstrate consistent elevation dependence patterns and geomagnetic stratification with systematic λ ranges (CODE: 3,225–7,499 km, IGS: 2,616–5,453 km, ESA: 1,600–3,914 km). Primary pooled fits show R² = 0.92–0.97 (distance-bin means, N_eff ≈ 25–28 bins used from 40 attempted). Sensitivity subsets (elevation/geomagnetic): R² = 0.70–0.91. This range (1,600–7,500 km) is a result of sensitivity analysis and is distinct from the primary finding.

Interpretation: TEP correlations are robust across different statistical sampling densities, ruling out sample-size-dependent artifacts. The consistent elevation dependence and geomagnetic stratification patterns across all centers demonstrate signal authenticity. ESA Final achieves excellent fits (R² = 0.81–0.91) despite having the smallest dataset (10.8M pairs), while CODE shows systematic elevation trends across the largest dataset (39.0M pairs).

Validation Score: 11/11 criteria passed with 100% validation achievement. The comprehensive multiple comparison correction analysis across 231 statistical tests supports statistical robustness, with 59-73% of findings surviving ultra-conservative corrections across 18 validation families including complete data quality, null hypothesis testing, bootstrap validation, and band diagnostic verification.

Bootstrap Confidence Intervals

Bin-level bootstrap method: 1000 bootstrap iterations with replacement at distance-bin level, preserving pair count weights and spatial structure. Station-block bootstrap method: 50 iterations systematically removing ~30% of stations per resample to test robustness to network structure.

Complete λ and CI values: See Table in §3.1.1 for comprehensive correlation length estimates and confidence intervals across all validation methods and analysis centers.

Interpretation: Bin-level bootstrap CIs reflect correlation function uncertainty; station-block bootstrap CIs test robustness to station removal and are narrower because they assess stability under station subset removal rather than full correlation function uncertainty.

Jackknife Resampling Validation

Method: Systematic leave-one-out resampling to assess parameter stability across geographic subsets.

• Jackknife estimates: See Table in §3.1.1 for complete values with CV of λ across subsets metrics

Convergence: Bootstrap and jackknife estimates converge with <5% difference, confirming statistical robustness.

Cross-Validation Results

Leave-One-Station-Out (LOSO): Removes all pairs involving each station to test robustness to network structure and station-specific artifacts.

• Parameter stability: λ estimates remain within bootstrap confidence intervals across LOSO iterations
• Model fit preservation: R² remains > 0.90 for all LOSO subsets
• Maximum deviation: <8% from full-network estimates

Leave-One-Day-Out (LODO): Tests temporal stability by systematically removing each day from the 912-day dataset.

• Temporal consistency: λ stable across all LODO iterations (CV of λ over time < 5%)
• No day-specific artifacts: No single day drives the observed correlations

Block-wise Validation: Leave-N-stations-out blocks provide additional independence testing for regional biases.

• Monthly blocks: Consistent λ across all monthly subsets
• Spatial blocks: Continental-scale resampling preserves correlation structure
• 100% cross-validation convergence: All validation approaches yield consistent parameters

Statistical Independence Considerations

Pair-level dependencies: Station pairs sharing common stations create covariance structures that could inflate precision estimates. This is addressed through:

• Distance-bin aggregation: Primary analysis operates on binned means rather than individual pairs, reducing dependency effects
• LOSO validation: Removes all pairs involving each station, directly testing impact of shared-station dependencies
• Effective N estimation: Bootstrap confidence intervals reflect ~25-28 independent bins, not 62.8M individual pairs
• Conservative uncertainty: CIs properly account for spatial correlation structure in station network

Statistical Independence and Effective Degrees of Freedom

A primary statistical challenge in network analysis is addressing the inherent covariance of station pairs that share a common station, which can potentially inflate precision estimates when weighting by pair counts. This study's validation framework was designed to systematically address this through a combination of data aggregation and rigorous empirical testing, rather than explicit model-based covariance fitting.

• Distance-Bin Aggregation: The primary analysis operates on binned means rather than individual pairs. This standard geostatistical technique reduces pair-level dependency effects by averaging thousands of pairs within each spatial bin.
• Leave-One-Station-Out (LOSO) Validation: This is the critical empirical test of pair-dependence. By systematically removing *all pairs* associated with a given station and refitting the model, LOSO directly probes the influence of high-degree stations. The observed parameter stability (λ deviation <8%, R² > 0.90 across iterations) confirms that the result is not driven by the over-representation of any single station or its associated non-independent pairs.
• Station-Block Bootstrap Validation: Beyond single-station removal, systematic removal of multiple stations simultaneously (50 bootstrap iterations removing ~30% of stations per resample) demonstrates robust λ stability within bootstrap confidence intervals across all processing methodologies. This addresses the specific concern about high-degree station bias by showing the signal persists even when removing entire subnetworks.
• Bootstrap on Bins, Not Pairs: The bootstrap resampling is performed on distance bins, not the 62.7M individual pairs. This correctly estimates the uncertainty of the *correlation function itself*, yielding conservative confidence intervals. **All statistical claims in this work are based on bin-level analysis, not the total number of station pairs.**

Assessment: While formal geostatistical or mixed-effects modeling represents a valid alternative for future work on refining parameter precision, the comprehensive empirical framework employed here—particularly the stability under LOSO cross-validation—provides high confidence that the detected physical effect is robust and not an artifact of pair-level covariance. The methods used provide a powerful, data-driven assessment of the effective degrees of freedom.

Geographic Consistency Assessment

Methodology: Statistical resampling validation (jackknife) using the complete dataset to assess geographic consistency, hemisphere balance, elevation independence, and ocean vs. land baseline characteristics across three independent analysis centers.

Key Findings:

• Jackknife cross-validation: See Table in §3.1.1 for complete values with CV metrics
• Hemisphere balance: 69.4% North, 30.6% South (272/120 stations, ratio = 2.27, moderate bias addressed through subset validation)
• Balanced subset validation: 10 subsets with perfect 50:50 hemisphere balance show consistent correlation parameters, proving hemispheric imbalance does not create spurious correlations
• Elevation independence: Consistent correlation strength (0.74-0.88) across all topographic bands
• Ocean vs. land: Ocean baselines show slightly lower coherence (-5.8%), indicating no artificial enhancement from geographic bias
• Validation confidence: 0.813/1.0 (HIGH confidence with 62.7M measurements)

Eclipse Scale Consistency

• Eclipse shadow scale: Direct solar blockage spans ~2,000–3,000 km diameter
• Observed effect extent: Coherence modulations observed to distances matching TEP λ = 3,330–4,549 km
• Cross-center consistency: Eclipse type hierarchy (Partial > Annular > Total) observed across independent centers; a geometry-matched control yields the same ordering
• Limitations: Five eclipses with 1-2 events per type provide preliminary evidence requiring independent replication

Opposition Event Scale Consistency

• Temporal scope: Measures short-term coherence changes (±30-60 day windows), distinct from long-term correlations (Section 3.4.1)
• Center-specific responses: Varying sensitivities (CODE: -14.79% to +0.24%, IGS: up to +27.71%, ESA: up to +31.86%) reflect different processing approaches to short-term gravitational perturbations
• Physical interpretation: Short-term event enhancements can coexist with long-term anti-phase coupling, representing different physical timescales

Acknowledged Shared Systematic Dependencies

Methodological Context: True independence is inherently limited by shared systematic components across all analysis centers:

• Common data source: All centers process the same underlying IGS global tracking network observations
• Shared physical models: Similar satellite orbit determination using JPL ephemeris (DE440/DE441), common Earth rotation parameters (IERS), and shared reference frame realizations (ITRF2014)
• Fundamental processing constraints: All centers apply similar GNSS processing principles, multipath mitigation strategies, and environmental correction models
• Common calibration standards: Shared atomic frequency standards, similar relativistic corrections, and coordinated UTC realizations

Critical assessment: This shared infrastructure could theoretically produce correlated systematic errors across centers, representing a fundamental challenge in GNSS-based validation studies.

Evidence for Processing Independence Despite Shared Infrastructure

The analysis was designed to provide multiple lines of evidence for genuine signal detection despite these systematic dependencies:

1. Fundamentally Different Processing Philosophies

• CODE (Network Solution): Uses double-differenced observations with global network constraints, inherently coupling all stations through relative measurements and unified network adjustment
• ESA (Precise Point Positioning): Processes each station independently using undifferenced observations, treating each receiver as isolated without explicit inter-station connectivity
• IGS (Multi-Center Combination): Weighted meta-analysis combining diverse processing approaches, including both network and PPP solutions

Key Insight: Network solutions and PPP represent philosophically opposite approaches to the systematic dependencies that could create artifacts. Network solutions would amplify inter-station artifacts through explicit connectivity, while PPP would suppress them through station isolation. The convergence on λ = 3,330-4,549 km (CV of λ across centers = 18.2%) across these opposite systematic vulnerabilities provides strong evidence for genuine physical phenomena.

Robust Control Test: To directly test the independence of PPP vs network processing, station-block bootstrap analysis was performed that systematically removes random subsets of stations (up to ~30% per resample). Despite removing potential high-connectivity hub stations and entire geographic regions, λ values remain stable within their respective bin-level bootstrap confidence intervals (see Table in §3.1.1). The station-block bootstrap CIs are narrower and shifted compared to the primary bin-level bootstrap CIs because they test robustness to station removal rather than correlation function uncertainty. This demonstrates that neither network connectivity effects nor geographic clustering biases can account for the observed TEP signature.

2. Signal Survival Through Processing Suppression

A key indicator of authenticity: GNSS processing is explicitly designed to remove spatially correlated signals through network constraints, environmental corrections, and quality control filtering. The persistence of strong correlations (R² = 0.920-0.970) through processing designed to eliminate such signals provides substantial evidence for signal robustness rather than processing artifacts.

3. Comprehensive Systematic Validation

• Null hypothesis testing: Substantial signal separation (ΔR² = 0.86-0.94) over randomized controls provides strong evidence of spatial and temporal structure beyond processing artifacts
• Cross-validation robustness: LOSO/LODO validation with consistent parameters across geographic and temporal subsets
• Physical dependencies: Systematic variations with elevation, orbital velocity, and gravitational fields represent characteristics unexpected from shared processing algorithms

4. Systematic Environmental Dependencies

A key physical discriminator: A processing artifact arising from shared models would not be expected to show systematic dependencies on the local physical environment of individual stations. However, the advanced analysis reveals clear, physically plausible patterns:

• Elevation Stratification: The correlation length (λ) shows a monotonic increase with station elevation quintiles, from λ ≈ 1,600–3,200 km at the lowest elevations to λ ≈ 3,900–7,500 km at the highest. This is consistent with a screened field where atmospheric density modulates the coupling strength.
• Geomagnetic Stratification: The analysis reveals distinct correlation regimes based on geomagnetic latitude, with equatorial regions (λ ≈ 1,760–2,080 km) showing different characteristics from polar regions (λ ≈ 1,300–3,400 km).

This coupling to the local physical and electromagnetic environment is a signature of a genuine physical phenomenon interacting with the Earth system, not an artifact of data processing algorithms.

5. Coherence with External Physical Phenomena

The key discriminator: Perhaps the most compelling evidence against processing artifacts is the signal's coherence with independent, external physical phenomena that are not inputs to the GNSS processing chain. A self-contained systematic error would be blind to these external drivers. The analysis reveals multiple such correlations:

• Planetary Ephemeris: The detection of 11 significant coherence modulations (σ > 3.0) time-locked to planetary events—such as the 5.98σ Saturn opposition—links the signal to the gravitational dynamics of the solar system.
• Geophysical Oscillations: The significant detection of the 433-day Chandler wobble (coefficient of determination R² = 0.377–0.471, all p<0.01) demonstrates that the signal is coupled to the physical motion of the Earth itself.
• Orbital Dynamics: The robust negative correlation between anisotropy ratios and Earth's orbital velocity (Fisher's meta-analysis: combined p = 2.7×10⁻¹⁴ unadjusted or 1.4×10⁻¹⁰ dependence-adjusted, Q = 1.85, I² = 0%; p-values autocorrelation-corrected via Bretherton; N_eff ≈ 47) provides evidence of coupling to the motion of the entire Earth system through the solar system.
• Gravitational coupling: Planetary correlations with timescale strengthening, Moon-Chandler coupling, and dynamic astronomical responses
• Temporal structure: Diurnal and seasonal modulation with systematic patterns

The signal's phase-locking with these grand, independent "clocks" of the solar system provides a final, powerful line of evidence that the observed phenomenon is physical in origin, not an artifact of the measurement system.

Residual Systematic Risk Assessment

While no analysis can completely eliminate the possibility of sophisticated shared systematic errors, several observations from this study's design support genuine signal detection:

1. Opposing systematic vulnerabilities converge: Network and PPP processing should amplify different artifact types, yet yield consistent results
2. Processing suppression paradox: Signal persistence despite suppression designed to remove correlated signals
3. Physical validation: Earth motion coupling, planetary correlations, and astronomical event responses validate signal authenticity through independent physical mechanisms
4. Multi-modal confirmation: Convergent evidence across spatial correlations, temporal dynamics, and gravitational coupling provides systematic validation beyond processing effects

Assessment: While shared systematic dependencies represent a legitimate and acknowledged limitation, the convergence of opposing processing philosophies on consistent physical parameters, combined with comprehensive validation through independent physical mechanisms, provides substantial evidence for genuine field coupling phenomena transcending processing artifacts.

Predicted Spectral Signatures

Classical tidal contamination and TEP universal coupling make distinct, testable predictions for the frequency dependence of correlations:

Observational Evidence

1. Post-tidal band persistence: The 30-40 µHz frequency band—immediately beyond primary tidal forcing frequencies—exhibits mean R² = 0.946 across three centers, ranking among the strongest bands (CODE: 1st of 13, IGS: 4th of 13, ESA: 3rd of 14). This value equals or exceeds tidal band correlations (mean R² = 0.941), appearing inconsistent with classical tidal contamination which would predict weakest correlations in this transition region.

2. Smooth spectral response: Analysis across 13 bands reveals gradual decline in correlation strength from 10-1500 µHz without sharp transitions. The CODE center demonstrates strong spectral uniformity with CV of R² across bands = 2.9% across 10-200 µHz, maintaining R² > 0.85 throughout this 20-fold frequency range. This smooth response appears incompatible with frequency-selective tidal mechanisms.

3. Enhancement ratio analysis with theoretical justification: Observed ratios of 1.26-1.90× (mean: 1.58×) between TEP and control bands fall well below the >3× threshold expected for classical tidal dominance. This threshold is theoretically grounded in several principles:

The observed modest enhancement ratios (1.58× mean) represent universal characteristics: tidal (1.52×), TEP (1.54×), and post-tidal (1.53×) bands show statistically indistinguishable enhancement levels. This pattern contradicts frequency-selective tidal contamination, which would exhibit sharp spectral features, and supports universal coupling mechanisms with gravitational modulation operating across all frequency scales.

4. Spatial scale transitions: While correlation lengths decrease from λ = 4,677 km (tidal bands) to λ = 1,502 km (post-tidal), the fit quality (R²) remains consistently high (>0.85). This decoupling of spatial scale from correlation strength indicates that tidal frequencies couple to larger-scale gravitational gradients while maintaining the same underlying exponential decay mechanism—a pattern more consistent with scalar field response to varying gravitational forcing than with classical tidal residuals.

Theoretical Interpretation

Within the TEP framework, tidal effects represent gravitational field gradients that modulate the φ-field through disformal coupling B(φ)∇_μφ∇_νφ. The observed pattern—longest correlation lengths at tidal frequencies but persistent correlations across all bands—suggests the φ-field responds to gravitational forcing at multiple scales rather than being contaminated by classical tidal residuals. The universal conformal term A(φ) provides broadband coupling, while the disformal term enhances response to large-scale gravitational gradients at tidal frequencies.

Conclusion: The multi-band analysis appears to exclude classical tidal contamination as the dominant mechanism through four independent lines of evidence: post-tidal signal persistence, smooth spectral response, modest enhancement ratios, and spatial scale decoupling from correlation strength. While tidal gravitational effects clearly influence the correlations (evidenced by longest λ at tidal frequencies), the broadband nature and persistent post-tidal signals suggest these effects manifest through universal field coupling rather than frequency-selective contamination.

Frequency Band Overlap: Alternative Discriminators Required

Important limitation: TIDs exhibit periods of 10-180 minutes (92-1,667 µHz), directly overlapping the analysis band (10-500 µHz). Temporal/frequency separation cannot exclude TID contamination. The analysis therefore relies on spatial structure analysis, processing pipeline evidence, and ionospheric independence validation to distinguish these phenomena.

Spatial Structure Incompatibility

• TIDs: Coherent plane-wave propagation with defined k-vectors and wavelengths of 100–3,000 km, creating directional propagation patterns
• Observed signals: Isotropic exponential correlation decay with screening length λ = 3,330-4,549 km, inconsistent with directional wave propagation
• Model discrimination: Plane-wave models would produce correlation patterns dependent on station pair azimuth; observed patterns show azimuthal independence (Section 3.2)
• Different physics: Wave propagation (phase-coherent traveling disturbances) vs field screening (exponential correlation decay)

Processing Pipeline Mitigation

GNSS analysis centers apply comprehensive ionospheric corrections including dual-frequency combinations, global ionospheric maps (GIMs), and common-mode signal removal specifically designed to mitigate TID effects. The persistence of strong correlations (R² = 0.920-0.970) after these corrections indicates the signal originates below the ionosphere, in the atomic clock frequencies themselves. Multi-center λ consistency (CV of λ across centers = 18.2%) across different correction approaches further supports non-ionospheric origin.

Ionospheric Independence Validation

Direct correlation with space weather indices (Kp, F10.7 flux) shows weak relationships (r = 0.12-0.13, p > 0.29), well below thresholds for ionospheric contamination. TIDs are strongly modulated by solar activity and geomagnetic conditions; the observed signal's independence from these drivers provides additional evidence against TID origin.

Conclusion: While TIDs operate in an overlapping frequency band, the observed signals are distinguished through (1) spatial structure incompatibility (plane-wave vs exponential decay), (2) survival through aggressive ionospheric processing, (3) independence from space weather drivers, and (4) multi-center convergence on identical correlation lengths. These discriminators provide strong evidence against TID contamination.

Frequency Band Incompatibility

• TEP signals: L-band GNSS frequencies (1.2–1.6 GHz)
• TEQ: VHF/UHF amateur bands (30–300 MHz)
• Frequency separation: order-of-magnitude frequency mismatch (GNSS L-band ~1.2–1.6 GHz vs TEQ VHF/UHF ~30–300 MHz)

Temporal Characteristics Mismatch

• TEP signals: Continuous over months/years (persistent correlations)
• TEQ: Hours post-sunset duration (transient propagation)
• Geographic scope: Global vs regional propagation patterns

Conclusion: Trans-equatorial propagation (TEQ) is excluded as an explanation for the observed Global Time Echo correlations. The frequency band incompatibility (order-of-magnitude mismatch: GNSS L-band 1.2–1.6 GHz vs TEQ VHF/UHF 30–300 MHz) and temporal persistence mismatch (continuous multi-year TEP signals vs transient hours-duration TEQ propagation) provide fundamental physical exclusion criteria that are independent of statistical analysis.

1. Atmospheric Loading Effects: Systematic Analysis

Mechanism hypothesis: Large-scale atmospheric pressure variations could create correlated timing effects through gravitational loading, where atmospheric mass redistribution affects local gravity and influences atomic clock rates through gravitational redshift (δν/ν ≈ δg/g ≈ 10⁻⁹ for mb-scale pressure changes).

Critical discriminators against atmospheric loading:

• Processing immunity: GNSS analysis centers already apply comprehensive atmospheric loading corrections (Global Pressure and Temperature models, Vienna Mapping Functions), yet strong correlations persist (R² = 0.920-0.970)
• Temporal persistence: Atmospheric loading effects operate on synoptic timescales (3-10 days), incompatible with the multi-year correlation persistence observed
• Orbital velocity coupling: The observed negative correlation with Earth's orbital speed (r = -0.701 to -0.796) represents a signature absent from atmospheric loading mechanisms
• Planetary correlations: Detection of 8 Bonferroni-significant astronomical events with mass-scaled enhancement patterns contradicts terrestrial atmospheric origin

2. Other Geophysical Phenomena: Scale and Physics Analysis

• Planetary-scale atmospheric waves: Rossby waves have wavelengths of 6,000–10,000 km (Holton & Hakim 2012), significantly longer than observed λ ≈ 3,330-4,549 km, and exhibit strong seasonal patterns inconsistent with observed temporal stability
• Ionospheric traveling disturbances: Large-scale TIDs propagate at 400–1000 km/h with wavelengths of 1,000–3,000 km (Hunsucker & Hargreaves 2003), but show strong diurnal and solar cycle dependencies systematically excluded through comprehensive TID analysis (Section 4.2.5)
• Magnetospheric current systems: Ring current and field-aligned currents create magnetic field variations at 2,000–5,000 km scales (Kivelson & Russell 1995), but these primarily couple to magnetic rather than timing systems, and lack the exponential spatial decay characteristic of screened fields
• Tropospheric delay correlations: Water vapor patterns show correlations up to 1,000–2,000 km (Bevis et al. 1994), insufficient to explain the 3,330-4,549 km scale, and are systematically removed by zenith delay modeling in standard GNSS processing

Systematic conclusion: Comprehensive analysis of major geophysical phenomena reveals fundamental incompatibilities in spatial scale, temporal signature, processing immunity, and coupling mechanisms. The systematic exclusion of atmospheric loading through temporal, processing, and correlation pattern analysis, combined with scale mismatches for other geophysical processes, supports the interpretation as novel field coupling transcending conventional Earth system processes.

Conformal vs. Disformal Coupling Evidence

The observed spectral pattern—broadband correlations with gravitational enhancement—suggests contributions from both conformal and disformal metric modifications:

• Conformal coupling A(φ): Evidenced by broadband response (R² > 0.85 from 10-200 µHz, CV of R² across bands = 2.9%) indicating universal coupling to all frequency modes without frequency selectivity
• Disformal coupling B(φ): Evidenced by enhanced correlation lengths at tidal frequencies (λ = 4,677 km vs 1,502 km post-tidal), indicating preferential response to gravitational gradients
• Combined metric structure: g̃_μν = A(φ)g_μν + B(φ)∇_μφ∇_νφ appears consistent with observed frequency-dependent spatial scales but frequency-independent correlation strength

Frequency-Scale Relationship

The systematic variation of correlation length with frequency provides insight into the physical mechanisms:

• Tidal frequencies (10-30 µHz): λ = 4,677 km — planetary-scale gravitational forcing (Moon/Sun)
• Post-tidal (30-100 µHz): λ = 1,502 km — transition to regional-scale coupling
• Intermediate (100-500 µHz): λ = 1,459 km — stabilized local-scale response
• Control (>1000 µHz): λ = 2,149 km — systematic instrumental effects

The factor of 3× decrease in correlation length from tidal to post-tidal frequencies, while correlation strength (R²) remains high, suggests the same coupling mechanism operates at different spatial scales depending on the frequency of gravitational forcing.

Systematic Effects and Signal Purity

Control band analysis reveals that systematic instrumental effects are non-negligible but clearly discriminated from physical correlations:

• Systematic contribution: R² = 0.618 in control bands indicates reduced model fit quality compared to TEP band
• Model fit difference: ΔR² ≈ 0.33 between TEP band (R² = 0.952) and control bands (R² = 0.618)
• Physical interpretation: Observed correlations represent genuine physical signal superimposed on systematic background, not pure signal with zero systematics
• Validation strength: Clear quantitative separation enables robust signal discrimination despite presence of instrumental effects

Key Insight: Signal Persistence Despite Processing

The robustness of TEP correlations to GNSS processing provides critical evidence for signal authenticity through three key observations:

Conclusion: The persistence of physically validated correlations despite processing designed to remove them strengthens the case for genuine field coupling.

Systematic Signal Suppression Mechanisms

Standard GNSS analysis center processing applies multiple corrections that would specifically attenuate TEP signals:

1. Common-Mode Removal and Network Constraints

• Network datum constraints: All centers apply network-wide reference frame constraints that force the sum of clock corrections to zero, explicitly removing any globally coherent signal components
• Common-mode filtering: Systematic removal of signals common across multiple stations, precisely the signature predicted by TEP theory
• Reference clock stabilization: Centers typically stabilize against ensemble averages, further suppressing spatially correlated variations
• Impact: Network constraints are designed to eliminate globally coherent variations

2. Sidereal and Environmental Corrections

• Sidereal filtering: Routine removal of 24-hour and harmonically related periodicities would eliminate TEP signals coupled to Earth's rotation
• Environmental modeling: Tropospheric and ionospheric corrections remove spatially correlated atmospheric effects, potentially including genuine field-atmosphere coupling
• Multipath mitigation: Advanced multipath modeling may inadvertently remove coherent field effects that manifest as apparent signal path variations
• Effect: Environmental corrections remove spatially correlated atmospheric variations

3. Outlier Detection and Quality Control

• Statistical outlier removal: Automated detection of "anomalous" correlations would systematically exclude the strongest TEP events
• Inter-station consistency checks: Quality control algorithms designed to ensure station independence would flag and remove genuine field correlations
• Temporal smoothing: Many centers apply temporal smoothing that would blur rapid field variations during astronomical events
• Result: Quality control processes are designed to flag and remove inter-station correlations

Evidence for Systematic Underestimation

Multiple lines of evidence suggest the measurements represent lower bounds on true field coupling strength:

Processing-Dependent Signal Strength

• Center-specific variations: The 36% variation in λ across centers (3,330-4,549 km) likely reflects different degrees of signal suppression rather than measurement noise
• ESA vs CODE comparison: ESA's shorter λ (3,330 km) and higher amplitude (A=0.250) compared to CODE's longer λ (4,549 km) and lower amplitude (A=0.114) suggests different processing approaches preserve different aspects of the TEP signal
• IGS intermediate values: IGS shows intermediate characteristics (λ=3,764 km, A=0.194, exponential model), consistent with processing-dependent signal recovery

Residual Signal Persistence

• Survival of strong correlations: The fact that R² = 0.920-0.970 correlations persist despite aggressive processing suggests the underlying field coupling may be substantially stronger than measured
• Elevation dependence: The systematic increase in λ with station elevation (Q1 to Q5: 1,785 to 4,549 km) indicates atmospheric screening effects that would be amplified in raw data
• Astronomical modulations: Detection of eclipse and supermoon effects despite processing suggests stronger signatures may exist in unprocessed data

Phase Information Preservation

• Why phase survives: While amplitude corrections are applied per-station, relative phase information between stations remains largely intact, explaining why the phase-based method succeeds where amplitude methods fail
• Incomplete phase suppression: The phase-alignment index approach captures residual phase relationships that survive processing, but this represents only a fraction of the original signal
• Raw data potential: Access to raw pseudorange and carrier phase measurements could provide additional insights into the field coupling mechanism

Future Research Directions

The robustness of the findings to GNSS processing suggests several promising avenues for future investigation:

Immediate Research Priorities

• Raw pseudorange analysis: Direct analysis of unprocessed GPS pseudorange measurements before any corrections or constraints
• Carrier phase coherence: High-precision analysis of raw carrier phase data to detect field-induced timing variations at the wavelength level
• Multi-constellation validation: Extension to GLONASS, Galileo, and BeiDou raw data to confirm field universality
• Real-time monitoring: Development of dedicated TEP detection networks bypassing standard GNSS processing chains

Expected Raw Data Advantages

• Signal characterization: Direct analysis of unprocessed measurements could provide additional insights
• Temporal resolution: Sub-second detection of rapid field variations during astronomical events
• Spatial precision: Millimeter-level coherence detection through carrier phase analysis
• Dynamic range: Full access to field variations across all timescales and amplitudes

Scientific Impact: Raw data access would provide critical validation of the signal strength hypothesis and could reveal substantially stronger coupling than currently measured, advancing understanding of potential modifications to spacetime structure and enabling more rigorous tests of fundamental physics.

Key Insight: Signal Survival Strengthens Authenticity

The persistence of strong correlations despite GNSS processing designed to remove spatially correlated signals provides evidence for robust underlying field coupling. Processing artifacts would likely be eliminated by these standard procedures; the signal persistence combined with multi-modal physical validation is consistent with genuine phenomena.

Implication: The robustness of detected signals to standard processing procedures, combined with their physical validation through Earth motion and astronomical correlations, supports the interpretation as genuine field coupling phenomena rather than processing artifacts.

Summary of Evidence

1. Signal authenticity supported (Section 4.1): Eleven independent validation criteria passed, including comprehensive multiple comparison correction analysis (137-168 of 231 statistical tests maintain significance after ultra-stringent Bonferroni α = 0.000216, FDR q = 0.05, and family-wise corrections), complete validation spanning data quality verification through astronomical period analysis, and systematic statistical robustness through bootstrap, LOSO, LODO, and cross-validation methods
2. Alternative explanations addressed (Section 4.2): Systematic analysis rules out classical tidal contamination (post-tidal 30-40 µHz shows R² = 0.946, enhancement ratio only 1.58× not >3×), tested ionospheric effects (TIDs, TEQ), processing artifacts, and examined conventional geophysical phenomena through scale, temporal, and spectral mismatches
3. Spectral coupling characterized (Section 4.3): Cross-center multiband validation demonstrates remarkable consistency across independent processing chains. Multi-band analysis reveals conformal (broadband, CV of R² across bands = 2.9%) and disformal (gravitational enhancement, λ = 4,627 km mean at tidal frequencies) coupling contributions, with systematic effects quantified through control bands showing reduced model fit quality (R² = 0.618, ΔR² ≈ 0.33 from TEP band)
4. Theoretical consistency observed (Section 4.4): Observed primary correlation lengths (λ = 3,330–4,549 km) fall within predicted range for screened scalar fields, with derived field mass m_φ ≈ 4.34–5.93 × 10^-14 eV (same conversion as §1.1) and potential implications for fundamental physics including dark matter connections and fifth force constraints. Bootstrap validation shows center-specific ranges (see Table in §3.1.1), corresponding to a conservative union mass range of m_φ ≈ 3.65–6.53 × 10^-14 eV.
5. Robustness to processing effects (Section 4.5): GNSS processing is designed to remove spatially correlated signals; the survival of strong correlations (R² = 0.920-0.970) through such filtering, combined with multi-modal physical validation, provides evidence for genuine field coupling rather than processing artifacts
6. Dynamical time validation (Section 4.6): Comprehensive diurnal analysis of 72.4 million hourly records reveals systematic temporal variations consistent with TEP predictions for variable proper time flow rates, with synchronized early morning peaks, seasonal modulation, and time rate variations of ~10⁻¹⁷-10⁻¹⁶ fractional amplitude. The patterns are consistent with the central TEP tenet that "when" is as dynamical as "where," though alternative slow environmental covariates remain plausible without raw-data confirmation

Open Questions: Next-Generation Alternative Testing

Building upon the systematic exclusion of primary alternatives, the following hypotheses represent the natural research frontier. Each offers specific testable predictions that could further strengthen the TEP interpretation or reveal alternative physics:

1. Systematic Processing Correlations

• Hypothesis: GNSS analysis centers use similar reference models, creating artificial correlations that survive current null tests
• Mechanism: Shared satellite orbit models, common ionospheric corrections, or similar tropospheric delay models could induce distance-structured correlations
• Required Test: Analysis of centers using fundamentally different processing approaches (e.g., PPP vs. network solutions)

2. Atmospheric Loading Effects

• Hypothesis: Large-scale atmospheric pressure variations create correlated timing effects through gravitational loading
• Mechanism: Atmospheric mass redistribution affects local gravity, influencing atomic clock rates through gravitational redshift
• Predicted Scale: Continental-scale correlations (1,000-5,000 km) matching observed λ values
• Required Test: Correlation with atmospheric reanalysis data (ECMWF, NCEP)

3. Solid Earth Tidal Correlations

• Hypothesis: Imperfect solid Earth tide corrections create residual correlations with distance-dependent structure
• Mechanism: Earth tide models may not perfectly capture local geological variations, leaving systematic residuals
• Required Test: Analysis with enhanced tidal corrections, geological substrate correlations

4. Electromagnetic Field Coupling

• Hypothesis: Large-scale electromagnetic fields (magnetospheric, atmospheric electrical) couple to atomic clocks
• Mechanism: AC Stark shifts from ambient electromagnetic fields could modulate transition frequencies
• Required Test: Correlation with magnetometer networks, atmospheric electrical measurements

5. Thermal Environment Correlations

• Hypothesis: Large-scale temperature patterns create correlated thermal effects on clock performance
• Mechanism: Continental-scale weather patterns could induce systematic thermal effects on station infrastructure
• Required Test: Correlation with meteorological reanalysis, seasonal decomposition

Research Strategy: Systematic testing using independent datasets (atmospheric reanalysis, geological surveys, electromagnetic monitoring, meteorological data) represents the natural next phase of investigation. Success in excluding these alternatives would further strengthen the TEP interpretation through comprehensive process of elimination, while positive findings would redirect theoretical development—either outcome advances fundamental understanding.

Critical Priorities and Path Forward