Global Time Echoes: Raw RINEX Consistency Test

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 $\lt$ 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; the analyzed 539-station set itself contains 285 NH and 254 SH stations) 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 $\lt$ 50 ns AND
       max_jump $\lt$ 500 ns AND
       total_range $\lt$ 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 Temporal Topology correlation lengths (λ_T). If TEP were an artifact of a specific processing method, different modes would produce different signatures. Outcome: modes do not return similar λ_T — fitted scales range from ~0.6×10³ km (MSC) to ~4.8×10³ km (ionofree phase), so this criterion is not met as stated. The mode dependence is attributed to noise-floor projection in the single-frequency and MSC channels. What stays consistent across modes is the sign of the E-W > N-S ordering and of the orbital-velocity coupling, not the scale.

Outcome: modes do not return similar λ_T. Fitted scales range from ~0.6×10³ km (MSC) to ~4.8×10³ km (ionofree phase), so this criterion is not met as stated. The mode dependence is attributed to noise-floor projection in the single-frequency and MSC channels (§3.7). What stays consistent across modes is the sign of the E-W > N-S ordering and of the orbital-velocity coupling, not the scale.

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 λ_T.
• 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 Temporal Topology correlation lengths are observed (λ_T = 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 testing because:

• Directly measures temporal fluctuations at the receiver
• Shows strongest directional anisotropy (E-W/N-S ≈ 1.22 at short distances)
• Temporal Topology correlation length λ_T ≈ 700–1,100 km (lower bound; see §3.7)

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 ($\lt$ 500 km)

Metric 2: Mean Phase-Offset Alignment Index

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

The resultant magnitude $R = |\sum_f w_f e^{i\phi_f} / \sum_f w_f|$ quantifies the concentration of phase differences. Near-zero resultants make the mean phase angle unstable; a diagnostic of one representative station-day found that 1.5% of pairs had R $\lt$ 0.10. The reported PA statistic retains its established definition and is interpreted as mean phase-offset alignment, not phase-locking strength. 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 indicates a larger directional contrast in the mean phase-offset score than in MSC. Because the archived calculation does not retain the resultant magnitude, it does not by itself establish stronger phase locking.

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.
• Mean phase-offset alignment = cos(arg(Σw·e^iφ)) — Sensitive to the direction of the weighted mean phase vector. Its concentration is R = |Σw·e^iφ|/Σw and is not encoded in PA alone.

MSC measures the strength of a shared signal. Mean phase-offset alignment measures whether the average timing offset points near zero; the resultant magnitude is needed to tell whether that average is sharply defined or produced by cancelling phase values.

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: the PA result reports a directional difference in mean phase offset. A concentration analysis using R can determine whether that difference is accompanied by phase locking.
• 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: MSC and the mean phase-offset score provide complementary descriptive information. PA must not be interpreted as a phase-locking-strength statistic without its resultant magnitude.

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]
• Temporal Topology Correlation length (λ_T): 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 $\lt$ 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² $\lt$ 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 $\lt$ 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 globally; the analyzed set is nearer balance at 285 NH vs 254 SH stations, though NH stations contribute substantially more station-days), making NH seasons the natural stratification choice.

Analysis Methodology

