Time Transfer From GPS for Designing a SmallSat-Based Lunar Navigation Satellite System

With the advent of new frontiers in space exploration, traveling to the Moon serves as a crucial stepping stone for the success of future deep space missions. Many international space agencies, such as the National Aeronautics and Space Administration (NASA), European Space Agency (ESA), China National Space Administration (CNSA), Japan Aerospace Exploration Agency (JAXA), and Roscosmos State Corporation for Space Activities (Roscosmos), have demonstrated a growing interest in building a sustainable human presence on the Moon (Laurini & Gerstenmaier, 2014). More than 40 lunar missions are being planned within the next decade by 10 international space agencies (Tai et al., 2020). After more than 50 years since the Apollo program ended, NASA’s Artemis mission will land humans on the Moon in the mid-2020s, including the first woman and first person of color (Smith et al., 2020). Furthermore, active efforts on traveling to the Moon are also being invested in by private sector companies like SpaceX and Blue Origin. User position, navigation, and timing (PNT) serves as the basis for various lunar applications that include mission activity planning, search-and-rescue operations, and geotagging calibrated scientific payload data samples (Israel et al., 2020). Stanford University Abstract For the future lunar navigation satellite system (LNSS), a GPS-like constellation for the Moon, NASA has expressed interest in using a smallsat platform. However, many design decisions have yet to be finalized, including the orbit, constellation, and onboard clock. Furthermore, designing a dedicated smallsat-based LNSS limits the payload capacity, thus restricting the size, weight, and power (SWaP) of the onboard clock. We propose an LNSS design with time transfer from intermittently available GPS signals wherein we design a Kalman filter to estimate timing corrections for a low-SWaP clock. Additionally, we devise a lunar User Equivalent Range Error (UERE) metric to characterize the ranging accuracy of signals transmitted by an LNSS satellite. We further validate our time-transfer technique with simulations of an LNSS satellite in an elliptical lunar frozen orbit (ELFO) with an onboard chip scale atomic clock (CSAC) while analyzing the visibility effects of GPS, timing errors, and lunar UERE across ELFOs.

