© ESO 2019

1. Introduction

Earth precession-nutation models describe the long-term and long-period changes of the celestial intermediate pole (CIP) direction in the geocentric celestial reference system (GCRS). Before 2000, the IAU 1976/1980 precession-nutation model (Lieske et al. 1977; Wahr 1981; Seidelmann 1982) was adopted, which is based on a celestial reference system based on the bright star catalog FK5. With the help of very long baseline interferometry (VLBI), these models were developed and thus finally replaced previous models with the IAU 2006 precession (Capitaine et al. 2003) and IAU 2000A nutation (Herring et al. 2002; Mathews et al. 2002) models.

The CIP axis is used as a rotation axis to orient the international terrestrial reference system with respect to a kinematically nonrotating celestial coordinate system, the latter being given by a group of distant radio source positions (Fey et al. 2015). In this work, we focus on one part of the Earth orientation parameters (EOP), namely the celestial pole offsets (CPO), which represent the differences between observations and theoretical predictions of the CIP locations (i.e., dX, dY).

By far, VLBI has been playing the most crucial role in this domain, given its microarcsecond- level accuracy. However, VLBI observations have been available for a short time so far (about 38 years, 1979–2017), which can be a disadvantage in revealing effects with long periods such as the deficiencies in the Earth precession model, which has a period of about 26 000 years.

The satellite space-geodetic techniques, such as the Global Positioning System (GPS), lunar laser ranging (LLR), or satellite laser ranging (SLR), can also theoretically define an inertial reference frame in a dynamical sense. LLR, which has been operated since 1969, is the only space-geodetic technique now capable to offer a stable dynamical reference system, with sufficient determination accuracy of the lunar orbit (Zerhouni & Capitaine 2009, denoted as ZC09 hereafter). ZC09 provided a basic method of determining CPO from LLR observations using data in the interval of 1969–2008. They have also estimated the corrections to nutation terms, compared the results with those of VLBI observations, and presented joint analyses of VLBI and LLR. The typical standard deviation of the fitted coefficients is about 0.1–0.2 mas for the nutation terms, and 3 mas cy⁻¹ for the secular terms. Hofmann et al. (2018, denoted as H18 hereafter) have estimated the reflector coordinates, station coordinates, and velocities and EOP at the same time with LLR data during 1969–2016. The typical standard deviation for corrections to nutation terms in H18 is about 0.1 mas. It was concluded in both works that the LLR technique lacks accuracy as well as sufficient observations to determine the Earth-related parameters with the same level of accuracy as VLBI.

The history of optical observations, traced back to the nineteenth century, is much longer than that of either VLBI or LLR observations. Vondrák & Štefka (2010) have constructed a series of Earth orientation catalogs (EOC) based on a combination of HIPPARCOS and Tycho catalogs and long lasting ground-based astrometric observations. With robust determinations of the proper motions of a binary star system from 4.5 million individual observations, the latest version, the EOC-4 catalog (Vondrák 2011, priv. comm.), contains 4418 celestial objects (including stars, double star components, and photocenters). These catalogs were used to derive the CPO during the time interval 1899.7–1992.0, in the HIPPARCOS celestial reference frame (HCRF), which can be regarded as the optical counterpart of the ICRS. The current EOP solution, denoted OA10 (Vondrák 2011, priv. comm.), was composed of a series of polar motion, UT1, but the CPO were only presented as quadratic functions of time. Thus, we used the CPO series of the penultimate version (OA00) in this work.

We here aim to improve the determination of precession-nutation corrections based on analyses of VLBI (Sect. 3), LLR (Sect. 4), and optical (Sect. 5) observations, and we also present joint analyses of VLBI and LLR data (Sect. 6). We follow the basic method of ZC09 for the analyses of LLR and the joint analyses of LLR and VLBI, and make a comparison with the results of H18. Further discussions and the potential of LLR are presented in Sect. 7. The time distribution of CPO series derived from LLR, VLBI (opa2018a), and optical observations are plotted in Fig. 1 to show the differences in the number of observations and in the time span.

Fig. 1. Time distributions of LLR, VLBI (opa2018a), and optical CPO series. The time span of every CPO series is LLR 1970.6–2017.7, VLBI (opa2018a) 1979.9–2018.5, and optical (OA00) 1899.7–1991.9.

2. Models used in the CPO analyses

We considered several models that are commonly used in CPO analyses (e.g., Capitaine et al. 2009) to estimate corrections to the long-period part of the precession-nutation models. All models are functions of t (centuries from the basic epoch). The basic epoch was set to be J2000.0 for the VLBI and LLR series but to be J1956.0 for optical series, which is almost the central epoch of the time coverage of the optical observations.

