EDP Sciences
Free Access
Issue
A&A
Volume 584, December 2015
Article Number A127
Number of page(s) 8
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201527005
Published online 04 December 2015

© ESO, 2015

1. Introduction

We revisit the hydrostatic equilibrium modeling of the global shape, gravity field, and moment of inertia for rapidly rotating bodies. This study is motivated by the arrival of the Dawn spacecraft at Ceres (Raymond et al. 2011), which is the main application of the proposed modeling approach. The global shape properties contain information on the density profile, which is critical to understanding a body’s current internal structure and past evolution (e.g., McCord & Sotin 2005; Castillo-Rogez & McCord 2010; Castillo-Rogez 2011). In addition, knowledge of the theoretical hydrostatic shape, gravity coefficients, and moments of inertia are necessary to assess non-hydrostatic contributions to the shape, which contain information on a body’s thermal current state and evolution.

The theoretical approach to modeling the hydrostatic equilibrium of a planet or satellite typically are first-order developments that are sufficient for most known bodies because these rotate slowly (i.e., tens of hours to a few days) (e.g., Dermott 1979; Zharkov et al. 1985) or simply because available observations are not accurate enough to justify more complex forward modeling. Earth is one of several exceptions, requiring a second-order development, because its parameters are more accurately known. Ceres also needs a high-order shape expansion because of its fast rotation and low density. Computing the hydrostatic shape for a multi-layered solid body is not straightforward and involves theoretical and numerical techniques to reach high-order accuracy. Here, we follow the currently best technique for Earth studies, which is based on Kopal (1960), Lanzano (1974), and Chambat et al. (2010).

In the first part of this paper (Sect. 2), we introduce the shape-modeling framework inherited from Earth geodetic studies by Chambat et al. (2010) and extended to third order by using the formalism of Lanzano (1974). Then we introduce the density profiles we used to apply our modeling to Ceres. We highlight the importance of using density profiles that preserve chemical balance, which would otherwise lead to erroneous interpretations (Sect. 3). Finally, we apply our third-order approach to three representative interior models of Ceres and compare these results with previous shape models. In addition, we compute the difference in equatorial and polar radii, a signature of internal differentiation, for a large parametric space.

2. Hydrostatic equilibrium models

The dwarf planet Ceres rotates at a period of 9.074170 ± 0.000001 h (Chamberlain et al. 2007). Under the assumption of hydrostatic equilibrium, its shape is spheroidal in response to the centrifugal and the self-gravitational potentials. Hence the hydrostatic shape is defined by the angular spin velocity and the radial density profile.

2.1. Homogeneous MacLaurin ellipsoid

It has long been known that an ellipsoid can be a hydrostatic self-gravitating equilibrium figure of a rotating homogeneous planet. For this Maclaurin’s ellipsoid, the eccentricity e (e.g., Chandrasekhar 1969) is given after (1)where Ω is the angular spin velocity, G the gravitational constant, ρ the density, and e the eccentricity defined by (2)with a and c the equatorial and polar radii, respectively.

Equation (1) is either solved numerically or expanded in series of the small parameter . In the latter case, the eccentricity is expressed as (3)up to the third order in ε. This yields the expansion of (ac): (4)where is the equivolumetric radius of the ellipsoid. We use this expansion to quantify the uncertainties in the computation of the shape parameters at each order.

2.2. Shape of a heterogeneous planet

The Maclaurin equation, Eq. (1), is valid for any eccentricity, but does not apply to stratified bodies. The shape modeling of radially heterogeneous planets originates in Clairaut’s work. It consists of solving the hydrostatic equilibrium equations to first order with respect to the height of equipotential surfaces, which is of first order with respect to the (small) geodetic parameter (M is the mass and R the external radius) (e.g., Kopal 1960; Moritz 1990; Murray & Dermott 1999). For the Earth, the method has been extended to the second order in m (Lanzano 1974; Chambat et al. 2010). For a fast rotator like Ceres, it is necessary to develop Clairaut’s approach up to third order to reach a hydrostatic shape estimate with an uncertainty of a few tens of meters consistently with the measurement accuracy expected from Dawn (Polanskey et al. 2014). We solve this problem following the equations provided in Lanzano (1974).

