Phoebe’s differentiated interior from reﬁned shape analysis (cid:63)

Context. Phoebe is an irregular satellite of Saturn, and its origin, from either between the orbits of the giant planets or the Kuiper Belt, is still uncertain. The extent of di ﬀ erentiation of its interior can potentially help inform its formation location because it is mainly determined by heat from 26-aluminum. The internal structure is reﬂected in the shape, assuming the body is relaxed to hydrostatic equilibrium. Although previous data analysis indicates Phoebe is close to hydrostatic equilibrium, its heavily cratered surface makes it di ﬃ cult to tease out its low-order shape characteristics. Aims. This paper aims to extract Phoebe’s global shape from the observations returned by the Cassini mission for comparison with uniform and stratiﬁed interior models under the assumption of hydrostatic equilibrium. Methods. The global shape is derived from ﬁtting spherical harmonics and keeping only the low-degree harmonics that represent the shape underneath the heavily cratered surface. The hydrostatic theoretical model for shape interpretation is based on the Clairaut equation developed to the third order (although the second order is su ﬃ cient in this case). Results. We show that Phoebe is di ﬀ erentiated with a mantle density between 1900 and 2400kgm − 3 . The presence of a porous surface layer further restricts the ﬁt with the observed shape. This result conﬁrms the earlier suggestion that Phoebe accreted with su ﬃ cient 26-aluminium to drive at least partial di ﬀ erentiation, favoring an origin with C-type asteroids.


Introduction
Phoebe is the largest of Saturn's irregular moons and presents unique features, such as a dark regolith surface related to hydrated silicates (Clark et al. 2005) and a retrograde orbit, which indicate that Phoebe could be a captured body. The origin of Phoebe is debated; it could be a trans-Neptunian object (TNO) (Johnson & Lunine 2005) or a C-type body based on the study of its spectrum (Clark et al. 2005) as well as on comparisons between geophysical properties expected for TNOs and C-type asteroids (Castillo-Rogez et al. 2019). Indeed, Castillo-Rogez et al. (2019) modeled the thermal evolution of TNOs assuming they formed after the extinction of short-lived radioactive elements. In the absence of a significant heat source, these bodies could retain high porosity, consistent with the low densities measured at these bodies, which are about half the densities measured at C-type asteroids (see database in Castillo-Rogez et al. 2019).
The images of Phoebe acquired by the joint NASA/ESA Cassini-Huygens mission in 2004 have revealed a heavily cratered surface. However, the overall figure of Phoebe is a flattened spheroid (Castillo-Rogez et al. 2012). Phoebe is too far away for Saturn's tidal friction to drive spin-orbit synchronization. Therefore, Phoebe has a rotation period of about 9.2735 h (Bauer et al. 2004), while its orbital period is 551 days.  Castillo-Rogez et al. (2012) analyzed the topography of Phoebe, and, by comparing these data with a theoretical model of hydrostatic equilibrium developed to the first order, they argued that Phoebe could be close to hydrostatic equilibrium within the 2 km error bars on each radius of the Cassini-derived shape radii.
We revisit the interpretation of Phoebe's shape in order to better constrain its interior structure and its possible porosity by extracting Phoebe's average form from a development in spherical harmonics. The concept is to use the spherical harmonics to filter out high-degree features (e.g., craters). In addition, we updated the modeling of Phoebe's hydrostatic shape with a development of the Clairaut equation to the third order (although the second order is sufficient considering the rotation period and the precision of the observations). As shown in Rambaux et al. (2015) for the case of rapid rotators (e.g., most asteroids and some dwarf planets), the Clairaut equation generally used to describe the shape needs to be developed to at least the second order, and potentially a third order depending on the sought precision.
The purpose of this paper is therefore to extract information on Phoebe's internal structure from the shape derived from the observations. A summary of the Cassini-Huygens observations and the fit to spherical harmonics are presented in Sect. 2. Then, the theory on hydrostatic equilibrium computation is summarized in Sect. 3. Section 4 presents the geophysical models that we explored, and the comparison between the theoretical models and the observed shape is presented in Sect. 5. The density profile inferred from this new shape inversion is consistent with A&A 643, L10 (2020) 0.9 ± 1.5 0.9 ± 1.00 -Notes. The uncertainties in Castillo-Rogez et al. (2012) and in Thomas et al. (2018) did not take into account the heavily crater surface and as argued in Castillo-Rogez et al. (2012) the estimated uncertainties is around 2 km. The quantityā corresponds to (a + b)/2. In the column Thomas et al. (2018), we invert the a and b values.
a partially differentiated interior, reinforcing the suggestion by Castillo-Rogez et al. (2012) that Phoebe went through a phase of partial melting driven by aluminum-26 decay.

