Sensitivity analysis of polar orbiter motion to lunar viscoelastic tidal deformation

Abstract

We investigate the impact of viscoelastic tidal deformation of the Moon on the motion of a polar orbiter. The dissipative effects in the Moon’s interior, i.e., tidal phase lags, are modeled as Fourier series sampled at given frequencies associated with linear combinations of Delaunay arguments, the fundamental parameters describing the lunar motion around the Earth and the Sun. We implement the tidal model to evaluate the temporal lunar gravity field and the induced perturbation on the orbiter. We validate the numerical scheme via a frequency analysis of the perturbed orbital motion. We show that, in the case of the Lunar Reconnaissance Orbiter at a low altitude of less than 200 km, the main lunar tides and hence the potential Love numbers around the monthly and some multiple frequencies are dynamically separable. The omission of those effects in practice introduces a position error at the level of a few decimeters within 10 days.

1 Introduction

The Moon deforms as a result of tidal forces raised by the Earth, and to a lesser extent, the Sun (Harrison 1963; Williams and Boggs 2015). The deformation manifests in surface displacement as well as a concurrent change in the body’s gravitational field, which are dependent on the structure and material properties of the interior. The susceptibility or resistance of the body to deformation is described by the Love numbers, h, l for the vertical or horizontal surface variations and k for the induced gravitational potential. The external forcing occurs over a wide spectrum of frequencies governed by the relative motions of the Moon-Earth-Sun system. The Love numbers vary across the dynamic spectrum in amplitude and phase with respect to the excitation. While a perfectly elastic body deforms instantaneously, i.e., the phases of its Love numbers are zero, the viscosity delays the body response and causes energy dissipation via internal friction.

The variability of the lunar gravitational field in turn exerts perturbations on an orbiting spacecraft. The sensitivity of the orbiter motion to the lunar tides depends on the specific orbit configuration. The practical consequence of these perturbations is also a question of tracking capability and precision. Nowadays, the orbit determination of lunar spacecraft has achieved a customary decameter-level accuracy or better and requires elaborate force models of commensurate accuracy (Mazarico et al. 2012, 2018; Löcher and Kusche 2018). For instance, it is essential to account for the perturbations stemming from the elastic tidal response in gravity recovery campaigns (Konopliv et al. 2001; Goossens and Matsumoto 2008; Matsumoto et al. 2010; Lemoine et al. 2014; Yan et al. 2020). The viscoelastic tidal deformation of the Moon, well understood in the literature, is an even smaller effect compared with the elastic response. However, as the precision of spacecraft tracking improves, it may soon become necessary to take into account viscoelastic tidal effects in orbit modeling, which motivates the present study. Furthermore, the determination of the Moon’s tidal response, including phase lags at different frequencies, would have broad implication on our understanding of the lunar interior and evolution. Complementary to seismic measurements which are planned for future landed missions, the phase lags would in particular provide information on the low-viscosity zones and other dissipative layers in the deep interior (Briaud et al. 2023).

We remind ourselves of a rich literature over decades on the theory and determination of perturbed motions of Earth’s satellites, in which the role of (solid) Earth’s tidal deformation has been considered in detail and with rigor. Love numbers and phase lags at prominent frequencies have been independently derived from satellite tracking and altimetric data (Kozai 1967; Cazenave et al. 1977; Ray et al. 1996). Following such investigations as by Kaula (1964); Kozai (1965); Lambeck et al. (1974); Felsentreger et al. (1976), etc., we aim to assess the impact of lunar viscoelasticity in the motion of a polar orbiter, representing a typical class of survey satellites. We use the lunar tidal deformation model by Williams and Boggs (2015) (Sect. 2). We employ both numerical schemes and analytical expressions to investigate the long-term orbiter motion and as a means to validate model implementation (Sect. 3), followed by an interpretation of the orbital spectra (Sect. 4); finally, the errors in the orbit solutions due to omission of the viscoelastic effects are preliminarily evaluated (Sect. 5).

2 Lunar viscoelastic tidal deformation: model formulation

The external gravitational potential V of the Moon can be expressed as a spherical harmonic series in the body-fixed frame (Heiskanen and Moritz 1967),

$$\begin{aligned} V = \frac{\mu }{r} \displaystyle \sum ^{l_\textrm{max}}_{l=0} \sum ^l_{m=0} \left( \frac{R}{r} \right) ^l P_{lm}(\sin \varphi ) \left( C_{lm} \cos m\lambda + S_{lm} \sin m\lambda \right) \end{aligned}$$

(1)

where $\mu =GM$ and R are the gravitational constant and reference radius of the body, respectively, $r,\lambda ,\varphi$ are the spherical coordinates of the field point, $P_{lm}$ are the associated Legendre functions, and $C_{lm},S_{lm}$ are Stokes’ coefficients of the field variations at degree l and order m. Equation (1) is expressed in the lunar Principal Axes (PA) frame (Lemoine et al. 2014).

Referring to Williams and Boggs (2015), the time-varying components of the Stokes’ coefficients due to lunar tidal deformation at degree 2 can be expressed as a series expansion at designated frequencies associated with the motions of the Moon, Earth, and Sun,

$$\begin{aligned} \begin{aligned} \Delta C_{2m}(t)&= \displaystyle \sum _j \textrm{Re} \left\{ \textsf{k}_{2(j)} \ (-\textrm{i})^{\delta _{m1}} C_{2m(j)} e^{\textrm{i} \zeta _{(j)}} \right\} , \\ \Delta S_{2m}(t)&= \displaystyle \sum _j \textrm{Im} \left\{ \textsf{k}_{2(j)} \ \textrm{i}^{\,\delta _{m1}} S_{2m(j)} e^{\textrm{i} \zeta _{(j)}} \right\} . \end{aligned} \end{aligned}$$

(2)

The angular arguments, $\zeta _{(j)}(t)$, are linear combinations of the Delaunay arguments, whose linear change rates then specify the respective frequencies (Simon et al. 1994; Petit and Luzum 2010). The real constants $C_{2m(j)}, S_{2m(j)}$ measure the frequency-dependent contribution of the potential raised by Earth and Sun with respect to the lunar reference sphere of radius R, more specifically, through

$$\begin{aligned} \begin{aligned} \textrm{Re} \left\{ (-\textrm{i})^{\delta _{m1}} C_{2m(j)} e^{\textrm{i} \zeta _{(j)}}\right\}&= C_{2m(j)} \left[ \begin{matrix} \cos \zeta _{(j)} \\ \sin \zeta _{(j)} \end{matrix} \right] ^{m\ne 1}_{m=1} \\ \textrm{Im} \left\{ \ \textrm{i}^{\,\delta _{m1}} S_{2m(j)} e^{\textrm{i} \zeta _{(j)}} \right\}&= S_{2m(j)} \left[ \begin{matrix} \sin \zeta _{(j)} \\ \cos \zeta _{(j)} \end{matrix} \right] ^{m\ne 1}_{m=1} \end{aligned} \end{aligned}$$

(3)

The complex Love numbers $\textsf{k}_{2(j)}$, which describe the body’s response to the degree 2 tidal potentials, can be expressed as,

$$\begin{aligned} \textsf{k}_{2\,(j)} = \vert \textsf{k}_{2\,(j)} \vert e^{-\textrm{i} \Delta \zeta _{(j)}} = k_{2(j)} - \textrm{i} k^{\Im }_{2(j)} \end{aligned}$$

(4)

where the imaginary component, $k^\Im _{2(j)}$, causes a phase shift (lag) with respect to $\zeta _{(j)}$. It follows that $\Delta \zeta _{(j)} = \textrm{arctan}(k^\Im _{2(j)}/k_{2(j)})$. The phase shift is related to the quality factor of energy dissipation by $\textrm{sin} \vert \Delta \zeta _{(j)} \vert$ = $1/Q_{(j)}$ (Efroimsky 2012). For numerical evaluation, we rewrite Eq. (2) as follows,