Three empirical models were fit to the CPO series in our analyses:

1. A parabola, that is, a quadratic function of t.

2. A linear term and 18.6-year nutation term:

   $$\begin{array}{c}\hfill \mathrm{d}X,\mathrm{d}Y=A_0+A_{1}t+A_{\mathrm{s}}sin{\mathrm{\Omega }}_1+A_{\mathrm{c}}cos{\mathrm{\Omega }}_1,\end{array}$$(1)

3. A linear term, and 18.6-year and 9.3-year nutation terms:

   $$\begin{array}{cc}\hfill \mathrm{d}X,\mathrm{d}Y& =A_0+A_{1}t+A_{\mathrm{s}1}sin{\mathrm{\Omega }}_1+A_{\mathrm{c}1}cos{\mathrm{\Omega }}_1\hfill \\ \hfill & +A_{\mathrm{s}2}sin{\mathrm{\Omega }}_2+A_{\mathrm{c}2}cos{\mathrm{\Omega }}_2,\hfill \end{array}$$(2)

   where Ω₁ = 2π/0.186 and Ω₂ = 2π/0.093, representing the arguments of the two principle nutation terms. The 18.6-year term is related to the motion of the lunar node.

All models were used in VLBI analyses for comparison with the other two techniques. Models 2 and 3 are used in LLR analyses to estimate corrections to the IAU 2006/2000 precession-nutation models. Models 1 and 2 are used in optical analyses to estimate effects with relatively longer periods.

3. CPO analyses with VLBI observations

We used the CPO series from 1979 to 2018 derived by the Calc/Solve software at the Paris Observatory analysis center (OPA)¹ (Fig. 2) and also the series derived by the Goddard Space Flight Center (GSF)² (not plotted because it is very similar to the OPA series) in the framework of the international VLBI service for geodesy and astronomy (IVS, Nothnagel et al. 2017). The quasi-periodic free core nutation (FCN), with a period of 430.23 days, stands out in the structures of the series, and needs to be removed before the residuals are analyzed. We removed this effect according to the model recommended by the IERS conventions 2010 (Petit & Luzum 2010; updated yearly³), and the OPA series without FCN are shown in Fig. 3.

Fig. 2. VLBI celestial pole offsets from the OPA analysis center with respect to the IAU 2006/2000 models.

Fig. 3. dX and dY residuals of the OPA VLBI time series with respect to the IAU 2006/2000 models. FCN is removed in the updated IERS model.

We fit all of the three models to the residuals and show the results in Table 1; the two data sources are presented for comparison. The uncertainties of the fit parameters are smaller than those presented in Capitaine et al. (2009) because the data have accumulated over ten more years and also because they have a higher quality. A parabola clearly is not an effective model to fit the curvature because the secular and quadratic term are strongly correlated (−0.9). For the second model, the weighted root mean square (WRMS) is reduced by about 15%. The correlation coefficient between the secular term and the sine term of the 18.6-year nutation (0.4) of model 3 stands out among others. The short time-span of the VLBI observations (38 years) leads to an evident correlation between the 18.6-year nutation term and long-period terms (t¹ and t²). All the correlation coefficients can be found in Appendix A.

By comparing the results of models 2 and 3, we can see that there is no significant difference between the fitted coefficients of the 18.6-year term with and without the 9.3-year term. The correlation coefficients between the two nutation terms are quite small (∼0.2), indicating that they are well separated. The results of CPO series obtained from the two data centers are highly consistent in terms of nutation terms and uncertainties of all fit coefficients, but differences can be seen especially in the secular terms of dY. This is probably related to the strong correlation between constant and secular terms (−0.8 in both solutions) and to the significant uncertainties compared to those of other estimated terms.

The fitted coefficients of the secular terms show an underestimation in the IAU model of the precession rate in X about 0.3 mas cy⁻¹. Capitaine et al. (2009) have pointed out that the 18.6-year nutation term is the most sensitive to the error in the precession-nutation model, with VLBI data up to 2008. With the accumulated high-precision data, we here obtain more consistent results of the corrections of the 18.6-year nutation term. All fit coefficients of the 18.6-year nutation term reveal that the amplitude is underestimated by about 35 μas. However, the formal errors of the CPO derived from VLBI data are probably underestimated by about a factor 2 according to Herring et al. (2002).

4. LLR observations and CPO series