Observed shape
Phoebe was observed by the Cassini-Huygens space mission during its arrival at Saturn in 2004 (Porco et al. 2005). These observations revealed a heavily cratered surface that masks its global equilibrium shape. However, Castillo-Rogez et al. (2012) inferred its global shape by using statistical analysis of craterization. The triaxial global shape from that study is given in Table 1, which lists the semi-axes a b c. Phoebe's global shape is close to an oblate body: The departure of b to a is below 1%, or 0.9 km, but it is within the error bar (see Table 1). The difference (a − c) used as a key constraint on interior models is equal to 7.1 ± 1.0 km. The error bar is on the order of 1 km, but Castillo-Rogez et al. (2012) show that the effect of cratering on altering the original shape could be up to 2 km. Thomas et al. (2018) recently provided a revised estimate of Phoebe's shape that is close to the Castillo- Rogez et al. (2012) results. In order to reduce the uncertainty on the extraction of Phoebe's equilibrium shape, we decomposed its surface in spherical functions as r(λ, ϕ) = N max n=0 n m=0 A n,m cos mλ + B n,m sin mλ P n,m (sin ϕ). (1) where r is the radius, λ is the longitude, ϕ is the latitude, P n,m is the associated Legendre polynomial, and A n,m and B n,m are the coefficients of the spherical functions. The long semi-axis a and (a − c) equilibrium shape for an oblate body was then deduced from (Rambaux et al. 2015): where A 0,0 corresponds to R, the mean surface radius. We started from the digital terrain models (DTM), called phoebe_quad64q.tab, of Gaskell (2013). From this model, we fit the surface to the 13th degree and order (see Table A.1). During the adjustment, we corrected the center of Phoebe's figure as well as the rotation of these axes so that the coefficients A 1,0 , A 1,1 , B 1,1 , A 2,1 , B 2,1 , and B 2,2 are minima. The correction on the figure center is then (−0.4, 0.4, 0.4) km, and the correction of the orientation angles is (0.62, 1.34, 91.98) degrees with the rotation matrices in angular sequence 1,2,3. The uncertainty column corresponds to the diagonal elements obtained in the covariance matrix. By applying Eq.
(2), we find a − c = 8.4605 km, with a very small formal uncertainty about meters. We estimated the residuals by computing the root mean square (RMS) as where r i is the radial distance deduced directly from the DTM and r(λ, ϕ) is the radius computed from the spherical functions. Figure 1 shows the residuals as a function of the maximum degree. At the 13th degree, the residuals are about 1 km, and we infer that the accuracy is about ±0.5 km. The equilibrium shape of a self-gravitating and rotating body is a spheroid with a nonzero degree for A 0,0 , A 2,0 , A 4,0 , A 6,0 , . . . (see, e.g., Chambat et al. 2010;Rambaux et al. 2015, and references therein). In Table A.1, these terms dominate, and, by comparing their values with the theoretical values of hydrostatic equilibrium (see Sect. 5 below), the A 4,0 and A 6,0 present large departures that can be related to non-hydrostatic contributions. Hence, we used the value of A 2,0 to obtain information on the average interior structure. The surface of Phoebe is heavily cratered, and the largest crater has a diameter of 101 km (Jason crater). A priori, this value is below the surface wavelength of 2πR/2 = 333 km, and A 2,0 should reflect only the global shape. However, as Phoebe is heavily cratered, the cumulative craters could yield a topography profile that mimics a second degree. In order to estimate this potential effect, we followed the statistical approach from Castillo-Rogez et al. (2012) and generated a synthetic surface for Phoebe using a size-frequency distribution with a slope of −2.26. For this statistical test, the depth-todiameter ratio was fixed to 0.2 for craters below 10 km, based on the topographic values of Phoebe's crater derived by Giese et al. (2006). For larger craters, we used a power-law function such that a large (130 km diameter) crater on Phoebe has a depth of 14.4 km similar to the value featured in Castillo-Rogez et al. (2012). For the rim and ejecta blanket, we used a parametrization developed for the Moon, where the height rim is 0.092 r crater and the ejecta at r > r crater follow a power law in −3.64. This parametrization led to larger mean values in the final shape, but the standard deviation of (a − c) is in agreement with the results from Castillo-Rogez et al. (2012). We obtain, on average, a difference equal to 0.55 km on A 2,0 , below the 2 km, for (a − c), found by Castillo-Rogez et al. (2012) This uncertainty is equivalent to the uncertainty determined by using the RMS. Therefore, in the following, we focus mainly on A 2,0 and use an uncertainty of 0.55 km. We also computed the RMS spectra for the topography, as defined in Ermakov et al. (2018). Figure 2 shows the variation of the RMS spectra in blue as well as the difference with the RMS spectra deduced by Ermakov et al. (2018). The two estimations are very close to each other, and the difference is below 10 −4 . As the DTM of Phoebe is spatially dense, we performed a powerlaw fit following Kaula's rule, and we obtain a slope of −1.4 ± 0.1 for the RMS spectra and −2.2 ± 0.2 for the variance spectra. The latter slope corresponds to the constant slope deduced for the other icy satellites (Nimmo et al. 2011). The interpretation of the slope's value involves physical mechanisms, such as cratering, as well as isostasy. These two processes can be described by studying high topographic degrees that are outside the scope of this study, which focuses on the interpretation of the low degree.

