N Rambaux

J Castillo-Rogez

V Dehant
email: vdehant@oma.be

P Kuchynka
email: kuchynka@imcce.fr

Constraining Ceres' interior from its Rotational Motion

Keywords: rotation -Ceres -Precession-Nutation -l, o, d

Context. Ceres is the most massive body of the asteroid belt and contains about 25 wt.% (weight percent) of water. Understanding its thermal evolution and assessing its current state are major goals of the Dawn Mission. Constraints on internal structure can be inferred from various observations. Especially, detailed knowledge of the rotational motion can help constrain the mass distribution inside the body, which in turn can lead to information on its geophysical history. Aims. We investigate the signature of the interior on the rotational motion of Ceres and discuss possible future measurements performed by the spacecraft Dawn that will help to constrain Ceres' internal structure. Methods. We compute the polar motion, precession-nutation, and length-of-day variations. We estimate the amplitudes of the rigid and non-rigid response for these various motions for models of Ceres interior constrained by recent shape data and surface properties. Results. As a general result, the amplitudes of oscillations in the rotation appear to be small, and their determination from spaceborne techniques will be challenging. For example, the amplitudes of the semi-annual and annual nutations are around ∼ 364 and ∼ 140 milli-arcseconds, and they show little variation within the parametric space of interior models envisioned for Ceres. This, combined with the very long-period of the precession motion, requires very precise measurements. We also estimate the timescale for Ceres' orientation to relax to a generalized Cassini State, and we find that the tidal dissipation within that object was probably too small to drive any significant damping of its obliquity since formation. However, combining the shape and gravity observations by Dawn offers the prospect to identify departures of non-hydrostaticity at the global and regional scale, which will be instrumental in constraining Ceres' past and current thermal state. We also discuss the existence of a possible Chandler mode in the rotational motion of Ceres, whose potential excitation by endogenic and/or exogenic processes may help detect the presence of liquid reservoirs within the asteroid.

Conclusions.

Introduction

The largest asteroid, Ceres, will be rendez-vous'd by the Dawn spacecraft in 2015. During the nominal mission, the asteroid will be mapped with high-resolution imaging in the visual and infra-red wavelengths (VIR instrument), and its gravity field will be determined to the tenth degree and order [START_REF] Mccord | Ceres: Its Origin, Evolution and Structure and Dawns Potential Contribution[END_REF]. A major goal of the Dawn Mission is to constrain the internal structure of Ceres, and especially quantify the extent of differentiation, and thus of internal evolution, achieved by this dwarf planet [START_REF] Russell | Dawn Mission to Vesta and Ceres[END_REF][START_REF] Mccord | Ceres: Its Origin, Evolution and Structure and Dawns Potential Contribution[END_REF].

While the primary objective of high-resolution camera VIR is to map the surface composition and surface feature, an interesting contribution of these observations will be to characterize Ceres' rotational motion. Rotation properties have proved to bring important constraints on the interiors of planetary bodies, (e.g., [START_REF] Mathews | Modeling of nutation and precession: new nutation series for nonrigid Earth and insights into the Earths interior[END_REF][START_REF] Koot | Constraints on the coupling at the core-mantle and inner core boundaries inferred from nutation observations[END_REF][START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF][START_REF] Margot | Large Longitude Libration of Mercury Reveals a Molten Core[END_REF][START_REF] Dehant | Lander radioscience for obtaining the rotation and orientation of Mars[END_REF] because they depend on mass distribution and viscoelastic properties. Indeed, the departures from a uniform rotation and changes in the orientation of a body are responses to an external forcing such as the gravitational force of another celestial body of the solar system or of the Sun. These responses depend on the structure and composition of the interior. In particular the possible presence of liquid layers inside a body and the elastic or inelastic properties of the solid parts drive the object's response to the external forcing. In particular when there are liquid layers within a body, such as Mercury or the icy satellites, the librational response of the solid part is amplified. A similar effect is expected for the nutation when considering rapid orientation changes such as for the Earth and Mars (resonance to the Free Core Nutation). The material composing Ceres and its physical properties thus determine the rotational behavior of Ceres. The observation of the response of Ceres to any forcing may thus provide information on its interior.

The purpose of this paper is to investigate the temporal variations of Ceres' rotation in order to identify potential observations to be performed by the Dawn Mission. So far, observations of Ceres' surface have been obtained with the Hubble Space Telescope (HST), adaptive optics (Keck Telescope), and occultations. However, it has proved difficult to track surface landmarks with enough accuracy to constrain the polar orientation and rotation of the asteroid [START_REF] Thomas | Differentiation of the asteroid Ceres as revealed by its shape[END_REF][START_REF] Carry | Near-infrared mapping and physical properties of the dwarf-planet Ceres[END_REF][START_REF] Drummond | Triaxial ellipsoid dimensions and rotational poles of seven asteroids from Lick Observatory adaptive optics images, and of Ceres[END_REF]. Due to the non-spherical shape of Ceres related to the centrifugal potential, the gravitational potential of the Sun exerts a non-zero torque on Ceres dynamical figure. Consequently, Ceres' axis of figure is expected to exhibit a precessional motion around its normal to the orbital plane and the periodic part of the torque generates a nutational motion along the precessional cone. In addition, in the reference frame of Ceres, the spin axis describes a polar motion or wobble around the figure axis, and, finally, tidal deformations arising from the Sun induce perturbations in its rotational velocity, leading to length-of-day variations (l.o.d.).

The first part of this article describes the shape and interior structure parameters used for modeling Ceres. Then, we review the observed pole position and discuss briefly its secular evolution. In Section 4, we introduce the main equations used to determine the rotational parameters of the body. The geophysical and rotational models are then combined to compute the rotational motion of Ceres as a function of polar motion and precession-nutation (Section 5). Based on these results, we discuss the prospects for characterizing Ceres' rotation with the Dawn Mission, and the constraints these observations will provide on Ceres' interior (Section 6).

The shape and interior structure of Ceres

Available constraints on Ceres' interior come from ground-based and space telescope observations and best estimates on the density and mean radius are gathered and discussed in [START_REF] Mccord | Ceres: Evolution and current state[END_REF] and Castillo-Rogez and McCord (2010) (Table 1). Key information on the interior structure comes from shape data, which led to constraints on the mean moment of inertia, assuming that the object is in hydrostatic equilibrium. Evidence for Ceres' shape hydrostaticity was first suggested by [START_REF] Millis | The size, shape, density, and albedo of Ceres from its occultation of BD+8 deg 471[END_REF] based on 13 ground-based occultation observations, who concluded that the asteroid is an oblate spheroid. This configuration was confirmed by [START_REF] Thomas | Differentiation of the asteroid Ceres as revealed by its shape[END_REF] from a dataset of ∼380 images obtained with the HST over 80% of Ceres' rotation period, as well as from other observation campaigns (e.g., [START_REF] Carry | Near-infrared mapping and physical properties of the dwarf-planet Ceres[END_REF][START_REF] Drummond | Triaxial ellipsoid dimensions and rotational poles of seven asteroids from Lick Observatory adaptive optics images, and of Ceres[END_REF]. While it is not possible to rule out the possibility that Ceres bears non-hydrostatic anomalies with amplitudes of the order of the current uncertainty on the shape data (∼ 2 km), we take as a working hypothesis that the asteroid is in hydrostatic equilibrium. We will discuss the validity of this assumption in Section 6.

Several shape models have been suggested for the past decade, inferred from different measurement techniques. Although the data are globally consistent, they show some discrepancy, as summarized in [START_REF] Zolotov | On the composition and differentiation of Ceres[END_REF] and [START_REF] Castillo-Rogez | Ceres' evolution and present state constrained by shape data[END_REF]. The difference is in part due to the surface coverage enabled by the various techniques. As pointed out by Rivkin and Volquardsen (2009), longitudinal variations in the surface composition are likely to induce a bias in the interpretation of optical images obtained over a short longitudinal range. The difference between the (ac) radii difference (where (a) and (c) are the equatorial and polar radii, respectively) inferred by [START_REF] Thomas | Differentiation of the asteroid Ceres as revealed by its shape[END_REF] and [START_REF] Carry | Near-infrared mapping and physical properties of the dwarf-planet Ceres[END_REF] is significant, of the order of 8 km, i.e., beyond the error bars of 2 km estimated in both cases. The difference between the equatorial (a) and the polar (c) radii varies from 31.5 to 35.5 km. The upper bound indicates that Ceres interior does not present large contrasts in density with depth. [START_REF] Zolotov | On the composition and differentiation of Ceres[END_REF] inferred from that observation that Ceres is not chemically differentiated and that a small density gradient is due to the variation of porosity with depth, while [START_REF] Castillo-Rogez | Ceres' evolution and present state constrained by shape data[END_REF] showed that even a warm icy satellite model whose core is dominated by hydrated silicate is consistent with that upper bound. A smaller value of (ac) is the signature of an increasing concentration with depth, for example due to the presence of an inner core composed of dry silicates (ordinary chondrite-like composition) that did not evolve since accretion or result from the dehydration of hydrated silicates [START_REF] Castillo-Rogez | Ceres' evolution and present state constrained by shape data[END_REF]. In order to compute the rigid and non-rigid response of Ceres to external perturbations, we assume that the asteroid is stratified in a rocky and icy shell, after [START_REF] Castillo-Rogez | Ceres' evolution and present state constrained by shape data[END_REF]. [START_REF] Zolotov | On the composition and differentiation of Ceres[END_REF] suggested that Ceres is a porous assemblage of hydrated minerals. However, Castillo-Rogez (submitted) demonstrated that hydrated minerals dehydrate in response to the moderate temperature increase undergone by Ceres in the course of its evolution. The main characteristics of this interior model are summarized in Fig. 1, and the geophysical parameters tested in this study are presented in Table 1. 1. Main parameters ranged in the models tested for this study for the stratified model (Fig. 1). Ceres' mean radius and density are taken as 476.2 km and 2078 kg/m 3 , respectively.

