Ceres internal structure from geophysical constraints

Thermal evolution modeling has yielded a variety of interior structures for Ceres, ranging from a modestly differentiated interior to more advanced evolution with a dry silicate core, a hydrated silicate mantle, and a volatile‐rich crust. Here we compute the mass and hydrostatic flattening from more than one hundred billion three‐layer density models for Ceres and describe the characteristics of the population of density structures that are consistent with the Dawn observations. We show that the mass and hydrostatic flattening constraints from Ceres indicate the presence of a high‐density core with greater than a 1σ probability, but provide little constraint on the density, allowing for core compositions that range from hydrous and/or anhydrous silicates to a mixture of metal and silicates. The crustal densities are consistent with surface observations of salts, water ice, carbonates, and ammoniated clays, which indicate hydrothermal alteration, partial fractionation, and the possible settling of heavy sulfide and metallic particles, which provide a potential process for increasing mass with depth.


INTRODUCTION
Ceres, the largest object in the asteroid belt at approximately 470 km in radius, is the most accessible volatile-rich proto-planetary body (McCord et al. 2011;Russell et al. 2016). Occupying a size intermediate between asteroids and planets, it is unclear to what extent differentiation processes, including the possible formation of a dense core, played a role in the evolution of Ceres. Ceres has a bulk density of 2162 AE 8 kg m À3 , which translates to a bulk composition of~17-27 wt% water ice (McCord et al. 2011) depending on the amount of water trapped in minerals. Thermal evolution modeling shows that radiogenic heating would be sufficient to melt the ice, fractionating the ice and rock components, and leading to the production of salts and hydrated silicates (McCord and Sotin 2005; see also McSween et al. 2018). Assuming a carbonaceous chondrite composition, thermal modeling prior to the arrival of Dawn predicted that Ceres would undergo some degree of thermal metamorphism at depth (McCord and Sotin 2005;Castillo-Rogez and McCord 2010;McCord and Castillo-Rogez 2018), accompanied by the migration of salt species toward the surface (Castillo-Rogez 2011). Thermal history calculations showed that the degree of heating, therefore the degree of alteration and differentiation, is dependent on both the long-and shortlived radiogenic element concentrations (see discussion in McCord and Castillo-Rogez 2018).
In spite of the absence of potential meteorite analogs (McSween et al. 2018), there is information regarding the composition of Ceres. Prior to Dawn, remote sensing with terrestrial and space-based telescopes indicated the possible presence of ammonium-bearing minerals (King et al. 1992) and brucite mixed with carbonate (Milliken and Rivkin 2009). Dawn's Visible and Infrared Spectrometer (VIR) has confirmed the presence of ammonia-bearing hydrated minerals , as well as water ice, various types of carbonates, salts, serpentines, and organic material on Ceres' surface (De Sanctis et al. 2015Ammannito et al. 2016), requiring extensive aqueous alteration (e.g., Castillo-Rogez et al. Forthcoming). Aqueous alteration results in the separation of the material into a briny liquid and lessmobile solid residue, and iron would be largely concentrated in the solid phase. Because of Ceres's large size, separation of the brine and residue by gravitational settling would be expected with the less dense brine forming the water, carbonate, salt, and organic-rich crust and more metallic phases sinking to the core.
Crater depths (Bland 2013;Bland et al. 2016) and topographic power spectra Fu et al. 2017) constrain the rheology of the outer 40-100 km shell of Ceres. This in turn places constraints on the allowable composition of the outer 40-100 km shell. Crater relaxation limits the ice content of the outer 40 km of Ceres to be no more than 35%, with a stronger (more viscous) component or components making up the bulk of the crust Fu et al. 2017). Ahuna Mons, the only volcanic construction identified on Ceres's surface, indicates that brines at depth existed <200 My ago in at least one location , raising the possibility that a briny liquid subsurface layer may still be present, if not globally at least in local regions, as also suggested by Fu et al. (2017).
Gravity and topography from the Dawn spacecraft indicate that Ceres is a partially differentiated body . The gravity and topography constraints allow for a range of density structures. For example, assuming a partially differentiated body with a CM chondrite composition, Park et al. inferred a 70-190 km thick outer shell with a density of 1680-1950 kg m À3 and an interior with a density of 2460-2900 kg m À3 . They did not explore the possibility of further refining the interior into an intermediate density "mantle" and a higher density "core" as the problem becomes underconstrained. Additional insights come from combining the information contained in the gravity and topography data as done in the Ermakov et al. (2017) admittance analysis and Fu et al. (2017) topography relaxation modeling. In the former study, the outer shell density is matched by a mixture of volatiles and denser materials including silicates and salts, consistent with both the surface observations from the Dawn VIR instrument De Sanctis et al. 2016 and the aforementioned crater ) and geophysical constraints Fu et al. 2017). The latter studies indicate a low core density (~2400 kg m À3 ) that is interpreted as evidence for a cool thermal evolution.
Whether Ceres's differentiation proceeded as far as to produce a partially dehydrated silicate core, or even small iron core, remains an open question. Analysis of Dawn's GRaND data is consistent with several weight percent iron depletion in the equatorial average concentration of elemental iron relative to cosmochemical iron abundance or an average of CI/CM chondrites (Prettyman et al. 2017). The joint measurements of iron and hydrogen concentrations imply that Ceres underwent some degree of ice-rock fractionation, if the starting composition was similar to carbonaceous chondrites (Prettyman et al. 2017).
In this work we address this question of whether Ceres may have a high-density core using the Dawn shape data and hydrostatic flattening constraints. Although  and Ermakov et al. (2017) restricted their focus to two-layer models for the radial density structure of Ceres, we consider a suite of threelayer density models varying the densities and radii of the layers (shells). After computing a large ensemble of density models that are consistent with the mass and shape observations, we compare the results with several surface compositions and discuss how these assumed compositions would impact the interior density structure.

