Evidence for differentiation of the most primitive small bodies

Dynamical models of Solar System evolution have suggested that P-/D-type volatile-rich asteroids formed in the outer Solar System and may be genetically related to the Jupiter Trojans, the comets and small KBOs. Indeed, their spectral properties resemble that of anhydrous cometary dust. High-angular-resolution images of P-type asteroid (87) Sylvia with VLT/SPHERE were used to reconstruct its 3D shape, and to study the dynamics of its two satellites. We also model Sylvia's thermal evolution. The shape of Sylvia appears flattened and elongated. We derive a volume-equivalent diameter of 271 +/- 5 km, and a low density of 1378 +/- 45 kg.m-3. The two satellites orbit Sylvia on circular, equatorial orbits. The oblateness of Sylvia should imply a detectable nodal precession which contrasts with the fully-Keplerian dynamics of the satellites. This reveals an inhomogeneous internal structure, suggesting that Sylvia is differentiated. Sylvia's low density and differentiated interior can be explained by partial melting and mass redistribution through water percolation. The outer shell would be composed of material similar to interplanetary dust particles (IDPs) and the core similar to aqueously altered IDPs or carbonaceous chondrite meteorites such as the Tagish Lake meteorite. Numerical simulations of the thermal evolution of Sylvia show that for a body of such size, partial melting was unavoidable due to the decay of long-lived radionuclides. In addition, we show that bodies as small as 130-150 km in diameter should have followed a similar thermal evolution, while smaller objects, such as comets and the KBO Arrokoth, must have remained pristine, in agreement with in situ observations of these bodies. NASA Lucy mission target (617) Patroclus (diameter~140 km) may, however, be differentiated.