$$\begin{aligned} \begin{aligned} \Delta C_{2m}(t) =&\displaystyle \sum _j \vert \textsf{k}_{2(j)} \vert C_{2m(j)} \left[ \begin{matrix} \cos \tilde{ \zeta }_{(j)} \\ \sin \tilde{ \zeta }_{(j)} \end{matrix} \right] ^{m\ne 1}_{m=1}, \\ \Delta S_{2m}(t) =&\displaystyle \sum _j \vert \textsf{k}_{2(j)} \vert S_{2m(j)} \left[ \begin{matrix} \sin \tilde{ \zeta }_{(j)} \\ \cos \tilde{ \zeta }_{(j)} \end{matrix} \right] ^{m\ne 1}_{m=1}, \end{aligned} \end{aligned}$$

(5)

where $\tilde{\zeta }_{(j)} = \zeta _{(j)} - \Delta \zeta _{(j)}$.

Table 1 is an excerpt from Williams and Boggs (2015) whose series approximation of Earth’s and solar tides consists of 21 frequencies with a truncation error of about 1% (Fig. 1a). We will refer to the frequencies (periods) near 27 days, namely, $j=1,2,3,14$ (and 20,21), as (lunar) monthly components. Those between 13 and 16 days ($j=4,5,6,7,11,18$) are the semimonthly components. Accordingly, we designate frequencies around 9 days ($j=9,12,13,16,19$) as tritmonthly, or thrice per month. In a similar way, the frequencies at 365.3 and 188.2 days are the annual and semiannual components ($j=10,15$), respectively.

If the Moon were purely elastic, its response would not only be in scale (i.e., with $k_{2(j)}=\vert \, \textsf{k}_{2(j)}\, \vert$) but also in phase with the external, tide-raising potential with $\zeta _{(j)}$ at all frequencies. A plausible set of Love numbers was derived based on the lunar interior simulation by (Briaud et al. 2023). The Love numbers were sampled at the same frequencies as in Table 1. Note that $k^\Im _{2\,(0)}=0$ or $\Delta \zeta _{(0)}=0$. We note that the tidal dissipation of a librating body, as is the case of the Moon, is a complicated problem that requires delicate treatment (Efroimsky 2018).

An illustration of the viscoelastic tidal response of the gravitational potential at the latitude of 45 degrees and 15 degrees in longitude on the lunar surface is given in Fig. 1b. The viscoelastic curve lags slightly behind the elastic one (panel b), due to $k^\Im _2$ being mostly between 30 and 40 times smaller than the real $k_2$, or $\Delta \zeta$ between 1 and 2 degrees, except at the end frequencies. Note that the phase lags are dependent on the interior model in the simulation. As a check, the spectrum of the tide-raising potential series exhibits precisely the 21 frequencies, as depicted in Fig. 1c. The gravitational potential and its spectrum at other locations, e.g., along a moving platform at $r(t),\lambda (t),\varphi (t)$, will be different. The formulation and evaluation of the tidal effects via Eq. (5) remain the same, which is subsequently implemented for the computation and analysis of the spacecraft orbit.

Table 1 Frequency-dependent constants of tide-raising potentials by Earth and Sun (excerpt from Williams and Boggs 2015) and modeled Love numbers of the Moon (Briaud et al. 2023)

Fig. 1

Temporal variability of degree 2 lunar gravitational field due to tidal deformation. a Potential at latitude of 45 deg and longitude of 15 deg on the lunar surface between the years 2010 and 2020. The black curve indicates truncation errors of the series approximation according to Williams and Boggs (2015) with respect to the exact expression. b Comparison of elastic and viscoelastic deformations over a period of 30 days, with $\textsf{k}_{2\,(j)}$ given in Table 1. c Spectrum of potential variability. Dashed vertical lines mark five periods belonging, respectively, to the tritmonthly, semimonthly, monthly, semiannual, and annual groups. The annotated periods are in days

3 Perturbed spacecraft motion due to viscoelastic lunar tides

The motion of a spacecraft is affected by and, in turn, can be used to estimate the body’s tidal response or potential Love numbers. On the one hand, the possibility of estimation depends on the sensitivity of the motion to perturbations that vary with orbit configurations. On the other hand, there comes the question if the induced motions are observable for given measurement types and precisions. This work addresses the former question, i.e., how the motion of a low-altitude polar orbiter is influenced by, and thus sensitive to, the viscoelastic effects expressed by Eq. (5).

3.1 Force model

The perturbed orbiter motion is evaluated numerically as follows. Let the state vector of the orbiter be

$$\begin{aligned} \textbf{X} = \begin{pmatrix} \textbf{r} \\ \dot{ \textbf{r} } \end{pmatrix}, \end{aligned}$$

(6)

where $\textbf{r}$, $\dot{ \textbf{r} }$ denote, respectively, the position and velocity in the Moon-centered, inertial frame. The state evolution over time is given by

$$\begin{aligned} \textbf{X}(t) = \textbf{X}(t_0) + \displaystyle \int ^{t}_{t_0} \dot{ \textbf{X} }(t,\textbf{X})\, \textrm{d}t \end{aligned}$$

(7)

with $\dot{ \textbf{X} }^\textrm{T} = ( \dot{ \textbf{r} }^\textrm{T}\ \ \ddot{ \textbf{r} }^\textrm{T} )$. If we assume the spacecraft motion is influenced only by the time-varying lunar gravitational field, the acceleration of an orbiting spacecraft can be expressed as

$$\begin{aligned} \ddot{ \textbf{r} } = (\textbf{R}^\textrm{B}_\mathrm {\,I})^\textrm{T} \frac{\partial V}{ \partial (\textbf{R}^\textrm{B}_\mathrm {\,I} \textbf{r} ) }, \end{aligned}$$

(8)

where $\textbf{R}^\textrm{B}_\mathrm {\,I}$ stands for the rotation matrix of coordinate system transformation, such that $\textbf{R}^\textrm{B}_\mathrm {\,I} \textbf{r}$ is expressed in the lunar body-fixed frame (in which V is defined). In the case of simple, uniform body rotation, $\textbf{R}^\textrm{B}_\textrm{I}$, becomes $\textbf{R}_z(\theta )$ where $\theta$ measures positive counterclockwise the angle of rotation with respect to the z-axis.

Of main interest is the perturbed orbit motion due to viscoelastic tidal response of the Moon with respect to the “nominal” case in which viscoelasticity, or the phase lag of deformation, is neglected, e.g.,

$$\begin{aligned} \Delta \textbf{X}_{(j)} = \textbf{X}_{(j)} - \textbf{X}. \end{aligned}$$

(9)

In the following, the “(j)” subscripted quantities are associated with the viscoelastic case, i.e., in the presence of $k_{2(j)}$ or $\Delta \zeta _{(j)}$ as in Eq. (4) with a consistent notation. The indices j are as given in Table 1. We may treat and denote accordingly more than one component, for instance, “$(j_1,j_2,\cdots )$” indicates some particular frequencies at $j_1,j_2,\cdots$, and “$(\sum j)$” includes all 21 frequencies in Table 1 (Williams and Boggs 2015). For matrix operations, multiple $k_{2\,(j)}$ are arranged into a column vector, e.g., $\textbf{k}_{2\,(\sum j)} = ( k_{2\,(1)}\ \cdots \, k_{2\,(21)})^\textrm{T}$.

$\Delta \textbf{X}_{(j)}$ can be evaluated via a direct comparison of simulated orbits with and without the viscoelastic effects modeled by $k_{2\,(j)}$. It can be approximated to the first order of $k_{2(j)}$ as,

