A Framework for Earth–Moon Reference-Clock Modeling, Compensation, and Timing-Requirement Assessment in Cislunar PNT

1. Introduction

Renewed scientific and engineering activity in the Earth–Moon system has increased interest in positioning, navigation, and timing (PNT) architectures that can operate beyond the immediate vicinity of Earth. Recent studies of lunar and cislunar navigation have considered architectures that combine Earth-based infrastructure, lunar-orbiting sources, surface assets, and Earth–Moon or inter-satellite communication links [1,2,3,4]. The definition of a lunar reference frame and the realization of a corresponding lunar time scale are also recognized as essential elements of future lunar navigation architectures [5,6]. Within such architectures, time is required not only for clock synchronization, but also for signal transmission, observation time tagging, orbit determination, and pseudorange formation. The relationship between terrestrial and lunar time references is therefore relevant both to the physical definition of reference times and to their implementation in navigation measurements [5,7,8,9,10].

Clocks associated with Earth-based and Moon-based reference systems do not accumulate proper time at identical rates. In the weak-field and slow-motion approximation, the rate of a clock depends primarily on the gravitational potential and the clock’s kinematic state [7,9,11]. Consequently, clocks located in the terrestrial and lunar environments experience different mean rates and time-dependent variations. The dominant contribution is associated with the different gravitational environments of the Earth and Moon, whereas external-body potentials, orbital motion, rotational motion, and higher-order gravitational structure introduce additional corrections [12,13,14]. Even a small difference in fractional frequency can accumulate into a measurable time offset over sufficiently long intervals.

The relativistic treatment of these effects requires a consistent definition of reference systems and coordinate time scales. The International Astronomical Union framework provides the basis for barycentric and geocentric relativistic reference systems and their associated time transformations [7,8,9]. The practical realization of these systems, including the treatment of terrestrial time scales and Earth-related reference quantities, is further specified by the IERS conventions [10]. These foundations have been extended in recent lunar studies to Earth–Moon time transformations, lunar reference systems and frames, and possible realizations of lunar coordinate time [5,6,15].

Recent work has also investigated the numerical realization of lunar time and the construction of lunar time products. Turyshev et al. derived relativistic transformations among barycentric, geocentric, and lunicentric time coordinates and quantified both secular and periodic components of the terrestrial–lunar time difference [15]. Fienga et al. discussed lunar reference systems, lunar reference frames, and possible implementations of a lunicentric time scale for orbiting and surface-based clocks [5]. In addition, lunar time-ephemeris studies have considered the generation and dissemination of numerical correction products for practical applications [16,17]. These studies demonstrate that an operational treatment of lunar timing must distinguish the underlying physical relation between reference times from the residual error that remains after the relation has been modeled, predicted, or disseminated as a correction.

In parallel, relativistic clock modeling and time-frequency transfer have been extensively studied in the terrestrial and solar-system contexts. Relativistic corrections are required for high-accuracy time transfer because gravitational and kinematic effects influence both clock rates and the propagation of timing information [12,13,14]. Satellite-navigation analyses generally represent receiver-clock errors as range-equivalent terms and estimate receiver clock bias jointly with the user position [18,19,20,21]. These results provide the theoretical, numerical, and estimation foundations needed for a cislunar timing analysis. However, the connection between a modeled Earth–Moon clock-rate difference and a system-level PNT timing requirement is not obtained by a single universal conversion.

In particular, multiplying a clock offset by the speed of light gives a one-way range-equivalent quantity,

$$\Delta {\rho }_{\mathrm{eq}}=c\Delta t,$$

(1)

but this quantity should not be interpreted directly as a three-dimensional position error. Consider a set of simultaneous source observations whose timing residual can be written as

$$\delta t_i=\delta t_{\mathrm{com}}+\delta t_{\mathrm{sd},i},$$

(2)

where $\delta t_{\mathrm{com}}$ is common to all sources and $\delta t_{\mathrm{sd},i}$ is the source-dependent component. The common component has, to first order, the same measurement signature as receiver clock bias and can therefore be absorbed predominantly by the estimated clock state. The source-dependent component, in contrast, projects into the position-sensitive observation subspace. Its contribution to the position solution depends on source geometry, observation weighting, estimator conditioning, and the remaining non-timing error budget [18,19,21].

The treatment of compensation introduces a second distinction. A model estimated from past observations and extrapolated into future epochs is a forecast, whereas a tabulated product containing correction nodes over the evaluation interval provides interpolation or reconstruction. These two cases rely on different information conditions and should not be evaluated as though they were equivalent. A continuously available correction product may produce a small nominal interpolation residual, while its performance can degrade during delivery outages if local propagation does not adequately represent the time-dependent clock rate. Consequently, chronological validation, correction cadence, update continuity, and outage behavior must be evaluated separately. This distinction is particularly important when a numerical correction product is generated from the same modeled series that is subsequently used for validation.

A third difficulty concerns the derivation of timing requirements from positioning targets. Timing accuracy is only one component of the complete observation error budget. Orbit uncertainty, source-clock noise, propagation delay errors, measurement noise, hardware delays, and estimator geometry may impose a position-error floor. Once the source-dependent timing contribution is reduced below these other components, the total position error approaches a geometry-dependent limit. Some positioning targets may therefore be unattainable under a given geometry even if the source-dependent timing error is reduced arbitrarily. A timing requirement obtained without checking this floor can assign an unnecessarily strict or physically irrelevant specification to the clock subsystem [19,20,21].

Accurate numerical evaluation of the Earth–Moon clock-rate difference also depends on the adopted planetary and lunar ephemerides and lunar gravity field. The DE430/DE431 and DE440/DE441 ephemerides provide successive high-precision frameworks for planetary and lunar state-vector calculations [22,23]. Likewise, lunar gravity models derived from Clementine and GRAIL tracking data provide the basis for representing the lunar gravitational environment at increasing spherical-harmonic resolution [24,25]. The present study uses DE440s barycentric states as the input for an idealized reference-clock calculation. The selected ephemeris therefore serves as the dynamical input for the numerical experiment, rather than as an attempt to establish a new lunar ephemeris or a complete lunar gravity-field realization.

This study examines the above connections through three linked numerical experiments. The first constructs an idealized Earth–Moon reference-clock series from DE440s barycentric states using a first-order weak-field proper-time model. The calculation separates the contributions from self-gravity, external gravity, orbital motion, and equatorial rotation, and quantifies both the secular accumulation and the remaining time-dependent variation. The resulting series is used as a controlled physical input for the subsequent compensation and PNT experiments.

The second experiment evaluates a sequence of compensation strategies. These include no compensation, a training-only linear-drift model, prescribed harmonic forecast models, and a full-interval tabulated correction product. Chronological training and validation intervals are used to distinguish out-of-sample forecasting from in-sample fitting. Multi-window validation and harmonic ablation are used to assess the stability of the parametric models, while cadence and update-outage tests characterize the behavior of the idealized correction product under different delivery conditions.

The third experiment propagates the resulting compensation residuals through a linearized four-state PNT estimator that jointly solves for three-dimensional position and receiver clock range bias. The residual is decomposed into common and source-dependent components so that their different effects on the estimator can be examined. Controlled source geometries and a parameterized non-timing error budget are then used to evaluate position RMSE and to invert the maximum allowable source-dependent timing error for selected positioning targets. The geometry-sensitivity analysis is motivated by prior lunar PNT and lunar-augmented navigation studies [1,2,3]. The configurations used here, however, are deliberately controlled cases rather than proposed operational constellations.

The analysis addresses three questions:

1.

What secular and time-dependent structures arise in the adopted first-order model of the Earth–Moon reference-clock difference?

2.

How do training-only forecast models and a continuously updated correction product differ in residual accuracy, robustness, cadence sensitivity, and outage behavior?

3.

How do common and source-dependent timing residuals map into position error, and how do source geometry and non-timing errors constrain the timing accuracy associated with a specified positioning target?

The resulting framework links the physical clock model, compensation residual, observation-error decomposition, and conditional timing allocation within a consistent sequence. Its purpose is not to define a formal lunar coordinate time, to reproduce a complete TCG–TCL or TT–TL transformation, or to predict the performance of a particular flight constellation. Rather, it provides a controlled numerical analysis of how assumptions made at each stage affect the subsequent timing and PNT results.

The primary contributions of this work are threefold:We establish an architecture-level framework that functionally maps first-order Earth-Moon weak-field clock differentials to explicit deep-space PNT budgets; Through rigorous harmonic ablation and product-outage simulations, we expose the architectural dichotomy between the stability of autonomous analytical forecasting and the precision of high-cadence correction architectures; We mathematically demystify the scalar range-equivalent metric under practical lunar geometries, mandating a coupled geometric-stochastic paradigm that redefines how timing allocations must be constrained strictly by the non-timing error floor.

These contributions are intended as a methodological and architecture-sensitivity study. The reference clocks are idealized, and the relativistic model retains leading weak-field and slow-motion terms rather than implementing a formal lunar coordinate-time realization. Non-spherical lunar gravity, tidal deformation, detailed body orientation, site-specific clock locations, orbit–rotation cross terms, higher-order relativistic terms, and mission-specific time-transfer procedures are outside the present scope. In addition, the PNT geometries are controlled configurations rather than candidate flight constellations, and the non-timing error budget is represented parametrically. The numerical values should therefore be interpreted together with the assumptions under which they are obtained. A mission-specific assessment would require time-varying source trajectories and visibility, explicit time-transfer links, correlated error models, realistic observation scheduling, and a nonlinear navigation filter.

2. Theoretical Framework

2.1. Relativistic Time-Reference Difference

