# Conditions for Accretion Favoring an Unmelted Callisto and a Differentiated Ganymede

Yannis Bennacer, Olivier Mousis, Marc Monnereau, Vincent Hue, Antoine Schneeberger

► To cite this version:

Yannis Bennacer, Olivier Mousis, Marc Monnereau, Vincent Hue, Antoine Schneeberger. Conditions for Accretion Favoring an Unmelted Callisto and a Differentiated Ganymede. The Planetary Science Journal, 2025, 6, pp.138. 10.3847/psj/add719. hal-05390740

HAL Id: hal-05390740

https://hal.science/hal-05390740v1

Submitted on 1 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Conditions for Accretion Favoring an Unmelted Callisto and a Differentiated Ganymede

Yannis Bennacer<fnref n="1" />,<fnref n="4" /> Ⓓ, Olivier Mousis<fnref n="1" />,² Ⓓ, Marc Monnereau<fnref n="3" /> Ⓓ, Vincent Hue<fnref n="1" /> Ⓓ, and Antoine Schneeberger<fnref n="1" /> Ⓓ

Received 2024 November 21; revised 2025 May 6; accepted 2025 May 7; published 2025 June 10

## Abstract

Analysis of Callisto's moments of inertia, derived from Galileo's gravity data, suggests that its structure is not fully differentiated. This possibly undifferentiated state contrasts sharply with the globally molten state inferred in its counterpart, Ganymede, and poses unique challenges to theories of the formation and evolution of the Galilean moons. During their formation, both moons experienced multiple heating mechanisms, including tidal heating, radiogenic heating from short-lived radionuclides, accretional heating from impacts, and heat from the surrounding circumplanetary disk (CPD). Our study investigates the optimal conditions required to account for Callisto's partially differentiated state in contrast to Ganymede's complete differentiation. We investigate crucial accretion parameters, such as the timing of accretion onset, the duration of accretion, and the impactor size distribution. We find that the observed dichotomy between Ganymede and Callisto can be attributed to similar formation conditions, assuming an identical impactor size distribution and composition in the Jovian CPD. The key differences in the formation of Ganymede and Callisto are the disk temperature at their respective formation locations and their final radii. Our results indicate that both moons accreted gradually over more than 2 Myr, concluding at least 5.5 Myr after the formation of calcium-aluminum-rich inclusions in the protosolar nebula. Our model demonstrates that Callisto can remain undifferentiated despite accreting a substantial influx of kilometer-sized impactors—potentially contributing up to 30% of the total mass inflow—while still allowing for the complete differentiation of Ganymede.

Unified Astronomy Thesaurus concepts: Solar system (1528); Natural satellite formation (1425); Ganymede (2188); Callisto (2279); Planetesimals (1259)

## 1. Introduction

The remarkable diversity of giant planet moons in our solar system provides a unique opportunity to investigate the origin and evolution of ocean worlds, both within these systems and in the broader context of exoplanetary science. A particularly fascinating feature of the Jovian system is the disparity in the internal structures of its moons. Specifically, the outermost moon, Callisto, shows no sign of global melting, in stark contrast to the fully differentiated interiors of the inner moons —Io, Europa, and Ganymede.

Analysis of Callisto's moments of inertia, based on Galileo gravity data ($C/MR^2 = 0.3549$; J. Anderson et al. 2001), indicates that its structure is not fully differentiated. This incomplete separation of ice and rock suggests that Callisto has remained cold from its accretion phase to the present day. The undifferentiated state of Callisto, in contrast to the globally molten state observed in its counterpart Ganymede (M. G. Kivelson et al. 1998; J. Anderson et al. 2001), imposes unique constraints on theories regarding the formation and evolution of the Galilean moons.

Radiogenic heating, especially from short-lived radionuclides such as ²⁶Al, is thought to be a primary driver of ice-rock separation, possibly triggering thermal runaway early in the history of giant planet satellites (J. S. Lewis 1971). Accretional heating, the heat generated by impacts during

accretion, has also been proposed as a mechanism for melting the ice within the ice-rock mixture (G. Schubert et al. 1981). The contribution of accretional heating depends on factors such as the size and velocity of the impactors within the Jovian circumplanetary disk (CPD) and the growth timescale of the protosatellite. The undifferentiated nature of Callisto is often attributed to late formation, a long accretion timescale, and the predominance of small accreting particles (less than ~100 m in size; A. Barr & R. Canup 2008; A. C. Barr et al. 2010). However, larger impactors may have been present in Callisto's formation region (P. R. Estrada et al. 2009; A. Barr & R. Canup 2008; K. Batygin & A. Morbidelli 2020). The size distribution of impactors during accretion is uncertain and varies widely between formation models. Two main scenarios are considered for the formation of satellites: pebble accretion and oligarchic growth. In pebble accretion, gas drag favors the accumulation of centimeter-to-meter-sized particles, leading to satellite formation (T. Ronnet et al. 2017). In contrast, the oligarchic growth scenario assumes kilometer-sized particles in the CPD that grow through collisions, eventually reaching the size of the Galilean satellites (P. R. Estrada et al. 2009; K. Batygin & A. Morbidelli 2020). Since oligarchic growth requires larger impactors, the collision energy is buried deeper in the moon. J. Monteux et al. (2014) investigated Callisto's accretional heating and found that a substantial population of impactors larger than 1 km within the Jovian CPD would be incompatible with Callisto remaining undifferentiated. While pebble accretion seems plausible, there is evidence for the presence of large particles in the CPD: kilometer-sized satellites exist both within Io's orbit and beyond Callisto's orbit in the current Jovian system. These bodies may represent collisional debris from the past evolution of the CPD, or

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Aix- Marseille Université, CNRS, CNES, Institut Origines, LAM, Marseille, France 38 Rue Frédéric Joliot Curie, 13013 Marseille, France
²Institut Universitaire de France (IUF), France</footnote>

<footnote n="3">³IRAP, University of Toulouse, CNRS, Toulouse, France</footnote>

<footnote n="4">⁴ Corresponding author.</footnote>
planetesimals trapped by the CPD that have undergone partial ablation (T. Ronnet & A. Johansen 2020).

Given the presence of both small and large particles in the Jovian CPD, a key question arises: could a substantial population of kilometer-sized impactors during accretion be consistent with a possibly undifferentiated state of Callisto, while still allowing for the full differentiation of Ganymede? This question is the focus of our study.

In this analysis, we present a thermal evolution model designed to identify the specific accretion conditions that could result in an unmelted Callisto while allowing for a fully differentiated Ganymede. Our model simulates the growth of a satellite embryo within the CPD, accounting for impacts from diverse particle populations (cisplanetary impactors) originating within the disk. Although high-velocity heliocentric particles are often cited as the primary contributors to the Galilean satellites' cratered surfaces (K. Zahnle et al. 1998), the energy flux from cisplanetary impactors is more influential during the rapid accretion phase owing to their higher collision frequency (S. W. Squyres et al. 1988). In addition to impact heating, our model includes radiogenic heating from short-lived isotopes, thermal energy from the surrounding Jovian CPD, and tidal dissipation. Key variables in our analysis include the onset time and duration of accretion, as well as the size distribution of impactors.

Section 2 outlines the thermal evolution model and the transport equation, incorporating heat contributions from radionuclides, impacts, tidal dissipation, and interactions with the surrounding CPD. Section 3 describes the criteria assumed for moon differentiation and provides estimates of the characteristic timescales for ice-rock separation and transport within their interiors. In Section 4, we present the results of a nominal case, followed by an analysis of how varying parameters influence the outcomes. Section 5 is dedicated to interpreting our findings and discussing their implications for the primordial dichotomy between Callisto and Ganymede.

## 2. Model

In this section, we provide an overview of the thermal evolution model used in our analysis, highlighting the modifications and enhancements made to address the specific objectives of our research.

### 2.1. Thermal Evolution Model

We use the one-dimensional (1D) model developed by M. Monnereau et al. (2013) to simulate the thermal evolution of a growing moon, incorporating radiogenic heating as a key heat source. This heating is essential during the early history of satellites and can potentially drive global melting within their interiors. Specifically, short-lived radionuclides such as $^{26}$Al and $^{60}$Fe provide nearly all the radiogenic energy during accretion, while the contribution from long-lived radionuclides is negligible during this phase. We opt to neglect $^{60}$Fe, as its initial $^{60}$Fe/$^{56}$Fe ratio is two orders of magnitude lower than the initial $^{26}$Al/$^{27}$Al ratio (R. K. Mishra & M. Chaussidon 2014), resulting in a significantly lower initial heating rate. The time-dependent heating rate $Q(t)_{\text{rad}}$ of $^{26}$Al is given by

$$Q(t)_{\text{rad}} = m_r \rho X_{\text{Al}} q_0 \exp[-\lambda_{26}(t + t_{\text{start}})], \quad (1)$$

where $\rho$ is the uniform density of the body, $m_r$ is the rock mass fraction, and $\lambda_{26} = 9.68 \times 10^{-7} \text{yr}^{-1}$ is the decay constant of $^{26}$Al. The initial heating rate per kg of Al is