$$\begin{aligned} \Delta \textbf{X}_{(j)}(t) \approx \left( \displaystyle \int ^t_{t_0} \frac{ \partial \dot{ \textbf{X} }_{(j)} }{ \partial \textbf{k}_{2\,(j)} } \textrm{d}t \right) \textbf{k}_{2\,(j)} = \Phi _{(j)} \textbf{k}_{2\,(j)}, \end{aligned}$$

(10)

where the matrix $\Phi _{(j)}$ evaluates the sensitivity of the orbit (change) to $\textbf{k}_{2\,(j)}$ from $t_0$ through t, and $\textbf{k}_{2\,(j)}$ can be a column vector or a scalar.

3.2 Coordinate systems and reference frames

The numerical orbit propagation is evaluated by Eq. (7) in the Earth Mean Equator and Equinox of 2000 (EME2000) inertial frame in terms of Cartesian coordinates. The transformation matrix (i.e., $\textbf{R}^\textrm{B}_\textrm{I}$ in Eq. 8) between the lunar PA body-fixed frame and the EME2000 is retrieved from the DE421 ephemeris via SPICE of NASA’s Navigation and Ancillary Information Facility (Folkner et al. 2008).

The orbital elements, i.e., the semimajor axis a, the eccentricity e, the inclination i, the longitude of the ascending node $\Omega$, the argument of the periapsis $\omega$, and the mean (or true) anomaly M (or f), are more convenient for describing and interpreting the orbital behavior. They are referred to a Moon-centered, equatorial frame (MEq).¹ In order for the frame to be inertial, we let it coincide with the lunar PA frame at a given epoch, say $t_0$, and assume it is nonrotating so that the transformation matrix into and from the EME2000, $\textbf{R}^\textrm{B}_\textrm{I}(t_0)$, is constant. The epoch $t_0$ is selected as UTC 2:00:00 on April 5, 2014.

Table 2 Orbit parameters of LRO in the MEq frame on 2014 Apr. 5 at 2:00:00 (UTC)

3.3 The case of the Lunar Reconnaissance Orbiter

The analysis hereafter is based on the orbit of the Lunar Reconnaissance Orbiter (LRO) from NASA (Mazarico et al. 2018). Over a majority of its mission phase, the LRO is operating in a low, near-circular, polar orbit. The orbital elements at $t_0$ in the MEq frame are listed in Table 2.

We make use of the software Technical University Delft Astrodynamics Toolbox (TUDAT) for orbit propagation, taking advantage of, among others, its generic, multi-body simulation setting, which is particularly suited for modeling gravitational interactions (Dirkx et al. 2019).² The tidal deformation is incorporated by expanding the existing module based on the elastic case to include amplitudes and phase lags at arbitrary frequencies, i.e., per user specification.

3.3.1 Spacecraft tidal acceleration and comparison with other perturbations

In addition to nonspherical primary gravitation, other sources of perturbations on a free-orbiting spacecraft include solar radiation pressure and gravity of external bodies. A comparison of the perturbations acting on the LRO over two orbits is shown in Fig. 2. Thermal radiation and sunlight reflected off the lunar surface may need to be considered (Smith et al. 2008; Löcher and Kusche 2018); these effects are not assessed here. While elastic tidal deformation plays a significant role in affecting spacecraft dynamics, the viscoelastic effects are at the low end of all considered mechanisms but might become relevant for high-precision applications.

Fig. 2

Comparison of perturbations on LRO over two orbit periods. The viscoelastic effects are the weakest among the presented mechanisms. The static gravitational field is up to degree 300. Drop-offs of the solar radiation pressure occur when the spacecraft is entering eclipse phases

3.3.2 Orbit difference

Figure 3 shows $\Delta \textbf{X}_{(\sum j)}$, the difference between orbits with and without considering the viscoelastic perturbations, evolving from the same initial condition (Table 2). The duration of 4 days roughly corresponds to the single-arc length in orbit determination. Simulations were performed for two contrasting resolutions of the gravitational field model, $l_\textrm{max}=2$ or 300. The former corresponds to the dominant contribution of the tidal response at degree 2. The latter is roughly the minimum resolution needed for precise orbit determination of the LRO (Mazarico et al. 2018).

The impact is least significant in the radial direction or on the orbit size, as is typical with conservative forcing in a near-circular orbit. The oscillation at the short, orbital period of about 2 h for the LRO amplifies over time but does not exceed 2 cm in magnitude. The transverse component exhibits a systematic trend that grows steadily to reach about 0.5 m after four days. The normal component is, as the radial, dominated by the amplifying, oscillatory pattern with no clear systematic trend over time.

The nonshort-periodic, or long-term, variations are more clearly revealed in the angular arguments, whereas the semimajor axis and eccentricity exhibit far less or no long-term variations (Kaula 1966, see also Sect. 3.4). The inclination, longitude of the ascending node, and argument of latitude, $\omega + f$, all display periodicities on the order of days (Fig. 4). The differences correspond to those observed in Fig. 3 in magnitude, i.e., at the level of $R \Delta i \approx 0.1$ m. The stronger deviation in $\omega + f$ is associated with the transverse acceleration under $k^\Im _{(\sum j)}$(Fig. 3). The orbit differences are only slightly influenced by the resolution of the gravitational field.

Fig. 3

Orbit difference $\Delta \textbf{X}_{(\sum j)}$ in the radial (R), transverse (T), and normal (N) directions over four days in the case of LRO. The solid and the dot-dashed curves correspond, respectively, to the gravitational field resolutions of degree 2 and 300

Fig. 4

Evolution of inclination (top), longitude of the ascending node (middle), and argument of latitude (bottom) due to $k^\Im _{(\sum j)}$. The solid and the dot-dashed curves correspond, respectively, to the gravitational field resolutions of degree 2 and 300

3.3.3 Note on model simplification, limitation

The orbit difference in Figs. 3 and 4 no greater than 50 cm can be easily overwhelmed by other mechanisms, some of which are illustrated in Fig. 2. Another phenomenon of note is the physical librations of the Moon, which concern the orientation of the body (frame) with respect to the inertial space and hence the gravitational perturbation acting on the LRO. The monthly libration amplitudes alone (i.e., at $\mathcal {F}$, $\ell$) reach about 100 arcseconds in both latitude and longitude (Rambaux and Williams 2011). Figure 5a shows the evolution of the three Euler angles for the transformation from the EME2000 into the lunar PA frame relative to the case of uniform rotation from the same initial epoch (Table 2) (Folkner et al. 2008; Archinal et al. 2018).³ The gravitation of a librating body differs from that of a uniformly rotating one by as large as $10^{-5}$ $\mathrm {m^2 s^{-2}}$ (Fig. 5b). Assuming a simple, uniform lunar rotation thus incurs an error already larger than all but lunar gravitation itself (Fig. 2).

In practice, e.g., precise orbit computation and determination, such a simple model would be of no use. Nonetheless, it is legitimate to simplify the dynamic model in the present analysis, which focuses solely on the role of viscoelasticity. Because the perturbation is extremely small, the induced orbit variation relative to some (viscoelastically) unperturbed reference orbit depends little on the latter. Namely, the accuracy of the reference orbit, or that of the propagating model, does not affect the evaluated orbit difference, at least insofar as the first-order effect is concerned. Figures 3 and 4 show, for instance, that the viscoelastic signatures are quite similar for (static) gravitational fields up to degree 2 and degree 300. The same is true for the impact of physical librations. Figure 5c shows the discrepancies between the evaluated angular arguments, $\Delta i$, $\Delta \Omega$, and $\Delta (\omega +f)$, with or without librations, which reach only $10^{-5}$ of the respective variations shown in Fig. 4. Hence, the curves coincide with those in Fig. 4, irrespective of the rotation model.

Consequently, we can employ analytic approaches taking into account only low-degree static gravitational field and uniform rotation of the Moon. It should be reemphasized, however, that this is only feasible when analyzing the viscoelastic effects in isolation. For accurate orbit modeling, care is needed to ensure not only completeness of perturbation mechanisms but also consistency of the coordinate frames.