Hydrostatic equilibrium model
Phoebe's hydrostatic shape was used to constrain Phoebe's interior model. We used a forward approach in which we computed hydrostatic shapes for various interior models built on geophysical and formation arguments (see the next section). There are no analytical solutions for hydrostatic shapes of bodies with arbitrary density profiles, hence we used a numerical approach. Various approaches to this problem have recently been suggested (e.g., Hubbard 2013), and we used a direct numerical integration of the Clairaut equations developed for studying hydrostatic shapes of the Earth (Chambat et al. 2010), Ceres (Rambaux et al. 2015;Park et al. 2016), and large TNOs (Rambaux et al. 2017). These equations have been developed to an order that is a function of the geodetic parameter m 1 (Chambat et al. 2010), defined as the ratio between the centrifugal acceleration and the gravitational acceleration: where Ω is the spin angular velocity, G is the gravitational constant, R is the mean radius, and M is the mass of the body. At the first order, the Clairaut equation is written as (Kopal 1960;Lanzano 1974): where f 2 corresponds to the coefficient of the Legendre polynomial to the second order in the development of the equipotential surface s s(r, θ) = r(1 + f 2 (r)P 2 (cos θ)), γ = ρ(r)/ρ(r) is the ratio between the density of the layer r and the mean density at r, and θ is the co-latitude. Lanzano (1974) has developed these equations at order 3 by introducing the following coefficients: s(r, θ) = r(1+ f 2 (r)P 2 (cos θ)+ f 4 (r)P 4 (cos θ)+ f 6 (r)P 6 (cos θ)). (7) We checked the expressions of the equations and the boundary conditions from Lanzano (1974) by using Maple software. These equations can be solved in the general case by using an arbitrary density profile through the numerical scheme introduced by Chambat et al. (2010) and Rambaux et al. (2015) for the third order.