For each season, 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 seasonal estimates across processing configurations (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 reflects a stable terrestrial screening structure: λ should be stable across seasons in high-reliability multi-constellation channels (Δ $\lt$ 15%), with any seasonal excursions in noise-sensitive single-frequency channels bounding propagation-medium artifacts.

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. It diagnoses how long-distance exponential fits behave under global network geometry and varying seasonal noise conditions.
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 directional signature
• Clock bias + MSC: Moderate anisotropy (E-W/N-S ≈ 1.02–1.05) — consistent but weaker
• Clock drift (any coherence): Weak anisotropy (E-W/N-S ≈ 1.07) — derivative preserves structure
• Planetary events (all metrics): Year-specific windows yield 2.8× higher detection than the date-shuffled null (59–68% vs. 20–26%, p $\lt$ 0.001 for all 6 metrics); season-preserving controls (Step 2.6b) attribute the excess to calendar-position sampling of the autocorrelated seasonal series rather than event-locking, and neither the GM/r² acceleration proxy nor the GM/r³ tidal-gradient scaling shows dependence.

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.

The direction-only template is compared with two named fixed-frame candidates and with competitive in-ecliptic orbital controls:

• 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 a fixed 20 km/s phenomenological background term. Different background directions produce different annual phases in the velocity-declination predictor. The scan fixes this magnitude and therefore estimates direction/phase only; it does not estimate a physical background speed.

The CMB dipole remains a physically motivated fixed-frame control, but it is not privileged over the aphelion and perihelion velocity tangents. The corrected results identify the aphelion tangent and reject the CMB-apex phase.

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

Total combinations: 3 × 4 × 3 × 2 = 72 analysis configurations

This exhaustive approach allows identification of which combinations recover the CMB signal most cleanly and assessment of 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.

Orbital-velocity convention: The orbital term implements Earth's physical heliocentric velocity (prograde motion), verified against independent ephemeris checks: at the March equinox the velocity points to ecliptic longitude λ = 270° (equatorial RA ≈ 270°, Dec ≈ −23.4°), at perihelion to RA ≈ 192°, Dec ≈ −5.1°, and at aphelion to RA ≈ 12°, Dec ≈ +5.1°. An earlier pipeline version carried the antiparallel convention; because the predictor cos(dec_net) is invariant under V_net → −V_net, that convention returned the identical correlation landscape mirrored through the origin (verified to 10⁻¹⁵), so the directional results of earlier versions read at the antipode. All reported directions use the corrected convention.

Grid Search Parameters

The 20 km/s background speed is the fixed phenomenological parameter used to construct the directional predictor, not the physical 370 km/s CMB-dipole speed. Because the predictor cos(declination of V_orbital + V_bg) is phase-dominated by the rotating orbital term, the inferred direction re-parametrizes the annual phase of the modulation; the dependence on the assumed speed is therefore measured directly by a |V_bg| scan from 5 to 370 km/s (§3.13.3a), and the best-fit direction is compared against competitive controls that live in the ecliptic plane — where any orbit-driven modulation must lie — rather than only against the Solar Apex 53° off the ecliptic.

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. Configuration-grid interpretation: the Monte Carlo global p-value is reported for the selected configuration; no experiment-wide aggregate significance across the 54 correlated configurations is claimed.

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 directional recovery: 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 $\lt$ 0.05 after look-elsewhere correction
• Different station filters produce inconsistent directions (high variance)
• Baseline and Multi-GNSS modes find different preferred directions

Conversely, support for a genuine frame detection — as distinct from an annual-phase measurement — would additionally require:

• Majority of clean (non-Ionofree) combinations find RA within 20° of CMB
• At least one combination achieves global p $\lt$ 0.05
• Zero variance across station filters (all converge to same RA)
• The CMB dipole outperforms competitive in-ecliptic controls (perihelion/aphelion velocity tangents), not merely the off-ecliptic Solar Apex
• The inferred direction is stable under a |V_bg| scan up to the physical CMB speed of 370 km/s

As reported in §3.13.3a–3.13.9, under the corrected orbital-velocity convention the CMB hypothesis fails outright: the best-fit direction is ~160° from the dipole apex (farther than the Solar Apex at ~93.5°), the CMB-direction template anti-correlates (r ≈ −0.39), the aphelion-velocity tangent outperforms it (r ≈ +0.49, 3.9° from the best fit), and the inferred direction migrates tens of degrees across the speed scan. The result is therefore reported as a well-defined annual phase expressible as an in-ecliptic preferred direction — coincident with Earth's velocity tangent at aphelion — and not as a cosmological-frame detection.

Under the same corrected convention, CODE's 25-year analysis recovers best-fit at RA = 6°, Dec = +4° (162° from the CMB dipole, ~6° from the aphelion velocity tangent). With only 3 years of data, weaker Dec constraints are expected due to limited seasonal sampling, but RA should converge to the orbital-tangent direction if the effect is real.

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

The analysis achieves consistent detection of distance-structured correlations across all 72 correlated metric configurations:

Key insight: Mean phase-offset alignment (labelled "Phase Alignment" in tables and figures throughout) 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 the mean phase-offset score captures directional consistency in average timing offset that persists over longer distances (~1,700–3,500 km); as defined in §2.4.4, it is not a phase-locking-strength statistic without the resultant magnitude.

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 $\lt$ 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 $\lt$ 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 $\lt$ 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 $\lt$ 20% (most $\lt$ 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.

Processing Independence

Raw SPP results (λ = 727–1,072 km) are comparable to precise-product analyses (λ ≈ 3,330–4,549 km, Papers 1, 2, 5), demonstrating that the observed correlation structure is present directly in raw pseudorange observations. The baseline GPS-only mode exhibits a shorter apparent length (727 km) because uncorrected ionospheric fluctuations introduce short-range spatial noise. When removed via dual-frequency linear combination, the longer-range correlation (1,072 km) is recovered. The Multi-GNSS mode (815 km) provides a cross-constellation verification across GPS, GLONASS, Galileo, and BeiDou.

Seasonal Modulation of the Anisotropy Ratio

A static global average of the E-W/N-S anisotropy ratio in raw baseline SPP yields values near 0.95 (suggesting apparent N-S dominance), differing from the CODE longspan finding of E-W dominance (Ratio ~2.16). However, analyzing the directional ratio on a monthly basis reveals a systematic annual oscillation. In baseline single-frequency GPS, the monthly series exhibits an annual cycle superposed on a semiannual modulation with peaks near the equinoxes. In Multi-GNSS phase alignment, this semiannual power is suppressed by a factor of 4 to 6, isolating an overwhelmingly annual harmonic ($A_1/A_2 \approx 3.98$) phase-locked to Earth's orbital velocity and perihelion.

The Three Signatures of TEP

1. The "Summer Enhancement" (OPTIMAL_100/Ionofree): λ ≈ 6060 km — upper-bound seasonal estimate in the noise-amplified ionosphere-free combination, bounding propagation-channel sensitivity
2. The "Core Baseline" (DYNAMIC_50/Multi-GNSS): λ = 1700–1900 km — low-variation baseline across seasons (Δ ≈ 7–13%), recovering a stable processing-conditioned clock-covariance scale
3. The "All-stations Baseline" (ALL_STATIONS/Baseline): λ = 1750–1890 km (Δ $\lt$ 8%) — baseline robustly detectable in the full network across all meteorological seasons

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

Finding: When ionospheric delay 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 consistent in scale with CODE's 25-year directional family (cross-sector mean and dispersion 4,201 ± 1,967 km).

Methodological interpretation and propagation bounds:

• Propagation channel sensitivity: In the dual-frequency Ionofree combination, eliminating first-order dispersive delay inflates measurement noise by a factor of ~3.2, making long-distance exponential fits sensitive to seasonal ionospheric and tropospheric gradient variations
• Kinetic screening vs. atmospheric plasma: Under TEP foundations (Paper 0 §2.2), the environmental screening operator $S_\Sigma(\mathcal{E})$ is governed by the scalar-gradient state $X$ in high-acceleration regimes, while atmospheric plasma ($\rho \sim 10^{-18}\ {\rm g\,cm^{-3}}$) carries negligible scalar charge and gradient; atmospheric weather does not source or screen the gravitational scalar field
• Propagation bounds: The seasonal variation in this specific channel bounds the maximum possible propagation-medium contribution to apparent correlation lengths, defining the observational sensitivity envelope
• Core baseline stability: In contrast to this noise-sensitive channel, the multi-constellation core network (DYNAMIC_50 Multi-GNSS, §3.7.3) exhibits only 7–13% seasonal variation, isolating a stable underlying spatial baseline ($\lambda \approx 1700\text{--}1900\ {\rm km}$) across all seasons

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 demonstrate that the underlying correlation structure ($\lambda \approx 1700\text{--}1920\ {\rm km}$) is robust against seasonal propagation variations, isolating a stable terrestrial baseline independent of atmospheric fluctuations.

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.

Observed Seasonal Stability and Propagation Bounds

The multi-constellation core (DYNAMIC_50 Multi-GNSS) shows a seasonally stable scale of ~1.7–2.0×10³ km (7–13% seasonal variation). The ionofree/OPTIMAL_100 channel extends to ~6×10³ km in summer, which bounds the propagation-driven distortion. No screening-derived value for the terrestrial scale is claimed here.

Seasonal Stratification Summary

Key findings:

• DYNAMIC_50 seasonal variation: Δ $\lt$ 13% across all seasons (Multi-GNSS: Δ $\lt$ 8%)
• OPTIMAL_100 summer enhancement: λ = 6,060 km is consistent in scale with the CODE directional family
• ALL_STATIONS baseline: Signal detectable across network configurations (Δ $\lt$ 12%)
• Physical consistency: Seasonal stability in the core network (Δ $\lt$ 13%) demonstrates resilience against propagation disturbances, while single-channel variations bound the propagation-medium envelope

Conclusion: The seasonal stratification confirms that the multi-constellation core network isolates a stable underlying correlation scale ($\lambda \approx 1700\text{--}1900\ {\rm km}$) across all seasons (Δ ≈ 7–13%), while the larger apparent scale in specialized single-channel configurations (such as OPTIMAL_100 Ionofree in summer) bounds observational noise and geometry variations rather than field-level screening.

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 (Δλ $\lt$ 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 correlated analysis configurations at the primary threshold (3 filters × 4 modes × 3 metrics × 2 coherence types):

• Primary result (Kp $\lt$ 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.

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 observed association between station-pair distance and coherence:

• Discrimination: Real data maintains R² > 0.70 in all 72 tests (mean 0.95); shuffled data yields R² $\lt$ 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² $\lt$ 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 in fit quality after shuffling demonstrates that the exponential structure depends on the observed association between station-pair distance and coherence.

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.8) 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 seasonal analysis (§3.7) and geomagnetic stratification (§3.8), 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

The CODE-referenced sector comparison yields corrected E-W/N-S ratios of 1.80–1.86. This is an empirical comparison, not an independent confirmation against CODE because its correction factor is constructed using the CODE reference.

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 empirical geometry diagnostic provides interpretive context for the full-distance results. The primary evidence (E-W > N-S at short distances, §3.10.1) is independent of this analysis and does not use CODE calibration values. The separately normalized GPS-geometry simulation is reported in §6.3.3.

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 factor is calculated from the sector comparison rather than freely fitted, but CODE enters its construction, so agreement with CODE afterward is not an independent validation:

• 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 is compatible with a stable direction-dependent suppression in this empirical comparison, rather than a mode-specific adjustment.

Suppression Factor Consistency

The suppression factor ranges from 2.42× to 3.16× across modes. The factor is computed directly from the sector-by-sector λ comparison and is not introduced as a free parameter; its use of CODE means it is an empirical, CODE-referenced diagnostic.

Hemispheric Consistency

In the ALL_STATIONS filter, 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). The ALL_STATIONS result is consistent with a heliocentric rather than purely local origin, while recognizing that subset selection can change hemispheric behavior.

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)
• Preferred direction (corrected convention): RA = 6°, Dec = +4° — in-ecliptic, near-equatorial

The hemispheric orbital-coupling asymmetry (Southern stronger than Northern) is corroborated by the CODE longspan analysis. The recovered preferred direction, by contrast, is near-equatorial and in-ecliptic under the corrected orbital convention, so it is neutral on the hemispheric question; it confirms the shared annual phase rather than a southern celestial origin. The directional grid search shares the same annual-phase template as the coupling test, so the pattern is a consistency check across datasets rather than fully independent confirmation; it supports further investigation of regional asymmetries.

Note: The unbalanced pair-epoch counts (53.0M NH vs 8.9M SH) arise from differing observation volume per hemisphere, not station count alone: the analyzed set contains 285 NH and 254 SH stations, but NH stations contribute several times more station-days, so within-hemisphere pair-epochs are dominated by the North. The underlying IGS network bias (238 NH vs 106 SH stations globally) remains a relevant population-level caveat.

Diagnostic: In the higher-quality DYNAMIC_50 subset, hemisphere-stratified short-distance ratios exhibit a Southern Hemisphere inversion (E-W/N-S $\lt$ 1) for the short-distance estimator, motivating the additional falsification tests proposed in §4.

Summary Statistics

• Nominal ≥3σ detections: 5/18 (28%) — 3 MSC, 2 Phase cells (all pos_jitter or ALL_STATIONS clock_bias)
• Nominal 2.5–3σ: 2/18 — pos_jitter+MSC (OPTIMAL_100) and clock_drift+MSC (ALL_STATIONS)
• Nominal values assume 36 independent monthly points. After Bretherton autocorrelation correction (N_eff ≈ 10.5–38.5 across cells), 6/18 remain ≥2.5σ, with best corrected significance ≈3.0σ (pos_jitter/phase, ALL_STATIONS); the balanced DYNAMIC_50 pos_jitter/phase cell reaches r = −0.864, ≈3.3σ corrected
• Direction consistency: 17/18 (94%) show negative correlation, matching CODE's sign; the sole exception (ALL_STATIONS clock_drift/phase) is a null-magnitude cell (r = +0.013)
• Best result: Position Jitter + Phase: r = −0.763, nominal 5.4σ, ≈3.0σ autocorrelation-corrected (N_eff ≈ 12.5)

Position Jitter and Clock Bias

Position jitter and clock bias show comparable orbital coupling on the MSC metric:

• Clock Bias + MSC: r = −0.581, 3.7σ nominal
• Position Jitter + MSC: r = −0.610, 4.0σ nominal

This symmetric response between spatial and temporal observables directly discriminates against downstream clock-internal hardware artifacts: if the signal were an oscillator-specific thermal or circuit effect, it would alter receiver clock bias without perturbing pseudorange geometry. In Single Point Positioning, the state vector is $[X, Y, Z, c\cdot\Delta t]^T$, where $c\cdot\Delta t$ carries units of metres. Because all four parameters are inverted jointly from the pseudorange observation equation $\rho = |\mathbf{r}_{\rm sat} - \mathbf{r}_{\rm rec}| + c\cdot\Delta t + \epsilon$, unmodeled pseudorange perturbations upstream of the navigation solve project with comparable geometric dilution into both coordinate errors and receiver clock bias.

The match holds for MSC in baseline GPS (Δ ≈ 1–9%) and in the Multi-GNSS ALL_STATIONS filter (Δ ≈ 5%), but not for phase alignment or in the other Multi-GNSS filters (e.g. −0.763 vs −0.100 in ALL_STATIONS). Because position and clock are estimated jointly in the navigation solve, correlated projection is expected; the comparison is therefore not used as a test of conformal scaling.

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 $\lt$ 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\cdot\Delta t]^T$ estimates all four unknowns simultaneously, where $c\cdot\Delta t$ carries units of metres. Position jitter and clock bias show similar orbital coupling on the MSC metric in baseline GPS (Δ ≈ 1–9%) and in the Multi-GNSS ALL_STATIONS filter (Δ ≈ 5%), though not for phase alignment or in the other Multi-GNSS filters. Because all four parameters are inverted jointly from the pseudorange observation equation, any unmodeled upstream pseudorange error projects through the geometry matrix into both coordinates and clock bias — correlated projection is expected from the navigation solve itself, so the comparison is not used as a test of conformal scaling. Its retained discriminating content is against downstream clock-internal hardware systematics (such as oscillator temperature drift, which does not alter pseudorange geometry); discrimination against classical upstream systematics is established by the in-phase hemispheric alignment (§3.11.5), multi-constellation resilience, and ionospheric-free persistence.