Fig. 5

Impact of lunar physical librations. a Evolution of Euler angles for the coordinate transformation from EME2000 frame into lunar Principal Axes frame with respect to the case of uniform body rotation. b Discrepancy of the gravitational perturbations in the two frames. c Discrepancy of the evaluated $\Delta i$, $\Delta \Omega$, and $\Delta (\omega +f)$ using the two frames

3.4 Frequency analysis and long-term effects

Further insight into the orbiter motion can be gained by an inspection of the orbital elements in the frequency domain, due to distinct periodicities of the tide-raising potentials and the induced temporal variability of the gravitational field. The purpose is to average out the short-periodic variations associated with M so that the long-term trend is distinguished by the contributions from the tidal components and the slow-varying arguments of $i,\Omega ,\omega$. The analytic approach below has been routinely applied in previous investigations, above all, on satellite orbit subject to Earth’s solid and ocean tides (Kaula 1964; Kozai 1965; Lambeck et al. 1974). The analysis is likewise instructive in the case of lunar orbiters. The approach is theoretically justified, since the error accumulation in the linear approximation of orbit variation is tolerable over the short timespan of days. The comparison of the numerically propagated orbits with the analytic results offers a means to validate the implementation of the tidal deformation model.

The gravitational potential of degree l and order m is expressed in terms of orbit elements as (Kaula 1966),

$$\begin{aligned} V_{lm} = \frac{ \mu }{ a } \left( \frac{R}{a} \right) ^l \displaystyle \sum _{p=0}^{l} \sum _{q=-\infty }^{\infty } F_{lmp}(i)\,G_{lpq}(e)\, S_{lmpq}( \omega , M, \Omega -\theta ), \end{aligned}$$

(11)

with

$$\begin{aligned} \begin{aligned} S_{lmpq} = \left[ \begin{array}{l} \textrm{Re} \\ \textrm{Im} \end{array} \right] ^{l-m\ \textrm{even}}_{l-m\ \textrm{odd}}&\Big \{ \left( C_{lm} - \textrm{i} S_{lm} \right) \\&\times \exp \textrm{i}\Big [ (l-2p)\omega + (l-2p+q) M + m(\Omega -\theta ) \Big ] \Big \}. \end{aligned} \end{aligned}$$

(12)

$\theta$ measures the rotation angle since $t_0$ in the MEq frame, and $\dot{\theta }=2\pi /27.322$ $\mathrm {rad/day}$ is assumed to be constant.

The second-degree Stokes’ coefficients are time-varying according to Eq. (5). We augment the notation of the periodic function as $S_{2mpq}(\omega , M, \Omega -\theta , \tilde{\zeta }_{(j)})$ to explicate the variability through the dependence on $\tilde{\zeta }_{(j)}=\zeta _{(j)}-\Delta \zeta _{(j)}$. Long-term variations are associated with $l-2p+q=0$. The eccentricity function, $G_{2\,p\,(2p-2)}$, only exists for $p=1$, whereby the argument $\omega$ vanishes. The disturbing potential is therefore simplified as,

$$\begin{aligned} \Delta V_{2m\,(j)} = \frac{ n^2 R^2 }{ (1-e^2)^{3/2} }\, F_{2m1}(i)\, S_{2m10}( \Omega , \theta , \tilde{\zeta }_{(j)} ), \end{aligned}$$

(13)

in which the inclination function is one of the following,

$$\begin{aligned} F_{201} = \frac{3}{4} \sin ^2i - \frac{1}{2}, \quad F_{211} =-\frac{3}{2} \sin i \cos i, \quad F_{221} = \frac{3}{2} \sin ^2 i. \end{aligned}$$

(14)

The equations of perturbed orbital motion are (dropping the subscripts) (Brouwer and Clemence 1961)

$$\begin{aligned} \begin{aligned} \frac{ \textrm{d} (\Delta i) }{ \textrm{d} t }&=-\frac{ 1 }{ n a^2 (1-e^2)^{1/2} \sin i } \frac{ \partial (\Delta V) }{ \partial \Omega }, \\ \frac{ \textrm{d} (\Delta \Omega ) }{ \textrm{d} t }&= \frac{ 1 }{ n a^2 (1-e^2)^{1/2} \sin i } \frac{ \partial (\Delta V) }{ \partial i }, \\ \frac{ \textrm{d} (\Delta \omega ) }{ \textrm{d} t }&=-\frac{ \cos i }{ n a^2 (1-e^2)^{1/2} \sin i } \frac{ \partial (\Delta V) }{ \partial i } + \frac{ (1-e^2)^{1/2} }{ n a^2 e } \frac{ \partial (\Delta V) }{ \partial e }, \\ \frac{ \textrm{d} (\Delta M) }{ \textrm{d} t }&= -\frac{ 1-e^2 }{ n a^2 e } \frac{ \partial (\Delta V) }{ \partial e } - \frac{ 2 }{ na } \frac{ \partial (\Delta V) }{ \partial a }. \end{aligned} \end{aligned}$$

(15)

Equation (13) indicates that the disturbing potentials and induced perturbations are distinguished by the order m. Consequently, the components associated with tesserals and nontesserals in Table 1 should influence the spacecraft orbit motion differently.

3.4.1 Inclination

Although tidal deformation gives rise to a time variability of the zonal coefficient $\Delta C_{20}(t)$, the inclination is not directly influenced by the zonal harmonics since $\partial \Delta V_{l0}/\partial \Omega =0$. The first expression of Eq. (15) can be rewritten as

$$\begin{aligned} \begin{aligned} \frac{ \textrm{d} (\Delta i ) }{\textrm{d}t}&= -\frac{3n R^2}{4 a^2 (1-e^2)^2} \times \\&\quad \Big [\ \cos i \sum _j \sum ^-_+ k^\Im _{2\,(j)} (C_{21\,(j)} \pm S_{21\,(j)})\cos (\Omega - \theta \pm \zeta _{(j)} ) \\&\quad + 2 \sin i \sum _j \sum ^-_+ k^\Im _{2\,(j)} (-S_{22\,(j)} \pm C_{22\,(j)})\cos (2\Omega - 2\theta \pm \zeta _{(j)} ) \ \Big ], \end{aligned} \end{aligned}$$

(16)

where $\sum _{+}^{-} f(x_1 \pm y_1, \cdots , x_i \pm y_i, \cdots ) = f(x_1 + y_1, \cdots , x_i + y_i, \cdots ) + f(x_1 - y_1, \cdots , x_i - y_i, \cdots )$ as in Lambeck et al. (1974). For a linear solution, the sinusoidal functions $\cos x$ in Eq. (16) can be integrated directly as $\dot{x}^{-1} \sin x$. The change rates, $\dot{\Omega }$, $\dot{\theta }$, and $\dot{\zeta }_{(j)}$ are treated as constants. While they may assume small values, the divisors, $m(\dot{\Omega }-\dot{\theta }) \pm \dot{\zeta }_{(j)}$, are never zero so that there exists no strictly secular variation.

Meanwhile, Eq. (14) indicates $F_{221}(i) \gg F_{211}(i)$ for a polar orbiter, such as the LRO, since $\cos (i\approx \pi /2) \approx 0$. Hence, the trend of $\Delta i$ is governed by the sectoral component. The tide-raising potential and the body response are the strongest at the monthly period(s). As a result, $\Delta i$ is associated foremost with $j=2$:

$$\begin{aligned} \Delta i_{(2)}= & {} \, -\frac{3n R^2}{2a^2 (1-e^2)^2}\,\sin i \nonumber \\{} & {} \times \sum ^{-}_{+} k^\Im _{2\,(2)} (-S_{22\,(2)} \pm C_{22\,(2)}) \frac{\sin (2\Omega - 2\theta \pm \zeta _{(2)}) }{2\dot{\Omega }-2\dot{\theta } \pm \dot{\zeta }_{(2)}} - i_0, \end{aligned}$$

