J G Williams

Keywords: Moon, Physical Librations, spin-orbit synchronous resonance. N. Rambaux 1, 2

The Lunar Laser Ranging experiment has been active since 1969 when Apollo astronauts placed the first retroreflector on the Moon. The data accuracy of a few centimeters over recent decades, joined to a new numerically integrated ephemeris, DE421, encourages a new analysis of the lunar physical librations of that ephemeris, and especially the detection of three modes of free physical librations (longitude, latitude, and wobble modes). This analysis was performed by iterating a frequency analysis and linear least-squares fit of the wide spectrum of DE421 lunar physical librations. From this analysis we identified and estimated about 130 ∼ 140 terms in the angular series of latitude librations and polar coordinates, and 89 terms in the longitude angle. In this determination, we found the non-negligible amplitudes of the three modes of free physical libration. The determined amplitudes reach 1.296" in longitude (after correction of two close forcing terms), 0.032" in latitude and 8.183" × 3.306" for the wobble, with the respective periods of 1056.13 days, 8822.88 days (referred to the moving node), and 27257.27 days. The presence of such terms despite damping suggests the existence of some source of stimulation acting in geologically recent times.

Introduction

At first appearance, the Moon seems to be unchanging and rotating close to its minimum energy state. Analyses of Lunar Laser Ranging data reveal a lunar rotation not quite at minimum, even without the influence of oceans or a significant atmosphere. It has been recognized for three decades (Calame 1976a(Calame , 1976b[START_REF] Calame | Free librations of the Moon from Lunar Laser Ranging[END_REF] that the amplitudes of two of the free libration modes exceed 1 arcsec. The damping times for these two modes are estimated as 2 10 4 and 2 10 6 years [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF] for the longitude and wobble modes, respectively, so there must be some source of stimulation acting in geologically recent times. A third mode is much smaller and consequently more elusive. A clear determination is needed. As a first step toward a better understanding of lunar free librations, this paper makes an accurate determination of amplitudes, periods, and phases of the free libration modes along with the large amplitude forced librations.

There are three modes of free librations for a solid Moon. They correspond to the dynamical normal modes of the Moon in spin-orbit resonance when its spin is displaced from its dynamical equilibrium position. The period of the free rotational modes of the solid Moon can be calculated in linear first order theory (see e.g. [START_REF] Eckhardt | Theory of the libration of the moon[END_REF][START_REF] Moons | Analytical theory of the libration of the Moon[END_REF]; [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF][START_REF] Varadi | Normal modes of synchronous rotation[END_REF]. One is in longitude, with rotation parallel to the equatorial plane of the Moon, and has a period of 2.9 years. For the second mode, the latitude mode, the axis normal to the lunar equatorial plane traces out a small cone. That motion in space is retrograde with a period near 81 years. The third mode of free libration is related to the motion of the axis of figure about the rotation axis. As seen from the lunar body reference frame, the axis of rotation traces out an elliptical path with a period around 75 years. This motion is analogous to the Earth's Chandler wobble, but as seen from the inertial reference frame the period is around 27 days [START_REF] Newhall | Estimation of the Lunar Physical Librations[END_REF]. As emphasized in [START_REF] Bois | Proposed terminology for a general classification of rotational swing motions of the celestial solid bodies[END_REF], the free librations at 2.9 and 81 years are the result of the spin-orbit resonance problem of the Moon around the Earth and do not correspond to the free Eulerian modes of rotation of a body alone in space. While the wobble mode is analogous to the Chandler wobble, the expressions for the frequency and ratio of axes of the elliptical path are different by a factor of 2 so that mode is also modified by the synchronous spin-orbit motion of the Moon.

While the periods of the free modes can be calculated analytically, the determination of their amplitudes requires data analysis. The first valid determination of the free librations is due to Calame (1976aCalame ( , 1976b[START_REF] Calame | Free librations of the Moon from Lunar Laser Ranging[END_REF]. More recent determinations have been made by [START_REF] Jin | Determination of some physical parameters of the Moon with lunar laser ranging data[END_REF], [START_REF] Newhall | Estimation of the Lunar Physical Librations[END_REF] and [START_REF] Chapront | Complements to Moon's Lunar Libration Theory[END_REF]. Calame (1976a[START_REF] Calame | Free librations of the Moon from Lunar Laser Ranging[END_REF] fit Lunar Laser Ranging data using both an analytical series and a JPL numerical integration. [START_REF] Jin | Determination of some physical parameters of the Moon with lunar laser ranging data[END_REF] fit Lunar Laser Ranging data. [START_REF] Chapront | Complements to Moon's Lunar Libration Theory[END_REF] and [START_REF] Newhall | Estimation of the Lunar Physical Librations[END_REF] determined the amplitudes and phases by comparing an analytical representation to the Euler angles in a JPL ephemeris. Therefore, the free librations depend on the orbital and lunar physical model and the method of fit. The models and fits have improved with time and we expect the free librations embedded in the Euler angles to also improve.

Several later determinations of the amplitude, period and phase of the librations from a numerical ephemeris [START_REF] Newhall | Estimation of the Lunar Physical Librations[END_REF][START_REF] Chapront | Complements to Moon's Lunar Libration Theory[END_REF][START_REF] Chapront | The lunar libration: comparisons between various models -a model fitted to LLR observations[END_REF] are done using DE403. The DE403 integration was done in 1995 and used a tidally flexing solid Moon. The present paper extracts the physical librations from the modern lunar Euler angles of DE421. The lunar orbit and physical librations are significantly improved in DE421 relative to DE403 (Williams et al. 2008a;[START_REF] Folkner | The Planetary and Lunar Ephemeris DE 421[END_REF]. DE421 numerically integrated the rotation of both lunar solid mantle and fluid core. Lunar geophysical parameters for tides, fluid core, core flattening, dissipation and gravity field were fit or input. To generate the lunar ephemeris of DE421, 38 yr of Lunar Laser Ranging data were fit with a weighted rms residual of less than 2 cm at the end (Williams et al. 2008a;[START_REF] Folkner | The Planetary and Lunar Ephemeris DE 421[END_REF]). As a first step toward understanding the origin and evolution of lunar librations, this paper focuses on extracting the free librations from the recent DE421 Euler angles.

In this paper, we follow the method of Newhall & Williams (1997) improved by some refinenements described in Section 4. Before that Section, we describe briefly the geometry of the physical librations and then the calculation of the free periods. In Section 5, we present the physical libration series resulting from our new analysis focusing on the free librations and finally we conclude.

Geometry of the physical librations

The physical librations are small and are extracted from the three Euler angles describing the orientation and spin rotation of the Moon. The Euler angles come from the recent JPL ephemeris called DE421 [START_REF] Folkner | The Planetary and Lunar Ephemeris DE 421[END_REF]Williams et al. 2008a). This accurate ephemeris is obtained by jointly fitting a numerical model of the lunar rotation along with orbits of the Moon, the Sun, the Earth and planets. The dynamical model contains relativistic Earth-Moon interactions, gravitational harmonics for the Moon (up to degree 4), Earth (zonal), and Sun (J 2 ), the tides on Earth and Moon, and an oblate fluid lunar core. The lunar orbit and rotation are fit to Lunar Laser Ranging data from 1970 to 2007 with adjustable parameters for moments of inertia, gravity field, tides, dissipation, interaction with a fluid core, and initial conditions for both solid mantle and fluid core [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF](Williams et al. , 2008b)). For recent years the lunar range residuals are less than 2 centimeters and the Moon's Euler angles are good to a few milliarcseconds. Periodic terms well separated in frequency can typically be determined to about 1 milliarcsecond [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF].

The three Euler angles (ψ, θ, ϕ) express the transformation between the mean Earth equatorial reference frame with equinox at J2000 (OXY Z) and the selenographic reference frame (Oxyz), where O is the lunar center of mass. The Euler angles are ψ the angle along the J2000 mean Earth equatorial plane, θ the nutation angle that expresses the inclination of the figure axis Oz with respect to the normal to the mean Earth equator OZ and ϕ the rotation angle expressed along the lunar equatorial plane.

The libration angles are designated by Iσ, ρ, and τ and are defined as small corrections to the lunar Cassini laws [START_REF] Eckhardt | Theory of the libration of the moon[END_REF][START_REF] Newhall | Estimation of the Lunar Physical Librations[END_REF]). However, the libration angles are defined with respect to the mean ecliptic of date and not the mean Earth equatorial reference frame with equinox at J2000 as in the numerical ephemeris. Thus, the Euler angles (ψ, θ, ϕ) expressed in the J2000 equatorial Earth reference system were transformed to the Euler angles (ψ ′ , θ ′ , ϕ ′ ) expressed with respect to the ecliptic of date by applying the four successive rotation matrices

R 3 (-Λ A )R 1 (π A )R 3 (Π A )R 1 (ǫ)
, where R i denotes the rotation matrix around an x-axis if i = 1, and around a z-axis if i = 3. The angles Λ A , π A , Π A , ǫ represent the ecliptic longitude from the mean equinox of date along the moving ecliptic to the intersection with the fixed ecliptic, the inclination of the moving ecliptic on the fixed ecliptic, the longitude of the ascending node of the moving ecliptic of date upon fixed ecliptic, and the Earth's obliquity at J2000, respectively [START_REF] Lieske | Expressions for the precession quantities based upon the IAU /1976/ system of astronomical constants[END_REF][START_REF] Simon | Numerical expressions for precession formulae and mean elements for the Moon and the planets[END_REF][START_REF] Hilton | Report of the International Astronomical Union division I Working Group on precession and the ecliptic[END_REF]). Thus, the libration angles are defined as

θ ′ = I + ρ ψ ′ = Ω + σ ψ ′ + ϕ ′ = π + λ + τ, (1)
where Ω is the ascending node of the mean orbit with respect to the mean ecliptic plane of date and λ is the mean longitude. I represents the tilt of the mean equator of the Moon to the mean ecliptic equal to about 1.543 degrees (see Figure 1). The tilt θ ′ of the equator to the ecliptic is in the opposite direction of the orbit inclination. So the angle θ ′ is small implying that the sum of ψ ′ and ϕ ′ is better determined than either angle alone. Consequently, σ is multiplied by I to get the smaller Iσ. τ is the libration in longitude and ρ and Iσ are referred to as the latitude librations. The physical librations in latitude are also well described by using the two components (lunar body frame) of the unit vector pointing towards the mean pole of the ecliptic of date on the lunar equatorial plane (p 1 , p 2 ) defined as

p 1 = -sin θ ′ sin ϕ ′ (2)
p 2 = -sin θ ′ cos ϕ ′ (3)
These two variables are coupled together (see Eq.7) and are often described in the analytical theories (see e.g. [START_REF] Eckhardt | Theory of the libration of the moon[END_REF][START_REF] Moons | Analytical theory of the libration of the Moon[END_REF][START_REF] Petrova | Analytical extension of lunar libration tables[END_REF]).

Determination of the free periods

The lunar physical librations have three free rotational modes: one in longitude and two for pole direction. The free librations come from the solutions of the rotational differential equations reduced by using the first Cassini law and linearized with respect to small quantities. The equations governing the rotation of a solid Moon are the classical Euler-Liouville equations

d[I]ω dt + ω × [I]ω = T (4)
where [I] is the variable tensor of inertia composed of a rigid part, a centrifugal part and a tidal part (see [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF]. ω is the rotation vector and T is the sum of the external torques acting on the Moon. The main torque is created by the pointmass Earth acting on the Moon's dynamical figure developed in spherical harmonics of second degree. It is expressed as

T = 3Gm r 3 u × [I]u (5)
with G the gravitational constant, m the mass of the Earth, and u the unit vector toward the Earth in the Moon's reference system. At first order, whithout fluid core, tidal distortion or dissipation, the homogeneous linearized equations governing the rotation of the spin-orbit resonant Moon are

p1 -(1 -β)n ṗ2 + 4βn 2 p 1 = 0 p2 + (1 -α)n ṗ1 + αn 2 p 2 = 0 (6) τ + 3γn 2 τ = 0
where p 1 , p 2 represent the pole position normal to the ecliptic plane in the Moon reference system and α

= (C -B)/A, β = (C -A)/B, γ = (B -A)/C, with A < B < C
the principal moments of inertia. n is the lunar mean motion. For the linearized case, the last differential equation is decoupled from the first two. The eigenfrequencies of this system are √ 3γn for the libration in longitude, 2 √ αβn for wobble and (1 + 3 2 β)n for the Moon referenced libration in latitude (equivalent to 3 2 βn for the space referenced frequency). The elastic deformation and presence of the core alter these periods by small amounts [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF]. The free periods correspond to the eigenfrequencies of the linear system.

The fitting method

The lunar rotation exhibits combined forced and free librations and we determine in the same fit the two types of librations. The libration angles and polar coordinates are shown in Figures 2 and3 for 1070 yr. The libration angles and polar coordinates are fit with a combination of polynomial and periodic terms. The periods of the forced libration components are known because they arise from periods present in the lunar orbit from (i) Earth-Moon-Sun effects, with periods related to the Delaunay arguments of lunar theory [START_REF] Chapront-Touzé | The lunar ephemeris ELP 2000[END_REF][START_REF] Chapront | A new determination of lunar orbital parameters, precession constant and tidal acceleration from LLR measurements[END_REF], and (ii) planetary effects [START_REF] Bretagnon | Theory for the motion of all the planets -The VSOP82 solution[END_REF][START_REF] Eckhardt | Planetary and earth figure perturbations in the librations of the Moon[END_REF]. The arguments of the physical libration terms φ i (t) are composed of Delaunay and planetary angles, expressed as time polynomials, and the frequencies are their derivatives.

The amplitudes of the periodic forced libration terms depend on both the strength of the torque applied to the aspherical Moon and how close the forcing period is to the resonant periods of the rotation, the free periods. The strength of a torque term can be variable in time, as for example due to the eccentricity of the Earth-Moon orbit around the Sun, which varies with time and alters the amplitude of the librations. We represent these variations by introducing Poisson terms, linear terms in time in the amplitudes. The analytical form of the librations is thus assumed to be composed of Fourier and Poisson series (with varying frequencies through the Delaunay arguments) and polynomials. The polynomial coefficients are used to represent small corrections to the uniform rotation of the Moon and very long period effects. Therefore the fit function takes the following form

f (t) = m j=1 a j t j + n i=1 [(C i + ǫ c i t) cos (φ i (t)) + (S i + ǫ s i t) sin (φ i (t))] (7)
where the adjustable coefficients are a j , the polynomial coefficients, (ǫ c i , ǫ s i ) the Poisson coefficients, and (C i , S i ) the Fourier series coefficients. The arguments φ i (t) are combination of Delaunay arguments, planetary angles, free libration terms of the form ω i t + α i , or simply ω i t for terms with unidentified origin. The polynomial expression of the Delaunay arguments come from the recent determination of [START_REF] Chapront | A new determination of lunar orbital parameters, precession constant and tidal acceleration from LLR measurements[END_REF], while for the planetary arguments we used [START_REF] Bretagnon | Theory for the motion of all the planets -The VSOP82 solution[END_REF]. The origin of the time t is 1 st January 2000 at 12h.

The fitting method is based on the sequential use of frequency analysis and leastsquares reduction methods. The frequency analysis had been developed by [START_REF] Laskar | Secular evolution of the Solar System over 10 million years[END_REF][START_REF] Laskar | Frequency Map analysis and quasi periodic decompositions[END_REF] for the use of Celestial Mechanics and it allows accurate determination of the frequencies and associated amplitudes contained in the problem. The least-squares method is required to take into account varying frequencies and amplitudes in the fit. Between the frequency analysis and the least-squares fit, the empirically derived frequencies are identified with the frequencies of the φ i (t) expressions. The combination of the two methods is iterative in order to converge toward robust values and to allow for non-linearity introduced by the fit for the free frequencies.

In the first step, we performed a frequency analysis directly on the three libration angles and two polar components. The initial frequency analysis provided a set of 30 frequencies for Iσ and ρ, and only 5 frequencies for τ , p 1 , and p 2 . The identification for the angle τ is more difficult because there are many planetary terms, including long period terms which are strongly correlated with the polynomial terms. We identified the empirical frequencies with the well-known frequencies of the lunar theory by using [START_REF] Eckhardt | Planetary and earth figure perturbations in the librations of the Moon[END_REF][START_REF] Chapront | A new determination of lunar orbital parameters, precession constant and tidal acceleration from LLR measurements[END_REF]. The angle I = 5553.6 arcseconds is large, so the leading libration term for both p 1 and p 2 presents a very large amplitude [START_REF] Eckhardt | Theory of the libration of the moon[END_REF]; Moon's 1984) and a non-negligible Poisson term.

The Poisson term generates spurious frequencies in the frequency analysis performed over the entire time-span of 1070 years. In order to identify these spurious terms, we analyze the time-series by using a slide frequency analysis over 260 years. A linear regression for the amplitude allows a first estimate of the Poisson coefficients and indicates if the Poisson term is necessary or not in the fit process. The second step consists of using the least-squares method to fit the formal expression [START_REF] Chapront | Complements to Moon's Lunar Libration Theory[END_REF] to the libration angles with adequate Poisson terms. The frequencies used in the function [START_REF] Chapront | Complements to Moon's Lunar Libration Theory[END_REF] have been determined in the first step. Then we compute the residuals for this solution.

In the third step, we apply the frequency analysis to the residuals of the previous step. We obtain 143 frequencies for Iσ, 134 frequencies for ρ, 89 frequencies for τ , 143 frequencies for p 1 , and 130 frequencies for p 2 . The improvement of the analysis by taking into account the polynomial and Poisson terms is obvious. The large number of frequencies improve the fit significantly over the previous determination of Newhall [START_REF] Newhall | Estimation of the Lunar Physical Librations[END_REF], which was based on a set of 30 frequencies. The longer time span (1070 years against 718 years) gives better frequency resolution.

In the fourth step, we search for the major Poisson terms, corresponding to an increase superior to 1 mas over one thousand years, and finally the last step consists in fitting by the least-squares method the formal solution composed of Fourier, polynomial, and Poisson terms to the residuals from step four. Then, we compute the residuals of this solution, shown in Figures 4 and5.

Physical libration series

In this section, we analyze the series obtained in the previous section that are available in supplementary material in electronic format where, we have listed the first terms of the series in Tables 1 to 10. The series are formed by polynomial (shown in Table 2), Fourier, and Poisson terms and the arguments are polynomials (see function 7). The physical librations could be classified in three different groups based on the origin of their periods. The first group is related to the main lunar problem, Sun-Earth-Moon, and the frequencies are combination of purely Delaunay arguments. For the second group, the precessing node and planetary interaction present arguments dependent on Ω, the precessing node, and the mean longitudes of the planets. The first two groups . The third group depends on lunar free frequencies, about 1056 days, 81 years and 75 years. The last group is the main objective of this paper and they will be described separately in following sections.

The physical librations can also be arranged in different groups based on the origin of the librations [START_REF] Bois | Proposed terminology for a general classification of rotational swing motions of the celestial solid bodies[END_REF]. Such classification is more convenient for studies dedicated to the physical mechanisms and for comparison studies of the rotation of celestial bodies.

Forced librations

The number of periodic terms, Poisson terms, the residual values and the minimum amplitude extracted from our iterative process are listed in Table 1. According to the number of terms isolated and identified and the residual values, the fits for Iσ, ρ, p 1 , and p 2 are better than the fit for τ . The main limitation in τ is due to the presence of many terms involving the mean longitudes of the planets. Some, the 9 century great inequality and 18 and 19 century terms [START_REF] Eckhardt | Planetary and earth figure perturbations in the librations of the Moon[END_REF], have such long periods that they are approximated with the polynomials over the 1070 years time-span of the ephemeris. The polynomials associated with each libration angle or polar component are listed in Table 2. During the fitting procedure a preliminary value of I=5553.60965 arcseconds was used. With the correction in the table a final value of 5553.595 arcseconds is derived for year 2000. The reference polynomials for orbital mean node and mean longitude come from [START_REF] Chapront | A new determination of lunar orbital parameters, precession constant and tidal acceleration from LLR measurements[END_REF], so any incompatibilities between those polynomials and DE421 will show up as polynomial corrections for Iσ and τ . Long period planetary terms will also be absorbed by the polynomial corrections. The constant for Iσ is associated with dissipation in the Moon [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF]). The constant terms for p 1 and p 2 arise from gravity field coefficients above second degree. Along with the τ constant, they are of interest for the transformation between the DE421 lunar principal axes of this paper and the mean Earth/mean rotation axes typically used for cartography. See [START_REF] Williams | Lunar moments, tides, orientation, and coordinate frames[END_REF]Williams et al. ( , 2008a) ) for examples of transformations.

The slowly changing planetary eccentricities and inclinations and the slow motion of their nodes and perihelia require Poisson terms. For p 1 and p 2 , the main limitation is related to the large difference of amplitude between the first (largest) term and the following libration terms. Indeed, the ratio of the amplitude of the largest over the next largest amplitude is 45 and 73 for p 1 and p 2 respectively. As a consequence, a small error in the determination of the amplitude of the first libration term and especially in its Poisson coefficient causes interference in the fit process. The largest Poisson terms for Iσ and ρ in electronic supplementary materials Tables 1 to 4 are associated with the lunar node Ω. There are multiple terms with similar periods that cannot be separated with the 1070 year span, so the Poisson terms account for the slightly different frequencies. Perturbations from the Earth's J 2 depend on the node referred to the moving equinox, the angle used in the Tables, while the moving ecliptic plane causes terms with a slightly different period and a different phase. Small effects from the planets would involve the difference between the lunar and planetary nodes. For τ , the large Poisson term for the annual term is because the eccentricity of the Earth-Moon orbit about the Sun is shrinking. The large 273 year Venus-Earth-Moon term depends on planetary perturbations that have changing amplitudes and phases, which can be represented with multiple terms for very long time spans. Finally, we compare our series with the forced physical librations computed analytically or numerically by [START_REF] Eckhardt | Theory of the libration of the moon[END_REF], [START_REF] Chapront | Complements to Moon's Lunar Libration Theory[END_REF] and revised version of [START_REF] Chapront | The lunar libration: comparisons between various models -a model fitted to LLR observations[END_REF], and [START_REF] Newhall | Estimation of the Lunar Physical Librations[END_REF].

The pole direction p 1 , p 2 In Eckhardt (1981) the lunar libration theory is solved for the main lunar problem (Earth-Sun-Moon system). So the Eckhardt solution coincides with our first group composed of only Delaunay arguments. Allowing for different constants, we obtain good agreement for the amplitude of each libration coefficient, but we find one discrepancy. Above the truncation level of the Eckhardt series, we find only one additional term D + l ′ -2F (27.10 days).

Chapront and Francou (2005) produced polar coordinates series of the Moon based on the analytical theory of [START_REF] Moons | Analytical theory of the libration of the Moon[END_REF] and improved by [START_REF] Chapront | Complements to Moon's Lunar Libration Theory[END_REF]. These series contain terms related to the main lunar problem, planetary terms, and free librations. At the same truncation level (1.2 mas), our series have fewer terms because the frequency analysis cannot separate frequencies spaced too close (roughly at 1.6 10 -5 rad.days -1 separation). Nevertheless the comparison gives very good agreement except about ten terms that are present only in our series or in Chapront and Francou series. Although we did not clearly identify the origin of these discrepancies, we suspect that these terms are related to the presence of the Poisson coefficients.

The orientation of the Moon Iσ, ρ The spectral behavior of Iσ and ρ has been studied in [START_REF] Newhall | Estimation of the Lunar Physical Librations[END_REF]. They determined the first 30 terms in the series. The agreement between the series is very good and the slight differences are due to the arguments of the test function that depend on polynomials in the present study, a larger number of Poisson coefficients, and obviously a different ephemeris.

The libration in longitude τ The fit of the libration in longitude τ is made difficult by the presence of many planetary terms. We compare with the Newhall & Williams (1997) series and we find a very good agreement at the same truncation level, as in the Iσ and ρ comparison.

Free libration in longitude

Figure 6 shows the resulting spectrum composed of forced and free librations for the lunar libration in longitude for the main terms (the complete spectrum with a log-log scale is shown in Fig. 7(c)). The three largest libration amplitudes result from the eccentric Earth-Moon orbit around the Sun (365.26 days), the anomalistic orbit period of the Moon around the Earth (27.55 days), and half the argument of perigee period (1095.17 days). The 10th largest term, 1.808 arcseconds amplitude with a 223.5 • phase, is a candidate for the free libration in longitude with a period measured at 1056.21 days. Indeed, this term is a combination of the free libration in longitude and two close forced terms identified by [START_REF] Eckhardt | Planetary and earth figure perturbations in the librations of the Moon[END_REF] at 1056.416 (-21V e + 23Ea + 2D -l) and 1056.342 days (V e -2Ea -D + 2l -F ). A difficulty in removing the contribution from these two forced terms results from the fact that the amplitudes of the two forced terms depend on the value of the free period. Thus, we untangle the free and the forced terms by fitting the combined period, amplitude and phase with a free period amplitude and phase while keeping the dependence of the forced amplitudes on the free frequency.

The resulting free libration in longitude has a period equal to 1056.13 days with an amplitude of 1.296" and a phase of 207.01 • at the JD of 2000.

τ free = 1.296" sin U (t) [START_REF] Chapront | A new determination of lunar orbital parameters, precession constant and tidal acceleration from LLR measurements[END_REF] where the angle U in degrees is

U (t) = 207.01 • + 360 • (t -2451545.0) 1056.13 (9)
Fig. 6 Spectrum of the libration in longitude. The presence of the free libration candidate is indicated by the arrow and the black box. The frequency of the candidate is 1056.2 days.

Libration in latitude and wobble

The free latitude mode of the Moon is very small with an amplitude around 0.03 arcseconds. This mode is detected thanks to both the 11 centuries long time series of DE421 and the quality of the fit. This small mode has a 24.156 yr period prograde motion when referred to the moving node using ρ and Iσ ρ free = 0.032" cos V (t) [START_REF] Dickey | Lunar Laser Ranging: A continuing legacy of the Apollo program[END_REF] Iσ free = 0.032" sin V (t)

V (t) = 160.81 • + 360 • (t -2451545.0) 8822.88 [START_REF] Eckhardt | Planetary and earth figure perturbations in the librations of the Moon[END_REF] but the motion of the rotation axis is retrograde with an 80.86 yr period as viewed from the inertial frame. The wobble mode is large enough to cause many terms in ρ and Iσ and is more easily expressed in the selenographic reference system p 1 , p 2 . In this case, the rotation axis of the Moon describes a small cone of semi-axes 8.18 arcseconds by 3.31 arcseconds with respect to the body polar axis with a period of 74.626 years. p 1wobble = -3.306" sin W (t) [START_REF] Eckhardt | Passing through resonance: The excitation and dissipation of the lunar free libration in longitude[END_REF] p 2wobble = 8.183" cos W (t) ( 14)

W (t) = 161.60 • + 360 • (t -2451545.0) 27257.27 (15)
In the Tables 7 and9 of the electronic supplementary materials, the wobble presents out-of-phase terms of ∼ 0.001" that reflect an increased noise level due to high correlation (0.7 ∼ 0.8) with the librational combination of -D + l -3Ea + 4M a or 26561.273 days. The beat period of the two frequencies is 2899 years that's about three times longer than the time span of the ephemeris.

Free Core Nutation

Table 3 lists the periods and the amplitudes of the free libration modes that were found. The DE421 integration included a fluid core with an oblate or flattened core/mantle boundary (CMB). This adds a fourth possible free mode, a retrograde precession of the pole in space, commonly called the free core nutation by analogy with the Earth's rotation. This mode would mainly affect the fluid core, but there would also be a small response of the mantle. With the DE421 CMB flattening the period in space would be about 197 years. The FCN period depends on the orbit period and CMB flattening. Expressed in days P F CN ∼ 27.32/fc [START_REF] Gastineau | Manuel de référence TRIP[END_REF] where fc is the DE421 CMB flattening of 3.8 10 -4 . With the ρ and Iσ angles from the DE421 integration, it should show up as a prograde period near 7500 days. A libration term at 7367 days was identified with a forced term predicted by [START_REF] Eckhardt | Planetary and earth figure perturbations in the librations of the Moon[END_REF]. A candidate pair of terms for the free core nutation occur at 7481 days in Iσ and 7468 days in ρ with amplitudes of 0.045 and 0.054 arcseconds, respectively. The equivalent terms in p 1 and p 2 have periods of 27.312 days with amplitudes of 16 mas. Therefore, the amplitudes are different in the two representation. A careful check shows that the 7481 days and 7468 days candidates are strongly (0.95) correlated with the 7367 days term. Indeed, the beat period between these periods is around one thousand years, which corresponds to the time span of the ephemeris. The amplitudes of the FCN candidates in Iσ and ρ are overestimated due to a strong correlation and the amplitude related to p 1 and p 2 seems more realistic. These numerical results apply to DE421, but the real values of the flattening, FCN amplitude and period are very uncertain.

Libration Spectra

Figures 7(a,b,c,d,e) present plots of physical libration and polar coordinates amplitudes vs periods corresponding to the Tables available in electronic forms. The free libration periods are resonances for the forced physical librations. For τ , forced terms at short periods are much reduced compared to their orbital counterparts, long periods have amplitudes similar to the terms in orbital longitude, and terms near the resonant period of 1056 days are magnified. Consequently, electronic Tables and Figure 7 show many terms from 1 yr to 273 yr. While there are also many terms near a month, these angular amplitudes are three orders-of-magnitude smaller than their orbital longitude counterparts. Periods shorter than a month are even more strongly suppressed. The motion of the lunar pole is described by the Iσ and ρ pair and the p

Discussion and Conclusion

In this paper, we have determined the amplitudes, phases and periods of the three free libration modes. Table 3 shows the periods and amplitudes for each free libration mode. The results of this paper are compared to two other studies. The second row gives the results of Newhall and Williams (1997) who analyzed DE403; we used unpublished material from that study to get phases and add digits. The bottom row shows the results obtained by Chapront et al. (1999) using DE403. The comparison is in very good agreement for the wobble mode. The three studies also agree for the blend of forced and free longitude terms. This paper and Newhall and Williams (1997) agree after extracting the free libration in longitude mode from the blend, but Chapront et al.

(1999) did not untangle the forced and free terms of the larger amplitude blend. We note also that we used a longer time series so our periods are expected to be more accurate than in [START_REF] Chapront | Complements to Moon's Lunar Libration Theory[END_REF], although some difference could result from different numerical values in the calculation of the free periods. The amplitude of the latitude mode is small and the phase is different for this paper using DE421 and the two older studies using DE403. The determination of the latitude mode in this paper is consistent for the sine and cosine components of the ρ and Iσ librations and we conclude that the latitude mode in DE421 has been recovered. Some possible excitation mechanisms have been explored in the past without satisfactory explanation. It has been shown that a recent meteoroid impact is an unlikely source of such excitation [START_REF] Peale | Dynamical consequences of meteorite impacts on the moon[END_REF]. [START_REF] Eckhardt | Passing through resonance: The excitation and dissipation of the lunar free libration in longitude[END_REF] proposed an excitation process related to a resonance crossing of the longitude normal mode (of 2.9 years) and a close forced frequency. During the evolution of the lunar orbit, the free and forced frequencies change slowly and can cross. However, the mechanism excites only the libration in longitude mode. [START_REF] Yoder | The free librations of a dissipative moon[END_REF] proposed an alternative mechanism, based on turbulent fluid core interaction, to excite the wobble mode. The new determination of the amplitudes of the free librations invites new investigation of their excitation mechanisms.

Fig. 1

 1 Fig.1Representation of the Euler angles (ψ ′ , θ ′ , ϕ ′ ) with respect to the ecliptic of date (OX ′ Y ′ Z ′ ).

Fig. 2

 2 Fig. 2 Temporal evolution of the libration angles over 1070 years.

Fig. 3

 3 Fig. 3 Ecliptic pole precesional cone over 1070 years. Units are arcseconds.

Fig. 4

 4 Fig. 4 Residuals of the librational angles over 1070 years.

Fig. 5

 5 Fig. 5 Residuals for the two polar components p 1 , p 2 over 1070 years.

  The[START_REF] Newhall | Estimation of the Lunar Physical Librations[END_REF] study also shows consistency of components for DE403, and that study agrees with[START_REF] Chapront | Complements to Moon's Lunar Libration Theory[END_REF]. The accurate determination of the physical libration initial conditions during the Lunar Laser Ranging data fit is a separate question. The integration model for DE403 had a solid Moon with tidal deformation while the DE421 model included an oblate fluid core. A fluid core will cause a free mode and small changes in the slow motion of the lunar pole in space (see Sidereal Terms in[START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF]) that may change the fit of the 81 year free latitude mode using the shorter span of Lunar Laser Ranging data (38 years for DE421, 25 years for DE403). While the detection of the latitude mode in this paper is valid for the DE421 physical librations, the amplitude and phase may be sensitive to the lunar interior model used to fit the Lunar Laser Ranging data so that different physical libration integrations have different values. Caution is also advised for the candidate free core nutation (FCN) terms. While the amplitude of about 16 milliarcseconds is substantial and those terms occur in DE421, the real value of the core/mantle boundary flattening, and consequently the FCN frequency, is still very uncertain[START_REF] Williams | Lunar Core and Mantle. What Does LLR See?[END_REF]. Analysis of another ephemeris with a different flattening would test the candidate terms.Each normal mode is characterized by a damping time scale. Based on the DE421 tide and core dissipation values given inWilliams et al. (2008b) and the expressions in[START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF], the damping times are 2 10 4 yr for the longitude mode, 1.5 10 5 yr for the 81 yr latitude mode, and 2 10 6 yr for the wobble mode. As a consequence, the observational detection of free librations requires recent excitation mechanisms compared to the damping times.

Fig. 7

 7 Fig. 7 Figures (a,b,c,d,e) show the spectrum of librations for different angles, Iσ, ρ, τ and polar orientation p 1 and p 2 . The corresponding Tables are available in the electronic supplementary materials.

Table 1

 1 Statistical values in the fit process of the online Tables in electronic supplementary materials.

		Iσ	ρ	τ	p 1	p 2
	Number of frequencies	143	134	89	143	130
	Number of Poisson terms	28	33	38	22	27
	rms (mas)	9.7	7.4	35.6	4.9	4.1
	Minimum amplitude (mas)	1.7	1.9	7.2	1.2	1.2

Table 2

 2 Polynomial coefficients a j for the libration angles. The first coefficient is for the constant term, followed by linear, quadratic and finally cubic terms.

		constant	linear	quadratic	cubic
		(as)	(as/cy)	(as/cy) 2	(as/cy) 3
	Iσ	-0.249	0.0189	-0.0036	2.45 10 -4
	ρ	-0.014	3.201 10 -4	-1.504 10 -4	1.375 10 -5
	τ	67.753	0.3162	-0.1034	7.434 10 -3
	p 1	-78.513	1.075 10 -5	-6.604 10 -6	6.663 10 -7
	p 2	0.290	2.411 10 -4	-6.364 10 -5	4.586 10 -6

Table 3

 3 Periods (days), amplitudes and phases (degrees) of the free librations determined from the ephemeris DE421 and comparison with the values obtained in the Newhall and Williams (1997) and[START_REF] Chapront | Complements to Moon's Lunar Libration Theory[END_REF] analyses.

  1 and p 2 pair. The resonances correspond to the two free modes at 8822.88 d prograde for Iσ and ρ and the wobble mode at 27257.27 d for p 1 and p 2 . When the wobble mode is converted to Iσ and ρ it gives multiple periods including 27.185 d and 27.239 d.

Acknowledgments

We thank Dale H. Boggs for integrating DE421 for a millenium. The research described in this paper was carried out at the Jet Propulsion Laboratory of the California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Government sponsorship acknowledged. Nicolas Rambaux started this research under the Oak Ridge NASA Postdoctoral Program and finished it as an Independent Research Advisor. We thank S. Bouquillon and G. Francou for providing the series of Chapront et al. (1999) and J. Laskar and M. Gastineau for providing the frequency analysis software TRIP (Gastineau & Laskar (2008)).