In this theory, the equipotential surfaces are described by their distance s to the center of the body for each co-latitude θ, and s is developed in the form of Legendre polynomials P2(cosθ),P4(cosθ),andP6(cosθ). Each surface is referenced by its mean spherical radius r in such a way that we write (5)where f2(r),f4(r),andf6(r) are functions to be determined. The condition that the surfaces must be equipotentials leads to differential relations that allow the numerical determination of these three functions for a given rotation velocity Ω and density profile ρ(r). We upgrade the free code provided by Chambat et al. (2010) with the third-order terms given by Lanzano (1974), Eqs. (26), (34), and (46). The shape parameters f2(r),f4(r),f6(r) are of magnitudes m, m2, and m3, respectively. Then the (ac) radii difference is obtained by (6)

2.3. Gravity, mass, and moment of inertia

In the same way as for the shape, the external gravity potential can be expanded in the form of Legendre polynomials up to third order by (7)where the gravitational coefficients J2n are non-dimensionalized using the length scale that we chose as the mean radius R. We here used the radius R to normalize all geodetical quantities since there is no reference value for a. When the Dawn measurements become available, these expressions can be normalized with the observed a, as is customary for Earth studies. Writing the condition that the external surface s(R,θ) is a gravity equipotential and expanding the quantities up to third order, we can express the gravitational coefficients as functions of the f2n and of m: (8)where the equatorial radius is given by .

Other useful geodetic parameters are the mass and mean inertia (one third of the trace of the inertia tensor), which are defined as (9)where V is the volume of the body. We can write it with variables r,θ instead of s,θ using the Jacobian of the transformation. By defining the angular average of a function f over the unit sphere as (10)we can express the mass and inertia for a stratified planet as (11)Using Eq. (5) and after an expansion correct up to third order, we use the properties of Legendre polynomials (e.g., Lanzano 1974, Eqs. (12), (13)) to infer

This is different from the mass and inertia of the spherical reference model defined as (14)and the differences between M and M0, as well as and I0, are of second order.

Notably, based on Eq. (5) and Pn ⟩ = 0, the mean radius r is implicitly defined in this paper as r = ⟨ s (Chambat et al. 2010). This is a practical definition for geophysics even if it does not conserve the mass. The equivolumetric radius preserves the mass but has a nonlinear definition: . Up to third order, the two are related after .

2.4. Radau-Darwin approximation

The Radau-Darwin approximation (RDA) establishes a relationship between the flattening, angular spin velocity, and the moment of inertia of a body. This approximation is valid only at first order because it is applied to solve Clairaut’s equations. It gives the following expression of the inertia (Chambat et al. 2010): (15)where k is the degree-two Love number defined by (16)which is computed from(17)We assess the inertia from the Radau-Darwin approximation to illustrate to which extend this first-order approximation is different from the third-order computation presented earlier.

3. Ceres’ differentiation state

3.1. Chemical balance

Two types of interior models have been suggested for large ice-rock bodies like Ceres: (1) variable fractions of rock and ice distributed throughout the interior; (2) differentiation of a core of hydrated silicates and an ice-rich shell. Hydrated materials can store up to 10wt.% of water. A large proportion of hydrated materials are expected in Ceres by analogy with carbonaceous chondrite parent bodies (CM) and icy satellites (see Castillo-Rogez 2015, for a review). Ground-based observations also point to the direct detection of carbonates, and possibly brucite, two products of hydrothermal activity (Milliken & Rivkin 2009). Aqueous alteration is generally accompanied by the global scale redistribution of material as a consequence of the leaching of many elements under aqueous alteration. Some of these elements (especially alkali and earth alkaline) are soluble in water and eventually freeze in the form of hydrated salts (e.g., Kargel et al. 2000). Iron-rich compounds (e.g., iron sulfide) may sink and separate in a small seed, or at least lead to densification of the hydrated core with depth (e.g., Scott et al. 2002). Hence it is important that density profiles account for the chemical balance of hydrothermal systems so that the actual extent of differentiation of the interior is properly reflected.

However, previous models (e.g., Thomas et al. 2005; McCord & Sotin 2005) have assumed models of Ceres that were stratified with a core of density ~2700 kg/m3 and an icy shell at ~1000 kg/m3. These models did not properly account for the exchange of material from the rock to the icy shell as a consequence of leaching in hydrothermal conditions and also appear unrealistically devoid of metals. This problem was studied at length by Kargel et al. (2000), for example, who explored (in their Fig. 8) multiple chemical evolution pathways for Europa. This study showed that a model of Europa highly evolved from a geochemistry standpoint is characterized by a high moment of inertia as a result of the shell enrichment in a variety of salts and hydrates. Furthermore, Kirk & Stevenson (1987) pointed out that the oceanic layer of an icy body may be enriched in phyllosilicate particles that do not sink, but remain in suspension; these would eventually become trapped in Ceres’ frozen ocean layer at the interface between the shell and the rock-dominated core.