METHOD
The triaxial shape of Ceres has been determined using Dawn spacecraft images obtained during the High Altitude Mapping Orbit phase of the Dawn encounter with Ceres. The images cover Ceres's surface with a resolution of about 140 m ). Stereophotogrammetric reconstruction of these images enabled the production of a three-dimensional shape model of Ceres (Gaskell et al. 2008;Preusker et al. 2015). The global shape is fit by an ellipsoid with principal axes a = 4.831 9 10 5 m and b = 4.810 9 10 5 m in the equatorial plane and c = 4.459 9 10 5 m along the spin axis. The uncertainty in radius at each point is AE200 m (Park et al. 2016). Ceres's shape departs from hydrostatic equilibrium by almost 2000 m between the equatorial principal axes, or 0.4%. The inferred topographic flattening, f = (a À c)/a, is assumed to be dominated by rotation; however, for a body in hydrostatic equilibrium, the equatorial principal axes are equal (i.e., a = b). In the absence of internal density structures, the degree 2 order 1 and degree 2 order 2 coefficients of the spherical harmonic expansion of the gravitational potential are zero for a rotating body in hydrostatic equilibrium. The degree 2 coefficients of the Ceres gravity model from the LAMO orbit (Konopliv et al. 2017) are given in Table 1. The ratio of the degree 2 order 2 component of the gravitational potential to the degree 2 order 0 component of the gravitational potential is 0.031. Ermakov et al. (2017) assumed that the nonhydrostatic component of the degree 2 order 0 component of the gravitational potential is of the same order as the degree 2 order 2 component of the gravitational potential (i.e., 3%). Because we found that a AE3% range yielded too few successful models for statistical analysis, we expanded the uncertainty range of J 2 to AE5% of the inferred value in the evaluation of density models. Tests using the topographic flattening, f, produced nearly identical results to those presented here. Rambaux et al. (2015) showed that first-order hydrostatic equilibrium theory underestimates the difference between the equatorial and polar radii (i.e., a À c) by as much as 1800 m for Ceres, which in turn leads to an underestimate of the extent of mass concentration. They found that a second-order expansion of Clairaut's equations is accurate to about 200 m, consistent with uncertainty associated with the present Dawn shape and gravity models . In this work, we solve Clairaut's equations to second order as described in Chambat et al. (2010).
The three-layer oblate spheroid density structure can be described by five parameters (three densities and the radii of the two density boundaries). Here for sake of familiarity, we adopt the terms "crust" for the outer shell, "mantle" for the intermediate shell, and "core" for the inner sphere. There are only two geophysical constraints used in this study, the effect of rotation on the gravitational potential, J 2 , and the total mass of the satellite. We use a Monte Carlo approach, considering a more than 1,000,000 density surface, mantle, and core densities, and for each density triplet, we allow the core-mantle boundary and crust-mantle boundary radii to vary in 2 km increments ( Table 2). The range of crust-mantle boundary radii is consistent with Ermakov et al. (2017). For the crust (outer shell), we randomly chose density values that range from 1150 to 1450 kg m À3 , slightly expanding the range of crustal density models (1200-1357 kg m À3 ) that is consistent with admittance studies Mitri et al. Forthcoming). For the mantle (middle shell), we randomly chose mantle densities ranging from 1650 to 2450 kg m À3 , and for the core we randomly chose densities ranging from 1950 to 5150 kg m À3 . The results are nearly identical to a fully grid-based approach where we considered the densities in 10 and 20 kg m À3 uniformly spaced increments spanning the ranges above (see also Table 2).
We reject density structures unless the resulting mass is within the uncertainty bounds of the observed mass (Park et al. 2016) (Table 3) and the normalized J 2 coefficient is within 5% of 0.0119 (Konopliv et al. 2017). Varying the misfit of the mass and J 2 constraints had only a minor impact on the distribution of density structures. The theory here assumes that each shell has a uniform density; therefore, the densities described here are best thought of as representing average densities within each of the three regions.

