Issue 
A&A
Volume 627, July 2019



Article Number  A81  
Number of page(s)  10  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/201834785  
Published online  03 July 2019 
Analyses of celestial pole offsets with VLBI, LLR, and optical observations
School of Astronomy and Space Science, Key Laboratory of Modern Astronomy and Astrophysics (Ministry of Education), Nanjing University, 163 Xianlin Avenue, 210023 Nanjing, PR China
email: yuting_cheng@smail.nju.edu.cn, jcliu@nju.edu.cn
Received:
6
December
2018
Accepted:
8
May
2019
Aims. This work aims to explore the possibilities of determining the longperiod part of the precessionnutation of the Earth with techniques other than very long baseline interferometry (VLBI). Lunar laser ranging (LLR) is chosen for its relatively high accuracy and long period. Results of previous studies could be updated using the latest data with generally higher quality, which would also add ten years to the total time span. Historical optical data are also analyzed for their rather long timecoverage to determine whether it is possible to improve the current Earth precessionnutation model.
Methods. Celestial pole offsets (CPO) series were obtained from LLR and optical observations and were analyzed separately by weighted leastsquare fits of three empirical models, including a quadratic model, a linear term plus an 18.6year nutation term, and a linear term plus two nutation terms with 18.6year and 9.3year periods. Joint analyses of VLBI and LLR data is also presented for further discussion.
Results. We improved th determination of the nutation terms with both VLBI and LLR data. The VLBI data present a most reliable feature of the CPO series with the highest accuracy, and they are most important for determining the precessionnutation of the Earth. The standard errors of CPO obtained from the LLR technique have reached a level of several tens of microarcseconds after 2007, but they are probably underestimated because the models used in the calculation procedure are not perfect. Nevertheless, the poor time resolution of LLR CPO series is also a disadvantage. However, this indicates that LLR has the potential to determine celestial pole offsets with a comparably high accuracy with VLBI in the future and to serve as an independent check for the VLBI results. The current situation of LLR observations is also analyzed to provide suggestions of future improvement. The typical standard error of CPO series from historic optical observations is about two hundred times larger than that of the VLBI series and can therefore hardly contribute to the contemporary precessionnutation theory.
Key words: astrometry / reference systems / Moon / ephemerides
© ESO 2019
1. Introduction
Earth precessionnutation models describe the longterm and longperiod changes of the celestial intermediate pole (CIP) direction in the geocentric celestial reference system (GCRS). Before 2000, the IAU 1976/1980 precessionnutation 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 spacegeodetic 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 spacegeodetic 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^{−1} 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 Earthrelated 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 groundbased astrometric observations. With robust determinations of the proper motions of a binary star system from 4.5 million individual observations, the latest version, the EOC4 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 precessionnutation 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 longperiod part of the precessionnutation 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:

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

A linear term and 18.6year nutation term:

A linear term, and 18.6year and 9.3year nutation terms:
where Ω_{1} = 2π/0.186 and Ω_{2} = 2π/0.093, representing the arguments of the two principle nutation terms. The 18.6year 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 precessionnutation 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)^{1} (Fig. 2) and also the series derived by the Goddard Space Flight Center (GSF)^{2} (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 quasiperiodic 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^{3}), 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.6year nutation (0.4) of model 3 stands out among others. The short timespan of the VLBI observations (38 years) leads to an evident correlation between the 18.6year nutation term and longperiod terms (t^{1} and t^{2}). All the correlation coefficients can be found in Appendix A.
Weighted fits of three models to VLBI residuals (1979.6–2018.6), corresponding to the IAU 2006/2000 model.
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.6year term with and without the 9.3year 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^{−1}. Capitaine et al. (2009) have pointed out that the 18.6year nutation term is the most sensitive to the error in the precessionnutation model, with VLBI data up to 2008. With the accumulated highprecision data, we here obtain more consistent results of the corrections of the 18.6year nutation term. All fit coefficients of the 18.6year 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 socalled 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 roundtrip times and then the residuals of the roundtrip time [observation minus calculation (O−C)], which can be converted into residuals in oneway 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 oneway distance in centimeters) provided by Pavlov et al. (2016) ^{4}. 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 groundbased 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.
Basic information of LLR data of each station.
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 roundtrip time may cause errors of hundreds of meters in residuals of the oneway 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 roundtrip time as Δt = (D_{1} + D_{2})/c instead of Δt = 2D/c, in which D_{1} and D_{2} are the stationreflector vectors of the forward and backward ways of light travel in the barycentric celestial reference system (BCRS).
The stationreflector vector can be written as the summation of three vectors, as presented in ZC09,
where −SE is the geocentric position of the groundbased 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:
In the expression above, only 𝓠, the transformation matrix for precession and nutation, is related to X and Y:
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
in which t_{1} is the emission time, t_{2} is the reflection time, and t_{3} is the reception time. We follow Sect. 4 in Pavlov et al. (2016) to calculate t_{2} and t_{3}, 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 oneway 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 oneway distance with the speed of light in vacuum. Then we obtained the CPO by weighted leastsquare fits using
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 oneway distance. The weighted mean of the resulting differences in the oneway 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 oneway distance with a factor three before we obtained the CPO with weighted leastsquare 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 longperiod 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 inphase and outofphase terms of dψ and dϵ, and we converted them into dψ sin ϵ and dϵ. The decrease in WRMS after the weighted leastsquare fits is clear (about 20%), and the largest deviations of all nutation terms occur for the 18.6year 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.6year 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.
Weighted fits of models 2 and 3 to LLR residuals.
H18 presented an overall estimation with LLR data from 1970 to 2016. Five nutation terms were estimated along with the station coordinates, velocities, and retroreflector 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 ydirection, up to 0.4 between the sine term of dψ and cosine term of dϵ with the zcomponent of the retroreflector coordinates, and up to 0.7 between the cosine term of dψ and sine term of dϵ with the initial zcomponents of the lunar position. In our method, however, these values are fixed by applying the ephemeris DE430 and the coordinates of the stations and retroreflectors estimated by Pavlov (2018, priv. comm.).
In our results, the fit coefficients of the 18.6year nutation term with and without the 9.3year 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.6year and 9.3year 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 leastsquares 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^{−5} 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 longterm 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 longlasting groundbased 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 precessionnutation 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 precessionnutation 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,
Then we calculated the nutation matrix with Δψ, Δϵ, and the mean obliquity of date (ϵ_{A}) by
We applied the precession matrix (IAU 1976) of date to form the precessionnutation 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:
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 longperiodic 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 longperiod terms in precessionnutation 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 longterm 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.6year nutation term, whereas optical data cannot help in improving the precessionnutation model today because of the poor accuracy.
Weighted fits of models 1 and 2 to optical residuals (1899.71991.9) corresponding to the IAU 2006/2000 model.
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 rigidbody rotation (ϵ_{x}, ϵ_{y}, ϵ_{z}) for two reference frames and the difference of EOPs (ΔX, ΔY) expressed in these two frames:
in which ϵ_{x}, ϵ_{y}, and ϵ_{z} are the rotation angles about the three axes. In the case of short timeintervals (several hundred years), the CIP coordinates X and Y are small quantities, and thus , so that the effect on EOP can be limited to the first order in X, Y, and ϵ_{x, y, z}:
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
representing the accuracy of LLR and VLBI observations. The results are
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 70day 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 70day 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.
Weighted fits of the joint series (1970.7–2018.5) to three models, corresponding to the IAU 2006/2000 model.
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 slidingwindow method reduces the WRMS of the series, and the estimated deviations of all coefficients are smaller than those in Table 3. The changingwindow 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 slidingwindow method and 0.112 mas for the changingwindow method, respectively. These results are within expectations. Observations are both more frequent and more accurate after 1985, so that with the slidingwindow method we obtain a larger portion of the CPO of smaller uncertainties and thus reduce the weighted average. As for the changingwindow 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.
Weighted fits of model 3 to LLR residuals obtained with slidingwindow method and changingwindow method, corresponding to the IAU 2006/2000 models.
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 precessionnutation models in terms of longperiod 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 precessionnutation 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^{−1}. All fit coefficients of the 18.6year nutation term revealed an underestimate of the amplitude by about 35 μas.
The correlation between 18.6year and 9.3year 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.6year 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 longperiod part of precessionnutation, 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 nonnegligible correlations between the fit terms remind us that LLR is not yet dependable enough to contribute significantly to the improvement of the current precessionnutation model.
In terms of the optical series, different from the results of VLBI series, the 18.6year 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 precessionnutation 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^{5}.
Acknowledgments
The 13th version of software routines from the IAU SOFA Collection (Hohenkerk 2012) was used to calculate the precessionnutation 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 PingJie 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
 Capitaine, N., & Wallace, P. T. 2006, A&A, 450, 855 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Capitaine, N., Wallace, P. T., & Chapront, J. 2003, A&A, 412, 567 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Capitaine, N., Mathews, P. M., Dehant, V., Wallace, P. T., & Lambert, S. B. 2009, Celest. Mech. Dyn. Astron., 103, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Chapront, J., ChaprontTouzé, M., & Francou, G. 1999, A&A, 343, 624 [NASA ADS] [Google Scholar]
 Fey, A. L., Gordon, D., Jacobs, C. S., et al. 2015, AJ, 150, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., & Kuchynka, P. 2014, Interplanet. Netw. Progr. Rep., 196, 1 [Google Scholar]
 Herring, T. A., Mathews, P. M., & Buffett, B. A. 2002, J. Geophys. Res. (Solid Earth), 107, 2069 [NASA ADS] [CrossRef] [Google Scholar]
 Hofmann, F., Biskupek, L., & Müller, J. 2018, J. Geodesy, 92, 975 [NASA ADS] [CrossRef] [Google Scholar]
 Hohenkerk, C. Y. 2012, in Journées Systèmes de Référence Spatiotemporels 2011, eds. H. Schuh, S. Boehm, T. Nilsson, & N. Capitaine, 21 [Google Scholar]
 Kopeikin, S. M. 1990, Sov. Astron., 34, 5 [NASA ADS] [Google Scholar]
 Lieske, J. H., Lederle, T., Fricke, W., & Morando, B. 1977, A&A, 58, 1 [NASA ADS] [Google Scholar]
 Liu, J.C., Capitaine, N., Lambert, S. B., Malkin, Z., & Zhu, Z. 2012, A&A, 548, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marini, J. W., & Murray, C. W. 1973, Correction of laser range tracking data for atmospheric refraction at elevation above 10 degrees, http://hdl.handle.net/2060/19740007037 [Google Scholar]
 Mathews, P. M., Herring, T. A., & Buffett, B. A. 2002, J. Geophys. Res. (Solid Earth), 107, 2068 [NASA ADS] [CrossRef] [Google Scholar]
 Mendes, V. B., & Pavlis, E. C. 2003, in Proceedings of the 13th International Laser Ranging Workshop, Washington D.C., eds. R. Noomen, S. Klosko, C. Noll, & M. Pearlman, NASA/CP2003212248 [Google Scholar]
 Mendes, V. B., & Pavlis, E. C. 2004, Geophys. Res. Lett., 31, L14602 [NASA ADS] [CrossRef] [Google Scholar]
 Mendes, V. B., Prates, G., Pavlis, E. C., Pavlis, D. E., & Langley, R. B. 2002, Geophys. Res. Lett., 29, 1414 [NASA ADS] [CrossRef] [Google Scholar]
 Nothnagel, A., Artz, T., Behrend, D., & Malkin, Z. 2017, J. Geodesy, 91, 711 [Google Scholar]
 Pavlov, D. A., Williams, J. G., & Suvorkin, V. V. 2016, Celest. Mech. Dyn. Astron., 126, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Petit, G., & Luzum, B. 2010, IERS Technical Note, 36 [Google Scholar]
 Seidelmann, P. K. 1982, Celest. Mech., 27, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Viswanathan, V., Fienga, A., Gastineau, M., & Laskar, J. 2017, Notes Scientifiques et Techniques de l’Institut de Mecanique Celeste, 108 [Google Scholar]
 Vondrák, J., & Ron, C. 2000, in IAU Colloq. 180: Towards Models and Constants for SubMicroarcsecond Astrometry, eds. K. J. Johnston, D. D. McCarthy, B. J. Luzum, & G. H. Kaplan, 248 [Google Scholar]
 Vondrák, J., & Štefka, V. 2010, A&A, 509, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wahr, J. M. 1981, Geophys. J., 64, 705 [NASA ADS] [CrossRef] [Google Scholar]
 Zerhouni, W., & Capitaine, N. 2009, A&A, 507, 1687 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Correlation coefficients
Correlation coefficients of all leastsquares fits are listed for reference.
Correlation coefficients of the VLBI fitting results in Sect. 3.
Correlation coefficients of the LLR fitting results in Sect. 4.
Correlation coefficients of the optical fitting results in Sect. 5.
Correlation coefficients of the fitting results of the combined series in Sect. 6.
Correlation coefficients of the LLR fitting results in Sect. 7, of the slidingwindow (the first part) and changingwindow (the second part) methods.
All Tables
Weighted fits of three models to VLBI residuals (1979.6–2018.6), corresponding to the IAU 2006/2000 model.
Weighted fits of models 1 and 2 to optical residuals (1899.71991.9) corresponding to the IAU 2006/2000 model.
Weighted fits of the joint series (1970.7–2018.5) to three models, corresponding to the IAU 2006/2000 model.
Weighted fits of model 3 to LLR residuals obtained with slidingwindow method and changingwindow method, corresponding to the IAU 2006/2000 models.
Correlation coefficients of the fitting results of the combined series in Sect. 6.
Correlation coefficients of the LLR fitting results in Sect. 7, of the slidingwindow (the first part) and changingwindow (the second part) methods.
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 70day 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.