For a clock moving with coordinate velocity $\mathbf{v}$ in a weak gravitational field, the relation between proper time $\tau$ and a chosen coordinate time t can be written, to first order in the weak-field and slow-motion approximation, as

$${\displaystyle \frac{\mathrm{d}\tau }{\mathrm{d}t}}=1-{\displaystyle \frac{1}{c^2}}\left(U+{\displaystyle \frac{v^2}{2}}\right)+\mathcal{O}\left(c^{-4}\right),$$

(3)

where U denotes the adopted positive gravitational potential, v is the coordinate speed of the clock, and c is the speed of light.

For two reference clocks, denoted by E and L, the differential rate with respect to the same coordinate time is therefore

$${\displaystyle \frac{\mathrm{d}\Delta t}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}({\tau }_L-{\tau }_E)}{\mathrm{d}t}}={\displaystyle \frac{1}{c^2}}\left[U_E-U_L+{\displaystyle \frac{v_E^2-v_L^2}{2}}\right]+\mathcal{O}\left(c^{-4}\right).$$

(4)

The potential and velocity terms can be separated according to their physical origin:

$${\displaystyle \frac{\mathrm{d}\Delta t}{\mathrm{d}t}}={\dot{\Delta t}}_{\mathrm{grav}}+{\dot{\Delta t}}_{\mathrm{kin}},$$

(5)

with

$$\begin{array}{cc}\hfill {\dot{\Delta t}}_{\mathrm{grav}}& ={\displaystyle \frac{U_E-U_L}{c^2}},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\dot{\Delta t}}_{\mathrm{kin}}& ={\displaystyle \frac{v_E^2-v_L^2}{2c^2}}.\hfill \end{array}$$

(7)

In the cislunar case, the gravitational contribution may contain the self-potentials of the Earth and Moon as well as the potentials generated by external bodies. The kinematic contribution may include the orbital and rotational motion terms used in the adopted reference model. The decomposition used in the numerical experiments is given in Section 3.1.

The accumulated time-reference difference is obtained by integration:

$$\Delta t\left(t\right)=\Delta t\left(t_0\right)+{\int }_{t_0}^t{\displaystyle \frac{\mathrm{d}\Delta t\left(\xi \right)}{\mathrm{d}\xi }}\mathrm{d}\xi .$$

(8)

Equation (8) shows why a small persistent rate difference can result in a large accumulated time difference over a long interval. It also explains why a compensation model must represent both the secular component and the time-dependent variations of the rate difference.

The time difference can be expressed as a one-way range-equivalent quantity:

$$\Delta {\rho }_{\mathrm{eq}}\left(t\right)=c\Delta t\left(t\right).$$

(9)

This conversion is used to compare the magnitude of the timing effect with pseudorange errors. It is not interpreted as a direct three-dimensional position error because the receiver clock is estimated jointly with the position states.

2.2. Observation Model and Error Decomposition

For $N_s$ observed signal sources, the linearized pseudorange observation model is written as

$$\delta \rho =\mathbf{H}\delta \mathbf{x}+\epsilon ,$$

(10)

where

$$\delta \mathbf{x}={\left[\begin{array}{cccc}\delta r_x& \delta r_y& \delta r_z& \delta b\end{array}\right]}^{\mathsf{T}}$$

(11)

contains the three position errors and the receiver-clock range-equivalent error. The design matrix is

$$\mathbf{H}=\left[\begin{array}{cc}-{\mathbf{u}}_1^{\mathsf{T}}& 1\\ -{\mathbf{u}}_2^{\mathsf{T}}& 1\\ ⋮& ⋮\\ -{\mathbf{u}}_{N_s}^{\mathsf{T}}& 1\end{array}\right],$$

(12)

where ${\mathbf{u}}_i$ is the unit line-of-sight vector from the receiver to source i.

The observation error is decomposed into a timing component and a non-timing component:

$$\epsilon ={\epsilon }_{\mathrm{time}}+{\epsilon }_{\mathrm{non}-\mathrm{time}}.$$

(13)

The timing component is further separated into common and differential parts:

$${\epsilon }_{\mathrm{time}}=c{r}_{\mathrm{com}}\mathbf{1}+c{\mathbf{r}}_{\mathrm{diff}},$$

(14)

where $\mathbf{1}$ is an $N_s$-dimensional vector of ones.

The common component has the same value in every observation. Its observation signature is parallel to the receiver-clock column of $\mathbf{H}$ and can therefore be estimated predominantly through $\delta b$. The differential component has source-dependent values and is not generally parallel to that column. It is consequently capable of producing a position error.

For a weighted least-squares solution, let

$$\mathbf{K}={\left({\mathbf{H}}^{\mathsf{T}}\mathbf{W}\mathbf{H}\right)}^{-1}{\mathbf{H}}^{\mathsf{T}}\mathbf{W},$$

(15)

where $\mathbf{W}$ is the observation-weight matrix. The estimated state error caused by an observation perturbation is

$$\delta \widehat{\mathbf{x}}=\mathbf{K}\epsilon .$$

(16)

Let ${\mathbf{S}}_r$ extract the three position components from the state. The corresponding position error is

$$\delta \widehat{\mathbf{r}}={\mathbf{S}}_r\mathbf{K}\epsilon .$$

(17)

For the common timing term,

$${\epsilon }_{\mathrm{com}}=c{r}_{\mathrm{com}}\mathbf{1},$$

(18)

the resulting position contribution is ideally zero or small when the receiver-clock state is correctly included and the observation model is consistent. For the differential term,

$${\epsilon }_{\mathrm{diff}}=c{\mathbf{r}}_{\mathrm{diff}},$$

(19)

the position contribution is generally non-zero and depends on $\mathbf{H}$.

This provides the theoretical basis for separating the total compensation residual from the source-dependent timing residual in Experiment 03.

2.3. Compensation Residuals and Temporal Prediction

Let $\Delta t\left(t\right)$ denote the reference time-difference series and ${\widehat{\Delta t}}_m\left(t\right)$ its estimate under compensation model m. The residual is

$$r_m\left(t\right)=\Delta t\left(t\right)-{\widehat{\Delta t}}_m\left(t\right).$$

(20)

For a parametric forecast model, the estimate is determined by a training data set:

$${\widehat{\Delta t}}_m\left(t\right)={\mathcal{M}}_m\left[\Delta t\left(\xi \right);\xi \in {\mathcal{T}}_{\mathrm{train}}\right],t\in {\mathcal{T}}_{\mathrm{valid}}.$$

(21)

The validation residual therefore measures out-of-sample prediction performance. The model complexity determines which temporal components can be represented, but increasing complexity does not guarantee improved validation performance when the fitted components are not stable beyond the training interval.

For a full-interval correction product, the estimate is instead obtained from correction nodes available throughout the evaluation period:

$${\widehat{\Delta t}}_{\mathrm{prod}}\left(t\right)=\mathcal{I}\left[\Delta t\left(t_k\right);t_k\in {\mathcal{T}}_{\mathrm{full}}\right],$$

(22)

where $\mathcal{I}$ denotes the selected interpolation or reconstruction operator.

Equations (21) and (22) represent different information conditions. The former describes prediction without future values; the latter describes reconstruction while the correction product remains available. The two cases are therefore reported separately in Experiment 02.

2.4. PNT Error Propagation and Timing Requirement Inversion

Let the differential timing-error covariance be represented as

$${\mathbf{C}}_{t,\mathrm{diff}}={\sigma }_t^2{\mathbf{R}}_{\mathrm{diff}},$$

(23)

where ${\sigma }_t$ is the characteristic source-dependent timing error and ${\mathbf{R}}_{\mathrm{diff}}$ describes its normalized spatial correlation or source pattern.

The corresponding pseudorange-error covariance is

$${\mathbf{C}}_{\rho ,\mathrm{time}}=c^2{\mathbf{C}}_{t,\mathrm{diff}}.$$

(24)

Let ${\mathbf{C}}_{\rho ,\mathrm{non}-\mathrm{time}}$ denote the covariance of the remaining error budget. Under the independent-error assumption used in the primary analysis,

$${\mathbf{C}}_{\rho }=c^2{\sigma }_t^2{\mathbf{R}}_{\mathrm{diff}}+{\mathbf{C}}_{\rho ,\mathrm{non}-\mathrm{time}}.$$

(25)

The state-error covariance is then

$${\mathbf{C}}_{\widehat{x}}=\mathbf{K}{\mathbf{C}}_{\rho }{\mathbf{K}}^{\mathsf{T}}.$$

(26)

The three-dimensional position RMS is defined by

$${\sigma }_{\mathrm{pos}}={\left[tr\left({\mathbf{S}}_r{\mathbf{C}}_{\widehat{x}}{\mathbf{S}}_r^{\mathsf{T}}\right)\right]}^{1/2}.$$

(27)

Substitution of Eq. (25) into Eq. (26) gives

$${\sigma }_{\mathrm{pos}}^2=c^2{\sigma }_t^2{\kappa }_t+{\sigma }_{\mathrm{floor}}^2,$$

(28)

where

$$\begin{array}{cc}\hfill {\kappa }_t& =tr\left[{\mathbf{S}}_r\mathbf{K}{\mathbf{R}}_{\mathrm{diff}}{\mathbf{K}}^{\mathsf{T}}{\mathbf{S}}_r^{\mathsf{T}}\right],\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{floor}}^2& =tr\left[{\mathbf{S}}_r\mathbf{K}{\mathbf{C}}_{\rho ,\mathrm{non}-\mathrm{time}}{\mathbf{K}}^{\mathsf{T}}{\mathbf{S}}_r^{\mathsf{T}}\right].\hfill \end{array}$$

(30)

The factor ${\kappa }_t$ contains the influence of source geometry and differential-error structure. The quantity ${\sigma }_{\mathrm{floor}}$ is the position-error floor caused by non-timing errors.