The mean moment of inertia is computed from the density profile after: where V body is the volume of the body, ρ(r) is the density inside the body as function to the radius of the point r. The geophysical information contained in the global shape and degree-two gravity coefficients is a function of the secular tidal Love numbers k s [START_REF] Munk | The Rotation of the Earth: a Geophysical Discussion[END_REF], and defined in Eq. ( 23))

I = 8π 3 V body ρ(r)r 4 dr (1)
J 2 = 1 3 mk s R a 2 (2)
where m is the rotational parameter equal to Ω 2 R 3 /GM, with Ω is the angular rotation rate, G the gravitational constant, M the asteroid mass, and R the mean radius of Ceres. In order to relate the internal structure to the observables, we can used the Radau-Darwin relationship written, in its approximated form, as [START_REF] Van Hoolst | The librations, shape, and icy shell of Europa[END_REF]:

I MR 2 = 2 3          1 - 2 5 4 -k s 4 + k s          (3)
where the body is assumed to be in hydrostatic equilibrium. We can then calculate the values of the equatorial and polar moments of inertia A and C from:

C MR 2 = I MR 2 + 2 3 J 2 (4) A MR 2 = -J 2 + C MR 2 (5)
where J 2 is the degree-two gravity coefficients. For the model presented in Figure 1 we obtain a mean moment of inertia I/MR 2 equal to 0.347 and A/MR 2 = B/MR 2 = 0.3394, C/MR 2 = 0.3623. We compute the complex tidal Love number k 2 , from the integration of the equations of motion by, e.g., [START_REF] Takeuchi | Seismic surface waves[END_REF] (see [START_REF] Tobie | Tidal dissipation within large icy satellites: Applications to Europa and Titan[END_REF][START_REF] Castillo-Rogez | The tidal history of Iapetus. Dissipative spin dynamics in the light of a rened geophysical model[END_REF] for details about the computational approach.) Ceres' dissipation factor is inferred from the imaginary part of k 2 . Mechanical attenuation is computed after the composite dissipation law introduced by [START_REF] Castillo-Rogez | The tidal history of Iapetus. Dissipative spin dynamics in the light of a rened geophysical model[END_REF]. That model is based on the observation that the attenuation spectrum of planetary materials shows a major shift in frequency-dependence as a function of the Maxwell time τ M characterizing these materials (ratio of the viscosity to the shear modulus). At forcing frequencies greater than 2π/τ M , the dependence of the dissipation factor Q on the angular frequency χ is such that Q -1 ∼ χ -γ with γ = 0.2 -0.4. At low frequencies, the dissipation factor follows a Maxwellian behavior such that γ = 1. This change reflects an evolution in the microstructural mechanisms driving dissipation: anelasticity-driven at high frequency and viscosity-driven at low frequency. 1 Castillo-Rogez et al. (2011) parameterized anelasticity using the Andrade model, with application to Iapetus, an icy satellite subject to a tidal stressing of a few kPa. The tidal stress amplitude in Ceres is of the order of 100 Pa, thus a priori we expect that the response to that stress involves the same physical mechanisms as described by [START_REF] Castillo-Rogez | The tidal history of Iapetus. Dissipative spin dynamics in the light of a rened geophysical model[END_REF] for Iapetus. For the parameters displayed in Table 1, we find k 2 of the order of 10 -3 . Considering the absence of robust constraints on Ceres' temperature profile, a detailed calculation of the dissipation factor is meaningless. However, that parameter can be roughly quantified as a function of frequency. Indeed, [START_REF] Castillo-Rogez | The tidal history of Iapetus. Dissipative spin dynamics in the light of a rened geophysical model[END_REF] demonstrated that the dissipation factor of a water-rich object tends toward 1 at Ceres' orbital period (1681 days), but can significantly become greater than 100 at forcing periods as short as Ceres' spin period. [START_REF] Drummond | Triaxial ellipsoid dimensions and rotational poles of seven asteroids from Lick Observatory adaptive optics images, and of Ceres[END_REF], green shaded line, from space telescope and ground-based measurements. The black dot represents the position of the orbital pole [START_REF] Giorgini | JPL's On-Line Solar System Data Service[END_REF]. The cross represents the mean intersection of the three measurements and the ellipses the uncertainties.

Ceres' Pole Position

Polar orientation

The precession-nutation theory of Ceres is defined for a given pole position of Ceres in space. However, the determination of the pole position of Ceres is difficult because of the size of the object and the lack of outstanding spectral features at its surface. Therefore, in this first section, we use the available observations based on adaptive optics and HST [START_REF] Thomas | Differentiation of the asteroid Ceres as revealed by its shape[END_REF][START_REF] Drummond | Triaxial ellipsoid dimensions and rotational poles of seven asteroids from Lick Observatory adaptive optics images, and of Ceres[END_REF][START_REF] Carry | Near-infrared mapping and physical properties of the dwarf-planet Ceres[END_REF] to address Ceres' pole position. The orientation data for Ceres' pole are provided in the International Celestial Reference Frame (ICRF) and are listed in the Table 2. The first column represents the right ascension α s , the second column the declination δ s , and the last column the error bars. The Figure 2 displays the three spin pole determinations of Table 2 projected onto the XY plane of the ICRF, i.e. X ICRF and Y ICRF :

X ICRF = cos δ cos α, (6)
Y ICRF = cos δ sin α. (7)
The three determinations overlap and their intersection is centered around α s =289.658 deg, δ s =63.189 deg, i.e. (X ICRF = 0.1517, Y ICRF = -0.4248), indicated by a cross in the Figure 2. Orbital orientation α n (deg) δ n (deg) Horizons [START_REF] Giorgini | JPL's On-Line Solar System Data Service[END_REF] 293.39 62.85

Obliquity

The Figure 2 shows the position of the orbital pole (black point) computed from the Horizons ephemerides [START_REF] Giorgini | JPL's On-Line Solar System Data Service[END_REF]) and listed at the last line of the Table 2. The orbital pole location coincides with the mean value of the pole position from [START_REF] Drummond | Triaxial ellipsoid dimensions and rotational poles of seven asteroids from Lick Observatory adaptive optics images, and of Ceres[END_REF]. It is then interesting to investigate the information contained in the obliquity value, as previously suggested by [START_REF] Bills | Forced obliquities and moments of inertia of Ceres and Vesta[END_REF]. The obliquity ε is defined as the angle between the normal to the orbital plane and the figure axis of Ceres. If Ceres' obliquity has reached its equilibrium position as a consequence of internal dissipation, it is possible to obtain a relationship between the obliquity and the moment of inertia known as a Cassini state (e.g. Yoder 1995):

ν n sin (ε -I) = - 3 4 sin 2ε C -A C ( 8
)
where ν = -50.48 kyr and I= 10.6 deg are the precession period and inclination of the orbit of Ceres with respect to the ecliptic plane (which is taken coincident with the Laplace plane), n is the mean motion, and A, C are the moments of inertia of Ceres (here, A = B < C). This formulae comes from the generalized Cassini states that result from an equilibrium position of the spin axis by taking into account the precessing orbit of the body [START_REF] Colombo | Cassini's second and third laws[END_REF][START_REF] Peale | Generalized Cassini's Laws[END_REF][START_REF] Henrard | Colombo's top[END_REF][START_REF] Lemaitre | The 3:2 spin-orbit resonant motion of Mercury[END_REF]. By using a simple uniform precessional orbital period at 50.48 ky (Bills & Nimmo 2010) and our model of Ceres differentiated into a rocky core and icy shell (Section 2), we obtain an equilibrium obliquity of ∼0.01 deg. This small value is mainly due to the long precession period (∼-50.48 kyr) with respect to the orbital period (1681 days). The secular motion of Ceres is influenced by the oscillation at -22 kyr (see Bills & Nimmo 2010) and in this case the equilibrium obliquity is equal to 0.02 deg.

Relaxation to the Cassini state is achieved when the obliquity ε meets this equilibrium criterion.

From the right ascension α n and declination δ n of the orbit pole, we could express the obliquity as