$q_0 = 2.06 \times 10^{-5} \text{ W kg}^{-1}$ using an initial $^{26}$Al/$^{27}$Al ratio of $5.85 \times 10^{-5}$ (K. Thrane et al. 2006), which is relevant for carbonaceous chondrites.

The heating rate value $q_0$ is consistent with those derived in the studies by A. Barr & R. Canup (2008) and K. T. Trinh et al. (2023), although slight variations are found in the literature. For instance, J. Castillo-Rogez et al. (2007) adopted a value of $q_0 = 7.3 \times 10^{-6} \text{ W kg}^{-1}$. The mass fraction of Al in CI chondrites is denoted by $X_{\text{Al}}$ with CI abundances taken from K. Lodders (2003). The newly accreted material layers produce heating associated with their formation time $t + t_{\text{start}}$, where $t_{\text{start}}$ is the time when satellite accretion begins, following a chosen time span after the formation of calcium-aluminum inclusions (CAIs) in the protosolar nebula (PSN).

During accretion, the heat received by protosatellites is primarily transported by conduction within the growing moon. Protosatellites acquire heat primarily from their surface owing to accretion energy, resulting in a temperature gradient that opposes the conditions required for convection. While radiogenic heating in the satellite's interior can drive core convection, this process is unlikely to occur on a timescale shorter than the satellite's accretion period (G. Schubert et al. 1981). The onset time for convection in Ganymede and Callisto has been estimated to be on the order of $\sim 10^8 \text{ yr}$ (A. Barr & R. Canup 2008), based on the melting point viscosity of Ice V and Ice VI ($\eta \sim 10^{17} \text{ Pa s}$). Since its timescale, $\tau_{\text{conv}}$, is longer than the accretion timescale, $\tau_{\text{acc}} \sim 10^6 \text{ yr}$ (G. Schubert et al. 1981), we do not expect subsolidus ice convection to significantly affect the accretion temperature profiles.

The thermal evolution of the growing satellite is described by a 1D heat diffusion equation in spherical geometry with a time-dependent radius (M. Monnereau et al. 2013):

$$\rho c_p \left( \frac{\partial T}{\partial t} - \frac{r' \dot{R}_s}{R_s} \frac{\partial T}{\partial r'} \right) = \frac{1}{(r'R_s)^2} \frac{\partial}{\partial r'} \left( k r'^2 \frac{\partial T}{\partial r'} \right) + Q_{\text{rad}}, \quad (2)$$

where $R_s$ is the radius of the growing moon, $\dot{R}_s$ is its time derivative (the growth rate), and $r'$ is the nondimensional radial coordinate. The mean specific heat is denoted by $c_p$, and $k$ represents the thermal conductivity described in Section 2.3. The heat diffusion equation is numerically solved using the implicit Euler method on a grid with 10,000 radial meshes, where the spatial step size $\Delta r$ increases over time to account for the satellite's growth. The time step in our simulations is adaptive, determined by the characteristic timescale of impact dynamics.

### 2.2. Circum-Jovian Disk

In this study we investigate the interaction between accretion heating and the temporal evolution of the CPD temperature in order to understand how cooling of the disk during the later stages of accretion affects satellite formation. The CPD model used in this study is described in A. Schneeberger & O. Mousis (2025) and builds on the work of A. B. Makalkin & V. A. Dorofeeva (1995, 2014) and R. Heller et al. (2015). This model assumes that the CPD actively accretes material from the PSN with a protosolar metallicity $Z/H = 2.45 \times 10^{-2}$ (K. Lodders 2019) and forms shortly after the end of Jupiter's gas runaway accretion phase. Over time,
the CPD undergoes significant transformations, evolving from a massive, hot disk to a late-stage, depleted and cold state. Early in its evolution and close to the planet, the CPD is optically thick and vertically adiabatic. However, at larger distances from the planet and in later evolutionary stages, the CPD becomes optically thin and vertically isothermal. Since the model assumes that the CPD forms immediately after Jupiter's formation, it incorporates the planet's early high surface temperature, which can reach up to ~2000 K owing to accretion energy (J. Szulágyi et al. 2016). As a result, Jupiter's radiation heats the surface of the disk when it is optically thick and the midplane when it becomes optically thin. The model also takes into account the geometry of the CPD, which can cast shadows on itself, creating cold regions. In addition to Jupiter's radiative heating, the model takes into account viscous stress heating, accretion heating, and the temperature of the surrounding PSN.

The surface temperature (or photosurface temperature) of the CPD is defined by

$$T_s^4(r) = \frac{1 + (2\kappa_p(r, z_s)\Sigma_g(r))^{-1}}{\sigma_{SB}} \times (F_{vis} + F_{acc} + k_s F_p) + T_{neb}^4, \quad (3)$$

where $\Sigma_g(r)$ represents the gas surface density, $\sigma_{SB}$ the Stefan-Boltzmann constant, and $z_s$ the photospheric height. $T_{neb}$, set to 40 K (A. Schneeberger & O. Mousis 2025), denotes the temperature of the surrounding PSN, and $\kappa_p$ is the radius-dependent opacity. $F_{vis}$, $F_{acc}$, and $F_p$ are the heat fluxes produced by viscous stress, accretion of material onto the CPD, and Jupiter's radiative heating, respectively. The factor $k_s$ is the fraction of radiation absorbed by the disk, the rest being reflected and scattered away. The midplane effective temperature, $T_e$, is derived from the surface temperature, $T_s$, under the assumption that the inner region of the CPD is optically thick and adiabatic. Accordingly, $T_e$ is determined using a gray atmosphere radiative model (A. B. Makalkin & V. A. Dorofeeva 2014; A. Schneeberger & O. Mousis 2025). Figure 1 shows the evolution of the CPD midplane temperature $T_e(t)$ at the present-day positions of Ganymede and Callisto, situated at 15 and 26.3 Jupiter radii, respectively. Accretion of the moons is assumed to begin when the snowline crosses their orbits, which are fixed in our study, typically occurring several tens of thousands of years after the CPD formation for Ganymede and Callisto. Since Ganymede would form in a more heated environment compared to Callisto, the condensation front migrates from Callisto's location to Ganymede's in only a few thousand years, a time span significantly shorter than the timescale of the accretion phase, $\tau_{acc}$. The resolution method and detailed structure of the CPD are described in A. Schneeberger & O. Mousis (2025, and references therein).

## 2.3. Growing Conditions

The accretion conditions of the Galilean moons are primarily determined by three factors: the lifetime of the CPD, the mass inflow rate, and the size of the bodies orbiting around the gas giant. The growth of the satellites occurs as they accrete particles from the CPD surrounding Jupiter, with a

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Evolution of the midplane effective temperature $T_e(t)$ in the Jovian CPD at the current orbits of Ganymede (solid line) and Callisto (dashed line). The CPD model assumes a metallicity $Z/H = 2.45 \times 10^{-2}$, a disk viscosity parameter of $10^{-3}$, a centrifugal radius of $50R_J$, an ambient temperature $T_{neb} = 40$ K, and a planetary temperature $T_p = 2000$ K.</figcaption></figure>

mass inflow rate depicted by the following relation:

$$\dot{M}_{Sat} = \frac{M_{Sat}}{\tau_{acc}}, \quad (4)$$

where $M_{Sat}$ is the mass of the satellite and $\tau_{acc}$ is the duration of the accretion phase. As accretion timescales can differ widely between formation scenarios, we treat $\tau_{acc}$ as a free parameter to explore a broad range of accretion conditions. This parameter ranges from 10 kyr, which aligns with the solids-enhanced minimum mass model (I. Mosqueira & P. R. Estrada 2003a, 2003b), to a few Myr, a value suitable for the longer accretion timescales derived from the classic gas-starved model (R. Canup & W. Ward 2002).

We assume that Ganymede and Callisto accreted from the same source of pristine material within the Jovian CPD. The particles impacting the embryos have a radius $r_{imp}$ and the same mean density ($\rho \sim 1830 \text{ kg m}^{-3}$) as Callisto, consisting of a uniform mixture of ice and rock, with a rock mass fraction $m_r = 0.44$ (W. B. McKinnon 1997).