Moreover, with an increase in human and robotic exploration, PNT services will be needed everywhere on the Moon. In particular, the Global Exploration community has expressed the need to achieve a positioning accuracy of less than 50 m on the Moon within the next 10 years (InsideGNSS, 2021).
State-of-the-art methods on PNT estimation for deep space navigation include vision-based terrain relative navigation (TRN) and radio navigation via NASA's deep space network (DSN). However, these techniques do not provide reliable global coverage across the entire lunar landscape. Specifically, TRN techniques provide a map-relative position fix by matching the onboard camera images with a pre-compiled terrain database of distinctive landmarks in space gathered from previous missions. These approaches have demonstrated great success for interplanetary missions, including the recent safe landing of the Mars 2020 Perseverance rover (Bolles, 2020). However, these methods work poorly in the sparsely surveyed regions of the Moon such as the poles and not at all in the darkness of its night and permanently shadowed craters (Christensen & Geller, 2011). Also, TRN techniques require computationally expensive onboard processing and extensive storage (Clerc et al., 2010).
Some prior works applying TRN techniques to exploration missions to the Moon include Downes et al. (2020), Johnson and Montgomery (2008), and Steffes et al. (2019). DSN-based radio navigation downlinks measurements back to ground monitoring stations on Earth and, subsequently, receives uplink with the estimated PNT. However, the PNT information is estimated with a delay, thus limiting how accurately the users on or near the lunar surface, including spacecraft, rovers, and astronauts, can navigate in real time (DeLange et al., 2016). Another drawback is that these radio navigation approaches require a dedicated DSN station (i.e., antenna is to be reserved for downlink/uplink tasks, thus, restricting its scalability with the increase in number of lunar missions). Furthermore, radio navigation is not applicable for a lunar user on the far side of the Moon since the line-of-sight vectors with DSNs on Earth are obstructed. Some relevant works on the use of radio navigation for the lunar context include Christian and Lightsey (2009) and Jun et al. (2020).
Recently, there has been a growing interest in exploring the feasibility of an Earth-like GPS constellation around the Moon. For instance, the NASA Goddard Space Flight Center (GSFC) has conceptualized an initial framework of LunaNet to provide PNT and communication services for future lunar missions (Israel et al., 2020). The ESA has also announced its Moonlight initiative that aims to utilize dedicated lunar satellites to build a navigation and communication infrastructure (Cozzens, 2021). These initiatives aim to address the current gap towards the needs expressed by the Global Exploration community and are, thus, targeting to achieve a positioning accuracy of less than 50 m (InsideGNSS, 2021). However, many design decisions related to the lunar satellite navigation are not yet finalized, such as orbit, constellation, payload, clock, etc. NASA has expressed particular interest in utilizing a smallsat platform for a lunar PNT system due to the cost-effectiveness of this platform and its potential for rapid deployment around the Moon (Israel et al., 2020).
In contrast to the legacy GPS, designing a dedicated smallsat-based lunar navigation satellite system (LNSS) involves major additional challenges (Israel et al., 2020) as follows: (1) The potential inclination of lunar PNT design toward smallsat configuration limits the payload capacity of the LNSS satellites. Given this, the clock onboard the future LNSS constellation will also be limited in the size, weight, and power (SWaP), thus restricting the associated timing stability. Thus, instead of higher-grade atomic clocks based in Cesium and Rubidium, the feasibility of using low-SWaP clocks such as chip scale atomic clocks (CSACs) and micro atomic clocks (MACs) is being explored (Batista et al., 2012). The timing stability of the onboard clock plays a key role in lunar PNT design as it determines the ranging accuracy of transmitted signals from the LNSS satellite. (2) A limited number of ground monitoring stations can be deployed on the Moon, which implies that the LNSS satellites should require fewer station-keeping maneuvers. (3) The Moon has an anisotropic gravitational field that can potentially induce dynamical orbit perturbations. Given the limited lunar ground monitoring stations and the Moon's anisotropic gravitational field, a potential orbit being investigated for lunar PNT in recent literature (Ely & Lieb, 2006) is the elliptical lunar frozen orbit (ELFO). The ELFO minimizes sensitivity to external perturbations for providing persistent, stable coverage to either the North or South Pole with no requirement for station keeping. (4) As compared to that of the GPS, a lower financial budget is allocated toward developing an LNSS constellation.
There exists rich preliminary studies on the design of future lunar navigation constellations. However, none of the existing works account for the SWaP constraints imposed by smallsats in a standalone manner (i.e., without requiring a reliance on ground monitoring infrastructure of Earth or the Moon). Furthermore, they do not estimate a measure of lunar User Equivalent Ranging Error (UERE) and its variation based on the choice of lunar orbit and the grade of the onboard clock of an LNSS satellite.
Similar to the UERE metric defined for a GPS satellite, the lunar UERE provides key insights regarding the LNSS design (at the satellite level) as it characterizes the ranging accuracy of transmitted signals from an LNSS satellite and, thereby, affects the position accuracy achieved by a lunar user. The authors of Murata et al. (2022) designed an ELFO-based constellation whose navigation performance at the South Pole was assessed using dilution of precision (DOP) and satellite visibility. To provide global lunar coverage, prior works (Pereira & Selva, 2020 investigated the trade-offs among constellations based in near-circular polar orbits and frozen orbits by evaluating their DOP, ∆V, and overall cost. Prior work (Schönfeldt et al., 2020a) proposed high-level requirements related to signal frequency and antenna gain by analyzing the DOP and satellite visibility of various hybrid constellations based in an ELFO, near-rectilinear halo orbit, and distant retrograde orbit. Another work (Batista et al., 2012) proposed a Rider constellation of two orbital planes with eight lunar satellites per plane, wherein each lunar satellite was equipped with a low-SWaP CSAC that would be updated hourly by a ground station on the Moon.
Thus, it is worthwhile to investigate if we can harness the signals already broadcast by legacy GPS to alleviate the challenges in designing a smallsat-based LNSS and to address the research gaps in the preliminary LNSS designs. In particular, one can perform the time-transfer technique from GPS to correct the lower-grade clock onboard an LNSS satellite, thus mitigating the need to install and maintain ground station(s) on the Moon and/or the need to schedule one-on-one communication with Earth's ground stations. This is justified because GPS satellites maintain high onboard timing stability, given that they have higher-grade onboard atomic clocks and are aided by a well-established network of Earth ground monitoring stations. In fact, extensive research (Bhamidipati & Gao, 2019;Nie et al., 2007;Sivrikaya & Yener, 2004;Wouters & Marais, 2019) exists on GPS time transfer to maintain high-precision timekeeping at various terrestrial assets, such as low-Earth-orbit satellites, 4G/LTE communication networks, and electrical substations. Moreover, prior works (Chen et al., 2016;Hsu & Jan, 2014;Montenbruck & Ramos-Bosch, 2007;Sun et al., 2017) investigated various state estimation frameworks, such as multilateration, Kalman filters, and real-time kinematic positioning to leverage GPS signals to estimate timing corrections.
While the use of GPS for near-Earth space applications is well-researched, utilizing GPS signals for lunar satellites is not straightforward. This is because the GPS transmit antennas point toward Earth, as shown in Figure 1, thus causing a major part of its main lobe to be occluded by Earth. Thus, signals are received in lunar orbits from the sidelobes and small unoccluded parts of the main lobe of GPS satellites that are located on the far side of Earth, as shown in Figure 1. Furthermore, the GPS satellites are not always available (i.e., only intermittently) due to occultation from not only Earth but also the Moon. In addition to the transmit antenna gain pattern and occultation, the received carrier-to-noise density (C/N 0 ) depends on the transmit antenna power and free-space path loss. Given this, the strength of GPS signals is largely attenuated at lunar distances about 385,000 km from Earth.
Recently, progress has been made toward addressing the challenges in using GPS at lunar distances. In NASA's Antenna Characterization Experiment (ACE; Donaldson et al., 2020), extensive observations of GPS L1 signals were conducted in geostationary orbit to characterize the signal performance in the transmit antenna sidelobes. The reconstructed transmit antenna gain patterns were made publicly available for all GPS satellites in Block IIF, IIR, and IIR-M. In recent years, NASA has also utilized sidelobes of GPS signals in their Magnetospheric Multiscale Mission (MMS) and have successfully achieved a position fix at distances nearly halfway to the Moon (Winternitz et al., 2017).
While GPS signals are not always available in a lunar orbit due to occultations from Earth/Moon and signal attenuation, prior simulation works Schönfeldt et al., 2020a;Winternitz et al., 2019) have demonstrated the presence of continual time segments during which received C/N 0 values were greater than a prespecified threshold. This establishes that successful acquisition and tracking of GPS signals is possible in a lunar orbit in an intermittent manner.
Furthermore, the NASA GSFC has pioneered work in developing the Navigator, which is an autonomous, real-time, fully space-flight-qualified high-sensitivity GNSS receiver (Winternitz et al., 2009). The Navigator has faster acquisition times and the ability to track very weak signals, thus enabling the use of GPS at farther distances from Earth. This receiver is scheduled for testing on the Moon in 2023 during the Lunar GNSS Receiver Experiment (LuGRE), and if successful, it will be the first GNSS positioning fix on the lunar surface (Parker et al., 2022). On a similar note, the SpacePNT has developed a high-sensitivity spaceborne GNSS FIGURE 1 GPS signal at lunar distances is received from the sidelobes and small parts of the main lobe; the received C/N 0 depends on the gain pattern and power of the GPS satellite's transmit antenna, occultation from Earth and the Moon, and free-space path loss. receiver named NAVIMOON that will orbit the Moon onboard the ESA's lunar Pathfinder (Capuano et al., 2016). The ESA's lunar Pathfinder mission will perform the first-ever satellite navigation positioning fix in a lunar orbit by the end of 2023 (Cozzens, 2021). These advances prove the feasibility of using GPS for designing a smallsat-based LNSS.
We propose the design of a smallsat-based LNSS with time transfer from GPS, wherein the LNSS satellites will listen to signals already broadcast by the GPS and process those signals to estimate timing corrections. An illustration of our proposed time-transfer technique is shown in Figure 2. We designed a timing filter that corrects the lower-grade clock onboard an LNSS satellite when GPS signals are available and, when GPS signals are unavailable, we propagated these clock estimates forward in time. We also developed a GPS Continual Outage Period (ECOP) metric to analyze the visibility effects of GPS on the timing stability of the onboard clock.
Our proposed GPS-to-LNSS time-transfer technique alleviates requirements with regards to timing stability and, correspondingly, the SWaP of the onboard clocks. Additionally, our design leverages GPS signals, thereby eliminating the need for additional ground infrastructure on the Moon to maintain LNSS services. Furthermore, our design facilitates navigation across the complete lunar landscape with the LNSS constellation. This work is based on our recent ION GNSS+ 2021 conference paper, Design Considerations of a Lunar Navigation Satellite System with a Time-Transfer from Earth-GPS (Bhamidipati et al., 2021).

