Global Time Echoes: MGEX Multi-GNSS Clock Replication, 2025–2026

1. Introduction

Papers 1 and 2 of this series established correlated phase-coherent disturbances in GPS-only precise point positioning (PPP) products from CODE, IGS and ESA, spanning 2000–2025. The disturbances exhibit a spatial correlation length of order thousands of kilometres, an east–west anisotropy exceeding the north–south counterpart, coupling to Earth orbital velocity, and a preferred axis near the CMB dipole. Paper 3 moved to raw RINEX single-point positioning (SPP), demonstrating that the same signatures persist without reliance on analysis-centre orbit or clock products. Between them, those papers already provide multi-centre, multi-decade, raw-observation, and GPS-only product-family independence. The principal forms of robustness they do not yet supply are independence from the historical training epoch and independence from a combined multi-GNSS clock product family.

This paper (Paper 14) addresses those two gaps. The analysis uses the public MGEX combined multi-GNSS receiver-clock solution—the CODE COD0MGXFIN CLK product distributed by NASA CDDIS—for the window 2025-01-01 to 2026-05-01. The present MGEX analysis uses a combined multi-GNSS clock product distinct from the GPS-only PPP and raw RINEX products analysed earlier, and a largely held-out 2025–2026 window, though it is not fully temporally disjoint from all earlier analyses. The MGEX clock product is generated by combining GPS, GLONASS, Galileo and BeiDou observations into a single receiver-clock estimate per station, using software, orbit and combination strategies that differ from the GPS-only PPP of Papers 1–2. It therefore constitutes an independent combined-clock data product analysed over a largely held-out epoch. Receiver-clock offsets for 256 globally distributed stations are read directly from the daily 5-minute CLK files; no positioning is performed.

We emphasise at the outset what this product tests and what it does not. Because the MGEX clock files contain one combined multi-GNSS solution per station, they do not permit an independent per-constellation comparison; that form of independence is the province of raw per-system processing, which lies outside the scope of this combined-product analysis. The independence offered here is temporal (a largely held-out epoch) and product-type (a combined multi-GNSS clock from a different analysis pipeline), not per-constellation.

2. Data and Product Scope

2.1 MGEX Combined Clock Product

All data are public and freely available from NASA CDDIS (cddis.nasa.gov) and the International GNSS Service (IGS). The analysis uses the Multi-GNSS Experiment (MGEX) combined receiver-clock product COD0MGXFIN produced by the Center for Orbit Determination in Europe (CODE), with the Wuhan (WUM) and CNES (GRG) MGEX products as fallbacks on the few days CODE is unavailable (WUM 12 days, GRG 2 days out of 486 calendar days in the window, i.e. ≈3 % fallback usage). These are daily clock (CLK) files at a 5-minute sampling interval. Each file reports a single combined multi-GNSS receiver-clock estimate per station, derived by CODE from GPS, GLONASS, Galileo and BeiDou observations; the combined solution is the quantity analysed. No positioning is performed and no raw RINEX observation files are used.

2.2 Analysis Window and Station Set

The analysis period is 2025-01-01 to 2026-05-01 (473 days with usable data out of 486 calendar days in the window), yielding approximately 16 months of daily files. The window is chosen to provide a largely held-out 2025–2026 epoch, though it is not fully temporally disjoint from Papers 1–2, which span 2000–2025. Receiver-clock offsets are extracted from the AR (receiver-clock) records of each CLK file for 256 globally distributed stations, converted to nanoseconds, and assembled into per-station daily time series.

2.3 Product Limitation: No Per-Constellation Separation

Unlike Papers 1–3, which can isolate individual systems, the MGEX CLK product contains one combined receiver clock per station rather than per-constellation clocks. The four constellations therefore cannot be separated in this product, and the analysis is performed once on the combined solution. The independence offered here is temporal (a largely held-out epoch) and product-type (a combined multi-GNSS clock from a different analysis pipeline), not per-constellation.

2.4 Evaluated Signatures

The analysis evaluates the following signatures, consistent with the approach in earlier papers: correlation length, azimuthal anisotropy, orbital-velocity coupling, CMB-frame alignment, ionospheric independence, and geometric robustness.

3. Methods

3.1 Receiver-Clock Extraction