Geophysical models
Phoebe's low density, which was revised to 1642 ± 18 kg m −3 by Thomas et al. (2018), is indicative of a large fraction of volatiles, about 50-60 wt.%, assuming anhydrous rock. Phoebe's surface displays ice (Clark et al. 2005). The low (0.06) average surface albedo suggests that the ice is mixed with a dark component that may include silicates, carbon (graphite, amorphous), and metal compounds (e.g., de Sanctis et al. 2015). For forward modeling of the shape, we assumed a range of models with two and three layers, as illustrated in Fig. 3: a rocky mantle, an ice-rich crust, and a porous outer layer. The mantle density ranges from 1700 (ice mixed with rock) to 3200 kg m −3 . The lower bound represents an undifferentiated mixture of ice and rock, and the upper bound represents an anhydrous rock-metal mixture. The ice-rich layer density ranges from 1000 to 1300 kg m −3 to represent various impurity contents. The upper bound is based on the density of Ceres's crust (Ermakov et al. 2017). The megaregolith layer thickness ranges from 10 to 30 km, and its density is set to 800 and 900 kg m −3 to represent an up to 10-40% porosity created by large impacts and gardening over the long term (see Table 2).

Results
First, we computed the interior models assuming a homogeneous interior. Figure 4 shows the value of (a−c) over the range of possible mean density values within the error bars. The values computed from the Clairaut equation to the first order (black crosses) and the third order (black circle) are compared to the observed shape derived by Castillo-Rogez et al. (2012) (blue line) and that derived in the present paper (green line). The difference between the solutions at the first and third orders is about +0.6 km, that is, the first order underestimates the a − c shape by that much. The a − c value from the Clairaut equation for the uniform model is about 11 km, that is, more than two kilometers higher than the observed 8.46 ± 0.71 km. This indicates that Phoebe's interior is differentiated and/or it holds a fossil shape. Next, we explored a subset of five reference differentiated interior models, listed in Table 2, with one two-layer model and four three-layer models, where we varied the mantle density between 1700 and 3200 kg m −3 and the ice density between 1000 and 1300 kg m −3 . The crust layer for the three-layer model is fixed to 800 and 900 kg m −3 , and its thickness is fixed to 10 and 15 km (thicker porous crust is not consistent with the physical constraints available). The computed r f 2 (or, equivalently, A 2,0 ) is mapped against the ice and mantle densities in Figs. 5 and 6. These results show that Phoebe's interior is better Maps of r f 2 for two-layer assumption on Phoebe's interior structure. All cases include a rocky mantle and ice-rich layer whose densities are ranged (x-and y-axis, respectively). The observed value of A 2,0 to be matched is −6.11 ± 0.55 km and corresponds to the shaded band.
described as a stratified interior with a large contrast between the outer layer and the rocky interior, especially for a high-density ice-rich layer.
A two-layer model with an ice-rich crust and rocky mantle matches the observed A 2,0 = −6.11 ± 0.55 km for a realistic hydrated rock density based on a CM chondrite density (Fig. 5). The addition of an outer porous layer between 10-15 km, and assuming 10-40% porosity, restricts the space of models that can match the observations, as shown in Fig. 6. Solutions for mantle density lower than 2000 kg m −3 or a porous layer thicker than 20 km are not consistent with the observed mean density. Decreasing the thickness of that megaregolith increases the match between the computed r f 2 and the observed A 2,0 for a lower mantle density. The inferred rocky mantle density is Maps of r f 2 (or equivalently A 2,0 ) for different assumptions on Phoebe's interior structure. Crustal density is ranged on the y-axis and the rocky mantle density is ranged on the x-axis. The observed value of A 2,0 to be matched is −6.11 ± 0.55 km and represented by the shaded bands. All cases include a rocky mantle and ice-rich layer whose densities are ranged (x-and y-axis, respectively). Models differ by the characteristics of an upper low density layer assumed as a megaregolith formed as a consequence of large impacts and long-term gardening with a density of 800 kg m −3 (right column) and 900 kg m −3 (left column) and varying thickness, from top to bottom: 10 km and 15 km. The white part are not consistent with available constraints (Phoebe's density in particular).
between 2000 kg m −3 and 2400 kg m −3 . For Phoebe's A 2,0 value, which is close to 6.11 km, observations are better matched with a mantle density between 2000-2400 kg m −3 , if the crustal density is about 800 kg m −3 , or between 2000-2500 (2800 kg m −3 for few extreme cases) for a crustal density of 900 kg m −3 . The density of the icy crust itself has little effect on the calculated r f 2 , hence it cannot be constrained. Figures 5 and 6 show that the second degree shape of Phoebe can be explained by a rotating body at 9.27 h at hydrostatic equilibrium. However, it is not possible to rule out a fossil shape. Figure 7 shows the r f 2 as a function of the spin period. We note that for a faster body, at 8.5 h, the r f 2 drops by 1 km and the curves coinciding with A 2,0 are shifted toward the left in Figs. 5 and 6. On the other hand, for a slower spinner, around 10 h, A 2,0 increases by 1 km and the A 2,0 curves shift toward the right in Figs. 5 and 6. The quantity A 4,0 can fit interior models for the current spin rate at 7.27 h, but A 2,0 presents a large discrepancy. The quantity A 6,0 does not fit any models. Therefore, these two last quantities present large non-hydrostatic contributions.