(17)

where $i_0$ is a constant to fulfill the initial condition $\Delta i=0$. Because $\dot{\zeta }_{(2)} \approx \dot{\theta }$ with a deviation of about 1%, $\Delta i_{(2)}$ presents two distinct frequencies at $2\dot{\theta } \pm \dot{\zeta }_{(2)} \approx 3\dot{\theta }$ and $\dot{\theta }$, so long as $\dot{\Omega } \ll \dot{ \theta }$ (as is the case for the LRO).

Figure 6 shows $\Delta i_{(2)}$ in comparison with that under $k^\Im _{2\,(\sum j)}$ as also presented in Fig. 4. $\Delta i_{(2)}$ is a main constituent of the total amplitude as well as periodicity for the variation. Equation (17) offers a reasonable approximation and explanation for the long-term trend displayed in Fig. 6. The amplitude is obviously in proportion to $k^\Im _{2\,(2)}$. One evident periodicity is the tritmonthly, associated with $1/3\dot{\theta } \approx 9$ days. The curves comprise also a monthly component of $1/\dot{ \theta }$, which is only incompletely depicted and thus resembles a linear trend over 15 days.

Fig. 6

Evolution of orbit inclination over 15 days. The solid gray curve indicates the impact of $k^\Im _{(\sum j)}$ as displayed also in the top panel of Fig. 4. The dotted red curve is for a single component of $k^\Im _{(2)}$. The dashed black curve approximates analytically the long-term trend for the latter via Eq. (17). The analytic expression and numerical result are compared more closely over two and half orbit periods at the beginning (lower-left embedded panel)

3.4.2 Ascending node

The precession of the orbital plane differs from the behavior of the inclination. According to Eq. (15), we write the precessing rate of the ascending node as

$$\begin{aligned} \frac{\textrm{d} (\Delta \Omega ) }{ \textrm{d} t}= & {} \frac{n R^2}{2a^2 (1-e^2)^2} \sum _{m=0}^{2} \sum _j \sum ^-_+ \frac{ F^\prime _{2m1}(i) }{ \sin i }\, k^\Im _{2\,(j)} \nonumber \\{} & {} \left[ \begin{matrix} -S_{2m\,(j)} \pm C_{2m\,(j)} \\ -( C_{2m\,(j)} \pm S_{2m\,(j)}) \end{matrix} \right] ^{m\ne 1}_{m=1} \,\sin \big ( m(\Omega -\theta ) \pm \zeta _{(j)} \big ), \end{aligned}$$

(18)

where the outermost summation applies through all orders. $S_{20\,(j)}=0$ is present only for formality or symmetry.

Due to the absence of the term $\cos i$, $F'_{211}$ outweighs $F'_{201}$ and $F'_{221}$. $\Delta \Omega$ is predominantly induced by the monthly tidal response, i.e., as is the inclination but at $j=1$ and a slightly shorter period of 27.2 days (Table 1). Assuming $\dot{\zeta }_{(1)}=\dot{\theta }$ and $\dot{\Omega } = 0$, we integrate the two periodic components of $m=1$ separately (not by reason of singularity, as Eq. 18 can be integrated directly),

$$\begin{aligned} \begin{aligned} \Delta \Omega _{(1)}&=\,-\frac{3 n R^2}{4a^2 (1-e^2)^2}\, k^\Im _{2\,(1)} C_{21\,(1)}\, \times \\&\qquad \left[ \, \sin ( \Omega -\theta +\zeta _{(1)}) \Delta t - \frac{\cos ( \Omega -\theta -\zeta _{(1)})}{ 2 \dot{\theta } } \,\right] - \Omega _0, \end{aligned} \end{aligned}$$

(19)

where the contribution from $S_{21\,(1)}$($\ll C_{21\,(1)}$) is neglected and where $\Omega _0$ takes such a value that $\Delta \Omega =0$ initially.

As shown in Fig. 7a, $\Delta \Omega$ exhibits a linear trend superimposed by a semimonthly variation governed by $\Delta \Omega _{(1)}$. The contribution from other frequencies is perceptible, e.g., yielding a better approximation of the numerical result, but secondary. A closer look shows that the analytic approximation deviates from the numerical result in an apparent secular trend but does not grow substantially (gray curve in Fig. 7b). We infer that this is due to the errors of the linear approximation, i.e., orbit changes evaluated with respect to fixed reference orbit elements.

Fig. 7

Evolution of longitude of right ascension. a Impact of $k^\Im _{2\,(j)}$ evaluated by numerical integration (solid gray curve) and long-term approximation (dashed black and dotted red curves). b Indirect effect caused by $\Delta i_{(2)}$ marked by dashed red curve accounts for the difference between numerical and long-term approximation (solid gray curve)

3.4.3 Indirect effect

One reason for the discrepancies between the numerical integration and analytical approximation in Figs. 6 and 7a is that the spacecraft in two diverging orbits, e.g., those separated by $\Delta i$ or $\Delta \Omega$, will experience increasingly different perturbations. The perturbations, in turn, widen the separation of the orbits (Lambeck et al. 1974). Such indirect effects are, above all, associated with the second-degree (static) field coefficients, $C_{20}$ and $C_{22}$. However, for low-altitude lunar orbiters, i.e., $a \approx R$, the contribution of higher degrees may still be significant (Felsentreger et al. 1976). The mechanism is implicitly addressed by numerical integration of differential orbit changes and is therefore of no concern in practice.

For a diagnosis of the numerical results, the magnitude of the contribution to the total orbit difference is roughly assessed as follows. The perturbations due to $C_{20}, C_{22}$ cause the ascending node of the spacecraft orbit to precess at a rate of $\dot{ \Omega } = \frac{3n R^2}{ 2a^2 (1-e^2)^2 } \left[ C_{20} + 2 C_{22} \cos (2\Omega -2\theta ) \right] \cos i$. Leaving out $R/a \approx 1$ and $(1-e^2)^{-2}$ and differentiating the expression with respect to i and $\Omega$, we approximate the additional variation caused by $\Delta i,\, \Delta \Omega$ as,

$$\begin{aligned} \begin{aligned} \frac{ \textrm{d} (\Delta \Omega ^+) }{ \textrm{d}t } = -&\frac{3n }{2} \Big \{ \left[ C_{20} + 2 C_{22} \cos (2\Omega -2\theta ) \right] \sin i \Delta i \\&+ 4C_{22} \cos i\, \sin (2\Omega -2\theta ) \Delta \Omega \Big \}. \end{aligned} \end{aligned}$$

(20)

Since the main component of $\Delta i$ comes from $\Delta i_{(2)}$ at the frequencies $\dot{ \theta }$ and $3\dot{ \theta }$, the induced $\Delta \Omega ^+$ will not grow secularly, for $\int \Delta i_{(2)} \textrm{d}t$ and $\int \cos 2(\Omega -\theta ) \Delta i_{(2)} \textrm{d}t$ remain periodic. It suffices to consider only $C_{20}\approx -2\times 10^{-4}$ to estimate the oscillatory amplitude. According to Eq. (17), the tritmonthly component is $\frac{3}{2} n k^\Im _{2\,(2)} (S_{22\,(2)}+C_{22\,(2)})/3\dot{\theta } \approx 4 \times 10^{-6}$ degrees (see also Fig. 6). The corresponding $\Delta \Omega ^+$ is $\frac{3}{2} n C_{20}/3\dot{\theta } \cdot (4 \times 10^{-6}) \approx 1.5 \times 10^{-7}$ degrees; the magnitude of the monthly component is similar at about $2 \times 10^{-7}$ degrees. Both correspond to the difference between the numerical result and the long-term approximation observed in Fig. 7.