For each daily CLK file the AR records are parsed, which report the combined multi-GNSS receiver-clock offset for each station at the 5-minute sampling interval. The offsets are converted from seconds to nanoseconds and ordered in time, giving one clock time series per station per day. No positioning is performed and no coordinate residuals are computed; the receiver clock is the sole observable.

3.2 Spectral Pre-processing

Each daily series is linearly detrended to remove the dominant clock drift. Cross- and auto-spectra are estimated with Welch's method (Hann window, segment length up to 96 samples). Spectral estimates are restricted to the nominal TEP-sensitive band [10 µHz, 500 µHz], but with the daily Welch segmentation used here the effective lowest non-zero resolved frequency is set by the segment length, approximately 35 µHz for 96 five-minute samples. The practical analysed band is therefore approximately 35–500 µHz, corresponding to periods from about 8 hours down to 33 minutes, which lies safely within the 5-minute Nyquist limit.

3.3 Spectral Phase-Alignment Metric

For every station pair the magnitude-weighted circular mean is computed of the cross-spectral phase across the band. The phase-alignment metric is the cosine of this weighted mean phase, ranging from −1 (anti-phase) through 0 (random) to +1 (in-phase). This metric replaces the spectral-magnitude coherence used in earlier papers; it isolates the sign and stability of the inter-station phase relationship and is the quantity carried through all subsequent tests. Pairs are binned by great-circle distance into 40 logarithmically spaced bins from 50 km to 13,000 km, with a minimum of 10 pairs per bin.

3.4 Correlation-Length Model

The distance-binned phase alignment is modelled as

where A is the amplitude, λ is the correlation length, and C₀ is an incoherent offset. Fitting uses bounded non-linear least squares (Trust Region Reflective), weighting each bin by the inverse of its standard error of the mean (1/SEM, where SEM = σ_bin / sqrt(count)). Fit quality is quantified by R², and the uncertainty on λ is propagated from the covariance matrix.

3.5 Anisotropy, Orbital, CMB, Ionospheric, and Null Tests

The EW > NS test fits λ separately in eight 45°-wide azimuth sectors and compares the mean of the east and west sectors with the mean of the north and south sectors. Because ionospheric decorrelation inverts the anisotropy at long baselines, a longitude-matched subset is also reported (station-pair longitude difference below 30°), which, following Paper 3, isolates pairs at similar local solar time. Significance and confidence intervals for the ratio are obtained from 500 bootstrap resamples. For anisotropy, pair-bootstrap results are reported as sensitivity tests, while station/spatially-clustered resampling is treated as the more conservative inference because station pairs are not independent. The ratio is additionally reported in three distance strata (0–500, 500–1,000, >1,000 km).

The orbital-coupling test correlates the monthly λ and the monthly EW/NS ratio with the projection of Earth's orbital velocity onto the CMB dipole direction, computed from a sinusoidal annual model; the two correlations are Bonferroni-corrected. Because exponential-fit noise on short baselines can degrade these metrics, a supplementary PA-difference metric is also computed: the monthly mean EW minus NS phase alignment, which avoids fitting entirely and directly measures anisotropy strength. The CMB-frame test performs a full-sky grid search for the anisotropy axis on a 37 × 19 grid in right ascension and declination (10° steps), correlating the daily EW/NS ratio with the projection of Earth's orbital velocity onto each candidate axis, with a permutation null and a look-elsewhere correction factor of 50. Ionospheric controls stratify the fit by geomagnetic Kp index and exclude storm days; geometry controls repeat the fit on hemisphere-balanced and distance-matched subsets; and four null tests (temporal shuffle, spatial shuffle, phase randomisation, and solar-rotation reassignment) verify that the recovered structure vanishes under label permutation.

4. Results

The analysis comprises 1,753,922 station pairs drawn from 256 stations over 473 days. Because the MGEX product is a single combined solution, the results are reported once rather than per constellation. The analysis evaluates the core correlation-length signature, directional anisotropy, orbital-velocity coupling, CMB-frame alignment, and ionospheric/geometric robustness controls. The correlation length, anisotropy in the longitude-matched subset, and ionospheric persistence are observed. Orbital-velocity modulation is not recovered by the primary monthly λ or EW/NS metrics, but is recovered by the supplementary PA-difference metric. The CMB-frame test detects a significant anisotropy axis that lies 92° from the CMB dipole.