For a prescribed position target $P_{\mathrm{target}}$, the allowable source-dependent timing error follows from

$${\sigma }_{t,max}={\left[{\displaystyle \frac{P_{\mathrm{target}}^2-{\sigma }_{\mathrm{floor}}^2}{c^2{\kappa }_t}}\right]}^{1/2},$$

(31)

provided that

$$P_{\mathrm{target}}>{\sigma }_{\mathrm{floor}}.$$

(32)

If Eq. (32) is not satisfied, the target is not attainable by improving timing accuracy alone. In the numerical implementation, the requirement is obtained through the forward Monte Carlo model and numerical threshold search described in Section 3.11. Equation (31) is retained as the analytical interpretation of that numerical inversion.

The theoretical result has two limiting behaviors. When $c^2{\sigma }_t^2{\kappa }_t$ dominates ${\sigma }_{\mathrm{floor}}^2$, the positioning error is timing limited and approximately proportional to ${\sigma }_t$. When the timing contribution is small, the position error approaches ${\sigma }_{\mathrm{floor}}$. This transition explains why stricter time-reference accuracy does not necessarily improve the positioning result once the non-timing error floor has become dominant.

2.5. Scope of the Framework

The theoretical framework is used for an architecture-level sensitivity analysis. The time-reference model retains the leading weak-field and slow-motion terms and is not intended to constitute a complete operational realization of a lunar coordinate time. Higher-order relativistic terms, non-spherical gravity, tidal effects, detailed station motion, and mission-specific time-transfer procedures are outside the present scope.

The observation model is a linearized four-state position–clock model. The source geometries are prescribed controlled configurations rather than time-varying geometries from a specific constellation. The source-dependent timing structure is represented through a parameterized covariance or source pattern, and is therefore used to establish sensitivity rather than to claim a universal empirical value.

The non-timing error budget is represented by the component model defined in Section 3.6. Correlations among orbit, propagation, hardware, and clock errors are not included in the primary analysis. These assumptions are appropriate for isolating the effect of timing compensation and geometry. A mission-specific assessment would require time-varying ephemerides, source visibility, correlated error models, realistic observation scheduling, and a nonlinear navigation filter.

3. Methodology

The numerical analysis consisted of three consecutive experiments. Experiment 01 generated the reference Earth–Moon clock-difference series. Experiment 02 evaluated the compensation models and correction-product outages. Experiment 03 propagated the validated compensation residuals through a four-state PNT solution and inverted the corresponding timing requirements. The output of each experiment was used directly as the input to the next experiment, thereby preserving a consistent data lineage throughout the analysis.

3.1. Generation of the Reference Clock Series

The relativistic rate relation introduced in Section 2.1 was evaluated over a simulation interval of

$$0\le t\le 365.25\mathrm{day}.$$

(33)

The calculation used a uniform sampling interval of

$$\Delta t_{\mathrm{samp}}=0.25\mathrm{day},$$

(34)

giving 1462 epochs including the initial and final epochs.

The positions and velocities required by the reference model were obtained from the DE440s planetary ephemeris. All contributions were evaluated at the same epochs and expressed with respect to the common coordinate-time argument used by the numerical model.

The rate difference was assembled from the gravitational and kinematic terms retained in the idealized model:

$$\dot{\Delta t}\left(t\right)={\dot{\Delta t}}_{\mathrm{self}}\left(t\right)+{\dot{\Delta t}}_{\mathrm{ext}}\left(t\right)+{\dot{\Delta t}}_{\mathrm{orb}}\left(t\right)+{\dot{\Delta t}}_{\mathrm{rot}}\left(t\right),$$

(35)

where the four terms denote the contributions from the adopted Earth–Moon self-potentials, external gravitational potentials, orbital motion, and rotational motion, respectively. The positive-potential convention of Eq. (3) was used throughout.

The accumulated time-reference difference was evaluated numerically from the sampled rate series:

$$\Delta t_j=\Delta t_{j-1}+{\displaystyle \frac{{\dot{\Delta t}}_{j-1}+{\dot{\Delta t}}_j}{2}}\left(t_j-t_{j-1}\right),$$

(36)

with

$$\Delta t_0=0.$$

(37)

Equation (36) applies the trapezoidal rule to the reference-rate sequence. The accumulated difference therefore represents the relative evolution from the initial epoch rather than an absolute synchronization offset.

The resulting data product contained, at each epoch, the elapsed time, total fractional rate difference, accumulated time difference, and the individual physical contributions in Eq. (35). This series was stored as the reference input for Experiments 02 and 03. No compensation model was used during its generation.

To examine the temporal structure independently of the secular drift, a linear trend was fitted to the accumulated series:

$$\Delta t_{\mathrm{lin}}\left(t\right)=a_0+a_1(t-t_{\mathrm{ref}}),$$

(38)

and the detrended component was calculated as

$$\Delta t_{\det }\left(t\right)=\Delta t\left(t\right)-\Delta t_{\mathrm{lin}}\left(t\right).$$

(39)

The fitted trend was used only to characterize the reference sequence in Experiment 01. Compensation performance was assessed separately using the training–validation procedure described below.

3.2. Training and Validation Design

The reference series was divided chronologically into a training interval and an independent validation interval:

$$\begin{array}{cc}\hfill {\mathcal{T}}_{\mathrm{train}}& =\left\{t:0\le t\le 180.0\mathrm{day}\right\},\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill {\mathcal{T}}_{\mathrm{valid}}& =\left\{t:180.0<t\le 365.25\mathrm{day}\right\}.\hfill \end{array}$$

(41)

Parameters of the forecast models M1–M3 were estimated exclusively from ${\mathcal{T}}_{\mathrm{train}}$. Their performance statistics were calculated over ${\mathcal{T}}_{\mathrm{valid}}$. This chronological separation was used instead of random cross-validation because the intended application is future extrapolation from previously available clock data.

For any compensation model m, the residual at epoch $t_j$ was defined as

$$r_{m,j}=\Delta t_j-{\widehat{\Delta t}}_{m,j}.$$

(42)

Unless otherwise stated, the reported M0–M3 residual statistics refer only to epochs in ${\mathcal{T}}_{\mathrm{valid}}$.

The time residual was converted to an equivalent one-way range residual using

$$r_{\rho ,m,j}=c{r}_{m,j}.$$

(43)

This quantity was used for scale comparison and as an input to the source-dependent PNT error model. It was not treated as a direct three-dimensional position error.

3.3. Compensation Models

Five compensation cases, denoted M0–M4, were evaluated. M0 provided the uncompensated reference case. M1–M3 were training-only forecast models of increasing temporal complexity. M4 represented a continuously available tabulated correction product.

3.3.1. M0: No Compensation

For the uncompensated case,

$${\widehat{\Delta t}}_0\left(t\right)=0,$$

(44)

and therefore

$$r_0\left(t\right)=\Delta t\left(t\right).$$

(45)

M0 retained the full accumulated reference-clock difference and provided the baseline for calculating relative improvement.

3.3.2. M1: Constant-Rate Compensation

M1 represented the reference series by an offset and a constant rate:

$${\widehat{\Delta t}}_1\left(t\right)={\beta }_{1,0}+{\beta }_{1,1}(t-t_{\mathrm{ref}}).$$

(46)

The coefficients were estimated by least squares over ${\mathcal{T}}_{\mathrm{train}}$. The reference epoch $t_{\mathrm{ref}}$ was fixed within the training interval to reduce correlation between the fitted intercept and slope. The fitted rate used in the experiment was

$${\widehat{\beta }}_{1,1}=56.426517818\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}.$$

(47)

The fitted model was extrapolated without re-estimation throughout ${\mathcal{T}}_{\mathrm{valid}}$.

3.3.3. M2: Linear and Single-Period Compensation

M2 augmented the linear model with one harmonic component:

$${\widehat{\Delta t}}_2\left(t\right)={\beta }_{2,0}+{\beta }_{2,1}(t-t_{\mathrm{ref}})+A_{2,1}cos\left({\omega }_{1}t\right)+B_{2,1}sin\left({\omega }_{1}t\right),$$

(48)

where ${\omega }_1=2\pi /P_1$ and $P_1$ was the adopted principal lunar period. Both sine and cosine terms were included so that the amplitude and phase were estimated from the training data rather than prescribed.

For a fixed period, Eq. (48) is linear in its coefficients. The parameter vector was therefore estimated by ordinary least squares over ${\mathcal{T}}_{\mathrm{train}}$ and applied unchanged to the validation interval.

3.3.4. M3: Linear and Multi-Period Compensation

M3 extended the harmonic representation to the set of astronomical periods retained in Experiment 02:

$${\widehat{\Delta t}}_3\left(t\right)={\beta }_{3,0}+{\beta }_{3,1}(t-t_{\mathrm{ref}})+{\displaystyle \sum _{k=1}^{N_p}}\left[A_{3,k}cos\left({\omega }_{k}t\right)+B_{3,k}sin\left({\omega }_{k}t\right)\right],$$

(49)

where ${\omega }_k=2\pi /P_k$. The periods $P_k$ were fixed before model fitting and were not selected from the validation residuals.

The coefficient vector was obtained from

$${\widehat{\theta }}_3=arg\underset{{\theta }_3}{min}\sum _{t_j\in {\mathcal{T}}_{\mathrm{train}}}{\left[\Delta t_j-{\mathbf{z}}_3^{\mathsf{T}}\left(t_j\right){\theta }_3\right]}^2,$$

(50)

where ${\mathbf{z}}_3\left(t\right)$ contains the intercept, linear term, and paired harmonic basis functions. The fitted model was then extrapolated over ${\mathcal{T}}_{\mathrm{valid}}$.