cos ε = sin δ s sin δ n + cos δ s cos δ n cos (α s -α n ) (9)
To compute the orbit pole coordinates, we use the Horizons ephemeries that provides the orbital coordinates in the ecliptic reference frame (the orbital inclination of Ceres and the ascending node are 10.6 deg and 80.5 deg, respectively). Then we express these coordinates in the ICRF by using the Earth's obliquity. The final coordinates of the orbital pole are reported in Table 2. The obliquity is equal to 4.01 degrees for the mean pole orientation of [START_REF] Thomas | Differentiation of the asteroid Ceres as revealed by its shape[END_REF], 0.23 deg for [START_REF] Drummond | Triaxial ellipsoid dimensions and rotational poles of seven asteroids from Lick Observatory adaptive optics images, and of Ceres[END_REF], and 3.91 deg for [START_REF] Carry | Near-infrared mapping and physical properties of the dwarf-planet Ceres[END_REF]. The observation of [START_REF] Drummond | Triaxial ellipsoid dimensions and rotational poles of seven asteroids from Lick Observatory adaptive optics images, and of Ceres[END_REF] seems to be close to that expected if Ceres is relaxed to a Cassini state. The uncertainty on the obliquity is represented in Figure 3 in the case of the measurement obtained by [START_REF] Thomas | Differentiation of the asteroid Ceres as revealed by its shape[END_REF]. The obliquity is between 0 and 10 degrees i.e. that contains the Cassini state but present very large variation.

For this paper, we use as a working reference for Ceres' orientation the upper value of the obliquity of 9.6 deg. The reason for using that large value, while data overlap for an obliquity value of ∼ 3 deg, is that it will yield upper bounds on our estimates of the rotational perturbations. This will help us to assess whether or not these perturbations can be measured with spaceborne techniques. [START_REF] Bills | Forced obliquities and moments of inertia of Ceres and Vesta[END_REF] predicted that the obliquity of Ceres is ∼12 deg based on a secular orbital model of Ceres. However, available observations and the location of the orbital pole yielded by the Horizons ephemeris indicate that the obliquity is most likely between 0 and 10 deg and consistent with the conclusion of [START_REF] Thomas | Differentiation of the asteroid Ceres as revealed by its shape[END_REF] that the obliquity of Ceres is around 3 ∼ 4 deg. The discrepancy seems to reside in the initial value used in Eq. ( 20) of [START_REF] Bills | Forced obliquities and moments of inertia of Ceres and Vesta[END_REF] and it could be resolved by using the initial values output by the Horizons ephemerides.

Damped obliquity

We now evaluate to what extent Ceres has evolved toward that equilibrium state. The obliquity damping rate ε may be computed from the following equations (e.g., Néron de Surgy et al. 1996;[START_REF] Levrard | Tidal dissipation within hot Jupiters: a new appraisal[END_REF] Ceres (for an orbit without planetary perturbations):

ε = Kn CΩ sin ε[cos ε f 1 (e) Ω 2n -f 2 (e)] ( 10
) Ω = - Kn C [ f 1 (e) 1 + cos ε 2 2 Ω n -f 2 (e) cos ǫ] (11)
where

f 1 (e) = 1 + 3e 2 + 3e 4 /8 (1 -e 2 ) 9/2 , f 2 (e) = 1 + 15e 2 /2 + 45e 4 /8 + 5e 6 /16 (1 -e 2 ) 6 ( 12
)
where e is the eccentricity and the constant K is defined as

K = 3 k 2 Q GM 2 R e M ⊙ M 2 R e D 6 n ( 13
)
with G is the gravitational constant, M ⊙ the mass of the Sun, D Ceres' semi-major axis, R e its equatorial radius, and Ω the angular rotation rate. The parameters k 2 and Q correspond to the tidal Love number and dissipation factor, respectively at the orbital frequency. For the sake of simplicity, we assume for this calculation that these parameters remain constant. The inverse dependence of that equation on D 6 suggests that Ceres' large distance to the Sun is a severe limitation to any tidal evolution of its dynamical properties. In addition, in this expression, we neglect the impact of orbital perturbations.

From Equations (11), as quoted by [START_REF] Correia | Secular Evolution of a Satellite by Tidal Effect: Application to Triton[END_REF] the timescale of evolution of the rotation rate is shorter than the timescale for the obliquity evolution. So it is expected that first the rotation reaches its equilibrium that will take a damping rate of the order of Ω ∼ 6.14

-11 k 2 Q , ǫ ∼ 1.00 -14 k 2 Q ( 14
)
where time is expressed in years, and then the obliquity will reach its equilibrium state. The parameters k 2 and Q are not known but may be approached from geophysical modeling. Let us note that Bills and Nimmo (2010) considered a situation where Ceres' material is in equilibrium at the dissipation peak. However, such a situation is less than likely. The contribution of tidal dissipation to the total heat budget of the object is negligible (below 0.1%) in comparison to insolation and, to a lesser extent, long-lived radioisotopes. This precludes that heat source to drive the geophysical state of the asteroid. For the reference model chosen for this study (detailed in the previous Section), differentiated in an icy shell and rocky core, k 2 is of the order of 10 -3 and Q is of the order of 10. In that case the dissipation time is of the order of 10 13 years for Ω and 10 17 years for ǫ. Hence we infer that Ceres' obliquity has probably not fully damped over its lifetime. As a consequence, our result contrasts with the conclusion of [START_REF] Bills | Forced obliquities and moments of inertia of Ceres and Vesta[END_REF] that Ceres' obliquity could have damped in a few hundred My, even if we assume that the asteroid is in a very dissipative state.

Rotational model of Ceres

Euler-Liouville Equations

If Ceres were perfectly spherical and rigid, then its rotation would be uniform. However, the Hubble Space Telescope measurements showed that the figure of Ceres is an oblate body (within the error bars), from which we inferred an equatorial oblateness α = (C -A)/A = 0.0675 (see section 2). Thus, the Sun exerts a non-zero torque on Ceres dynamical figure, which responds through precession and nutation of its orientation axes. The Sun also raises tides that deform its surface and perturb its rotational velocity.

Thus it is convenient to describe Ceres' rotation by using the approach developed for Earth, which is oblate at first order. The rotation of the body is described through the classical Euler-Liouville equation written as (see [START_REF] Moritz | Earth Rotation: Theory and Observations[END_REF][START_REF] Dehant | Earth Rotation Variations[END_REF])

dH dt + Ω ∧ H = Γ (15)
This describes the variations of the angular momentum H disturbed by an external torque Γ. This equation is expressed in the rotating frame tied to the body through the spin velocity Ω and written in the Tisserand frame [START_REF] Munk | The Rotation of the Earth: a Geophysical Discussion[END_REF]. As shown in the section (Section 6) the wobble damping time is some ten thousands of years and we could assume that the instantaneous axis of rotation is near the polar principal axis of the body. Thus,

Ω =         m 1 m 2 1 + m 3         Ω (16)
where Ω is the mean rotation of the body and the quantities m i are small and dimensionless. The pair m 1 , m 2 describes the polar motion of Ceres, i.e. the orientation of the rotational speed in the body reference frame, while m 3 corresponds the variation in the rotational speed as shown in the following linearized expression

||Ω|| = Ω(1 + m 3 ). ( 17
)
The angular momentum is H = IΩ with I the tensor of inertia of the body expressed as

I =         A 0 0 0 A 0 0 0 C         +         c 11 c 12 c 13 c 21 c 22 c 23 c 31 c 32 c 33         (18)
where the c i j are symmetric and represents the departure from the reference ellipsoid, i.e. the deformation of the body surface. By introducing the moment of inertia Eq. ( 18) into the dynamical equations Eq. ( 15) and developing at first order in m j and c i j , the linearized dynamical are

AΩ ṁ1 + (C -A)Ω 2 m 2 + Ωċ 13 -Ω 2 c 23 = L 1 AΩ ṁ2 -(C -A)Ω 2 m 1 + Ωċ 23 + Ω 2 c 13 = L 2 (19) CΩ ṁ3 + Ωċ 33 = L 3
or by introducing complex notations, as usual for Earth rotation studies, m = m 1 + im 2 , L = L 1 + iL 2 , c = c 13 + ic 23 , we obtain one complex equation for the polar motion

A ṁ -iαAΩm + ċ + iΩc = L Ω ( 20
)
and one equation for the l.o.d (length of day) variations

CΩ ṁ3 + Ωċ 33 = L 3 (21)
It is interesting to note that, at the order of approximation of axi-symmetric body, the polar motion and l.o.d variation are described by well separated dynamical equations.

Deformed tensor of inertia

The centrifugal and tidal potentials deform the body and that deformation may be expressed as [START_REF] Dehant | Atmospheric and oceanic excitation of the rotation of a three-layer Earth[END_REF])

c = αA k 2 k s m -3αA k 2 k s W 21 Ω 2 d 2 (22)
where the first term results from the centrifugal potential and the second one from the tidal potential and especially the tesseral potential W 21 (see Section 4.4). G is the gravitational constant, d is the mean Sun-Ceres distance, k 2 is the tidal Love number, and k s the secular Love number defined as [START_REF] Munk | The Rotation of the Earth: a Geophysical Discussion[END_REF]:

k s = 3 (C -A)G Ω 2 R 5 (23)
where R the radius of the surface.

The c 33 tensor varies as a function of the centrifugal and tidal potentials after (Greff-Lefftz et al. 2000)