The LLR observations are presented as so-called normal points. They refer to lines of data that contain the emission time of the laser, the observed round trip time in UTC, the telescope and reflector ID, and some atmospheric parameters of each observation. These data can be used to calculate the round-trip times and then the residuals of the round-trip time [observation minus calculation (O−C)], which can be converted into residuals in one-way distance in centimeters. Finally, we obtain CPO series based on these residuals.

4.1. LLR data

We used the residuals spanning 1970–2017 (O−C of the one-way distance in centimeters) provided by Pavlov et al. (2016) ⁴. The rejection procedure is taken in advance to exclude data with relatively poor quality. The rejection criterion is such that O−C values that are higher than three times the respective formal error and the total WRMS of all normal points of the same station are excluded. The numbers of the used and rejected normal points and the corresponding time coverages of each LLR ground-based station are listed in Table 2. Data of the McDonald Laser Ranging Station (MLRS2) between 2000 and 2015 have a relatively lower accuracy and accordingly cause lower accuracies in the CPO that were obtained between 2000 and 2005 (the WRMS of MLRS2 in this time span is 12.18 cm, while the overall WRMS of this time span is 4.54 cm) in the following sections.

4.2. General methods of calculating CPO from LLR residual

ZC09 offered a method of calculating the CPO from LLR residual. We generally followed their method in converting LLR residuals into CPO (i.e., Eqs. (3)–(6)). The methods of deriving LLR residuals from the observational data (also known as “normal points”) in ZC09 (Chapront et al. 1999) and Pavlov et al. (2016) are basically the same. It is noted in both methods that not separating the two ways of light travel in the calculation of the round-trip time may cause errors of hundreds of meters in residuals of the one-way distance. Therefore, different from ZC09, we also treated the forward and backward ways of light travel separately in this part for a more solid estimation (further discussion in Sect. 4.5), namely taking the round-trip time as Δt = (D₁ + D₂)/c instead of Δt = 2D/c, in which D₁ and D₂ are the station-reflector vectors of the forward and backward ways of light travel in the barycentric celestial reference system (BCRS).

The station-reflector vector can be written as the summation of three vectors, as presented in ZC09,

$$D=SE+EM+MR,$$(3)

where −SE is the geocentric position of the ground-based station, EM is the geocentric position of the moon, and MR is the selenocentric position of the reflector. In order to calculate the partial derivative of Δt with respect to the celestial pole coordinate X and Y, we need to use the transformation from ITRS to GCRS for the station coordinates:

$$\begin{array}{c}\hfill {\mathit{SE}}_{\mathrm{GCRS}}=\mathcal{Q}(t)\mathcal{R}(t)\mathcal{W}(t){\mathit{SE}}_{\mathrm{ITRS}}=\mathcal{M}(t){\mathit{SE}}_{\mathrm{ITRS}}.\end{array}$$(4)

In the expression above, only 𝓠, the transformation matrix for precession and nutation, is related to X and Y:

$$\begin{array}{c}\hfill \mathcal{Q}=\mathcal{PN}(X,Y){\mathcal{R}}_z(s),\end{array}$$(5)

where 𝒫𝒩 is the precession nutation matrix and s is the CIO locator. Thus the partial derivatives of the calculated time delay with respect to X and Y can be written as

$$\begin{array}{cc}& \frac{\partial \mathrm{\Delta }t}{\partial X}=\frac{1}{c}[\frac{{\mathit{D}}_1}{D_1}(-\frac{\partial \mathcal{M}}{\partial X}(t_1,t_2){\mathit{ES}}_{\mathrm{ITRS}})+\frac{{\mathit{D}}_2}{D_2}(-\frac{\partial \mathcal{M}}{\partial X}(t_3,t_2){\mathit{ES}}_{\mathrm{ITRS}})]\hfill \\ \hfill & \frac{\partial \mathrm{\Delta }t}{\partial Y}=\frac{1}{c}[\frac{{\mathit{D}}_1}{D_1}(-\frac{\partial \mathcal{M}}{\partial Y}(t_1,t_2){\mathit{ES}}_{\mathrm{ITRS}})+\frac{{\mathit{D}}_2}{D_2}(-\frac{\partial \mathcal{M}}{\partial Y}(t_3,t_2){\mathit{ES}}_{\mathrm{ITRS}})],\hfill \end{array}$$(6)