3.3.5. M4: Tabulated Correction Product

M4 was formed from clock-correction nodes spanning the complete simulation interval. The nominal node spacing was

$$\Delta T_{\mathrm{prod}}=1.0\mathrm{day}.$$

(51)

At an observation epoch between adjacent nodes $t_k$ and $t_{k+1}$, the correction was reconstructed by the interpolation operator used in Experiment 02:

$${\widehat{\Delta t}}_4\left(t\right)=\mathcal{I}\left[\Delta t\left(t_k\right),\Delta t\left(t_{k+1}\right);t\right],t_k\le t\le t_{k+1}.$$

(52)

Because the correction nodes were available within the evaluation interval, M4 was classified as product reconstruction rather than out-of-sample forecasting. Its residual characterized the representation and interpolation error of the idealized correction product. M4 and the reference series were generated from the same underlying model; the M4 result therefore did not include independent product-generation, time-transfer, or dissemination errors.

3.4. Compensation Performance Metrics

The validation residuals were summarized using the mean error, root-mean- square error, and maximum absolute error:

$$\begin{array}{cc}\hfill {\overline{r}}_m& ={\displaystyle \frac{1}{N_v}}{\displaystyle \sum _{j=1}^{N_v}}{r}_{m,j},\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill {\mathrm{RMS}}_{t,m}& ={\left({\displaystyle \frac{1}{N_v}}{\displaystyle \sum _{j=1}^{N_v}}{r}_{m,j}^2\right)}^{1/2},\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill r_{m,max}& =\underset{1\le j\le N_v}{max}\left|r_{m,j}\right|,\hfill \end{array}$$

(55)

where $N_v$ is the number of validation epochs.

The equivalent range-domain metrics were obtained using Eq. (43). Relative improvement over M0 was calculated as

$$I_m=\left(1-{\displaystyle \frac{{\mathrm{RMS}}_{t,m}}{{\mathrm{RMS}}_{t,0}}}\right)\times 100\%.$$

(56)

The same epochs and reference values were used for all models. This ensured that the reported differences arose from the compensation strategy rather than from unequal validation samples.

3.5. Robustness, Harmonic Ablation, and Product Outage

Additional tests were used to examine whether the compensation ranking was dependent on the nominal fitting configuration.

For the parametric models, the training boundary and retained harmonic terms were varied within the prescribed test configurations. Each modified model was fitted using only the corresponding training data and was evaluated over data that were excluded from fitting. Harmonic ablation was performed by removing selected periodic basis-function pairs from M3 and repeating the full estimation and validation procedure. The change in validation RMS was used to assess the influence of each retained period.

Correction-product continuity was examined separately from nominal M4 reconstruction. Product updates were interrupted at

$$t_{\mathrm{out}}=180.0\mathrm{day}.$$

(57)

At the outage epoch, the most recent available clock correction and the estimated pre-outage rate were retained. The rate used by the holdover model was

$${\widehat{\dot{\Delta t}}}_{\mathrm{out}}=75.693078733\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}.$$

(58)

During the interruption, the correction was propagated as

$${\widehat{\Delta t}}_{\mathrm{out}}\left(t\right)={\widehat{\Delta t}}_4\left(t_{\mathrm{out}}\right)+{\widehat{\dot{\Delta t}}}_{\mathrm{out}}(t-t_{\mathrm{out}}).$$

(59)

The outage residual was

$$r_{\mathrm{out}}\left(t\right)=\Delta t\left(t\right)-{\widehat{\Delta t}}_{\mathrm{out}}\left(t\right).$$

(60)

Outage durations of 1, 7, 30, 90, and 180 days were evaluated. For each duration, the maximum absolute timing residual and its one-way range-equivalent value were reported. This test represented loss of future correction nodes and was not included in the nominal M4 interpolation statistics.

3.6. Source Geometry and Four-State PNT Model

Experiment 03 used the M0–M4 residual series from the validation interval. The PNT solution estimated three receiver-position corrections and one receiver-clock correction:

$$\delta \mathbf{x}={\left[\begin{array}{cccc}\delta r_x& \delta r_y& \delta r_z& \delta b\end{array}\right]}^{\mathsf{T}}.$$

(61)

For source i, the linearized pseudorange error equation was

$$\delta {\rho }_i=-{\mathbf{u}}_i^{\mathsf{T}}\delta \mathbf{r}+\delta b+{\epsilon }_i,$$

(62)

where ${\mathbf{u}}_i$ is the unit line-of-sight vector and ${\epsilon }_i$ contains the injected timing and non-timing errors.

The observations were assembled as

$$\delta \rho =\mathbf{H}\delta \mathbf{x}+\epsilon ,$$

(63)

using the design matrix defined in Eq. (12). Each tested geometry was required to satisfy

$$rank\left(\mathbf{H}\right)=4.$$

(64)

Four controlled geometry classes were evaluated: GOOD, MODERATE, DEGRADED, and NEAR_SINGULAR. The first three were used to study increasing geometric amplification. NEAR_SINGULAR was retained as a numerical stress test and was not interpreted as a representative operational constellation.

The condition number of each design matrix was calculated as

$$\kappa \left(\mathbf{H}\right)={\displaystyle \frac{{\sigma }_{max}\left(\mathbf{H}\right)}{{\sigma }_{min}\left(\mathbf{H}\right)}},$$

(65)

where ${\sigma }_{max}$ and ${\sigma }_{min}$ are the largest and smallest singular values. Rank and condition-number checks were performed before the Monte Carlo propagation.

3.7. Common and Differential Timing Injection

A compensation residual that is identical in all observations has the same signature as the receiver-clock column of the design matrix. Experiment 03 therefore propagated only the prescribed source-dependent fraction of the M0–M4 residual into the position-sensitive observation subspace.

For compensation model m at validation epoch $t_j$, the raw source-dependent timing pattern was written as

$${\tilde{\mathbf{r}}}_{m,j}=f_{\mathrm{diff}}{r}_{m,j}\mathbf{q},$$

(66)

where $f_{\mathrm{diff}}$ is the differential mapping fraction and $\mathbf{q}$ is a dimensionless source pattern.

To exclude a residual common mode, the source pattern was projected onto the differential subspace. For the weight matrix $\mathbf{W}$, this gives

$${\mathbf{q}}_{\mathrm{diff}}=\left[\mathbf{I}-\mathbf{1}{\left({\mathbf{1}}^{\mathsf{T}}\mathbf{W}\mathbf{1}\right)}^{-1}{\mathbf{1}}^{\mathsf{T}}\mathbf{W}\right]\mathbf{q}.$$

(67)

The resulting pattern satisfies

$${\mathbf{1}}^{\mathsf{T}}\mathbf{W}{\mathbf{q}}_{\mathrm{diff}}=0.$$

(68)

After common-mode removal, the pattern was normalized according to the convention used in Experiment 03:

$${\displaystyle \frac{1}{N_s}}{\mathbf{q}}_{\mathrm{diff}}^{\mathsf{T}}{\mathbf{q}}_{\mathrm{diff}}=1.$$

(69)

The injected differential range error was then

$${\epsilon }_{m,j,\mathrm{time}}=c{f}_{\mathrm{diff}}{r}_{m,j}{\mathbf{q}}_{\mathrm{diff}}.$$

(70)

The tested range included $f_{\mathrm{diff}}=0$, representing complete common-mode absorption, and $f_{\mathrm{diff}}=0.05$, representing a case in which 5% of the compensation residual was mapped into source-dependent errors. Intermediate fractions were retained where specified by the multifactor experiment.

The mapping in Eq. (70) is an architecture-level sensitivity model. It does not assign the differential component to a particular source-clock, signal-generation, or hardware mechanism.

3.8. Non-Timing Error Budget

The non-timing pseudorange error was decomposed into ranging noise, source-clock error, orbit or ephemeris line-of-sight error, propagation residual, and hardware delay:

$${\epsilon }_{\mathrm{nt}}={\epsilon }_{\mathrm{rng}}+{\epsilon }_{\mathrm{clk}}+{\epsilon }_{\mathrm{orb}}+{\epsilon }_{\mathrm{prop}}+{\epsilon }_{\mathrm{hw}}.$$

(71)

Three error-budget levels were used to represent progressively less favorable measurement conditions. For budget level b, each component was generated from its assigned zero-mean distribution:

$${\epsilon }_k^{\left(b\right)}\sim \mathcal{N}\left(\mathbf{0},{\mathbf{C}}_k^{\left(b\right)}\right),k\in \{\mathrm{rng},\mathrm{clk},\mathrm{orb},\mathrm{prop},\mathrm{hw}\}.$$

(72)

Under the independence assumption of the primary experiment, the non-timing covariance was

$${\mathbf{C}}_{\mathrm{nt}}^{\left(b\right)}=\sum _k{\mathbf{C}}_k^{\left(b\right)}.$$

(73)

The numerical standard deviations assigned to the three budgets are listed in Table ??. Orbit and ephemeris uncertainty was represented as an effective line-of-sight range error rather than by perturbing the complete source orbit state.

For each realization, the total injected observation error was

$${\epsilon }_{m,j}^{\left(b\right)}={\epsilon }_{m,j,\mathrm{time}}+{\epsilon }_{\mathrm{nt}}^{\left(b\right)}.$$

(74)

3.9. PNT Error Propagation

For each combination of compensation model, validation epoch, geometry, differential fraction, and non-timing budget, the state error was calculated using weighted least squares:

$$\delta \widehat{\mathbf{x}}={\left({\mathbf{H}}^{\mathsf{T}}\mathbf{W}\mathbf{H}\right)}^{-1}{\mathbf{H}}^{\mathsf{T}}\mathbf{W}\epsilon .$$

(75)