Hemisphere Phase Discriminator and Harmonic Decomposition (Step 2.8)

The decisive discriminator between a heliocentric driver (Earth's orbital velocity, shared globally) and a local-seasonal driver (insolation, ionospheric season, which is anti-phased between hemispheres) is the relative phase of the annual anisotropy modulation in the two hemispheres. The hemisphere-stratified monthly anisotropy series (DYNAMIC_50, multi-GNSS mode, 36 months) was fit with a constant-plus-annual-harmonic model in each hemisphere and channel. All four channels return annual phases in phase between the hemispheres: clock_bias/MSC Δφ = +52°, clock_bias/phase-alignment Δφ = −81°, pos_jitter/MSC Δφ = −14°, pos_jitter/phase-alignment Δφ = +9° — all within the in-phase sector (|Δφ| $\lt$ 90°) and far from the ~180° required by a local-seasonal driver. The modulation is therefore consistent with a shared heliocentric driver rather than hemisphere-local seasonal forcing, directly supporting the orbital-velocity interpretation of §3.11.

Furthermore, Fourier harmonic decomposition ($A_1$ annual vs. $A_2$ semiannual) reveals that Multi-GNSS phase alignment is overwhelmingly dominated by a pure annual harmonic ($A_1 = 0.159$, $A_1/A_2 = 3.98$, $R^2 = 0.61$). The annual minimum phase-locks to Earth's perihelion (DOY 344.3, offset by only 24.9 days, within the 30.4-day resolution of monthly binning), coinciding with Earth's maximum orbital velocity ($v \approx 30.29\ {\rm km/s}$). In contrast, semiannual equinoctial power is prominent only in single-frequency baseline GPS ($A_2 = 0.226$, peaking within 1.3 to 3.6 days of the equinoxes), where it is readily explained by Russell–McPherron geomagnetic coupling and GPS satellite eclipse seasons; this semiannual power is suppressed by over 5.5-fold under multi-constellation phase tracking.

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

1. Detection at nominal 5.4σ (p = 6×10⁻⁸); ≈ 3.0σ after Bretherton autocorrelation correction (N_eff ≈ 12.5) and ≈ 2.0σ under a conservative 18-cell Bonferroni bound
2. Direction consistency: 17/18 results show negative correlation, consistent with CODE's 25-year finding; the exception is a null-magnitude cell
3. Filter consistency: Consistent negative correlation across three station selection methods
4. Position jitter and clock bias show matched MSC orbital coupling in baseline GPS (Δ ≈ 1–9%) and in the Multi-GNSS ALL_STATIONS filter (Δ ≈ 5%), but not in the other Multi-GNSS filters or on phase alignment; because both are jointly estimated in the navigation solve, correlated projection is expected — this is a hardware-artifact discrimination check, not a conformal-scaling test
5. Metric complementarity: MSC excels at temporal modulation (orbital), phase alignment at spatial structure (anisotropy) and achieves strongest orbital coupling (nominal 5.4σ)
6. Hemisphere and harmonic consistency: the annual anisotropy modulation is in-phase between hemispheres in all four channels (Step 2.8, with pos_jitter/phase_alignment at Δφ = +9°), inconsistent with a local-seasonal driver; Multi-GNSS phase alignment isolates an overwhelmingly annual harmonic ($A_1/A_2 = 3.98$) phase-locked to perihelion (offset ~22–25 days), while equinoctial semiannual power is suppressed by >5.5× relative to single-frequency GPS baseline

Together with the directional anisotropy (Section 3.10), null tests (Section 3.9), 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).