c 33 = - 4 3 αA k 2 k s m 3 -2αC k 2 k s W 20 Ω 2 d 2 (24)
where W 20 is the zonal potential (see Section 4.4).

Gravitational torque

The tesseral degree 2 tidal potential W acts on the Ceres' equatorial bulge and thus involves an equatorial torque such that [START_REF] Dehant | Earth Rotation Variations[END_REF])

L 1 + iL 2 = - 3iαA d 2 W 21 ( 25
)
where the complex potential W 21 is developed in the next section. The torque L 3 is equal to zero due to the symmetry axis of Ceres.

The tidal potential

The gravitational tidal potential induced by the Sun may be expressed in both, a frame tied to Ceres (MBRF=Mean Body Reference Frame) and the celestial frame (MCRF=Mean Celestial Reference Frame). Choosing one reference frame, phenomena induced by the gravitational forcing have to be expressed in the same frame, with the particularities that the frequency in the frame tied to Ceres and the frequency in space are related by the rotational velocity of the body (see Eq. ( 29)). So a constant torque applied in the MCRF will appear periodical in the MBRF at the rotational frequency and vice versa.

Following, the method of [START_REF] Dehant | Earth Rotation Variations[END_REF], the degree-two potential exerted by the Sun on Ceres in the MBRF is developed as

W 2 = W 20 + W 21 W 20 = GM d 3 (z 2 - 1 3 ) ( 26
)
and

W 21 = GM d 3 (xz + iyz) (27)
where the W 20 leads to the zonal part and W 21 leads to the sectorial part in the tidal torque. Here, the sectorial part is zero due to the axi-symmetry shape of the body. The cosine directions (x, y, z) are the direction of the Sun in the MBRF. They are evaluated from the Horizons ephemeris [START_REF] Giorgini | JPL's On-Line Solar System Data Service[END_REF] and rotate from the ecliptic frame to the MBRF by using the polar direction of Ceres with a right ascension of 286 deg and declination of 54 deg consistent with an obliquity of 9.2 deg. The uncertainty related to the direction of the polar direction of Ceres is discussed in section 3.1. The zonal part is developed as a Fourier series W 20 = W j 20 e iω j t where W j 20 contains also the phases and the tesseral part is developed as Fourier series according prograde (index +) and retrograde (index -) components:

W 21 = W + 21 e iω j t + W - 21 e iω -j t (28)
where we used the same notations as [START_REF] Roosbeek | Potentiel de marées, nutations et précession d'une Terre rigide[END_REF][START_REF] Roosbeek | RDAN97: An Analytical Development of Rigid Earth Nutation Series Using the Torque Approach[END_REF]. The prograde and retrograde circular motions allow to express the elliptical motion in two symmetric motion components. The frequencies ω j of the tidal potential are expressed in the MBRF, and they are related to the prograde and retrograde frequencies (∆ω j , -∆ω j ) expressed in the space MCRF, through the following relation:

ω j = ∆ω j -Ω ω -j = -∆ω j -Ω ( 29
)
where Ω is the rotation period close to 9 hours. Henceforth, the periods appear long-period in space and short period in the frame tied to Ceres.

In the case of the Earth, one goes from the terrestrial reference frame tied to the planet to the celestial frame using several rotation matrices. These rotations first bring the terrestrial frame attached to the figure axis of the Earth to the intermediate pole, accounting for polar motion; then a rotation is performed along the true equator of date around this intermediate pole, accounting for the Earth's rotation (uniform part and length-of-day or UT1 variations). Then the precession and nutations are accounted for in order to bring the true equator of date to the celestial frame. Precession and nutations are thus those of the true equator of date. One has to keep in mind however that the choice concerning the intermediate frame is purely conventional. The more logical choice is of course related to the way UT1 or the length-of-day variations are expressed. In the recently adopted conventions for Earth, the intermediate frame is the equator of the CIP (Celestial Intermediate Pole), a conventional pole that has no retrograde diurnal motion in a reference frame tied to the Earth, and only long-period motions (precession, nutations) in space. The instantaneous rotation pole and the mean rotation pole are not identical; they differ by small changes in their direction due to atmosphere, ocean and hydrology excitation of polar motion at very short periods. When computing the precession and nutations, these axes have identical long-period motion in space and retrograde diurnal motion in a frame tied to the Earth. Similarly, for Ceres, we can ignore the differences at short period in space. We then work with the instantaneous rotation pole (m 1 , m 2 ) in a frame tied to Ceres. Long-term motion in that frame will be related to the Chandler Wobble, if excited. Long-term motion of the pole in space or retrograde diurnal motions in the frame tied to Ceres are representing precession and nutations. We have a one-to-one relation between the frequencies of these motions in a frame tied to Ceres and in space. First, we solve the polar motion by introducing the potential expressed in (28) and using a Fourier transform e iσt to express the budget equations at a given frequency in Equation ( 20). The complex polar motion is then:

Description of the rigid rotational motion

Polar motion

m = m c e iσ c t e -λt + (30) j -3αW + 21 Ω 2 d 2 Ω -(Ω + ω j )k 2 /k s ω j -Ωα + α(Ω + ω j )k 2 /k s e i(ω j t+φ j ) + -3αW - 21 Ω 2 d 2 Ω -(Ω + ω -j )k 2 /k s ω -j -Ωα + α(Ω + ω -j )k 2 /k s e i(ω -j t+φ j )
hence the polar motion is composed of a free mode (first term) and a sum of forced modes, prograde and retrograde.

The frequency of the free mode is called the Chandler frequency σ c by analogy with Earth rotation and is written as

σ c = αΩ k 2 s -|k 2 |k s cos δ + |k 2 |α(k s cos δ -|k 2 |) k 2 s + 2α|k 2 |k s cos δ + α 2 |k 2 | 2 (31)
Its period is about 5.48 days and the correction due to the deformation k 2 = |k 2 |e (-iδ) with δ the phase lag representing the dissipative part, is between 3 and 5% of its value. This contrast to the Earth's case, for which that deformation induces a difference of 100 days in the period. The difference of behavior comes mainly from the value of α that differs by about one order of magnitude between the two bodies. The amplitude m c and the phase φ c of the Chandler mode depend on the dynamical and geophysical history of the body (see discussion section 6.2). The amplitude of the Chandler mode is damped with a typical timescale 1/λ, a function of the imaginary part of the Love number and expressed as

λ = αΩ|k 2 |k s sin δ(α + 1) k 2 s + 2α|k 2 |k s cos δ + α 2 |k 2 | 2 (32)
Its value strongly depends on the interior model as discussed in section 6. The forced terms of the rigid Ceres' polar motion are shown in Table 3. The polar motion of Ceres oscillates at short periods close to 9 hours and its motion projected onto the surface of Ceres is very small (see Figure 4). Its amplitude multiplied by the mean radius of Ceres is ∼ 0.5 millimeter, and summing all the contributions regardless of the phase yields an amplitude no greater than 1 millimeter.

Precession-Nutation of Ceres

The rotational motion of Ceres' polar axis describes in the inertial reference frame, the MCRF, a precessional nutational motion. The rotation angles and their derivatives are easily computed by using the kinematic Euler equation allowing to express the instantaneous rotation pole components in terms of the nutation angles as

θ + i ψ sin θ = Ωme iΩt (33)
where Ω accounts for the expression of the pole in space due to the rotation around the Z-axis. So after integration of Eq (33), except for the case ω j = -Ω that leads to the precessional motion, we obtain the following nutation series

∆θ + i∆ψ sin θ = -i Ω Ω + σ c m c e i(Ω+σ c )t e -λt (34)
+i j 0

-Ω ∆ω j m + j e i(∆ω j t+φ j ) + Ω ∆ω j m - j e -i(∆ω j t+φ j )

The rigid nutations of Ceres are described by the periodic components of the last equation and they are listed in Table 4. The amplitudes of the long-period nutations are positively affected because the amplitude is inversely proportional to the forcing frequency.

The main term (in absolute amplitude) is the semi-annual nutation 2λ c related to the obliquity of Ceres, and then the terms related to harmonics. We also note the presence of a term related to Jupiter's mean longitude λ J . The amplitude of the annual nutation is around 364 mas, that is a 0.84 m surface displacement, for a mean radius of 476 km. The detection of such a small displacement requires tracking of Ceres' surface with a beacon for long periods of time.

The first term in equation ( 35) represents the Chandler mode observed from space. In this case, it has a period of 9h40 minutes i.e. an increase of 36 minutes with respect to the proper rotation of the body. The precessional motion of the figure axis is represented by the oscillation at ω j = -Ω in the body reference frame that is purely imaginary m 0 = 0.1576 mas, leading to a precession time of 226981.8 years that is longer than the 218654.2 years period calculated with the classical formulae:

ψ = - 3 2 n 2 Ω C -A C cos ǫ (35)
The discrepancy between the two results (4%) is essentially due to the ephemeris timescale used in the frequency analysis of the potential. In section 4 we have introduced the rotational equations for a non-rigid body. The application of these equations to the geophysical models shows that the non-rigid contributions to the shape deformation bear a negligible effect on polar motion and nutational motion. However, the variations of the moments of inertia in response to the tidal forcing exerted on the body generates a non-zero torque along the figure axis that would perturb the uniform rotational motion in the form of length of day (l.o.d) variations. By combining Eq. ( 21) and the inertia deformation (24) we deduce the variations of the l.o.d m 3