RESULTS
Using the approach described above, the 1,000,000 randomly selected density models resulted in the evaluation of more than a hundred billion density models. Of these, only 7078 or 0.000007078% were within observed mass and J 2 bounds. We plot histograms of the crustal density (Fig. 1a), mantle density (Fig. 1b), core density (Fig. 1c), core-mantle boundary radius (Fig. 1d), and crust-mantle boundary radius (Fig. 1e) in 5 kg m À3 density bins (25 kg m À3 for the core density) and 10 km radius bins. Of the five parameters, the mantle density histogram has a shape that more closely resembles a normal distribution with a mean of 2367 kg m À3 and a standard deviation of 26 kg m À3 (Fig. 1b). The crustal density histogram has no values below 1300 kg m À3 and the frequency of crustal densities increases with increasing density (Fig. 1a). The core density distribution peaks in the 0 À1.19 9 10 À2 -2 1 4.60 9 10 À9 3.63 9 10 À9 2 2 2.46 9 10 À4 À2.75 9 10 À4 a The gravitational potential field, V, as a function of latitude, h, and longitude, φ, can be written as Vðh; uÞ ¼ P lmax l¼2 P l m¼Àl P lm ðcosðhÞÞ½C lm cosðmuÞ þ S lm sinðmuÞ where the P lm ðcosðhÞÞ's are the associated Legendre polynomials (see Konopliv et al. 2017). C 20 = ÀJ 2 .  ) 1650-2450 kg m À3 Inner sphere (core) density (q c ) 1950-5150 kg m À3 Core-mantle boundary radius (r cmb ) 1 0 -434 km in 2 km steps Crust-mantle boundary radius (r crust ) 10-434 km in 2 km steps 2375-2425 kg m À3 range with a long tail (Fig. 1c). The core-mantle boundary radii histogram (Fig. 1d) is nearly flat from 10 to 50 km with frequencies of 10% in each bin, and the frequencies drop almost linearly from 50 to 150 km with a long, low-frequency tail. The crust-mantle boundary radii distribution is a maximum of 50% at 35 km (the minimum value considered) and decreases to 0 at 425 km (Fig. 1d). While the histograms in Fig. 1 indicate the frequency of occurrence of each of the five parameters, the parameters are not independent. To address this point, we create 2-D gridded distributions of density versus radius for the crust, mantle, and core. Considering the plot of density as a function of radius (Fig. 2a), for each density model, we create a function that counts the frequency of occurrence of the densityradii pairs by incrementing the value of the function the grid points corresponding to the crustal density (horizontal axis) for all radii (vertical axis) from the surface to the crust-mantle boundary. This produces a vertical line of grid points from the crust-mantle boundary to the surface at the appropriate crustal density. We repeat the same process for the mantle densities (from the core-mantle boundary to the crustmantle boundary, Fig. 2b) and for the core densities (from the center of the body to the core-mantle boundary, Fig. 2c). We then divide each value in the grid by the total number of acceptable models, giving a function between zero and one that shows the frequency of each density-radius pair. The vertical lines indicate the frequency of the density/thickness of the crust, mantle, and core, illustrating the trade-offs between the thickness and density of each layer. Considering the crust, the thickness of the crust increases as the density increases (Fig. 2a) and as was observed in the histogram plots, the higher crustal density models occur more frequently. For the mantle (Fig. 2b), once again, the mean density of 2367 kg m À3 stands out; however, we now see that there is a 20% probability that models have mantle densities in the range 2360-2380 kg m À3 over the radius range from 434 km to 80 km radius (tracing the 0.2 or 20% contour), while there is a 0.1 or 10% probability that the models have mantle densities in the range 2340-2385 kg m À3 over the radius range from 434 km to 25 km radius. For the core (Fig. 2c), the frequency plot has a number of vertical streaks (light gray), indicating that the number of successful is   insufficient to create a smooth frequency map. However, the trend is clear, the core radius decreases as the core density increases. The most frequent core densities are the range 2385-2405 kg m À3 . While this is larger than the mean mantle density (2367 kg m À3 ), indicating the likely presence of a core, it does partially overlap the 1r standard deviation of the mean mantle density (2367 + 26 = 2393 kg m À3 ). Based on admittance modeling, Ermakov et al. (2017) found the best-fitting, two-layer density structures have crustal densities between 1200 and 1357 kg m À3 depending on the isostatic mechanism assumed in their analysis, with a preferred crustal density of 1287 þ70 À87 kg m À3 and a crustal thickness of 41 þ3:2 À4:7 km. They estimate the mean interior density to be 2434 þ5 À8 kg m À3 . In another admittance study, Mitri et al. (Forthcoming) report a crustal density of 1183 þ141 À177 kg m À3 and a crustal thickness of 38 þ7 À6 km with a mean interior density of 2442 þ119 À94 kg m À3 . Both the Ermakov et al. and the Mitri et al. results overlap when taking into account the range of reported uncertainties. The three-layer models described above have an average mantle density value of 2367 kg m À3 , smaller than and outside of the uncertainty range of the interior densities reported in the Ermakov et al. study even when taking into account the 1 À r standard deviation (2341-2393 kg m À3 ), although it does overlap the more generous bounds of the Mitri et al. study. The difference in mean mantle density is not surprising because the crustal densities in the hydrostatic density models are always greater than 1300 kg m À3 , which is larger than the preferred value from the Ermakov et al. study. It is not surprising that the admittance and hydrostatic studies find slightly different density structures because the degree 2 and higher order degree admittance results do not agree (see fig. 11 in Ermakov et al. 2017). There is something unusual about Ceres at degree 2. Figure 3 shows the resulting gridded frequency maps after truncating the population of acceptable density models from this study to the range of crustal densities inferred from the Ermakov et al. study (1200-1357 kg m À3 ). In this truncated distribution, the mean mantle density rises slightly to 2379 kg m À3 and the standard deviation is reduced to 10 kg m À3 . The most frequent core densities decrease to 2385-2395 kg m À3 for a core that can be as much as 350 km in radius, with additional peaks in the density ranges 2510-2570 kg m À3 (<80 km in radius), 2600-2620 kg m À3 (<30 km in radius), 2645-2665 kg m À3 (<30 km in radius), and 2685-2705 kg m À3 (<30 km in radius).
While instructive, these figures still do not illustrate the trade-offs among the five parameters. For example, if the density of the crust increases, then in order to maintain the mass constraint, either the density of one of the other layers must decrease or the thickness of a dense shell must be reduced. The most direct way to address the possible presence of a core is to plot a histogram of the density difference between the core and mantle for each model in 50 kg m À3 bins. The results for the full population (Fig. 4a) and the populations truncated to the Ermakov et al. crustal density range (Fig. 4b) have similar trends with 2.5% of the models having equal mantle and core densities and Fig. 2. Plots of frequency of (a) crust, (b) mantle, and (c) core densities and radii (thickness) for the population of models from the Monte Carlo simulation. The vertical streaks are the result of the assumed uniform density within each shell. a peak density difference in the 50-200 kg m À3 range with long, noisy tails in the 1-2% frequency range. While the trend in histogram in Fig. 4b matches that in Fig. 4a, the distribution is not as smooth, likely because there are fewer models in the distribution.
While the frequency in each individual bin in the histogram is small, the probability that the density is within the range between A and B is the integral of the probability distribution function between A and B. As a histogram is a discrete approximation of a probability distribution function, we can approximate the integral by plotting the cumulative probability histogram (Fig. 5) for the two populations described in Fig. 4. While the curves differ in minor details, in each case, more than 67% of the models (1r) have core densities at least 300 kg m À3 greater than the corresponding mantle density.

