Impact of Atmospheric Delay on Equivalence Principle Tests Using Lunar Laser Ranging

1. Introduction

Lunar Laser Ranging (LLR) is a scientific mission that calculates the distance between Earth-based telescope stations and lunar retroreflectors by measuring the round-trip flight time of laser pulses. Since the deployment of the first retroreflector (A11) by astronauts during the Apollo 11 mission in 1969, four additional retroreflectors were placed on the lunar surface for international scientific purposes during the last century. These include two retroreflectors deployed during the Apollo 14 and Apollo 15 missions (A14, A15), and two retroreflectors carried by Soviet Lunokhod rovers (L1, L2). In March 2025, the sixth lunar retroreflector (NGLR-1, Next Generation Lunar Retroreflector-1) was deployed on the Moon’s surface [1]. Laser pulses traversing the Earth–Moon distance are subject to extreme attenuation: the received intensity falls with the inverse fourth power of the separation [2]. This geometric loss, compounded by atmospheric turbulence and diffraction at the lunar retro-reflectors [3], reduces the mean return to fewer than one photon per pulse [4,5]; additional losses arise in the transmit and receive optics. To compensate, observatories integrate returns over several minutes (typically 5–20) and compress them into a single “normal point” via specialized algorithms. During five decades of LLR, this procedure has yielded more than 30 000 normal points and has progressively refined Earth–Moon distance measurements from decimetre to millimetre precision [3].

LLR provides important support for multiple research fields, such as lunar ephemeris calculation [6,7], determination of Earth orientation parameters [8,9], and studies of the Moon’s internal structure [10]. With the gradual improvement in ranging accuracy, LLR has now become one of the most powerful tools in the solar system for testing general relativity, such as Lorentz symmetry [11,12,13,14,15], temporal variation of the gravitational constant [16], the inverse-square law [17], the equivalence principle [16,18,19]. This study is devoted to research on testing the equivalence principle. If the Weak Equivalence Principle (WEP) is violated, it would cause the Earth and Moon to experience different accelerations in the Sun’s gravitational field, resulting in periodic oscillations in the Earth-Moon distance. Since the Earth and Moon differ in composition and possess significant gravitational self-energy, LLR can test not only the WEP but also the Strong Equivalence Principle (SEP)——a capability beyond the reach of ground-based or satellite experiments. Currently, LLR achieves a precision of ${10}^{-14}$ in testing the WEP within the Sun’s gravitational field, and a precision of ${10}^{-4}$ in testing the SEP [16,18]. LLR precision is constrained by several factors [20,21,22,23], and Earth’s atmosphere effects is one of the most important limitations for testing the equivalence principle through LLR [24].

While atmospheric effects can be assessed at the millimeter level with the aid of additional measurements such as GNSS, the vast majority of researchers lack access to these supplementary data streams. Typically, only temperature, pressure, and humidity data are available from the LLR normal point data. Using these three measurements in conjunction with atmospheric models, the residual error in atmospheric delay correction remains on the order of 5-10 %, corresponding to decimeter-level residuals. This level of residual error significantly exceeds the currently millimeter-to-centimeter measurement precision of LLR. Therefore, to meet the demand for high-precision LLR data from a broad scientific researchers, alternative strategies to mitigate the impact of atmospheric delay are required. Recently, a novel differential lunar laser ranging (DLLR) is claimed to reduce atmospheric delay errors and improve the accuracy of testing equivalence principle parameters in LLR [25]. Its practical feasibility still requires further investigation. Here, by analyzing the inherent characteristics of atmospheric delay, we have proposed suggestions for reducing the impact of the atmosphere delay on equivalence principle tests through measurement schemes or data processing methods.

By investigating how atmospheric delay, and especially its dependence on observation elevation angle, limits equivalence principle tests with Earth–Moon laser-ranging data, this paper can offer a theoretical reference for future higher-precision tests of the equivalence principle. These findings are timely given the rapid expansion of lunar exploration: the Chang’e program, Blue Ghost’s delivery of the NGLR-1 retro-reflector [1], the maturation of hollow corner-cube retro-reflectors for next-generation LLR [21,26], and lunar time-frequency comparison experiments [27,28,29]. Insights into atmospheric delay presented herein will hold significant reference value for these missions. The paper is organized as follows: Section 2 introduces the atmospheric delay model. Section 3 describes the LLR methodology for testing the equivalence principle. Section 4 analyzes the influence of different elevation angle distributions on equivalence principle tests using several datasets. Finally, Section 5 summarizes our findings and conclusions.