Season-Preserving Null Controls (Step 2.6b)

The permutation null above shuffles coherence values across dates, which destroys the seasonal autocorrelation documented in §3.6–§3.8 while the real events retain their calendar positions on the seasonal curve. A season-destroying null cannot distinguish an event-locked anomaly from ordinary sampling of the seasonal modulation, and the 2022–2024 planetary events are seasonally clustered (outer-planet oppositions fall preferentially in late summer/autumn). Three complementary controls were therefore constructed on the identical daily coherence series (Step 2.6b):

• Uniform-date null: ensembles of 37 fake events drawn at uniformly random dates on the real (autocorrelated) series — the null now samples the true seasonal structure (200 ensembles)
• Calendar-twin control: each real event re-evaluated at the same day-of-year in the other years of the span — seasonal phase held fixed while the planetary configuration is broken
• Deseasonalized re-analysis: annual + semi-annual harmonics fitted and subtracted from the daily series, with real events and the uniform-date null re-scored on the residuals

In addition, the tidal-gradient scaling channel ($GM/r^3$, the leading classical tidal term) is evaluated on the stored per-event JPL distances, closing the untested scaling gap identified in §3.12.6.

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 $\lt$ 0.001, with the year-specific methodology achieving higher detection rates than the previous DOY-pooled approach. However, the shuffled null destroys the seasonal autocorrelation of the daily series; whether this excess is event-locked or a calendar-sampling artifact is adjudicated by the season-preserving controls of §3.12.2a.

Event-Locked Excess Is Not Separable from Seasonal Autocorrelation

Rebuilding the identical daily clock-bias coherence series (2,419,544 valid pair-days; 1,096 days; OPTIMAL_100 baseline) and re-scoring all 37 events under season-preserving controls yields:

Three independent diagnostics converge on the same conclusion:

1. The detector is non-specific on an autocorrelated series. Fake events at uniformly random dates reach a 65.5% "significant" rate on the real MSC series — at or above the real-event rate (59.5%; p = 0.805). The season-destroying permutation null had suppressed this baseline to ~22%, manufacturing the apparent 2.8× excess.
2. The significance pattern is a calendar-position function. Re-evaluating the identical day-of-year in the other two years reproduces the real detection rate exactly (59.5% vs 59.5% for MSC; 55.4% vs 59.5% for phase alignment). The pulse-significance map is therefore locked to calendar position, not to the planetary configuration.
3. No residual excess remains after deseasonalization. On the harmonic-residual series the real events sit at or below the uniform-date null (MSC: 56.8% vs 61.6%, p = 0.735; PA: 62.2% vs 59.5%, p = 0.420).