m 3 = 2α k 2 k s 1 -4 3 α A C k 2 k s W j 20 Ω 2 d 2 (36)
The resulting oscillations of the m 3 variations are very small under 0.001 mas (Table 5), largely below the expected accuracy for space-borne observational techniques. The term at 1374.3 days is related to Jupiter with the combination 2λ cλ J .

Table 5. m 3 variations of Ceres for a tidal love number k 2 = 0.017i 1.27 10 -7 with frequencies expressed in the inertial reference frame. The amplitudes are truncated at 10 -5 mas.

Freq

Per Amp in-phase Amp out-of-phase (rad/days) (days) (mas) (mas) 0.00374 1681.9 0.00072 0.0 0.00747 840.8 0.00020 0.0 0.01121 560.6 0.00004 0.0 0.00457 1374.3 0.00001 0.0

Wobble

The rotational motion of Ceres appears to be relatively uniform because all the nutational oscillations, polar motion, and l.o.d variations show very small amplitudes. Therefore, if a sizeable departure from a quiet rotation is detected by the Dawn Mission at a period of about of 9h40, then we could assign this motion to the Wobble. Indeed, the presence of a Chandler mode is expected as soon as any perturbation, exterior or interior to the body, shifts the figure axis from its equilibrium position. However, this mode is also damped due to internal dissipation. The Chandler mode expressed in the inertial reference frame has a period of 9h40 minutes, i. e., around 36 minutes longer than the proper rotation of the body. As the nutations in the inertial frame are of long periods (harmonics of the orbital period), then any observed departure of the uniform rotation at the short period of 9h40 might be attributed to the Chandler mode. The Chandler period is sensitive to the value of the Love number k 2 as shown in Figure 5, where the period is expressed in both reference frame MBRF (5.5 days) and MCRF (9h40). We consider a wide range of possible values for k 2 covering the spectrum of possible models envisioned for Ceres. This parameter is computed at the period of the Chandler mode at 5.5 days, i.e. the period in the body reference frame. The damping timescale T λ is proportional to Q and the damping timescale can be as long as 120,000 years for dissipative models with large Q and as short as as few decades if the object is very dissipative. Thus the damping timescale could be very short, hence a non-zero Chandler mode requires a continuous physical process or a recent impulse in order to be observable today.

In the Earth's case, the Chandler wobble is mainly excited by the atmosphere and the ocean. In the case of Ceres there is no atmosphere, but Ceres is in a rich dynamical environment, the asteroid belt, thus it is exposed to a constant meteoritic flux. Such meteoritic flux may involve impacts exciting the Chandler mode for Ceres. The wobble excitation may be expressed by using [START_REF] Peale | Dynamical consequences of meteorite impacts on the moon[END_REF] expression introduced in the case of the Moon: where m e represents the polar response of the impact, N represents the maximum angular momentum potentially induced on Ceres by a collision (N = mvR with m and v corresponding to the mass and velocity of the bolide, respectively, and R Ceres' radius ). The parameter H(tt 0 ) is the Heaviside function associated with an impact at time t 0 . This expression contains two components: the angular momentum transfer and the modification of the moment of inertia due to the ejected matter and formation of a crater. The subsequent response of the pole is composed of a constant offset due to the first term in Eq. ( 37) and excitation of the Chandler mode.

m e = - ic Aα H(t -t 0 ) + (- c A + ic Aα + N Aα )H(t -t 0 )e
We use the formulation of [START_REF] Gauchez | Simulation of post-impact rotational changes through multi-dimensional parametrization[END_REF] for crater modeling and the scaling law is borrowed from [START_REF] Holsapple | The scaling of impact processes in planetary sciences[END_REF].

We search for possible impact configurations leading to the excitation of the Chandler mode with an amplitude of 10 arcseconds, i.e., a displacement at the surface of 20 meters (the amplitude observed today would be damped due to the dissipation of the Chandler mode, so such events have to be recent). This may be achieved for a cometary projectile (heliocentric) with a diameter of 2.5 km, a density of 0.6 g/cm 3 , and a velocity of 20 km/s; or by a neighbor asteroid of 4 km diameter with a density of 1.3 g/cm 3 , colliding at 5 km/s [START_REF] Farinella | Collision rates and impact velocities in the Main Asteroid Belt[END_REF].

