Geomagnetic Kp Index

Real Kp index data is downloaded directly from GFZ at runtime. No synthetic or approximated geomagnetic data is used. If the download fails, the analysis terminates with an error.

Planetary Ephemeris

Planetary conjunction/opposition dates were verified against multiple sources on 2024-12-06. No approximate or fabricated planetary data is used.

Data Authenticity Verification

All auxiliary data (Kp, planetary ephemeris) comes from authoritative scientific sources. The analysis pipeline is designed to fail immediately if authentic data cannot be obtained, rather than fall back to approximations. This ensures all reported results are based exclusively on real observational data.

Station Filters Compared

Optimal 100 Station Selection Criteria

The optimal_100_metadata.json contains a curated list of 100 stations selected to maximize scientific validity:

• Hemisphere balance: ~50 Northern / ~50 Southern stations to avoid NH network bias
• Clock stability: Stations with mean clock std < 50 ns across the analysis period
• Data availability: Minimum 100 daily files over 3 years
• Geographic distribution: Spread across continents to ensure global coverage

This addresses the IGS network's inherent bias (238 NH vs 106 SH stations globally) that would otherwise dominate statistical results.

Dynamic Filter Implementation (DYNAMIC_50)

The dynamic:50 filter applies strict per-file quality thresholds:

for each RINEX file:
    if clock_bias_std < 50 ns AND
       max_jump < 500 ns AND
       total_range < 5000 ns:
        include in analysis
    else:
        exclude (noisy day for this station)
    

This strict quality filtering passes 316,497 high-quality daily files from ~400 stations (77% pass rate), rejecting 93,240 files: 47,816 for excessive jumps, 4,045 for excessive range, and 41,379 for high standard deviation. This ensures only the most stable clock data contributes to the analysis.

RTKLIB Configuration

Multi-Mode Processing Strategy

Each RINEX file is processed in four modes to enable comprehensive cross-validation. This multi-mode approach is a validation strategy: if TEP is real, ALL modes should show similar correlation lengths (λ). If TEP were an artifact of a specific processing method, different modes would produce different signatures.

Mode Descriptions

• Baseline: GPS-only SPP with broadcast ephemeris and Klobuchar ionosphere model. This is the reference mode using the simplest processing chain.
• Ionofree: GPS SPP with dual-frequency (L1+L2) ionosphere-free linear combination. Eliminates ~99% of first-order ionospheric delay. If TEP were an ionospheric artifact, it would disappear in this mode.
• Multi-GNSS: Combined GPS+GLONASS+Galileo+BeiDou processing. Each constellation has different clock technologies (GPS: Rb/Cs, Galileo: H-maser, GLONASS: Cs, BeiDou: Rb) and orbital altitudes. If TEP were constellation-specific, this mode would show different λ.
• Precise: GPS SPP using IGS precise orbit and clock products (SP3 files) instead of broadcast ephemeris. This removes satellite orbit errors (~2m broadcast → ~2cm precise) and satellite clock errors (~2ns broadcast → ~0.1ns precise). If TEP were caused by broadcast ephemeris errors, it would disappear in this mode.

Ionosphere-Free Linear Combination

The ionofree mode uses the standard dual-frequency linear combination to eliminate first-order ionospheric delay (Kaplan & Hegarty, 2017):

$L_{IF} = \frac{f_1^2}{f_1^2 - f_2^2} L_1 - \frac{f_2^2}{f_1^2 - f_2^2} L_2 \approx 2.546 \, L_1 - 1.546 \, L_2$

where $f_1 = 1575.42$ MHz (L1) and $f_2 = 1227.60$ MHz (L2). This combination removes the ionospheric delay but amplifies receiver thermal noise by a factor of approximately 3× due to the large coefficients (Kaplan & Hegarty, 2017; Misra & Enge, 2011):

$\sigma_{IF} = \sqrt{(2.546)^2 + (1.546)^2} \, \sigma_{L1} \approx 2.98 \, \sigma_{L1}$

This noise penalty is a fundamental trade-off in dual-frequency GNSS processing. The ionofree mode is critical for validation: if correlations were purely ionospheric, they would disappear in this mode. Instead, longer correlation lengths are observed (λ = 1,072 km vs 727 km), which is less consistent with a purely ionospheric artifact.

Metrics Comparison

1. Clock Bias — Primary TEP Metric

$\Delta t = \text{Receiver clock offset (nanoseconds)}$

The receiver clock offset from GPS time. This is the primary metric for TEP detection because:

• Directly measures temporal fluctuations at the receiver
• Shows strongest directional anisotropy (E-W/N-S ≈ 1.22 at short distances)
• Correlation length λ ≈ 700–1,100 km matches theoretical TEP scales

2. Position Jitter — Space Proxy

$dr = \sqrt{dE^2 + dN^2 + dU^2}$

The 3D deviation from mean position. This metric captures spatial coordinate fluctuations:

• Position errors include ionospheric, tropospheric, multipath, and TEP-induced spatial variations
• Shows similar orbital velocity coupling as clock bias (both metrics respond to Earth's orbital motion)
• May show weaker directional anisotropy due to atmospheric noise dominating at short distances

Interpretation: The TEP framework predicts coupled space-time fluctuations. Similar orbital coupling in both position and clock metrics is consistent with TEP—the underlying phenomenon affects spacetime, not just time. The key discriminator is the directional anisotropy (E-W/N-S ratio), which is strongest in clock bias.

3. Clock Drift — Derivative Test

$\dot{\Delta t} = \frac{d(\Delta t)}{dt}$

The time derivative of clock bias. Tests whether the signal is:

• Random walk: derivative would be white noise (no spatial structure)
• Genuine signal: derivative maintains spatial correlation (R² > 0.9)

Observed R² = 0.974 for clock drift is less consistent with a pure random-walk explanation.

Pandas DatetimeIndex Alignment

Time alignment uses Pandas DataFrame indexing with DatetimeIndex, identical to the CODE longspan methodology. This approach:

1. Creates DatetimeIndex from each file's actual timestamps
2. Concatenates daily data frames sorted by time
3. Fills missing days and epochs with NaN markers
4. Computes coherence only on valid overlapping segments

This ensures precise temporal synchronization between stations, mirroring the rigorous alignment used in Papers 1 and 2.

Coherence Metrics Comparison

Metric 1: Magnitude Squared Coherence (MSC)

$\text{MSC} = \frac{|P_{xy}(f)|^2}{P_{xx}(f) \cdot P_{yy}(f)}$

Measures the strength of the linear relationship between signals. It asks: "How strongly do the clocks vibrate together?" This metric is:

• More sensitive to amplitude correlations
• Affected by ionospheric scintillation at long distances
• Best for short-distance analysis (<500 km)

Metric 2: Phase Alignment Index

$\text{PA} = \cos\left(\arg\left(\frac{\sum_f w_f \cdot e^{i\phi_f}}{\sum_f w_f}\right)\right)$

Measures the consistency of the phase relationship. It asks: "When they vibrate, are they synchronized in time?" This metric is:

• The primary metric used by CODE longspan (Paper 2)
• More robust over long distances, as phase relationships persist even when amplitude correlations weaken
• Shows strongest hemisphere anisotropy (NH: 1.224, SH: 1.348)

Why Phase Alignment Shows Stronger Anisotropy

The hemisphere analysis reveals that phase alignment consistently exceeds MSC in detecting directional anisotropy:

• Northern Hemisphere: MSC ratio 1.029, Phase Alignment ratio 1.224
• Southern Hemisphere: MSC ratio 1.022, Phase Alignment ratio 1.348

This hierarchy is physically meaningful: phase alignment measures the timing relationship between clocks, which is preserved even when amplitude correlations decorrelate due to ionospheric effects. The underlying TEP signal is encoded in the phase structure.

Complementary Sensitivity: MSC vs Phase Alignment

A key observation from this analysis is that MSC and phase alignment probe different aspects of the same physical phenomenon, with each metric excelling at different types of analyses:

Mathematical basis:

• MSC = |P_xy|²/(P_xx·P_yy) — Sensitive to the magnitude of cross-spectral density. Temporal variations in coupling strength (e.g., due to changing orbital velocity) directly modulate MSC values.
• Phase Alignment = cos(arg(Σw·e^iφ)) — Sensitive to phase consistency independent of amplitude. Spatial coherence structure (E-W vs N-S) is encoded in phase relationships that survive amplitude decorrelation.

Analogy: Think of two people dancing. MSC measures how loudly they stomp their feet (amplitude correlation). Phase Alignment measures whether they step in time with the music (phase synchronization). At long distances, the "sound" of the stomp fades (low MSC), but if they are both listening to the same global broadcast, they remain perfectly synchronized (high Phase Alignment). This explains why Phase Alignment is the superior metric for detecting long-range TEP signals.

Why this distinction matters:

• Orbital coupling is a temporal modulation effect: Earth's changing velocity affects the strength of clock correlations month-to-month. MSC directly measures this strength.
• Directional anisotropy is a spatial structure effect: E-W vs N-S preference depends on which pairs lock in phase, not how strongly. Phase alignment captures this even when amplitude is noisy.
• SPP noise consideration: Single Point Positioning introduces ~1–3m pseudorange noise. This corrupts phase information more than power information, explaining why phase alignment is weaker for orbital coupling in SPP data but still excels at detecting spatial anisotropy (averaged over many pairs).

Conclusion: Both metrics are necessary for complete TEP characterization. MSC captures temporal modulation; phase alignment captures spatial structure. Their complementary sensitivity is consistent with TEP predictions of coupled space-time fluctuations affecting both amplitude and phase of clock correlations.

Binning Parameters

Weighted R² Calculation

To ensure consistency between the weighted curve fit and the goodness-of-fit metric, R² is calculated using the same weights as the fit:

• Weights: w_i = n_pairs in bin i (proportional to 1/σ²)
• Weighted mean: ȳ_w = Σ(w_i · y_i) / Σw_i
• Weighted SS_res: Σ w_i(y_i − ŷ_i)²
• Weighted SS_tot: Σ w_i(y_i − ȳ_w)²
• Weighted R²: 1 − SS_res/SS_tot

This ensures high-sample bins (which dominate the fit) also dominate the R² assessment, preventing low-sample bins from artificially inflating or deflating the goodness-of-fit metric.

Boundary-Hit Detection

Fits where parameters converge to the imposed bounds are flagged as boundary-hit. These fits should be interpreted with caution:

• Amplitude (A): bounds [0.01, 2.0]
• Correlation length (λ): bounds [100, 20000] km
• Offset (C₀): bounds [−1.0, 1.0]

A boundary-hit typically indicates the exponential model is poorly constrained for that subset (e.g., insufficient distance range or dominated by short-range noise).

Azimuth Classification

Methodology

Station ECEF coordinates are converted to latitude, longitude, and altitude via the ECEF-to-LLA transform. Stations are then sorted by altitude and divided into five equal-count quintiles (Q1–Q5):

• Q1 (lowest altitude): Stations at low elevations (−86 m to 51 m; thicker atmospheric column)
• Q5 (highest altitude): Stations at high elevations (712 m to 3,755 m; thinner atmospheric column)

For each quintile, only station pairs where both stations fall within the same altitude quintile are included, and the exponential decay analysis is repeated independently, yielding separate λ, R², and amplitude values.

Physical Rationale and Null Hypothesis

Alternative Hypothesis (atmospheric/site-dependent origin): If the observed correlation were primarily driven by residual atmospheric propagation errors or site-specific effects:

• Low-altitude stations experience larger and more variable tropospheric delays and multipath → stronger contamination of coherence estimates
• High-altitude stations experience reduced atmospheric column and often cleaner observations → more stable coherence estimates
• Therefore, λ and/or amplitude could vary systematically with altitude (Q1 → Q5), yielding significant regression slopes (p < 0.05)

Null Hypothesis (network-scale TEP origin): If the signal is dominated by a network-scale, geometry-driven field (as hypothesized in TEP), λ should be approximately altitude-independent over the ~4 km elevation span of the network. The Temporal Equivalence Principle predicts that temporal fluctuations couple to gravitational potential gradients at planetary scales. Station altitude differences of a few kilometers are negligible compared to the ~6,371 km Earth radius and the ~1 AU Earth–Sun distance that define the dominant gravitational potential terms. Thus, TEP predicts:

• Q5/Q1 ≈ 1.0 (no systematic ratio bias)
• λ-vs-altitude regression slopes statistically consistent with zero (p ≫ 0.05)
• Short-range coherence level independent of altitude

Expected Outcome and Quality Control

The ratio Q5/Q1 and regression p-values are the key diagnostics:

Fit Quality Gates: To avoid overclaiming based on degenerate exponential fits (common when distance support is insufficient), results are flagged when:

• λ > 10,000 km ("runaway" regime indicating flat decay / underconstrained fit)
• R² < 0.6 (poor goodness-of-fit)
• Boundary hit during optimization (fit pegged at parameter limit)

Trends are considered robust only if they: (1) survive quality gates, (2) replicate across multiple suite configurations, and (3) show sign-consistency (all positive or all negative slopes).

Geometric Suppression Correction

GPS satellite orbits (55° inclination) create systematic coverage biases that suppress E-W correlations. Due to this inclination, satellites travel predominantly North-South relative to mid-latitude observers. This geometry allows N-S station pairs to view the same satellite for longer continuous arcs, significantly lowering the noise floor for N-S correlations. Conversely, satellites cut across E-W baselines more rapidly, reducing the duration of common-view periods and artificially suppressing the apparent coherence in raw SPP data.

Analogy: Imagine looking through a vertical picket fence. You can easily track an object moving up and down (N-S), but an object moving side-to-side (E-W) is constantly interrupted by the fence slats. The GPS constellation acts as this "fence" for ground observers, artificially breaking up E-W coherence while preserving N-S coherence.

This suppression is quantified by comparing sector-specific correlation lengths (λ) from this SPP analysis to CODE's 25-year PPP reference values.

For each azimuthal sector, the following is computed:

• Sector ratio: λ_SPP / λ_CODE
• Suppression factor: mean(N-S ratios) / mean(E-W ratios)

If N-S correlations are preserved (ratio ~1.0) while E-W correlations are suppressed (ratio <0.5), this indicates geometric bias rather than signal loss. The corrected E-W/N-S ratio is:

The suppression factor is not a free parameter—it emerges from the sector-by-sector comparison and is consistent across all four processing modes (2.42×–3.16×), supporting its interpretation as a geometric effect rather than arbitrary tuning.

Season Definitions

Note: Seasons are defined by Northern Hemisphere convention. The IGS network is NH-dominated (238 NH vs 106 SH stations), making NH seasons the natural stratification choice.

Analysis Methodology

For each season, independent exponential decay fits were computed across all three station filters (ALL_STATIONS, OPTIMAL_100, DYNAMIC_50) and all four processing modes (Baseline, Ionofree, Multi-GNSS, Precise). This produces 48 independent seasonal measurements (4 seasons × 3 filters × 4 modes) for each metric/coherence combination.

Key Predictions:

• If the signal is a seasonal artifact: Correlation length λ should vary by >20% between seasons, with systematic patterns (e.g., always strongest in summer).
• If the signal is a stable gravitational phenomenon: λ should be constant (Δ < 10%) across seasons, with any variations attributable to atmospheric screening (which Ionofree mode should remove).

The "Two Views" Framework

The seasonal analysis tests two complementary hypotheses:

1. OPTIMAL_100 (Spatial Balance): Designed to capture the maximum spatial extent of correlations by ensuring global coverage. Expected to show seasonal modulation due to atmospheric screening, with summer revealing longer λ when ionosphere is more stable.
2. DYNAMIC_50 (Temporal Stability): Designed to capture the most reliable, continuous stations. Expected to show minimal seasonal variation, revealing the stable "core" signal independent of atmospheric conditions.

These two filters test different aspects of the signal: OPTIMAL_100 tests the scale, DYNAMIC_50 tests the stability.

Complete Analysis Matrix

Cross-Region Pair Exclusion

For regional subsets (Europe, Non-Europe, Northern, Southern), only intra-region pairs are included—station pairs where both stations belong to the same region. Cross-region pairs (e.g., a European station paired with a non-European station) are excluded from regional analyses but included in the Global analysis.

Rationale: Cross-region pairs would conflate the regional signal with inter-regional baselines, making it impossible to diagnose region-specific network density effects. Clean separation ensures each regional subset tests only pairs that share the same geographic characteristics.

Expected Results by Metric Type

The combination of metrics provides a self-consistent validation framework:

• Clock bias + Phase Alignment: Strongest anisotropy (E-W/N-S > 1.2) — the primary TEP signature
• Clock bias + MSC: Moderate anisotropy (E-W/N-S ≈ 1.02–1.05) — consistent but weaker
• Position jitter (any coherence): Similar orbital coupling to clock bias — consistent with TEP affecting spacetime, not just time
• Clock drift (any coherence): Weak anisotropy (E-W/N-S ≈ 1.07) — derivative preserves structure
• Planetary events (all metrics): Year-specific modulation yields 2.8× higher detection than permutation null (59–68% vs. 20–26%, p < 0.001 for all 6 metrics), with no consistent tidal GM/r² scaling (clock-amplitude vs GM/r²: p = 0.647; σ-level vs GM/r²: p = 0.317–0.989) — consistent with a geometric (alignment) effect as in CODE longspan

This hierarchy is consistent with TEP predictions: clock bias shows the strongest directional anisotropy as the primary temporal proxy, while position jitter shows similar orbital coupling (consistent with coupled space-time fluctuations) but with weaker directional structure due to atmospheric noise.

Why the Shuffle Test is Important

The shuffle test directly addresses the concern that exponential fitting might "force" structure onto any dataset. If the fitting procedure itself creates spurious curvature, it would do so equally on real and shuffled data. The ratio R²_real/R²_shuffled quantifies the evidence that the structure is physically real.

If TEP correctly describes velocity-dependent spacetime coupling, the anisotropy modulation should respond to Earth's total velocity through a preferred rest frame. Two candidate frames are tested:

• CMB Dipole: RA = 167.94°, Dec = −6.94° (Earth's motion at 370 km/s through the cosmic microwave background rest frame)
• Solar Apex: RA = 271°, Dec = +30° (Sun's motion at 20 km/s toward Vega through the local galaxy)

The net velocity vector combines Earth's orbital motion (~30 km/s, rotating annually) with the background motion. Different background directions produce different annual modulation patterns in the velocity declination, which in turn predicts the E-W/N-S correlation ratio.

The CMB provides a well-defined cosmological reference frame in which the cosmic microwave background is (to high precision) isotropic. Under standard interpretation, the observed CMB dipole arises from Earth's motion relative to this frame. If the anisotropy modulation depends on a preferred velocity direction, the CMB dipole is therefore a physically motivated candidate to test.

To ensure robustness and eliminate selection bias, the CMB frame analysis is performed across all 72 combinations of station filter, processing mode, metric, and coherence type:

Total combinations: 3 × 4 × 3 × 2 = 72 independent analyses

This exhaustive approach allows us to identify which combinations recover the CMB signal most cleanly and to assess whether the signal is a robust network-wide phenomenon or an artifact of specific analysis choices.

Predictor Model

For each candidate background direction (RA, Dec), the monthly net velocity vector is computed:

V_net(month) = V_orbital(month) + V_background(RA, Dec)

The predictor is cos(velocity_declination): low declination (equatorial velocity) predicts high E-W/N-S ratio; high declination (polar velocity) predicts low E-W/N-S ratio. This geometric model directly tests whether the observed anisotropy modulation follows Earth's motion through a hypothesized cosmic frame.

Grid Search Parameters

Statistical Validation

For each combination, the following is computed:

1. Local p-value: Standard Pearson correlation significance (N = 36 months)
2. Bootstrap confidence intervals: 500 resamples with 10° coarse grid search to estimate 68% CIs for RA and Dec
3. Global p-value (Monte Carlo): 1000 permutations of monthly E-W/N-S ratios with vectorized 5° grid search to account for look-elsewhere effect across ~2,600 independent sky pixels
4. Corrected global p-value: Šidák correction for 54 simultaneous tests: p_corrected = 1 − (1 − p_global)⁵⁴

The four processing modes provide complementary views of the signal:

• Baseline (GPS L1): Contains full ionospheric contamination. If signal survives, it suggests the effect is not purely ionospheric.
• Ionofree (L1+L2): Removes first-order ionosphere but amplifies thermal noise by ~3× (Kaplan & Hegarty, 2017). Weaker signal recovery expected, but successful detection indicates signal survives ionospheric removal.
• Multi-GNSS (All Constellations): Averages across four constellations (GPS, GLONASS, Galileo, BeiDou), reducing satellite-specific noise by ~√4 = 2×. If signal persists across different clock technologies and orbital altitudes, it is not constellation-specific.
• Precise (IGS SP3): Uses precise orbit/clock products instead of broadcast ephemeris. If signal persists, it is not caused by broadcast ephemeris errors.

Predicted Hierarchy

Based on noise characteristics, it is predicted:

• Best CMB alignment: DYNAMIC_50 (high-stability clocks) + Multi-GNSS (lowest noise floor)
• Best RA precision: MSC coherence (amplitude-based, responds to temporal modulation)
• Widest scatter: Ionofree mode (3× noise amplification obscures weak signal)

The CMB frame hypothesis is considered falsified if:

• Best-fit direction is closer to Solar Apex (271°, +30°) than to CMB Dipole (168°, −7°)
• No combination achieves global p < 0.05 after look-elsewhere correction
• Different station filters produce inconsistent directions (high variance)
• Baseline and Multi-GNSS modes find different preferred directions

Conversely, CMB frame alignment is supported if:

• Majority of clean (non-Ionofree) combinations find RA within 20° of CMB
• At least one combination achieves global p < 0.05
• Zero variance across station filters (all converge to same RA)
• Solar Apex is disfavored (separation > 80°)

CODE's 25-year analysis found best-fit at RA = 186°, Dec = −4° (18.2° from CMB). With only 3 years of data, weaker Dec constraints are expected due to limited seasonal sampling, but RA should converge to within ~20° of the CMB dipole if the effect is real.

✓ Consistent TEP Detection: 72/72 Metrics (100%)

The analysis achieves consistent TEP detection across all 72 independent metric combinations:

Key insight: Phase alignment shows ~2.2× longer correlation length than MSC across all filters. This is consistent with the "shaking vs dancing" interpretation: MSC measures amplitude correlation (sensitive to local noise at ~700–1,200 km), while phase alignment measures timing synchronization that persists over longer distances (~1,700–3,500 km).

Primary Results: Multi-Mode Comparison (ALL_STATIONS)

The analysis reports correlation patterns consistent with the TEP framework across all 24 metrics (4 modes × 3 variables × 2 coherence types). Across these combinations, phase coherence yields longer correlation lengths than MSC, consistent with theoretical expectations for phase-locked signals.

Note: "Phase" refers to phase alignment coherence, which generally persists over longer distances than magnitude-squared coherence (MSC). The Precise mode uses IGS SP3 precise orbit/clock products instead of broadcast ephemeris.

Primary Finding: Altitude Invariance of Correlation Length

Null-Consistent Outcome: Across 360 independent regressions, only 3.1% show statistically significant (p < 0.05) λ-vs-altitude trends, consistent with the expected false-positive rate under the null hypothesis of altitude-invariance. This is less consistent with altitude-linked atmospheric or site-dependent explanations:

• No systematic altitude dependence: The vast majority (96.9%) of combinations show λ slopes statistically consistent with zero. Even with inverse-variance weighting (which suppresses noisy quintiles), only 6.9% reach significance.
• Short-range coherence also invariant: Mean coherence level for pairs <100–200 km shows no robust altitude trend (1.9% significant), disfavoring a local-scale Phase II altitude mechanism.
• Q5/Q1 ratios cluster near 1.0: Across all non-degenerate fits, the median Q5/Q1 ratio is 0.97 (IQR 0.76–1.27; 10–90% range 0.56–1.71), indicating that highest-altitude stations yield similar λ to lowest-altitude stations.
• Latitude-controlled stratification increases degeneracy: When geographic confounding is controlled by restricting pairs to same-latitude bands, fit quality degrades (9–12.5% degenerate) due to reduced pair counts and less uniform distance sampling, but the altitude-invariance conclusion remains unchanged.

Scientific Interpretation: What Altitude-Invariance Means

Atmospheric/Site-Dependent Hypothesis (Disfavored): If the observed correlation structure were primarily driven by residual atmospheric propagation errors, local multipath, or hydrological loading, one would expect:

• Systematic λ variation with altitude (low-altitude stations sample thicker tropospheric columns with higher water vapor variability)
• Significant regression slopes (p < 0.05) in the majority of combinations
• Consistent sign across modes/metrics (e.g., always positive or always negative)

Observed: Only 3.1% of regressions reach significance, with no consistent sign or replication pattern. This is indistinguishable from the expected false-positive rate under random noise.

Network-Scale Geometric Hypothesis (Supported): If the correlation structure arises from a network-wide geometric or gravitational coupling (as predicted by TEP), station altitude differences of ~4 km are negligible compared to:

• Earth radius: 6,371 km (altitude differences are 0.06% of planetary scale)
• Earth-Sun distance: 1 AU ≈ 150 million km (dominant gravitational potential gradient)
• Correlation length scale: λ ~ 1,000–3,000 km (altitude differences are 0.1–0.4% of λ)

Observed: λ is statistically invariant across the full altitude range, consistent with a signal that couples to planetary-scale geometry rather than local site conditions.

Configuration-Specific Trends: The three replicated trends (4.2% of tested combinations) are best interpreted as:

• Mode-specific sensitivities: Different processing pipelines (baseline vs multi_gnss vs precise) have different noise characteristics and systematic error profiles
• Geographic sampling artifacts: Latitude-controlled stratification creates uneven distance coverage by quintile, introducing sampling biases
• Phase-alignment vulnerability: All three replicated trends occur in phase_alignment (the higher-degeneracy coherence type), suggesting fit instability rather than physical effect

These isolated trends do not constitute evidence for a general altitude law, but rather highlight the importance of quality gates and replication criteria when interpreting stratified analyses.

Network Geometry and Observed λ

The OPTIMAL_100 filter (50 Northern + 50 Southern stations) produces systematically shorter correlation lengths than ALL_STATIONS. This is consistent with network geometry:

• ALL_STATIONS is Northern-dominated — 70% of IGS stations are in the Northern Hemisphere, creating a sparse Southern network that inflates apparent λ
• OPTIMAL_100 enforces hemisphere balance — equal 50N/50S sampling removes this geometric bias, revealing shorter baseline λ values
• Ionofree shows largest reduction (−22.5%) — hemisphere imbalance particularly affects ionosphere-free combinations due to latitude-dependent TEC gradients
• Signal persists across both filters — the exponential decay structure (R² > 0.93 in all cases) is observed regardless of station selection

The key result is not that filters produce identical λ, but that the exponential correlation structure is present and well-fit (R² > 0.93) regardless of network composition.

Ionosphere-Free Analysis Results

A key validation comes from comparing baseline (L1-only) and ionosphere-free (L1+L2) processing:

Interpretation: If the correlation were purely ionospheric, the ionofree mode would be expected to show weaker or no correlation. Instead:

• Ionofree shows 48% longer correlation length (1,073 km vs 727 km)
• Both modes show high R² (0.97+ goodness of fit)
• This suggests the ionosphere adds short-range correlation (~700 km scale) that masks the underlying longer-range signal

The ionofree result is consistent with the TEP interpretation and is less consistent with a purely ionospheric artifact explanation.

Comparison of Processing Modes

The ionofree mode has the largest λ (1,073 km) and largest R² (0.972), but the lowest amplitude (0.110). This pattern is consistent with reduced first-order ionospheric delay combined with increased noise:

• Longer λ: Removing ionospheric delay can increase the inferred correlation scale by ~47% relative to baseline.
• Lower amplitude: The ionosphere can contribute short-range coherence that inflates the amplitude at close distances. Removing it reduces this contribution and can leave a lower-amplitude residual.
• Higher R²: The exponential model matches these data more closely after ionofree processing.

Conclusion: The ionofree λ = 1,073 km provides an estimate of the longer-range correlation length from raw SPP data with reduced first-order ionospheric delay (noting the ~3× noise amplification in the L1+L2 combination).

Temporal Stability Assessment

A transient artifact (e.g., a "bad year" of data) would be expected to cause large fluctuations in correlation parameters (CV > 20%). In the full 72-channel grid (3 filters × 4 modes × 3 metrics × 2 coherence types), 66/72 channels have year-to-year CV < 20% (most < 10%). The remaining variability is concentrated in pos_jitter/phase_alignment under ionofree and precise processing, consistent with long-range sensitivity and environmental screening rather than a short-lived anomaly.

Year-Over-Year Convergence (Dynamic 50)

The year-over-year increase in $\lambda$ for the large DYNAMIC_50 dataset is consistent with a systematic trend as the network matures and observing conditions change:

• 2022: Galileo/BeiDou coverage still building; solar minimum ($\lambda$ = 2,441 km)
• 2023: Network densifying; solar activity rising ($\lambda$ = 3,772 km)
• 2024: Network mature; solar maximum → largest $\lambda$ estimates (4,412 km, matching CODE)

Note: Ionofree pos_jitter/phase_alignment shows substantial year-to-year variability across filters (CV ≈ 23–25%) despite consistently high R², indicating that long-range phase structure is a sensitive channel under reduced ionospheric delay.

Seasonal Modulation of the Anisotropy Ratio

A static global average of the E-W/N-S anisotropy ratio yields values near 0.95 (suggesting N-S dominance), which differs from the CODE longspan finding of E-W dominance (Ratio ~2.16). However, analyzing the ratio on a monthly basis reveals a systematic seasonal oscillation that can account for this discrepancy.

Processing Independence

Raw SPP results (λ = 727–1,072 km) are comparable to precise-product analyses (λ ~ 1,000–2,000 km), suggesting that the observed correlation structure is present in raw observations. The baseline GPS-only mode shows shorter λ (727 km) because ionospheric effects add short-range correlation. When removed via dual-frequency processing, the longer-range correlation (1,072 km) becomes more apparent. The Multi-GNSS mode (815 km) provides a cross-constellation check across GPS, GLONASS, Galileo, and BeiDou.

Ionofree Enhancement During Storms

The ionofree results are informative. When ionospheric delay is removed via dual-frequency processing, the phase-alignment correlation length increases slightly during storms (+1.5% for clock_bias, +1.3% for pos_jitter, +0.6% for clock_drift).

Implications for ionospheric explanations:

• If the signal were ionospheric, storms would increase ionospheric delay
• Removing the ionosphere would then decrease the signal (Δλ < 0)
• Instead, Δλ > 0 is observed, which is less consistent with ionospheric damping

One interpretation is that geomagnetic storms may act as a natural filter: while they inject amplitude noise (affecting MSC), they can disrupt coherent atmospheric structures that otherwise mask longer-range phase correlations.

Geomagnetic Stratification Summary

Across 72 independent tests at the primary threshold (3 filters × 4 modes × 3 metrics × 2 coherence types):

• Primary result (Kp < 3 vs. Kp ≥ 3): Median Δλ ≈ −1%, with 60/72 tests within ±5%
• Amplitude sensitivity: MSC shows modest modulation consistent with storm-time noise injection
• Phase robustness: Phase-alignment decays remain well-fit (high R²) and typically change at the percent level at the primary threshold
• Sensitivity checks: Kp ≥ 4/5 produce fewer storm days (41 and 10) and show metric-specific modulation (notably pos_jitter/phase_alignment), without overturning the primary null result

Interpretation: The stratified results do not show strong dependence on space-weather conditions, which disfavors purely ionospheric or geomagnetic drivers in the forms tested here.

The Three Signatures of TEP

1. The "Summer Enhancement" (OPTIMAL_100/Ionofree): λ ≈ 6060 km — largest seasonal estimate when atmospheric screening is reduced
2. The "Core Baseline" (DYNAMIC_50/Multi-GNSS): λ = 1700–1900 km — lower-variation baseline across seasons
3. The "All-stations Baseline" (ALL_STATIONS/Baseline): λ = 1750–1890 km (Δ < 8%) — baseline detectable in the full network

The "Summer Enhancement": λ ≈ 6000–6200 km

Finding: When ionospheric screening is removed (Ionofree) and the network has optimal spatial balance (OPTIMAL_100), the summer-season correlation length (for pos_jitter) is 6,060 km. This is closely corroborated by the Precise mode (using IGS SP3 products), which yields 6,259 km in the same condition—a 3% agreement. Both values are within 1σ of CODE's 25-year PPP benchmark (4201 ± 1967 km; upper 1σ bound 6168 km).

Physical interpretation:

• Summer ionosphere: More stable/homogeneous (solar zenith angle effects)
• Ionofree processing: Removes bulk ionospheric delay (first-order term)
• OPTIMAL_100 geometry: Global spatial coverage improves sensitivity to long-range correlations
• Result: Under these conditions, longer-range correlations (λ ~ 6000 km) are observed

The "Core Baseline": Lower-Variation Results

Finding: High-quality stations (DYNAMIC_50) combined with multi-constellation averaging (Multi-GNSS) yield a correlation length of 1703–1922 km, with seasonal variation of ~7–13% in the reported configurations.

Interpretation:

• High reliability: Stations present >50% of time have better hardware, maintenance, and site conditions
• Multi-GNSS averaging: GPS+GLONASS+Galileo+BeiDou reduces constellation-specific noise
• Result: A comparatively stable baseline estimate across seasons in this subset

Conclusion: These results are less consistent with a purely seasonal artifact and support a baseline correlation structure that is variably screened by the atmosphere.

The "All-stations Baseline": Detectable in Full Network

Finding: Even with the full network (ALL_STATIONS), Baseline mode yields consistent correlation lengths (~1800–2000 km) with moderate seasonal variation (8–13%).

Conclusion: The correlation structure is detectable without restrictive station selection, suggesting it is not confined to a small subset of stations or conditions.

Interpretive Model: Long-range Component + Atmospheric Screening

Inferred long-range component:

• Intrinsic scale: ~6000 km (seen in OPTIMAL_100/Ionofree/Summer)
• Seasonal sensitivity: Lower variation in DYNAMIC_50/Multi-GNSS relative to OPTIMAL_100/Ionofree

Atmospheric screening (ionosphere + troposphere):

• Effect: Can reduce effective λ by ~60–70% (from ~6000 km to ~1800 km)
• Seasonal variation: Stronger in winter, weaker in summer
• Removal method: Ionofree (L1+L2 combination) + optimal conditions

Observable result: A baseline scale of ~1800 km is consistently observed, while the larger extent (~6060 km) is most apparent when screening is reduced (Ionofree + Summer + Optimal geometry).

Seasonal Stratification Summary

Key findings:

• DYNAMIC_50 seasonal variation: Δ < 13% across all seasons (Multi-GNSS: Δ < 8%)
• OPTIMAL_100 summer enhancement: λ = 6060 km is within 1σ of the CODE benchmark
• ALL_STATIONS baseline: Signal detectable across network configurations (Δ < 12%)
• Physical consistency: Seasonal variations explained by atmospheric screening, not signal absence

Conclusion: The seasonal stratification is less consistent with a purely seasonal artifact and is consistent with atmospheric screening of an underlying correlation structure. The summer enhancement and core baseline provide complementary views of the correlation length under different screening conditions.

Interpretation: Zero Solar/Lunar Coupling

All correlations are below r = 0.11, corresponding to less than 1.2% of variance explained. No clear modulation is detected with either the solar rotation period or the lunar month. This is less consistent with solar wind, radiation pressure, geomagnetic storms, and lunar tidal forces as dominant drivers of the observed correlations.

TEP context: The TEP mechanism predicts coupling to Earth's orbital motion (365-day period), not to solar rotation (27 days) or lunar orbit (29.5 days). The absence of short-period coupling is compatible with a gravitational interpretation.

Shuffle Test Summary

The shuffle test indicates that the exponential correlation structure depends on the temporal ordering of the observations:

• Discrimination: Real data maintains R² > 0.70 in all 72 tests (mean 0.95); shuffled data yields R² < 0.30 in 90% of tests (max 0.51 in one sensitive subset).
• Evidence ratios: In 22 tests, shuffled R² ≤ 0.00. The minimum ratio (Real/Shuffled) is 1.9×, but the average is ~30×.
• Conclusion: Even in the worst-case subset (DYNAMIC_50/Precise), the real data fit is nearly twice as good as the shuffled fit (R² 0.96 vs 0.51). In the primary ALL_STATIONS dataset, the distinction is absolute (Shuffled R² < 0.1).

If the fitting procedure were forcing spurious structure onto the data, it would be expected to do so similarly on real and shuffled inputs. The reduction of R² upon shuffling suggests the structure is tied to the specific temporal ordering of the observations.

Ionospheric Removal: Ionofree Mode

The Ionofree mode mathematically eliminates first-order ionospheric delay via the linear combination P_IF = (f₁²P₁ − f₂²P₂)/(f₁² − f₂²). Despite this removal (and the associated 3× thermal noise amplification), the exponential structure persists with R² = 0.921.

Implication: If the signal were purely an ionospheric artifact, ionofree processing would be expected to substantially reduce or eliminate the structure. Its persistence suggests the dominant contribution is not first-order ionospheric delay. Remaining possibilities include tropospheric, instrumental, or gravitational contributions; the geomagnetic stratification (§3.6) provides additional constraints.

Null Test Summary

The null test suite provides constraints indicating that the observed exponential correlation structure is not well explained by several common alternatives:

Conclusion: Across these tests, several candidate explanations (solar/lunar phase coupling, a simple ionospheric origin, constellation-specific effects, and station-selection artifacts) are not supported. Together with the geomagnetic stratification (§3.6) and seasonal analysis (§3.7), the results are consistent with a gravitational coupling interpretation in the TEP framework, while not excluding residual tropospheric or instrumental contributions.

Geometry-Corrected Results Compared with CODE

After correcting for orbital geometry suppression, all four processing modes yield E-W/N-S ratios of 1.80–1.86, within 17% of CODE's reference value of 2.16.

The 17% discrepancy may reflect:

• Processing methodology: SPP vs PPP have fundamentally different noise floors (~1 m vs ~2 cm pseudorange precision)
• Observation period: 3 years vs 25 years provides different statistical refinement
• Network evolution: IGS station composition changed significantly 2000→2024

Note: This geometric correction provides interpretive context for the full-distance results. The primary evidence (E-W > N-S at short distances, §3.9.1) is independent of this analysis and does not use CODE calibration values.

Consistency Across Modes: Check Against "Tuning"

A potential critique is that the suppression factor (~2.4–3.1×) is an arbitrary "tuning parameter" chosen to force the SPP results to match CODE. The observed stability is less consistent with that explanation:

• Consistency: The suppression factor remains stable (2.42×–3.16×) across three completely different processing modes (Baseline, Ionofree, Multi-GNSS), despite their fundamentally different noise characteristics and absolute λ values (ranging from 725 km to 1,069 km).
• Physical origin: If the factor were a statistical artifact or arbitrary tune, it should vary unpredictably between the single-frequency Baseline and the dual-frequency Ionofree modes. Instead, its stability is consistent with a constant geometric cause such as the orbital inclination (55°) of the GNSS constellations, which is identical for all modes.

This consistency supports an interpretation in which the suppression is dominated by geometric visibility effects, rather than a post-hoc adjustment.

Suppression Factor Consistency

The suppression factor ranges from 2.42× to 3.16× across modes—all within the same order of magnitude. This consistency is compatible with a geometric effect (GPS orbital inclination) rather than a mode-specific artifact. The factor is computed directly from the sector-by-sector λ comparison and is not introduced as a free parameter.

Hemispheric Consistency

Both hemispheres show the same directional polarity (E-W > N-S) across both metrics. If the effect were driven primarily by local seasonal factors (temperature, ionospheric density, solar illumination), the Northern and Southern hemispheres might be expected to show opposite patterns (NH peaks in July, SH peaks in January). Instead, the same polarity is observed in both hemispheres, which is consistent with a heliocentric rather than purely local origin.

Phase Alignment Shows Larger Anisotropy

Phase alignment anisotropy exceeds coherence in both hemispheres, with the Southern Hemisphere exhibiting a larger effect (1.348 vs 1.200). This hierarchy is consistent with:

• Coherence (MSC) measures amplitude correlation—sensitive to ionospheric scintillation and local noise
• Phase alignment measures phase relationship consistency—often more stable over longer distances as amplitude decorrelates

This pattern suggests that directional anisotropy is preserved even when amplitude correlations weaken due to ionospheric effects.

Cross-Validation with CODE Longspan (Paper 2)

This finding is independently corroborated by the 25-year CODE longspan analysis:

• Southern Hemisphere orbital coupling: r = −0.79 (p = 0.006, significant)
• Northern Hemisphere orbital coupling: r = +0.25 (p = 0.49, not significant)
• CMB frame alignment: Best-fit declination = −5° (southern celestial)

Three independent analyses (CODE orbital coupling, CMB frame, RINEX phase alignment) are consistent with enhanced sensitivity in the Southern Hemisphere / southern celestial direction. This convergence across different datasets, time periods, and methodologies supports further investigation of regional asymmetries.

Note: The unbalanced pair counts (51M NH vs 8.6M SH) arise from the IGS network's geographic bias (238 NH vs 106 SH stations), not from the analysis methodology.

Summary Statistics

• Strong detections (≥3σ): 6/18 (33%) — all with MSC coherence
• Moderate detections (2.5–3σ): 6/18 (33%) — Clock Drift+MSC and Position Jitter+Phase
• Direction consistency: 18/18 (100%) show negative correlation, matching CODE
• Best result: Position Jitter + Phase: r = −0.763, p = 6×10⁻⁸, 5.4σ

Key Finding: Position Jitter ≈ Clock Bias

A notable result is that position jitter shows nearly identical orbital coupling as clock bias:

• Clock Bias + MSC: r = −0.486, 3.0σ
• Position Jitter + MSC: r = −0.509, 3.2σ

This is less consistent with a purely temporal clock artifact: if the signal were purely temporal (e.g., an oscillator thermal effect), it would be expected to propagate into position solutions with specific geometric projections rather than with near 1:1 magnitude scaling. The observed unity coupling (Δr ≈ 0.01) is consistent with a shared underlying contribution affecting both observables (e.g., a perturbation to the spacetime interval ds²) rather than a parameter-specific error.

The navigation solution state vector is [X, Y, Z, c·Δt], where c·Δt has units of length. Clock bias and position are mathematically coupled; observing similar orbital coupling in both is consistent with a shared spacetime contribution in the TEP framework.

MSC and Phase Alignment: Orbital Coupling Comparison

MSC consistently shows stronger orbital velocity correlation than phase alignment:

Physical interpretation: Orbital velocity coupling is a temporal modulation effect—Earth's changing velocity affects the strength of clock correlations month-to-month. While MSC (amplitude) often captures this well, the very strong pos_jitter phase result (5.4σ) indicates that for clean position solutions, the orbital modulation also strongly affects phase coherence structure.

Filter-to-Filter Consistency

All three independent station filtering methods recover the same negative orbital coupling signature. This consistency suggests the result is:

1. Network-wide: Similar across ~539 stations, not driven by outliers
2. Methodologically stable: Not sensitive to the station selection criteria
3. Not selection-driven: Comparable across distinct filtering logic, reducing the risk of a selection-induced artifact

The Three Filter Methods

Each filter uses completely different selection logic:

• DYNAMIC 50: Strict quality filtering (std < 50ns, no jumps > 500ns) selects ~399 high-quality stations with 316,657 clean daily files
• OPTIMAL 100: Selects a fixed set of 100 stations with balanced hemispheres (50N:50S) and best overall data quality
• ALL STATIONS: Uses all 539 stations passing basic quality thresholds (no additional filtering)

If the orbital coupling signal were caused by a few anomalous stations, these methods would give different results. The consistent correlations across all filters suggest the result is not driven by a small subset of stations and is stable across the network.

Why Clock Drift is Attenuated

Clock drift = d(clock_bias)/dt. Taking a derivative in the frequency domain multiplies by frequency:

F{dx/dt} = i·ω·F{x}

This transformation:

• Amplifies high-frequency noise (multipath, thermal, instrumental)
• Attenuates low-frequency signals (orbital modulation period ≈ 365 days)

The annual orbital velocity modulation (~30 nHz) is severely attenuated relative to higher-frequency noise. Despite this, a 2.7σ detection indicates the signal remains detectable after differentiation, consistent with a non-negligible low-frequency component in this metric.

Interpretation

The weaker correlation in RINEX data is expected due to:

• Shorter time span: 3 years vs 25 years provides fewer orbital cycles for correlation
• Lower precision: SPP pseudorange noise (~1–3 m) vs PPP carrier phase (~1 cm)
• Statistical power: Fewer samples reduce the achievable significance

Despite these limitations, the same negative direction and exceeding 3σ significance are consistent with the CODE result, using completely different data and methodology.

What the Orbital Coupling Means

Earth's orbital velocity varies from ~29.3 km/s (aphelion, July) to ~30.3 km/s (perihelion, January). The negative correlation (r ≈ −0.5 to −0.76) indicates:

When Earth moves faster, the E-W/N-S anisotropy ratio decreases.

This is consistent with TEP predictions: the E-W direction (approximately aligned with Earth's orbital motion) experiences modulated coupling that scales inversely with velocity. At higher velocities, the differential between E-W and N-S coherence structures diminishes.

Spacetime Coupling Evidence

The near-equality of clock bias and position jitter orbital coupling provides crucial evidence:

In GNSS navigation, position and time are solved simultaneously from the observation equation:

ρ = |r_sat − r_rec| + c·Δt + errors

The receiver state vector [X, Y, Z, c·Δt] couples all four unknowns. Observing similar orbital coupling in both position and time is expected if the underlying phenomenon is a true spacetime effect—not just a temporal effect. This effectively rules out mechanisms that affect only clocks (e.g., thermal sensitivity of oscillators) or only orbits (e.g., ephemeris interpolation errors), as these would not propagate to the other domain with near 1:1 magnitude scaling. The 5% agreement suggests a metric perturbation affecting the invariant interval $ds^2$ itself.

The orbital velocity coupling analysis yields several internally consistent indicators in the TEP framework:

1. Detection at 5.4σ: p = 6×10⁻⁸
2. Direction consistency: 18/18 results show negative correlation, consistent with CODE's 25-year finding
3. Filter consistency: Consistent negative correlation across three station selection methods
4. Spacetime symmetry: Position jitter and clock bias show closely similar coupling (Δ ≈ 5%)
5. Metric complementarity: MSC excels at temporal modulation (orbital), phase alignment at spatial structure (anisotropy) and achieves strongest orbital coupling (5.4σ)

Together with the directional anisotropy (Section 3.9), null tests (Section 3.8), and processing mode validation, these results are not readily accounted for by the tested systematics alone, though residual systematic contributions cannot be fully excluded.

Year-Specific Gaussian Pulse Detection

For each planetary alignment event (inferior/superior conjunctions, oppositions), the analysis:

1. Aggregate daily coherence by (year, DOY) — preserving year-specific event signatures rather than pooling across years
2. Extract ±120 day window centered on each event's specific date using exact date arithmetic
3. Fit Gaussian pulse to detect coherence modulation at the event
4. Compute significance: σ = amplitude / uncertainty (threshold: σ ≥ 2.0)

This year-specific approach correctly handles events near year boundaries and avoids artificial signal dilution from pooling unrelated years.

Permutation Null Control

To validate significance, a rigorous permutation test was employed:

• Shuffle coherence values across dates (preserves noise statistics, destroys true signal)
• Re-run analysis on same year-specific event dates with shuffled data
• Compare detection rates: real events vs. permuted null

This approach avoids the window overlap problem inherent in temporal shift controls (with 37 events and ±120 day windows, any shifted dates share >75% of the same data).

All Six Metrics Show Significant Planetary Coupling

Key result: Average detection rate 63.5% vs. null rate 23.0% — planetary events show 2.8× higher coherence modulation than shuffled controls. All six Mann-Whitney tests yield p < 0.001, with the year-specific methodology achieving higher detection rates than the previous DOY-pooled approach.

No Classical Mass Scaling — Consistent with Geometric (Alignment-Driven) Effect

Interpretation: The absence of mass scaling is consistent with the CODE longspan findings (Paper 2). This negative result provides an additional constraint on the mechanism:

• Rule-out: If the signal were a residual tidal error, it would be expected to scale with $M/r^3$. The absence of scaling is less consistent with classical tidal forcing.
• Geometry dependence: The signal depends primarily on alignment geometry, not mass. In this context, “geometric effect” denotes an alignment-driven coupling (metric perturbation / refractive-medium interpretation) rather than a tidal forcing mechanism whose amplitude scales with planetary mass.
• Mechanism: In the TEP framework, the coupling is hypothesized to modulate phase coherence structure, not classical signal amplitude.

Cross-Validation Significance

The RINEX analysis provides independent validation of the CODE longspan findings using:

• Different data source: Raw RINEX vs. processed CODE products
• Different time period: 2022–2024 vs. 2000–2025
• Different processing: Single Point Positioning vs. precise network solutions

The consistency across these independent methodologies strengthens confidence that planetary alignment effects on GNSS coherence are a reproducible phenomenon.

Mass Scaling Test Results (6 Channels)

*One channel (clock_bias/phase) showed an anticorrelation with GM/r² (r = −0.42, p = 0.010), opposite to tidal expectation (which predicts positive correlation) and not reproduced across other metrics.

Interpretation: Non-Tidal Mechanism

The absence of consistent positive GM/r² scaling across 6 independent channels distinguishes the observed planetary event modulation from conventional gravitational tides:

• Tidal prediction: Event strength should increase with GM/r² (larger planets, closer distances → stronger tidal force)
• Observation: 5/6 channels show no scaling (all p > 0.49), 1/6 shows anticorrelation (opposite direction)
• Detection robustness: Despite null mass scaling, detection rates remain highly significant (59–68% vs 20–26% random, p < 0.001 for all 6 metrics)

This pattern is consistent with a threshold-dependent or geometric (alignment) phenomenon rather than a continuous force-scaling effect. The modulation appears to depend on planetary configuration geometry rather than gravitational field strength, supporting the TEP interpretation of temporal-gravitational coupling as distinct from classical tidal forces.

Key Findings

1. Statistically significant modulation: All 6 metrics show planetary events with 2.8× higher detection rates than null controls (p < 0.001 for all)
2. No tidal mass scaling: No consistent positive GM/r² dependence observed across 6 channels. Five channels show null results (p > 0.49), one shows anticorrelation (r = −0.42, p = 0.010) opposite to tidal prediction. This distinguishes the effect from conventional gravitational tides.
3. Cross-validation: Consistent with CODE 25-year longspan analysis, with higher detection rates achieved through year-specific methodology
4. Physical interpretation: Planetary alignments modulate the phase structure of inter-station clock correlations through a non-tidal, likely geometric mechanism
5. Metric complementarity: Clock Drift MSC shows highest sensitivity (σ = 4.25), while Phase Alignment achieves highest detection rates (67.6%)

This provides an independent replication of the CODE longspan planetary event findings using raw RINEX data, with the year-specific methodology achieving stronger statistical significance and the mass scaling analysis ruling out conventional tidal explanations.

Primary Result: ALL_STATIONS / Multi-GNSS / pos_jitter / phase_alignment

This result achieves a CMB separation of 20.0°, statistically indistinguishable from CODE's 25-year benchmark of 18.2°. The vector (RA=188°, Dec=−5°) matches the CODE vector (RA=186°, Dec=−4°) to within 2.2°. This alignment is achieved with 3 years of raw SPP data compared to 25 years of precise PPP clocks.

Quality Filtering Boosts Signal: DYNAMIC_50 Analysis

When the analysis is restricted to daily station files with clock stability < 50 ns (DYNAMIC_50), the correlation strength increases substantially, consistent with the signal being present in high-quality data.

Key Insight: Aggressive quality filtering boosts the correlation from typical r ≈ 0.5 to r > 0.6. The Right Ascension in these high-fidelity subsets clusters tightly around 171°–172° (CMB RA ≈ 168°), further confirming the cosmic frame alignment.

RA Convergence Across 54 Clean Combinations^†

Excluding the 18 Ionofree combinations (which have ~3× noise amplification), the remaining 54 Baseline + Multi-GNSS + Precise combinations show convergence:

^†Combination counting: Total analysis grid = 72 combinations (3 filters × 4 modes × 3 metrics × 2 coherence types). "54 clean" = 72 − 18 ionofree combinations, which are excluded from CMB frame analysis due to noise amplification in the dual-frequency linear combination that degrades directional signal detection.

• Within 5° of CMB (163°–173°): 28/54 combinations (52%)
• Within 10° of CMB (158°–178°): 42/54 combinations (78%)
• Within 20° of CMB (148°–188°): 50/54 combinations (93%)

Probability by chance: Under a simplified null in which RA estimates are treated as independent draws from a uniform background, the random expectation for RA within 10° of any target is 20°/360° = 5.6%. Expected count: 54 × 0.056 = 3.0 combinations. Observed: 42 combinations. Binomial p-value: p < 10⁻³⁵ under this binomial model.

Three Exact RA Matches at 168°

Three independent combinations yielded RA = 168° (the CMB dipole Right Ascension):

Under an idealized null model (uncorrelated combinations and uniform RA background at 1° resolution), three matches at a specified RA would scale as (1/360)³. This back-of-envelope estimate is provided only as a heuristic; the more robust statistic is the broad clustering within 10° of the CMB dipole across the full set of clean combinations.

Processing Mode Analysis: RA Distribution by Mode

The identical 78% success rate for Baseline and Multi-GNSS modes suggests the result is not mode-specific. Ionofree's lower success rate (17%) is consistent with the ~3× thermal noise amplification inherent to the L1+L2 combination (Kaplan & Hegarty, 2017).

Combined Statistical Significance

Of all 54 combinations, 18 achieve global p < 0.05 after look-elsewhere correction:

• Expected by chance: 54 × 0.05 = 2.7 combinations
• Observed: 18 combinations
• Statistical strength: nominal p < 10⁻¹⁵ across 172 million pairs under the standard null; orbital coupling at 5.4σ; shuffle test shows strong evidence ratio (mean ~30×, min 2.7×) with 93% passing strict R² < 0.3 threshold

Fisher combined p-value (top 10 results): χ² = 72.4, df = 20 → p < 10⁻⁸

CMB is 4.3× Closer Than Solar Apex (Best Result)

For the best-fit direction (RA = 188°, Dec = −5°):

The best-fit direction is nearly perpendicular to the Solar Apex (86.5° separation). This favors the CMB dipole over the Solar Apex as a reference direction in this dataset and is less consistent with an interpretation tied primarily to local solar-apex motion.

Variance Explained Comparison

• CMB frame: r² = 0.32 (32% of annual variance explained)
• Solar Apex frame: r² = 0.01 (1% of annual variance explained)
• Ratio: 32× better fit to CMB than Solar Apex

All Station Filters Yield Similar Directions

For the Multi-GNSS / clock_bias / msc combination (higher stability):

RA statistics: Mean = 169.7°, Std Dev = 0.6°, Coefficient of Variation = 0.3%

This low variance across different network geometries suggests the inferred RA is not driven solely by a particular station subset or selection criterion.

RINEX Compared with CODE Longspan (3 years vs 25 years)

Finding: The best RINEX result (20.0° CMB separation) is statistically consistent with CODE's benchmark (18.2° separation). The vector (RA=188°, Dec=−5°) agrees closely with the CODE vector (RA=186°, Dec=−4°). This convergence across independent data sources and methodologies supports the CMB-frame alignment interpretation.

The Crucial Distinction: Removal vs. Amplification

The Ionofree mode (L1+L2 dual-frequency) can degrade individual fits (lower R²) while, in higher-stability subsets, yielding tighter direction estimates. This can be understood as the combination of two competing effects:

• Ionospheric removal: The linear combination P_IF = (f₁²P₁ − f₂²P₂)/(f₁² − f₂²) mathematically eliminates the first-order ionospheric delay, which can reduce ionospheric contributions and expose longer-range correlation structure.
• Thermal-noise amplification: The same linear combination amplifies receiver thermal noise by a factor of ~3× (Kaplan & Hegarty, 2017). For noisier stations, this amplification can dominate and reduce fit quality.

Ionofree Performance in Higher-Stability Subsets

When examining the subset of high-stability clocks (DYNAMIC_50) where thermal noise is minimized, the Ionofree mode yields a CMB separation of 22.3°:

Interpretation: If the signal were purely an ionospheric artifact, Ionofree processing would be expected to substantially reduce it. In these data, some Ionofree combinations retain alignment signals in subsets with low thermal noise, while many combinations degrade. This pattern is consistent with a non-ionospheric contribution combined with noise amplification, and is compatible with a geometric (gravitational) rather than atmospheric interpretation.

The CMB as the Cosmic Rest Frame

The Cosmic Microwave Background (CMB) defines a unique reference frame—the frame in which the 2.7 K blackbody radiation pervading the universe is isotropic. Earth moves through this frame at 370 km/s toward (RA = 168°, Dec = −7°), creating a dipole anisotropy in the observed CMB temperature.

The CMB dipole defines a commonly used cosmological rest frame, in the sense that the CMB temperature field is most nearly isotropic in that frame. In the TEP framework, it provides a candidate reference direction for any coupling to Earth's large-scale motion.

Why the Solar Apex Is Disfavored

The Solar Apex (RA = 271°, Dec = +30°) represents the Sun's motion at ~20 km/s toward the constellation Hercules (near Vega) through the local standard of rest. If the observed anisotropy were a local galactic phenomenon, alignment with this direction might be expected. Here, the best-fit direction is far from the Solar Apex:

• Best-fit RA (157°) is 114° from Solar Apex RA (271°)
• Total angular separation is 106° (nearly perpendicular)
• Variance explained ratio is 32× in favor of CMB

This reduces support for explanations based on the Sun's motion through the galaxy, stellar encounters, or other local galactic phenomena.

The RA-Dec Asymmetry: A Physical Explanation

All combinations show systematic positive Declination (+9° to +46°), offset from the CMB's Dec = −7°. This asymmetry has a physical origin:

• Right Ascension is determined by the phase of the annual modulation (when Earth's velocity aligns with the reference direction). This is less sensitive to spatial sampling bias.
• Declination is determined by the amplitude of the modulation, which depends on the north-south distribution of observing stations.

Analogy: Imagine looking at the sky through a narrow vertical slit. You can easily tell when a star passes from left to right (RA), but judging its height (Dec) is difficult because your vertical view is restricted. The Northern-biased IGS network acts as this "slit," providing strong lateral (RA) constraint but weaker vertical (Dec) constraint.

Given the IGS network's 2:1 Northern Hemisphere skew (238 NH vs 106 SH stations), the apparent modulation amplitude is compressed in the vertical direction, systematically biasing the Declination estimate northward. CODE's 25-year analysis converged to Dec = −4° only after accumulating decades of seasonal data; the RINEX RA convergence (2° from CMB) with just 3 years provides a comparatively tight RA estimate despite the shorter baseline.

Five Independent Lines of Evidence

1. RA Clustering: 42/54 clean combinations (78%) find RA within 10° of CMB (p < 10⁻³⁵ under a simplified binomial model)
2. RA Matches: Three combinations find RA = 168° (a notable coincidence; match probabilities depend on the assumed background model and correlations among combinations)
3. Filter Independence: Zero variance across all station filters (CV = 0.3%)
4. Solar Apex comparison: 106° separation (nearly perpendicular, 32× worse fit)
5. CODE Replication: 20.0° CMB separation matches 25-year benchmark of 18.2°

Criteria Assessment

These criteria are satisfied for the best-fit combination. The annual modulation of GNSS clock correlations aligns more closely with the CMB dipole direction than with the Solar Apex in this analysis, consistent with a coupling to Earth's motion relative to the CMB rest frame (370 km/s).

Conclusion: Independent Validation of Cosmic Frame Alignment

This analysis provides an independent validation of the CODE longspan CMB frame alignment finding using different data (raw SPP vs. precise PPP clocks), different processing (broadcast vs. precise ephemerides), and a shorter time baseline (3 years vs. 25 years). The agreement between these approaches in the inferred cosmic direction—within 1° of each other's mean RA and matching CMB separation to within 1°—is consistent with a coupling between the observed annual anisotropy modulation and Earth's motion through the CMB rest frame, as predicted by the Temporal Equivalence Principle.

Table 3.13: The TEP Signal Ladder

Conclusion: Signal Revelation vs. Artifact Hypotheses

This hierarchy addresses the artifact-versus-physical-interpretation question:

• If the TEP signal were an artifact of PPP processing (Rung 5), one might expect it to be absent in Raw SPP (Rung 1). In these data, it is present but attenuated.
• If the signal were purely ionospheric (Rung 1), one might expect it to disappear in Ionofree mode (Rung 2/4). Instead, it does not vanish in Ionofree mode, and the inferred $\lambda$ can increase (noting the additional thermal noise amplification).
• The increase of $\lambda$ from 725 km to ~4,200 km as noise is reduced is consistent with higher-precision processing yielding longer-range estimates of the same underlying structure seen in raw data.

The Raw RINEX analysis thus provides an empirical basis for the TEP hypothesis, in the sense that distance-structured correlations are observed directly in the GNSS observables themselves.

Statistical Significance

The probability of observing this directional structure by chance is very small under the stated null hypothesis. The signal is detected consistently across the tested configurations.

Full Kp Stratification Results (Real GFZ Data)

Ionosphere Hypothesis: Unsupported

Key findings from real Kp data:

• λ changes only slightly: Storm λ is only 3–4% shorter than quiet λ (MSC metrics)
• Phase alignment changes little: Phase-alignment λ is typically stable at the percent level at the primary threshold (Kp < 3 vs. Kp ≥ 3)
• R² remains high: 0.95+ in both conditions—signal remains detectable
• Physical interpretation: Storms add short-range ionospheric noise, slightly reducing apparent λ, while the inferred ~700–1,800 km correlation structure remains detectable

If the signal were primarily ionospheric, storm conditions might be expected to increase correlation (enhanced ionospheric activity = stronger ionospheric correlations). Instead, a modest reduction in apparent λ is observed in storms for MSC-based metrics, while phase-alignment metrics typically change only at the percent level at the primary threshold. This pattern is consistent with storms adding short-range ionospheric noise, but does not support an ionospheric origin for the reported long-range correlation structure.

The clock drift metric (time derivative of clock bias) serves as an internal consistency check. Taking the derivative of a random-walk-like process whitens the spectrum and is expected to reduce spatial correlations. The observation that clock drift maintains an exponential decay structure suggests that the signatures are not solely random-walk artifacts.

Conclusion: The preservation of spatial correlation in the derivative suggests that the signal is not due to static offsets or random walks. The signal represents dynamic, spatially-correlated fluctuations.

A key validation test addresses the question: "Does the methodology itself create artificial exponential decay?" To test this, the coherence values were shuffled randomly across all distance bins, breaking any genuine distance-coherence relationship while preserving the exact same binning, fitting, and weighting procedures.

Methodology

The shuffled null test randomly permutes the coherence values across all station pairs, destroying any physical distance-dependence while preserving:

• Identical log-spaced distance binning (40 bins, 50–13,000 km)
• Identical bin mean calculation
• Identical weighted exponential fitting
• Identical R² assessment

Results

Implications

If log-binning or exponential fitting created artificial decay, the shuffled data would also show exponential decay. It does not. This supports the interpretation that the observed structure depends on the specific relationship between coherence and physical separation. Destroying the spatial topology destroys the signal, which is consistent with a genuine geometric property of the GNSS network.

• Log-spaced binning does not transform Y-values — it only groups data
• The exponential model cannot fit flat data with high R²
• The decay observed in real data is a property of the data, not the analysis

This supports the interpretation that the exponential decay is a property of the data, not a methodological artifact.

If the methodology artificially created exponential decay, all metrics would show identical correlation lengths. Instead, different metrics show distinctly different λ values:

Key observation: The ionofree mode shows λ = 1,072 km, which is 47% longer than baseline (727 km). Removing ionospheric effects increases the correlation length because ionospheric decorrelation was shortening the apparent λ. An artificial methodology effect would not "know" to respond this way to ionospheric correction. This physically meaningful response is consistent with a real signal rather than a fitting artifact.

The "Time Slip" Problem and Its Solution

A critical validation step was verifying that the Pandas DatetimeIndex alignment method correctly handles missing data. Without proper alignment:

• Missing days cause cumulative "time slips"
• Station pairs become desynchronized
• Coherence artificially degrades
• Exponential decay becomes spuriously steep

Validation: The agreement between raw SPP results (using Pandas DatetimeIndex alignment) and precise-product results (using the same methodology in Papers 1 and 2) supports the alignment approach.

Dataset Independence Statement

This analysis (Paper 3) was designed for maximum independence from Papers 1 and 2:

The two analyses share no common processing. Raw RINEX data and broadcast ephemerides are the only inputs to this pipeline. Any agreement between results therefore constitutes independent confirmation, not circular validation.

Cross-Validation Results

Comparison with CODE longspan reference values (from Paper 2):

Interpretation: The agreement between completely independent processing pipelines provides strong evidence that the TEP signal is a physical phenomenon, not a processing artifact. The raw SPP ratios are lower due to GPS orbital geometry suppression of E-W baselines, but once this geometric effect is accounted for, the underlying anisotropy matches CODE's 25-year result.

Why This Is Not Circular Reasoning

A common concern with cross-validation is circularity: is this just comparing results to themselves? This is not the case here:

• Primary evidence requires no correction: The core finding—E-W > N-S at short distances (<500 km)—uses raw, uncorrected values and matches CODE's prediction directly. The geometric suppression analysis (§3.9.4) is secondary validation only, explaining long-distance λ inversions
• Different data: CODE products are network-adjusted; RINEX/SPP is purely local
• Different processing: CODE uses IGS analysis centers; RTKLIB is independent
• Different time periods: 25-year vs 3-year (mostly non-overlapping stations)
• Results are not tuned: Raw RINEX results are reported as computed; comparison is post-hoc

The CODE reference values are used only for interpreting secondary full-distance results, never to adjust or tune the primary short-distance evidence. The short-distance finding (E-W > N-S at <500 km) is fully independent and requires no external calibration.

Table 4.8.1: Summary of validation tests for exponential decay signature

The Processing Artifact Hypothesis — Addressed

Critics of Papers 1 and 2 could reasonably argue that the sophisticated algorithms used by CODE, ESA, and IGS to generate precise products might inadvertently create correlated residuals. These algorithms include:

• Network adjustment with inter-station constraints
• Reference frame alignment procedures
• Common ionospheric and tropospheric models
• Clock constraint strategies

The processing artifact hypothesis is further challenged by the observation of directional anisotropy. The short-distance E-W/N-S ratios of 1.033 (MSC) and 1.224 (Phase Alignment) are highly significant (nominal p < 10⁻¹⁵ under the standard null) and consistent with CODE's directional signature. A critical audit reveals this signal is not a distance artifact: E-W pairs are actually 13 km longer than N-S pairs (a bias against the signal), and robust distance-matching strengthens the ratio (1.033 → 1.041).

Mechanism 1: Geometric Suppression (PDOP Anisotropy)

GPS orbital inclination (55°) creates anisotropic satellite visibility. Pure geometry simulation with PDOP-weighted synthetic clocks yields:

• E-W mean λ: 265 km
• N-S mean λ: 3,994 km
• E-W/N-S Ratio: 0.066 (15× suppression factor)

Key Point: Derived from GPS constellation geometry without any reference to CODE or empirical TEP results. This breaks the circularity argument.

Mechanism 2: Ionospheric Local-Time Decorrelation

E-W station pairs span different time zones, experiencing different ionospheric phases. Simulation with diurnal TEC model (TEC = TEC₀ × [1 + 0.5 × cos(2π × (LST − 14)/24)]) yields:

• E-W mean λ: 1,959 km
• N-S mean λ: 100 km
• E-W/N-S Ratio: 19.6 (20× enhancement, opposite direction)

Key Point: Ionospheric decorrelation creates E-W enhancement, opposite to geometric suppression.

Distance-Dependent Behavior

These mechanisms operate at different distance scales, explaining the observed pattern:

Resolution: The short-distance E-W/N-S ratio (1.20–1.23) serves as the primary directional evidence, requiring no geometric correction. At long distances (>1000 km), both mechanisms are active with ionosphere dominating, creating the observed sign reversal. After correcting for the 15× geometric suppression factor (derived from first principles), the long-distance ratio recovers to 1.46, within 32% of CODE's 25-year benchmark (2.16). This provides secondary validation while maintaining independence of the primary evidence.

Physical Implications of CMB Alignment

• Best-fit RA = 188°, Dec = −5°, only 20.0° from CMB dipole (168°, −7°)—matching CODE's 25-year benchmark of 18.2°
• 78% RA clustering: Of 54 clean (non-Ionofree) combinations, 42 find RA within 10° of CMB (p < 10⁻³⁵ under a simplified binomial model)
• Signal Booster: Aggressive quality filtering (Dynamic-50: daily files with clock std < 50 ns) boosts correlation to r = 0.660 (vs. typical r ≈ 0.51), confirming the signal is physical and high-fidelity.
• Solar Apex disfavored at 86.5° separation (4.3× farther than CMB, 32× worse variance explained)
• Zero-variance filter independence: All three station filters converge to same RA (CV = 0.3%)

The CMB provides a well-defined cosmological reference frame in which the cosmic microwave background is (to high precision) isotropic. Under standard interpretation, the observed CMB dipole arises from Earth's motion relative to this frame (Burde, 2016; Consoli & Pluchino, 2021). If the anisotropy modulation depends on a preferred velocity direction, the CMB dipole is therefore a physically motivated direction to test. The large separation from the Solar Apex direction suggests that any preferred direction inferred from the data is not aligned with the Sun's local galactic motion.

Theoretical Resolution: Bi-Metric Geometry & Local Invariance

The apparent conflict with standard Lorentz invariance can be addressed within the Bi-Metric Geometry framework detailed in the companion theory paper (Smawfield, 2025, TEP theory preprint). The theory postulates that while matter couples to a causal metric $\tilde{g}_{\mu\nu}$ (preserving exact local Lorentz invariance and a locally invariant $c$ in freely falling frames), the global time field $\phi$ induces path-dependent synchronization non-integrability. The "preferred" CMB frame is not an ether, but the natural cosmological rest frame of the scalar time field $\phi$, consistent with the background evolution of the universe. Thus, the signal represents a breakdown of global simultaneity, not local covariance.

Orbital Velocity Coupling — Detected

The orbital velocity correlation (as in Paper 2) has been tested and is detected:

• Best result: r = −0.763, 5.4σ (Multi-GNSS pos_jitter/phase); MSC yields r = −0.610, 4.0σ
• Baseline GPS: r = −0.509, 3.2σ
• Ionospheric independence: Signal persists under ionofree (best: r = −0.416, 2.5σ)
• Seasonal breathing: Equinox/Solstice ratio 1.33–1.58 across all modes

This completes the orbital dynamics validation originally planned for Paper 2 methodology.

TEP Predictions vs. Observations

Primary Results

1. Directional anisotropy detected: E-W correlations are 2–5% (MSC) to 22% (Phase Alignment) stronger than N-S at short distances (<500 km), matching CODE's directional signature (nominal p < 10⁻¹⁵ under the standard null). This finding is robust to distance bias: audit confirms E-W pairs are 13 km longer than N-S (bias against signal), and distance-matched analysis strengthens the ratio (1.033 → 1.041). This supports the "Vanishing Bias" principle: as baseline length approaches zero, distance-dependent atmospheric biases fade away, revealing the true E-W > N-S signal without need for geometric correction
2. Monthly temporal stability: E-W > N-S detected in 94–100% of all 36 months across all processing modes (p = 1.5 × 10⁻¹¹ for 36/36). Multi-GNSS shows strongest effect (phase alignment ratio 1.279). Short-distance ratios show CV = 0.7–1.0% (coherence), 3–6% (phase alignment). The constant short-distance signal combined with the orbitally-modulated full-distance λ ratio (r = −0.509 to −0.763) supports the "Screened Signal Model"—a constant gravitational signal masked by variable atmospheric screening
3. Multi-mode validation: Signal detected in GPS L1 (ratio 1.033), ionofree L1+L2 (1.019), multi-GNSS (1.050), and precise (IGS SP3), suggesting it is not ionospheric, constellation-specific, or caused by broadcast ephemeris errors
4. Geometry-corrected match: After correcting for GPS orbital suppression, E-W/N-S ratios converge to 1.80–1.86, within 17% of CODE's 25-year PPP reference (2.16)
5. Geomagnetic independence (comprehensive): Kp stratification using real GFZ data (primary split Kp<3 vs Kp≥3; 72 tests) shows near-invariance (median Δλ ≈ −1%, with 60/72 tests within ±5%). Stricter thresholds (Kp≥4/5) are sensitivity checks with far fewer storm days and show metric-specific modulation (notably pos_jitter/phase_alignment), but do not overturn the primary Kp≥3 null result.
6. Hemisphere consistency: Both Northern and Southern hemispheres show E-W > N-S across both metrics (coherence: 1.029/1.022; phase alignment: 1.200/1.348), consistent with heliocentric rather than local/seasonal origin
7. Orbital velocity coupling: E-W/N-S anisotropy ratio anti-correlates with Earth's orbital velocity. Multi-GNSS yields the strongest detection (r = −0.763, 5.4σ for pos_jitter/phase; r = −0.610, 4.0σ for MSC), with baseline GPS at r = −0.509, 3.2σ. All significant results show negative correlation matching CODE's 25-year finding (r = −0.888, 5.1σ). Signal persists under ionospheric removal (best ionofree: r = −0.416, 2.5σ), disfavoring an ionospheric origin.
8. Filter consistency: All three station filtering methods (DYNAMIC, OPTIMAL, ALL STATIONS) produce consistent negative correlations across all 6 metric/coherence combinations, suggesting the signal is network-wide and methodologically robust
9. Metric complementarity: MSC detects temporal modulation (orbital coupling: 3.0–4.2σ), phase alignment detects spatial structure (directional anisotropy: 1.35 ratio) and achieves strongest orbital coupling (5.4σ)—different aspects of the same underlying phenomenon
10. Regional control tests (Step 2.1a): Exponential coherence decay is reproduced in Global, Non-Europe, and hemisphere-specific subsets with MSC correlation lengths of order 700–900 km and phase-alignment lengths ≈2–3× larger. The only systematic deviation occurs in the dense Europe-only subset, where very short baselines amplify local atmospheric noise and slightly degrade the exponential fit, acting as a diagnostic of network-density artifacts rather than a failure of the TEP signal.
11. Planetary event modulation: Year-specific coherence modulation detected around 37 planetary conjunction/opposition events with 2.8× higher detection rates than permutation null controls (p < 0.001 for all 6 metrics). Detection rates of 59–68% vs. 20–26% null rate. Mass scaling analysis rules out tidal mechanism: No consistent positive GM/r² scaling observed across 6 channels (5/6 show p > 0.49, 1/6 shows anticorrelation r = −0.42, p = 0.010 opposite to tidal prediction). Despite null mass scaling, detection remains highly significant, indicating a non-tidal, threshold-dependent or geometric mechanism consistent with TEP's prediction of temporal-gravitational coupling distinct from classical tidal forces. Clock Drift MSC shows highest sensitivity (mean σ = 4.25). This independently replicates and strengthens CODE 25-year longspan planetary event findings while ruling out conventional tidal explanations
12. CMB frame alignment: Comprehensive 72-combination full-sky grid search yields results consistent with the CODE 25-year benchmark. The Multi-GNSS/Pos_Jitter/Phase combination produces a best-fit vector (RA=188°, Dec=−5°) that is statistically indistinguishable from the CODE reference (RA=186°, Dec=−4°), with a separation of just 20.0° from the CMB dipole. Aggressive quality filtering (Dynamic-50: daily files with clock std < 50 ns) boosts the correlation to r = 0.660 (vs. typical r ≈ 0.51), confirming the signal is an underlying high-fidelity feature of the data. Of 54 clean combinations, 78% find RA within 10° of CMB (p < 10⁻³⁵). Solar Apex is disfavored (86.5° separation, 4.3× farther, 32× worse variance explained). This provides independent support for CODE's finding that the annual anisotropy modulation is coupled to Earth's motion relative to the cosmic rest frame
13. Seasonal stability (comprehensive): Seasonal stratification analysis reveals three complementary signatures: (1) "Summer Breakthrough" (OPTIMAL_100/Ionofree: λ = 6060 km; confirmed by Precise mode: λ = 6259 km), (2) "Invariant Core" (DYNAMIC_50/Multi-GNSS: λ = 1700–1900 km, Δ < 13% across seasons), and (3) "Network-wide Baseline" (ALL_STATIONS/Baseline: Δ < 8%). The signal is not a seasonal artifact—it is a stable gravitational phenomenon variably screened by the atmosphere. The CODE result (4,201 ± 1,967 km) is statistically consistent with the Annual Ionofree average (~4,170 km) and encompasses the "Summer Breakthrough" within its uncertainty range.
14. Null tests passed (comprehensive): Rigorous validation across 72 independent tests constrains several non-gravitational origins: (1) Solar rotation shows zero correlation (all r < 0.09, 72/72 tests pass), (2) Lunar tides show zero correlation (all r < 0.11, 71/72 tests pass), (3) Shuffle test confirms genuine structure (Real R² = 0.945 vs. Shuffled R² = 0.029 mean, min ratio 1.9×, 90% pass strict R² < 0.3). The signal survives ionospheric removal (Ionofree R² = 0.921), persists across four constellations (Multi-GNSS R² = 0.956), and shows direction-level filter convergence in the CMB frame analysis (CV = 0.3%). A non-gravitational explanation consistent with this complete evidence suite has not yet been identified; within the TEP framework, a gravitational coupling remains a plausible interpretation
15. Statistical significance: t-statistics up to 112.13, Cohen's d up to 0.304, 95% CI excludes unity

Significance of Mass-Independence in Raw Data

The absence of consistent positive GM/r² mass scaling in raw SPP data is a notable finding that supports the TEP interpretation beyond previous PPP-based results:

• No "Filtering" Argument: Unlike PPP, SPP processing does not rigorously model and remove solid Earth tides. If the observed modulation were a residual tidal effect, it should scale with planetary mass and distance ($M/r^3$).
• Distinct Phenomenon: The fact that the signal is detectable but shows no mass dependence is consistent with it being physically distinct from tidal forces.
• Geometric Origin: The modulation depends only on alignment geometry, consistent with TEP's prediction of spacetime metric variations rather than Newtonian force perturbations.

License: This work is licensed under a Creative Commons Attribution 4.0 International License.

Version: v0.3 (Kathmandu) · Last updated: 17 December 2025