2. Atmospheric Delay

The observation accuracy of laser ranging technology is primarily affected by propagation model errors when signals pass through the troposphere and stratosphere, particularly the neutral atmospheric refraction delay. The lasers used for ranging mostly employ green light (532 nm) or near-infrared light (1064 nm). For these wavelength bands, the atmospheric delay mainly considers the neutral atmosphere component, while the ionospheric delay can be neglected for millimeter-level ranging accuracy.

In LLR, atmospheric delay effects manifest in two aspects: firstly, atmospheric refraction bends the laser path, resulting in increased propagation distance; secondly, the laser’s propagation speed through the atmosphere is slower than the speed of light in vacuum, causing temporal delay. As shown in Figure 1.

The atmospheric delay value is obtained by multiplying the tropospheric zenith delay (i.e., the atmospheric delay when the station’s elevation angle is 90 degrees) by a dimensionless mapping function (which projects the zenith delay to other directions and is normalized to 1 at an elevation angle of 90 degrees). Currently, there are two main models commonly used for applying atmospheric delay corrections to satellite laser ranging: the Marini-Murray model [30] and the Mendes-Pavlis model [31]. The Mendes-Pavlis model demonstrates higher accuracy, with errors at the sub-millimeter level in the zenith direction [31], gradually increasing to the millimeter level as the elevation angle decreases [32]. (However, if only the pressure, humidity, and temperature parameters of standard points are utilized without additional measurements, the accuracy is only approximately 5% to 10% [33,34].) Therefore, it is used here as the atmospheric delay model to correct the measurement data. The Mendes-Pavlis model includes a zenith delay and a mapping function, adopting the ray-tracing method to fit the constructed zenith delay and mapping function for coefficient determination. The ray-tracing technique is currently the most accurate method available. This approach computes the total delay by integrating contributions from individually modeled atmospheric layers.

2.1. Zenith Propagation Delay

The fundamental formula for neutral atmospheric refraction delay of electromagnetic waves in the zenith direction is:

$$\begin{array}{c}\hfill d_{\mathrm{atm}}^z={10}^{-6}\times {\int }_{r_s}^{r_a}Ndz={\int }_{r_s}^{r_a}(n-1)dz\end{array}$$

(1)

where $r_s$ is the geocentric radius of the laser station, $r_a$ is the geocentric radius at the top of the neutral atmospheric layer, N is called the group refractivity, $dz$ is the path integral element along the zenith direction, having length units, n is the atmospheric group refractive index along the electromagnetic wave propagation path.

The calculation formula for group refractivity is given by [35]:

$$\begin{array}{c}\hfill N=N_h+N_{nh}=\left({\displaystyle \frac{{\rho }_a}{{\rho }_{axs}}}\right)N_{gaxs}+\left({\displaystyle \frac{{\rho }_w}{{\rho }_{ws}}}\right)N_{gws}\end{array}$$

(2)

Here, the hydrostatic component of group refractivity $N_h$ is given by $N_h={\displaystyle \frac{{\rho }_a}{{\rho }_{axs}}}{N}_{gaxs}$, where ${\rho }_a$ is the dry air density (kg/m³), ${\rho }_{axs}$ is the standard dry air density at 15°C and 101325 Pa, and $N_{gaxs}$ is the dry air group refractivity. The non-hydrostatic component of group refractivity $N_{nh}$ is expressed as $N_{nh}={\displaystyle \frac{{\rho }_w}{{\rho }_{ws}}}{N}_{gws}$, with ${\rho }_w$ being the water vapor density (kg/m³), ${\rho }_{ws}$ the pure water vapor density at 20°C and 1333 Pa, and $N_{gws}$ the water vapor group refractivity. Subsequently, all six physical quantities on the right-hand side of equation (2) are represented using surface parameters. $N_h$ and $N_{nh}$ are then integrated separately and summed to obtain the zenith atmospheric delay. The integration of $N_h$ yields the zenith hydrostatic delay $d_h^z$, while the integration of $N_{nh}$ yields the zenith non-hydrostatic delay $d_{nh}^z$. The results after integration are as follows:

$$\begin{array}{c}\hfill d_h^z=0.00002416579{\displaystyle \frac{f_h\left(\lambda \right)}{f(\phi ,H)}}{P}_s\end{array}$$

(3)

$$\begin{array}{c}\hfill d_{nh}^z={10}^{-6}\left(5.316f_{nh}\left(\lambda \right)-3.759f_h\left(\lambda \right)\right){\displaystyle \frac{e_s}{f(\phi ,H)}}\end{array}$$

(4)

with

$$\begin{array}{c}\hfill f_h\left(\lambda \right)={10}^{-2}\left[k_1^{*}{\displaystyle \frac{(k_0+{\sigma }^2)}{{(k_0-{\sigma }^2)}^2}}+k_3^{*}{\displaystyle \frac{(k_2+{\sigma }^2)}{{(k_2-{\sigma }^2)}^2}}\right]C_{{\mathrm{CO}}_2}\end{array}$$

(5)

$$\begin{array}{c}\hfill f(\phi ,H)=1-0.00266cos2\phi -0.00028H\end{array}$$

(6)

$$\begin{array}{c}\hfill f_{nh}\left(\lambda \right)=0.003101\left({\omega }_0+3{\omega }_1{\sigma }^2+5{\omega }_2{\sigma }^4+7{\omega }_3{\sigma }^6\right)\end{array}$$

(7)

Here, $\lambda$ is the vacuum wavelength in $\mathrm{\mu }\mathrm{m}$, $\phi$ is the station latitude, H is the station height in units of $\mathrm{k}\mathrm{m}$, $P_s$ is the surface barometric pressure ($\mathrm{Pa}$), $e_s$ is the surface water vapor pressure ($\mathrm{Pa}$), which can be calculated using relative humidity and temperature. The coefficients are determined as $k_1^{*}=19990.975 \mathrm{\mu }{\mathrm{m}}^{-2}$, $k_3^{*}=579.55174 \mathrm{\mu }{\mathrm{m}}^{-2}$, $k_0=238.0185 \mathrm{\mu }{\mathrm{m}}^{-2}$, $k_2=57.362 \mathrm{\mu }{\mathrm{m}}^{-2}$; ${\omega }_0=295.235$, ${\omega }_1=2.6422 \mathrm{\mu }{\mathrm{m}}^2$, ${\omega }_2=-0.032380 \mathrm{\mu }{\mathrm{m}}^4$, ${\omega }_3=0.004028 \mathrm{\mu }{\mathrm{m}}^6$. $\sigma ={\lambda }^{-1}$, $C_{{\mathrm{CO}}_2}=1+0.534\times {10}^{-6}(x_c-450)$, where $x_c$ is the CO2 content, in ppm. The influence of ${\mathrm{CO}}_2$ concentration on atmospheric delay is negligible. Normal point data typically does not include ${\mathrm{CO}}_2$ concentration measurements, and it can therefore be assumed to be 375 ppm.

The final calculation of the zenith delay is given by [31]:

$$\begin{array}{c}\hfill d_{\mathrm{atm}}^z=d_h^z+d_{nh}^z={10}^{-6}\left[4.16579{\displaystyle \frac{f_h\left(\lambda \right)}{f(\phi ,H)}}{P}_s+\left(5.316f_{nh}\left(\lambda \right)-3.759f_h\left(\lambda \right)\right){\displaystyle \frac{e_s}{f(\phi ,H)}}\right]\end{array}$$

(8)

2.2. Mapping Function

All combinations of zenith delays and mapping functions are based on the theoretical assumption of a spherically symmetric Earth atmosphere. With the mapping function, we can project the tropospheric delay from the zenith direction to any arbitrary direction, thereby obtaining the atmospheric delay along the slant path. Consequently, the accuracy of the mapping function affects the precision of the slant delay, which in turn influences the estimation of atmospheric delay.

The total delay can be expressed as the product of the zenith atmospheric delay and the mapping function [36]:

$$\begin{array}{c}\hfill d_{\mathrm{atm}}=d_h^z{m}_h\left(ϵ\right)+d_w^z{m}_w\left(ϵ\right)\end{array}$$

(9)

where $ϵ$ is the laser elevation angle. $m_h\left(ϵ\right)$ and $m_w\left(ϵ\right)$ are the hydrostatic and wet mapping functions, respectively. Since water vapor’s contribution to atmospheric refraction in the visible spectrum is negligible, under these conditions we have:

$$\begin{array}{c}\hfill d_{\mathrm{atm}}=d_{\mathrm{atm}}^{z}m\left(ϵ\right)\end{array}$$

(10)

The mapping function $m\left(ϵ\right)$ is usually based on a truncated continued fraction expansion of ${\displaystyle \frac{1}{sinϵ}}$ [37]. This empirical formula is normalized to unity at the zenith. The mapping function is given by:

$$\begin{array}{c}\hfill m\left(ϵ\right)={\displaystyle \frac{1+{\displaystyle \frac{a_1}{1+{\displaystyle \frac{a_2}{1+a_3}}}}}{sinϵ+{\displaystyle \frac{a_1}{sinϵ+{\displaystyle \frac{a_2}{sinϵ+a_3}}}}}}\end{array}$$

(11)

The coefficients $a_i$ in equation (11) can be calculated using only the station position and surface temperature. The coefficients of the mapping function have the following mathematical expressions:

$$\begin{array}{c}\hfill a_i=a_{i0}+a_{i1}{t}_s+a_{i2}cos\phi +a_{i3}h(i=1,2,3)\end{array}$$

(12)

where $a_{ij}$ are known quantities, $t_s$ is the station temperature (°C), h is the sorthometric height of the station (m), and the remaining parameters can be obtained from reference [38].

3. Lunar Laser Ranging and Equivalence Principle Tests

WEP is equivalent to the statement that “gravitational mass equals inertial mass”. Gravitational mass $m_g$ is defined via gravitational effects, characterizing an object’s capacity to generate and respond to gravity, and inertial mass $m_i$ is introduced through Newton’s second law as a measure of an object’s inertia. If the $m_g/m_i$ ratio differs between two objects—indicating a violation of WEP—their accelerations in the same gravitational field would diverge. The Earth-Moon system in the Sun’s gravitational field serves as an exemplary testbed for this principle:

$$\begin{array}{c}\hfill {\mathit{a}}_{\mathrm{E}}={\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{E}}\mathit{g}\end{array}$$

(13)

$$\begin{array}{c}\hfill {\mathit{a}}_{\mathrm{M}}={\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{M}}\mathit{g}\end{array}$$

(14)

where $\mathit{a}$ denotes the acceleration of an object in the Sun’s gravitational field, where the subscripts E and M represent the Earth and Moon respectively, and $\mathit{g}$ is the acceleration in the Sun’s gravitational field. If the WEP is violated:

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{E}}\ne {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{M}}\end{array}$$

(15)

it will result in:

$$\begin{array}{c}\hfill {\displaystyle \frac{\Delta \mathit{a}}{\mathit{a}}}={\displaystyle \frac{{\mathit{a}}_{\mathrm{E}}-{\mathit{a}}_{\mathrm{M}}}{\frac{1}{2}({\mathit{a}}_{\mathrm{E}}+{\mathit{a}}_{\mathrm{M}})}}\approx {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{E}}-{\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{M}}=\Delta {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{EM}}\end{array}$$

(16)

This implies that the Earth-Moon distance would exhibit oscillations. Nordtvedt treated the equivalence principle-related tidal term as a perturbation, revealing a polarization phenomenon of the Moon’s orbit along the solar direction accompanied by radial perturbations [39], with a signal period of 29.53 days (distinct from the Moon’s orbital period of 27 days).

$$\begin{array}{c}\hfill d{r}_{\mathrm{EM}}=\mathrm{S}\Delta {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{EM}}cosD\end{array}$$

(17)

The additional acceleration resulting from the violation of the equivalence principle leads to an additional range term $r_{\mathrm{EM}}$ between Earth and Moon. $\mathrm{S}=-2.943\times {10}^{10}\mathrm{m}$ [40]. D is the lunar phase angle. Biskupek et al. used LLR data to obtain $\Delta {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{EM}}=(-2.1\pm 2.4)\times {10}^{-14}$ [18].