Key Contributions
The key contributions of our smallsat-based GPS-to-LNSS time-transfer technique are as follows: 1. We developed a timing filter based on a Kalman filter framework to update the lower-grade clock onboard an LNSS satellite using intermittently available GPS signals. 2. We devised a lunar UERE metric to characterize the ranging accuracy of an LNSS satellite. For a given orbit and an onboard clock, the lunar UERE depends on the root-mean-square (RMS) timing error estimated using our proposed time-transfer technique. This is because any timing error is directly reflected in the ranging accuracy. 3. We validate our time-transfer technique using high-fidelity simulations of an LNSS satellite in ELFO with an onboard CSAC. We also performed a comparison analysis across different ELFO configurations by varying the Keplarian parameters in terms of various metrics, such as maximum ECOP, timing errors, and lunar UERE. Additionally, we performed a sensitivity analysis of lunar UERE as the broadcast ephemeris errors are varied.
The rest of the paper is organized as follows. Section 2 describes our proposed time-transfer technique from GPS to LNSS. Section 3 describes the details related to the modeling of our simulation setup, including the GPS constellation, attenuation of GPS signals, and LNSS satellites in various ELFOs. Section 4 discusses the implementation details and results of the simulated experiment setup. Section 5 provides concluding remarks.

PROPOSED TIME TRANSFER FROM GPS TO LNSS
We first describe the details of our timing filter and, later, explain the formulation of the lunar UERE metric.

Clock Estimation via Timing Filter
We designed a timing filter using a Kalman-filter-based state estimation method (Coleman & Beard, 2020;Krawinkel & Schön, 2015) to continuously maintain precise timing onboard the LNSS satellite. The high-level architecture of our timing filter that leverages intermittently available GPS signals is shown in Figure 3. We considered the LNSS satellite to be equipped with a spaceborne GPS receiver and a lower-grade clock. Particularly, we designed a Kalman filter framework that performs time and measurement updates to compute the clock estimates. Additionally, we aided the timing filter with position and velocity of the LNSS satellite that is obtained from known/pre-computed ephemeris.
We, first, check if GPS satellites are visible by analyzing the ECOP metric. During ECOP, we only perform a time update to propagate the current clock estimate using the propagation model described in Section 2.1.1. When GPS signals are available, we perform both a time and measurement update to correct the clock estimate, where the measurement update is described in Section 2.1.2.
At any time step t, the state vector of our timing filter comprised of LNSS clock estimates (bias and drift) is defined as follows: where b t is the clock bias in m and ḃ t is the clock drift in m/s. In particular, b t and ḃ t are converted from timing bias and timing drift errors to range and range rate errors via multiplication by the speed of light c = 299, 792, 458 m/s.

Time Update Step
At each time step, our timing filter propagates the current clock estimates, which is comprised of clock bias and clock drift, forward according to the following integrated random walk process: where x t is the propagated, or predicted, state vector at time step t, P t is the propagated state error covariance matrix at time step t, and Δt is the update time interval (i.e., the time difference between consecutive time steps of the propagation model).
The process noise covariance matrix Q in Equation (3) is defined as: where h 0 , h −1 , and h −2 are the power spectral density coefficients obtained from the Allan deviation plot (Krawinkel & Schön, 2015). These coefficients provide insights regarding the short-term and long-term stability of the onboard clock and are independent of the LNSS satellite orbit being considered (Krawinkel & Schön, 2015). Note that the squared speed of light is added as a scaling factor in the process noise covariance matrix to maintain consistency with the units of state vector in Equation (1).

FIGURE 3
High-level architecture of our timing filter; the LNSS satellite is equipped with a lower-grade clock and utilizes a GPS receiver with position and velocity aiding to correct for clock bias and drift errors.