in which t₁ is the emission time, t₂ is the reflection time, and t₃ is the reception time. We follow Sect. 4 in Pavlov et al. (2016) to calculate t₂ and t₃, but ignore some minor effects such as solid tide deformations of the Earth and the Moon because we do not need an accuracy of the moment of time as high as that in the calculation residuals of the one-way distance. The relativistic transformations between terrestrial and selenocentric reference systems to the BCRS (Petit & Luzum 2010) were applied to transform the vectors in Eq. (3) into BCRS. The relativistic gravitational delay (Kopeikin 1990) and the tropospheric delay (Mendes & Pavlis 2004; Mendes et al. 2002) were also taken into account. We used Jet Propulsion Laboratory (JPL) planetary ephemerides DE430 (Folkner et al. 2014) for the geocentric and barycentric positions of the Sun, Earth, Moon, and major planets, and also for the lunar libration angles and TT−TDB transformations.

After obtaining the partial derivatives in Eq. (6), we derived dΔt by dividing twice the residuals in one-way distance with the speed of light in vacuum. Then we obtained the CPO by weighted least-square fits using

$$\begin{array}{c}\hfill \mathrm{d}\mathrm{\Delta }t=\frac{\partial \mathrm{\Delta }t}{\partial X}\mathrm{d}X+\frac{\partial \mathrm{\Delta }t}{\partial Y}\mathrm{d}Y.\end{array}$$(7)

4.3. Uncertainties of the determined CPO

Uncertainties are presented as two and three times the formal error in ZC09 and H18, respectively, after fitting the CPO series to the models to account for unmodeled effects and further model deficiencies. The factor three was checked by analyses of subsets of used LLR normal points. Here we present an estimation of uncertainties in advance of obtaining CPO series from LLR normal points.

We investigated the uncertainty of the CPO derived from LLR observations, considering that many aspects would add to the uncertainty, except for the observational formal errors. According to IERS conventions (Petit & Luzum 2010), the accuracy of satellite and lunar laser ranging (SLR and LLR) is greatly affected by the residual errors in modeling the effect of signal propagation through the troposphere and stratosphere. The traditional approach in laser ranging data analyses used a model developed in the 1970s (Marini & Murray 1973), until new mapping functions (Mendes et al. 2002) and the zenith delay model (Mendes & Pavlis 2004) were used after January 1, 2007. Mendes & Pavlis (2003) assessed the “model minus ray tracing (cm)” for several wavelengths used in LLR and SLR techniques with the elevation angle set to be 10° for both models mentioned above. For the typical wavelength used in the LLR technique (532 nm), the assessed “model minus ray tracing (cm)” was 0.82 cm for the currently used models.

Moreover, uncertainties of orbit and libration determinations of the Moon will also add to that of the calculated CPO. We estimated the influence of these aspects by using different ephemerides in our calculation of the CPO and comparing the results. We applied two different ephemerides (DE430 and INPOP17a (Viswanathan et al. 2017) from IMCCE of Paris Observatory) to the calculation of the CPO series to estimate the level of differences in residuals in the one-way distance. The weighted mean of the resulting differences in the one-way distance is 0.11 cm (DE430−INPOP17a).

Uncertainty estimates of the calculation process can only be obtained in a statistical way as presented above because the process of determining residuals of LLR observations is extremely complicated and prevents us from obtaining actual uncertainties for each normal point. The sum of the two estimates above is approximately twice the formal error of LLR observations (0.61 cm).

Because of the complexity and according to the results of statistical tests above, we multiplied the formal errors of residual in the one-way distance with a factor three before we obtained the CPO with weighted least-square fits to compensate for the mixture of systematic and random errors. In this way, the weighted average uncertainty of CPO determined from LLR observations is 0.061 mas and is almost comparable with VLBI, but we stress that the quality of the CPO series also depends on the resolution (∼70 days for LLR instead of 1 day for VLBI), which appears to be another shortcoming of LLR that is due to the poor time distribution of the observation.

The fitting interval of each CPO was chosen as 70 days as in ZC09 for the whole time span of LLR observations. Finally, we obtained a series of 220 corrections (dX and dY) out of 23619 normal points. The obtained CPO series is plotted in Fig. 4, and details after 1985 are presented in Fig. 5.

Fig. 4. dX and dY residuals of LLR observations in mas, using separated windows of 70 days.

Fig. 5. Details after 1985 of dX and dY residuals from LLR observations, using separated windows of 70 days.

4.4. Analyses of LLR residuals

We present estimates using the empirical models 2 and 3 for the long-period components of nutation after the FCN was removed (cf. Sect. 3 for VLBI data). Results of the previous works are listed for comparison in Table 3. H18 presented the results as in-phase and out-of-phase terms of dψ and dϵ, and we converted them into dψ sin ϵ and dϵ. The decrease in WRMS after the weighted least-square fits is clear (about 20%), and the largest deviations of all nutation terms occur for the 18.6-year sine term, in accordance with the results of ZC09 and H18. The fit coefficients show no evident correction to the precession rate, but a deviation of about −0.3 mas for all of the coefficients of the 18.6-year nutation term. The determined accuracies are better than those in ZC09 because of ten more years of observational time coverage. Furthermore, Pavlov et al. (2016) treated the LLR data meticulously, including a discussion of the limits of the widely used IERS C04 series (see details in Sect. 4 in their paper), providing us with a O−C series of high quality, which also contributes to the better accuracy.