4.1 Correlation Length

The isotropic fit gives a well-constrained correlation length.

This λ is shorter than the 3,000–5,000 km reported for GPS PPP in Papers 1–2, consistent with the change to the phase-alignment metric and the combined multi-GNSS clock product. The quoted λ uncertainty is the formal covariance-based fit uncertainty for the binned exponential model and does not include all sources of network, product, or model-selection uncertainty.

4.2 EW / NS Anisotropy

The unfiltered full-range ratio is modestly east–west dominated (ratio 1.23) but does not reach significance under the station-clustered bootstrap (p = 0.48). A suggestive east–west excess emerges once pairs are matched in longitude (ratio 2.28, pair-bootstrap p = 0.002), though the significance weakens under spatial-clustered resampling (95% CI [0.16, 14.12], p = 0.244). The 500–1,000 km stratum also shows the expected excess.

4.3 Orbital-Velocity Coupling

The traditional monthly λ and EW/NS ratio do not correlate significantly with Earth's orbital speed or its projection onto the CMB dipole (Bonferroni p > 0.5 for both). These metrics are degraded by the exponential-fit noise on short baselines. A supplementary PA-difference metric (monthly mean EW − NS phase alignment, which avoids exponential fitting) recovers a negative correlation with orbital speed (r = −0.670, p = 0.017) and with the velocity projection onto the CMB dipole (r = −0.770, p = 0.003).

4.4 CMB-Frame Alignment

The best-fit anisotropy axis is RA = 60°, Dec = −60° (r = 0.161), with its antipode at RA = 240°, Dec = 60°. The closer pole lies 92° from the CMB dipole. The axis is significant (LEE p < 0.0001), so the data do contain a preferred anisotropy axis, but it points away from the CMB direction. Widening the azimuth sectors from 45° to 60° increases day inclusion from 298 to 300 and shifts the best-fit to RA = 140°, Dec = −60°, leaving the CMB separation at 57° (or 123° from the antipodal direction). This is more consistent with ionospheric or product-geometry contamination than with a stable TEP/CMB-frame effect. The result is sensitive to sector width: 56% of pairs fall outside the 45° EW/NS sectors and are discarded from the daily ratio computation.

4.5 Ionospheric and Geometry Controls

The correlation length persists on geomagnetically quiet days ($\lambda = 1356$ km, Kp ≤ 2, 688,537 pairs). Active days (Kp ≥ 5, 39,000 pairs) give a longer $\lambda = 1943$ km; the smaller sample increases the variance and the value awaits confirmation with additional active-period data. The distance-matched subset yields a comparable $\lambda = 1847$ km (R² = 0.677), indicating the signal is not a product of network geometry. Hemispheric $\lambda$ differs north to south, as expected from the uneven station distribution; the distance-matched control is the geometry-robust comparison.

4.6 Null Tests

Null tests were run using the same MGEX combined-clock phase-alignment pipeline. In all four cases the fitted distance-decay structure collapsed to R² ≈ 0, confirming that the recovered correlation length depends on the true temporal, spatial and phase information. Fitted λ values under null conditions are not interpreted because the exponential model has no explanatory power when R² ≈ 0.

5. Discussion

The MGEX combined-clock product shows a correlation length of λ = 1396 ± 90 km, which survives the ionospheric and geometry controls and all four null tests, while collapsing to negligible structure under label permutation. That a CODE MGEX combined-clock product, distinct from the GPS PPP products analysed earlier, returns a correlation length of the same order as earlier papers is the central positive result of this analysis. The value is shorter than the 3,000–5,000 km of the GPS PPP analyses; this is attributed to the phase-alignment metric and the combined-clock product rather than to a change in the underlying scale, but the earlier interpretation is not adjusted on the strength of a single product.