Introduction
The Cybele region, at the outer rim of the asteroid belt (3.27-3.7 au), is essentially populated by P-and D-type asteroids and to a lesser extent by C-type bodies (DeMeo & Carry, 2013. P-and D-type asteroids are thought to have formed in the outer Solar System (beyond 10 au), among the progenitors of the current Kuiper Belt, and to have been implanted in the inner Solar System (asteroid belt, Lagrangian Points of Jupiter) following giant planet migrations (see Levison et al., 2009;Vokrouhlický et al., 2016). This implies that the P/D-type main belt asteroids and the Jupiter Trojans could be compositionally related to outer small bodies such as Centaurs, short-period comets, and small (D≤300km) Kuiper Belt Objects (KBOs). This dynamical scenario is currently supported by the Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under program 073.C-0851 (PI Merline), 073.C-0062 (PI Marchis), 085.C-0480 (PI Nitschelm), 088.C-0528 (PI Rojo), 199.C-0074 (PI Vernazza).
The reduced and deconvolved SPHERE images are available at the CDS via anonymous ftp to http://cdsarc.u-strasbg. fr/ or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+ A/xxx/Axxx Deceased similarity in size distributions between the Jupiter Trojans and small KBOs (Fraser et al., 2014) as well as the similarity in terms of spectral properties between P/D-type main belt asteroids, the Trojans of Jupiter and comets (Vernazza et al., 2015;Vernazza & Beck, 2016).
Overall, the outer Solar system is of tremendous interest, as it is recognized as the least processed since the dawn of the Solar System and thus the closest to the primordial materials from which the Solar System formed (e.g., McKinnon et al., 2020). This is currently supported by the analysis of the spectral properties of P/D-type main belt asteroids, Jupiter Trojans, and comets that reveal a surface composition compatible with that of anhydrous chondritic porous interplanetary dust particles (CP IDP, see Vernazza et al., 2015). The CP IDPs are currently seen among the available extra-terrestrial materials as the closest to the starting ones (Bradley, 1999). In particular, based on its albedo and visible, near-and mid-infrared spectrum, it is now well established (see Vernazza et al., 2013) that the aqueously altered Tagish Lake meteorite cannot be representative of the surface composition of D-type asteroids nor that of most P-type asteroids as suggested earlier (Hiroi et al., 2001). As such, CP IDPs are currently the most likely analogs of the refractory materials present at the surface of P/D type asteroids.
With an estimated diameter of nearly 280 km (Carry, 2012, and reference therein), (87) Sylvia is the largest body in the Cybele region and more generally the largest P/D-type asteroid in the inner (≤5.5 au) Solar System. Its surface composition is fully consistent with that of anhydrous chondritic porous IDPs (Vernazza et al., 2015;Usui et al., 2019). Furthermore, it is the first asteroid around which two satellites were discovered (Marchis et al., 2005). Because of its rather large angular size at opposition (≈0.14 ) and its two moons, Sylvia is an ideal target for high angular-resolution adaptive-optics observations as these allow an accurate characterization of its 3D shape and of the mass of the system, hence of its bulk density. Furthermore, the two moons allow probing the internal structure and in particular the harmonics of the gravity field (at least the lowest-order gravitational moment, the quadrupole J 2 ). This information, in turn, probes the distribution of material inside Sylvia, i.e., whether its interior is homogeneous or differentiated.
Several authors have studied Sylvia's dynamical system (Marchis et al., 2005;Winter et al., 2009;Frouard & Compère, 2012;Fang et al., 2012;Beauvalet & Marchis, 2014;Berthier et al., 2014;Drummond et al., 2016). While these studies agreed on the mass of Sylvia and its low density (around 1,300 kg·m −3 ), there is no consensus regarding its gravitational potential (of which only the quadrupole J 2 has been studied, Table 1). The latter is, however, a direct consequence of the internal structure of Sylvia.
To constrain the bulk density of Sylvia and its internal structure, we observed it as part of our imaging survey of D≥100 km main-belt asteroids (ID 199.C-0074, PI P. Vernazza, see Vernazza et al., 2018Vernazza et al., , 2020Viikinkoski et al., 2018;Fétick et al., 2019;Carry et al., 2019;Hanuš et al., 2019Hanuš et al., , 2020Marsset et al., 2020;Yang et al., 2020;Ferrais et al., 2020). We imaged Sylvia over two apparitions, separated by a year, throughout its rotation at high angular resolution with the SPHERE extreme adaptive-optics (AO) IRDIS and ZIMPOL cameras (Thalmann et al., 2008;Dohlen et al., 2008;Beuzit et al., 2019) mounted on the European Southern Observatory (ESO) Very Large Telescope (VLT). Beyond contributing to our understanding of the formation and evolution of the Solar System's most primitive bodies, the present study provides the context for the future insitu exploration of primitive P/D-type bodies by NASA's Lucy mission to the Jupiter Trojans.

High-angular-resolution imaging
The data used in the present study to extract the position of the moons comprises all the high angular resolution images of Sylvia taken with the Hubble Space Telescope (HST) and ground-based telescopes equipped with adaptive-optics (AO) cameras: Gemini North, ESO VLT, W. M. Keck, and SOR (Drummond et al., 2016). The data span 51 different epochs, with multiple images each, over 17 years from February 2001 to November 2018. For the reconstruction of Sylvia's 3D shape, however, only the images with the highest resolution were used, i.e., those acquired with VLT/SPHERE/ZIMPOL (Section 3).
The images from the HST were obtained with the second Wide Field and Planetary Camera (WFPC2, Holtzman et al., 1995). The images from the VLT were acquired with both the first generation instrument NACO (NAOS-CONICA, Lenzen et al., 2003;Rousset et al., 2003) and SPHERE (Spectro-Polarimetric High-contrast Exoplanet REsearch, Fusco et al., 2006;Beuzit et al., 2019), the second generation extreme-AO instrument designed for exoplanet detection and characterization. The images acquired with SPHERE were taken with both the IRDIS (InfraRed Dual-band Imager and Spectrograph, Dohlen et al., 2008) andZIMPOL (Zurich Imaging Polarimeter, Thalmann et al., 2008) sub-systems. Images taken at the Gemini North used the NIRI camera (Near InfraRed Imager, Hodapp et al., 2003), fed by the ALTAIR AO system (Herriot et al., 2000). Finally, observations at Keck were acquired with the guiding camera of NIRSPEC in 2001(McLean et al., 1998) and NIRC2 (Near-InfraRed Camera 2, van Dam et al., 2004Wizinowich et al., 2000) later on.
To reduce the AO-imaging data, a standard data processing protocol (sky subtraction, bad-pixel identification and correction, and flat-field correction) was followed using in-house routines developed in Interactive Data Language (IDL) (see Carry et al., 2008). The images were then deconvolved with the Mistral algorithm (Fusco, 2000;Mugnier et al., 2004) to restore their optimal angular resolution (see Fétick et al., 2019, for details). Separately, the reduced images were processed to subtract the bright halo surrounding Sylvia to enhance the detectability of the satellites (see Pajuelo et al., 2018; 2019, for details).

Optical lightcurves
We used the 40 lightcurves from Kaasalainen et al. (2002) to create a convex 3-D shape model of Sylvia 1 , compiled from the Uppsala Asteroid Photometric Catalog 2 (Lagerkvist & Magnusson, 2011). We also compiled 11 lightcurves acquired by amateur astronomers within the Courbes de rotation d'astéroïdes et de comètes database (CdR 3 ).
In addition to these data, we acquired 12 lightcurves using the 60 cm André Peyrot telescope mounted at Les Makes observatory on La Réunion Island, operated as a partnership among Les Makes Observatory and the IMCCE, Paris Observatory, and 7 lightcurves with the 60 cm TRAPPIST telescopes located at La Silla Observatory in Chile and the Oukaïmeden observatory in Morocco (Jehin et al., 2011). Finally, we extracted 51 lightcurves from the data archive of the SuperWASP survey (Pollacco et al., (Parley et al., 2005;Grice et al., 2017). In summary, a total of 121 lightcurves observed between 1978 and 2017 (Table A.1) were used in this work and are presented in Figure A.1.

Stellar occultations
Nine stellar occultations by Sylvia have been recorded since 1984, mostly by amateur astronomers during the last decade (see Mousis et al., 2014;Dunham et al., 2016;Herald et al., 2020). We converted the timings of disappearance and reappearance of the occulted stars 4 into segments (called chords) on the plane of the sky, using the location of the observers on Earth and the apparent motion of Sylvia following the recipes listed in Berthier (1999). Only five stellar occultations had multiple chords that could be used to constrain the size and apparent shape of Sylvia ( Figure B.1). For the January 2013 and October 2019 occultations, several observers reported secondary events due to the occultation of the stars by either Romulus or Remus (Berthier et al., 2014;Vachier et al., 2019). We thus used the relative positions between Sylvia and its satellites at the time of the occultations to constrain their mutual orbits. We list the observers of the occultations in Table B.1.

Sylvia's 3D shape
We fed the All-Data Asteroid Modeling (ADAM) algorithm with all SPHERE/ZIMPOL images, lightcurves, and stellar occultations to determine the spin and 3D shape of Sylvia (Viikinkoski et al., 2015a). Optical lightcurves are often required for ADAM to constrain the regions not imaged and to stabilize the shape optimization. The procedure is similar to that published in our previous studies with SPHERE, and we refer to these for more details (e.g., Viikinkoski et al., 2015b;Marsset et al., 2017;Vernazza et al., 2018). We used the sidereal rotation period and the spin-axis coordinates of Sylvia from the literature (Kaasalainen et al., 2002;Hanuš et al., 2013Hanuš et al., , 2017 as input values to ADAM. We further model the shape by using the Multi-resolution PhotoClinometry by Deformation (MPCD) method (Jorda et al., 2016), following the procedure of our previous works (e.g., Ferrais et al., 2020). MPCD gradually deforms the vertices of a previous mesh (here ADAM model) to minimize the difference between the observed images and synthetic images of the model (Jorda et al., 2010). Both models only present marginal differences (Table 2), and in what follows, we report on the MPCD model. The derived shape is in essence similar to that based on lower angular resolution images (Berthier et al., 2014;Hanuš et al., 2017). A striking feature of Sylvia is its remarkable elongated shape (Table 2, Figure 2). To put its peculiar shape into context, we measured the tri-axial diameters of 103 asteroids larger than 100 km from their shape models 5 and compiled their rotation periods from the Planetary Data System (Harris et al., 2017). Sylvia appears to be more elongated and to spin faster than most asteroids larger than 100 km ( Figure 2). In particular, the population of asteroids with satellites stands out from the population of singletons, with a median ratio of equatorial diameters a/b of 1.37 (vs 1.16 for the background) and a median rotation period of 5.5 h (vs 7.9 h for the background). We compared the a/b ratio and rotation period distributions of asteroids with and without satellites using the Kolmogorov-Smirnov test. The p-values for both are 5 · 10 −4 and 7 · 10 −5 respectively. We thus conclude that the distributions are different above the 99.5% confidence level. While asteroids with a diameter smaller than about 15 km are subject to YORP spin-up and surface re-arrangment (Walsh et al., 2008;Vokrouhlický et al., 2015), Sylvia is too large to have been affected. The origin of this difference is thus unclear. It may result from the impact at the origin of the satellite formation (Margot et al., 2015), or alternatively, satellites may be more stable around elongated bodies (see Winter et al., 2009).

Dynamics of the system
The two satellites orbit Sylvia on equatorial, circular, and prograde orbits ( Table 3). The root mean square (RMS) residual between the observations and the computed positions is only 9.3 mas, that is within the pixel size of most observations. The positive occultations by both satellites in October 2019 provide a practical estimate of the reliability of the orbital solution, as the two satellites were detected at only 5 mas from the positions we predicted . It also highlights the importance of accurate ephemerides to prepare the observation of the occultation by placing observers on the path of the satellites.
The mass of Sylvia is constrained with an uncertainty of less than 1% : (1.44 ± 0.01) × 10 19 kg, thanks to the long baseline of observations. Combined with our volume-equivalent diameter estimate (271 ± 5 km, see above), the density of Sylvia is found to be 1378 ± 45 kg·m −3 , reminiscent of that of other large asteroids with a surface composition consistent with that of IDPs (C/P/D types, see Carry, 2012;Vernazza et al., 2015). We present in Figure 3 all the possible bulk compositions of Sylvia, consid-  ering a mixture of rocks and ices, with voids. We use a range of density from 2200 to 3000 kg·m −3 for the rock phase (from rocks with organics to the density of the silicate phase reported by the Stardust mission, Brownlee et al., 2006). The density of Sylvia implies the presence of both ices and macroporosity in its interior.
We estimate the masses of Romulus and Remus at (1.4 ± 1.2) × 10 15 kg and (7.8 ± 7.3) × 10 14 kg (i.e., effectively upper limits), very close to the masses reported by Fang et al. (2012). Assuming a similar albedo for Sylvia and its two moons, their magnitude differences with Sylvia (Table 3) imply diameters of 15 +10 −6 km and 10 +17 −6 km for Romulus and Remus (consistent with occultation chords, Berthier et al., 2014), hence densities of 790 ± 680 and 1480 ± 1400 kg·m −3 . The density of both satellites is loosely constrained and similar to that of Sylvia.
Finally, we note that the best orbital solutions are obtained for the smallest quadrupole J 2 (Figure 4), tending toward J 2 =0 (i.e., fully Keplerian orbits over the 19 years baseline). Although there are orbits fitting the data within 1 σ of the observations, their residuals are systematically larger.

Implication for the internal structure
Under the assumption of a homogeneous density in the interior, the shape of Sylvia implies a J 2 of 0.082 ± 0.005 (computed with SHTOOLS 6 , see Wieczorek & Meschede, 2018). This values contrasts with the null J 2 determined dynamically (Section 4). This discrepancy reveals an inhomogeneous density distribution inside Sylvia and hints at a more spherical mass concentration than suggested by Sylvia's oblate and elongated 3D shape. This implies a denser, more spherical core, surrounded by a less dense envelope. Based on similar considerations, similar internal structures have recently been proposed for other large P-types such as the Cybele (107) Camilla (Pajuelo et al., 2018) and the Jupiter Trojan (624) Hektor .
This differentiated structure is at odds with the IDP-like spectral properties, evidence for an absence of both thermal metamorphism and aqueous alteration. This suggests that partial differentiation occurred, limited by the insufficient amount of heat generated by radionuclides which did not propagate to the Table 3: Orbital elements of the satellites of Sylvia, expressed in EQJ2000, obtained with Genoid: orbital period P, semi-major axis a, eccentricity e, inclination i, longitude of the ascending node Ω, argument of pericenter ω, time of pericenter t p . The number of observations and RMS between predicted and observed positions are also provided. Finally, we report the mass of Sylvia M Sylvia , the mass of Romulus M Romulus , the mass of Remus M Remus , their apparent magnitude difference ∆m with Sylvia, the ecliptic J2000 coordinates of the orbital pole (λ p , β p ), the equatorial J2000 coordinates of the orbital pole (α p , δ p ), and the orbital inclination (Λ) with respect to the equator of Sylvia. Uncertainties are given at 3-σ.

Romulus
Remus 10.3 ± 2.1  The difference between asteroids with and without satellites is striking.
Building upon the work of Neveu & Vernazza (2019), we model the thermal and internal structure history of Sylvia. The evolution of internal temperatures and structure is computed nu-merically using a one-dimensional code (Desch et al., 2009). Sylvia is assumed to be made of rock (idealizing a mixture of refractory materials such as silicates, metals, and organic material), water ice, and voids (the macroporosity). The mass is distributed assuming spherical symmetry over 200 grid zones initially evenly spaced in radius. The internal energy in each grid zone is computed from the initial temperature using equations of state for rock and ice. Material is never hot enough in our simulations for rock-metal differentiation, which is neglected. Initial radionuclide abundances are provided in Table 1 of Neveu & Vernazza (2019). Simulations start once Sylvia is fully formed, neglecting the progressive accretion of material over time. Because of this and the near-absence of short-lived radionuclide heating given the assumed formation time, Sylvia's simulated early evolution is cold. The implementation of instantaneous differentiation in the central regions that warm above 273 K rests on the assumption that sufficiently large rock grains settle via Stokes flow on timescales smaller than one time step.
Sylvia is assumed to accrete homogeneously at 60 K, 6 million years (My) after the formation of Ca-Al-rich inclusions (consistent with a surface without aqueous alteration; Neveu & Vernazza, 2019). It is the equilibrium temperature for an albedo 0.05 at a distance of 17-18 au (i.e., the postulated accretion distance of KBOs, Morbidelli & Nesvornỳ, 2020) from the Sun with 70% of the present-day luminosity. The surface temperature is instantaneously raised to 148 K (the present-day equilib-Article number, page 5 of 27 A&A proofs: manuscript no. main rium temperature) at the time of heliocentric migration. We test different timings, from a late planet migration (hundreds of My) such as the one described by the Nice model Morbidelli et al., 2005;Gomes et al., 2005;Levison et al., 2009) to an early dynamical instability occurring a few My after the dissipation of the gas disk (Nesvorný et al., 2018;Clement et al., 2019), and find that the timing of implantation into the asteroid belt seldom affects its structure or peak temperature. In the results below, the time of migration is set to 6 My after formation (i.e., 12 My after Ca-Al-rich inclusions).
The thermal structure is determined by balancing conductive heat transfer with primarily radiogenic heating by 26 Al, 40 K, 232 Th, 235 U, and 238 U, using a finite-difference method and a 50yr time step, for 5 billion years (Gy). Thermal conductivities depend mainly on porosities ( Figure 6), but also on composition and temperature (Desch et al., 2009). Porosity is allowed to compact at rates determined from material viscosities as described in Neveu & Rhoden (2017). Sylvia's bulk density constrains the void porosity and rock volume fractions to about 40-55% and 20-35%, respectively, mainly depending on the rock density ( Figure 3).
Convection is generally neglected as it is assumed that the postulated porous, rock-rich internal structure for Sylvia is not prone to fluid or ice advection. In the ice-rich case ( Figure 5, right panel), solid-state convection is allowed to occur but does not, because the critical Rayleigh number is never exceeded. Although this was not simulated, in the percolation case below we expect convection to be possible in the central region rich in liquid water, until this region refreezes. In all simulations, volume changes due to water melting or freezing are neglected.
The viscosity of ice-rock mixtures, used to compute the Rayleigh number and pore compaction, is calculated following Roberts (2015). Above 30 vol.% ice, it is equated to the viscosity of ice. Below 30% ice volume fraction, as a first approximation, it is set to the geometric mean of the rock and ice viscosities at given grain size, stress (equated with hydrostatic pressure), and temperature. Roberts (2015) noted that this approximation tends to underestimate mixture viscosities relative to extrapolations of laboratory measurements. The ice and rock flow laws adopted in the model are the composite rheology of Goldsby & Kohlstedt (2001) and the dry diffusion creep flow law for olivine of Korenaga & Karato (2008), respectively.
The key factor governing thermal evolution in these simulations is porosity Φ, which decreases thermal conductivities k of rock-ice mixtures, of order 1 W m −1 K −1 (Desch et al., 2009, and references therein), by up to two orders of magnitude. The adopted thermal conductivity-porosity relationship (Shoshany et al., 2002) is derived from Monte-Carlo modeling of porous cometary ice: k is decreased with increasing Φ via multiplication by a factor (1 -Φ /0.7) nΦ+0.22 . This relationship is shown in Figure 6. We adopt n=4.1 (grey curve) for the canonical simulation with 52% porosity, following the discussion in Shoshany et al. (2002). We set n=3.5 (black curve) for the simulation with 60% porosity; choosing n=4.1 would result in more heating and compaction than shown in Figure 5. There is a wide spread in the Monte Carlo results of Shoshany et al. (2002) in how Φ affects k, with some of their models suggesting a lesser effect of porosity. Conversely, independent measurements of porous silica aggregates (Krause et al., 2011) and extrapolated results from models of highly porous aggregates (Arakawa et al., 2017) both suggest similar or slightly higher decrease factors due to porosity ( Figure 6) than the relationships we have assumed. Thus, the assumed relationships seem to adequately represent the state of understanding of how porosities of 50 to 60 vol.% decrease thermal conductivities.
Crucially, to be compatible with Sylvia's observed anhydrous surface and low J 2 despite an oblate shape, time-evolution simulations with this model constrain the volume fraction of water to be low, less than 40 vol% relative to rock or 15 vol% overall. For higher volume fractions, ice grains tend to become adjoined and control the mechanical properties of the interior. In that case, the material viscosity is assumed equal to that of water ice (Goldsby & Kohlstedt, 2001), and any void porosity rapidly decreases as the interior warms due to radiogenic heating and porosity insulation. Instantaneous ice-rock differentiation happens first once the central (warmest) regions warm above 273 K. It then proceeds outward if the interior keeps warming. This yields a gravitationally unstable structure: the topmost undifferentiated layers are denser than underlying layers, which are poorer in rock. Once Sylvia is differentiated out to more than   The orange curve is a thermal conductivity in W m −1 K −1 (rather than a decrease factor). Data points of a given color are from the same source as the fit curve of that color.
half its radius, differentiation is assumed to proceed by gravitational (Rayleigh-Taylor) instabilities: layers overturn if their viscosity is below a threshold that corresponds to T≈150 K (Rubin et al., 2014). Since Sylvia's post-migration surface temperature is warmer, 148 K, differentiation out to the surface is essentially inevitable. This ought to result in evidence of surface water, either as ice or as mineral hydration, as observed on aster-oids linked to carbonaceous chondrites (Rivkin & Emery, 2010;Campins et al., 2010). This is inconsistent with Sylvia's anhydrous, IDP-like surface composition. Ice-rock differentiation can be prevented if ice dominates the volume fraction ( Figure 5, right column), since in that case the combination of low rock (i.e., radionuclide) content and low insulating void porosity results in a cold interior in which ice never melts. However, in such a homogeneous interior, the mass distribution should result in a higher J 2 than observed given Sylvia's oblate shape.
It follows that to retain a pristine anhydrous external envelope, Sylvia's water volume fraction must be low (consistently with observations of the comet 67P nucleus by the Rosetta spacecraft; Pätzold et al., 2019;Choukroun et al., 2020) so that interior solids are less prone to deformation, inhibiting both porosity compaction and instability-driven differentiation. As a canonical case, we assume an interior comprised of 52 vol.% void porosity, 35% rock of density 3075 kg m −3 , and 13% water ( Figure 5, central column). The void porosity decreases the interior's thermal conductivity by a factor of ≈15 relative to a compact rock-ice mixture. This favors accumulation of radiogenic heat, melting ice in the central regions after ≈0.15-0.2 Gy. At such low water volume fraction and high void porosity, it is sensible for liquid water to percolate downward in pores without significantly disturbing the remaining rock-void porosity structure, since rock grains already tend to be adjoined (see Fig. 3a,b of Neumann et al., 2020, for a pictorial description). This would result in a three-layer structure ( Figure 5, bottom central panel): a central region where the porosity has been filled with percolated water, surrounded by a porous layer free of water, and a primordial outer layer that remains too cold for ice to melt. Our thermal evolution simulations do not explicitly track percolation. The example interior structure of Figure 5, bottom central panel is obtained by manually moving the mass of water that is liquid at 0.2 Gy to fill porosity at the center, and assuming that in the middle layer, the empty volume left behind by the displaced water is compacted (so that the void porosity remains 52% in this layer). This compaction leads to a volume-averaged diameter decrease from 283.4 km to 271.4 km, the observed value. The spherical central mass distribution in this three-layer model implies a J 2 of 5 · 10 −5 only, consistent with the observed dynamics of both satellites.
Although this is not a unique solution, assuming lower rock densities (Figure 3) requires increasing the rock volume fraction at the expense of ice or void porosity so as to keep matching Sylvia's bulk density. Neither the bulk ice volume fractions nor the bulk void porosity are likely to be much lower than assumed given the need to invoke the migration of melted ice, enabled by the insulating effect of porosity, to explain a more spherical mass distribution (low J 2 ) than suggested by Sylvia's oblate shape.
Another, less likely explanation for a spherical mass distribution inside Sylvia is central compaction of the rock. Although we assume a rather low viscosity for rock-ice mixtures, Sylvia's relatively low gravity (lithostatic pressures) precludes compaction below 900 K. The required thermal insulation could be achieved with a bulk void porosity as low as 60% ( Figure 5; left column). However, such a hot evolution would result in advection and, likely, outward outgassing of water (Prialnik & Podolak, 1999;Young et al., 2003), which are not captured in these simulations and would cool the interior.
We thus deem percolation of water in the deep interior as being the likelier explanation for Sylvia's low J 2 despite its oblate shape. This implies that, unlike the pristine outer layers comparable to the CP IDPs, the innermost region may be instead analogous to hydrated material exemplified by chondritic smooth IDPs or perhaps the Tagish Lake meteorite (Fujiya et al., 2019). The minimum body diameter for such percolation to take place (holding all other quantities constant) is between 130 and 150 km, implying that even objects as small as Patroclus (diameter ≈140 km, see Hanuš et al., 2017), target of NASA's Lucy mission, may have experienced a low degree of central liquid water percolation.

Conclusions
We used newly acquired high-angular resolution imaging observations of (87) Sylvia with the SPHERE instrument on the ESO VLT, along with archival images, lightcurves, and stellar occultations to reconstruct its 3D shape and to constrain the orbital properties of its two moons.
We find that Sylvia possesses a low density of 1378 ± 45 kg·m −3 , similar to that of other large C/P/D asteroids whose surface composition is mostly consistent with that of anhydrous interplanetary dust particles. Sylvia spins fast and is oblate and elongated, a property shared by most 100+ km multiple asteroids, contrasting with the physical properties of large asteroids without satellites.
The orbits of the two satellites is in apparent contradiction with the oblate shape of Sylvia: the two orbits do not show the nodal precession expected from the shape. We interpret it as an evidence for a central spherical mass concentration, due to water percolation over millions of years triggered by long-lived radionuclides. This long lasting heating episode allowed for partial differentiation, the outer shell of Sylvia remaining pristine. It follows that even the most primitive small bodies with diameters larger than 150 km did not avoid thermal processing leaving only their outermost layers intact. Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation.
This research has made use of the Keck Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration.
The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.
We thank the AGORA association which administrates the 60 cm telescope at Les Makes observatory, La Reunion island, under a financial agreement with Paris Observatory. Thanks to A. Peyrot, J.-P. Teng for local support, and A. Klotz for helping with the robotizing.
B Our colleague and co-author M. Kaasalainen passed away while this work was carried out. Mikko's influence in asteroid 3D shape modeling has been enormous. This study is dedicated to his memory.
This paper makes use of data from the DR1 of the WASP data (Butters et al., 2010) as provided by the WASP consortium, and the computing and storage facilities at the CERIT Scientific Cloud, reg. no. CZ.1.05/3.2.00/08.       (