The SEP extends the equivalence principle to encompass gravitational properties arising from gravitational self-energy, with the subject of discussion being “self-gravitating systems”. Testing the SEP requires selecting celestial systems where gravitational self-energy cannot be neglected, and the Earth-Moon-Sun system serves as an ideal experimental platform for verifying the SEP. Therefore, for the Earth-Moon system within our solar system, both WEP and SEP can be tested concurrently. The mass ratio ${\displaystyle \frac{m_g}{m_i}}$ for an object of mass M can be expressed using the Nordtvedt parameter $\eta$ [16]:

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{SEP}}=1+\eta {\displaystyle \frac{U}{M{c}^2}}\end{array}$$

(18)

where the gravitational self-energy ${\displaystyle U=-\frac{3G{M}^2}{5R}}$, with R being the radius. In General Relativity, $\eta =0$. If only the SEP is violated:

$$\begin{array}{c}\hfill d{r}_{\mathrm{EM}}=\mathrm{S}\eta \left[{\displaystyle \frac{U_{\mathrm{E}}}{M_{\mathrm{E}}{c}^2}}-{\displaystyle \frac{U_{\mathrm{M}}}{M_{\mathrm{M}}{c}^2}}\right]cosD\end{array}$$

(19)

The difference of the gravitational self energy between Earth and Moon is:

$$\begin{array}{c}\hfill {\displaystyle \frac{U_{\mathrm{E}}}{M_{\mathrm{E}}{c}^2}}-{\displaystyle \frac{U_{\mathrm{M}}}{M_{\mathrm{M}}{c}^2}}=-4.45\times {10}^{-10}\end{array}$$

(20)

Substituting this value into Eq. (19) yields:

$$\begin{array}{c}\hfill d{r}_{\mathrm{EM}}={\mathrm{C}}_0\eta cosD\end{array}$$

(21)

where ${\mathrm{C}}_0=13.1\mathrm{m}$. The Nordtvedt parameter can be estimated using Earth-Moon laser ranging:

$$\begin{array}{c}\hfill \eta =\left(-0.2\pm 1.1\right)\times {10}^{-4}\end{array}$$

(22)

4. Impact of Atmospheric Delay on Equivalence Principle Tests

The schematic diagram of LLR is shown in Figure 2, with its fundamental principle being:

$$\begin{array}{c}\hfill d=c{\displaystyle \frac{\tau }{2}}=\left|{\mathit{r}}_{\mathrm{EM}}-{\mathit{r}}_{\mathrm{station}}+{\mathit{r}}_{\mathrm{reflector}}\right|+c(\Delta {\tau }_{\mathrm{atm}}+\Delta {\tau }_{\mathrm{GR}})\end{array}$$

(23)

Here, d is the distance between the station and the reflector, c is the speed of light, $\tau$ is the pulse round-trip propagation time, ${\mathit{r}}_{\mathrm{EM}}$ is the vector connecting the Earth’s center and the Moon’s center, ${\mathit{r}}_{\mathrm{station}}$ is the geocentric position vector of the observatory, and ${\mathit{r}}_{\mathrm{reflector}}$ is the selenocentric position vector of the reflector array. The initial coordinate systems for the ground station and lunar reflector are International Terrestrial Reference System (ITRS) and Principal Axis System respectively, while the coordinate system for laser ranging data processing is Barycentric Celestial Reference System (BCRS). Therefore, conversion to BCRS is required for both, during which coordinate corrections due to tidal effects must be considered, typically at the level of several tens of centimeters. $\Delta {\tau }_{\mathrm{atm}}$ describes the propagation time correction caused by atmospheric effects. This term varies significantly, generally ranging from 2 to 10 m, but achieves sub-centimeter accuracy when the elevation angle exceeds 20° [41]. $\Delta {\tau }_{\mathrm{GR}}$ represents the general relativistic delay, also known as the Shapiro effect, with corrections typically around 7 to 8 m.