The contribution of $\Delta \Omega$ to $\Delta \Omega ^+$ is weaker as a result of the small factor $\cos i$ and $C_{22}$ is about one-tenth of $C_{20}$ in magnitude. In the long(er) term, the main variation, $\Delta \Omega _{(1)}$, composed of secular and $2\dot{\theta }$ trends seems more likely to become significant. According to Eq. (19), the (absolute) linear rate is $\alpha = \frac{3}{4} n k^\Im _{2\,(1)} C_{21\,(1)} \approx 2 \times 10^{-6}$ degrees/day; the semimonthly amplitude is $\alpha /2 \dot{\theta } \approx 4 \times 10^{-6}$ degrees. Then, the variation induced by $\Delta \Omega _{(1)}$ in Eq. (20), namely, $6 n C_{22}\cos i\,\int \sin 2(\Omega - \theta ) \Delta \Omega _{(1)}\, \textrm{d}t$ produces two amplifications in scale with $(\alpha /2\dot{\theta }) \Delta t \cos 2(\Omega -\theta )$ and $(\alpha /2\dot{\theta }) \Delta t$, respectively. Both are on the order of $10^{-9}$ degrees/day and far slower than $\Delta \Omega _{(1)}$. Thus, $\Delta \Omega$ will be dominated by $\Delta \Omega _{(1)}$, and the indirect contribution mainly comes from $\Delta i_{(2)}$ within about 100 days. Figure 7b confirms that the discrepancy between the numerical and the long-term approximation is due to the presence of $\Delta \Omega ^+$. The apparent linear trend belongs in part to a monthly variation induced indirectly by $\Delta i_{(2)}$.

At last, we briefly assess the indirect effect of $\Delta \Omega$ on $\Delta i$. Following the same procedure above, we find

$$\begin{aligned} \frac{ \textrm{d} (\Delta i^+)}{ \textrm{d} t} = 3 n \sin iC_{22} \sin 2(\Omega -\theta ) \Delta \Omega . \end{aligned}$$

(21)

The quasi-secular and semimonthly components of $\Delta \Omega _{(1)}$ will cross-generate a periodic (i.e., semimonthly) amplification and a secular trend in $\Delta i^+$. The magnitude of the former varies with $3n \sin i\, C_{22}\, (\alpha /2\dot{\theta }) \Delta t \approx 2 \times 10^{-8}$ degrees/day. The secular variation is of the same order of magnitude.

4 Distinction of frequencies and practical implication

The viscoelastic signature is most distinct around 27 days at which the excitation is most prominent, which has been revealed by the analysis of the Gravity Recovery and Interior Laboratory (GRAIL) and Lunar Laser Ranging (LLR) observations (Williams et al. 2014). The variations at other frequencies will provide further constraints for interior modeling and data analysis (Williams and Boggs 2015). With the tidal deformation model implemented and the basic perturbation mechanism identified, we next examine how the tidal frequencies are reflected and their distinction in the spectra of the orbiter motion.

Figure 8 shows the spectra of $\Delta i$ (panel a) and $\Delta \Omega$ (panel b) of the LRO. The orbit spectra were obtained by integrating Eq. (15) over three years. The potential has signals only from about 9 days and beyond in period, as indicated in Table 1 and depicted in Fig. 1. Orbital variations consist additionally of short-period components that are multiples of the orbit frequency. These signals do not exceed 0.1 days and are scattered toward the lower end of the periodicity.

4.1 Monthly frequencies and multiples

The strong amplitudes around the monthly frequencies and those in the vicinity of their multiples result from the prominence of similar frequencies in the tide-raising potentials, e.g., at $j \le 7$, likewise clustered about the monthly and semimonthly (Fig. 1c). As discussed, the distinguished amplitudes of orbital variations do not correspond one-to-one to those of tidal deformation.

4.1.1 Inclination

Only temporal variability of the sectoral harmonics affects inclination. The two peaks near 27 days both arise from $2\dot{ \theta } - \dot{ \zeta }_{(j)}$ ($\dot{ \Omega }$ omitted hereafter), namely,

$$\begin{aligned} \begin{aligned} 2\dot{ \theta } - \dot{ \zeta }_{(2)}&= \left( \frac{2}{27.322} - \frac{1}{27.555} \right) \quad \textrm{day}^{-1} = \frac{1}{27.091} \quad \textrm{day}^{-1}, \\ 2\dot{ \theta } - \dot{ \zeta }_{(3)}&= \left( \frac{2}{27.322} - \frac{1}{31.812} \right) \quad \textrm{day}^{-1} = \frac{1}{23.943} \quad \textrm{day}^{-1}, \end{aligned} \end{aligned}$$

(22)

with counterparts mirrored to around 9 days, e.g.,

$$\begin{aligned} \begin{aligned} 2\dot{ \theta } + \dot{ \zeta }_{(2)}&= \left( \frac{2}{27.322} + \frac{1}{27.555} \right) \quad \textrm{day}^{-1} = \frac{1}{9.133} \quad \textrm{day}^{-1}, \\ 2\dot{ \theta } + \dot{ \zeta }_{(3)}&= \left( \frac{2}{27.322} + \frac{1}{31.812} \right) \quad \textrm{day}^{-1} = \frac{1}{9.557} \quad \textrm{day}^{-1}. \end{aligned} \end{aligned}$$

(23)

Similarly, the neighboring groups at 7 and 5 days are associated with $2\dot{ \theta } + \dot{ \zeta }_{(j)}=1/7.096$ and 1/5.643 $\textrm{day}^{-1}$ for $j=4$ and 9, respectively (and in their spectral proximity).

The strongest peak of inclination change at 13.7 days differs from those above in that it indicates a distinct indirect effect, i.e., $\Delta i^+$ as in Eq. (21) (Fig. 8a). Note that $2\dot{ \theta } \pm \dot{ \zeta }_{(j)}$ invariably shifts the frequency away from the semimonthly. Only for (very) slow-varying components, e.g., $\dot{ \zeta }_{(10)}=$1/365 $\textrm{day}^{-1}$, will the separation be small and the resulting frequency remains near 13 days. Therefore, the semimonthly variation of the inclination must be indirectly induced by the secular $\Delta \Omega _{(1)}$. It has been shown that the corresponding semimonthly variation amplifies at a rate of about $2 \times 10^{-8}$ degrees/day. The magnitude of $10^{-5}$ degrees after 1000 days thus explains the semimonthly amplitude in Fig. 8a.

4.1.2 Ascending node

The spectral pattern of $\Delta \Omega$ as shown in Fig. 8b is explainable with the same approach. The dominant semimonthly component of $\Delta \Omega _{(1)}$ accounts for the strongest peak around 13.7 days. The impact of the sectoral harmonics is reduced by the factor $\cos i$. This is illustrated by the comparison of $\Delta \Omega$ and $\Delta i$ at 9.13 days, i.e., $2\dot{ \theta }-\dot{ \zeta }_{(2)}$, where the former is smaller by an order of magnitude.

Unlike inclination, $\Delta \Omega$ is also affected by the zonal harmonic, in which case the frequencies at which $C_{20\,(j)}$ are nonzero are directly introduced in the orbit motion. This concerns, for instance, the amplitudes around 27 days, of which the strongest is due to $C_{20\,(2)}=-605.4 \times 10^{-9}$. Though the impact of the zonal harmonic is reduced by the factor $\cos i$ as of the sectoral harmonics, the magnitude of $C_{20\,(2)}$ still produces a strong signal. Another telling observation about the role of the zonal harmonics is the peak at 31.8 days introduced via $C_{20\,(3)}$ at its original frequency, which is absent in inclination (Fig. 8b).

4.2 Ambiguity at semiannual and indistinctness of annual frequencies