For numerical evaluation, the corresponding system was solved using a rank-revealing matrix method rather than by explicitly forming a matrix inverse. Singular values and condition numbers were retained as solution diagnostics.

The three-dimensional position-error magnitude for realization ℓ was

$$e_{\mathrm{pos}}^{(\ell )}={∥{\mathbf{S}}_r\delta {\widehat{\mathbf{x}}}^{(\ell )}∥}_2,$$

(76)

where

$${\mathbf{S}}_r=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right].$$

(77)

The position RMS for a given test condition was calculated from the ensemble of Monte Carlo realizations:

$${\mathrm{RMS}}_{\mathrm{pos}}={\left[{\displaystyle \frac{1}{N_{\mathrm{MC}}}}{\displaystyle \sum _{\ell =1}^{N_{\mathrm{MC}}}}{\left(e_{\mathrm{pos}}^{(\ell )}\right)}^2\right]}^{1/2}.$$

(78)

The receiver-clock estimate was retained separately and was not included in the three-dimensional position metric. Common timing offsets were therefore not converted directly into position error.

3.10. Monte Carlo Design and Contribution Analysis

Monte Carlo propagation was used because the analysis combined time-varying compensation residuals, multiple error components, controlled geometries, and several differential mapping fractions. All compared compensation models used the same set of random non-timing error realizations for a given geometry and budget. This common-random-number design reduced sampling variation in model-to-model differences.

A fixed pseudorandom seed was used for reproducibility. Before aggregation, all required numerical outputs were checked for finite values. The input time series was required to be strictly monotonic, and only validation epochs were admitted to Experiment 03.

The contribution of error component k was evaluated by repeating the propagation with that component active and the remaining stochastic components suppressed. The contribution fraction was calculated from the component position-error variance:

$${\eta }_k={\displaystyle \frac{{\sigma }_{\mathrm{pos},k}^2}{{\sum }_j{\sigma }_{\mathrm{pos},j}^2}}\times 100\%.$$

(79)

Variance fractions, rather than RMS ratios, were used because independent error variances are additive under the linearized model. The contribution percentages were checked for closure after numerical rounding.

The complete multifactor design crossed the compensation mode, source geometry, differential mapping fraction, and non-timing budget. The same estimator and performance metric were used for all combinations.

3.11. Timing-Requirement Inversion

Timing requirements were obtained independently of the M0–M4 model labels by injecting a prescribed source-dependent timing standard deviation ${\sigma }_t$ into the same geometry and non-timing error model. For each tested value, the differential timing vector was generated as

$$\delta {\mathbf{t}}_{\mathrm{diff}}^{(\ell )}={\sigma }_t{\mathbf{q}}_{\mathrm{diff}}^{(\ell )},$$

(80)

and converted to the range domain:

$${\epsilon }_{\mathrm{time}}^{(\ell )}=c\delta {\mathbf{t}}_{\mathrm{diff}}^{(\ell )}.$$

(81)

The resulting position RMS was evaluated over the timing-error search grid for each geometry and non-timing budget. The positioning targets were

$$P_{\mathrm{target}}\in \left\{10\mathrm{m},100\mathrm{m},1\mathrm{k}\mathrm{m}\right\}.$$

(82)

The maximum allowable differential timing error was defined as

$${\sigma }_{t,max}=sup\left\{{\sigma }_t:{\mathrm{RMS}}_{\mathrm{pos}}\left({\sigma }_t\right)\le P_{\mathrm{target}}\right\}.$$

(83)

The numerical position-RMS curve was checked for effective monotonicity before the threshold was reported. Where the target crossing occurred between adjacent timing-grid values, the threshold was refined using the same forward model and acceptance criterion.

Attainability was tested first at

$${\sigma }_t=0.$$

(84)

If

$${\mathrm{RMS}}_{\mathrm{pos}}\left(0\right)>P_{\mathrm{target}},$$

(85)

the target was classified as unattainable for that geometry and non-timing budget. No timing threshold was assigned in such cases because the positioning target was already excluded by the non-timing error floor.

The numerically inverted thresholds were compared with the analytical behavior given by Eq. (31). The analytical expression was used to interpret the roles of geometry and the non-timing floor; the reported requirement values were taken from the numerical propagation.

3.12. Numerical Quality Control

Quality-control tests were applied at the input, model, and output levels. The Experiment 01 and Experiment 02 data versions and file checksums were verified before Experiment 03 was executed. The time axis was checked for monotonicity, and Experiment 03 was restricted to the declared validation interval.

All geometry matrices were required to have rank four. Operational geometries were checked against the adopted condition-number criterion, while the NEAR_SINGULAR case was required to trigger the intended stress condition. The zero and 5% differential-fraction cases were explicitly retained in the multifactor output.

The Monte Carlo implementation verified use of common random numbers across comparable cases. Required outputs were checked for finite values, error contribution fractions were checked for closure, and timing-requirement curves were checked for effective monotonicity. The output metadata also retained the interpretation that M1–M3 were training-only forecasts and M4 was correction-product reconstruction.

4. Results

The results are presented according to the information flow established in Section 2 and Section 3. Experiment 01 quantified the modeled Earth–Moon reference-clock difference. Experiment 02 evaluated the residual timing errors after parametric compensation and correction-product reconstruction. Experiment 03 mapped the source-dependent fraction of those residuals into a four-state PNT solution and inverted the timing accuracy required for specified positioning targets.

4.1. Earth–Moon Reference-Clock Difference

Figure 1 shows the instantaneous Earth–Moon clock-rate difference, the accumulated clock difference, and the residual remaining after removal of a linear trend. The accumulated series is dominated by a secular increase, whereas the instantaneous rate and detrended series contain substantial time-dependent variations.

Starting from the imposed initial condition $\Delta t\left(t_0\right)=0$, the accumulated difference at the end of the 365.25-day interval is

$$\Delta t_{\mathrm{final}}=20483.311346\mathrm{\mu }\mathrm{s}=20.483311\mathrm{m}\mathrm{s}.$$

(86)

The equivalent one-way range scale is approximately $6140.74$ $\mathrm{k}$$\mathrm{m}$. This value characterizes the magnitude of the uncompensated timing offset and is not interpreted as a direct three-dimensional position error.

The linear-fit rate is

$${\widehat{\dot{\Delta t}}}_{\mathrm{lin}}=56.017058\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}.$$

(87)

The arithmetic and time-weighted mean rates are $56.080533\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$  and $56.080250\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$, respectively. Their difference is only $0.000283\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$, indicating that the six-hour sampling interval does not materially bias the annual mean.

The instantaneous rate varies between $25.925973\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$ and $85.130911\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$. Its standard deviation is $19.204330\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$, and its peak-to-peak variation is $59.204938\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$. Thus, although the accumulated series is approximately linear at the annual scale, the instantaneous rate cannot be treated as constant at shorter timescales.

After the fitted linear trend is removed, the residual clock difference has an RMS of

$$\mathrm{RMS}\left(\Delta t_{\det }\right)=90.758824\mathrm{\mu }\mathrm{s},$$

(88)

and a maximum absolute value of $140.694069$ $\mathrm{\mu }$$\mathrm{s}$. These values establish the scale of the periodic and non-secular structure that remains after constant-rate compensation.

Table 1. Summary of the modeled Earth–Moon reference-clock difference.

Quantity	Value
Final accumulated difference	$20483.311346$$$$\mathrm{s}$
Linear-fit rate	$56.017058\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Arithmetic mean rate	$56.080533\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Time-weighted mean rate	$56.080250\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Rate standard deviation	$19.204330\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Minimum instantaneous rate	$25.925973\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Maximum instantaneous rate	$85.130911\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Detrended RMS	$90.758824$$\mathrm{\mu }$$\mathrm{s}$
Maximum absolute detrended residual	$140.694069$$\mathrm{\mu }$$\mathrm{s}$

Table 1. Summary of the modeled Earth–Moon reference-clock difference.

Quantity	Value
Final accumulated difference	$20483.311346$$$$\mathrm{s}$
Linear-fit rate	$56.017058\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Arithmetic mean rate	$56.080533\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Time-weighted mean rate	$56.080250\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Rate standard deviation	$19.204330\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Minimum instantaneous rate	$25.925973\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Maximum instantaneous rate	$85.130911\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$
Detrended RMS	$90.758824$$\mathrm{\mu }$$\mathrm{s}$
Maximum absolute detrended residual	$140.694069$$\mathrm{\mu }$$\mathrm{s}$

4.2. Physical Contributions to the Clock-Rate Difference

The mean contributions of the physical terms retained in the model are listed in Table 2. Their sum closes exactly to the arithmetic mean of the total rate.

The self-gravity term is the dominant positive contribution, with a mean rate of $57.365317\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$. External gravity and orbital motion reduce the rate by $0.980659\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$ and $0.408092\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$, respectively. The adopted equatorial rotation term adds $0.103967\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$.

The combined correction to the self-gravity contribution is therefore approximately

$$-0.980659-0.408092+0.103967=-1.284784\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}.$$

(89)

Although the self-gravity term determines the sign and principal magnitude of the mean rate, the external-potential and kinematic terms are required to recover the total result.

Numerical closure was achieved to $4.1\times {10}^{-10}\mathrm{\mu }\mathrm{s}$ for the accumulated decomposition and to $2.2\times {10}^{-11}\mathrm{\mu }\mathrm{s}$ between rate integration and the endpoint difference. The reported temporal structure is therefore not attributable to decomposition or numerical-integration inconsistency.

4.3. Primary Compensation Performance

Figure 2 compares the M0–M4 compensation cases. M0 retains the complete accumulated clock difference. M1–M3 are evaluated as out-of-sample forecasts after fitting to the first 180 days, whereas M4 is evaluated as full-interval reconstruction from a tabulated correction product.