Conclusions
We determined a spherical harmonic decomposition of the topography of Phoebe to the 13th degree and order. From this decomposition, we extracted the zonal even degrees describing the global shape for comparison with forward modeling.
This study emphasized the importance of using the Clairaut equations to at least the second order when modeling the shape of fast rotators (<15 h) for comparison with observations of bodies large enough to be hydrostatically relaxed. This definitely rules out a uniform interior for Phoebe. Instead, it indicates that it has an interior stratified in a low-density crust (<1000 kg m −3 ) that is likely porous, with a thickness smaller than 20 km, and a mantle with a density between 2000 and 2400 kg m −3 with a bounding case at 2800 kg m −3 for a 10 km crust at 900 kg m −3 . These results support the earlier suggestion by Castillo-Rogez et al. (2012) that Phoebe was subject to at least partial melting, with the separation of an ice-rich shell and a period during which rock interacted with ice, an evolution that might be similar to carbonaceous chondrite parent bodies (Castillo-Rogez & Young 2016). The interior models inferred for Phoebe, with an ice-silicate mixture mantle density, support that type of evolution model. The low-density shell inferred from this modeling may indicate an enrichment in ice and/or the presence of a significant fraction (10-40%) of macroporosity resulting from Phoebe's intense collisional history. Phoebe is believed to have accreted other irregular asteroids present in its neighborhood, as illustrated by the large density of craters found across its surface (Porco et al. 2005). Values of r f 2 as function of the spin period for one interior models inside the hach band of Fig. 6a. The estimated A 2,0 is represented in vertical band and the horizontal line represents the actual spin period of Phoebe. This plot shows that the A 2,0 value is consistent with rotation period between 9 to 9.75 h red band. The lightsalmon band represents variation of ±1 km used in the discussion (see main text).
A major difference from C-type bodies, though, is that Phoebe exhibits evidence for ice on its surface (Clark et al. 2005). Ice is not stable in the warmer conditions in the main belt of asteroids. However, ice could be buried a few meters below the surface (Schorghofer 2008(Schorghofer , 2016 or could have been eroded by impacts, as suggested by Marsset et al. (2020) for Pallas based on its density of about 2700 kg m −3 . Hence, Phoebe might represent a mostly pristine example of a C-type body and, in turn, provide constraints on the volatile composition of these bodies based on the results of the Cassini mission (e.g., Clark et al. 2005).