Normal points are the final observational data provided by LLR. The content of these observations mainly includes the laser emission time from the ground station, using Coordinated Universal Time (UTC), the round-trip flight time of the laser and its accuracy, meteorological data from the ground station (pressure, temperature, relative humidity, etc.), laser azimuth and elevation angles, laser wavelength, and station and retroreflector identifiers. After obtaining these data, each station will upload the data to the International Laser Ranging Service (ILRS) website for use.

After more than 50 years of LLR experiments, various stations have collectively contributed over 30,000 normal point data. During 1975-1985, the precision of LLR data obtained by LLR stations improved to 15 cm. Subsequently, through enhancements in laser and detector performance, the ranging precision of LLR data was elevated to the centimeter level. From the 1990s to the early 21st century, the precision of LLR data showed no significant improvement, reaching a root mean square value of 2 cm by 2015. LLR has significantly enhanced the precision of Earth-Moon distance measurements. Currently, owing to advancements in high-precision measurement equipment (e.g., the development of picosecond lasers), improved accuracy of time systems, and refined modeling approaches, the generated normal points have achieved millimeter-level precision [3], while the random errors in ranging have improved from the initial decimeter level to the present millimeter level.

The temporal distribution of these normal points is non-uniform, with significantly fewer data points near full moon and new moon (as shown in Figure 3). This occurs because: during new moon, the Moon’s proximity to the Sun makes operations difficult; during full moon, solar illumination noise is too strong [42]. These new- and full-Moon normal points correspond to $cosD\approx 1$ and are the most sensitive for equivalence principle tests. Although the number of normal points is lower at these phases, practical analyses use the entire data set—more points yield tighter constraints—and do not discard phases with smaller counts. Comparatively, infrared light demonstrates more uniform distribution than green light, representing one of the current advantages of using infrared light for LLR [43].

4.1. Current Impact Level of Atmosphere Delay on Equivalence Principle Test

Among the five lunar retroreflectors deployed in the last century, A15 has the largest reflective area, with the number of normal points obtained from A15 observations accounting for nearly two-thirds of the total from all retroreflectors. We therefore adopt A15’s data. Subsequently, we fit the cosine signal of the equivalence principle to several sets of atmospheric delay data to observe their amplitude intensities. The data processing flowchart is shown in Figure 4, with each set of atmospheric delay data undergoing the complete procedure.

The Mendes-Pavlis atmospheric delay model itself has inherent errors. According to Ref. [44], when the elevation angle is above 45°, the residual error compared to the ray-tracing method is at the millimeter level. As this model error eludes rigorous parameterisation, we approximate it as white noise of millimetre amplitude in the subsequent analysis, thereby obtaining a conservative estimate of its influence on equivalence principle tests.

From recent A15 observations, we extracted atmospheric-delay data provided by Grasse—the station contributing the largest number of normal points in recent years. A cosine

$$\begin{array}{c}\hfill y\left(x\right)=Acos(0.2128x+0.6064)+C\mathrm{m}\end{array}$$

(24)

was fitted to these delays with the fixed synodic period T=29.53 days, and initial phase locked to the lunar phase of the first record which can be computed using publicly available source code; only amplitude A and a constant offset C were allowed to vary (C is irrelevant to the test of the equivalence principle). The temporal distribution of normal points covers the period between October 2022 and January 2025. The best-fit values are $A=0.3345\pm 0.1710\mathrm{m}$. Interpreting this modulation as an equivalence principle violation yields $\Delta {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{EM}}=(-1.14\pm 0.58)\times {10}^{-11}$ and $\eta <\left(2.55\pm 1.31\right)\times {10}^{-2}$. Although atmospheric delay itself does not directly test the equivalence principle, it contributes an uncertainty at the ${10}^{-12}$ level in $\Delta {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{EM}}$—about two orders of magnitude above current experimental precision. The fit is illustrated in Figure 5.

4.2. Impact of Temporal and Elevation Angle Inhomogeneities in LLR Data on Equivalence Principle Tests

To systematically evaluate how temporal and elevation-angle inhomogeneities in LLR data influence equivalence principle tests, we simulated atmospheric-delay signatures at the Grasse station and generated the following four cases for discussion.