A somewhat deceptive feature in the orbital spectra is the semiannual peak(s) at 182.7 days, which is the lowest identifiable frequency (Fig. 8). It can be shown, however, that it cannot be caused by the semiannual component of the tidal potential (at $j=15$), of which the tesseral harmonics are nonzero. By definition, $\Delta i$ does not admit the semiannual frequency via the zonal harmonic (for $\partial V_{l0}/\partial \Omega = 0$). On the other hand, referring to Eq. (16), we see that the resulting orbital frequency in inclination via the time-varying tesseral harmonics is shifted from the monthly as $\dot{\theta } \pm \zeta _{(15)} \approx \frac{7}{6} \dot{\theta }$ and $\frac{5}{6} \dot{\theta }$, both of which still lie close to $\dot{\theta }$. Therefore, $\zeta _{(15)}$ is not responsible for the semiannual peak.

Subsequently, we can rule out the possibility that it is related to the (double) annual frequency of 1/365.3 $\textrm{day}^{-1}$ at $j=10$, which affects the zonal and sectoral harmonics. The sectoral does not affect annual orbital variations in inclination for the same reason as with tesseral harmonics, since the resulting orbit frequencies are shifted to $2\dot{\theta } \pm \dot{ \zeta }_{(10)} \approx 2\dot{\theta }$. The zonal harmonic does produce annual variations of $\pm \dot{ \zeta }_{(10)}= \pm \dot{\theta }/13.4$, yet this does not result automatically in semiannual signals. The only viable mechanism to double the frequency seems to be through a second-order effect, namely, $\cos ^2 \zeta = \frac{1}{2}( 1+\cos 2\zeta )$. If this were the case, Fig. 8b would have exhibited a prominent annual peak (i.e., of the first-order), as well. The absence of the annual signal is not surprising, however, as $C_{20\,(10)}$=$-0.1\times 10^{-9}$, far smaller than $S_{22\,(10)} =-13.6\times 10^{-9}$ and four orders of magnitude below the monthly coefficients ($j=1,2$, Table 1).

Therefore, the last strong, semiannual peak in Fig. 8a is induced by neither semiannual ($j=15$) nor annual ($j=10$) tides. Instead, a scrutiny of Table 1 led to the following condition related to a semimonthly term $j=4$:

$$\begin{aligned} 2 \dot{\theta } - \dot{ \zeta }_{(4)} = \left( \frac{2}{27.322} - \frac{1}{14.765} \right) \quad \textrm{day}^{-1} \approx \frac{1}{182.7} \quad \textrm{day}^{-1}. \end{aligned}$$

(24)

We refer to Eq. (17) and find the amplitude of the “composite” semiannual variation of the inclination as

$$\begin{aligned} \left| \Delta i_{(4)} \right| _\mathrm {182.7~d} = \frac{3n}{2( 2 \dot{\theta } - \dot{ \zeta }_{(4)} )} k^\Im _{2\,(4)} \left| C_{22\,(4)}-S_{22\,(4)} \right| \approx 2\times 10^{-6} \quad \textrm{degrees}. \end{aligned}$$

(25)

The peak at 182.7 days with a magnitude of $1.47 \times 10^{-6}$ degrees in Fig. 8b is thus fully attributable to the quasi-cancellation of the semimonthly tidal deformation by the rotation rate of the Moon.

In comparison, the semiannual peak of $\Delta \Omega$ has another component (dashed line in Fig. 8b). The direct constituent, i.e., of $\Delta \Omega _{(4)}$, is small because the sectoral harmonics are factored by $\cos i$. The dominant contribution comes indirectly from $\left| \Delta i_{(4)} \right| _\mathrm {182.7~d}$, i.e., via static $C_{20}$ in Eq. (20)

$$\begin{aligned} \left| \Delta \Omega ^+ \right| _\mathrm {182.7~d} = \frac{3n}{2(2\dot{ \theta }-\dot{ \zeta }_{(4)})} C_{20} \left| \Delta i_{(4)} \right| _\mathrm {182.7~d} \approx 1\times 10^{-6}\ \ \textrm{degree}. \end{aligned}$$

(26)

The most significant direct contribution is rather from $j=14$ and slightly offset at $1/177.85=1/27.322-1/32.281$ $\textrm{day}^{-1}$:

$$\begin{aligned} \left| \Delta \Omega _{(14)} \right| _\mathrm {177.85~d} = \frac{3n}{4(\dot{ \theta } - \dot{ \zeta }_{(14)})} k^\Im _{2\,(14)} \left| C_{21\,(14)} + S_{21\,(14)}\right| \approx 6 \times 10^{-7}\ \ \textrm{degree}. \end{aligned}$$

(27)

Such coupling occurs with other nearby semimonthly components, all producing even lower frequencies than the semiannual. For instance, in the same mechanism as (24), $\zeta _{(5)}$ with a period of 13.777 days accounts for a slow orbital change of about 1620 days.

Fig. 8

Spectra of spacecraft orbit evolution due to viscoelastic effect. a Inclination with several distinguished frequencies marked in red. b Longitude of ascending node. Dashed lines indicate indirect effects. Units are days unless annotated otherwise

5 Error budget in orbit adjustment

The errors of neglecting the viscoelastic effect in the least-squares orbit adjustment can be roughly assessed as follows. The linearized dynamic equation is

$$\begin{aligned} \begin{aligned} \textbf{X}(t)&= \bar{ \textbf{X} }(t) + \textbf{x}(t), \\ \textbf{x}(t)&= \Phi (t,t_0)\, \textbf{x}(t_0), \quad \dot{ \Phi }(t,t_0) = \frac{ \partial \dot{ \textbf{X} } }{ \partial \textbf{X} } \Phi (t,t_0), \end{aligned} \end{aligned}$$

(28)

where $\bar{ \textbf{X} }$ is a reference state from which the true state deviates by $\textbf{x}$ and where the state transition matrix $\Phi (t,t_0)$ propagates the initial deviation to the epoch t. The linearized observation equation is

$$\begin{aligned} \textbf{y} = \textbf{H}\, \textbf{x}(t_0) + \varvec{\varepsilon }, \end{aligned}$$

(29)

with $\varvec{ \varepsilon }$ representing the noise. For the sake of simplicity, we assume that the orbit positions and thus their errors are directly measurable, so that the observation matrix is $\textbf{H} = \textbf{I}_{3} \Phi$. For a given weight matrix $\textbf{W}$, the estimate and its covariance matrix are

$$\begin{aligned} \begin{aligned} \hat{ \textbf{x} }(t_0)&= ( \textbf{H}^\textrm{T} \textbf{W}\, \textbf{H} )^{-1} \textbf{H}^\textrm{T} \textbf{W y}, \\ \textbf{P}&= ( \textbf{H}^\textrm{T} \mathbf {W\, H} )^{-1}. \end{aligned} \end{aligned}$$

(30)

Unless otherwise noted, hereafter $\textbf{x}$ refers to $\textbf{x}(t_0)$.

5.1 Consider analysis

We treated the modeled viscoelastic effects as consider parameters in the adjustment. The effects were included in the orbit and uncertainty propagation, but not estimated with the spacecraft orbit. The observation equation is augmented as,

$$\begin{aligned} \textbf{y} = \left( \, \textbf{H}\ \ \textbf{H}_\textrm{c} \right) \left( \begin{array}{l} \textbf{x} \\ \textbf{c} \end{array} \right) + \varvec{\varepsilon }, \end{aligned}$$

(31)

where $\textbf{c} = \textbf{k}_{2\,(j)}$ and $\textbf{H}_\textrm{c} = \Phi _{(j)}$, which evaluates the orbit change due to $\textbf{k}_{2\,(j)}$, as in Eq. (10). Assuming there is no correlation between $\textbf{x}$ and $\textbf{c}$, the estimated initial orbit and the associated uncertainties become (Montenbruck and Gill 2000),