The anisotropy result requires candour. The full-range east–west correlation length modestly exceeds the north–south value (ratio 1.23), but the difference is not significant under the station-clustered bootstrap (p = 0.48). The longitude-matched subset shows a large EW/NS ratio under pair bootstrap (ratio 2.28, p = 0.002, 95% CI [1.20, 2.69]), but the effect is not significant under spatial-clustered resampling (95% CI [0.16, 14.12], p = 0.244). We therefore treat the MGEX anisotropy result as suggestive rather than confirmatory. The non-significant full-range ratio is reported alongside the longitude-matched result for completeness. The orbital-velocity coupling is not recovered by the traditional monthly λ or EW/NS ratio metrics when tested against either scalar speed or the velocity-vector projection onto the CMB dipole (Bonferroni p > 0.5). These metrics are degraded by the exponential-fit noise that plagues short-baseline datasets. The PA-difference metric — monthly mean EW minus NS phase alignment, which avoids exponential fitting — is used because the ~16-month data window provides too few monthly points for stable λ-versus-speed regression. It recovers a negative correlation with orbital speed (r = −0.670, p = 0.017) and with the velocity projection onto the CMB dipole (r = −0.770, p = 0.003). The traditional monthly λ and EW/NS ratio remain the standard for longer time series; here they are inconclusive (Bonferroni p > 0.5), consistent with the limited temporal leverage. Because the PA-difference metric was adopted only after the primary λ and EW/NS tests returned null, the coupling it recovers is suggestive and requires independent replication before it can be treated as confirmatory.

The CMB-frame test detects a significant anisotropy axis (LEE p < 0.0001) at RA = 60°, Dec = −60° (antipode RA = 240°, Dec = 60°), but the closer pole lies 92° from the CMB dipole. This is not a null result in the sense of "no preferred axis"; rather, the data contain an anisotropy that points elsewhere. The data favour ionospheric contamination of the daily EW/NS ratio: the axis direction and sector-width instability are more consistent with ionospheric or product-geometry contamination than with a stable CMB-frame alignment. Widening the azimuth sectors to 60° shifts the best-fit to RA = 140°, Dec = −60° and reverses the correlation sign (r = −0.283 versus +0.161 for 45°), leaving the CMB separation at 57° (or 123° from the antipodal direction). The sign reversal indicates that the anisotropy pattern is unstable with respect to sector definition. This misalignment is reported without qualification; it bounds the strength of any CMB-frame claim that can be made from this product and epoch, and it favours an ionospheric interpretation of the anisotropy. The permutation null evaluates significance at the direction that maximises observed correlation, with a look-elsewhere factor of 50 as an approximation; it does not re-run the full grid search on each permutation, so the reported LEE p-value is a lower bound on the true family-wise error.

A satellite-clock analogue was also attempted using SP3 orbit geometry, but the MGEX combined solution estimates all satellite clocks relative to a single reference time scale; after per-epoch common-mode removal the residuals are dominated by estimation noise and the exponential model is actively rejected (R² < 0, e.g. GPS R² = −0.38), which is stronger evidence of no correlation than R² ≈ 0. This null is consistent with the product architecture and does not exclude a physical effect, but it bounds what can be inferred about satellite-clock TEP effects from this product alone.

The structural constraint is that the MGEX clock files provide one combined multi-GNSS solution per station. Per-constellation independence, which the earlier framing envisaged, is not achievable with this product because the four constellations are not separable in the MGEX combined solution. Genuine per-system replication remains the domain of raw observation processing, which the data volume places outside the present scope. What this paper adds is product-type independence and a largely held-out epoch, and within those bounds the core correlation-length signature replicates.

6. Conclusions

An analysis of Temporal Equivalence Principle signatures has been carried out using the public MGEX combined multi-GNSS receiver-clock product (CODE COD0MGXFIN) over 2025-01-01 to 2026-05-01, using a data product independent of the GPS PPP family and a largely held-out epoch. The primary positive result is a reproducible receiver-clock correlation length, λ = 1396 ± 90 km, R² = 0.486, which persists under ionospheric and geometry controls and collapses under four null tests.

Directional signatures are more mixed. The longitude-matched anisotropy is suggestive but not significant under spatial-clustered resampling. Orbital-velocity modulation is recovered only post-hoc with a short-baseline-robust PA-difference metric and requires independent replication. The CMB-frame test detects a significant anisotropy axis, but it lies 92° from the CMB dipole, favouring ionospheric or product-geometry contamination over a stable CMB-frame alignment. A satellite-clock analysis finds the exponential model actively rejected (R² < 0), consistent with the single-reference-time nature of the MGEX combined solution.