Consistent with this picture, the event dates do not preferentially sample high-coherence days (event-date placement z = +0.09); the apparent modulation arises because the Gaussian-pulse criterion fires on the majority of an autocorrelated seasonal series rather than because planetary configurations produce anomalous windows. Within this three-year RINEX channel, an event-locked planetary component is therefore not statistically supported: the observed pattern is consistent with calendar-position sampling of the seasonal coherence modulation already documented in §3.6–§3.8 (itself phase-locked to perihelion and in-phase across hemispheres, §3.11).

This negative control outcome does not bear on the independent CODE 25-year planetary analysis (Paper 2), which uses different data, a longer baseline with denser calendar coverage, and its own null construction; application of the same season-preserving control there is required before the corpus-level planetary claim can be assessed.

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). Read together with the season-preserving controls of §3.12.2a, however, the constraints are reinterpreted:

• Tidal rule-out: If any genuine event-locked signal were present and were a residual tidal error, it would be expected to scale with $GM/r^3$. No such scaling is detected (§3.12.6) — though the event-locked component itself is not established in this channel.
• No alignment-locked component detected: Because the season-preserving controls show the detection pattern is a calendar-position function, the data do not establish an alignment-driven (geometric) coupling in this channel. Any mass-independent alignment coupling would constitute new phenomenology beyond the scalar-sector mechanism and would require explicit derivation; none is claimed here.
• Consistent residual reading: the seasonal coherence modulation that the event windows sample is itself the annual, perihelion-locked, hemisphere-synchronous structure attributed in §3.11 to orbital-velocity coupling — the event test is non-specific precisely because it re-samples that documented modulation.

Cross-Validation Significance

The RINEX analysis tests the CODE longspan planetary 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

Under the season-preserving controls, the RINEX channel does not independently confirm an event-locked planetary effect: the elevated detection rate is accounted for by calendar-position sampling of the autocorrelated seasonal series. This constrains, rather than strengthens, the cross-dataset planetary claim pending an equivalent control on the CODE series — while also demonstrating that the null suite can falsify event-locked interpretations where they are absent.

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 a positive acceleration-proxy relationship and not reproduced across other metrics.

Interpretation

The planetary-event analysis finds no consistent positive dependence on the tested GM/r² gravitational-acceleration proxy. Five of six channels show no significant relationship, while one exhibits a negative correlation. This disfavors a simple acceleration-amplitude explanation. The distinct tidal-gradient channel has now been tested directly (Step 2.6b): event σ-level vs $GM/r^3$ yields r = +0.011 (p = 0.948) and |coherence modulation| vs $GM/r^3$ yields r = −0.044 (p = 0.797) across the 37 events — no tidal-gradient scaling is detected, closing the leading classical-amplitude channel.

Key Findings

1. Flat-null excess: All 6 metrics show planetary events with 2.8× higher detection rates than the date-shuffled permutation null (p $\lt$ 0.001 for all)
2. Season-preserving falsification: under a uniform-date null on the autocorrelated series, same-DOY calendar twins, and harmonic-deseasonalized rescoring (Step 2.6b), the event-rate excess does not persist — real 59.5% vs uniform null 65.5% (p = 0.805), twins 59.5%, deseasonalized 56.8% vs 61.6% (p = 0.735) for MSC; the significance pattern is a calendar-position function rather than an event-locked anomaly in this channel
3. No mass or tidal scaling: No consistent positive GM/r² dependence was observed across 6 channels (five null at p > 0.49, one anticorrelation r = −0.42, p = 0.010), and the GM/r³ tidal-gradient channel is null (σ-level r = +0.011, p = 0.948; |modulation| r = −0.044, p = 0.797)
4. Cross-validation status: the RINEX channel does not independently confirm the CODE 25-year planetary-event excess pending an equivalent season-preserving control on the CODE series
5. Interpretation: no alignment-locked (geometric) coupling is established; the event windows re-sample the documented perihelion-locked seasonal modulation (§3.6–§3.8, §3.11), and any mass-independent alignment phenomenology would require explicit derivation before being claimed

This provides a calibrated falsification exercise for the planetary-event channel: the null suite is shown to be capable of distinguishing calendar-position sampling from event-locking, and in this three-year RINEX sample the former accounts for the observations.

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

Under the corrected orbital-velocity convention, the vector (RA = 8°, Dec = +5°) matches the corrected CODE vector (RA = 6°, Dec = +4°) to within ~2.2° — the cross-dataset replication survives the sign correction intact, because both pipelines previously shared the same antiparallel convention. The best-fit direction lies at ecliptic latitude β = +1.4° — in the ecliptic plane, where any modulation driven by Earth's orbital velocity must reside — and sits 3.9° from Earth's orbital-velocity tangent at aphelion (RA = 11.9°, Dec = +5.1°). This is the physically consistent reading: §3.11 established that the directional ratio troughs at perihelion and peaks near aphelion (the r = −0.763 kinematic anti-correlation), so the annual-phase extremum points along the velocity tangent at aphelion. The result 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 $\lt$ 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 highest-correlation filtered configurations have right ascensions ranging from 336° to 352°, clustering near the CMB-antipode RA (348°) — the antipodal image of the previously reported cluster — ~35°–45° from the aphelion velocity tangent and ~140°–145° from the CMB dipole itself.

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. Under the corrected orbital-velocity convention the recovered RA cluster sits antipodal to the direction reported in earlier versions — i.e., near the CMB antipode RA of 348° rather than the dipole apex RA of 168°:

^†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 frame-direction analysis due to noise amplification in the dual-frequency linear combination that degrades directional signal detection.

• Within 5° of the CMB-antipode RA (343°–353°): 24/54 combinations (44%)
• Within 10° of the CMB-antipode RA (338°–358°): 28/54 combinations (52%)
• Within 20° of the CMB-antipode RA (328°–8°): 40/54 combinations (74%)

The 54 combinations share station data, processing modes, and time windows; they are consistency checks, not independent Bernoulli trials. The 52% agreement rate (28/54 within 10° of the cluster RA, versus 5.6% expected by chance for a single draw) demonstrates that the recovered direction is robust to processing choices, though with the scatter quantified below rather than a single tight point. Because every configuration re-samples the same 36 monthly points under the same geometric template, the clustering primarily re-expresses the stability of the shared annual phase: the phase of the aphelion-peaking modulation maps, under the 20 km/s template, to an in-ecliptic direction lying ~20° from the CMB antipode and ~160° from the dipole apex (§3.13.3a). This consistency is reported as a robustness check on the phase measurement, not as an independent combined significance; a binomial p-value assuming independence would be invalid here because the combinations are correlated.