H18 presented an overall estimation with LLR data from 1970 to 2016. Five nutation terms were estimated along with the station coordinates, velocities, and retro-reflector coordinates, and thus uncertainties were affected by the correlations. For instance, correlations were up to 0.5 between nutation coefficients and the perturbation rotation rate in y-direction, up to 0.4 between the sine term of dψ and cosine term of dϵ with the z-component of the retro-reflector coordinates, and up to 0.7 between the cosine term of dψ and sine term of dϵ with the initial z-components of the lunar position. In our method, however, these values are fixed by applying the ephemeris DE430 and the coordinates of the stations and retro-reflectors estimated by Pavlov (2018, priv. comm.).

In our results, the fit coefficients of the 18.6-year nutation term with and without the 9.3-year term are not consistent. This feature is different from the results obtained by VLBI analyses. Meanwhile, the correlation coefficients between two nutation terms are over 0.7, revealing the incapability of LLR data to separate the components effectively. This is probably because that LLR observations are directly related with the motions of the Moon, which is also the most important excitation of the 18.6-year and 9.3-year nutations. Nevertheless, the correlation coefficients between the secular term and the nutation terms are generally smaller than those of VLBI, probably benefiting from the longer time span of ten years. Furthermore, the correlation coefficient between the secular term and the sine term of 18.6 yr nutation remains larger than its counterparts, which may reveal a common problem shared by the VLBI and LLR techniques.

4.5. Influence of not separating the two ways of light travel in LLR calculation

As explained in Sect. 4.2, we separated the two ways of light travel when we converted LLR residuals into CPO. We tested the effect of this factor by using only the emission time when we calculated the partial derivatives (this equals the first iteration in the calculation process used in this work), and placed the external residuals (Pavlov et al. 2016) into the least-squares fit to obtain dX and dY. Then we used this CPO series to estimate the same corrections for the precession and nutation terms. The results only differed by 10⁻⁵ of their formal errors compared to the results in Sect. 4.4. This indicates that the Earth rotation parameters are not as sensitive to the accuracies of reflection and reception times as the laser ranging itself.

5. CPO analyses with optical observations

The history of optical observations can be traced back to the nineteenth century: this can be considered as an advantage in determining long-term effects in the Earth rotation. We used the OA00 series constructed by Vondrák & Štefka (2010), which is based on a combination of the HIPPARCOS and Tycho catalogs and long-lasting ground-based astrometric observations. It was derived from the Earth orientation catalogs (EOC) containing 4418 celestial objects. The catalogs were used to derive the CPO during the time interval 1899.7–1992.0, referred to as the HIPPARCOS celestial reference frame (HCRF).

5.1. Transformation of optical data

The OA00 series provides time series of nutation offsets dψ and dϵ with respect to the IAU 1976/1980 precession-nutation models (Lieske et al. 1977; Wahr 1981; Seidelmann 1982) between 1899 and 1972. Therefore, it is necessary to transform them into the same system as the VLBI series, that is, dX and dY with respect to the IAU 2006/2000 precession-nutation models.

We followed the method described by Capitaine & Wallace (2006). First we evaluated the observed nutation angles of date by adding the optical nutation offsets to the IAU 1980 nutation model,

$$\begin{array}{c}\hfill \mathrm{\Delta }\psi =\mathrm{\Delta }{\psi }_{\mathrm{IAU}1980}+\mathrm{d}\psi ,\mathrm{\Delta }ϵ=\mathrm{\Delta }{ϵ}_{\mathrm{IAU}1980}+\mathrm{d}ϵ.\end{array}$$(8)

Then we calculated the nutation matrix with Δψ, Δϵ, and the mean obliquity of date (ϵ_A) by

$$\begin{array}{c}\hfill \mathcal{N}={\mathcal{R}}_1(-{ϵ}_A){\mathcal{R}}_3(\mathrm{\Delta }\psi ){\mathcal{R}}_1({ϵ}_A+\mathrm{\Delta }ϵ).\end{array}$$(9)

We applied the precession matrix (IAU 1976) of date to form the precession-nutation matrix 𝓜_obs, of which the elements (3,1) and (3,2) are the observed X and Y coordinates of the CIP in the GCRS. Finally, we obtained the celestial pole offsets with respect to the IAU 2006/2000 model:

$$\begin{array}{c}\hfill \mathrm{d}X={\mathcal{M}}_{\mathrm{obs}}(3,1)-X_{\mathrm{IAU}2006/2000},\\ \hfill \mathrm{d}Y={\mathcal{M}}_{\mathrm{obs}}(3,2)-Y_{\mathrm{IAU}2006/2000}.\end{array}$$(10)

The obtained CPO in arcseconds is plotted in Fig. 6.

Fig. 6. dX and dY residuals from optical observations.

5.2. Analyses of the optical residuals

Because the OA00 series has a resolution of 5 days (Vondrák & Ron 2000), it is theoretically able to determine the long-periodic nutation terms. Vondrák & Ron (2000) have performed certain analyses of the residual series OA97 and OA99, the predecessors of OA00 and OA10, at the end of the 1990s, and compared the results to VLBI observations back then. According to their results, the fit results of OA99 were closer to those of the VLBI series than OA97. The accuracies of fit coefficients of OA99 and VLBI were on the same order owing to the low quality of VLBI data before 1990. Since then, VLBI observations have lasted 20 more years and have been greatly improved. The corresponding accuracy of the fit coefficients of corrections to long-period terms in precession-nutation models has reached μas level, leaving the optical data an only advantage of the long history. The uncertainties of the CPO derived from OA00 series are 11.33 and 7.66 mas (weighted average) in dX and dY, respectively, about two hundred times that of VLBI data.

Models 1 and 2 were fit to the transformed CPO to estimate the possible long-term corrections. The results are presented in Table 4. In both models t refers to the number of centuries from a reference epoch J1956.0, which is almost the central epoch of the time coverage of optical observations. In contrast to the VLBI results, there is no significant difference between the constant and secular terms in these two tables, and the correlation coefficients between secular term and nutation terms are negligible. This indicates that a time span long enough (over ninety years) may be sufficient to separate the quadratic term and the 18.6-year nutation term, whereas optical data cannot help in improving the precession-nutation model today because of the poor accuracy.

6. Joint analyses of VLBI and LLR residuals

Because the observing history of VLBI is short and the optical and LLR data lack high accuracy, a combined series would benefit from the advantages of each technique while avoiding shortcomings. However, the quality of the optical data is far too poor compared with the other two techniques and therefore can barely contribute anything. Therefore, we only performed joint analyses of VLBI (OPA series) and LLR observations.

6.1. Systematic deviations between LLR and VLBI systems

It should be noted that the reference frame for LLR observations is dynamically defined and the reference frame for VLBI is kinematic, so that it is necessary to take into account the small systematic deviations between these two reference systems before the two series are combined. Liu et al. (2012) has developed the relationships between relative orientation offset with a rigid-body rotation (ϵ_x, ϵ_y, ϵ_z) for two reference frames and the difference of EOPs (ΔX, ΔY) expressed in these two frames:

$$\begin{array}{cc}& \mathrm{\Delta }X=+{ϵ}_{\mathrm{z}}Y-{ϵ}_{\mathrm{y}}Z\hfill \\ \hfill & \mathrm{\Delta }Y=-{ϵ}_{\mathrm{z}}X+{ϵ}_{\mathrm{x}}Z,\hfill \end{array}$$(11)

in which ϵ_x, ϵ_y, and ϵ_z are the rotation angles about the three axes. In the case of short time-intervals (several hundred years), the CIP coordinates X and Y are small quantities, and thus $Z=\sqrt{1-X^2-Y^2}\approx 1$, so that the effect on EOP can be limited to the first order in X, Y, and ϵ_x, y, z:

$$\begin{array}{cc}& \mathrm{\Delta }X\approx -{ϵ}_y\hfill \\ \hfill & \mathrm{\Delta }Y\approx +{ϵ}_x.\hfill \end{array}$$(12)

6.2. Possibilities of obtaining more CPO corrections out of LLR observations

With this method, we obtained the upper limit of the bias between the dynamical and kinematic reference frames by calculating weighted means of the differences between CPO series obtained from VLBI and LLR observations. The uncertainties were derived as the average of

$$\begin{array}{c}\hfill \sigma =\sqrt{{\sigma }_{\mathrm{LLR}}^2+{\sigma }_{\mathrm{VLBI}}^2},\end{array}$$(13)

representing the accuracy of LLR and VLBI observations. The results are

$$\begin{array}{cc}& {ϵ}_x=16.42\pm 210.50\mu \mathrm{as}\hfill \\ \hfill & {ϵ}_y=-148.37\pm 195.02\mu \mathrm{as}.\hfill \end{array}$$(14)