Measurement Update Step
Unlike the time update, which is executed at each time step, we only perform the measurement update to correct the LNSS clock estimates when GPS signals are available (i.e., no ECOP). Furthermore, we utilize the LNSS satellite position and velocity from known/pre-computed ephemeris to aid our timing filter. Indeed, less frequent updates, depending on the orbit, are needed to correct the lunar orbit perturbations and are, thus, independently handled outside of our proposed time-transfer filter, such as via an onboard orbit determination module (Erdogan & Karslioglu, 2009) or via infrequent updates from Earth ground monitoring stations. Specifically, we utilize the known satellite position and velocity of the LNSS satellite to compute the expected range and range rate measurements (no clock effects incorporated here), which are later compared with received GPS measurements to formulate the input for the measurement update step.
At time step t, our measurement vector comprised of clock bias and clock drift observables from GPS is as follows: where N is the total number of visible satellites at time t and  ρ t k ( ) and ρ t k ( ) are the clock bias and clock drift observables, respectively, from satellite k. The clock bias observable is defined as the difference between the received GPS pseudorange and the expected range measurement. Similarly, the clock drift observable is defined as the difference between the received GPS pseudorange rate and the expected range rate measurement. These clock bias and clock drift observables for any satellite k can be mathematically represented as: where ρ t k ( ) and  ρ t k ( ) are the received pseudorange and pseudorange rate measurements, respectively, from the received signal of GPS satellite k at time step t. Similarly, ρ t k ( ) and  ρ t k ( ) are the expected range and range rate measurements, respectively, from GPS satellite k at time step t. Here, the expected range from GPS satellite k is calculated as the 2-norm of the position difference between GPS satellite k and LNSS satellite. Similarly, the expected range rate from GPS satellite k is calculated as the projection of the expected velocity difference between the GPS satellite k and the LNSS satellite along the line-of-sight vector from the GPS satellite k to the LNSS satellite. The measurement update step of the timing filter, thus, utilizes the received measurement vector z t to compute the corrected clock estimates ˆt x and covariance ˆt P via the following standard expressions: where K t is the Kalman gain matrix, R t denotes the measurement noise covariance matrix, and C N ∈ ×  2 2 is the measurement model associated with clock bias and clock drift observables such that C a matrix with size a b × of ones and zeros, respectively. I a denotes the identity matrix of size a a × . We adaptively formulate R t based on the tracking errors in the delay lock loop (DLL) and phase lock loop (PLL), which depend on the received C/N 0 of each GPS satellite signal. We model the measurement noise covariance matrix R t as: where σ ρ t k 2 ( ) and σ ρ  t k 2 ( ) correspond to the variance of the clock bias and clock drift observables, respectively, from the GPS satellite k.
We model the measurement variance of clock bias observable as the following, in m 2 : where the first term σ DLL t k 2 ( ) is the variance in the GPS satellite k due to the receiver DLL tracking error, the second term σ UERE, GPS 2 is the variance due to the space and control segment errors of the GPS system, including the broadcast clock and ephemeris errors (Kaplan & Hegarty, 2017). σ eph, LNSS 2 is the variance due to the uncertainty in the LNSS satellite ephemeris, which is used for aiding the timing filter. The multiplicative factor ( ) cT c 2 is added to scale the receiver DLL tracking error denoted in chips 2 to an error in clock bias observable denoted in m 2 . Note that for any GPS satellite k and time step t, we consider the second term σ UERE, GPS 2 and third term σ eph, LNSS 2 to be pre-defined values that are set during initialization, which is provided later in Section 4.1. Additionally, while errors in LNSS satellite ephemeris contribute to time-varying biases in clock bias and clock drift observables, we overbound these errors using an zero-mean Gaussian noise whose covariance is represented by σ eph, LNSS 2 .
We model the uncertainty due to the DLL tracking loop errors, represented as the variance σ DLL t k 2 ( ) in Equation (13), according to the following formulation from Kaplan and Hegarty (2017) and utilized in Capuano et al. (2015aCapuano et al. ( , 2015b for designing a non-coherent power DLL discriminator, in chips 2 : where B DLL is the code tracking loop bandwidth, d is the early-late correlator spacing, T is the coherent integration time, T c is the chipping period of the GPS L1 C/A signal at a chipping rate of 1 023 10 6 . ⋅ chips/s, such that T c = ( ) ≈ ⋅ 1 1 023 10 6 9 775 . . , Hz ns and B fe is the double-sided front-end bandwidth. We correspondingly model the uncertainty in the clock drift observable using the variance associated with the receiver PLL tracking error from GPS satellite k (Borio et al., 2011;Capuano et al., 2015a) as the following: where

Formulation of the Lunar UERE Metric
We formulated a lunar UERE metric to characterize the accuracy of the transmitted ranging signals from the LNSS satellite to users on or near the lunar surface, including spacecrafts, rovers, and other satellites, that assist in lunar missions. Table 1 lists the relevant error components to be used for the lunar UERE formulation while considering the error components of the GPS UERE from Kaplan and Hegarty (2017) as a reference. In contrast to the GPS constellation and terrestrial applications on Earth, the lunar environment has minimal atmospheric effects. As a result, the tropospheric and ionospheric delay error components in the lunar UERE can be considered to be negligible. Furthermore, we consider the multipath error component of the UERE to be negligible due to the lack of building infrastructure and foliage on the Moon. We additionally assume that the antenna gain pattern of the LNSS receiver will mitigate reflections from the ground, similar to terrestrial GPS receivers and, thus, consider ground multipath to also be negligible.
As shown in Table 1, the four error components that we consider in the formulation of lunar UERE are due to the broadcast clock, differential group delay, broadcast ephemeris of the LNSS satellite, and user LNSS receiver noise. Note that the errors due to the differential group delay σ gd, LNSS and receiver noise σ rec, LNSS greatly depend on the final signal structure for the LNSS constellation and tracking loops chosen for the user LNSS receiver, respectively.
Given that the broadcast ephemeris of LNSS satellites is considered to be known/ pre-computed, the associated error component σ eph, LNSS is modeled accordingly. In general, the broadcast ephemeris component for an LNSS satellite in a lunar orbit should account for errors induced via dynamic orbit perturbations from the anisotropic gravitational field of the Moon and the lower frequency of correction updates from ground monitoring station(s). In this work, we simulate and analyze the error due to LNSS broadcast clock inaccuracies, σ clk, LNSS , based on our proposed time-transfer technique from GPS, which, in turn, depends on the orbit and clock considerations for the LNSS satellite as well as the associated errors in the LNSS satellite ephemeris. Correspondingly, note that since the LNSS satellite ephemeris is used to aid our timing filter, which corrects the onboard clock, its error will be reflected twice in the formulation of the lunar UERE: Once explicitly in the broadcast ephemeris error and once implicitly in broadcast clock error.
Thus, the lunar UERE, which depends on the RMS of the four error components described above, is formulated as the following:

MODELING THE LUNAR SIMULATION SETUP
We designed a simulation experiment with the start time as Nov 9, 2025 at 00 : 00 : 00.000 UTC for a duration of 2 months. As described in Section 1, we consider the LNSS satellite to be equipped with an onboard CSAC and a spaceborne GPS receiver. We conducted our experiments for an LNSS satellite in ELFO and with an onboard CSAC, but in general, our proposed time-transfer technique would apply to any orbit and clock type.
We developed a high-fidelity simulation setup using the Systems Tool Kit (STK) software by Analytical Graphics, Inc. (AGI) and the MATLAB software developed by MathWorks. In particular, the STK software provides a user-friendly graphic interface for modeling complex systems and interactions inside a realistic and time-dynamic 3D simulation that includes high-resolution terrain, imagery, radio frequency environments, gravitational effects, solar radiation pressure, orbital propagation, and so on. STK allows the users to import a variety of generic objects of ground, sea, air, and space assets such as satellites, vehicles, facilities, transmitters, and receivers, that can combine to represent existing or proposed systems.
For instance, one can create a simulation of a satellite in orbit by, first, setting a desired simulation time and then inserting an in-built satellite object with desired properties related to central planetary body, orbital parameters, mass, propagator type, and altitude. This software can simulate the entire system-of-systems in action at any location and at any time that can later be easily interfaced with the MATLAB software to extract the relevant time-series of data outputs. Furthermore, the STK software allows users to determine the times one object can access (or see) another object while providing the flexibility to impose desired constraints on accesses between objects. For example, one can create an access link between a transmitter object and a receiver object to retrieve the link budget report that is comprised of time-series of fields related to signal quality, such as effective isotropic radiated power, power flux density at the receiver antenna, C/N 0 , and receiver radio frequency bandwidth. Our simulation environment in STK models the transmission and reception of GPS signals to simulate the visibility and C/N 0 at an LNSS satellite, which is explained in Sections 3.1 and 3.3. In particular, we modeled the orbits of the GPS and LNSS satellites, transmit power, and antenna gain patterns for the main lobe and sidelobes of each GPS satellite, attenuation of GPS signals due to free-space path loss, occultation due to Earth and the Moon, and spaceborne GPS receiver characteristics, such as antenna gain and tracking time.
To validate our proposed time-transfer technique, we generated simulations with an LNSS satellite in different ELFOs, which is explained in detail in Section 3.2. To model these various aspects, we employed four object types in STK, namely satellite, transmitter, receiver, and sensor objects, of which the satellite and transmitter objects were used to model GPS satellites while the satellite, receiver, and sensor objects were used to model LNSS satellites. Note that, for an LNSS satellite, the sensor object acts as an interface between the satellite and the receiver objects to provide steering properties to the receiver antenna such that it is always pointed toward Earth. We create access links between the transmitter object of each GPS satellite to the receiver object of the LNSS satellite to extract the received C/N 0 of GPS signals at the LNSS satellite.
We extracted the generated data from STK simulations into MATLAB where the rest of the necessary modules were simulated. In particular, for an LNSS satellite, we propagated the clock dynamics model to generate true clock estimates (explained in detail in Section 3.4). Additionally, we utilized the received C/N 0 and the pre-computed LNSS ephemeris errors to induce simulated errors that mimic realistic GPS measurements of clock bias and clock drift observables. We provide detailed explanations of these simulation steps in Section 3.5. These simulated measurements were given to our timing filter that was previously described in Section 2.
We formulated a conservative approach for demonstrating our algorithm's performance. In our simulation, we mitigated the terrestrial ionospheric and tropospheric effects in GPS signals by considering an Earth mask. We executed this by disregarding GPS signals that passed by Earth at altitudes less than or equal to a pre-specified value of Earth mask. Since the troposphere stretches up to 10 km and the ionosphere extends up to 965 km above Earth's surface, we chose an Earth mask of 965 km. Note that, instead of discarding these GPS signals, modeling atmosphere effects or using a second frequency to form ionospheric-free GPS measurements are other viable options. Exploring these options is left as an extension for future work. While our current paper only considers GPS, note that our proposed time-transfer technique can be extended to include other GNSS constellations as well.

Modeling the GPS Constellation in STK
To create each of the GPS satellites in STK, we inserted a satellite object with system properties loaded from a pre-compiled AGI catalog of GPS almanac data. We utilized the almanac obtained from the online AGI server to propagate the position and velocity states of the GPS satellites for the desired experiment duration. Specifically, we considered the date of the almanac to be October 25, 2020, at 17 : 04 : 00.000 UTC with week reference as April 7, 2019.
We considered our simulated GPS constellation to be made up of 31 operational satellites allocated in six orbital planes as seen in Figure 4(a). Our GPS constellation was comprised of eight satellites from Block IIR, seven satellites from IIRM, 12 satellites from IIF, and four satellites from Block III. We utilized the transmit power and antenna gain patterns of the L1 C/A signals for Blocks IIR, IIRM, and IIF from the NASA GPS ACE (Donaldson et al., 2020) to model the transmit antenna onboard the GPS satellites in our simulation setup.
For the sake of simulation simplicity, we assigned the same antenna gain pattern for satellites in each block of the GPS constellation, which is equivalent to the average experimental antenna gain pattern files made publicly available by Donaldson et al. (2020). Given that no data is publicly available for Block III satellites, we approximated the transmit antennas of Block III using the transmit power and average antenna gain pattern available in Donaldson et al. (2020) for Block IIF. To simulate a GPS transmit antenna, we attached a transmitter to the face side of the GPS satellite object and directed it toward the center of the Earth.
Given that no pre-compiled patterns exist in STK that could model the sidelobes of the GPS transmit antenna, we loaded a custom gain pattern file that we created from the results of NASA's ACE study into the transmitter object. The illustrations of these antenna gain patterns are provided in Figures 4(b), 4(c), and 4(d), for Block IIR, IIRM, and IIF, respectively. We additionally assigned the Block III satellites the same antenna gain pattern as that of Block IIF. The transmit power for the antenna gain pattern of the Block IIR satellites is 15 dBW, while for the Block IIRM, IIF, and III satellites, we used a transmit power of 14.3 dBW for the antenna gain pattern.
The antenna gain pattern files in Donaldson et al. (2020) only contained data for angles ranging between 16°−90° for off-boresight and 0°−360° for azimuth. Similar to the notation in Delépaut et al. (2020), an off-boresight angle of 0° represents the nadir direction of the GPS satellites, while 90° points towards the GPS satellite velocity/anti-velocity direction. Given that the NASA ACE does not include data for off-boresight angles of less than 16°, we assigned the antenna gain for any GPS satellite between 0°−16° with a constant value that equals the average value obtained for the 16° off-boresight angle across the azimuth values ranging between 0°−360°. This approximation had minimal effects on the simulation results, since off-boresight angles of less than 16° primarily make up the part of main lobe that is largely occluded by Earth.

Modeling the LNSS Satellite in STK: ELFO
We created realistic simulations for an LNSS satellite in ELFO by utilizing the High-Precision Orbit Propagator (HPOP) setting in the STK software (AGI, 2021). The HPOP builds upon precise force models of Earth, the Sun, and Moon for generating and propagating accurate orbit solutions of the LNSS satellite position and velocity (Liu et al., 2019). As mentioned earlier in Section 1, ELFOs have been the subject of great interest for LNSS development as they are highly stable and require significantly fewer station-keeping maneuvers compared to other types of lunar orbits.
We modeled ELFOs in STK by creating a satellite object and then defining its orbit using classical orbit mechanics, which states that any orbit requires six elements (i.e., Keplerian parameters) to propagate the position and motion of the satellite fully. Specifically, we referred to prior literature (Ely & Lieb, 2006;Schönfeldt et al., 2020b) for the Keplerian parameters of each ELFO, which include semi-major axis, eccentricity, inclination, argument of perigee, right ascension of the ascending node (RAAN), and mean anomaly, and are listed in Table 2. Figure 5 shows an illustration of the three ELFOs in the Moon-inertial frame. Particularly, we chose these three ELFOs given their persistent coverage of the Moon's South Pole, which is a key region of interest for future lunar missions due to the presence of water ice deposits (Li et al., 2018).

Modeling the LNSS Satellite in STK: Spaceborne GPS Receiver
We simulated a spaceborne GPS receiver onboard an LNSS satellite by modeling the GPS signal reception in STK software up to the computation of C/N 0 values. To maximize the visibility of GPS signals at the LNSS satellite, we designed a steering receiver antenna. In particular, we modeled steering characteristics in STK by creating a sensor object that was mounted on the LNSS satellite object and then constraining it to point toward Earth. We, then, attached an STK receiver object on the sensor object, which ensured that the GPS receiver onboard the LNSS satellite would always be Earth-pointed.
We selected a parabolic antenna type in the properties of the receiver object and, thereby, modeled a high-gain antenna with 14 dBi at 0° off-boresight angle and a 3-dB beamwidth of 12.2°. An illustration of the associated antenna gain pattern is shown in Figure 6. As explained in Delépaut et al. (2020), this antenna pattern is representative of the current design of the Pretty CubeSat mission to be developed for GNSS reflectometry applications (Fragner et al., 2020). Furthermore, in regard to the system noise temperature, we set a 30-dB gain for the low-noise FIGURE 5 Illustration of three ELFOs considered for validating our proposed time-transfer technique, where ELFO #1 is at an altitude of 6,541.4 km, ELFO #2 is at 7,500 km, and ELFO #3 is at 9,750.5 km FIGURE 6 Our simulated spaceborne GPS receiver is equipped with a high-gain steering receiver antenna with 14 dBi at 0° off-boresight angle and a 3-dB beamwidth of 12.2° amplifier (LNA) with a noise figure of 2 dB and an antenna noise temperature of 113 K. For simplicity, we did not consider implementation losses and cable losses in our simulation.

Modeling the LNSS Satellite in MATLAB: CSAC
Given the smallsat platform for the LNSS, the limited payload capacity restricts the SWaP of the onboard clock, which makes a CSAC a potentially useful option for a low-SWaP clock type. For our work, we simulated the onboard clock in accordance with the Microchip CSAC (Calero & Fernandez, 2015), which has a time deviation (TDEV) of 1.5 s µ per day and a SWaP of (17 cm 3 ·0.0035 kg·0.12 W) (Schmittberger & Scherer, 2020). For the onboard CSAC, we modeled the true clock model as having a constant drift in the clock bias, with the propagation model defined in Equation (2). In particular, we modeled the constant drift to be equal to the TDEV in order to observe the average drift consistent with the Microchip CSAC clock error model.
However, for our timing filter, we utilized the covariance matrix model Q defined in Equation (5), which has been derived in previous works. Since our work focuses on proposing an LNSS design that utilizes a time-transfer technique from GPS, we employed a simple two-state clock error model as a proof-of-concept. However, note that there exists rich literature on advanced clock modeling that can simulate additional effects like temperature variations and stability effects based on the duration for which the clock is steered. When using advanced true clock simulations, alternate state estimation techniques based on a three-state clock model (Zucca & Tavella, 2005) and four-state clock model (Shin et al., 2008;Van Buren et al., 2019) can be explored for our timing filter. These aspects will be investigated in our future works.

Modeling the LNSS Satellite in MATLAB: Received GPS Measurements
We analyzed the visibility statistics extracted from STK as explained earlier in Section 3.3 and induced simulated errors in pseudoranges and pseudorange rates in MATLAB, where the rest of our proposed time-transfer technique has been implemented.
Prior work (Musumeci et al., 2016) demonstrated that acquisition and tracking could be achieved for a low C/N 0 threshold of 10 dB-Hz while other works (Capuano et al., 2015a(Capuano et al., , 2015bDelépaut et al., 2020) utilized a threshold of 15 dB-Hz. To ensure a conservative margin, we set the acquisition and tracking threshold for our simulated spaceborne GPS receiver as 15 dB-Hz. Figure 7 shows the received C/N 0 at ELFO #1 from PRN 7 of the GPS constellation in orange dotted lines, where PRN stands for pseudorandom noise and is used as an identifier for a GPS satellite vehicle. The threshold of 15 dB-Hz is indicated by the black dashed line in the figure. We conservatively set the minimum required time duration of consistent signal visibility for acquisition and tracking to be 40 s (i.e., a GPS satellite is considered to be visible when the received C/N 0 value is greater than 15 dB-Hz for a continuous time duration of at least 40 s).
To simulate the received GPS measurements utilized in the measurement update step described earlier in Section 2.1.2, we first obtained the true range and true range rate between the GPS and LNSS satellites from STK and incorporated the simulated true clock bias and clock drift, as explained previously in Section 3.4. To generate a pseudorange-based clock bias observable, we then subtracted the true ranges with that of the expected ones computed using the GPS and LNSS satellite ephemeris from STK, while introducing stochastic errors according to the measurement error model described in Equations (13) and (14).
For this simulation, we sampled the position and velocity of GPS and LNSS satellites at a pre-specified sampling period of 60 s. For the LNSS satellite ephemeris error, we sampled this value only once per time step and incorporated the same measurement bias across all simulated clock bias observables. As a result, the simulated measurement vector z t had correlated errors; however, while formulating the measurement covariance matrix R t of our timing filter, we modeled the measurements as having uncorrelated errors.
For the pseudorange-rate-based clock drift observables, we subtracted the true range rate measurements with that of the expected ones and, thereafter, induced stochastic errors based on only the simulated uncertainties from the PLL tracking loops. Modeling other sources of uncertainties in the clock drift observable, including lunar velocity estimation errors, is left as an extension for our future work. For the sake of simulation simplicity, we did not consider light time delay effects (i.e., change in position and velocity of the GPS satellite between the transmit time and received time at the LNSS satellite). Incorporating a light time delay to further improve the simulation setup has been left as an extension for future work.

EXPERIMENTAL RESULTS AND ANALYSIS
We validated our proposed time-transfer technique using the simulated cases of an LNSS satellite with an onboard CSAC orbiting in various ELFOs, which were described earlier in Section 3.

Implementation Details
We designed our timing filter with an update period of ∆t = 60 s. To formulate the process noise covariance matrix described earlier in Equation (1), we extracted the power spectral density coefficients from the Allan deviation plot of CSACs 3 74 10 . . In the formulation of the measurement noise covariance matrix described earlier in Equation (12), the uncertainty due to the DLL tracking loop errors from Equation (14) depended on the tracking loop parameters, which were set to the following: B DLL = 0.5 Hz, d = 0.3 chips, T = 20 ms, and B fe = 26 MHz (Capuano et al., 2015a). To characterize the clock drift observable covariance described in Equations (15) and (16), which depend on the PLL tracking loop errors, we set B PLL = 0.5 Hz. Note that the same parameters were considered for simulating the received pseudorange and pseudorange rate measurements, as explained previously in Section 3.5.
To characterize the lunar UERE from Section 2.2, we considered the same group delay and receiver noise error magnitudes as with the GPS system (i.e., σ gd, LNSS = 0.15 m and σ rec, LNSS = 0.1 m). Compared to GPS, which is a terrestrial navigation system with extensive military applications, the LNSS proposed would have greater limitations on ground infrastructure and financial investment to monitor the lunar system. Given that no prior lunar missions have been conducted in ELFOs yet, we perform a sensitivity analysis of lunar UERE later in Section 4.6 by modeling the RMS error in LNSS satellite position to range from σ eph, LNSS = 0.3 m to σ eph, LNSS = 300 m in multiplicative increments of 10.
Our reasoning behind the choice of these values is multi-fold: a) The lower bound of σ eph, LNSS = 0.3 m denotes the broadcast ephemeris error component considered for legacy GPS currently; b) The values ranging between σ eph, LNSS = 30 m to σ eph, LNSS = 300 m have been demonstrated by past lunar missions (Liu et al., 2020;Mazarico et al., 2018), some of which include Chang'e 1, Chang'e 2, Chang'e 3, Chang'e 4, Chang'e 5T1, and the missions of the Lunar Reconnaissance Orbiter (LRO). Note that, for the analysis conducted later in Sections 4.4 and 4.5, we formulated a baseline case of σ eph, LNSS = 3 m which equals the error due to the broadcast ephemeris of GPS scaled by an order of magnitude.

Validation Metrics
We chose validation metrics to analyze the top considerations for designing an LNSS with our proposed time-transfer technique from GPS. We analyzed our algorithm's performance using the following four validation metrics: a) Satellite visibility, which indicates the percentage of time in the entire experiment's duration for which the number of visible GPS satellites were greater than a pre-specified threshold; b) Maximum ECOP to identify the region of maximum continuous time when no GPS satellites were visible; c) Maximum errors in clock estimates to analyze the worst-case performance of our proposed time transfer; and d) The lunar UERE metric to characterize the ranging measurement accuracy of signals transmitted by an LNSS satellite.
To characterize the satellite visibility for a given ELFO, we examined the percentage of time in which a certain minimum number of GPS satellites are visible. In particular, we examined the percentage of time that at least one GPS satellite was visible, since that is the minimum number required to correct the clock estimates, as well as the percentage of time that at least four GPS satellites were visible, since that is the minimum number to perform full state estimation involving position, clock bias, velocity, and clock drift. As previously explained in Section 2, the lunar UERE depends on the RMS timing error associated with our proposed time-transfer technique from GPS. Among the three ELFOs described in Table 2, we considered a particular ELFO to be desirable if it achieved some or all of the following: greatest satellite visibility, least maximum ECOP, least maximum timing error, and least lunar UERE.