Modal RA Bins at 352°–353°

Twelve configurations land on the adjacent modal bins RA = 352°–353° (Dec ≈ −31°), near the CMB-antipode RA of 348°. The three highest-correlation examples:

The twelve-fold pile-up on two adjacent RA bins would be improbable under an uncorrelated uniform background, but the RA values are deterministic phase readouts of the same annual modulation, so coalescence on a modal bin is expected under the template whenever the phase estimate lands in the corresponding range. The more robust statistic is the broad clustering within 10°–20° of the recovered RA across the full set of clean combinations, read as stability of the annual-phase measurement.

Processing Mode Analysis: RA Distribution by Mode

Convergence is concentrated in the Baseline (78%) and Multi-GNSS (61%) modes; Precise-mode combinations scatter (17%), indicating the directional recovery is mode-dependent rather than uniform across processing choices. Ionofree's low rate (11%) is consistent with the ~3× thermal noise amplification inherent to the L1+L2 combination (Kaplan & Hegarty, 2017).

Statistical Interpretation

The 54 configurations share observations and processing choices and cannot be combined as independent tests. Their agreement demonstrates robustness across analysis settings rather than an independent combined significance. The selected best-fit configuration reports a Monte Carlo global p = 0.026. Nominal pair-level probabilities elsewhere in the paper are retained with their dependence limitation stated explicitly.

The directional template is tested against control directions that are geometrically competitive under an annual-phase predictor — the ecliptic poles, the solar apex and antapex, the Galactic center, Earth's orbital-velocity tangents at perihelion and aphelion, and the CMB dipole and antipode (Step 2.7 corrected convention; step_2_7_cmb_direction_audit.json). Evaluating the same cos(declination) predictor at each control under the fixed 20 km/s template:

Among the competitive controls, the aphelion velocity tangent — the direction of Earth's orbital motion at the aphelion-peaking extremum of the EW/NS annual cycle (the ratio troughs at perihelion, §3.11) — lies 3.9° from the best-fit direction and attains r = +0.494. The CMB antipode attains r = +0.414 while the CMB dipole itself anti-correlates at r = −0.387: the annual phase that earlier versions attributed to the CMB apex is, under the corrected orbital convention, anti-aligned with it. The perihelion tangent returns r = −0.681 with the opposite sign, as required if the template is reading the annual phase rather than a fixed frame direction. All genuinely off-ecliptic controls fail. Under this template the recovered direction is therefore identified with the phase of the orbital-velocity-coupled annual modulation, expressed in celestial coordinates.

Dependence on the assumed background speed

Repeating the grid search over |V_bg| = 5–370 km/s under the corrected orbital convention (step_2_7_cmb_direction_audit.json):

The inferred direction is not stable under the speed scan: it sits near the aphelion velocity tangent (and hence ~20° from the CMB antipode) only in the ~10–20 km/s window, and migrates to ~130° from the dipole apex at the physical CMB speed of 370 km/s — the residual correlation at large |V_bg| (r ≈ 0.8) comes from the transverse projection of the orbital term onto the fixed vector, so the best-fit direction re-parametrizes the annual phase rather than locating a physical frame. The template therefore measures the phase of the annual modulation — expressible as an in-ecliptic preferred direction — rather than detecting the CMB rest frame. The proximity of that direction to the CMB antipode is a geometric consequence of the dipole lying at β ≈ −11°, close to the ecliptic plane.

Fixed-Frame Templates Fail; the Direction is Orbit-Kinematic

This comparison is geometrically asymmetric by construction: the Solar Apex sits at ecliptic latitude β = +53°, far outside the plane in which any orbit-driven annual modulation must live. For any in-ecliptic best-fit direction the Apex is disfavored a priori. Under the corrected convention the same asymmetry also disfavors the CMB dipole differently than earlier versions reported: the best fit is farther from the CMB apex than from the Solar Apex, and the CMB-direction template anti-correlates (r = −0.39). The competitive controls that discriminate — the perihelion/aphelion velocity tangents — are evaluated in §3.13.3a.

For the best-fit direction (RA = 8°, Dec = +5°):

The best-fit direction is nearly perpendicular to the Solar Apex (93.5° separation), excluding a local galactic-motion origin for the modulation phase, and 160° from the CMB dipole — the largest separation of any control except the perihelion tangent itself. The only nearby directions are the aphelion velocity tangent (3.9°) and the CMB antipode (20.0°), both near-ecliptic annual-phase encodings; the orbital reading is preferred because it follows directly from the measured modulation phase rather than from a fixed external frame.

Variance Explained Comparison

• Aphelion-tangent template: r² = 0.24 (positive sign, in phase)
• CMB-dipole template: r = −0.387 — anti-correlated; the dipole template is out of phase with the data
• Solar Apex template: r² = 0.015 (1.5% of annual variance)

All Station Filters Yield Similar Directions

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

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

This low variance across different network geometries shows the inferred RA — equivalently, the annual phase — is not driven solely by a particular station subset or selection criterion.

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

Finding: after the orbital-velocity sign correction applied to both pipelines, the best RINEX vector (RA = 8°, Dec = +5°) still agrees with the corrected CODE vector (RA = 6°, Dec = +4°) to within ~2.2°. This convergence is a genuine cross-dataset consistency check — different data (raw SPP vs precise PPP clocks) and a different time baseline (3 vs 25 years) recover the same in-ecliptic direction — and it survives the sign fix because both pipelines previously shared the same antiparallel convention, so both results mirrored through the origin together. It is not an independent validation of a CMB-frame origin: both pipelines apply the same cos(declination) template to data modulated by the same annual orbital geometry, so both are expected to recover the same annual phase and hence the same in-ecliptic direction regardless of which frame hypothesis motivates the template.

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 still recovers a direction within ~10° of the aphelion-tangent region (antipodal to the CMB-adjacent direction of earlier versions):

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 directional 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. Under the corrected orbital-velocity convention, however, the CMB hypothesis fails at the directional level: the best-fit direction is 160° from the dipole apex, and the CMB-direction template anti-correlates with the data (r = −0.39). The recovered direction instead coincides with Earth's orbital-velocity tangent at aphelion — the heliocentric frame the rest of this analysis already implicates (§3.11).

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 direction (RA = 8°, Dec = +5°) is 93.5° from the Solar Apex
• The Solar Apex template returns r = −0.121 (r² ≈ 0.015) against the best fit's r² = 0.25
• The CMB-dipole template fares worse still on sign: r = −0.387 (anti-phase)

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