In order to estimate the probability of such an impact on Ceres, we survey main-belt asteroids with absolute magnitude lower than 14. The population contains approximately 25000 objects (see [START_REF] Jedicke | Observational Selection Effects in Asteroid Surveys[END_REF]. Trajectories of all the considered asteroids were calculated for a 100-year time interval assuming Keplerian orbits. A fictional object evolving on the same orbit as Ceres, but with a cross-section 100,00 times greater would experience 200 collisions with other asteroids. Scaling this value to Ceres' size and a time span of 150,000 years, we obtain 0.003 impacts on Ceres during that timeframe. An asteroid diameter of 4 km corresponds approximately to an absolute magnitude of 15. According to [START_REF] Jedicke | Observational Selection Effects in Asteroid Surveys[END_REF], these objects are 2 to 3 times more abundant than the population considered here. As a consequence, the corresponding number of impacts onto Ceres amounts to approximately 0.007 per 150,000 yrs. The probability that Ceres experienced in the last 150,000 years a collision with an object greater than 4 km appears to be small (less than 1%). Such an estimate is more difficult to calculate in the case of cometary collisions due to the lack of constraints on the possible reservoir of comets.

Another consequence of collision with large objects is the alteration of the moments of inertia of Ceres that may lead to shift of its figure axis (first term in Eq. 37). The long-term consequence of that effect needs to be studied in details. It requires to properly model the respective timescales for the relaxation of the crater and of the equatorial bulge. As noted by [START_REF] Nimmo | Reorientation of icy satellites by impact basins[END_REF], both processes rely on the mechanical properties of the icy shell and thus should proceed over the same timeframe, which increases the complexity of the problem. For Ceres, the low subsurface viscosity should promote rapid crater relaxation preventing the re-orientation.

We also checked for the possible occurrence of close encounters during the Dawn mission lifetime, which could excite Ceres' spin axis by an impulse of its gravitational torque. From realistic (non keplerian) asteroid orbits, we found no close encounter between 2010 and 2020 that could modify the rotational dynamics of Ceres. The encounters are not sufficiently close or the bodies involved are not massive enough. This estimate accounts only for main-belt asteroids with absolute magnitudes lower than 14.

Another source of excitation of the Chandler wobble may be due to the presence of an equatorial sea inside Ceres. The existence of such a water reservoir has been suggested by [START_REF] Castillo-Rogez | Ceres' evolution and present state constrained by shape data[END_REF] based on the observation that Ceres' surface temperature at the equator is close to the eutectic temperature of salt impurities expected in the asteroid. Several recent astrophysical models also suggest that Ceres accreted a significant fraction of ammonia hydrates (up to 7wt.% of the ice phase, Dodson-Robinson et al. 2009), and possibly also methanol hydrates [START_REF] Mousis | Origin of volatiles in the main belt[END_REF]. The presence of these compounds would help preserve a deep liquid layer over extended periods of time, as, e.g., the ammonia hydrate peritectic temperature in Ceres' pressure conditions is ∼176 K [START_REF] Hogenboom | The Ammonia-Water System and the Chemical Differentiation of Icy Satellites[END_REF], i.e., similar to Ceres' surface temperature. The excitation process is then related to possible current circulation and loading due to the presence of the fluid reservoir. Also, the presence of a deep liquid layer may result in increased dissipation, as recently suggested by [START_REF] Tyler | Strong ocean tidal flow and heating on moons of the outer planets[END_REF] in the case of outer planet icy satellites. Although the modeling of this process is beyond the scope of this paper, circulation in closed ocean systems and its signature on the rotation is a recent topic of interest to planetary sciences (e.g., [START_REF] Tyler | Strong ocean tidal flow and heating on moons of the outer planets[END_REF][START_REF] Noir | An experimental and numerical study of librationally driven flow in planetary cores and subsurface oceans[END_REF]).

Hydrostatic State

There are multiple sources of departure from hydrostatic equilibrium at the large scale, starting with the large contrast in temperature between the equator and the poles, of at least 50 K (Fanale and Salvail 1988). [START_REF] Castillo-Rogez | Ceres' evolution and present state constrained by shape data[END_REF] suggested that Ceres' equatorial temperature may promote the preservation of a regional deep liquid layer while the polar regions would be entirely frozen. The contrast in density between water ice and liquid water saturated in brines may be up to 60% (e.g., Prieto-Ballesteros and Kargel 2005), which would increase the difference between A and C by about 5%. Another source of density anomalies are mascons (mass concentrations), for example due to topography anomalies at the silicate core, as inferred for Ganymede from Galileo measurements (e.g., [START_REF] Palguta | Constraints on the location, magnitude, and dimensions of Ganymede's mass anomalies[END_REF]. [START_REF] Schenk | The Lumpy Shape of Enceladus and Implications for the Interior[END_REF] have suggested in the case of Enceladus that an unrelaxed core is responsible for the departure of the satellite's shape from hydrostaticity, by ∼1.5 km, even if the outer shell of the satellite is likely to have relaxed. Topographic features, for example unrelaxed craters are another source of density anomalies.

The knowledge of the principal axis moment of inertia is key to estimating the departure from hydrostatic equilibrium that is generally assumed in order to interpret degree-two gravity and oblate shape data in terms of interior properties through simple relationships (Equation 2, [START_REF] Zharkov | Models, figures, and gravitational moments of the Galilean satellites of Jupiter and icy satellites of Saturn[END_REF]). A determination of the mean moment of inertia I from (A+ B+C)/3 independently from the former equation by using the rotational motion of the body (see [START_REF] Ferrari | Geophysical Parameters of the Earth-Moon System[END_REF][START_REF] Konopliv | A global solution for the Mars static and seasonal gravity, Mars orientation, Phobos and Deimos masses, and Mars ephemeris[END_REF] would enable the detection of large variations in internal structure. Unfortunately, the amplitudes of the nutation and the precessional motion of Ceres are very small. Their measurement requires the tracking of a landmark at the surface of Ceres with an accuracy better than 10 cm, and this for at least six months. Also, since it is unlikely that Ceres' obliquity is fully damped (section 3), we will not be able to rely on the assumption that Ceres is in the Generalized Cassini state as a means to determine the principal axis moments of inertia. Therefore, the comparison of gravity and shape data appears the best approach to infer the presence of non-hydrostatic anomalies (with the Chandler mode, if detected) in the case of Ceres.

The extent of global relaxation can also be assessed from the comparison of the secular Love numbers inferred independently from the shape h s and from the degree-two gravity field k s (e.g. [START_REF] Dermott | The Shape and Internal Structure of Mimas[END_REF].

The Love number k s can be inferred from the degree-two gravity field such as Eq. 2. If the object is in hydrostatic equilibrium, then the Love numbers are related by (e.g., [START_REF] Zharkov | Models, figures, and gravitational moments of the Galilean satellites of Jupiter and icy satellites of Saturn[END_REF]:

h s = k s + 1. ( 38
)
Departure from this relationship informs on the non-hydrostaticity of Ceres. The Dawn Mission is likely to yield the gravity field of Ceres to degree 10 as an outcome of the nominal mission [START_REF] Mccord | Ceres: Its Origin, Evolution and Structure and Dawns Potential Contribution[END_REF]. The ratio of the gravity data to the topography (admittance) is generally used to constrain the degree of isostatic relaxation achieved by geological features (e.g., [START_REF] Simons | Global Variations in the Geoid/Topography Admittance of Venus[END_REF]. For example, [START_REF] Nimmo | Geophysical implications of the long-wavelength topography of Rhea[END_REF] interpreted Rhea's degree-three gravity coefficient inferred by [START_REF] Iess | Gravity field and interior of Rhea from Cassini data analysis[END_REF] from the Cassini Orbiter, as the signature of unrelaxed impact craters. Line-of-sight gravity measurements are also most appropriate for detecting lateral variations in density, to be compared against the topography measurements to be inferred from high-resolution imaging.

Conclusion

We have characterized the main components of Ceres' rotation and quantified them assuming Ceres is differentiated into a rocky core and icy shell. First, our modeling predicts that Ceres' obliquity is not constrained by the dissipative history of the asteroid. However, multiple determinations of Ceres' pole agree that its obliquity lies between 0 and 4 deg. The lower bound suggests that Ceres could be relaxed to a generalized Cassini state. However, considering the very long damping timescale, such a situation is unlikely. This uncertainty will be resolved by the Dawn . In any case, important constraints can also be inferred from combining shape and gravity data. These will yield independent determinations of the secular Love number that will be used to constrain Ceres' hydrostatic state, from which the mean moment of inertia of the asteroid can be inferred.

Then, for the stratified, solid model considered in this study, we established upper bounds on the rigid and non-rigid components of the nutations, polar motion and l.o.d. These appear too small to be inferred from space measurement techniques. Then we identified that a detectable perturbation of Ceres' spin state (wobble) may be the signature of a Chandler mode. This mode would have to be excited by recent large impacts or currents in local liquid reservoirs at depth in order to yield a sizeable signature. This aspect needs to be quantified in detail as it offers the prospect to constrain Ceres' thermal state and geophysical evolution.

Fig. 1 .

 1 Fig. 1. Interior model of Ceres used in this study. The panels show, from left to right, the petrological structure, the corresponding density profile, and the viscosities assumed in the different layers.

Fig. 2 .

 2 Fig.2. Projection of the positions of the pole of rotation in the ICRF plane inferred by[START_REF] Thomas | Differentiation of the asteroid Ceres as revealed by its shape[END_REF], red line;[START_REF] Carry | Near-infrared mapping and physical properties of the dwarf-planet Ceres[END_REF], blue dotted line;[START_REF] Drummond | Triaxial ellipsoid dimensions and rotational poles of seven asteroids from Lick Observatory adaptive optics images, and of Ceres[END_REF], green shaded line, from space telescope and ground-based measurements. The black dot represents the position of the orbital pole[START_REF] Giorgini | JPL's On-Line Solar System Data Service[END_REF]. The cross represents the mean intersection of the three measurements and the ellipses the uncertainties.

Fig. 3 .

 3 Fig. 3. Obliquity of Ceres as a function of the orientation of its spin orientation in the equatorial reference frame. The nominal value detected by Thomas et al. (2005) is α s = 291 ± 5 deg and β s = 59 ± 5 deg. The curves represent isocontours every 2 degrees.

Fig. 4 .

 4 Fig. 4. The polar motion of Ceres is circular and present a main term oscillating at 840 days.

Table 2 .

 2 Polar orientation of Ceres in right ascension (α s ), declination (δ s ) and uncertainties (∆).

	References	α s (deg) δ s (deg) ∆ (deg)
	Thomas et al. 2005	291	59	5
	Drummond and Christou 2008	293	63	4
	Carry et al. 2008	288	66	5

Table 3 .

 3 Ceres rigid Polar motion (k 2 = 0). The prograde motions are shown on the left columns of theTable whereas the right part of the Table contains the retrograde motions.

	Freq ω j	Per	Amp	Freq ω -j	Per	Amp
	(rad/days)	(days)	(mas) (rad/days)	(days)	(mas)
	-16.617 -0.37813 0.1683
	-16.613 -0.37821 0.0299	-16.620 -0.3780 -0.0101
	-16.609 -0.37830 0.1638	-16.624 -0.3780	0.0011
	-16.605 -0.37838 0.0449	-16.628 -0.3779 -0.0003
	-16.605 -0.37840 0.0009	-16.629 -0.3779	0.0000
	-16.602 -0.37847 0.0085	-16.632 -0.3778 -0.0001
	-16.598 -0.37855 0.0014	-16.635 -0.3777	0.0000

Table 4 .

 4 Ceres rigid Nutations (k 2 = 0) and corresponding argument with λ c and λ J mean longitudes of Ceres and jupiter. The amplitudes are truncated at 10 -4 mas.

	Arg	Freq	Per Amp Nut prog Amp Nut retro
		(rad/days)	(days)	(mas)	(mas)
	2λ c	0.00747	840.8	-364.296	2.504
	λ c	0.00374 1681.7	-133.129	-44.976
	3λ c	0.01121	560.6	-66.627	-0.400
	4λ c	0.01495	420.4	-9.449	-0.054
	4λ c -λ J	0.01204	521.7	-1.256	0.009
	5λ c	0.01868	336.3	-1.217	0.007
	6. Geophysical constraints from space observations
	6.1. Non-rigid contributions: l.o.d

  Period of the Chandler model as a function of the value of the Love number k 2 . The lower bound of k 2 corresponds to a fully frozen model while the upper bound is expected if Ceres contains a global ocean.

		5.72	-9.639
		5.71	-9.638
	P c MCRF(days)	5.67 5.68 5.69 5.7	-9.637 -9.636 -9.635 -9.634	P c MBRF (hours)
		5.66	-9.633
		5.65	-9.632
		0.008 0.01 0.012 0.014 0.016 0.018 0.02
		abs(k 2 )
	Fig. 5.
		iαΩ(t-t0)		(37)

Anelastic strain is recoverable, but is a source of internal friction as it involves the motion of lattice defects. Viscoelastic involves the same defects, but is not recoverable.The anelasticity of planetary materials has been much studied, and a review can be found inMcCarthy and Castillo- Rogez (2011).

Acknowledgements

The authors wish to thank Jim Williams (JPL) and Richard Gross (JPL) for valuable discussions on the secular orbits of asteroids and excitation of Chandler wobble. The authors are also thankful to the anonymous reviewers who helped improve the quality of this manuscript. Part of this work was conducted at the Jet Propulsion Laboratory, California Institute of Technology, under contract to NASA. Government sponsorship acknowledged. All rights reserved. This research benefits of a financial support from Paris Observatory (2010).


Anthony Millot
email: millot@atcom.fr

Rym Feliachi

Rodolphe Weber
email: rodolphe.weber@univ-orleans.fr

Christophe Léger
email: christophe.leger@univ-orleans.fr

C Leger

Wireless water quality monitoring in the ism frequency band: a software and cognitive radio approach

INTRODUCTION

In recent years, water monitoring has become a great concern for all environmental decision-makers. However, classical o-line water monitoring (i.e. manual sampling or data loggers requiring manual downloading) are no longer suitable for current applications since it would be cost-prohibitive to collect data. Besides, new environmental regulations have dened more stringent requirements in term of eciency, quality and costs. Among the monitoring alternatives, Wireless Sensor Networks (WSN) are very promising approaches [START_REF] Yick | Wireless sensor network survey[END_REF]. For example, in [START_REF] Capella | A wireless sensor network approach for distributed inline chemical analysis of water[END_REF], a fresh water monitoring system is described. The WSN is implemented through proprietary wireless networks, and secondary data transmission services (GSM or GPRS) provided by a mobile telephone operator. It is also possible to use free access bands. This is the case in [START_REF] Riquelme | Wireless sensor networks for precision horticulture in Southern Spain[END_REF] where a WSN based on IEEE 802.15.4 radio components at 2.4 GHz is used for precision agriculture or in [START_REF] Willis | Radio propagation model for long-range wireless sensor networks[END_REF] where a Berkeley Mica 2 Mote with a specic transceiver for low frequency transmissions (40 MHz) has been experimented. In both cases, the WSNs were designed to be used in open spaces with little human activity. In particular, no or few radio frequency interferences (RFI) occur during transmissions.

In this paper, the concept of a wireless water quality monitoring system in a highly corruptive radio environment is described. It is based on a low-rate one-way data stream over the 433 MHz ISM band. This free access band is widely used for remote control and the amount of RFI is very high. Moreover, power transmission is limited to 10 mW in this band. The challenge is to achieve long range transmission for our low data rate sensor array despite the presence of numerous RFI.

To achieve low-cost and autonomous power objectives, the system complexity has been shifted to the reception level where a multi-antenna reception system will allow spatial signal processing. Since the transmission will be based on periodic (daily) short bursts of data, an o-line processing receiver is considered. In particular, demodulation and WSN data analysis will be done o-line.

In gure 1 the principles of our WSN are described.

Detailed information is also available in [START_REF] Millot | Transmission de données environnementales en bande ISM. Application à la surveillance à distance des masses d'eau[END_REF]. In quiet radio environments, our prototype achieves a sensitivity which allows a long transmission range ( 10 km).

Thus, in this paper we focus on the RFI issue. In We have developed a dedicated software which al- Table 1 provides some parameters extracted from some of these RFI survey campaigns. The rst conclusion that can be drawn is that a better RFI environment can be found when transmitting during the night and aside from the ISM center frequency (i.e.

f 0 =433.92 MHz). The main reason is that most remote devices transmit at 433.92 MHz and that human activity is obviously higher during the day. In the low occupancy 10 kHz band, we can consider that there is at most 1 RFI at a given instant. Another point is that 98% of the RFIs are shorter than 1 s with a bandwidth less than a few kHz.

CYCLOSTATIONARY DETECTION IN THE PRESENCE OF INTERFERENCE

Our relay station frames are short, typically less than a few seconds and there are only a few transmissions per day, typically one during the night. To enhance the transmission performance in the presence of RFI, spatial ltering will be considered in a forthcoming evolution of the prototype. However, to limit both the hardware and the rmware costs, the objective is to accomplish it oline. So, by considering the central station as a phased array, the rst step is to trigger the acquisition of the array output vector at the right moment. Obviously, this can be done by scheduling the transmissions but we would also like to consider the case of emergency T and F correspond respectively to the RFI average duration and bandwidth.

transmissions.

In this section, we propose an automatic triggering system based on the cyclic correlation matrices of the phased array. In subsection 3.1, we will dene our array model. Then, in subsection 3.2, the cyclostationary approach will be described and nally in subsection 3.3, performances will be analyzed through simulations. In particular, a comparison with a classical power approach will be proposed.

Model description

Consider a phased array consisting of M antennas, each having a received signal z k (t), k = 1, • • • M . From the previous discussions, it is clear that the narrow-band condition holds and that geometric delays for each antenna and each impinging source can be represented by phase shifts. In this case, the central station output vector can be modeled in complex base-band form as: Box n°1 is a narrow band RFI, the others are impulsive RFIs.

z(t) = a s s(t)+A r r(t) + n(t) (1)
where:

z(t) = [z 1 (t)...z M (t)]
r(t) = [r 1 (t)...r K (t)] T is the K × 1 vector of the K interferers and A r = [a r1 , a r2 , ..., a r K ] is a M × K matrix where each a r k = [a 1,r k ...a M,r k ]
T is the spatial signature of the corresponding k th RFI. n(t) is the system noise vector with independent Gaussian entries. There is no specic assumption made concerning the whitening of the covariance matrix and calibration of the data (i.e. n(t)n H (t) ∞ = σ 2 I, where . ∞ is the average operator and I is the M xM identity matrix).

At this stage, the phased array is not calibrated. Thus, all spatial signatures are modeled as random phase vectors.

Cyclostationary approach

Most telecommunication signals present a hidden periodicity due to the periodic characteristics involved in the signals construction (carrier frequency, baud rate, coding scheme...). These parameters are usually scrambled and hidden by the randomness of the message to be transmitted. However, by using a cyclostationary approach, this hidden periodicity can be recovered thus making identication of the telecommunications signal possible. An exhaustive overview of cyclostationarity theory and applications can found in [START_REF] Serpedin | Bibliography on cyclostationarity[END_REF][START_REF] Gardner | Cyclostationarity: Half a century of research[END_REF]. For array signal processing, cyclostationarity was introduced by Gardner [START_REF] Gardner | Simplication of MUSIC and ESPRIT by exploitation of cyclostationarity[END_REF]. In his analysis, the correlation matrix,

R z (τ ) = z(t + τ 2 )z H (t -τ 2 )
∞ , is replaced either by the cyclic correlation matrix:

R α z (τ ) = z(t + τ 2 )z H (t - τ 2 ) exp(-j2παt) ∞ (2)
or the cyclic conjugated correlation matrix:

R α z (τ ) = z(t + τ 2 )z T (t - τ 2 ) exp(-j2παt) ∞ ( 3
)
where α is the cyclic frequency. This parameter is linked to the above-mentioned periodic characteristics.

Any signal with the right non-zero cyclic frequency will generate non-zero cyclic or cyclic conjugated correlation matrices. Inversely, any stationary signal or cyclostationary signal with dierent cyclic frequencies will generate zero in the above expressions. We will consider next that we do not receive cyclostationary interferences with the same cyclic frequencies as our WSN. Thus, the cyclic correlation matrix for our model dened by Eq. 1 becomes:

R α0 z (τ ) = a s a H s R α0 s (τ ) + A r R α0 r (τ )A H r →0 + R α0 n (τ ) →0 (4)
where α 0 is the WSN cyclic frequency and R α0 s (τ )is the mono-dimensional version of Eq. 2. Similarly, the cyclic conjugated correlation matrix becomes:

R α0 z (τ ) = a r a T r R α0 s (τ ) + A r R α0 r (τ )A T r →0 + R α0 n (τ ) →0 (5)
where α 0 is the WSN cyclic conjugated frequency and R α0 s (τ )is the mono-dimensional version of Eq. 3.

In our application, the transmission between relay stations and the central station is based on a BPSK modulation scheme with a baud rate, 1 T0 , and a carrier frequency, f 0 :

s(t) = exp(-j2πf 0 t + jϕ 0 ). k Z c k h(t -kT 0 -t 0 )
where c k is the binary random sequence

( |c k | 2 ∞ = σ 2 c
), t 0 and ϕ 0 are unknown synchronization parameters and h(t) is the pulse shaping function.

From [START_REF] Gardner | Exploitation of spectral redundancy in cyclostationary signals[END_REF], the expression of R α0 s and R α0 s can be de- rived:

R α0 s (τ ) = σ 2 c T0 exp(-j2πα 0 t 0 ). exp(j2πf 0 τ ).r α0 h (τ ) R α0 s (τ ) = σ 2 c T0 . exp(-j2π(α 0 -2f 0 )t 0 ). exp(j2 ϕ 0 ).r (α0-2f0) h (τ ) where α 0 = k T0 k Z , α 0 = 2f 0 + k T0 k Z and r α h (τ )
is given by:

r α h (τ ) = +∞ -∞ h(t + τ 2 )h(t - τ 2 ) exp(-j2παt)dt
Except for a rectangular shaping function, ∃k = 0, r k T 0 h (τ = 0) = 0. Thus, we will only consider next the case τ = 0. T λ (z) computes the maximum singular value of

R α0 z or R α0 z T F (z) computes the Frobenius norm of R α0 z or R α0 z .
The Frobenius norm is computationally less costly

(O(M 2
)) than the singular value decomposition

(O(M 3
)). T diag (z) computes the Frobenius norm of the diagonal of R α0 z or R α0 z . This is our very low computation complexity approach (O(M )).

In the following section, the performances of the different cyclic detection algorithms will be studied, and compared to classic detectors based on R z . In particular, the dierent correlation matrices will be estimated over a nite set of samples (( . ∞ ⇒ . N = 1 N N (.)).

Performance analysis

For the simulations, the following signal model is considered:

H 0 : z(t) = √ ρa r r(t) + √ 1 -ρn(t) H 1 : z(t) = σ s a s s(t) + √ ρa r r(t) + √ 1 -ρn(t)
The spatial signatures a s and a r are chosen randomly ( a s = a r = 1 and M = 8), s(t) is BPSK signal with f 0 = 0.1, T 0 = 8, h(t) is a root-raised cosine lter (50 coecients, roll-o factor=0.5). The RFI, r(t), is also BPSK modulated but with f 1 = 0.318, T 0 = 13. The system noise, n(t), is assumed to be independent Gaussian entries. The interference to noise ratio is dened by IN R = ρ 1-ρ with 0 ≤ ρ < 1. The signal to noise plus interference ratio is dened by SN R = σ 2 s .

Monte Carlo simulations are performed on dierent reception scenarios. The Fisher criterion is used to evaluate the performances of the detectors:

F isher = (E H1 [T k (z)] -E H0 [T k (z)]) 2 var H1 (T k (z)) + var H0 (T k (z)) k = λ, F, diag
where E H0,1 [.] and var H0.1 (.) are respectively the mean and the variance under the H 0 or H 1 hypothesis.

section 2 ,

 2 analyses of some recent RFI surveys in the 433 MHZ band are presented. Based on that knowledge, we propose to implement a spatial cyclostationary detector to trigger the automatic storage of the WSN waveforms. Section 3 describes this cyclostationary approach. In particular, the performances of several implementations are compared through simulations. In section 4, we conclude with some perspectives. 2. RFI ANALYSIS IN THE 433 MHZ ISM BAND An experimental set-up has been developed for RFI surveys. It is based on a digital board from Xilinx (Xtreme DSP IV) which acquires the signal from a dedicated RF front-end centered in the 433 MHz ISM band. A digital down-converter is implemented in this board so that 10 kHz bandwidth waveforms can be acquired and stored in a computer for o-line statistical analysis. Several RFI survey campaigns have been conducted in a semi-urban area, for dierent frequencies and schedules.

Figure 1 .

 1 Figure 1.c shows the implantation of this receiver in suburb of Orléans (France). The location of the experimental relay station is also given for information. Figure 1.b compares the level of the relay station signal with the level of surrounding RFIs.

Figure 1 :

 1 Figure 1: (a) Description of our wireless water monitoring network demonstrator. A six-parameter probe (Aqua50 from NeoTex-Ponsel) is buried in ground water. It communicates to a surface emitter through an SDI12 wire protocol. Within a radius of 1km, several ( 20) probe/emitter couples transmit their data to a relay station (pictures on the bottom right). Finally, all the relay nodes transmit probe information to a central receiver. This is a one-way star topology. The transmission is based on periodic (daily) short bursts of data (a few seconds at 1200 bds). In quiet radio environments, our prototype achieves a sensitivity which allows a long transmission range ( 10 km). The main issue is the RFI present in the 433 MHz band. (b) Spectrum of the 433 MHz ISM band. It shows the maximum power over frequencies after 80 minutes of observation. The small black box corresponds to the received power level of our 10 mW relay station signal located 10 km away. Other spectral lines correspond to RFI. (c) Experiment location in the suburb of Orléans (France).

Figure 2 :

 2 Figure 2: RFI analysis in the 433 MHz ISM band. The antenna signal is acquired with the Xilinx board and then, digitally down-converted to a 10 kHz bandwidth complex waveform (here, from 433.910 MHz to 433.920 MHz). O-line, a lter bank with 256 channels provides a time-frequency representation for automatic segmentation. In the above example, 11 RFIs were detected (black boxes with a number on the bottom right).

T

  is the M × 1 vector of central station signals at time t. a s = [a 1,s ...a M,s ] T is the spatial signature of the relay station signal, s(t).

For

  the detection of our WSN signal, we propose to use the cyclic correlation matrix R α0 z at α 0 = 1 T0 or the cyclic conjugated correlation matrix R α0 z at α 0 = 2f 0 . Equations 4 and 5 have shown that asymptotically these matrices depend on the WSN signal only. Therefore, these matrices are rank one and asymptotically, all their singular values are zeros except one. This maximum singular value is directly linked with R α0 s and R α0 s . We propose to test 3 detectors based on R α0 z or R α0 z :

Figure 4

 4 draws the performances for dierent SN R and dierent percentages of noise uctuations between antennas. Once again, the robustness of cyclostationary approaches can be observed.3. scenario 3: a realistic simulation with 10% noise uctuations and one RFI with SN R = 8dB . Figure5draws the performances for dierent IN R and different numbers of samples. Better performances are obtained with α 0 = 2f 0 than with α 0 = 1 T0 . The reason is that r (α0-2f0) h > r α0 h . However, in practice it could be dicult to have exact knowledge of f 0 before any carrier synchronization. Note that T diag (z) detectors are less ecient than T F (z) detectors. Indeed, less cyclic information is used in the second case. In fact, T diag (z) can be seen as the average of M mono-dimensional cyclostationary detectors. From these simulations and from measurements on site ( IN R = -15dB and SN R = 8dB at 40 km), to ensure correct detection (Fisher >10) of our WSN signal at 40 km, N = 2048 samples (respectively N = 8192) are necessary for T F (z) (resp. T diag (z)) with α 0 = 1 T0 . These values go down to N = 512 with α 0 = 2f 0 . 4. PERSPECTIVES In the framework of wireless environmental monitoring, a radio frequency interference (RFI) survey has been conducted in the 433 MHz ISM free-access band. From analysis of the RFI characteristics, three new detectors based on cyclostationarity have been proposed. The objective is to detect the signal of interest (SOI) with a phased array in a highly corruptive environment. These algorithms have demonstrated by simulation their robustness against the presence of interferences or uncalibrated noise. These results are based on the fact that the cyclic correlation or cyclic conjugated correlation matrices depend asymptotically only on the SOI cyclostationary property. The next step will be to use one of these detectors to trigger the acquisition of SOI frames in real time. Then, triggered data are stored for o-line spatial ltering technique (MVDR), demodulation and analysis.

Figure 3 :

 3 Figure 3: Detector performances for scenario 1. The RFI is present and the noise is calibrated (i.e. n(t)n H (t) ∞ = I). Correlation matrices are estimated with N = 8192 samples. The 2 sets of curves are associated respectively to the 2 SNR given on the right.

Figure 4 :

 4 Figure 4: Detector performances for scenario 2. No RFI but noise is uncalibrated (i.e. n(t)n H (t) ∞ = I). Noise uctuations vary from 0% to 20%. Correlation matrices are estimated with N = 8192 samples. The 2 sets of curves are associated respectively to the 2 SNR given on the right.

Figure 5 :

 5 Figure 5: Cyclostationary detector performances for scenario 3. It corresponds to a realistic simulation 10% noise uctuations and one RFI with SN R = 8 dB (this value corresponds to a 40 km transmission). For continuous lines, R α0 z with α 0 = 1 T0 was used. For dashed

Table 1 :

 1

	n°,	f 0	RFI	T	F
	when	(MHz)	per hour	(s)	(kHz)
	1, day	433.92	5277	0.3	1.6
	2, night	433.92	4549	0.36	1.58
	3, day	434.2	2864	0.28	1.31
	4, night	434.2	567	0.32	1.33

Average RFI characteristics from 4 RFI survey campaigns. f 0 is the center frequency of the 10 kHz bandwidth. For each f 0 , night and day measurements have been conducted to observe the eect of human activity. RFI per hour corresponds to the average number of RFI detected per hour.

  1. scenario 1: RFI is present and the noise is calibrated (i.e. n(t)n H (t) ∞ = I). Figure3draws the performances for dierent ρ and SN R. As expected, the sensitivity of the classic approach (based on R z ) can be clearly seen. Indeed, such power based detectors cannot distinguish between WSN and RFI signals.

	On the contrary, the cyclostationary approaches are
	very robust. For example, whatever the RFI level,
	the cyclic T λ (z) behaves as a classic T λ (z) with
	no RFI. Note that T F (z) detectors are less ecient
	than T λ (z) detectors. Indeed, cyclic correlation ma-
	trices are no longer rank one since they are estimated
	over a nite set of samples. Therefore the M -1
	weaker singular values are no longer equal to zero.
	Since the T F (z) detector is the sum of the singular
	values, uncertainty is added.
	2. scenario 2:	no RFI but uncalibrated noise (i.e.
	n(t)n H (t)

∞ = I).

Acknowledgments

The authors would like to thank the Region Centre, the French agency ANRT, the French funding committee (FRED), the Ministère de la Recherche and the European committee (ERDF) for funding part of this work. This work was undertaken in collaboration with the Pôle Capteurs Automatismes de Bourges, France.