Variation in Satellite Visibility and Maximum ECOP
Figures 8(a) through 8(f) demonstrate the number of visible GPS satellites for the three simulated cases of ELFOs, which is shown in blue. The highlighted red vertical bars indicate the regions of ECOP. Based on Table 3, we observe that ELFO #3 achieved the maximum duration of at least one satellite visibility, which was 99.2% of the total time, as well as the maximum duration of at least four satellite visibility, which was 92.1% of the total time. Intuitively, this observation seems reasonable given that ELFO #3 has a maximum altitude of 9,750.5 km and, thus, experiences fewer occultations from the Earth and the Moon. Furthermore, ELFO #1 exhibited the least maximum ECOP of 3,060 s. Intuitively, this observation seems reasonable since the velocity of an LNSS satellite is the highest in the ELFO with the lowest altitude; therefore, this ELFO could pass through regions of occultation quicker. From Figures 8(a), 8(c), and 8(e), we observed that ELFO #1, which had the lowest altitude and, hence, the lowest orbital period, exhibited five occurrences of maximum ECOPs (duration of 3,060 s) in an experiment duration of 2 months, while only one occurrence of maximum ECOP was observed for ELFO #2 and ELFO #3. On the other hand, from Figures 8(b), 8(d), and 8(f), we observed that ELFO #3, which had the highest altitude, also exhibited the largest time spacing between ECOPs as compared to ELFO #1 and ELFO #2.