Because the derived uncertainties are larger than the biases themselves, we consider the two systems to be consistent.

6.3. Details of the combined series

The method of combination is almost the same as that applied in Zerhouni & Capitaine (2009) (Sect. 7): obtaining weighted means within sliding windows of 70 days and sliding every 10 data points, after removing the FCN from the opa2018a and LLR series. We obtained a combined series containing 655 points with this method. We also obtained a weighted mean series of VLBI (also with 70-day sliding windows) for comparison, and it contains 633 points. The combined series is plotted together with the series of mean CPO from VLBI in each window (see Fig. 7).

Fig. 7. Combined series of VLBI and LLR observations (blue) and the mean VLBI series obtained with 70-day sliding windows (red).

The features of the combined series are almost the same as that of the smoothed VLBI series, considering the better quality of the VLBI data. The differences in the figure indicate that the LLR CPO series has a comparable accuracy in certain periods (e.g., 1993–1995 and 2007–2010) that possibly reveals different details. The fitting results to all three models are listed in Table. 5. The coefficients are slightly different from those in Table. 1. The WRMS of the combined series reduces due to data smoothing in the combining method.

In Sect. 4 we obtained only 220 corrections from 23 619 normal points. We adoptedd two methods to explore the possibility of obtaining more corrections from the LLR normal points.

7. Discussion

First, because the observations are much more frequent after 1985 (see Fig. 1), it is possible to obtain more CPO data points out of the LLR normal points from this period of time by shortening the interval within which one point of the CPO (dX and dY) is derived. We decided to choose the length of the window according to the frequency of observation. We basically fixed the number of normal points to be 50 in each window, and the longest time span of the windows was set to be 70 days. By doing so, we obtained 67 windows that contained fewer than 50 normal points in periods without observations. The total number of windows is 483, more than twice the series presented in Sect. 4 (see in Fig. 8. Details after 1985 are shown in Fig. 10).

Another way to obtain more CPO data points is using a sliding window instead of separated windows. We took 70 days as the window length as in Sect. 4 and caused the window to slide at every 10 normal points. We obtadin 2513 corrections with this method, and the accuracy appears to be as good as that of the method using separated windows in Sect. 4 (see in Fig. 9, and details after 1985 are shown in Fig. 11).

Fig. 8. dX and dY residuals from LLR observations, using changing windows of 50 counts and within 70 days.

Fig. 9. dX and dY residuals from LLR observations, using sliding windows of 70 days.

Fig. 10. Details after 1985 of dX and dY residuals from LLR observations, using changing windows of 50 counts and within 70 days. (While it is stated that setting the number of normal points to 50, it appears that there may be 67 windows (out of 483, i.e., more than 10%) that contain fewer than 50 points, in which case the window length is about 70 days.)

Fig. 11. Details after 1985 of dX and dY residuals from LLR observations, using sliding windows of 70 days.

8. Concluding remarks

The fitting results using model 3 are presented in Table 6. The sliding-window method reduces the WRMS of the series, and the estimated deviations of all coefficients are smaller than those in Table 3. The changing-window method results in a lower accuracy in the CPO, but not in the fit coefficients. Moreover, the weighted means of the uncertainties of the two CPO series are 0.051 mas for the sliding-window method and 0.112 mas for the changing-window method, respectively. These results are within expectations. Observations are both more frequent and more accurate after 1985, so that with the sliding-window method we obtain a larger portion of the CPO of smaller uncertainties and thus reduce the weighted average. As for the changing-window method, the window durations are different while many parameters change with time in the calculation process. According to these two tests, the limited quantity of obtained CPO corrections is a cause of the estimate uncertainties of the nutation terms. Moreover, the time distribution of the observations and the uncertainties of the models that were used in the calculation process are clearly not perfect.

8.1. Future improvements of LLR observations for a better determination of the CPO

The accuracy of the CPO determined from LLR observations can be affected by many aspects. Of these, the observational error and frequency are most directly related to the observation itself. We show their changes with time in Fig. 12. The observational accuracy, although it suffers from instability, is improving. By comparing Figs. 4 and 12, especially between 2000 and 2010, we can conclude that the dispersion and larger uncertainties of the CPO in this period are the consequences of the lower frequency of the observations. Therefore, making LLR observations regular and sufficiently frequent to achieve a more uniform time distribution of normal points is quite essential in the future, before other necessary developments of related theories are possible.

Fig. 12. Time series of observational errors and number of used normal points; the second figure shows details of the first one. Each point or bin represents the weighted average of the total within a year.