Under the corrected convention the clustered combinations show systematic negative Declination (−9° to −46°), while the best primary-combination fit lands at Dec = +5°. Declination is the weakly identified coordinate for two reasons:

• 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 — and the predictor cos(dec_net) is nearly even in the background declination, so the hemisphere of the fixed vector is intrinsically weakly constrained (median |Δr| = 0.09 between north/south-mirrored trial directions).

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 NH weighting of the analyzed pair sample (285 NH vs 254 SH stations; 53.0M vs 8.9M pair-epochs), the apparent modulation amplitude is compressed in the vertical direction. The corrected CODE vector (Dec = +4°) shows the recovered direction lying in the ecliptic plane on the 25-year baseline; the RINEX RA convergence (~2.2° from the corrected CODE RA) with just 3 years provides a comparatively tight RA estimate despite the shorter baseline.

What the Directional Analysis Establishes

1. RA Clustering: 28/54 clean combinations (52%) find RA within 10° of the cluster center (RA ≈ 350°), demonstrating robustness to processing choices
2. RA Matches: Twelve configurations land on the adjacent modal bins RA = 352°–353° (Dec ≈ −31°), near the CMB-antipode RA of 348° (match probabilities depend on the assumed background model and correlations among configurations)
3. Filter Independence: Zero variance across all station filters (CV = 0.3%)
4. Frame identification: the recovered direction lies 3.9° from Earth's aphelion velocity tangent; fixed-frame candidates fail — Solar Apex 93.5° (r = −0.12), CMB dipole 160.0° (r = −0.39, anti-phase)
5. CODE Replication: within ~2.2° of the corrected CODE vector (RA = 6°, Dec = +4°) under the same corrected convention

Criteria Assessment

The CMB-direction criterion fails under the corrected orbital-velocity convention: the best-fit direction is nearly antipodal to the CMB dipole (160.0°) and its template anti-correlates with the data. What survives is a stable annual-phase measurement — an in-ecliptic preferred direction coincident with Earth's aphelion velocity tangent — together with a significant exceedance of the iid-shuffle look-elsewhere null. The grid search uses a fixed 20 km/s phenomenological background-speed template; it does not estimate Earth's physical 370 km/s CMB-dipole speed.

Conclusion: Cross-Dataset Confirmation of the Orbital-Phase Direction

This analysis provides a cross-dataset consistency check on the CODE longspan directional result 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). Under the corrected orbital-velocity convention, the agreement between the two approaches in the inferred direction—within ~2° of each other's recovered vectors—establishes that the same annual-phase structure is recovered in both datasets. Because the two pipelines share the same cos(declination) template and the same orbital geometry, the agreement does not independently validate any particular frame interpretation; it confirms that the measured annual phase is stable across data sources. Under the corrected convention the phase maps to Earth's aphelion velocity tangent — consistent with the kinematic picture of §3.11, in which the EW/NS ratio troughs at perihelion and peaks near aphelion — and is inconsistent with the CMB dipole direction.

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:

• λ_T changes only slightly: Storm λ_T is only 3–4% shorter than quiet λ_T (MSC metrics)
• Phase alignment changes little: Phase-alignment λ_T is typically stable at the percent level at the primary threshold (Kp $\lt$ 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 λ_T = 1,072 km, which is 47% longer than baseline (727 km). Removing ionospheric effects increases the Temporal Topology correlation length because ionospheric decorrelation was shortening the apparent λ_T. 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 (both observe the same physical tracking network). Raw RINEX data and broadcast ephemerides are the only inputs to this pipeline. Any agreement between results therefore constitutes confirmation on an estimation-chain-independent axis, 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 ($\lt$ 500 km)—uses raw, uncorrected values and matches CODE's prediction directly. The geometric suppression analysis (§3.10.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 $\lt$ 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 $\lt$ 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), the geometry simulation gives a 15× suppression estimate and a corresponding corrected ratio of 1.46. This simulation is independent of CODE, whereas the 2.42–3.16 empirical factors and 1.80–1.86 ratios are CODE-referenced. They are complementary diagnostics, not interchangeable estimates of a unique correction.

Physical Implications of the In-Ecliptic Direction

• Best-fit RA = 8°, Dec = +5°, in the ecliptic plane (β = +1.4°), 3.9° from Earth's aphelion velocity tangent (11.9°, +5.1°) and 160.0° from the CMB dipole — within ~2.2° of the corrected CODE 25-year vector (RA = 6°, Dec = +4°) under the shared template
• 52% RA clustering: Of 54 clean (non-Ionofree) combinations, 28 find RA within 10° of the cluster (RA ≈ 350°), demonstrating robustness to processing choices (not an independent combined significance)
• Signal Booster: Aggressive quality filtering (Dynamic-50: daily files with clock std $\lt$ 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 93.5° separation — a geometrically predetermined outcome for any in-ecliptic annual signal, since the Apex sits 53° off the ecliptic (see §3.13.3a for the discriminating controls)
• Zero-variance filter independence: All three station filters converge to same RA (CV = 0.3%)

The CMB dipole and Solar Apex are physically motivated fixed-frame controls. Neither reproduces the measured phase under the corrected convention: the CMB-apex template anti-correlates and the Solar-Apex template is off-ecliptic and weak. The aphelion-velocity tangent is the competitive control selected by the data.

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 Paper 0 (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. A cosmological rest frame may exist in the broader theory, but this directional analysis does not detect it. The measured phase is tied to heliocentric orbital motion, and the CMB-apex template is rejected under the corrected convention. Any future claim of ambient-field rest-frame coupling requires a template that responds to the field gradient rather than the velocity declination, or a baseline long enough to separate the annual phase from fixed in-ecliptic directions.

Orbital Velocity Coupling — Detected

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

• Best result: r = −0.763, nominal 5.4σ, ≈ 3.0σ corrected (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σ)
• Harmonic decomposition and perihelion phase-locking: Multi-GNSS phase alignment isolates an overwhelmingly annual harmonic ($A_1/A_2 \approx 3.98$, $R^2 = 0.61$) whose minimum phase-locks to Earth's perihelion (offset ~22–25 days), while semiannual equinoctial power is suppressed by >5.5× relative to single-frequency baseline GPS (which reflects Russell–McPherron and eclipse dynamics)
• Hemispheric in-phase consistency: All four tested channels are in-phase between hemispheres (Step 2.8, with pos_jitter/phase_alignment exhibiting Δφ = +9°), decisively refuting local insolation and instrumental geometric blind spots

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% (mean phase-offset alignment) stronger than N-S at short distances ($\lt$ 500 km), matching CODE's directional signature (nominal p $\lt$ 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). At short baselines, distance-dependent atmospheric and geometric confounds are substantially reduced, enabling a less biased estimate of local anisotropy without applying geometric correction in the primary estimator.
2. Monthly temporal stability: Month-by-month short-distance anisotropy shows stable polarity (E-W > N-S) at the 94–100% level across modes and metrics (worst case 34/36 months). In the ALL_STATIONS monthly analysis, Multi-GNSS shows the strongest mean effect (phase alignment ratio 1.314). Short-distance ratios show low CV (~0.7–1.0% for coherence; ~3–6% for phase alignment). The stable short-distance signal combined with the orbitally-modulated full-distance λ_T ratio (r = −0.509 to −0.763) supports a stable covariance baseline with orbital modulation; the transfer from the solved scalar field to that covariance and its processing-dependent response remains to be derived.
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 diagnostics: The CODE-referenced sector comparison yields corrected E-W/N-S ratios of 1.80–1.86, while an independent GPS-geometry simulation gives a larger suppression (~15×) and a corrected ratio of 1.46; both are supplementary comparisons, as the primary short-distance anisotropy uses no correction
5. Geomagnetic independence (comprehensive): Kp stratification using real GFZ data (primary split Kp$\lt$ 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 stratification (interpretive diagnostic): In the ALL_STATIONS anisotropy analysis, both Northern and Southern hemispheres show E-W > N-S (coherence: 1.029/1.022; phase alignment: 1.200/1.348). Higher-quality subsets motivate additional hemisphere-controlled falsification tests to assess sensitivity to network geometry and local-time decorrelation.
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, nominal 5.4σ for pos_jitter/phase, ≈ 3.0σ autocorrelation-corrected; r = −0.610, nominal 4.0σ for MSC, ≈ 2.5σ corrected), with baseline GPS at r = −0.509, nominal 3.2σ (nominal significances assume 36 independent monthly points). All significant results show negative correlation matching the sign of CODE's 25-year finding (r = −0.888; 3.65σ autocorrelation-corrected, N_eff ≈ 11). Signal persists under ionospheric removal (best ionofree: r = −0.416, 2.5σ), disfavoring an ionospheric origin. Fourier harmonic decomposition confirms that Multi-GNSS phase alignment is overwhelmingly dominated by a pure annual harmonic ($A_1/A_2 \approx 3.98$, $R^2 = 0.61$) phase-locked to perihelion (offset ~22–25 days), while semiannual equinoctial power is suppressed by >5.5× relative to single-frequency baseline GPS (which reflects Russell–McPherron and eclipse dynamics); the modulation is verified to be in-phase across hemispheres (Step 2.8, Δφ = +9°), refuting local insolation or instrumental geometric blind spots.
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: nominal 2.5–4.1σ; ≤ 2.9σ autocorrelation-corrected), phase alignment detects spatial structure (directional anisotropy: 1.35 ratio) and achieves strongest orbital coupling (nominal 5.4σ; ≈ 3.0σ corrected)—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 Temporal Topology 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 channel (falsified by season-preserving controls): Event windows show 2.8× higher detection rates than the date-shuffled permutation null (p $\lt$ 0.001 for all 6 metrics), but the season-preserving controls of Step 2.6b reassign the excess to calendar-position sampling of the autocorrelated seasonal series — uniformly random dates score at or above the real rate (65.5% vs 59.5%, p = 0.805), same-DOY calendar twins reproduce the real rate exactly (59.5%), and deseasonalized rescoring leaves no excess (56.8% vs 61.6%, p = 0.735). An event-locked planetary component is therefore not established in this channel. No consistent positive dependence is found on either the GM/r² acceleration proxy (five of six channels null, one anticorrelation r = −0.42, p = 0.010) or the GM/r³ tidal-gradient scaling (r = +0.011, p = 0.948).
12. In-ecliptic preferred direction: Under the corrected orbital-velocity convention, the full-sky analysis identifies a best-fit direction at RA = 8°, Dec = +5° (β = +1.4°), 3.9° from Earth's aphelion velocity tangent, 160° from the CMB dipole, and within ~2° of the corrected CODE longspan vector. Across the 54 non-ionofree configurations, 28 recover right ascension within 10° of the cluster center (RA ≈ 350°, near the CMB antipode rather than the dipole), with convergence concentrated in the Baseline and Multi-GNSS modes. These configurations share observations and processing choices; their agreement demonstrates robustness of the annual-phase measurement across analysis settings rather than independent combined significance. Competitive controls (§3.13.3a) show the aphelion tangent out-correlates the CMB dipole under the template (r = +0.49 vs −0.39) and that the inferred direction migrates tens of degrees across the 5–370 km/s speed scan, so the result is reported as an annual-phase measurement expressible as an in-ecliptic preferred direction — identified with Earth's heliocentric-orbital motion, not a cosmological rest frame. The selected best-fit configuration reports an iid-shuffle global p = 0.026.
13. Seasonal stability (comprehensive): Seasonal stratification provides 48 seasonal estimates across processing configurations (4 seasons × 3 filters × 4 modes) for each metric/coherence combination.
14. Null tests passed (comprehensive): Rigorous validation across 72 correlated analysis configurations constrains several non-gravitational origins. Pair-level and configuration-level results are reported as consistency checks, not independent aggregate tests.
15. Statistical significance: t-statistics up to 112.13, Cohen's d up to 0.304, 95% CI excludes unity

Status of the Planetary-Event Channel

The season-preserving controls resolve the interpretation of the event-window detections in this channel:

• Null specificity: On the autocorrelated daily series the Gaussian-pulse criterion is non-specific — uniformly random dates are scored "significant" at 65.5%, at or above the real-event rate. The 2.8× flat-null ratio reflected the season-destroying null, not event-locking.
• Calendar lock: Identical day-of-year placement in other years reproduces the real detection rate exactly, so the significance pattern is a function of calendar position on the seasonal modulation rather than of the planetary configuration.
• Scaling bounds: Neither the GM/r² acceleration proxy nor the leading classical tidal channel GM/r³ shows scaling; no mass-dependent or alignment-locked mechanism is established by these data.

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

Version: v0.8 (Kathmandu) · Last updated: 30 September 2026 · Originally published: 9 December 2025