Variation in Maximum Timing Errors
For a given CSAC onboard an LNSS satellite and σ eph, LNSS = 3 m, we compared the maximum errors in clock bias and clock drift across three ELFOs in Figure 9 with ELFO #1 indicated in cyan, ELFO #2 in magenta, and ELFO #3 in green. In the plots, we indicate the ECOP regions of ELFO #1 as vertical red bars, during which we correspondingly observed a steady increase in the satellite's clock bias errors. The slope of the clock bias error during ECOP is governed by the estimation error in clock drift right before the LNSS satellite enters the ECOP time duration, where no GPS satellites are visible.
Intuitively, maximum timing error indicates the worst-case performance of our proposed time transfer and, thus, the ELFO exhibiting the least maximum timing error is more favorable for the LNSS design in question. We observed that ELFO #2, which was at an altitude of 7,500 km, exhibited the least maximum timing error in clock bias of 0.26 μs and clock drift of 0.24 ns/s. Based on this, we observed that the maximum timing error not only depends on the maximum ECOP and altitude, but also other orbital parameters that govern the geometric configuration between GPS and LNSS, including eccentricity and inclination. As an alternative to this error analysis conducted in the time domain, Allan deviation plots of our time-transfer technique can be analyzed by executing multiple Monte Carlo runs, which is to be explored in our future works.