The numerical validation results are given in Table 3. Without compensation, M0 has a validation RMS of $15542.356500$ $\mathrm{\mu }$$\mathrm{s}$ and a maximum absolute residual of $20483.311346$ $\mathrm{\mu }$$\mathrm{s}$. Its one-way range-equivalent RMS is approximately $4659.48$ $\mathrm{k}$$\mathrm{m}$.

M1 reduces the validation RMS to $94.105805$ $\mathrm{\mu }$$\mathrm{s}$, corresponding to an improvement of $99.3945\%$ relative to M0. This result is consistent with the dominant secular component identified in Experiment 01.

M2 does not improve on M1. Its validation RMS increases to $138.475356$ $\mathrm{\mu }$$\mathrm{s}$, and its maximum absolute residual reaches $235.476273$ $\mathrm{\mu }$$\mathrm{s}$. Thus, adding the sidereal-month harmonic alone degrades validation performance by approximately $47.1\%$ relative to the linear model:

$${\displaystyle \frac{138.475356-94.105805}{94.105805}}\times 100\%\approx 47.1\%.$$

(90)

The additional periodic term therefore improves neither the RMS nor the maximum-error metric under the tested extrapolation interval.

M3 produces the best result among the autonomous forecast models. Its validation RMS is $14.420851$ $\mathrm{\mu }$$\mathrm{s}$, which is approximately $6.52$ times smaller than that of M1 and $9.60$ times smaller than that of M2. Its maximum absolute residual is reduced to $44.400852$ $\mathrm{\mu }$$\mathrm{s}$. The result shows that the clock series contains predictable multi-period structure, but that its representation requires more than a single sidereal-month term.

M4 gives the smallest nominal residual, with an RMS of $0.379690$ $\mathrm{\mu }$$\mathrm{s}$ and a maximum absolute residual of $0.780873$ $\mathrm{\mu }$$\mathrm{s}$. The corresponding range-equivalent RMS is $113.828$ $\mathrm{m}$. M4 is approximately 38 times smaller in timing RMS than M3. This comparison, however, measures the advantage of continued access to a tabulated correction product and must not be interpreted as a forecast-model comparison.

4.4. Forecast Robustness and Harmonic-Model Selection

The multi-window and periodic-basis tests are shown in Figure 3. They reveal two distinct features: M3 is consistently more accurate than M1 and M2 across the tested training windows, but the full four-period basis is not the validation-optimal harmonic combination.

Across training endpoints from day 90 to day 240, evaluated over 90-day validation horizons, M1 has a median RMS of $92.829977$ $\mathrm{\mu }$$\mathrm{s}$, with a range from $90.988392$ $\mathrm{\mu }$$\mathrm{s}$ to $115.824906$ $\mathrm{\mu }$$\mathrm{s}$. M2 remains less accurate, with a median RMS of $132.047979$ $\mathrm{\mu }$$\mathrm{s}$ and a range from $124.081772$ $\mathrm{\mu }$$\mathrm{s}$ to $135.776456$ $\mathrm{\mu }$$\mathrm{s}$.

M3 has a median RMS of only $7.853737$ $\mathrm{\mu }$$\mathrm{s}$. Its minimum and maximum values are $6.526654$ $\mathrm{\mu }$$\mathrm{s}$ and $17.171302$ $\mathrm{\mu }$$\mathrm{s}$, respectively. Even its maximum multi-window RMS remains approximately $5.3$ times smaller than the minimum M1 RMS. The superior performance of a multi-period representation is therefore not confined to the primary 180-day training split.

Table 4. Multi-window validation RMS of M1–M3.

Model	Minimum	Median	Maximum
	Validation RMS ($\mathrm{\mu }\mathrm{s}$)
M1	90.988392	92.829977	115.824906
M2	124.081772	132.047979	135.776456
M3	6.526654	7.853737	17.171302

Table 4. Multi-window validation RMS of M1–M3.

Model	Minimum	Median	Maximum
	Validation RMS ($\mathrm{\mu }\mathrm{s}$)
M1	90.988392	92.829977	115.824906
M2	124.081772	132.047979	135.776456
M3	6.526654	7.853737	17.171302

The ablation results in Table 5 provide a more detailed view of model specification. The synodic-plus-annual model, A7, achieves the lowest validation RMS:

$${\mathrm{RMS}}_{A7}=4.629502\mathrm{\mu }\mathrm{s}.$$

(91)

This is $67.9\%$ smaller than the $14.420851$ $\mathrm{\mu }$$\mathrm{s}$ RMS of the full M3/A8 basis.

Table 5. Validation results of the periodic-basis ablation experiment, ordered by increasing RMS.

Case	Retained periodic terms	Validation RMS ($\mathbf{\mu }\mathbf{s}$)	Design condition number
A7	Synodic and annual	4.629502	46.748490
A8	Sidereal, synodic, anomalistic, and annual	14.420851	70.800911
A2	Synodic	15.787814	4.422460
A5	Sidereal, synodic, and anomalistic	19.517915	62.949530
A0	Linear only	94.105805	4.385159
A1	Sidereal	138.475356	4.395363
A3	Anomalistic	144.670843	4.392953
A4	Annual	428.913918	46.079852
A6	Sidereal and annual	617.768056	46.330170

Table 5. Validation results of the periodic-basis ablation experiment, ordered by increasing RMS.

Case	Retained periodic terms	Validation RMS ($\mathbf{\mu }\mathbf{s}$)	Design condition number
A7	Synodic and annual	4.629502	46.748490
A8	Sidereal, synodic, anomalistic, and annual	14.420851	70.800911
A2	Synodic	15.787814	4.422460
A5	Sidereal, synodic, and anomalistic	19.517915	62.949530
A0	Linear only	94.105805	4.385159
A1	Sidereal	138.475356	4.395363
A3	Anomalistic	144.670843	4.392953
A4	Annual	428.913918	46.079852
A6	Sidereal and annual	617.768056	46.330170

Neither the sidereal term alone nor the anomalistic term alone improves on the linear baseline. Their validation RMS values are $138.475356$ $\mathrm{\mu }$$\mathrm{s}$ and $144.670843$ $\mathrm{\mu }$$\mathrm{s}$, respectively. The annual term alone and the sidereal-plus-annual combination perform considerably worse, with RMS values of $428.913918$ $\mathrm{\mu }$$\mathrm{s}$ and $617.768056$ $\mathrm{\mu }$$\mathrm{s}$.

The results demonstrate that model complexity is not monotonically related to forecast accuracy. In particular, the full four-period basis is more accurate than M1 and M2 but is not the validation-optimal representation. The synodic-plus-annual basis provides both a lower validation RMS and a smaller condition number than the full model.

4.5. Correction-Product Cadence and Update Outages

The correction-product results separate two operational effects: interpolation error while product updates remain available and prediction error after product delivery is interrupted. These effects are shown in Figure 4.

4.5.1. Correction-Product Cadence

At the native sampling interval of 0.25 day, the reconstructed values coincide with the reference samples and the reported residual is zero. Increasing the product spacing to 0.5 day produces an RMS of $0.092154$ $\mathrm{\mu }$$\mathrm{s}$, corresponding to $27.627$ $\mathrm{m}$ in one-way range.

The primary one-day product has an RMS of $0.379690$ $\mathrm{\mu }$$\mathrm{s}$, or $113.828$ $\mathrm{m}$. Increasing the spacing to two and three days raises the RMS to $1.513525$ $\mathrm{\mu }$$\mathrm{s}$ and $3.381062$ $\mathrm{\mu }$$\mathrm{s}$, respectively. At a seven-day spacing, the RMS reaches $17.764247$ $\mathrm{\mu }$$\mathrm{s}$, equivalent to $5325.587$ $\mathrm{m}$.

Table 6. Sensitivity of M4 reconstruction to correction-product node spacing.

Node spacing	Timing RMS	Range-equivalent RMS
(day)	($\mathrm{\mu }\mathrm{s}$)	($\mathrm{m}$)
0.25	0.000000	0.000
0.50	0.092154	27.627
1.00	0.379690	113.828
2.00	1.513525	453.743
3.00	3.381062	1013.617
7.00	17.764247	5325.587

Table 6. Sensitivity of M4 reconstruction to correction-product node spacing.

Node spacing	Timing RMS	Range-equivalent RMS
(day)	($\mathrm{\mu }\mathrm{s}$)	($\mathrm{m}$)
0.25	0.000000	0.000
0.50	0.092154	27.627
1.00	0.379690	113.828
2.00	1.513525	453.743
3.00	3.381062	1013.617
7.00	17.764247	5325.587

The interpolation residual increases approximately quadratically over the shorter tested spacings. For example, doubling the interval from 0.5 to one day increases the RMS by a factor of $4.12$, while doubling it from one to two days increases the RMS by a factor of $3.99$. This behavior is consistent with interpolation error driven by curvature in the clock series.

4.5.2. Update Outages

Correction-product interruption produces much larger errors than nominal interpolation. After a one-day outage, the maximum residual is already $5.016297$ $\mathrm{\mu }$$\mathrm{s}$, equivalent to $1503.848$ $\mathrm{m}$. After seven days, it increases to $139.480622$ $\mathrm{\mu }$$\mathrm{s}$, or approximately $41.815$ $\mathrm{k}$$\mathrm{m}$.

The 30-day outage produces a maximum residual of $569.752519$ $\mathrm{\mu }$$\mathrm{s}$, corresponding to approximately $170.808$ $\mathrm{k}$$\mathrm{m}$. After 90 and 180 days, the maximum residuals reach $1607.030210$ $\mathrm{\mu }$$\mathrm{s}$ and $3178.082888$ $\mathrm{\mu }$$\mathrm{s}$, respectively.