The thermal conductivity of the mixture is estimated as the mass-weighted harmonic average of the rock thermal conductivity, $k^{rock} = 4 \text{ W m}^{-1} \text{ K}^{-1}$ (K. Yomogida & T. Matsui 1983), and the water-ice thermal conductivity, $k^{ice} = 2.2 \text{ W m}^{-1} \text{ K}^{-1}$ (W. F. Waite et al. 2006), with $k = (m_r/k^{rock} + (1 - m_r)/k^{ice})^{-1} \sim 2.8 \text{ W m}^{-1} \text{ K}^{-1}$. It is worth noting that the thermal conductivity of ice, $k^{ice}$, is inversely proportional to temperature, following the relation $k^{ice} \approx 600/T$ with $T$ in K and $k^{ice}$ in $\text{W m}^{-1} \text{ K}^{-1}$ (J. Klinger 1980; E. M. Morris 2001; T. B. McCord & C. Sotin 2005). However, this temperature dependence has a minimal impact in our study owing to the relatively narrow temperature range considered (200–273 K) and the presence of rocks ($m_r = 0.44$), which reduce the sensitivity of the effective thermal conductivity to temperature. The mean specific heat, $c_p \sim 1480 \text{ J kg}^{-1} \text{ K}^{-1}$, is also inferred by
<table><caption>**Table 1**
Impactor Populations</caption><thead><tr><th>Population</th><th>Size Range</th><th>Energy Deposit</th></tr></thead><tbody><tr><td>Small</td><td>1 cm–100 m</td><td>Surface</td></tr><tr><td>Medium</td><td>100 m–1 km</td><td>Surface and subsurface</td></tr><tr><td>Large</td><td>1–100 km</td><td>Deep inside</td></tr></tbody></table>

the mass-weighted average, with $c_p^{ice} = 2100$ J kg$^{-1}$ K$^{-1}$ and $c_p^{rock} = 700$ J kg$^{-1}$ K$^{-1}$ (A. Barr & R. Canup 2008).

The size distribution of impactors is given by the power law

$$ \frac{dN}{dr_{\text{imp}}} \propto r_{\text{imp}}^{-\alpha}, \qquad (5) $$

with the slope $\alpha$ ranging from 1 to 6 (S. W. Squyres et al. 1988). The radius of cisplanetary impactors, $r_{\text{imp}}$, can vary from a few centimeters to several hundred kilometers, depending on the disk model and the formation mechanism considered (P. R. Estrada et al. 2009; T. Ronnet et al. 2017; K. Batygin & A. Morbidelli 2020). To account for the different sizes of bodies during accretion, we have considered three main particle populations, with each population contributing to the growth of the moon's radius $R_s(t)$ (see Table 1 for the considered size range of impactors).

Small impactors are assumed to be large enough to deposit energy on the surface but too small to penetrate significantly, leading to energy dissipation predominantly through radiation. Medium-sized impactors are large enough to deposit a fraction of their energy into the subsurface within a layer of thickness $\delta$, which exceeds the depth $\delta_{\text{cool}}$, subject to conductive cooling before accumulating a new layer (A. Barr & R. Canup 2008). To enhance computational efficiency, continuous accretion is assumed for particles with radii below 1 km. The accumulated mass from both medium and small particle populations is averaged and uniformly distributed across the surface using successive deposition layers. In contrast, larger impactors are assumed to generate significant shock waves that can penetrate deeply into the icy-rocky structure of protosatellites (J. Monteux et al. 2014). These impactors are treated explicitly, with their radius and impact frequency randomly determined.

The mass fraction $x_m^i$ of the ith population is obtained by the ratio between the mass of the ith population $M^i$ and the total mass. The mass $M^i$ is derived by integrating the number of particles (Equation (5)) multiplied by volume and density, yielding

$$ x_m^i = \frac{\int_{r_{\text{min}}}^{r_{\text{max}}} r_{\text{imp}}^{3-\alpha} dr_{\text{imp}}}{\int_{r_{\text{min}}}^{r_{\text{max}}} r_{\text{imp}}^{3-\alpha} dr_{\text{imp}}}, \quad \text{with } x_m^{\text{si}} + x_m^{\text{mi}} + x_m^{\text{li}} = 1, \quad (6) $$

where $r_{\text{min}}^i$ and $r_{\text{max}}^i$ represent the size range of the ith population (see Table 1) and superscripts si, mi, and li denote the populations of small, medium-sized, and large impactors, respectively. Large impactors are treated individually, with impact probabilities determined using a Monte Carlo approach. For each incoming particle, the radius of the impactor is randomly generated based on the probability density function $f(r_{\text{imp}}) = Cr_{\text{imp}}^{-\alpha}$, with C a constant. By reversing the cumulative distribution function $F_R(r)$, the radius

of the next impactor is given by

$$ \begin{gathered} r_{\text{imp}} = \left( \frac{F_R(r)(1-\alpha)}{C} + r_{\text{min}}^{(1-\alpha)} \right)^{\frac{1}{1-\alpha}}, \\ \text{with } C = \frac{1-\alpha}{r_{\text{max}}^{1-\alpha} - r_{\text{min}}^{1-\alpha}}, \qquad (7) \end{gathered} $$

where $F_R(r)$ is a random number drawn between 0 and 1. The latitude $\theta$ and longitude $\phi$ of each impact are also randomly generated to simulate random impact frequencies. We assume here an isotropic impact flux, as inclinations and eccentricities of large particles are randomized by frequent collisions, preventing any preferred impact direction. The latitude and longitude are given by