Comparison Analysis of Lunar UERE
For a given CSAC onboard an LNSS satellite and σ eph, LNSS = 3 m, we determined the lunar UERE metric for different ELFOs. As mentioned previously in Section 2.2, an important component of the lunar UERE is the broadcast clock error, which corresponds to the estimated RMS clock bias error associated with our proposed time-transfer technique. Table 4 reports the RMS timing errors in clock bias and clock drift according to which the ELFO #3 exhibited the least RMS error in clock bias of 7.9 m. With pre-specified RMS values of other error sources listed earlier in Section 4.1, we validated that our time-transfer technique could achieve a low lunar UERE of less than 10 m for an LNSS satellite in ELFO, even with a lower-grade clock onboard (i.e., a CSAC). Furthermore, we observed that the estimated lunar UERE for all three ELFOs was comparable in order of magnitude to that of GPS, thus validating our LNSS design which utilizes time transfer from GPS to lower the SWaP requirements of the onboard clock. Note: The estimated lunar UERE across all ELFOs had a comparable order of magnitude to that of the GPS UERE. Note that unlike Figure 9 that showcases maximum timing errors, this table reports the RMS errors with units of clock bias listed in m and that of clock drift in m/s.

Sensitivity Analysis of Lunar UERE as Broadcast Ephemeris Errors Varied
Considering the same low-SWaP onboard CSAC for all cases, Table 5 shows the lunar UERE variations for the other three cases of broadcast ephemeris errors with σ eph, LNSS = 0.3 m, σ eph, LNSS = 3 m, σ eph, LNSS = 30 m, and σ eph, LNSS = 300 m. Irrespective of the broadcast ephemeris error, we demonstrated that the variation in lunar UERE across different ELFOs would be minimal.
As expected, we observed a near-quadratic decrease in lunar UERE as the broadcast ephemeris errors decreased. This is because the lunar UERE was doubly dependent on the error in LNSS satellite ephemeris: Once explicitly in the broadcast ephemeris component σ eph, LNSS and once implicitly in broadcast clock component σ clk, LNSS (aids our timing filter).
Based on the lunar missions conducted so far, the ephemeris errors have been steadily decreasing from 230 m in the 2007 Chang'e 1 mission to 20 m in the recent LRO mission. Given this, we expect that the effect of broadcast ephemeris error on lunar UERE from our proposed time transfer from GPS to be less dominant.