Case 1 : Compared with the parameters used for the analysis in Figure 5, the elevation angle of the normal point was adjusted to a random value between 30° and 60°, while all other information remained unchanged, to simulate a scenario with a uniform distribution of elevation angles. Because atmospheric delay is strongly elevation-dependent, every Monte-Carlo realization produces a different amplitude. After 1000 simulations of the absolute amplitude, we obtain the amplitude as $A=0.0170\pm 0.0010\mathrm{m}$, corresponding to an equivalence principle bound $\Delta {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{EM}}=(-5.78\pm 0.34)\times {10}^{-13}$. The results of a representative random sampling and the corresponding fit are shown in Figure 6.

Case 2 : All parameters are identical to those in Case 1 except that the elevation angle is fixed at 40°, the corresponding results are presented in Figure 7. After 1000 fits, we obtain the amplitude as $A=0.001373\pm 0.000003\mathrm{m}$, corresponding to an equivalence principle bound $\Delta {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{EM}}=(-4.67\pm 0.01)\times {10}^{-14}$. When the elevation angle is set to fluctuate randomly within a small range between 40° and 45°, the amplitude is approximately twice that of the fixed elevation angle scenario. This case demonstrates a significant improvement over Case 1.

Case 3 : We adopted the following simulation protocol for normal-point generation: (1) Within each 6-day sampling window, one observation group was randomly selected; (2) Each group comprised 5 normal points with strictly controlled 15-minute inter-point intervals. The elevation angles were constrained by two conditions: (a) global uniform distribution between 30° and 60° across all points, and (b) intra-group angular coherence maintaining $\le 5^{\circ }$ variation within each 5-point cluster. The first time point in our simulation corresponds to the initial timestamp of the empirical normal point data, thereby allowing direct matching with the initial phase precomputed in Eq. (24). The meteorological data for each normal point were obtained from [45], with the air temperature assigned as a random value between the lowest and highest monthly average temperatures for the corresponding month, while the atmospheric pressure and relative humidity were set to the monthly average values. Using Grasse atmospheric data, we simulated two years of observations and performed 1000 Monte-Carlo fits. The resulting absolute amplitude is $A=0.0351\pm 0.0031\mathrm{m}$, yielding an equivalence principle bound $\Delta {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{EM}}=(-1.19\pm 0.11)\times {10}^{-12}$. A representative sample and its best-fit cosine are shown in Figure 8.

Case 4 : This is identical to Case 3, but with the elevation angle fixed at 40°, as shown in Figure 9. After 1000 iterations, we obtain the amplitude as $A=0.00013\pm 0.00001\mathrm{m}$, corresponding to an equivalence principle bound $\Delta {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{EM}}=(-4.42\pm 0.34)\times {10}^{-15}$. Figure 9 represents an idealized simulation, the practical implementation of which is highly unlikely in experimental settings. When the elevation angle is set to fluctuate randomly within a small range between 40° and 45°, the amplitude as $A=0.0071\pm 0.0006\mathrm{m}$, corresponding to an equivalence principle bound $\Delta {\left({\displaystyle \frac{m_g}{m_i}}\right)}_{\mathrm{EM}}=(-2.41\pm 0.20)\times {10}^{-13}$.

Comparing the four cases reveals that: a uniform distribution of elevation angles can effectively reduce the fitted amplitude; fixing the elevation angle can further decrease the amplitude and thereby improves the accuracy of the equivalence principle test; temporal non-uniformity in the distribution of normal points has little direct effect on the amplitude, yet it causes signal loss near the critical phases (90° and 270°), degrading the test precision. Therfore, once a large pool of data is available in the future, selecting a subset with fixed or nearly constant elevation angles will yield tighter constraints on the equivalence principle.

For the red fitting curves above, the achievable fitting accuracy is inherently limited. For Figure 5, the real atmospheric delay is intrinsically modulated at the lunar-phase frequency. This modulation arises because the mapping function in the delay model depends on telescope elevation, which in turn is set by the instantaneous Earth-Moon geometry. Consequently, when the delay is projected onto the equivalence principle test frequency, the resulting upper limit is large. The low fitting accuracy obtained for the red-line fit is mainly due to the non-uniform distribution of the data points, which does not fully sample the cyclic behaviour. For the other cases, the low fitting accuracy stems from two factors: (i) the data non-uniformity and (ii) the absence of lunar phase modulation when the elevation is held fixed or randomized. In the situations the atmospheric delay behaves like white noise: in the frequency domain its amplitude is nearly constant across all Fourier components. Since only the equivalence principle violation frequency is interested in, a single fixed-frequency cosine is fitted; hence the fitting accuracy is inevitably low.