$$\begin{aligned} \begin{aligned} \hat{ \textbf{x} }_{(j)}&= \hat{ \textbf{x} } + \textbf{S}\, \textbf{k}_{2\,(j)}, \\ \textbf{P}_{(j)}&= \textbf{P} + \textbf{S}\, \textbf{P}_\textrm{c} \textbf{S}^\textrm{T}, \quad \textbf{S}=-\textbf{P}\, \textbf{H}^\textrm{T} \textbf{W} \textbf{H}_\textrm{c}, \end{aligned} \end{aligned}$$

(32)

where $\textbf{S}$ is known as the sensitivity matrix describing the change of the estimate and its uncertainties with $\textbf{c}$. $\textbf{P}_\textrm{c}$ is the a priori covariance matrix of $\textbf{c}$, whose specification is discussed below.

5.2 Noiseless case

The impact of the considered parameters is largely independent of the measurement noise level, i.e., it does not decrease with noise reduction or increase of observation volume (Montenbruck and Gill 2000; Tapley et al. 2004).

We neglected noise in the present solution and therefore only assessed the impact of dynamical uncertainty (see, e.g., Dirkx et al., 2016). In this extreme case, the consider analysis has a special appeal. Without viscoelastic effects, the estimated $\hat{ \textbf{x} }$ converges to the true $\textbf{x}$ within machine precision or below a certain specified threshold, i.e., effectively free of errors. When considered in the orbit and error propagation but not solved for, the viscoelastic effects prevent the error-free solution. Rather, the estimate $\hat{ \textbf{x} }_{(j)}$ deviates from $\hat{ \textbf{x} }$ (and thus $\textbf{x}$ itself) by $\textbf{S}\, \textbf{k}_{2\,(j)}$ (Eq. 32). It follows that the change of $\textbf{P}_{(j)}$ is determinate in the noiseless case. Letting $\textbf{P}_\textrm{c} = \textbf{k}_{2\,(j)}\,\textbf{k}_{2\,(j)}^\textrm{T}$, we find

$$\begin{aligned} \textbf{S}\textbf{P}_\textrm{c} \textbf{S}^\textrm{T} = \left( \hat{ \textbf{x} }_{(j)} - \textbf{x} \right) \left( \hat{ \textbf{x} }_{(j)} - \textbf{x}\right) ^\textrm{T}, \end{aligned}$$

(33)

which indicates that the consider covariance matrix provides an exact measure of the estimation error in $\hat{ \textbf{x} }_{(j)}$. Equation (33) serves as a validating test for the consider covariance analysis (a confirmation provided later in Fig. 9).

5.3 Errors in estimated initial orbits

The estimated state vector consists of the following parameters,

$$\begin{aligned} \hat{ \textbf{X} }^\textrm{T} = \left( \hat{ \textbf{r} }_0^\textrm{T}\ \ \hat{ \textbf{v} }_0^\textrm{T}\ \ \hat{ C }_\textrm{SR}\ \ \hat{a}_\textrm{t} \right) , \end{aligned}$$

(34)

where $\hat{ C }_\textrm{SR}$ is the coefficient for the solar radiation pressure, and $\hat{a}_\textrm{t}$ is an empirical along-track acceleration that is commonly included in the orbit determination to absorb unmodeled perturbations. Note that $a_\textrm{t}$ is absent in the (simulated) true orbit, which however is influenced by the consider parameters, $\textbf{k}^\Im _{2\,(\sum j)}$ (Eqs. 10 and 31).

Orbit adjustment always involves iterations of least-squares solutions and corrections of the initial orbit $\hat{ \textbf{X} }=\bar{ \textbf{X} } + \hat{ \textbf{x} }$. Such refinement is not needed for the consider covariance analysis; instead, it suffices to obtain $\textbf{S}$ in just one solution, which changes little afterward.

Fig. 9

Errors of adjusted orbit position for the 3-day arc from 2014-04-05. The errors are projected into the radial (top), transverse (middle), and normal (bottom) directions. The errors derived from the consider covariance matrix are marked by dashed black curves. The true errors are solid curves

We analyzed the solutions of 3- and 9-day arcs. The lengths correspond roughly to the end choices in practice (Mazarico et al. 2018; Löcher and Kusche 2018). Table 3 shows the estimation errors of the initial position at six epochs in April 2014. The errors are evaluated as the square root of $\textrm{tr}[ (\textbf{S}\textbf{P}_\textrm{c} \textbf{S}^\textrm{T})_{3\times 3} ]$, where only the first (upper-left) $3\times 3$ block of the matrix is relevant.

As an illustration, Fig. 9 shows the estimation errors in the orbit position over three days from 2014-04-05 derived from consider covariance analysis in comparison with the true errors, i.e., $\hat{ \textbf{X} }-\textbf{X}$. The consider errors are nonnegative and coincide with the true errors when the latter is not negative; otherwise, they are of opposite signs. The condition of Eq. (33) is thus observed.

Over 3 days, the maximum error is between 1 and 2 decimeters. The error level oscillates and peaks after 4 arcs, i.e., 12 days. We infer that this is likely associated with $\Delta \Omega _{(1)}$ (Fig. 7a), and the effect was not (mis)absorbed by the adjusted orbits. As expected, the errors over 9 days are about 2 times more pronounced due to accumulation.

Table 3 Errors of adjusted initial orbit positions for six arcs in 2014 (units: m)

6 Discussion and conclusion

We have adopted the lunar tidal deformation model by Williams and Boggs (2015) to evaluate the motion of a polar orbiter, focusing on the impact of viscoelastic tidal perturbation, a mechanism of topical interest in lunar science. (Williams et al. 2014). The lunar tides are dominated by the monthly ($\sim$27 day) components associated with the orbital period around the Earth. The orbital spectra of a spacecraft are altered so that the frequencies of the long-term variations are offset order-wise from the multiples of the monthly, i.e., $m\dot{\theta }$. Tidal components close to such frequencies produce a wide range of composite periodicities, from a few to thousands of days. The impact of slower tides, of the annual and semiannual components for instance, is hardly separable from the respective monthly frequencies to the first order. Higher-order effects are too subtle to be identified in the analysis. Consequently, it is difficult (if not impossible) to identify slow-varying tidal components in the orbiter motion directly.

We evaluated the dynamical uncertainties in the least-squares orbit adjustment based on a simplified case neglecting noise. The errors at the level of decimeters over up to 9 days are present in the real-world precise orbit products. At present, they are overwhelmed by other, larger errors, e.g., due to additional dynamical uncertainties, limiting observation condition, and inevitably measurement errors, which amount to at least several meters (Mazarico et al. 2018; Löcher and Kusche 2018).

It can be anticipated that the viscous lag of the monthly tide, or the dissipation factor, and those at multiple frequencies can be realistically pursued in the near future, thus providing an independent probe into the viscoelastic structure of the Moon. Among the reasons are the prospective decimeter-level orbit accuracy, and also the distinct frequencies of orbit variation, many of which are even observable within a single arc. Moreover, several tidal components are likely separable in the orbit spectra. The monthly components, $j=1,2$, affect, respectively, the ascending node and inclination of the spacecraft orbit (Sect. 3.4); meanwhile, the semimonthly term, $j=4$, produces a cumulative effect over a semiannual cycle (Sect. 4). However, the same cannot be easily said for the lower frequencies, e.g., the annual one.

Measuring lunar tides at different frequencies would have broad implications for understanding the Moon’s interior and evolution. Whereas current data sets are not sufficient for a full evaluation, future lunar missions might allow for investigating the response of the lunar interior to different excitation frequencies. Those measurements would provide data sets complementary to seismic measurements at the surface. Tides and their phase lags would be in particular indicative of low-viscosity zones and the rheologic state and size of the lunar core.