DISCUSSION
It is not surprising that the mantle density in the three-layer models is smaller than the mantle density in the two-layer studies, because including a dense core necessitates reducing either the crust or mantle density. However, the mean density of the mantle plus core, appropriately weighted by the respective volumes of the mantle and core from this study is still slightly smaller than the mean interior density reported in the two-layer studies. The hydrostatic theory results here consistently produce higher density, thinner crustal layers than the admittance results. The difference is currently not understood and could be related to lateral variations in crustal structure that invalidate the assumption of isostasy in the admittance studies, despinning of Ceres (Mao and McKinnon 2016), or some other as yet unknown processes. As we restrict the crustal densities to the preferred range reported by Ermakov et al., we still find a lower mantle density than reported by Fig. 3. Plots of frequency of (a) crust, (b) mantle, and (c) core densities and radii (thickness) for the population of models restricted to the preferred crustal density range of Ermakov et al. (2017). The vertical streaks are the result of the assumed uniform density within each shell. Ermakov et al. This is likely because in two-layer modeling, the mass concentration with depth can only be accommodated by the position of the crust/mantle boundary and the density increase across that boundary. In three-layer modeling, the mass concentration can be accommodated by two density boundaries and two corresponding density increases, thus the first density increase in the three-layer model is almost certainly going to be smaller than in the corresponding two-layer case. Unfortunately, we are not able to place constraints on the core density beyond stating that the denser the core, the smaller it must be. A core with a composition ranging from serpentine to iron (e.g., McCord and Sotin 2005) has been proposed for Ceres based on thermal history modeling. Taking the mantle density of 2367 kg m À3 , core-mantle density differences up to 500 kg m À3 would be consistent with hydrated silicates, and this core density range includes close to 30% of the models. An anhydrous silicate core spans a core-mantle density difference range of 500-1500 kg m À3 and approximately 30% of the models have core densities from 2700 to 3700 kg m À3 . Finally, core densities greater than 3700 kg m À3 require a dense component in addition to silicates, possibly indicating an increase in Fe content and 40% of the models have core densities greater than 3700 kg m À3 .
The inferred crustal density range of 1200-1357 kg m À3 is consistent with low-density mineral phases including water, carbonates, salts, serpentine, ammonia-bearing hydrated minerals, and organic material De Sanctis et al. 2016. The crustal composition of water ice, carbonates, phyllosilicates, and salt and/or clathrate hydrate phases is most consistent with models of an ancient ocean layer that underwent progressive freezing, leading to the concentration of salts (Neveu and Desch 2015;Fu et al. 2017). It is possible that rocky particles enriched in magnetite and sulfides were concentrated at depth, per their greater densities, during the differentiation phase, when Ceres held a deep ocean. A similar model was proposed for icy satellites by Scott et al. (2002). Using data from Dawn's Gamma Ray and Neutron Detector (GRaND), Prettyman et al. (2017) determined the equatorial average concentration of elemental iron at the surface of Ceres is 16 AE 1 wt%. This is~10% lower than cosmochemical iron abundances (e.g., McDonough and Sun 1995) and 15-30% lower than CM and CI chondrite averages (Lodders and Fegley 1998), and may provide a possible mechanism and source for the highdensity component in the large core density models.