$$ \left\{ \begin{aligned} \theta &= \frac{180}{\pi} \arcsin(1 - 2r_1) \\ \phi &= 4 \times \frac{180}{\pi} |\arcsin(1 - 2r_2)|, \end{aligned} \right. \qquad (8) $$

where $r_1$ and $r_2$ are two independent random numbers between 0 and 1. An impact occurs when the incoming particle crosses a reference impact location. Neighboring impacts are considered by taking a margin of one impact diameter, i.e., the outer edges of the particle must at least touch the impact site, which is defined by a set of latitude–longitude coordinates. This 1D approach allows us to examine a random point on the protosatellite and study its internal structure.

When a large particle impacts the surface, a shock wave propagates from the impact site toward the center of the satellite. The initial kinetic energy of the impactor is converted into shock heating, internal energy produced by plastic work, and the kinetic energy of the ejecta (S. W. Squyres et al. 1988). We do not include the ejecta in the energy balance because the fraction of energy converted to kinetic energy is only between 7% and 10% (J. D. Okeefe & T. J. Ahrens 1977). We also neglect the post-shock impact heating caused by plastic deformation, as shock physics simulations have shown that the heat is primarily localized at the surface (K. Kurosawa et al. 2021).

When a shock wave propagates through a protosatellite, the peak shock pressure is concentrated within the isobaric core, a spherical region where the pressure is approximately uniform (S. K. Croft 1982). The temperature increase $\Delta T$ within the isobaric core due to each impact is given by the following relation (J. Monteux et al. 2007, 2014):

$$ \Delta T = \frac{4\pi \gamma_{li} \rho G R_s^2}{9 h_m c_p}, \qquad (9) $$

where $R_s$ is the growing satellite radius and $G$ is the gravitational constant. The term $h_m = \frac{V_{\text{eff}}}{V_{\text{ic}}} \sim 5.8$ (J. Monteux et al. 2014) represents the volume effectively heated normalized by the volume of the isobaric core, with $V_{\text{ic}} = 3V_{\text{imp}}$, where $V_{\text{imp}}$ is the volume of the impactor (H. Senshu et al. 2002). The parameter $\gamma_{li}$ denotes the fraction of the impactor's kinetic energy converted into shock heating, which ranges from 0.2 to 0.6 (J. D. Okeefe & T. J. Ahrens 1977) for low-velocity ($v_{\text{imp}} \lesssim 5$ km s$^{-1}$) and high-velocity ($v_{\text{imp}} > 10$ km s$^{-1}$) impacts, respectively. Since this study only considers cisplanetary particles, we neglect the impactor velocity at infinity, resulting
in $v_{\text{imp}} \approx v_{\text{esc}}$, where $v_{\text{esc}}$ is the escape velocity. For Ganymede and Callisto, where $v_{\text{imp}} < 3 \text{ km s}^{-1}$, we adopt $\gamma_{li} = 0.2$.

Outside the isobaric core, both the shock pressure and the thermal perturbations decay significantly via the following relation (J. Monteux et al. 2007):

$$T(\bar{r}) = \begin{cases} T + \Delta T & \bar{r} \le r_{\text{ic}} \\ T + \Delta T \left(\frac{r_{\text{ic}}}{\bar{r}}\right)^m & \bar{r} > r_{\text{ic}} \end{cases} \quad (10)$$

where $\bar{r}$ is the distance from the center of the isobaric core with radius $r_{\text{ic}}$ and $m$ is an exponent ranging from 2.8 to 4. In this study, we adopt the median value $m = 3.4$ (J. Monteux et al. 2014).

While large impactors are treated individually, small and medium impactors are uniformly added to the surface with an average mass inflow rate of $\dot{m}_{\text{imp}}^i = x_m^i \dot{M}_{\text{Sat}}$. At each time step, a thin uniform layer is deposited on the surface with a thickness given by

$$\delta^i = \left( \frac{3\dot{m}_{\text{imp}}^i \Delta t}{\rho 4\pi} + R_s^3 \right)^{1/3} - R_s. \quad (11)$$

Medium impactors heat the surface but also the subsurface, depositing a fraction of their energy uniformly to a depth of $\delta^{\text{mi}}$. Assuming $v_{\text{imp}} \approx v_{\text{esc}}$, the temperature profile of the additional layer is given by the relation (G. Schubert et al. 1981; J. Monteux et al. 2014)

$$T(R_s) = \frac{\gamma_{\text{mi}} GM_s}{c_p R_s} + T_e, \quad (12)$$

where $M_s$ is the mass of the growing moon and $T_e$ is the time-dependent temperature of the disk around the gas giant. The fraction of kinetic energy converted into internal energy from medium impactors is $\gamma_{\text{mi}} = 0.1$. The parameters $\gamma_{\text{mi}}$ and $\gamma_{li}$ differ from the parameter $h$ introduced by G. Schubert et al. (1981), which includes both the energy that goes into heating and the post-impact surface cooling (see discussion in J. Monteux et al. 2014).

The surface is continuously heated by small and medium impactors originating from the background disk around Jupiter. We assume that these impactors are at the same temperature as the CPD gas and are small enough to be approximated by continuous accretion. Besides the gravitational energy deposited during the impact, the surface is also heated or cooled, depending on whether these bodies are hotter or colder than the surface. The time-dependent surface temperature $T_{\text{surf}}$, used as the outer boundary condition for the heat diffusion equation, is determined by the balance between the gravitational energy deposited by impactors, radiation to space, and heat conducted to the surface or toward the center of the growing body:

$$\sigma (T_{\text{surf}}^4 - T_e^4) 4\pi R_s^2 = k \frac{\partial T}{\partial r} 4\pi R_s^2 + \frac{GM_s \dot{M}_{\text{sat}}}{R_s} (x_m^{\text{si}} + x_m^{\text{mi}} (1 - \gamma_{\text{mi}})) - c_p \dot{M}_{\text{sat}} (x_m^{\text{mi}} + x_m^{\text{si}}) (T_{\text{surf}} - T_e). \quad (13)$$

Here $\sigma$ is the Stefan-Boltzmann constant and $(1-\gamma_{\text{mi}})$ represents the fraction of the initial kinetic energy of medium impactors that is not buried inside the satellite but is instead

subject to radiative cooling. The left-hand side of the equation represents the difference between the radiation from the surface and the background disk. From left to right, the other terms of the equation describe heat conduction, gravitational potential energy delivered by small and medium impactors, and the heating (or cooling) induced by these particles. Larger impactors also heat the surface, but this energy is efficiently radiated away. After an elevation of temperature $\delta T_0$ caused by an impact, the surface temperature difference decreases over time according to $\delta T(t)/\delta T_0 \approx \sqrt{\kappa} \rho c_p / 4\sigma T_e^3 \sqrt{\pi t}$ (H. S. Carslaw & J. C. Jaeger 1959) for times longer than the cooling timescale $t_{\text{cool}} = 10\kappa (\rho c_p / 4\sigma T_e^3)^2$, where $\kappa \sim 10^{-6} \text{m}^2 \text{s}^{-1}$, $T_e = 200 \text{K}$, $\rho \sim 1830 \text{kg m}^{-3}$, and $c_p \sim 1480 \text{J kg}^{-1} \text{K}^{-1}$. After approximately $t_{\text{cool}} \sim 250$ days, the temperature difference drops to less than 20% of the initial temperature increase $\delta T_0$, and after 10 yr, our average time step $dt$, it falls below 5%. Since our study does not focus on the detailed short-term evolution of surface temperature, the fraction of the initial kinetic energy of large impactors lost through radiative cooling can be neglected.

As protosatellites grow, tidal forces increasingly dissipate energy within them, contributing significantly to their heat budget. The tidal dissipation in a homogeneous body in synchronous rotation with small eccentricity $e$ follows the prescription (J. Wisdom 2004)

$$\frac{dE}{dt} = \frac{21 k_2 GM_{\text{jup}}^2 n R_s^5 e^2}{2 Q a^6}, \quad (14)$$

where $k_2$ is the Love number, $Q$ is the tidal dissipation parameter, $n$ is the orbital mean motion, and $a$ is the semimajor axis. The factor $k_2/Q$ represents the deformability of the satellite.

For Ganymede and Callisto, we assume $k_2/Q \approx 0.05$ (R. Malhotra 1991; B. Downey et al. 2020) and adopt constant orbital parameters $(a, e)$ during accretion, assuming that the satellites formed in situ. Note that an additional expression for tidal dissipation in a body undergoing nonsynchronous rotation has been provided by J. Wisdom (2008), which shows that the dissipation in the case of synchronous rotation is greater than or equal to that in the case of asymptotic nonsynchronous rotation. The use of Equation (14) is justified by the fact that radiogenic heating dominates by several orders of magnitude during the accretion phase, as shown below, and because we assume in situ formation of the moons. However, due to the significant uncertainties associated with tidal dissipation, spin states, and orbital dynamics during this period, it is important to emphasize that Equation (14) represents an end-member scenario of minimal tidal heating. In this study, we have adopted several simplifying assumptions, including synchronous rotation, the neglect of dissipation in fluid layers, and the use of a constant dissipation factor that does not account for variations in temperature, composition, or pressure.

### 3. Melting and Differentiation

To determine whether melting occurs, we compare the temperature profiles within the newly formed satellite to the pressure-dependent melting temperature of water ice. To do so, we account for the different phases of ice as the pressure $P(r, t)$
<table><caption>**Table 2**
Ice Phases and Corresponding Coefficients (J. Leliwa-Kopystyński et al. 2002)</caption><thead><tr><th>Ice Phases</th><th>Pressure Range (MPa)</th><th>T<sub>0</sub> (K)</th><th>A (K Pa<sup>-1</sup>)</th><th>B (K Pa<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Ice I</td><td>0 &lt; P &le; 209</td><td>273.2</td><td>-7.95 &times; 10<sup>-8</sup></td><td>-9.6 &times; 10<sup>-17</sup></td></tr><tr><td>Ice III</td><td>209 &lt; P &lt; 344</td><td>247.7</td><td>2.38 &times; 10<sup>-8</sup></td><td>0</td></tr><tr><td>Ice V</td><td>344 &lt; P &lt; 626</td><td>242.5</td><td>4.9 &times; 10<sup>-8</sup></td><td>0</td></tr><tr><td>Ice VI</td><td>626 &lt; P &lt; 2150</td><td>190.3</td><td>1.54 &times; 10<sup>-7</sup></td><td>-3.43 &times; 10<sup>-17</sup></td></tr><tr><td>Ice VII</td><td>P &gt; 2240</td><td>149.8</td><td>9.14 &times; 10<sup>-8</sup></td><td>0</td></tr></tbody></table>

increases over time in a given layer of radius *r*:

$$P(r, t) = \frac{2\pi}{3} G \rho^2 (R_s^2(t) - r^2). \qquad (15)$$

We assume a uniform density and zero pressure at the surface.

The time-dependent liquidus temperature for a solid mixture of ammonia and water, including high-pressure phases, is determined using the formulation provided by J. Leliwa-Kopystyński et al. (2002):

$$T_m(P(r, t), X_{3\text{NH}}) \approx T_0 + AP + BP^2 - CX_{3\text{NH}} - DX_{3\text{NH}}^2, \qquad (16)$$

where $T_0$, $A$, and $B$ are the coefficients for the different phases of pure water ice given in Table 2 and taken from J. Leliwa-Kopystyński et al. (2002). $C$ and $D$ are coefficients with values of 53.8 and 650, respectively. Since other volatiles could condense in the CPD, ammonia and other salts could be present in the ice mixture, lowering the melting point of water and allowing a subsurface ocean to remain stable beneath a convecting ice shell (T. Spohn & G. Schubert 2003). However, the presence of ammonia raises the question of whether its inclusion is consistent with the formation of Callisto without significant melting. To investigate this effect, specifically the effect of NH₃ on an undifferentiated Callisto, we have included this species in the ice mixture, where $X_{NH3}$ represents the mass fraction. Since the phase diagram of NH₃ at high pressure remains poorly constrained (O. Mousis et al. 2002), we consider ammonia only at the liquidus temperature in the outer layers of the protosatellite.

If melting occurs during accretion, rock particles can settle to the bottom of the resulting localized pool of liquid water, leading to local differentiation. Specifically, an impact that melts the ice produces a transient pool of liquid water in which the denser rock particles may separate and settle at the base, provided that they can sink before the water refreezes. Modeling the rock particles as spheres of radius $r_{rock}$, their settling velocity can be estimated using Stokes's law. This yields a sinking timescale $t_{sink} = 9\eta_l d_{pool}/2g(\rho_r - \rho_l)r_{rock}^2$ (A. C. Barr & R. M. Canup 2010b), where $d_{pool}$ is the depth of the pool, $\eta_l = 10^{-3}$ Pa s is viscosity of liquid water, $\rho_r = 3000$ kg m⁻³ is the rock density, $\rho_l = 1000$ kg m⁻³ is the liquid water density, and $g = 1.24$ m s⁻² is the gravity on Callisto. For an impactor of radius $r_{imp} = 1$ km, which melts a volume corresponding to $V_{ic} = 3V_{imp}$ (H. Senshu et al. 2002), where $V_{ic}$ is the volume of the isobaric core, the pool depth is given by $d_{pool} = r_{ic}$, where $r_{ic}$ is the radius of the isobaric core. For small-grain particles $r_{rock} = 1$ cm, the estimated timescale is $t_{sink} \sim 30$ s. This must be compared to the timescale for the pool to freeze $t_{cryst} \sim L\rho_l r_{ic}^2/k^{\text{ice}}\Delta T$, where $L = 3 \times 10^5$ J kg⁻¹ is

the latent heat of water ice and $\Delta T$ is the temperature gradient between the cold surface and the depth of the pool. For $\Delta T = 150$ K, this yields $t_{cryst} \sim 2 \times 10^4$ yr, consistent with published estimates of freezing times ($10^3-10^5$ yr; C. Michaut & M. Manga 2014; M. A. Hesse & J. C. Castillo-Rogez 2019; C. J. Chivers et al. 2021; M. Naseem et al. 2023). Given that $t_{sink} \ll t_{cryst}$, even the smallest accreted and co-accreted particles would likely settle to the bottom of the melt pool well before it solidifies.

Once separated from the ice and settled at the base of the melt pool, the layer of dense rock particles becomes gravitationally unstable and may undergo convective overturn driven by Rayleigh-Taylor (RT) instabilities (M. E. Rubin et al. 2014; D. Shoji & K. Kurita 2014). These instabilities arise when the dense, rock-rich layer overlies a less dense, undifferentiated region, or when an undifferentiated layer is subsequently deposited atop a refrozen, rock-depleted melt pool. In either case, the instability facilitates the continued descent of rock material toward the interior, contributing to partial differentiation of the moon. Following the studies of M. E. Rubin et al. (2014) and D. Shoji & K. Kurita (2014), the timescale for the development of RT instabilities can be estimated by

$$\tau_{\text{RT}} = \eta^{\text{mix}} \left( (n-1)^{(1/n)} \frac{C_{L\Delta\rho}}{2n} \Delta\rho g L_{\text{RT}} \right)^{-1} \left( \frac{Z_0}{L_{\text{RT}}} \right)^{(1-n)/n}, \qquad (17)$$

where $\eta^{\text{mix}}$ is the viscosity of the ice-rock mantle, $C_{L\Delta\rho} \approx 0.76$ is a dimensionless constant that depends on the rheology and geometry of the system, and $n$ is the stress exponent, for which we adopt $n = 1.8$, assuming that grain boundary sliding dominates the creep behavior of the ice (M. E. Rubin et al. 2014). $Z_0$ represents the initial perturbation amplitude and $L_{RT} = nRT_0^2/E_a|dT/dz|$ the length scale over which the viscosity varies significantly, with $T_0$ the temperature at the interface, $E_a$ an activation energy, and $R = 8.314$ J mol⁻¹ K⁻¹ the gas constant. We assume $T_0 = 200$ K, $E_a = 49 \times 10^3$ J mol⁻¹, $|dT/dz| = 1.8$ K km⁻¹, and $Z_0 = 1$ km (M. E. Rubin et al. 2014). The viscosity of the ice-rock mantle is modeled as $\eta^{\text{mix}} = \eta^{\text{ice}}/f(\phi)$ (K. Nagel et al. 2004), where $\eta^{\text{ice}}$ is the ice viscosity, which depends exponentially on temperature, and $f(\phi)$ accounts for the influence of the rock fraction. For simplicity, we adopt an Arrhenius-type relation for the viscosity of ice, $\eta^{\text{ice}} = \eta^{\text{ref}}\exp(A(T_{\text{ref}}/T - 1)$, with $\eta^{\text{ref}} = 10^{14}$ Pa s, $T_{\text{ref}} = 273$ K, and $A = 20$ (D. Shoji & K. Kurita 2014). The influence of the rock volume fraction is given by $f(\phi) = (1 - \phi_0/\phi_{\text{CPL}})^3$, where $\phi_0 = m_r\rho/\rho_r = 0.27$ is the assumed volume fraction of rock, $\phi_{\text{CPL}} = 0.74$ is the close-packing limit, and $\beta$ is an empirical exponent ranging from 1.55 to 2.5 (K. Nagel et al. 2004). For our study, we adopt the median value, $\beta = 2.0$. Figure 2 illustrates the RT overturn timescale $\tau_{\text{RT}}$ as a function of typical mantle temperatures (150–273 K). The figure clearly illustrates the strong temperature dependence of the instability timescale, which spans from approximately 1 Gyr to just a few hundred years across the range of mantle temperatures examined in this study. Our analyses suggest that gravitationally unstable, rock-rich layers are likely to overturn on timescales comparable to the duration of satellite accretion, particularly when mantle temperatures exceed $\sim 180$ K, a threshold reached in all simulations.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Stokes sinking and overturn timescales vs. typical mantle temperature (150–273 K). Both timescales are inferred using the viscosity of the ice-rock mixture $\eta^{\text{mix}} = \eta^{\text{cc}}/f(\phi)$ and an Arrhenius-type relation for ice viscosity $\eta^{\text{ice}} = \eta^{\text{ref}}\exp(A(T_{\text{ref}}/T - 1))$. Timescales for Stokes sinking are displayed for different sizes of rock particles (100 m–1 km).</figcaption></figure>

Differentiation may also be facilitated by Stokes sinking, where ice-rock separation occurs via solid-state creep of ice around the denser rock. A rock fragment could cross the satellite by Stokes flow over a timescale given by

$$\tau_{\text{Stokes}} = \frac{9}{2} \frac{R_s \eta^{\text{mix}}}{(\rho_r - \rho_{\text{ice}}) g r_{\text{rock}}^2}, \qquad (18)$$

where $\rho_{\text{ice}} = 1400 \text{ kg m}^{-3}$ represents the compressed density of various high-pressure ice phases (A. Barr & R. Canup 2008). Figure 2 also presents the Stokes sinking timescale $\tau_{\text{Stokes}}$ as a function of typical mantle temperatures (150–273 K) and rock fragment sizes (100 m–1 km). This suggests that kilometer-scale particles (the typical size of the largest particles generated in our simulations) are unlikely to reach the base of the moon within accretion timescales unless the mantle approaches the melting temperature of water ice. Smaller particles (1 cm–100 m) are even less likely to settle before complete melting of the ice mantle occurs. However, it is important to note that Stokes sinking can also occur after rock has already separated from the ice within a melt pool. In this context, the size of a sinking fragment does not necessarily reflect that of the initially accreted particles. Once at the base of the melt pool, individual rock grains may accumulate and coalesce into larger fragments before the pool fully refreezes (A. C. Barr et al. 2010; A. C. Barr & R. M. Canup 2010b; W. B. Tonks et al. 2016). This process could lead to the formation of significantly larger rock bodies than those initially accreted, thereby substantially reducing their Stokes sinking timescale. However, because this mechanism is highly sensitive to the melt volume and remains poorly constrained, it is not quantified and is excluded from the scope of this study.

Overall, our calculations indicate that efficient rock-ice separation likely begins once the liquidus is reached and

melting initiates. Following decoupling from the ice, rock particles can migrate inward and contribute to core formation on timescales comparable to that of satellite accretion, primarily through convective overturn of gravitationally unstable layers. The gravitational energy released during this downward migration may significantly influence the satellite's thermal evolution, as it is substantial—approximately ~3 × 10⁵ J kg⁻¹, nearly two orders of magnitude greater than the energy required to melt all the ice, estimated at (1 − $m_r$)ML ~5 × 10³ J kg⁻¹. In this context, we assume that global melting and subsequent differentiation occur once the water liquidus is exceeded, an assumption also adopted by J. Monteux et al. (2014) and A. Barr & R. Canup (2008).

## 4. Results

Figure 3 shows an example temperature profile inside a satellite similar to Callisto and Ganymede, referred to as our nominal model, at the end of the accretion phase (left panel) and over time (right panel), including radiogenic and tidal heating. The nominal model is characterized by the critical accretion parameters $t_{\text{start}}$, $\tau_{\text{acc}}$, and $\alpha$, which are the onset time of accretion, the duration of the accretion phase, and the slope of the impactor size distribution, respectively. In our nominal model, we assume $\alpha = 4$, leading to a large impactor contribution of about $x_m^{\text{li}} = 30\%$ of total mass. The heat deposited deep within the growing moon by these large impactors is slow to dissipate, as evidenced by the temperature spikes shown in the left panel of Figure 3. Since conduction is a relatively slow heat transport mechanism compared to accretion, it is only effective on timescales of the order of ~10⁸ yr.

As the radius of the growing moon ($R_s$) increases, impactors accelerate ($v_{\text{imp}} \propto R_s^2$), leading to more energetic collisions. Consequently, the heat generated by these impacts intensifies toward the end of the accretion phase, potentially melting the ice in several outer layers, as shown by the high temperature rises in the right panel of Figure 3. Although impacts become more energetic during this stage, the peak temperature occurs primarily in the subsurface rather than at the surface. This is because radiative cooling, which dominates heat transport, prevents surface temperatures from rising significantly. When the temperature increases by $\delta T$, the surface rapidly cools through radiation within a few hundred days. As a result, large impactors are unable to efficiently heat the surface, and the surface temperature is instead determined by the balance of fluxes (see Equation (13)).

During accretion, localized melting can lead to the formation of transient melt pools, in which rock particles separate from ice and settle to the base, creating a dense, rock-rich layer above a less dense, undifferentiated region (see Section 3). As illustrated in Figure 3, the ice-rock mantles of both satellites reach temperatures above ~180 K, resulting in sufficiently low viscosities to enable gravitational instability. Consequently, the dense layer becomes prone to convective overturn on timescales comparable to that of satellite accretion. This overturn contributes both to the satellite's thermal evolution and to its partial differentiation. In this study, global melting and subsequent differentiation are assumed to occur once the water liquidus is exceeded. Since our simulations follow a Monte Carlo approach, the timing and frequency of impacts are randomly distributed, leading to slight run-to-run variations in the thermal evolution, even
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Simulation of the accretion of a protosatellite similar to Callisto and Ganymede, with parameters $\alpha = 4$, $\tau_{acc} = 2$ Myr, and $t_{start} = 4$ Myr. The body begins accreting material at $t_0 = t_{start}$, starting with an initial radius $R_0 = 100$ km. The simulation stops when the moon's radius reaches that of Callisto ($R_f = 2410$ km) and Ganymede ($R_f = 2630$ km) at $t_{end} = t_{start} + \tau_{acc}$. Left panel: temperature profile radius (black line) within the protosatellite at the end of the accretion phase. The red line represents the melting curve $T_m(P, X_{NH3})$ for water ice with 5% NH$_3$. The layer where the temperature profile exceeds the melting curve undergoes melting, where the percentage of material above the melting point of water ice (in volume) is 2.2% and 10.4% for Callisto and Ganymede, respectively. Right panel: time evolution of the temperature distribution during the accretion phase.</figcaption></figure>

under identical input parameters. Specifically, the volume of material exceeding the liquidus may vary by up to 5% between simulations. To account for this variability, we incorporate a margin of uncertainty into our differentiation criterion and define global melting as occurring when a small fraction—up to 5% by volume of the satellite's interior—surpasses the liquidus, consistent with the threshold adopted by J. Monteux et al. (2014). Under this criterion, the simulation shown in Figure 3, which yields only 2.2% melting by volume, suggests that Callisto, unlike Ganymede, likely avoided global melting during accretion for the considered set of parameters.

During the accretion phase, tidal forces can be neglected compared to radiogenic and accretional heating. Using Equation (14), we find that tidal dissipation is several orders of magnitude lower than radiogenic heating, with radiogenic

heating $Q_{rad} \approx 8.0 \times 10^{-6} \text{ W m}^{-3}$ at time $t_{start} = 3$ Myr after CAIs. In contrast, tidal dissipation $Q_{tid} \approx 9.3 \times 10^{-10} \text{ W m}^{-3}$ and $3.9 \times 10^{-9} \text{ W m}^{-3}$ for Callisto and Ganymede, respectively. However, during the post-accretion phase, which is beyond the scope of this study, the situation changes. Long-lived radionuclides like uranium and thorium provide heating sources that are comparable to tidal heating, making the latter a more significant contributor to the interior's warming. This additional tidal heating would complement both the heat retained from the accretion phase and the ongoing radiogenic heating.

Figure 4 illustrates the melting behavior of satellites as influenced by the key parameters $\alpha$, $\tau_{acc}$, and $t_{start}$. The top two panels, corresponding to a slope $\alpha = 3$, show that both Ganymede and Callisto undergo global ice melting when large
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Final states of Callisto and Ganymede as functions of $t_{\text{start}}$, $\tau_{\text{acc}}$, and $\alpha$. The white region indicates where Ganymede undergoes melting, while Callisto remains undifferentiated. From left to right, the panels display increasing values of $\alpha$, ranging from 3 to 5, with a total of 100 simulations per panel. Values equal to 6 are not shown, as the mass fraction distribution remains unchanged between $\alpha = 5$ and $\alpha = 6$, with ~100 wt% for small impactors ($x_m^{\text{si}}$) in both cases. To prevent melting during formation, Callisto must accrete a minimum amount of large impactors ($\alpha \gtrsim 4$), namely less than 30% of its final mass. For $\alpha = 4$ (or $\alpha \gtrsim 5$) Callisto's accretion likely occurred slowly, with $\tau_{\text{acc}} \gtrsim 2$ Myr (or $\tau_{\text{acc}} \gtrsim 0.6$ Myr), and started no earlier than $t_{\text{start}} = 3.5$ (or $t_{\text{start}} = 3$ Myr) after CAIs. Using the same parameters that prevented Callisto from melting, Ganymede could have differentiated early if $\alpha = 4$, suggesting that the observed dichotomy may be primordial.</figcaption></figure>

impactors predominantly dominate the size distribution, regardless of the accretion onset time or timescale. To prevent Callisto from melting during its accretion, the population of impactors must be skewed toward smaller bodies ($\alpha \gtrsim 4$). More precisely, the satellite must accumulate only a small fraction of large impactors, with less than $x_m^{\text{li}} \sim 30\%$ of its final mass. Otherwise, substantial energy will be deposited deep within the moon owing to shock heating, and neither prolonged accretion times nor delayed formation would be sufficient to inhibit melting.

Our calculations suggest that Callisto's accretion likely occurred slowly ($\tau_{\text{acc}} \gtrsim 2$ Myr) and concluded no earlier than 5.5 Myr for $\alpha = 4$. Unlike accretional heating, which effectively melts the ice in the outer layers of the satellites, radiogenic heating primarily affects their deep interiors, as the heating rate depends on the formation time of each layer. When considering both accretional and radiogenic heating, the thermal perturbations caused by impact heating add to the heat from radioisotopes. Consequently, if small particles dominate the impactor population, the period during which accretion occurs without melting could be shorter and begin earlier. Typically, for $\alpha \gtrsim 5$, we find that Callisto could have formed without melting over a timescale of $\tau_{\text{acc}} \gtrsim 0.6$ Myr and with $t_{\text{start}} \gtrsim 3$ Myr.

It is also noteworthy that an accretion onset time of less than $t_{\text{start}} = 3$ Myr results in global differentiation, regardless of the particle size distribution or accretion timescale. The substantial abundance of $^{26}$Al causes a significant rise in interior temperature, approaching the melting point. As the viscosity of the ice mantle decreases considerably, rock particles may

begin to sink toward the center of the protosatellite, creating conditions conducive to slow differentiation.

The main differences between the formation conditions of Ganymede and Callisto are the disk temperature and the final radius, which result in higher impact velocities for Ganymede ($v_{\text{imp}} \propto R_s^2$). Ganymede can only be preserved from melting if the power-law slope $\alpha \gtrsim 5$, indicating that only small bodies are accreted. Furthermore, longer accretion timescales or later accretion could not prevent melting if Ganymede is formed from objects larger than a few hundred meters.

Figure 4 further explores the parameter space ($\alpha, \tau_{\text{acc}}, t_{\text{start}}$) that could produce the primordial dichotomy observed in the classic gas-starved disk model (R. Canup & W. Ward 2002), where Ganymede and Callisto are thought to have formed around the same time. The middle panels show that $\alpha = 4$ is the most suitable slope to explain this dichotomy. With the same parameters $\tau_{\text{acc}}$ and $t_{\text{start}}$ that prevent Callisto from melting, we find that Ganymede could achieve early differentiation if both satellites accreted a significant number of medium- and large-sized bodies ($\alpha \sim 4$), over a timescale $\tau_{\text{acc}} \gtrsim 2$ Myr, and before $t_{\text{start}} \gtrsim 3.5$ Myr after CAIs.

## 5. Discussion

We have shown that Callisto can remain undifferentiated, despite forming through the accretion of a substantial influx of kilometer-sized impactors—up to 30% of its final mass—while simultaneously allowing the complete differentiation of Ganymede. This dichotomy between the two moons naturally emerges during their accretion phases, as both form under identical conditions governed by the same set of parameters
<table><caption>**Table 3**
Summary of Scenarios Leading to a Primordial Dichotomy</caption><thead><tr><th>Accretional Parameters<br>(τ<sub>acc</sub>, t<sub>start</sub>)</th><th>Callisto Undifferentiated</th><th>Ganymede Melts</th></tr></thead><tbody><tr><td>Same values</td><td>α ~ 4, τ<sub>acc</sub> > 2 Myr, t<sub>start</sub> ≳ 3.5 Myr<br>α ≥ 5, τ<sub>acc</sub> ∈ [0.6, 1] Myr, t<sub>start</sub> ≳ 3 Myr</td><td></td></tr><tr><td>Different values</td><td>α = 4: τ<sub>acc</sub> > 2 Myr, t<sub>start</sub> ≳ 3.5 Myr<br>α ≥ 5: τ<sub>acc</sub> > 0.6 Myr, t<sub>start</sub> ≳ 3 Myr</td><td>α = 4: no conditions on τ<sub>acc</sub> and t<sub>start</sub><br>α ≥ 5: τ<sub>acc</sub> < 1 Myr or t<sub>start</sub> ≲ 3 Myr</td></tr></tbody></table>

($τ_{acc}$, $t_{start}$, $α$). The main differences in their formation conditions are the disk temperature and the final radius, which result in higher impact velocities for Ganymede. These higher velocities favor more efficient heating and differentiation in Ganymede compared to Callisto. Additional simulations show that the dichotomy between the two moons can also be explained if they accrete with the same mass flux, $M_{Sat}$, meaning that Ganymede completes its accretion slightly later than Callisto.

For this dichotomy to be primordial, our results suggest that both moons likely accreted gradually ($τ_{acc} > 2$ Myr), starting after $t_{start} \gtrsim 3.5$ Myr, with $α \sim 4$. The preferred size distribution in this study lies between the two extremes, $α = 1$ (where the impactor mass is concentrated in large bodies) and $α = 6$ (where the impactors are smaller). Consequently, we propose that satellite bodies with radii $r_{imp} \gtrsim 1$ km were likely present in the Jovian CPD during satellite formation, contributing up to 30% of the total mass influx. These satellitesimals may represent collisional debris or trapped planetesimals that have lost a significant fraction of their mass through ablation (T. Ronnet & A. Johansen 2020), while the bulk of the mass is contained in smaller particles ranging from centimeters to tens of meters.

The values proposed here cannot constrain the size distribution of particles within the PSN at Jupiter's orbit, as the majority of planetesimals would have been ablated by the gas giant, contributing to a significant reservoir of pebbles within the Jovian subnebula. In this context, the size distribution of the PSN at Jupiter's orbit was likely modified and evolved toward steeper slopes during CPD accretion, where the mass influx was dominated by small particles. If the building blocks of satellites formed during the earliest stages of the CPD, our preferred value of $α \sim 4$ aligns with the pebble accretion scenario (T. Ronnet et al. 2017; T. Ronnet & A. Johansen 2020), challenging the view that collisions between satellitesimals dominate the growth process.

K. Batygin & A. Morbidelli (2020) did not include pebble accretion in their model, focusing instead on oligarchic growth, which resulted in a size-frequency distribution with a significantly flatter slope than $α = 3$. Even when considering debris from disk collisions, their model predicts an almost complete mass fraction of large impactors, leading to global melting of both moons. However, the size distribution inferred in this study may represent a second-generation distribution, in which the observed satellites accreted later from an evolved size-frequency distribution in which particles underwent both effective growth and collisional erosion. This interpretation is consistent with the measured value for main-belt asteroids, $α \sim 3.5$ (B. Zellner 1979), and suggests that oligarchic growth cannot be ruled out for the initial formation of large structures that were subsequently processed and eroded by collisional debris.

Ganymede is unlikely to avoid melting during accretion. Due to its larger radius and higher disk temperature compared to Callisto, resulting from its closer proximity to Jupiter, Ganymede can only prevent melting by accreting slowly ($τ_{acc} \gtrsim 1$ Myr) and from small particles with radii $r_{imp}$ around 100 m or less (see Table 3). Accreting even a small number of large impactors buries energy deep within the moon, where it cannot efficiently dissipate to the surface. Consequently, neither prolonged accretion times nor delayed formation would be sufficient to prevent melting, as removing the energy deposited at these depths would require a timescale of about $10^8$ yr. This finding contrasts with the traditional view that larger impactors require longer accretion times or delayed formation to avoid melting, as discussed in A. Barr & R. Canup (2008) and K. Batygin & A. Morbidelli (2020).

One could envision a scenario where Ganymede and Callisto formed exclusively from small bodies with radii of $r_{imp} \lesssim 100$ m and had an accretion timescale of $τ_{acc} \gtrsim 1$ Myr, which aligns well with the Canup & Ward model (R. Canup & W. Ward 2002; A. Barr & R. Canup 2008). However, the possibility of Ganymede and Callisto accreting only small particles is low, given the presence of kilometer-sized satellitesimals in the CPD. T. Ronnet & A. Johansen (2020) demonstrated that about 10% of planetesimals captured from the PSN avoided complete mass loss during ablation, retaining significant sizes of several kilometers. However, the current size-frequency distribution of craters on satellite surfaces may not provide additional constraints on impactor sizes during accretion, as no major satellites have surface ages that extend back to the CPD era (W. F. Bottke et al. 2024). The impact record was likely erased owing to resurfacing events and subsequent cratering by heliocentric impactors, with Jupiter-family comets accounting for over 90% of the craters observed on the Galilean satellites (K. Zahnle et al. 1998).

Nevertheless, if both Callisto and Ganymede accreted small bodies, the dichotomy between the two moons may eventually result from divergent evolutionary pathways. A. C. Barr et al. (2010) proposed that the late heavy bombardment (LHB) could explain the differences, with Ganymede possibly experiencing more energetic impacts than Callisto owing to probabilistic variations. However, a recent study by W. F. Bottke et al. (2024) questions this scenario, as their estimated mass flux of impactors during the LHB is five times lower than that assumed by A. C. Barr et al. (2010), reducing the likelihood of Ganymede melting owing to impacts. Similarly, F. Nimmo & D. Korycansky (2012) explain that the mass delivered to the satellites during the post-accretional bombardment is too high to be consistent with observations. To reconcile this discrepancy, the mass inflow from the destabilized population must be reduced by a factor of 10. Otherwise, Mimas, Miranda, and Enceladus would lose all of their volatiles, assuming that these icy bodies formed well before the bombardment.
Resonances could also explain the surface differences between Ganymede and Callisto. A. P. Showman & R. Malhotra (1997) and S. J. Peale & M. H. Lee (2002) proposed the existence of an evolutionary path to the Laplace resonance, where Ganymede's eccentricity increases to 0.07 before settling to its present value. Consequently, A. P. Showman et al. (1997) suggested that the significant tidal heating Ganymede experienced during this period would have led to resurfacing and the formation of the young, bright terrains observed today. However, the authors assume that Ganymede had a fully differentiated initial structure with a warm icy mantle before entering the Laplace resonance. This result may be consistent with our study, as the paper also suggests that differentiation may have occurred early, during the accretion phase, and before the formation of the dark terrains.

Given the uncertainty in the gas density distribution within the CPD, which directly influences the formation timescales, Callisto and Ganymede could also have fundamentally different accretion parameters, especially $\tau_{acc}$ and $t_{start}$. In this situation, the dichotomy could still be primordial, as illustrated by the various scenarios summarized in Table 3. In addition, Europa and Io probably accreted in the same satellitesimal feeding region as Ganymede (K. Batygin & A. Morbidelli 2020), implying similar or even shorter accretion timescales ($\tau_{acc} \sim 10^3-10^4$ yr; P. R. Estrada et al. 2009). In this scenario, melting of the two inner Galilean moons seems inevitable, probably leading to rock-iron differentiation (K. T. Trinh et al. 2023).

The onset time of Callisto's accretion, $t_{start}$, can also serve as a constraint on the timing of Jupiter's formation. In the scenarios outlined in Table 3, we find that $t_{start}$ for Callisto exceeds 3.5 Myr for $\alpha = 5$ and 3 Myr for $\alpha = 4$. Assuming that the formation of Jupiter ended with the formation of the Jovian CPD and that Callisto began accreting immediately thereafter, the formation of Jupiter probably ended no earlier than ~3 and 3.5 Myr after CAIs for $\alpha = 4$ and $\alpha = 5$, respectively. This result supports the late formation hypothesis for Jupiter, suggesting that the gas giant continued to accrete material from the PSN until about 3-5 Myr after the formation of the solar system (T. S. Kruijer et al. 2017; T. Kleine et al. 2020).

Here our temperature profiles are based on the disk model described in A. Schneeberger & O. Mousis (2025). The objective of this study was to couple accretional heating with the disk temperature evolution $T_c(t)$, to examine how the cooling of the disk at the end of the accretion phase impacts satellite formation, rather than introducing new free parameters or assuming constant temperature values throughout the satellite formation process. However, disk temperatures can exhibit significant variation across different models, with key influences from its opacity and the rate of material inflow. While these temperature variations could potentially impact the accretional temperature profiles, we find that this effect becomes less pronounced when the number of small particles decreases ($\alpha < 4$), as the surface would experience less heating compared to the interior, making it less sensitive to ambient temperature changes. Thus, if a high number of large impactors are present in the Jovian CPD, we expect that the temperature profiles inside protosatellites may be similar to the results presented in this work.

Furthermore, our model assumes in situ formation of both moons, although inward migration of the satellites may have occurred during their formation. According to R. M. Canup & W. R. Ward (2006), the satellites we observe today represent

the last generation that likely avoided Type I migration, since they formed after the material inflow ceased. Specifically, the simulations of R. M. Canup & W. R. Ward (2006) suggest that Callisto analogs migrated inward by less than 5 Jupiter radii, an effect that would have had only a minor impact on the disk temperature profiles. In contrast, I. Mosqueira & P. R. Estrada (2003a) predict that Callisto's building blocks originated in the cold outer regions of the disk, about 100 Jupiter radii. However, the satellite embedded model of I. Mosqueira & P. R. Estrada (2003a, 2003b) predicts a significant fraction of impactors with radii in the 1-100 km range. As noted earlier, within this typical size range, the accretion temperature profiles are not significantly affected by the disk temperature, since most of the energy is retained in the interior of the satellite.<fnref n="5" /> Therefore, even if Type I migration had occurred during satellite formation at a lower disk temperature, the results should not be significantly different from those presented here.

This study does not explicitly model water transport or water-rock separation, as we assume that differentiation is initiated once a few percent of the material exceeds the melting point of water ice. As discussed in Section 3, impacts can produce localized melt pools that rapidly refreeze over timescales of several tens of thousands of years, depending on the melt volume. Given that melting occurs near the end of accretion and is predominantly confined to the surface, we do not expect substantial melt migration before refreezing is complete.

Our model does not include plastic work heating—the conversion of mechanical energy to heat due to irreversible plastic deformation, where a fraction of the plastic strain energy is dissipated as heat. Recent studies have shown that plastic deformation during decompression can contribute to energy deposition in impacts, especially in low-velocity collisions (K. Kurosawa & H. Genda 2018; K. Kurosawa et al. 2021). This mechanism could lead to localized surface melting and devolatilization, potentially influencing the thermal evolution of icy bodies such as Ganymede and Callisto. However, incorporating plastic work heating into our model would introduce additional complexity by requiring the tracking of local deformation, which falls outside the scope of this study. Instead, we prioritize shock wave heating owing to its greater impact on thermal changes during collisions. Future work that integrates plastic work heating into accretion models could provide a more accurate representation of impact processes and provide deeper insights into the complex thermal dynamics of ice-rock environments.

## Acknowledgments

We thank the two anonymous reviewers for their helpful comments. O.M. acknowledges support from CNES. This research holds as part of the project FACOM (ANR-22-CE49-0005-01 ACT) and has benefited from funding provided by l'Agence Nationale de la Recherche (ANR) under the Generic Call for Proposals 2022. V.H. acknowledges support from the French government under the France 2030 investment plan, as part of the Initiative d'Excellence d'Aix-Marseille Université —A*MIDEX AMX-22-CPJ-04.

<footnote n="5">⁵ The assumption of an impactor size distribution containing only a few 100 km sized satellitesimals would also be sufficient to induce significant differentiation of Callisto, in contradiction with the gravity data.</footnote>
ORCID iDs

Yannis Bennacer https://orcid.org/0009-0000-5202-7107

Olivier Mousis https://orcid.org/0000-0001-5323-6453

Marc Monnereau https://orcid.org/0000-0002-4094-5581

Vincent Hue https://orcid.org/0000-0001-9275-0156

Antoine Schneeberger https://orcid.org/0000-0002-3289-2432

References

Anderson, J., Jacobson, R., McElrath, T., et al. 2001, Icar, 153, 157

Barr, A., & Canup, R. 2008, Icar, 198, 163

Barr, A. C., & Canup, R. M. 2010, NatGe, 3, 164

Barr, A. C., Citron, R. I., & Canup, R. M. 2010, Icar, 209, 858

Batygin, K., & Morbidelli, A. 2020, ApJ, 894, 143

Bottke, W. F., Vokrouhlick, D., Nesvorn, D., et al. 2024, PSJ, 5, 88

Canup, R., & Ward, W. 2002, ApJ, 124, 3404

Canup, R. M., & Ward, W. R. 2006, Natur, 441, 834

Carslaw, H. S., & Jaeger, J. C. 1959, Conduction of Heat in Solids (Oxford: Clarendon)

Castillo-Rogez, J., Matson, D., Sotin, C., et al. 2007, Icar, 190, 179

Chivers, C. J., Buffo, J. J., & Schmidt, B. E. 2021, JGRE, 126, e2020JE006692

Croft, S. K. 1982, Geological Implications of Impacts of Large Asteroids and Comets on the Earth (Boulder, CO: Geological Society of America)

Downey, B., Nimmo, F., & Matsuyama, I. 2020, MNRAS, 499, 40

Estrada, P. R., Mosqueira, I., Lissauer, J. J., D'Angelo, G., & Cruikshank, D. P. 2009, arXiv:0809.1418

Heller, R., Marleau, G. D., & Pudritz, R. E. 2015, A&A, 579, L4

Hesse, M. A., & Castillo-Rogez, J. C. 2019, GeoRL, 46, 1213

Kivelson, M. G., Warnecke, J., Bennett, L., et al. 1998, JGRE, 103, 19963

Kleine, T., Budde, G., Burkhardt, C., et al. 2020, ApJ, 216, 3404

Klinger, J. 1980, Sci, 209, 271

Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, PNAS, 114, 6712

Kurosawa, K., & Genda, H. 2018, GeoRL, 45, 620

Kurosawa, K., Genda, H., Azuma, S., & Okazaki, K. 2021, GeoRL, 48, e2020GL091130

Leliwa-Kopystyński, J., Maruyama, M., & Nakajima, T. 2002, Icar, 159, 518

Lewis, J. S. 1971, Icar, 15, 174

Lodders, K. 2003, ApJ, 591, 1220

Lodders, K. 2019, arXiv:1912.00844

Makalkin, A. B., & Dorofeeva, V. A. 1995, AVest, 29, 99

Makalkin, A. B., & Dorofeeva, V. A. 2014, SoSyR, 48, 62

Malhotra, R. 1991, Icar, 94, 399

McCord, T. B., & Sotin, C. 2005, JGRE, 110, E05009

McKinnon, W. B. 1997, Icar, 130, 540

Michaut, C., & Manga, M. 2014, JGRE, 119, 550

Mishra, R. K., & Chaussidon, M. 2014, E&PSL, 398, 90

Monnereau, M., Toplis, M. J., Baratoux, D., & Guignard, J. 2013, GeCoA, 119, 302

Monteux, J., Coltice, N., Dubuffet, F., & Ricard, Y. 2007, GeoRL, 34

Monteux, J., Tobie, G., Choblet, G., & Le Feuvre, M. 2014, Icar, 237, 377

Morris, E. M. 2001, AntSc, 13, 350

Mosqueira, I., & Estrada, P. R. 2003a, Icar, 163, 198

Mosqueira, I., & Estrada, P. R. 2003b, Icar, 163, 232

Mousis, O., Pargamin, J., Grasset, O., & Sotin, C. 2002, GeoRL, 29, 45

Nagel, K., Breuer, D., & Spohn, T. 2004, Icar, 169, 402

Naseem, M., Neveu, M., Howell, S., et al. 2023, PSJ, 4, 181

Nimmo, F., & Korycansky, D. 2012, Icar, 219, 508

Okeefe, J. D., & Ahrens, T. J. 1977, LPSC, 3, 3357

Peale, S. J., & Lee, M. H. 2002, Sci, 298, 593

Ronnet, T., & Johansen, A. 2020, A&A, 633, A93

Ronnet, T., Mousis, O., & Vernazza, P. 2017, ApJ, 845, 92

Rubin, M. E., Desch, S. J., & Neveu, M. 2014, Icar, 236, 122

Schneeberger, A., & Mousis, O. 2025, PSJ, 6, 23

Schubert, G., Stevenson, D., & Ellsworth, K. 1981, Icar, 47, 46

Senshu, H., Kuramoto, K., & Matsui, T. 2002, JGRE, 107, 5118

Shoji, D., & Kurita, K. 2014, JGRE, 119, 2457

Showman, A. P., & Malhotra, R. 1997, Icar, 127, 93

Showman, A. P., Stevenson, D. J., & Malhotra, R. 1997, Icar, 129, 367

Spohn, T., & Schubert, G. 2003, Icar, 161, 456

Squyres, S. W., Reynolds, R. T., Summers, A. L., & Shung, F. 1988, JGR, 93, 8779

Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853

Thrane, K., Bizzarro, M., & Baker, J. 2006, ApJ, 646L, 159

Tonks, W. B., Pierazzo, E., & Melosh, H. J. 2016, arXiv:1607.01282

Trinh, K. T., Bierson, C. J., & O'Rourke, J. G. 2023, SciA, 9, eadf3955

Waite, W. F., Gilbert, L. Y., Winters, W. J., & Mason, D. H. 2006, RScl, 77, 044904

Wisdom, J. 2004, AJ, 128, 484

Wisdom, J. 2008, Icar, 193, 637

Yomogida, K., & Matsui, T. 1983, JGR, 88, 9513

Zahnle, K., Dones, L., & Levison, H. F. 1998, Icar, 136, 202

Zellner, B. 1979, in Asteroids, ed. T. Gehrels & M. S. Matthews (Tucson, AZ: Univ. Arizona Press), 783