The average density of CM material (silicates mixed with salts metals Brearley 2006) taken as analog to Ceres’ core is about 2900 kg/m3 (e.g., McKinnon 1997). Following Scott et al. (2002), we have assumed that the metals separate from the silicates, a process that may not be warranted in the context of a low-gravity body like Ceres. The density of the shell is also difficult to estimate because it depends on the evolution of the salt mixture as the ice shell freezes and is a function of the temperature profile (Kargel 1991; Kargel et al. 2000). For this study we assumed an averaged salt and phyllosilicate density of between 1000 and 2000 kg/m3 based on Kargel et al. (2000) and Kirk & Stevenson (1987).

3.2. Test cases

Table 1

Ceres dimensions.

Table 2

Geophysical parameters of Ceres.

Table 3

Preferred model of Ceres based on observations of Thomas et al. (2005) and Drummond et al. (2014) with a 5 km frozen ocean layer.

Several shape models are available for Ceres, derived from different observation techniques. Although the results are generally consistent with each other, the differences encompass a broad range of (ac) values that overlap around 33–35 km, as shown in Table 1 (Millis et al. 1987; Thomas et al. 2005; Drummond & Christou 2008; Carry et al. 2008; Drummond et al. 2014). These studies are based on three different techniques: star occultation (Millis et al. 1987), visible imaging with the Hubble Space Telescope (HST; Thomas et al. 2005), and ground-based adaptive optics (AO) for the other references. Differences between the inferred shape models are in part due to the difference in the surface coverage enabled by the various techniques. As pointed out by Rivkin & Volquardsen (2008), longitudinal variations in the surface albedo are likely to induce a bias in the interpretation of optical images obtained over a short longitudinal range. Drummond et al. (2014) also noted the possible bias in limb fitting due to local darkening.

Table 4

(ac) radii difference for different interior models of Ceres assuming the mean radius derived from the Hubble Space Telescope observations by Thomas et al. (2005) (C1, C2, C3) and the mean radius derived from adaptive optics imaging Drummond et al. (2014) (C4).

We used the measurements of Thomas et al. (2005) and Drummond et al. (2014) because those of Thomas et al. (2005) are similar to the results of Drummond & Christou (2008) and the measurements of Drummond et al. (2014) are the most recent estimate of Ceres’ shape. The corresponding geophysical parameters are listed in Table 2. We used the same rotational periods as in these papers to facilitate the comparison of the computed shapes. The results are illustrated for four interior models. Three models use the Thomas et al. (2005) mass, radius and period: (C1) a homogeneous interior with a bulk density of 2077 kg m-3, (C2) a two-layer model with a core density of 2700 kg m-3 and an icy shell density of 1000 kg m-3, and (C3) a four-layer model whose stratification reflects a more realistic petrological model, with a core composed of Fe-FeS (5500 kg m-3), serpentine assemblage (2600 kg m-3), frozen ocean (1500 kg m-3), and ice layer (1000 kg m-3). The fourth model (C4) assumes the same structure as (C3) but using the Drummond et al. (2014) mass, radius and period. The parameters describing C3 and C4 are listed in Table 3.

4. Results

In this section we present the geodetic hydrostatic equilibrium computed for Ceres using the third-order approach described above for several petrological models. First, we assess the (ac) radii using various methods. Then we present the gravity coefficients and moments of inertia for Ceres.

4.1. Shape models

Table 4 presents the (ac) radii difference computed for several interior models and modeling methods. Column C1 shows the results for the homogeneous case. This case is used to provide a benchmark of the approximation intrinsic to the other approaches.

The first section of this table corresponds to the values numerically computed by the method described in Sect. 2.1. The results agree very well (below 0.001 m) with the Maclaurin solution displayed in the second section, which corresponds to the developments to first, second, and third order of the small parameter m. Moreoever, the last row of the Maclaurin section corresponds to the exact solution equal to 39.764 km. The departure of the first-order estimate from the exact value is significant (1.8 km). It is reduced to 200 meters and 25 meters at second and third order, respectively. Hence, development to third order is required for the forward modeling to be on the same level as Dawn’s observations. The numbers in parentheses correspond to the uncertainties computed at each order by using Eq. (4). These uncertainties are computed as the difference between the exact solution and the value obtained at each order. We assume that this uncertainty also applies to the heterogeneous models.

Column C1 also contains the analytical solution of the models of (c) Dermott (1979) and (d) Zharkov et al. (1985), which are developed to first order. The comparison of these values with the first order of this paper (a) agree well, as expected. Finally, the analytical method of (e) Thomas et al. (2005) agrees well with the exact value of the Maclaurin equation. The difference in (ac) is about 150 meters, as quoted by Thomas et al. (2005).