CONCLUSIONS
Solving Clairaut's equation to second order for J 2 , the gravitational potential due to rotation, we conclude that the Dawn shape and gravity models indicate the presence of a core within Ceres with a probability greater than 1r. There is little constraint on the density of the core with acceptable core densities, allowing for compositions of hydrous and/or anhydrous silicates as well as a mixture of a denser component with the silicates. Throughout the interior of Ceres, excluding approximately the innermost 100 km radius, the average density is 2367 AE 26 kg m À3 , slightly smaller than the average interior density from admittance modeling studies Mitri et al. Forthcoming). The presence of a dense core is consistent with laboratory simulations of heavy mineral (e.g., magnetite, sulfide-rich particles) separation under hydrothermal conditions and preferential concentration at depth when Ceres held an ocean. Denser material may also point to the presence of anhydrous silicates, either remnants from the originally accreted material or resulting from partial dehydration.
Acknowledgments-We thank Walter Kiefer, Anton Ermakov, and an anonymous reviewer for their constructive comments. We thank Fr ed eric Chambat for making the Matlab program based on Chambat et al. (2010) available. S. D. K. is supported by NASA award NNX15AI30G from the Dawn at Ceres Guest Investigator Program. We thank the Dawn team for the development, cruise, orbital insertion, and operations of the Dawn spacecraft at Ceres.