Table 7. Maximum residual after interruption of the tabulated correction product.

Outage duration	Maximum residual	Range equivalent
(day)	($\mathrm{\mu }\mathrm{s}$)	($\mathrm{m}$)
1	5.016297	1503.848
7	139.480622	41815.239
30	569.752519	170807.508
90	1607.030210	481775.537
180	3178.082888	952765.281

Table 7. Maximum residual after interruption of the tabulated correction product.

Outage duration	Maximum residual	Range equivalent
(day)	($\mathrm{\mu }\mathrm{s}$)	($\mathrm{m}$)
1	5.016297	1503.848
7	139.480622	41815.239
30	569.752519	170807.508
90	1607.030210	481775.537
180	3178.082888	952765.281

The one-day outage residual is more than 13 times the nominal one-day product RMS. Thus, nominal M4 interpolation performance cannot be used to represent performance during product unavailability. Correction cadence and update continuity are separate system-design variables.

4.6. Common and Source-Dependent Timing Effects

The preliminary PNT scenarios in Experiment 02 verify the common–differential distinction derived in Section 2.2. When the clock state is not included, the common compensation residual maps directly into large position errors. In the S1 common-error scenario, the no-clock-state RMSE ranges from $23.350$ $\mathrm{m}$ for M4 to $863219.357$ $\mathrm{m}$ for M0.

When receiver position and receiver clock are estimated jointly, the S1 position RMSE is

$${\mathrm{RMSE}}_{\mathrm{S}1,\mathrm{joint}}=10.751993\mathrm{m}$$

(92)

for every compensation model. The invariance across M0–M4 confirms that a strictly common residual is absorbed by the receiver-clock state rather than being mapped directly into position.

Adding independent transmitting-node clock noise in S2 increases the joint position RMSE only slightly, from $10.751993$ $\mathrm{m}$ to $10.765998$ $\mathrm{m}$. In contrast, the S3 mixed scenario introduces a source-dependent component equal to $5\%$ of the model residual. The joint-state position RMSE then becomes strongly model dependent: $407829.441$ $\mathrm{m}$ for M0, $2479.769$ $\mathrm{m}$ for M1, $3623.321$ $\mathrm{m}$ for M2, $379.073$ $\mathrm{m}$ for M3, and $14.665$ $\mathrm{m}$ for M4.

Table 8. Joint position–clock RMSE in the Experiment 02 PNT scenarios.

Model	S1 common	S2 node	S3 mixed
	Position RMSE ($\mathrm{m}$)
M0	10.751993	10.765998	407829.440752
M1	10.751993	10.765998	2479.769422
M2	10.751993	10.765998	3623.321151
M3	10.751993	10.765998	379.073356
M4	10.751993	10.765998	14.665242

Table 8. Joint position–clock RMSE in the Experiment 02 PNT scenarios.

Model	S1 common	S2 node	S3 mixed
	Position RMSE ($\mathrm{m}$)
M0	10.751993	10.765998	407829.440752
M1	10.751993	10.765998	2479.769422
M2	10.751993	10.765998	3623.321151
M3	10.751993	10.765998	379.073356
M4	10.751993	10.765998	14.665242

These results demonstrate that the range-equivalent residual $c{r}_m$ cannot be interpreted directly as position error. Position degradation is controlled by the source-dependent residual component, its pattern across the observations, and the source geometry.

4.7. Geometry Diagnostics

Experiment 03 extends the preliminary PNT analysis to a complete parameterized error budget. All three controlled geometries have rank four, but their dilution and conditioning differ substantially, as shown in Table 9.

GOOD geometry has a PDOP of 1.099 and a normal-matrix condition number of 4.17. MODERATE geometry remains well posed but increases the PDOP to 1.994 and the condition number to 44.57. POOR geometry is formally rank four but has a PDOP of 538.15 and a condition number of $5.63\times {10}^6$. It is therefore interpreted as a deliberately stressed, near-degenerate configuration rather than an operational geometry.

4.8. Full Error-Budget Positioning Performance

Figure 5 shows the three-dimensional position errors obtained after combining the compensation residuals with the range-noise, node-clock, orbit line-of-sight, propagation, and hardware terms. The timing residual was mapped into source-dependent errors using $f_{\mathrm{diff}}=0.05$.

The numerical RMSE and 95th-percentile values are listed in Table 10. Under GOOD geometry, M0 produces an RMSE of $274148.514$ $\mathrm{m}$. M1 and M2 reduce this to $1659.950$ $\mathrm{m}$ and $2444.507$ $\mathrm{m}$, respectively. The larger M2 position error is consistent with its larger validation timing residual.

M3 reduces the GOOD-geometry RMSE to $269.094$ $\mathrm{m}$, with a 95th-percentile error of $640.656$ $\mathrm{m}$. M4 reduces the RMSE further to $9.711$ $\mathrm{m}$ and the 95th percentile to $16.493$ $\mathrm{m}$. Thus, under the assumed $5\%$ differential mapping and complete error budget, only M4 reaches an RMSE below 10 $\mathrm{m}$ in the GOOD geometry.

Under MODERATE geometry, the RMSE values increase to $476933.047$ $\mathrm{m}$, $2940.736$ $\mathrm{m}$, $4405.874$ $\mathrm{m}$, $436.185$ $\mathrm{m}$, and $16.645$ $\mathrm{m}$ for M0–M4, respectively. M4 remains the only model below 100 $\mathrm{m}$; however, it no longer reaches the 10 $\mathrm{m}$ target.

POOR geometry causes extreme amplification. Even with M4, the position RMSE is $4700.631$ $\mathrm{m}$, and the 95th percentile reaches $9786.030$ $\mathrm{m}$. The POOR-geometry result is approximately 484 times the GOOD-geometry M4 RMSE. Improving compensation alone is therefore insufficient when the source geometry approaches degeneracy.

The model ranking remains consistent across all geometries:

$$\mathrm{M}4<\mathrm{M}3<\mathrm{M}1<\mathrm{M}2<\mathrm{M}0,$$

(93)

where the ordering is from smaller to larger position RMSE. This ranking matches the primary validation timing RMS and confirms that the source-dependent timing component remains an important contributor under the adopted $5\%$ mapping.

4.9. Geometry Sensitivity

Figure 6 isolates the dependence of position performance on source geometry. The effect is not a constant multiplicative scale across all cases because the total position error contains both source-dependent timing errors and independent non-timing components.

Relative to GOOD geometry, MODERATE geometry increases the M4 RMSE from $9.711$ $\mathrm{m}$ to $16.645$ $\mathrm{m}$, a factor of approximately 1.71. The corresponding factors for M3 and M1 are approximately 1.62 and 1.77, respectively.

The transition from MODERATE to POOR geometry is much more severe. For M4, the RMSE increases from $16.645$ $\mathrm{m}$ to $4700.631$ $\mathrm{m}$, a factor of approximately 282. For M3 it increases from $436.185$ $\mathrm{m}$ to $118651.953$ $\mathrm{m}$, a factor of approximately 272. The consistency of these amplification factors confirms that the POOR configuration is dominated by geometry rather than by the choice of compensation model.

4.10. Timing-Requirement Inversion

Figure 7 shows the inverted relation between source-dependent timing accuracy and three-dimensional position RMSE. The maximum timing-error values satisfying the 10 $\mathrm{m}$, 100 $\mathrm{m}$, and 1 $\mathrm{k}$$\mathrm{m}$ targets are reported in Table 11.

Under GOOD geometry, the maximum allowable source-dependent timing errors are $21.325$ $\mathrm{n}$$\mathrm{s}$, $280.675$ $\mathrm{n}$$\mathrm{s}$, and $2813.290$ $\mathrm{n}$$\mathrm{s}$ for the 10 $\mathrm{m}$, 100 $\mathrm{m}$, and 1 $\mathrm{k}$$\mathrm{m}$ targets, respectively.

Under MODERATE geometry, the 10 $\mathrm{m}$ target is not achievable within the tested timing grid because the non-timing position floor already exhausts or exceeds the target budget. The allowable timing errors for the 100 $\mathrm{m}$ and 1 $\mathrm{k}$$\mathrm{m}$ targets decrease to $154.490$ $\mathrm{n}$$\mathrm{s}$ and $1555.708$ $\mathrm{n}$$\mathrm{s}$, respectively. These values are approximately $45\%$ smaller than the corresponding GOOD-geometry thresholds.

None of the three targets is achievable under POOR geometry within the tested timing range. In this case, the limiting factor is the geometric and non-timing error floor rather than timing accuracy. The result directly demonstrates the attainability condition in Eq. (32): an arbitrarily small timing error cannot overcome an unfavorable observation geometry.

4.11. Compatibility of Compensation Models with Positioning Targets

To compare the compensation models with the inverted timing requirements, the Experiment 03 mapping gives implied source-dependent timing standard deviations of $149668.780$ $\mathrm{n}$$\mathrm{s}$, $4567.695$ $\mathrm{n}$$\mathrm{s}$, $6714.422$ $\mathrm{n}$$\mathrm{s}$, $582.535$ $\mathrm{n}$$\mathrm{s}$, and $18.965$ $\mathrm{n}$$\mathrm{s}$ for M0–M4, respectively.

Table 12. Implied source-dependent timing accuracy of the compensation models and compatibility with the GOOD-geometry requirements.

Model	Implied ${\mathit{\sigma }}_{\mathit{t}}$	10$\mathbf{m}$	100$\mathbf{m}$	1$\mathbf{k}$$\mathbf{m}$
	($\mathrm{n}\mathrm{s}$)