thumbnail Fig. 1

Comparison of the solutions (ac) for the models described in Table 3 with published observations. The blue boxes represent the values of (ac) with uncertainties listed in Table 4 for the four models. For each model, the solutions are given with their uncertainties at first (left), second (center), and third order (right). Reported observations include Millis et al. (1987, Mi87), Thomas et al. (2005, Th05), Drummond & Christou (2008, Dr08), Carry et al. (2008, Ca08), and Drummond et al. (2014, Dr14).

Open with DEXTER

thumbnail Fig. 2

Differences between the equatorial and polar radii as a function of the thickness of the hydrated silicate layer for the broad parametric space described in Table 5 . Each dot corresponds to one model, and the dot separation shows the discretization used in the calculation. The blue dots correspond to models with a salt and phyllosilicate density below 1500 kg/m3, the green points to a density between 1500 kg/m3 and 2000 kg/m3, and red points to a density equal to 2000 kg/m3.

Open with DEXTER

The heterogeneous models described in Sect. 3.2 are reported in Cols. C2 through C4. The corresponding (ac) values are lower than in the homogeneous case by 7 km at third order, including the uncertainties. This is a departure of 2 km from the first-order solution (accounting for uncertainties). However, it is outside the error bars of the HST shape model, which led Thomas et al. (2005) to conclude on the stratified nature of Ceres’ interior. The first-order computations with (c) and (d) agree very well, but the comparison with the determination of Thomas et al. (2005) shows a discrepancy of about 1 km for the C2 model.

Table 5

Parametric space ranged for the models presented in Fig. 2.

Recently, Kong et al. (2010), Schubert et al. (2011), and Tricarico (2014) used an alternative approach to solve the hydrostatic shape of Ceres in terms of concentric ellipsoids. The problem with this approach is that the equilibrium shape cannot be described by heterogeneous ellipsoids as stated in Hamy’s theorem (Hamy 1889; Tassoul 1978; Moritz 1990). This equilibrium is possible only at first order in ellipticity (Clairaut’s equations); at higher orders, the equidensity surfaces can no longer be described by ellipsoids. This theorem is also called the “no go theorem” by Moritz (1990). On the other hand, the method developed by Hubbard (2012, 2013), which consists of a numerically computed closed equation, the interior equipotential surfaces, and gravitational coefficients, also yields an accurate solution of the shape, as demonstrated by this author for Jupiter.

Figure 1 illustrates that a 25-m shape measurement accuracy is required to separate the four selected interior models we highlighted here. In addition to these four examples, we built 16 343 interior models covering the range of parameters listed in Table 5. These models include an ice layer thickness of density 1000–1200 kg/m3 varying between 0 and 80 km (step 20 km); it overlies a mixture of ice, salts, and phyllosilicates ranging from 0 to 220 km. We fixed the density of this layer at 1000–2000 kg/m3, 50 kg/m3 steps to reduce the computer time. The deep interior consists of a hydrated silicate layer with a density varying between 2610 and 2990 kg/m3, which may overlie a metallic core with a density of 5500 kg/m3 of up to 80 km radius. We solved the equation of mass conservation and deduced the hydrated silicate layer thickness. Figure 2 shows the results of (ac) to third order as a function of hydrated silicate layer thickness for core radius of 0 and 80 km, with color points as function of the densities chosen for the salt-phyllosilicates layer.

Table 6

Gravity coefficients J2, J4, and J6, mean moments of inertia , moment of inertia of the spherical model , moment of inertia computed with the RDA , M the mass of the flattened hydrostatic model, and M0 mass of the spherical model.

The first-order calculation would yield an offset of about −1.5 km for all models. Figure 2 shows that, as expected, the difference (ac) is especially sensitive to the mass of material in the shell relative to the mass of the silicate mantle. Parameters of secondary influence are the thickness of an outer shell dominated by ice and the density of the hydrated silicate. On the other hand, (ac) cannot distinguish the presence of a core. Upper limit models yield an (ac) of about 35 km and cannot reproduce the mean 37.5 km radii difference reported by Drummond et al. (2014). They correspond to models without an ice layer on the top and density of phyllosilicates equal to 2000 kg/m3. Future observations of the shape combined with gravity field determinations will solve this uncertainty.

4.2. Gravity models