The independence offered by this analysis is of epoch and of data product, not of constellation: the MGEX clock product is a single combined solution, so the four systems cannot be separated, and per-constellation replication remains the province of raw-observation processing addressed elsewhere in the series. Within those bounds, the robust observation of the core correlation length in an independent product and largely held-out epoch, together with suggestive longitude-matched anisotropy and supplementary PA-difference orbital modulation, provides additional context for the earlier results. The analysis pipeline and data-provenance records are open-source and reproducible at the project repository.

7. Data and Code Availability

All data are publicly available from NASA CDDIS (cddis.nasa.gov) and the IGS (igs.org) under their standard open-data policies. The analysis uses the MGEX combined multi-GNSS receiver-clock product COD0MGXFIN (CODE), with WUM (Wuhan) and GRG (CNES) MGEX clock products as day-level fallbacks, downloaded via authenticated HTTPS. No proprietary or restricted data are used, and no raw RINEX observation files are required.

Analysis code and pipeline outputs are available at github.com/matthewsmawfield/TEP-GNSS-MGEX.

8. References

8.1 TEP Series

• Smawfield, M. L. (2025). Temporal Equivalence Principle: Dynamic Time & Emergent Light Speed. v0.9 (Jakarta). DOI: 10.5281/zenodo.16921911.
• Smawfield, M. L. (2025). Global Time Echoes: Distance-Structured Correlations in GNSS Clocks (Paper 1, TEP-GNSS). DOI: 10.5281/zenodo.17127229.
• Smawfield, M. L. (2025). Global Time Echoes: 25-Year Temporal Evolution (Paper 2, TEP-GNSS-II). DOI: 10.5281/zenodo.17517141.
• Smawfield, M. L. (2025). Global Time Echoes: Raw RINEX / SPP Consistency Test (Paper 3, TEP-GNSS-RINEX).

8.2 GNSS Data and Products

• Montenbruck, O., Steigenberger, P., & Hauschild, A. 2014, GPS Solutions, 18, 49, "Multi-signal GNSS data from the IGS MGEX experiment — status and outlook".
• International GNSS Service (IGS). MGEX Multi-GNSS Experiment. https://igs.org/mgex/.
• NASA CDDIS. Crustal Dynamics Data Information System. https://cddis.nasa.gov.
• Center for Orbit Determination in Europe (CODE). MGEX combined products. https://www.aiub.unibe.ch.

8.3 Cosmology and CMB

• Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, "Planck 2018 results. VI. Cosmological parameters", doi:10.1051/0004-6361/201833910.

8.4 Methodological References

• Welch, P. D. 1967, "The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms", IEEE Trans. Audio Electroacoust., AU-15, 70–73.

Smawfield, M. L. 2026. Temporal Equivalence Principle series, Papers 0–25. Zenodo preprints and associated repositories.

9. Reproducibility Statement

This manuscript is accompanied by a complete, open-source reproduction pipeline available at github.com/matthewsmawfield/TEP-GNSS-MGEX. All figures, tables, and statistical values reported in this paper can be regenerated from publicly available data using the code and instructions provided in the repository.

Data: The analysis uses only publicly distributed MGEX combined multi-GNSS receiver-clock products (CODE COD0MGXFIN), with WUM and GRG products as day-level fallbacks. Data were downloaded from NASA CDDIS (cddis.nasa.gov) and the IGS (igs.org) via standard authenticated HTTPS endpoints. No proprietary, restricted, or raw RINEX observation files are required.

Software: The pipeline is written in Python (NumPy, SciPy, Astropy) and standard geodetic libraries. All dependencies are listed in requirements.txt in the repository root. No commercial or closed-source software is required.

Execution: Running python run_analysis.py from the repository root, after installing dependencies, reproduces the full analysis chain: data download, station-pair differencing, phase-coherent stacking, isotropic and anisotropic correlation fitting, orbital-velocity projection, and figure generation. Expected wall-clock time on a modern laptop is approximately 30 minutes for the full 2025-01-01 to 2026-05-01 dataset (256 stations).

Verification: The repository includes a CHECKSUMS.md file with SHA-256 hashes for all derived output files, permitting direct bit-level verification of reproduced results.

Environment: The pipeline has been tested on macOS 14 and Ubuntu 22.04 with Python 3.11. A Dockerfile is provided for fully containerised reproduction if environment consistency is required.