We presented comparative studies using VLBI, LLR, and optical data to estimate a possible improvement of the current precession-nutation models in terms of long-period nutation. VLBI, LLR, and optical data were first treated separately and then combined (LLR and VLBI only) for further discussion. The weighted average of the standard error in the LLR series is comparable with that of VLBI series, and the number rises to over two hundred in the optical series.

The VLBI data present a most reliable feature of the CPO series, and are still the most important data for determining the precession-nutation of Earth, whereas the formal errors may be underestimated. The fit coefficients of the secular terms with the opa2018a series show an underestimation of the precession rate in dX by about 0.3 mas cy⁻¹. All fit coefficients of the 18.6-year nutation term revealed an underestimate of the amplitude by about 35 μas.

The correlation between 18.6-year and 9.3-year nutation terms in LLR fitting results is very strong because of the direct relation between LLR observations and the motions of the Moon. The fit coefficients show no evident correction to the precession rate, but a deviation of about −0.3 mas for all of the coefficients of the 18.6-year nutation term.

The differences between the details of the combined series and the smoothed VLBI series indicate that LLR has a comparable accuracy with VLBI in its best state. The differences of the fit coefficients are almost the same as the deviations.

The potential of LLR for determining the CPO, and therefore the long-period part of precession-nutation, is confirmed, while the instability of the observational accuracy and the irregularity of observations, in particular, limit the consistency of the accuracies, as well as the time resolution of the derived CPO series. Furthermore, the non-negligible correlations between the fit terms remind us that LLR is not yet dependable enough to contribute significantly to the improvement of the current precession-nutation model.

In terms of the optical series, different from the results of VLBI series, the 18.6-year nutation term is not statistically significant. Considering the poor accuracy of optical observations, it is not reasonable to use optical observations to evaluate the contemporary precession-nutation models.

In short, corrections to nutation terms in current models are suggested by VLBI observations, whereas the formal errors may be underestimated. Meanwhile, the LLR technique can be a candidate to provide an independent check for the determination of EOP after major improvements in observation regularity. Historical optical data are no longer helpful in spite of their long history⁵.

Acknowledgments

The 13th version of software routines from the IAU SOFA Collection (Hohenkerk 2012) was used to calculate the precession-nutation matrix according to the IAU 1976/1980 and IAU2006/2000A models, and to convert timescales as well. We are grateful to D. A. Pavlov and his colleagues for their help and valuable advice in analyzing LLR data, and Niu Liu and Ping-Jie Ding for their suggestions to improve the manuscript. This research is funded by the National Natural Science Foundation of China (NSFC) No. 11473013 and No. 11833004. We also appreciate the detailed and valuable advice of the anonymous referee to improve this work.

References

Appendix A: Correlation coefficients

Correlation coefficients of all least-squares fits are listed for reference.

All Tables

All Figures

	Fig. 1. Time distributions of LLR, VLBI (opa2018a), and optical CPO series. The time span of every CPO series is LLR 1970.6–2017.7, VLBI (opa2018a) 1979.9–2018.5, and optical (OA00) 1899.7–1991.9.
In the text

	Fig. 2. VLBI celestial pole offsets from the OPA analysis center with respect to the IAU 2006/2000 models.
In the text

	Fig. 3. dX and dY residuals of the OPA VLBI time series with respect to the IAU 2006/2000 models. FCN is removed in the updated IERS model.
In the text

	Fig. 4. dX and dY residuals of LLR observations in mas, using separated windows of 70 days.
In the text

	Fig. 5. Details after 1985 of dX and dY residuals from LLR observations, using separated windows of 70 days.
In the text

	Fig. 6. dX and dY residuals from optical observations.
In the text

	Fig. 7. Combined series of VLBI and LLR observations (blue) and the mean VLBI series obtained with 70-day sliding windows (red).
In the text

	Fig. 8. dX and dY residuals from LLR observations, using changing windows of 50 counts and within 70 days.
In the text

	Fig. 9. dX and dY residuals from LLR observations, using sliding windows of 70 days.
In the text

	Fig. 10. Details after 1985 of dX and dY residuals from LLR observations, using changing windows of 50 counts and within 70 days. (While it is stated that setting the number of normal points to 50, it appears that there may be 67 windows (out of 483, i.e., more than 10%) that contain fewer than 50 points, in which case the window length is about 70 days.)
In the text

	Fig. 11. Details after 1985 of dX and dY residuals from LLR observations, using sliding windows of 70 days.
In the text

	Fig. 12. Time series of observational errors and number of used normal points; the second figure shows details of the first one. Each point or bin represents the weighted average of the total within a year.
In the text