CONCLUSION
We designed a smallsat-based LNSS with a time-transfer technique from GPS to alleviate the size, weight, and power (SWaP) as well as the timing stability requirements of the onboard clocks. We developed a timing filter that corrects the lower-grade clock when GPS signals are available and propagated these clock estimates forward in time when no GPS signals were available.
We analyzed our proposed time-transfer technique using high-fidelity simulations of an LNSS satellite with an onboard chip scale atomic clock (CSAC) for three cases of elliptical lunar frozen orbits (ELFOs). We validated that the least maximum GPS continual outage period (ECOP) of 3,060 s was observed for an ELFO with a semi-major axis of 6,541.5 km.
We demonstrated that a low lunar user equivalent range error (UERE) of less than 10 m, comparable in order of magnitude to that of the GPS UERE, can be achieved for all ELFOs, even with a lower-grade clock onboard (i.e., a CSAC). Additionally, by performing a sensitivity analysis, we demonstrated a near-quadratic decrease in lunar UERE across all ELFOs as the broadcast ephemeris error was reduced. Thus, our proposed time-transfer technique from GPS serves as a promising technique to maintain precise timing onboard the future lunar PNT constellation. Note: The estimated lunar UERE across all ELFOs exhibited a near-quadratic decrease in value as the broadcast ephemeris error decreased.