5. Conclusion

Atmospheric delay demonstrates high sensitivity to elevation angles. The temporal non-uniformity of ranging data not only leads to the absence of critical phase data but also results in an inhomogeneous elevation angle distribution, thereby inducing significant fluctuations in atmospheric delay that ultimately affect the testing accuracy of fundamental physical parameters. This paper focuses primarily on intrinsic influence of the atmospheric delay effect itself (rather than analysing the post-fit residuals that remain after subtracting an atmospheric model from the observations) on equivalence principle tests. Atmospheric delay introduces approximately 2-10 m correction in one-way propagation, with an error lower bound of about 5 mm, making it one of the primary factors limiting LLR precision. Fits to data from the Grasse station reveal a $0.3345\pm 0.1710\mathrm{m}$ cosine signature in the delay, which translates into an equivalence principle-violation bound of $(-1.14\pm 0.58)\times {10}^{-11}$—two orders of magnitude above current experimental precision and confirming that the uncertainty of atmospheric effect is a leading systematic limit in equivalence principle tests.

To further investigate how atmospheric delay and its inhomogeneous distribution with respect to elevation angle affect the precision of equivalence principle tests in LLR, we conducted a simulation focused on the atmospheric delay at the Grasse station. The results indicate that the highest improvement in test accuracy is achieved when the atmospheric delay is modeled at a fixed elevation angle. Allowing small angular variations—for instance, five degrees—yields the second-best performance. Both scenarios outperform those with non-uniform elevation angle distributions.

The conclusion is straightforward: the total atmospheric delay is the zenith atmospheric delay multiplied by a mapping function that scales with telescope elevation. Because elevation is dictated by the instantaneous Earth-Moon geometry, the mapping function, and hence the entire atmospheric correction, carries the lunar-phase periodicity. Consequently, the atmospheric delay term is intrinsically modulated at the Nordtvedt frequency and projects almost entirely onto the equivalence principle signal. When the atmospheric delay is projected onto the equivalence principle violation frequency, its amplitude is large. In contrast, simulations in which elevation is held fixed or randomized break the phase-lock to the lunar cycle. Projecting these realizations onto the equivalence principle violation frequency then produces a markedly smaller amplitude.

Although maintaining a fixed elevation angle is currently experimentally impractical, it will be feasible to acquire a large volume of normal point data in the future—especially after the NGLR-1 station becomes operational—and to selectively process datasets characterized by minimal fluctuations in elevation angle. This approach will effectively reduce errors caused by atmospheric delay. Increasing the proportion of infrared laser ranging and improving data coverage across lunar phases, particularly during full and new moons, would also enhance equivalence principle test precision.

Moreover, it is worth noting that our analysis of atmospheric delay model errors is based on the current literature, which indicates that the model residual error compared to the ray-tracing method is about at the millimeter level. Since this error is not straightforward to characterize, we treat it as white noise with millimetre-level amplitude to obtain a conservative estimate of its impact on equivalence principle tests. In fact, rigorously quantifying the relationship between elevation angle and equivalence principle test accuracy is significantly important and should be regarded as a priority for future related investigations.

In addition, any complete effort to constrain equivalence principle violation parameters from the post-fit residuals of lunar laser ranging must fully account for a broad set of error sources that originate primarily from the retro-reflectors, instrumentation, and environment. Lunar-libration-induced tilts of the corner-cube arrays remain the dominant random error because they temporally broaden the returned laser pulse. Instrumental noise is introduced by ground-based hardware such as avalanche-photodiode detectors, the laser, and the event-timing system. Signal propagation is perturbed by atmospheric delay, while geophysical effects—including ocean loading and geocentre motion—also imprint signatures on the ranges. All of these contributions combine to form the LLR error budget. Follow-up studies should therefore incorporate these additional error sources and exploit the residual data to set more stringent limits on equivalence principle violation parameters.