M0	149668.780	No	No	No
M1	4567.695	No	No	No
M2	6714.422	No	No	No
M3	582.535	No	No	Yes
M4	18.965	Yes	Yes	Yes

Table 12. Implied source-dependent timing accuracy of the compensation models and compatibility with the GOOD-geometry requirements.

Model	Implied ${\mathit{\sigma }}_{\mathit{t}}$	10$\mathbf{m}$	100$\mathbf{m}$	1$\mathbf{k}$$\mathbf{m}$
	($\mathrm{n}\mathrm{s}$)
M0	149668.780	No	No	No
M1	4567.695	No	No	No
M2	6714.422	No	No	No
M3	582.535	No	No	Yes
M4	18.965	Yes	Yes	Yes

Under GOOD geometry, M0–M2 exceed even the 1 $\mathrm{k}$$\mathrm{m}$ timing threshold of $2813.290$ $\mathrm{n}$$\mathrm{s}$. M3, with an implied timing error of $582.535$ $\mathrm{n}$$\mathrm{s}$, satisfies the 1 $\mathrm{k}$$\mathrm{m}$ timing requirement but not the 100 $\mathrm{m}$ or 10 $\mathrm{m}$ requirements.

M4 has an implied differential timing error of $18.965$ $\mathrm{n}$$\mathrm{s}$, below the GOOD-geometry $21.325$ $\mathrm{n}$$\mathrm{s}$ threshold for the 10 $\mathrm{m}$ target. This is consistent with its full-budget RMSE of $9.711$ $\mathrm{m}$. M4 therefore satisfies all three RMSE targets under GOOD geometry under the adopted assumptions.

Under MODERATE geometry, M4 satisfies the 100 $\mathrm{m}$ and 1 $\mathrm{k}$$\mathrm{m}$ timing requirements but cannot satisfy the 10 $\mathrm{m}$ objective because that target is excluded by the non-timing floor. M3 satisfies the MODERATE-geometry 1 $\mathrm{k}$$\mathrm{m}$ timing threshold but not the 100 $\mathrm{m}$ threshold. None of the compensation cases can overcome the POOR-geometry floor.

4.12. Summary of Principal Results

The three experiments establish the following sequence of results.

First, the idealized Earth–Moon reference-clock difference accumulates to $20.483311$ $\mathrm{m}$$\mathrm{s}$ over 365.25 days. Its mean rate is approximately $56.080\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$, but the instantaneous rate varies over a range of approximately $59.205\mathrm{\mu }\mathrm{s}{\mathrm{d}}^{-1}$. A constant rate therefore captures the secular growth but not the complete time-dependent structure.

Second, autonomous compensation performance depends on the choice of temporal basis. M1 reduces the timing RMS to $94.106$ $\mathrm{\mu }$$\mathrm{s}$, whereas the sidereal-only M2 model degrades it to $138.475$ $\mathrm{\mu }$$\mathrm{s}$. The multi-period M3 forecast reduces the RMS to $14.421$ $\mathrm{\mu }$$\mathrm{s}$. The ablation experiment further shows that a synodic-plus-annual basis achieves $4.630$ $\mathrm{\mu }$$\mathrm{s}$, outperforming the full four-period basis. Consequently, additional harmonics are beneficial only when their out-of-sample stability is supported by validation.

Third, the one-day M4 correction product reconstructs the reference series with an RMS of $0.380$ $\mathrm{\mu }$$\mathrm{s}$. This performance is conditional on continued product availability. Increasing the table spacing to seven days raises the residual to $17.764$ $\mathrm{\mu }$$\mathrm{s}$, while a 30-day update outage produces a maximum residual of $569.753$ $\mathrm{\mu }$$\mathrm{s}$.

Fourth, common timing residuals are absorbed by the jointly estimated receiver-clock state and do not directly determine the position error. Source-dependent timing residuals, in contrast, generate model-dependent position degradation. Under the complete error budget and GOOD geometry, the position RMSE decreases from $274.149$ $\mathrm{k}$$\mathrm{m}$ for M0 to $269.094$ $\mathrm{m}$ for M3 and $9.711$ $\mathrm{m}$ for M4.

Finally, timing requirements are conditional on source geometry and the non-timing error floor. Under GOOD geometry, the allowable source-dependent timing errors are $21.325$ $\mathrm{n}$$\mathrm{s}$, $280.675$ $\mathrm{n}$$\mathrm{s}$, and $2813.290$ $\mathrm{n}$$\mathrm{s}$ for the 10 $\mathrm{m}$, 100 $\mathrm{m}$, and 1 $\mathrm{k}$$\mathrm{m}$ targets. The 10 $\mathrm{m}$ target becomes unattainable under MODERATE geometry, and all three targets become unattainable under the deliberately stressed POOR geometry.

5. Discussion

The three experiments delineate a conditional chain spanning from the physical Earth–Moon reference-clock discrepancy to the allocation of precise PNT timing requirements. Our initial analysis corroborates that a lunar reference-clock difference cannot be adequately characterized by a solitary mean drift. Under the adopted first-order weak-field approximation, the accumulated offset exhibits a secular trend modulated by time-dependent structures, which originate from the contrasting gravitational and kinematic environments of the Earth and Moon. For compensation architectures, this implies that adequately modeling the secular rate is insufficient; the residual temporal structure—if left uncompensated—can undergo pronounced amplification over a navigation update interval.

Subsequently, the compensation evaluation explicitly decouples autonomous forecasting from disseminated correction products, treating them under fundamentally different operational constraints. A training-only drift or harmonic model measures autonomous sustainability. In contrast, the tabulated product serves as a reconstruction benchmark, conditional on uninterrupted connectivity. Its superior nominal accuracy should not be conflated with robust autonomous prediction. As demonstrated by the product-outage scenarios, local propagation based on outdated rates rapidly degrades, highlighting a distinct vulnerability in architectures heavily reliant on continuous low-latency links.

A subsequent and critical implication of this framework is the demystification of the one-way range-equivalent timing quantity, $c\Delta t$. Our findings demonstrate that $c\Delta t$ does not translate linearly into three-dimensional position error. Instead, common-mode timing residuals are largely absorbed by the receiver-clock state estimation. Position degradation is exclusively driven by the source-dependent timing residual, mapped through the observation geometry, and inextricably bound to the non-timing error floor. Consequently, arbitrarily tightening the timing allocation yields marginal positioning gains if the system operates near a geometry-dependent performance floor characterized by severe dilution (e.g., our simulated POOR geometry). This mathematical reality mandates a joint geometric-stochastic budgeting paradigm for future lunar constellation designs, rather than treating timing requirements in isolation.

Despite the architectural insights yielded by this framework, several theoretical and operational caveats necessitate rigorous qualification for future studies:

• Relativistic Boundary Limits: The current methodology relies on a first-order weak-field and slow-motion approximation. While sufficient for baseline accumulation estimates, long-term autonomous navigation warrants a formal error upper bound (Upper Bound) that quantifies the truncation errors. The long-term drift caused by ignoring higher-order relativistic terms, Solar System Barycenter (SSB) gravity potential variations, and lunar non-spherical gravity perturbations must be bounded before operational certification.
• Idealized Stochastic Assumptions: In Eq. (72) and Eq. (73), non-timing observation errors (ephemeris, ionosphere/plasma delays, and hardware biases) were modeled as independent, zero-mean Gaussian variables. In actual deep-space operational environments, these error sources are heavily correlated. The assumption of stochastic independence represents an idealized lower bound for the error budget. Future expansions of this framework must incorporate cross-correlation matrices to reflect the coupled nature of cislunar observation hazards.
• Empirical Mission Validation: The validation of the autonomous prediction models (multi-period/harmonic) currently relies on self-generated analytical baselines. To transition from a theoretical framework to an operative standard, these forecasting and compensation strategies must be cross-validated against empirical clock and trajectory data derived from ongoing benchmark missions, such as the NASA Artemis program or the Queqiao relay satellite series.

In brief, this study serves as an architecture-sensitivity methodology rather than a definitive lunar timescale realization or a predictive performance model for a specialized flight constellation. By transparently isolating the transition mechanics from reference-clock modeling to cislunar PNT requirements, this framework provides a stringent, theoretically consistent baseline for forthcoming lunar navigation and time-transfer design studies.

6. Conclusion

This study formalizes an integrated analytical framework that explicitly bridges Earth–Moon relativistic clock modeling with the allocation of cislunar PNT requirements. The fundamental contribution of this work is not merely a quantification of clock drift, but the establishment of a rigorous translation mechanism evaluating how unresolved timing dynamics propagate into deep-space navigation architectures.

Our findings compel a paradigm shift in how timing specifications are prescribed for lunar missions. We establish that defining timing requirements through an isolated, scalar range-equivalent metric ($c\Delta t$) is structurally inadequate. Because strict common-mode clock deviations are absorbed by state estimation, the degradation of user position is governed exclusively by the source-dependent residual fraction, acting in severe conjunction with the prevailing geometric dilution and non-timing error floor. Consequently, an arbitrarily stringent timing clock specification provides zero operational margin if the navigation configuration itself is geometrically degenerate.

Furthermore, the evaluation of distinct compensation strategies exposes a critical architectural dichotomy between autonomous endurance and nominal precision. Relying on high-cadence correction products dictates a system exceptionally vulnerable to data-link outages, whereas leveraging multi-harmonic autonomous forecasting guarantees stability at the cost of bounding the ultimate navigation accuracy.

Moving forward, optimizing cislunar PNT architectures requires abandoning single-variable time-budgeting. The framework presented herein dictates that future lunar constellation designs must adopt a coupled geometric-stochastic requirement allocation, ensuring that onboard clock hardware, compensation algorithms, and orbital geometries are jointly engineered to survive the unique kinematic and gravitational realities of the Earth–Moon system.