In addition to the equipotential shape, we computed and report in Table 6 the hydrostatic gravity coefficients of J2, J4, J6, the moments of inertia, and the masses for each interior models defined in Sect. 3.3. Table 6 shows that the difference in J2 amplitude between the homogeneous (C1) and the heterogeneous cases (C2, C3, C4) is about 30%. However, the variation of J2 across the set of heterogeneous cases (C2, C3, C4) is only ~3%. It is about 5% for J4 and less than 1% for J6. For the moment of inertia, whether or , the difference between the homogeneous (C1) and heterogeneous cases (C2, C3, C4) is about 11–13%, whereas the variation across the heterogeneous cases (C2, C3, C4) is below 3%. The values of the mean and mean spherical are similar, whereas the RDA overestimates that parameter by only 0.04%. This small offset agrees with Gao & Stevenson (2013), who found that the difference between the direct value of the moment of inertia and that computed via the RDA is small for fast-rotating bodies.

5. Discussion and conclusion

To infer the non-hydrostatic contributions to the shape and gravity, the hydrostatic shape needs to be accurately developed. We here investigated the hydrostatic shape and gravitational potential coefficients for Ceres. We expanded the Maclaurin solution to third order in small parameter m and compared test cases with the numerical integration of third-order Clairaut equations developed by Kopal (1960) and Lanzano (1974). We found that the third order is required to reach an accuracy below that expected for Dawn’s shape model (200 m/pixel) since the second order is inaccurate by about 200 meters. We also demonstrated that developments to first order underestimate the shape by up to 1.8 km. This in turn can lead to underestimating the extent of differentiation of Ceres’ interior. We also highlighted the importance of properly bookkeeping material redistribution for a body that is subject to early hydrothermal activity, as is believed to be the case for Ceres (Castillo-Rogez & McCord 2010; Castillo-Rogez 2015).

The Dawn spacecraft has been orbiting Ceres since March 2015 and will determine the shape, topography, and gravity field of Ceres with unprecedented accuracy (Raymond et al. 2011; Konopliv et al. 2011). Like for many other bodies, it is expected that the comparison of the measurements against the computed hydrostatic solution will reveal non-hydrostatic anomalies. The latter carry important information on past thermal evolution and large-scale events, such as large impact craters that may not have relaxed (e.g., Davison et al. 2015), fossil shape linked to a late stage of differentiation through silicate dehydration (Castillo-Rogez & McCord 2010), or large-scale material removal (Bowling et al. 2015). Supplementary information on the density profile can also be revealed by measuring the moment inertia by accurately measuring the precession-nutation motion from control point images and radio-science investigation as explained in Rambaux et al. (2011).

Acknowledgments

The authors are thankful to the anonymous referee for helpful comments. N. Rambaux is grateful to the CNU, Sect. 34, for supporting a six-month full-time research project through CRCT-2015 funding delivered by the MESR. Part of this work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract to NASA.

References

All Tables

Table 1

Ceres dimensions.

Table 2

Geophysical parameters of Ceres.

Table 3

Preferred model of Ceres based on observations of Thomas et al. (2005) and Drummond et al. (2014) with a 5 km frozen ocean layer.

Table 4

(ac) radii difference for different interior models of Ceres assuming the mean radius derived from the Hubble Space Telescope observations by Thomas et al. (2005) (C1, C2, C3) and the mean radius derived from adaptive optics imaging Drummond et al. (2014) (C4).

Table 5

Parametric space ranged for the models presented in Fig. 2.

Table 6

Gravity coefficients J2, J4, and J6, mean moments of inertia , moment of inertia of the spherical model , moment of inertia computed with the RDA , M the mass of the flattened hydrostatic model, and M0 mass of the spherical model.

All Figures

thumbnail Fig. 1

Comparison of the solutions (ac) for the models described in Table 3 with published observations. The blue boxes represent the values of (ac) with uncertainties listed in Table 4 for the four models. For each model, the solutions are given with their uncertainties at first (left), second (center), and third order (right). Reported observations include Millis et al. (1987, Mi87), Thomas et al. (2005, Th05), Drummond & Christou (2008, Dr08), Carry et al. (2008, Ca08), and Drummond et al. (2014, Dr14).

Open with DEXTER
In the text
thumbnail Fig. 2

Differences between the equatorial and polar radii as a function of the thickness of the hydrated silicate layer for the broad parametric space described in Table 5 . Each dot corresponds to one model, and the dot separation shows the discretization used in the calculation. The blue dots correspond to models with a salt and phyllosilicate density below 1500 kg/m3, the green points to a density between 1500 kg/m3 and 2000 kg/m3, and red points to a density equal to 2000 kg/m3.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.