© ESO 2020

1 Introduction

The formation of massive moons around gas giant planets is envisioned to take place in a gaseous disk surrounding the planet in the last stages of its accretion (see, e.g., Peale & Canup 2015, for a review). Their formation would thus be analogous to that of planets around stars and face the same theoretical challenges. These difficulties include, for example, the formation of satellitesimals to seed the accretion of more massive moons and the survival of these latter against their rapid inward migration. In addition, some issues, such as the delivery of solids to the circum-planetary disk (hereafter CPD), are specific to the formation of satellites around the giant planets. The observed properties of the Galilean moons orbiting Jupiter, as well as that of Saturn’s moons, provide additional constraints on the conditions under which they have formed.

The satellite systems of Jupiter and Saturn both represent a similar fraction of the mass of their parent planets (~10⁻⁴ M_p) and are quite compact (e.g., Canup & Ward 2006). In the case of the Galilean system, the three inner moons – Io, Europa and Ganymede – form a resonant chain (known as a Laplace resonant system). Additionally, the decreasing bulk densities of the Galilean satellites with respect to their distance from Jupiter is suggestive of a compositional gradient among the moons, that is, an increasingly volatile rich composition with respect to the distance from Jupiter (see, e.g., Hussmann et al. 2015). Finally, information regarding the internal structure of the moons has been inferred from the gravity measurements performed by the Galileo spacecraft at the Jovian system and the Cassini spacecraft in the case of Saturn’s moons. While the three inner Galilean moons (involved in the Laplace resonance system) are likely to be fully differentiated with the presence of an iron core, a silicate mantle and, in the case of Europa and Ganymede, an icy outer mantle, Callisto and Titan appear to be only partially differentiated (Anderson et al. 2001; Iess et al. 2010). The formation of undifferentiated satellites implies limited heating during accretion to prevent large scale ice melting and hence requires long formation timescales (≳1 Myr; Barr & Canup 2008). Even so, a later differentiation of the satellites could be difficult to avoid if compositional gradients prevent an efficient transport of radiogenic heating through convection (O’Rourke & Stevenson 2014). Also, nonhydrostatic effects could have an important impact on the derivation of the moment of inertia of slowly rotating satellites, such as Callisto and Titan (Gao & Stevenson 2013), so that the interpretation of their internal structure remains uncertain and these objects might in fact be differentiated. If this is the case, these satellites could have formed over much shorter timescales.

The consideration of the properties of the satellite systems, together with the available knowledge on processes relevant to planet formation, guided the development and refinement of satellite formation scenarios (e.g., Lunine & Stevenson 1982; Coradini et al. 1995; Canup & Ward 2002, 2006; Mosqueira & Estrada 2003a,b; Mousis et al. 2002; Sasaki et al. 2010). Although the gas-starved model developed by Canup & Ward (2002, 2006) has been well recognized as a plausible satellite formation scenario, in recent years, new paradigms have emerged regarding several key processes for the formation of planets, which challenge our current understanding of the formation of the giant planets’ moons.

These challenges are briefly discussed in Sect. 2, where we argue that the capture and ablation of planetesimals should be an important source of solids in giant planet’s CPD, as previously suggested (e.g., Estrada et al. 2009), but unlike the assumption of the gas-starved model. Such a mechanism is investigated in Sect. 4, which shows that most of the planetesimals captured in the CPD are efficiently ablated. Ultimately, the material ablated off of the surface of the planetesimals provides a source of dust in the CPD whose subsequent evolution is investigated in Sect. 5. These results provide the ground for the development of a revised formation model for the giant planets’ satellites (Sect. 6). Specifically, we propose that the seeds of the satellites initially form from the fraction of captured planetesimals that survived ablation in the CPD and subsequently grow through pebble accretion from the flux of dust supplied by the ablation of planetesimals. More aspects and implicationsare discussed in Sect. 7, before we summarize our results in Sect. 8.

2 Motivation

The so-called gas-starved model developed by Canup & Ward (2002, 2006) has been the leading scenario for the formation of satellites around gas giant planets because of its apparent ability to resolve long-standing issues which are primarily (i) the impossibility of satisfying the minimum mass requirement to grow the moons and the low temperature needed for ice stability simultaneously in a “self-consistent” viscous disk model (Makalkin et al. 1999; Canup & Ward 2002), (ii) the difficulty of forming undifferentiated massive satellites and (iii) the survival of the satellites against their inward migration.

In this scenario, the formation of the satellites takes place as the giant planet is still accreting material from the surrounding protoplanetary disk (PPD). The CPD is not a closed system then, but is instead constantly replenished by fresh material with an approximately solar composition (dust-to-gas mass ratio ε_d,⊙ ≈ 10⁻²). In this case, allthe mass required to grow the satellites need not be present at one time in the CPD, which allows to consider less dense environments which are more propitious to cold temperature and the formation of icy satellites. The accretion timescale of the moons would be regulated by the rate of inflow of fresh material onto the CPD, and the migration timescales would be lengthen due to the lower gas densities.

However, it should be noted that the validity of viscous disk models has been challenged in the recent years (see, e.g., Turner et al. 2014; Gressel et al. 2015; Bai 2017) due to the fact that non-ideal MHD effects tend to suppress the source of turbulent viscosity in the disk, and their evolution would then rather be driven by thermo-magnetic winds. This possibly results in inefficient viscous heating of the disks (Mori et al. 2019). The non-ideal MHD effects are expected to dominate even more in the case of CPDs (Fujii et al. 2014, 2017), so that CPDs could be denser than advocated by Canup & Ward (2002) while remaining cold. Moreover, the truncation of the giant planets’ CPD by an inner magnetic cavity, as seems to be required to explain their rotation rate (Takata & Stevenson 1996; Batygin 2018), would stop the migration of the satellites (Sasaki et al. 2010; Ogihara & Ida 2012). It thus seems that a gas-starved environment is not essential to allow for the formation of icy satellites and their survival.

On the other hand, one of the important issues of the gas-starved model is the formation of satellitesimals and satellites seeds. Whereas it has been customarily assumed that large objects would form out of the small dust grains in the CPD (e.g., Canup & Ward 2002, 2006; Sasaki et al. 2010; Ogihara & Ida 2012), it is now understood that the formation of planetesimals or satellitesimals likely requires some instability (e.g., streaming instabilities) or adequate environment (e.g., pressure bumps) to allow for the efficient concentration of dust that can then gravitationally collapse into 10–100 km sized objects (see Johansen et al. 2014, for a review). The rapid inward drift timescale of dust grains in a CPD precludes the formation of satellitesimals through any known mechanism (Shibaike et al. 2017). Even if satellitesimals were somehow able to form (see, for example, Drążkowska & Szulágyi 2018)¹, it is not expected that all the dust grains would be converted into large objects, contrary to the assumption in previous investigations. Rather, the satellitesimals should further grow by accreting the remaining (and possibly re-supplied) dust grains in a process known as pebble accretion (see the recent reviews by Johansen & Lambrechts 2017; Ormel 2017). The accretion of moons in CPDs thus needs to be reassessed considering this more likely growth channel.

Another important issue is that it is unlikely that the material accreted by a giant planet in the late stages of its formation has a solar dust-to-gas mass ratio, as advocated in the gas-starved models (see Ronnet et al. 2018, for a discussion). Above a mass of typically a few ~10 M_⊕, a protoplanet starts to significantly perturb the gas disk in its vicinity, opening a shallow gap whose outer edge acts as a barrier for drifting dust grains (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014). At a mass comparable to that of Jupiter, dust grains larger than a few 10–100 μm should be efficiently filtered (Paardekooper 2007; Zhu et al. 2012; Bitsch et al. 2018; Weber et al. 2018), likely resulting in a drop of the dust-to-gas mass ratio of the accreted material at the order of magnitude level. Combined with the fact that only a small fraction of the dust grains would be accreted by the forming moons (see Sect. 7.2 for an estimate), it becomes difficult to envision that the inflow of small dust grains entrained with the gas constitutes the main source of material for the formation of moons.

Estrada et al. (2009) proposed that gas-drag assisted capture and ablation of planetesimals on initially heliocentric orbits in the vicinity of the giant planets could have provided the solids necessary to form the moons. This mechanism has recently received some attention (Mosqueira et al. 2010; Fujita et al. 2013; Tanigawa et al. 2014; D’Angelo & Podolak 2015; Suetsugu & Ohtsuki 2017; Ronnet et al. 2018), although it has not been put in the perspective of constructing a more consistent scenario of the subsequent formation of satellites. Here we therefore investigate the capture and ablation of planetesimals in the CPD of a Jovian mass planet with the aim to parametrize the distribution of material (Sect. 4) so as it could be used as an input for models of the subsequent accretion of satellites (Sects. 5 and 6). We argue that the efficient ablation of planetesimals provides the ground for a pebble accretion scenario for the formation of giant planets’ moons, in many regards similar to the scenario proposed for the formation of compact super-Earth systems (Ormel et al. 2017; Izidoro et al. 2017, 2019; Lambrechts et al. 2019).

We note that a pebble accretion scenario for the origin of the Galilean satellites, which shares many similarities with the ideas presented in Sect. 6, has recently been investigated by Shibaike et al. (2019). However, similarly to Canup & Ward (2002), these authors assumed that dust grains were brought to the CPD with the gas accreted by Jupiter, which is unlikely to provide a sufficient source of solid material as discussed above. Moreover, the coupling between the accretion rate of gas and dust onto the CPD implies that the Galilean moons must have started to form as Jupiter was only 40% of its present mass, and that Jupiter accreted over a very long timescale of 30 Myr, whereas the typical lifetime of protoplanetary disks is only ~3 Myr (Hartmann et al. 2016), to reproduce their current properties. If on the other hand the building blocks of the moons are brought to the CPD through the capture and ablation of planetesimals, as proposed here, the formation of the satellites can take place much later in the accretion history of their parent planet, and these latter are not constrained to form on several tens of Myr.

3 Circum-planetary disk model

Recent findings suggest that the turbulence level inside a gaseous CPD is low (Fujii et al. 2014, 2017), and that the midplane of such disks are inefficiently heated by internal dissipation mechanisms (Mori et al. 2019). We thus consider here the case of a passively irradiated CPD.

We assume that the surface density of the disk follows a power-law such that $${\Sigma }_{\text{g}}={\Sigma }_{\text{out}}{\left(\frac{r}{r_{\text{out}}}\right)}^{-\gamma },$$(1)

where r is the radial distance to the planet in cylindrical coordinates, r_out is the outer radius of the CPD and $${\Sigma }_{\text{out}}=\frac{(2-\gamma )M_{\text{CPD}}}{2\pi r_{\text{out}}^2}$$(2)

is the surface density of gas at the outer edge of the disk and M_CPD is the total mass of the CPD. Here we adopt r_out = 0.2 R_H (with R_H the size of the Hill sphere of the planet) as reported by Tanigawa et al. (2012) from the analysis of their 3D hydrodynamic simulations. Although we note that there is no sharp boundary at the outer edge of the CPD allowing to precisely define its radial extent (e.g., Tanigawa et al. 2012; Szulágyi et al. 2014), both the capture and ablation of planetesimals and the subsequent growth of satellites take place in the inner part of the CPD in our scenario (e.g., mainly interior to r ≲ 0.05 R_H). We assume that the mass of the CPD represents some fraction of the mass of the planet, M_CPD = 1.5 × 10⁻³ M_p, in agreement with recent hydrodynamic simulations (Szulágyi 2017). We set the index of the power-law of the gas surface density at γ = 1.5, which is also in agreement with the simulation results of Tanigawa et al. (2012) and Szulágyi (2017)².

The midplane temperature T_d of a passively irradiated disk depends on the angle at which the light from the central radiating object strikes the disk, or grazing angle. It can be approximated by (Ruden & Pollack 1991) $$T_{\text{d}}=T_{\text{p}}{\left[\frac{2}{3\pi }{\left(\frac{R_{\text{p}}}{r}\right)}^3+\frac{1}{2}{\left(\frac{R_{\text{p}}}{r}\right)}^2\frac{H_{\text{ph}}}{r}\left(\frac{\text{d}\ln H_{\text{ph}}}{\text{d}\ln r}-1\right)\right]}^{1/4},$$(3)

where T_p and R_p are the temperature and radius of the central planet, respectively, and H_ph is the photospheric height of the disk at which radiation from the planet is intercepted. This latter quantity is usually assumed to be proportional to the pressure scale height of the disk, H_ph = χH_g, where H_g = c_g∕Ω_K (with c_g the isothermal speed of sound in the gas and Ω_K the orbital frequency), and the proportionality factor χ, which is of order unity, and depends on the opacity of the disk (e.g., Chiang & Goldreich 1997; Dullemond et al. 2001). Assuming that gas accretion onto the giant planet is still ongoing, T_p is related to the mass accretion rate onto the planet, since the planet’s luminosity, $L_{\text{p}}=4\pi R_{\text{pl}}^2{\sigma }_{\text{sb}}{T}_{\text{p}}^4$, is then given by the accretion luminosity L_p ≈ GM_pṀ_p∕R_p (e.g., Mordasini et al. 2017). To remain as general as possible, it is convenient to parameterize the accretion rate onto the planet through the accretion timescale τ_acc ≡ M_p∕Ṁ_p. The planet’s radiating temperature is then $T_{\text{p}}\propto M_{\text{p}}^{1/2}{R}_{\text{p}}^{-3/4}{\tau }_{\text{acc}}^{-1/4}$.

We can now use Eq. (3) to derive some scaling of the temperature of the CPD as a function of the distance from the planet. The first term inside the brackets of Eq. (3) accounts for the finite angular size of the radiating planet, whereas the second term arises from the flaring geometry of the disk. In the inner portions of the CPD, where H_ph ≲ R_p, the first term will dominate, while at larger distances, the grazing angle is instead determined by the flaring of the disk. There are thus two asymptotic solutions for the temperature of the CPD. In the inner, geometrically thin (with respect to the planet’s size), portion of the CPD, the temperature follows $$T_{\text{thin}}\approx 190{\left(\frac{M_{\text{p}}}{M_{\text{Jup}}}\right)}^{1/2}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{-1/4}{\left(\frac{r}{10R_{\text{Jup}}}\right)}^{-3/4}\text{K},$$(4)

which has the same dependency on distance (T ∝ r^−3∕4) than viscousdisk models (e.g., Ruden & Pollack 1991; Canup & Ward 2002). In the regions where flaring dominates instead, the temperaturebecomes $$\begin{array}{ll}{T}_{\text{flar}}\approx \hfill & 190{\left(\frac{r}{10R_{\text{Jup}}}\right)}^{-3/7}\hfill \\ \hfill & \times {\left(\frac{M_{\text{p}}}{M_{\text{Jup}}}\right)}^{3/7}{\left(\frac{R_{\text{p}}}{1.6R_{\text{Jup}}}\right)}^{-2/7}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{-2/7}{\left(\frac{\chi }{4}\right)}^{2/7}\text{K},\hfill \end{array}$$(5)

which follows the classical result derived by Chiang & Goldreich (1997) that T ∝ r^−3∕7. The temperature of the CPD can then be approximated by T_d = max(T_thin, T_flar). By equating the two asymptotic solutions, we find that the transition between the two regimes should occur at a distance $$r_{\text{tran}}\approx 10{\left(\frac{M_{\text{p}}}{M_{\text{Jup}}}\right)}^{2/9}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{1/9}{\left(\frac{R_{\text{p}}}{1.6R_{\text{Jup}}}\right)}^{8/9}{\left(\frac{\chi }{4}\right)}^{-8/9}{R}_{\text{Jup}}.$$(6)

The resulting aspect ratio of the CPD, h_g ≡ H_g∕r, is $$h_{\text{thin}}\approx 0.06{\left(\frac{M_{\text{p}}}{M_{\text{Jup}}}\right)}^{-1/4}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{-1/8}{\left(\frac{r}{10R_{\text{Jup}}}\right)}^{1/8},$$(7)

when r ≤ r_tran, and $$\begin{array}{lll}{h}_{\text{flar}}\approx \hfill & \hfill & 0.06{\left(\frac{r}{10R_{\text{Jup}}}\right)}^{2/7}\hfill \\ \hfill & \hfill & \times {\left(\frac{M_{\text{p}}}{M_{\text{Jup}}}\right)}^{-2/7}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{-1/7}{\left(\frac{R_{\text{p}}}{1.6R_{\text{Jup}}}\right)}^{-1/7}{\left(\frac{\chi }{4}\right)}^{1/7},\hfill \end{array}$$(8)

for r > r_tran.

In the above estimates, we have used a rather long accretion timescale, which is similar to the value adopted by Canup & Ward (2002, 2006) at the formation epoch of the Galilean moons, and considered a slightly inflated planet, as expected at the end of the runaway contraction phase (e.g., Lissauer et al. 2009; Mordasini et al. 2017). The resulting luminosity, for a Jovian mass planet, would be L_p ≈ 4 × 10⁻⁵ L_⊙, which is also consistent with numerical simulations (Lissauer et al. 2009; Mordasini et al. 2017). For simplicity, we have here adopted a constant photospheric to pressure scale height ratio χ = 4, corresponding to the typical value reported by Chiang & Goldreich (1997) for a disk with a solar metallicity. However, we elaborate in Sect. 5 on the fact that the dust-to-gas mass ratio in a CPD is likely to be much lower than expected from a solar composition mixture. A further improvement would thus be to self-consistently determine χ (as in Dullemond et al. 2001) based on the dust spatial and size distribution, but we leave such considerations to future investigations.

Applying this passive CPD model to the case of Saturn, we can note some differences in the resulting thermal structure in the region of satellites formation. The temperature scales to the Saturn system as $$T_{\text{thin}}\approx 120{\left(\frac{M_{\text{p}}}{M_{\text{Sat}}}\right)}^{1/2}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{-1/4}{\left(\frac{r}{10R_{\text{Sat}}}\right)}^{-3/4}\text{K},$$(9) $$\begin{array}{lll}{T}_{\text{flar}}\approx \hfill & \hfill & 130\text{\hspace{0.17em}K\hspace{0.17em}}{\left(\frac{r}{10R_{\text{Sat}}}\right)}^{-3/7}{\left(\frac{M_{\text{p}}}{M_{\text{Sat}}}\right)}^{3/7}\hfill \\ \hfill & \hfill & \times {\left(\frac{R_{\text{p}}}{1.6R_{\text{Sat}}}\right)}^{-2/7}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{-2/7}{\left(\frac{\chi }{4}\right)}^{2/7},\hfill \end{array}$$(10)

which is overall colder than in the case of the Jovian disk. The resulting aspect ratio of the CPD is on the contrary larger, due to the lower planet’s mass, where $$h_{\text{thin}}\approx 0.079{\left(\frac{M_{\text{p}}}{M_{\text{Sat}}}\right)}^{-1/4}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{-1/8}{\left(\frac{r}{10R_{\text{Sat}}}\right)}^{1/8},$$(11) $$\begin{array}{ll}{h}_{\text{flar}}\approx \hfill & 0.082{\left(\frac{r}{10R_{\text{Sat}}}\right)}^{2/7}\hfill \\ \hfill & \times {\left(\frac{M_{\text{p}}}{M_{\text{Sat}}}\right)}^{-2/7}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{-1/7}{\left(\frac{R_{\text{p}}}{1.6R_{\text{Sat}}}\right)}^{-1/7}{\left(\frac{\chi }{4}\right)}^{1/7}.\hfill \end{array}$$(12)

The transition between the two solutions occurs at a distance $$r_{\text{tran}}\approx 7.8{\left(\frac{M_{\text{p}}}{M_{\text{Sat}}}\right)}^{2/9}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{1/9}{\left(\frac{R_{\text{p}}}{1.6R_{\text{Sat}}}\right)}^{8/9}{\left(\frac{\chi }{4}\right)}^{-8/9}{R}_{\text{Sat}}.$$(13)

4 Capture and ablation of planetesimals

4.1 Methods

We model the evolution of a swarm of planetesimals embedded in a protoplanetary disk and orbiting close to a Jupiter-mass planet located at 5.5 au from a Sun-mass star. Figure 1 shows the surface density of gas in the vicinity of the planet (see below for details) along with the region where planetesimals were initially distributed, which is shown as a gray shaded area and corresponds to the so-called feeding zone of the planet. Following Shiraishi & Ida (2008), we define the feeding zone as the region where the Jacobi energy of the planetesimals E_J is positive. E_J can be expressed as (e.g., Nakazawa et al. 1989; Shiraishi & Ida 2008) $$E_{\text{J}}=\frac{1}{2}\left[{\left(\frac{e}{h}\right)}^2+{\left(\frac{i}{h}\right)}^2\right]-\frac{3}{8}{b}^2+\frac{9}{2}+O(h).$$(14)

In the above expression, $h=R_{\text{H}}/a_{\text{p}}={(M_{\text{p}}/3M_{\odot })}^{1/3}$ is the reduced Hill radius of the planet, with M_p its mass and a_p its orbital distance, e and i are the eccentricity and inclination of the planetesimal, and b is the non-dimensional orbital separation with the planet defined as $$b=\frac{a-a_{\text{p}}}{h{a}_{\text{p}}},$$(15)

where a is the (heliocentric) semimajor axis of the planetesimal. When the eccentricity and inclination distributions of the planetesimals are small (i.e., both e∕h and i∕h ≲ 1), the feeding zone of the planet extends to $b\approx \pm 2\sqrt{3}$. In our simulations, 3000 planetesimals are uniformly randomly distributed in the feeding zone of the planet with initial eccentricities and inclinations drawn from a Rayleigh distribution with a mean of 10⁻³ and 5 × 10⁻⁴, respectively. The planetesimals were considered as massless test particles and the orbits of all the bodies were integrated using the adaptive time stepping IAS15 integrator provided within the REBOUND³ package (Rein & Spiegel 2015). Planetesimals were assigned a physical size to account for the effect of gas drag on their orbital evolution (see next paragraph), and we performed runs for planetesimals with sizes of either 100, 10, or 1 km. The implementation of the effect of gas drag and ablation of planetesimals is described in the following paragraphs.

Fig. 1 Surface density profile of the protoplanetary disk obtained by a 2-dimensional hydrodynamic simulation of a Jupiter mass planet in a disk with a constant aspect ratio HPPD ∕R = 0.05.

4.1.1 Planetesimal dynamics

Since we arehere interested in the gas drag assisted capture of planetesimals, the effect of gas friction on the evolution of the objects was accounted for using the stopping time formalism (e.g., Whipple 1972), where the following acceleration was added to the equation of motion of the planetesimals, $$a_{\text{drag}}=-\frac{v-v_{\text{g}}}{t_{\text{s}}}.$$(16)

Here, v is the velocity vector of the planetesimal and v_g that of the gas. Thestopping time t_s is defined by (e.g., Perets & Murray-Clay 2011; Guillot et al. 2014) $$t_{\text{s}}={\left(\frac{{\rho }_{\text{g}}{v}_{\text{th}}}{{\rho }_{\text{pl}}{R}_{\text{pl}}}\min \left[1,\frac{3}{8}\frac{v_{\text{rel}}}{v_{\text{th}}}{C}_{\text{D}}\right]\right)}^{-1},$$(17)

where R_pl and ρ_pl (=1000 kg m⁻³) are the radius and internal density of the planetesimal, respectively, v_rel is the magnitude of the relative velocity between the planetesimal and the gas, ρ_g is the density of the gas disk (either the PPD or the CPD), $v_{\text{th}}=\sqrt{8/\pi }{c}_{\text{g}}$ is the thermal velocity of the gas, and C_D is the dimensionless drag coefficient. Although C_D is in principle a complicated function of the Reynolds and Mach numbers of the flow around an object (for example, see D’Angelo & Podolak 2015), it is of order unity at large Mach and Reynolds numbers (Tanigawa et al. 2014; Suetsugu & Ohtsuki 2017). Similarly to Suetsugu & Ohtsuki (2017), we therefore opted to set C_D = 1.

4.1.2 Planetesimal thermodynamics

To estimate the effect of the ablation of planetesimals on the delivery of solid material to the CPD, the evolution of the surface temperature of the objects must be followed. To do so, we assume that the surface temperature T_pl of the planetesimals is always at equilibrium with respect to heating from radiation at the ambient temperature of the CPD and friction with the gas, and cooling due to the release of the latent heat of vaporization of water (see Appendix A for details). This yields the following expression for the temperature $$T_{\text{pl}}^4=T_{\text{d}}^4+\frac{C_D{\rho }_{\text{g}}{v}_{\text{rel}}^3}{32{\sigma }_{\text{sb}}}-\frac{P_{\text{v}}(T_{\text{pl}})}{{\sigma }_{\text{sb}}}\sqrt{\frac{{\mu }_{{\text{H}}_{\text{2}}\text{O}}}{8\pi R_{\text{g}}{T}_{\text{pl}}}}{L}_{\text{w}},$$(18)

where T_d is the local temperature of the CPD (assumed to be isothermal in the vertical direction), σ_sb is the Stefan–Boltzmann constant, R_g is the ideal gas constant, ${\mu }_{{\text{H}}_{\text{2}}\text{O}}=0.018\text{\hspace{0.17em}kg\hspace{0.17em}}{\text{mol}}^{-1}$ is the molecular weight of water, L_w = 3 × 10⁶ J kg⁻¹ is the latent heat of vaporization of water and P_v is the saturated vapor pressure of water and is expressed as polynomials of the temperature (taken from Fray & Schmitt 2009). Here we focus on the ablation of water only, assuming that small silicate grains would be entrained with the water released from the surface of the planetesimals. This is only accurate if the planetesimals are not differentiated, and water ice and silicatesare well mixed. If the planetesimals are differentiated into an icy mantle and rocky core, then the ablation of the silicate core should be considered (Mosqueira et al. 2010). We show in Fig. A.1 that such a consideration should not have a great influence on the results presented here (see figure caption for details).

The cooling of the surface due to the ablation of water ice (third term on the r.h.s of Eq. (18)) depends upon the surface temperatureof the object. We must therefore solve Eq. (18) iteratively to determine the surface temperature of planetesimals. When T_pl is known, the mass ablation rate from the surface of the object can be expressed as (e.g., D’Angelo & Podolak 2015), $${\dot{m}}_{\text{abl}}=-4\pi R_{\text{pl}}^2{P}_{\text{v}}(T_{\text{pl}})\sqrt{\frac{{\mu }_{{\text{H}}_{\text{2}}\text{O}}}{2\pi R_{\text{g}}{T}_{\text{pl}}}}.$$(19)

The ablation of material also changes the size of planetesimals, which affects their dynamics through a change in their stopping time (Eq. (17)).

We solve Eq. (18) in between each timestep of the N-body integrator and evaluate the mass loss and radius change due to ablation for each planetesimal that interacts with the CPD (we ignore ablation of material in the PPD). Our approach is not fully numerically self-consistent in the sense that the timestep of the integrator changes according to the dynamical evolution of the objects only, regardless how much mass may be lost (equivalently, size may be reduced) during the timestep. However, as we show in the following, planetesimals are mostly ablated during their closest approach to the planet or, in the case of captured objects, around the pericenter of their orbit. This is also when the integrator timestep is the smallest, on the order of 10⁻⁴ –10⁻⁵ yr∕2π. Most strongly heated objects can reach surface temperatures of T_pl ~ 550 K, which, for a planetesimal with a size R_pl = 100 km, would imply a mass loss on the order of 10⁻² –10⁻³ of its total mass during a typical timestep.

4.1.3 Structure of the protoplanetary disk

When far from the planet (i.e, at distances >0.2 R_H), planetesimals interact with the gas of the PPD. In this case, the gas density and orbital velocity are taken from the results of a 2D hydrodynamic simulation including a Jupiter mass planet, similarly to that used in Ronnet et al. (2018), where an isothermal equation of state has been used for the gas and an α-viscosity (Shakura & Sunyaev 1973) with α = 2 × 10⁻³ was assumed. The surface density of the PPD (Σ_PPD) was normalized so that the unperturbed surface density (that is, in the absence of the planet) at 1 au is ~ 300 g cm⁻², corresponding to a moderately evolved disk (e.g., Bitsch et al. 2015). The resulting surface density is shown in Fig. 1. To obtain the volume density ρ_g of gas around a planetesimal, we assume vertical hydrostatic equilibrium, $${\rho }_{\text{g}}=\frac{{\Sigma }_{\text{PPD}}}{\sqrt{2\pi }{H}_{\text{PPD}}}\exp \left(-\frac{z^2}{2H_{\text{PPD}}^2}\right),$$(20)

where H_PPD is the scale height of the disk and is derived from the fact that the hydrodynamic simulations were performed assuming a constant aspect ratio of the disk H_PPD∕R = 0.05, R being the radial distance to the star in cylindrical coordinates.

4.1.4 Structure of the circum-planetary disk

Planetesimals experiencing a close encounter with the planet (i.e., at distances ≤ 0.2 R_H) interact with the CPD. The surface density of the CPD follows the power-law presented in the previous section with an index γ = 1.5 and a total mass M_CPD = 1.5 × 10⁻³M_p, yielding Σ_out ≈ 2 × 10² g cm⁻². This differs from the work of Suetsugu & Ohtsuki (2017) who scaled their analysis of the capture of planetesimals in CPDs to the gas starved model, assuming Σ_out = 1 g cm⁻². In such a low density disk, the largest planetesimals that can be efficiently captured have a size of R_pl ~ 3 km. However, primordial asteroids of the Solar System seem to record a birth size distribution of objects such that most of the mass was carried by large planetesimals (Morbidelli et al. 2009). This is also consistent with the results of numerical simulations of the formation of planetesimals through the streaming instability, showing that the mass of the planetesimal population is dominated by objects with a radius of ~100 km (Johansen et al. 2015; Simon et al. 2016). For our adopted CPD mass, which is consistent with recent numerical simulations (Szulágyi 2017), hundred-kilometer-sized planetesimals can readily be captured. This alleviates the need to invoke substantial collisional grinding of the initial planetesimal population to allow for their delivery to the giant planets’ CPDs (e.g., Estrada et al. 2009; Mosqueira et al. 2010). We assume vertical hydrostatic equilibrium to obtain the volume density as a function of height in the disk. Since the aspect ratio of the CPD is expected to vary only slightly with respect to distance from the planet in the inner portions of the disk, we here assume an average and constant aspect ratio h_g = 0.06, so that the results presented in this section are not particularly tied to a particular (viscous or passive) CPD model. We nevertheless verified that similar results are obtained when considering a broken power-law such as presented in Sect. 3.

4.2 Results

Here we present the results of the orbital integration of the planetesimals that were initially released on heliocentric orbits within the feeding zone of the giant planet. Planetesimals found on bound orbits around the planet with a (planetocentric) semimajor axis ≤0.05 R_H (~ 40 R_Jup) and an eccentricity <0.1 were removed from the simulation to save computation time. For the range of planetesimal sizes and disk densities explored here, we find that our eccentricity criterion is generally the more restrictive, that is, the semimajor axis of most captured planetesimals shrinks to values ≤ 0.05 R_H before their eccentricities are damped to values <0.1. As shown below, we find that objects that have reached our removal criterion should not contribute further significant amounts of dust through ablation as they have either already lost most of their initial mass, or are found on more distant orbits and experience little frictional heating. We also imposed a minimum cut-off radius of 10 m after which point further ablation was not considered. Since we investigated the evolution of large planetesimals, objects that have reached the cut-off radius have a negligible mass (m_cut-off∕m_init = 10⁻¹²–10⁻⁶ when R_pl = 100–1 km initially).

4.2.1 Distribution of ablated material

Planetesimals are efficiently ablated only when they experience a very close encounter with the planet. The ablation is mainly driven by the strong heating provided by friction with the gas, whereas the temperature of the CPD has a minor effect. If we consider for example a planetesimal approaching the planet at the local escape velocity, in the case of a prograde motion with respect to the rotation of gas in the CPD, the relative velocity between the gas and the planetesimal would be $v_{\text{rel}}=(\sqrt{2}-1)v_{\text{K}}$, or $v_{\text{rel}}=(\sqrt{2}+1)v_{\text{K}}$ in the case of a retrograde approach. For a closest approach at a distance r_CA = 10 R_Jup, ρ_g ≈ 10⁻⁶ g cm⁻³. Using these values to estimate the frictional heating (second term on the r.h.s of Eq. (18)), we find it to be about five orders of magnitude larger than the heating provided by radiation at the local temperature of the CPD (T_d ≈ 190 K). In fact, the temperature of the CPD should have been T_d ≳ 3000 K to provide a source of energy comparable to that of gas drag.

As a consequence, ablation mostly occurs in the densest inner parts of the CPD. This is illustrated in Fig. 2 which shows example trajectories of planetesimals with an initial radius R_pl = 100 km that are captured by the planet in both prograde and retrograde directions. Planetesimals captured on retrograde orbits have a larger relative velocity with respect to the gas than the planetesimals that are captured in the prograde directions. They are therefore more strongly heated and ablated.

Figure 2 reveals that planetesimals captured in the prograde direction have ablation timescales τ_abl≡ m∕ṁ_abl ≲ 10 h (where m = m(t) is the mass of the planetesimal) near their closest approach (top left panel of Fig. 2), whereas planetesimals captured in the retrograde direction have ablation timescales τ_abl ≲ hr (bottom left panel of Fig. 2). These ablation timescales are to be compared with the time “spent” at the closest approach, $$t_{\text{CA}}~\frac{r_{\text{CA}}}{v_{\text{peri}}}\approx 3.5{\left(\frac{r_{\text{CA}}}{5R_{\text{Jup}}}\right)}^{3/2}\text{\hspace{0.17em}h},$$(21)

where we have followed Tanigawa & Ohtsuki (2010) in assuming that the path length of the orbit can be approximated by the distance at closest approach r_CA, and we have used $v_{\text{peri}}\approx \sqrt{2}{v}_{\text{K}}$ (e.g., Fujita et al. 2013) which is a good approximation when ẽ ≈ 1 (here ẽ is the planetocentric eccentricity of the orbit of a planetesimal and is indeed close to unity right after capture).

This simple estimate suggests that a large fraction of the mass of a planetesimal can be ablated during a single passage at pericenter for either prograde or retrograde orbits. Figure 3 shows the time evolution of the distance to the planet and mass of a planetesimal captured in the prograde and retrograde direction, respectively, and confirm this suspicion. From close inspection at Fig. 3, we find that the planetesimal captured in the prograde direction has a fraction of about 40% of its initial mass ablated after its first close encounter with the planet (i.e., upon capture) whereas the planetesimal captured in the retrograde direction has been almost entirely ablated. It is evident from these figures that ablation occurs primarily near the closest approach to the planet and is otherwise negligible (which could also be inferred from Fig. 2).

It is interesting to note that by the time the two planetesimals shown on Fig. 3 are circularized (according to our criterion), they have been almost completely ablated. Their final masses are fractions of only ~10⁻⁵ and 10⁻⁹ that of theirinitial mass. We find that this is actually the fate of the vast majority of captured planetesimals. We report on Fig. 4 the final radii and semimajor axes (i.e., at the time of circularization) of planetesimals that have been captured by the planet. More than 60% of the captured objects have been ablated down to our imposed cut-off radius of 10 m. Those are mainly planetesimals that were captured on initially retrograde orbits. Objects that are captured in the prograde direction may “survive” ablation, with planetesimals captured at wider orbital distances experiencing less ablation (a typical example is the object shown on the top right panel of Fig. 2). We find that about 10% of the captured objects can remain as large objects in the CPD with sizes R_pl > 10 km and they all have a (planetocentric) semimajor axis ã ≳ 10 R_Jup. Although these objects cannot be deemed as completely ablated, they have lost a large fraction of their mass through ablation. Even the largest remaining planetesimal in our simulation (R_pl ≈ 65 km) has lost more than half of its initial mass.

Due to the fact that ablation is efficient only at very close distances from the planet (see Figs. 2 and 3), the main contributors to the supply of material are planetesimals that are captured by the planet or directly collide with it. Planetesimals that merely fly-by the planet and the CPD deliver only a negligible amount of mass through ablation.

We show in Fig. 5 the distribution of ablated material as a function of the distance from the planet in terms of the cumulative ablated mass for planetesimals with three different initial sizes, R_pl = 1, 10 and 100 km. To construct these curves, we divided the CPD in logarithmic bins and added the mass ablated off of planetesimals in each bin at any time. We then considered the (normalized) cumulative of this mass distribution as it appears much less noisy. We stress that these curves indicate where the mass is deposited on average by planetesimals but they do not represent the distribution of material in the CPD at any particular time (such considerations are presented in Sect. 5). We find that the cumulative ablated mass in the CPD can be approximately fit by a function of the form $$F(r;r_0)=1-e^{-r^2/(2r_0^2)},$$(22)

which actually corresponds to a cumulative distribution function of a Rayleigh distribution, where r is the distance to the planet and r₀ is a scale parameter. Light blue curves on Fig. 5 were obtained using Eq. (22) with different values of r₀ as indicated to the right of each curve. Although the fits are far from perfect, they reproduce the “S”-shaped distribution of material and we show in Sect. 5 how Eq. (22) can be related to some physical quantities, allowing us to construct a simple model of the evolution of solids in the CPD.

Fig. 2 Example of orbits of objects that are captured in the prograde (top panels) and retrograde (bottom panels) directions in the planetocentric frame. Each dot corresponds to the position of the planetesimal at a given timestep and the arrows mark the direction at which the objects approach the planet. The left row shows the “typical” orbital evolution during the gas drag assisted capture of planetesimals whereas the right row corresponds to the orbital evolution following lower probability events where planetesimals are captured on orbits with a large pericenter and their semimajor axis decays slowly (see Fig. 4; see also, Suetsugu & Ohtsuki 2017). Insets on the top and bottom left panels are close-up views of the trajectories near the pericenter of the orbits with X, Y ∈ [−10;10]. The color of the dots indicates the ablation rate from the surface of the planetesimal. Ablation occurs mainly around the point of closest approach to the planets where the density of gas is the highest. Due to their larger relative velocities with respect to the gas, and hence larger frictional heating, planetesimals on retrograde orbits (bottom) experience more sustained ablation with ablation timescales as short as τabl ≡ m∕ṁabl ≲hr. The planetesimals captured in the retrograde directions both reached our imposed cut-off radius of 10 m (at this small size, the objects couple with the gas and end up rotating in the prograde direction which can be more clearly seen in the bottom right panel). However, far from the pericenter of the orbit, ablation rapidly becomes negligible for planetesimals orbiting on both prograde and retrograde orbits (τabl ≳ 107–108 h $\gg {\Omega }_{\text{K}}^{-1}$).

4.2.2 Accretion rate onto the circum-planetary disk

Keeping track of the total mass being ablated within the circum-planetary disk, we find that after 100 orbital periods of the planet, it amounts to ~23% of the total mass of planetesimals initially in the feeding zone of the planet in the case of 100 km sized objects. If the initial surface density of planetesimals in the feeding zone would correspond to Σ_pl = 1 g cm⁻², the average accretion rate of solids onto the CPD due to ablation of planetesimals would be ≈ 4.5 × 10⁻⁵ M_⊕ yr⁻¹. At the same time, the feeding zone has been emptied by ≈65%, which means that it would be cleaned up in ~10⁴ yr, irrespective of the initial surface density of planetesimals.

This is problematic if the satellites of the giant planets must form on timescales of 10⁵ –10⁶ yr (e.g., Barr & Canup 2008). However, as shown in Ronnet et al. (2018), the feeding zones of giant planets could be fed themselves as neighboring planets form and destabilize nearby planetesimals, thus allowing the delivery of material to the CPDs on longer timescales.

Fig. 3 Top: evolution of the radial distance from the planet of a planetesimal captured in the prograde direction (panel a), and the retrograde direction (panel b), with an initial size Rpl = 100 km, as a function of time (normalized by the inverse of the orbital angular frequency of the planet ΩK). Bottom: corresponding evolution of the mass of the planetesimal due to ablation in the case of prograde (panel c) and retrograde (panel d) capture. Ablation essentially occurs at closest approach from the planet due to the strong frictional heating in the dense inner parts of the CPD.

Fig. 4 Final radius versus final semimajor axis of planetesimals in the circum-planetary disk (i.e., at the time of circularization of the orbit according to our criteria, ã ≤ 0.05 RH ~ 40 RJup and ẽ ≤ 0.1). Each planetesimal started with an initial radius Rpl = 100 km at the beginning of the simulation. On the right, the fraction of objects with a size Rpl > R0, defined as f> = N(Rpl > R0)∕Ntot, is reported for different threshold sizes R0 as indicated by the horizontal blue lines. The histogram at the top shows the corresponding distribution of the final semimajor axis of the planetesimals. The vast majority of planetesimals are captured at distances ≲ 10 RJup and are strongly, if not completely, ablated.

5 Evolution of pebbles in the circum-planetary disk

The ablation of planetesimals in the CPD provides an important source of small dust grains from which satellites can accrete (this process is discussed in Sect. 6). We thus construct here a simple parametrization of the dusty component of a CPD fed by planetesimals ablation.

5.1 Dust supply

In the previous section, we have found that the cumulative ablated mass within the CPD could be roughly approximated by a cumulative Rayleigh distribution (Eq. (22)), corresponding to the following integral $$F(r;r_0)=1-e^{-r^2/(2r_0^2)}={\displaystyle {\int }_0^r\frac{x}{r_0^2}}{e}^{-x^2/(2r_0^2)}\text{d}x.$$(23)

Using physical quantities, the cumulative ablated mass in the CPD can be expressed as $$\frac{M(r)}{M_{\text{tot}}}=\frac{{\displaystyle {\int }_0^{r}2}\pi r\Sigma (r)\text{d}r}{M_{\text{tot}}},$$(24)

where M(r) is the total mass that has been deposited through ablation up to the distance r from the planet, Σ(r) is then the surface density distribution of the ablated mass and M_tot is the total mass that has been deposited in the CPD through ablation. By identification of the terms in Eqs. (23) and (24), we find $$\Sigma (r)=\frac{M_{\text{tot}}}{2\pi r_0^2}{e}^{-r^2/(2r_0^2)}.$$(25)

The above expression yields the distribution of material deposited by ablation as a function of the distancefrom the planet. Assuming that mass is globally deposited onto the CPD at a rate Ṁ₀, the source of dust provided by ablation can then be expressed as $${\dot{\Sigma }}_{\text{d}}=\frac{{\dot{M}}_0}{2\pi r_0^2}{e}^{-r^2/(2r_0^2)}.$$(26)

The value of r₀, which controls where the ablated mass is deposited, should depend on the structure of the CPD and the planetesimal size distribution in the feeding zone of the giant planet (see Sect. 4.2.1 and Fig. 5) which are rather uncertain. We thus treat it as a free parameter, with plausible values on the order of 10 planetary radii considering our results and the current radial extent of the satellite systems of Jupiter and Saturn.

Fig. 5 Cumulative mass distribution of material ablated off of planetesimals with initial sizes of 100 km (solid black line), 10 km (dashed black line) and 1 km (dotted black line) as a function of the radial distance from the planet. Smaller planetesimals are captured at wider distances and deliver material outward compared to the larger objects. The light blue curves correspond to cumulative Rayleigh distributions, $F(r;r_0)=1-\exp \left[-r^2/(2r_0^2)\right]$, that we use as approximate fits to the distribution of ablated material. The number to the right of each blue curve corresponds to the value of the parameter r0 of the Rayleigh distribution.

5.2 Quasi-steady-state analytical solution

The dust released through planetesimal ablation is not bound to remain where it has been deposited. The dust grains should rapidly grow and start to radially drift inward. The dust growth timescale is proportional to the orbital timescale and depends on the dust-to-gas mass ratio in the disk (e.g., Birnstiel et al. 2012): $${\tau }_{\text{growth}}\equiv \frac{s_{•}}{{\dot{s}}_{•}}\approx \frac{1}{{\epsilon }_{\text{d}}{\Omega }_{\text{K}}}\approx 1.65\times {10}^{-3}{\epsilon }_{\text{d}}^{-1}{\left(\frac{r}{10R_{\text{jup}}}\right)}^{3/2}\text{yr},$$(27)

where s_∙ is the radius of the dust grains, ε_d = Σ_d∕Σ_g is the dust-to-gas mass ratio, and the numerical estimate is for a Jupiter-mass planet. The radial drift velocity of dust, on the other hand, is a function of the aerodynamic coupling of the dust particles with the gas, which is expressed through their Stokes number, St = t_sΩ_K, and is given by $$v_r=-\frac{2\text{St}}{1+{\text{St}}^2}\eta v_{\text{K}},$$(28)

where v_K is the Keplerian orbital velocity around the planet and η is a measure of the gas pressure support of the CPD, $$\eta =-\frac{1}{2}\frac{\text{d}\ln P}{\text{d}\ln r}{h}_{\text{g}}^2,$$(29)

with P the pressure in the midplane of the CPD. This yields the radial drift timescale of dust, $$\begin{array}{ll}{\tau }_{\text{drift}}\equiv \left|\frac{r}{v_r}\right|\hfill & \approx \frac{1}{2\eta \text{St}}{\Omega }_{\text{K}}^{-1}\hfill \\ \hfill & \approx 13\left(\frac{{10}^{-2}}{\text{St}}\right){\left(\frac{0.06}{h_{\text{g}}}\right)}^2{\left(\frac{r}{10R_{\text{Jup}}}\right)}^{3/2}\text{yr}.\hfill \end{array}$$(30)

Given these short growth and radial drift timescales, it is likely that the dust distribution in the CPD rapidly reaches a quasi-steady-state regulated by the rate at which dust is supplied by ablation. This should hold as long as Ṁ₀ does not significantly vary over timescales ≫ τ_drift. The steady-state solution can be obtained by considering that the flux of pebbles through the disk at any given distance r must be equal to the total flux deposited exterior to r by ablation, translating to $$2\pi r{v}_r{\Sigma }_{\text{d}}={\displaystyle {\int }_r^{+\infty }\text{}}2\pi r{\dot{\Sigma }}_{\text{d}}\text{d}r={\dot{M}}_0{e}^{-r^2/(2r_0^2)}.$$(31)

When St ≪ 1, the drift velocity is v_r ≈ 2Stηv_K, which, using Eq. (31), yields the following expression of the dust-to-gas mass ratio in the CPD, $${\epsilon }_{\text{d}}\approx \frac{{\dot{M}}_0{e}^{-r^2/(2r_0^2)}}{4\pi r\text{St}\eta v_{\text{K}}{\Sigma }_{\text{g}}}.$$(32)

In Fig. 6, we show a comparison of the dust-to-gas mass ratio and pebble flux through the CPD of a Jupiter-mass planet, obtained by numerically solving the radial advection-diffusion equation of the dust surface density (e.g., Birnstiel et al. 2012), with the analytical solutions given by Eqs. (31) and (32) for different values of the Stokes number of the dust particles. We have used here a characteristic deposition radius r₀ = 10 R_Jup and a global supply rate Ṁ₀ = 1 M_⊕ Myr⁻¹ (equivalently 3 × 10⁻⁹ M_p yr⁻¹, for a Jovian-mass planet). The numerical integration results were taken after 5 × 10³ yr of evolution and effectively show a steady-state which agrees well with the analytical estimates. It is interesting to note that the flux of pebbles through the disk is indeed being regulated by the ablation supply rate and thus independent of the Stokes number of the particles (bottom panel of Fig. 6). The quasi-steady-state solution for the pebble flux and surface density equally applies to the case of the CPD of a Saturn-mass planet as the dust growth and drift timescales are equivalently short.

The top panel of Fig. 6 also shows the physical size of the dust grains in the CPD corresponding to the assumed Stokes number. It is important to note that, in the CPD, the mean free path of the gas molecules is much shorter than the typical values encountered in protoplanetary disks. Thus, pebbles are found in the Stokes drag regime or even the nonlinear drag regime closer to the planet (as indicated by the change of slope of the gray curves in the top panel of Fig. 6), as opposed to the Epstein regime which is the most relevant under PPD conditions (see Johansen et al. 2014, for a review of the different drag regimes). The relation between the Stokes number of the particles and their physical sizes in the relevant regimes are presented in Appendix B.

5.3 Pebbles Stokes number

Although the flux of pebbles through the CPD does not depend on their Stokes number in the steady-state, the pebble accretion rate onto a growing moon does (e.g., Lambrechts & Johansen 2012, and Sect. 6). Dust particles can initially grow by coagulation, but this growth is limited by either the rapid radial drift of the pebbles or fragmenting collisions at high relative velocities (e.g., Brauer et al. 2008). In the drift limited regime, pebbles can grow up to the point when their drift timescale becomes comparable to their growth timescale, which occurs at a Stokes number of (Birnstiel et al. 2012) $${\text{St}}_{\text{drift}}=2{\epsilon }_{\text{d}}{\eta }^{-1}.$$(33)

Using the expression of the dust-to-gas mass ratio in the CPD given by Eq. (32) yields $${\text{St}}_{\text{drift}}\approx {\left(\frac{{\dot{M}}_0{e}^{-r^2/(2r_0^2)}}{2\pi r{\eta }^2{v}_{\text{K}}{\Sigma }_{\text{g}}}\right)}^{1/2}.$$(34)

On the other hand, fragmentation occurs when the relative velocity between colliding particles becomes higher than a threshold fragmentation speed v_f. The corresponding Stokes number of the particles is then (Birnstiel et al. 2012) $${\text{St}}_{\text{frag}}=\frac{1}{3}{\alpha }^{-1}{\left(\frac{v_{\text{f}}}{c_{\text{g}}}\right)}^2,$$(35)

where α is the turbulent viscosity parameter.

The maximum Stokes number achieved in the fragmentation and drift limited regimes in the CPDs of Jupiter and Saturn⁴ are presented in Fig. 7 along with the corresponding particle sizes. For both Jupiter and Saturn, we assume a dust mass supply of Ṁ₀ = 3 × 10⁻⁹ M_p yr⁻¹ and a characteristic deposition radius of 10 planetary radii when estimating the drift limited Stokes number of the particles. As for fragmentation, we consider a turbulent viscosity parameter of α = 10⁻⁴, in agreement with our assumption of low-viscosity, passively irradiated CPDs, and a threshold fragmentation speed v_f = 1 m s⁻¹ (e.g., Blum & Wurm 2008). The resulting Stokes number and particle sizes are very similar for both planets. The drift limit typically yields Stokes numbers on the order of 10⁻², corresponding to roughly ~cm sized pebbles, except at distances larger than a few r₀, where there is little dust deposition, which promotes particle drift over their growth. Fragmentation limits the Stokes number of the pebbles to values that are typically lower than in the case of the drift regime (except, again, when far from r₀), with St in the range 10⁻³–10⁻² (corresponding to physical sizes of ~1–10 mm). The growth of pebbles should thus mainly be limited by fragmentation in the regions where the massive regular satellites of Jupiter and Saturn are found and the moons would therefore have accreted from pebbles with small Stokes numbers (St ≲ 10⁻²).

Fig. 6 Top: dust-to-gas mass ratio in the circum-planetary disk of a Jupiter-mass planet for pebbles Stokes number of 10−3 (black solid lines) and 10−2 (black dashed lines), assuming a characteristic radius of deposition of ablated material r0 = 10 RJup and a global supply rate Ṁ0 = 1 M⊕ Myr−1, obtained through numerical integration of the advection-diffusion equation of the evolution of the dust surface density after 5 × 103 yr of integration. The red dotted curves show the corresponding analytical steady-state solutionsgiven by Eq. (32). The gray curves show the physical size of pebbles as a function of distance from the planet corresponding to Stokes number of 10−3 (solid) and 10−2 (dashed). Bottom: radial flux of pebbles through the CPD for the two considered Stokes number of the pebbles (overlapping black curves), along with the steady-state solution where ${\dot{M}}_{\text{peb}}={\dot{M}}_0{e}^{-r^2/(2r_0^2)}$ (red dotted curve).

Fig. 7 Comparison of the Stokes number of the dust particles in the CPDs of Jupiter (top) and Saturn (bottom) in the fragmentation limited regime (black solid lines), assuming vf = 1 m s−1 and α = 10−4, and the drift limited regime (black dashed lines), assuming Ṁ0 = 10−9 Mp yr−1 and r0 = 10 RJup∕Sat. The solid and dashed gray curves show the physical sizes of the dust grains in the fragmentation and drift limited regimes, respectively.

6 Accretion of moons

Now that we have parametrized an ablation supplied circum-planetary disk, we turn to the discussion of the accretion of moons in such CPDs. We derive some characteristics of the accretion process in our proposed framework, with the aim of paving the way for future investigations using dedicated numerical simulations.

6.1 Pebble accretion

Once a protosatellite grows massive enough, it is able to catch the radially drifting dust grains in a process known as pebble accretion (see, e.g., Lambrechts & Johansen 2012, 2014; Ormel 2017; Johansen & Lambrechts 2017). Here we use the formalism of Liu & Ormel (2018) and Ormel & Liu (2018), where the pebble accretion rate Ṁ_PA is given by $${\dot{M}}_{\text{PA}}={\epsilon }_{\text{PA}}{\dot{M}}_{\text{peb}}.$$(36)

In the above expression, Ṁ_peb represents the mass flux of pebbles in the CPD. As we have shown in Sect. 5, the flux of pebbles should be at equilibrium with the supply of material through ablation of planetesimals, implying $${\dot{M}}_{\text{peb}}={\dot{M}}_0{e}^{-r^2/2r_0^2}.$$(37)

The term ε_PA in Eq. (36) is the pebble accretion efficiency (i.e., the fraction of pebbles drifting past an accreting seed that is captured), for which the aforementioned authors provide formulae fit to numerical simulations that we report in Appendix C.

Maps of the pebble accretion efficiency as a function of protosatellite mass and distance from the planet are shown in Fig. 8 in the case of the Jovian and Saturnian CPDs, assuming that pebbles have a constant Stokes number St = 5 × 10⁻³. At mass ratios comparable to that of the Galilean moons and Titan (q_s ~ 10⁻⁴), the pebble accretion efficiency is on the order of ~25%, and it is lessthan ~1% when q_s ≲ 10⁻⁶. This impliesthat only a small fraction of the available pebble flux is accreted by the moons over the course of their formation and that more than the current mass of the satellite systems worth of pebbles must have been processed through the CPDs of the giant planets. We discuss this and some implications in Sect. 7.2.

6.1.1 Onset of pebble accretion

A necessary condition for pebble accretion to operate is that the timescale during which a pebble interacts with an accreting seed is comparable to (or shorter than) its stopping time (see, e.g., Ormel 2017). According to Liu & Ormel (2018), this condition is satisfied when $$q_{\text{s}}\gg {\eta }^3\text{St}\approx 1.25\times {10}^{-9}\left(\frac{\text{St}}{5\times {10}^{-3}}\right){\left(\frac{h_{\text{g}}}{0.06}\right)}^6.$$(38)

Around Jupiter, this corresponds to an object with a size of ≈83 km. Around Saturn, assuming a thicker aspect ratio h_g ~ 0.08, this would correspond to a size of ≈98 km. The largest captured objects within the CPD could thus readily accrete pebbles, although they might initially grow by accreting other large captured planetesimals as these latter would dominate at distances beyond the characteristic ablation deposition radius r₀ (see Fig. 4).

6.1.2 Pebble isolation mass

An important aspect of pebble accretion is that it has an end. As an object grows more and more massive, it starts to significantly perturb the distribution of gas around it and carve an initially shallow gap. The outer edge of the gap acts as a barrier to the drift of pebbles, which are then unable to reach the accreting seed, hence terminating pebble accretion (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Bitsch et al. 2018). Using 3D hydrodynamical simulations, Lambrechts et al. (2014) derived the mass at which pebble accretion ends, the so-called isolation mass. Using their scaling relation, we find that, in the CPD, pebble isolation occurs at $$q_{\text{iso}}=6\times {10}^{-5}{\left(\frac{h_{\text{g}}}{0.05}\right)}^3.$$(39)

If the final mass of satellites is set by the pebble isolation mass, Eq. (39) implies that the aspect ratio of the Jovian CPD was h_g ≈ 0.055 at the time and place Ganymede finished its growth, whereas it should have been h_g ≈ 0.08 around Saturn where and when Titan has formed. The pebble isolation mass as a function of distance from the planet in the CPDs of Jupiter and Saturn is drawn as a yellow dashed line in Fig. 8.

6.2 Type-I migration

As a protosatellite grows more and more massive, it starts to migrate inward due to tidal interactions with the gaseous disk (so-called type-I migration). The type-I migration rate is classically expressed as $$\dot{r}=-k_{\text{mig}}{q}_{\text{s}}\frac{{\Sigma }_{\text{g}}{r}^2}{M_{\text{p}}}{h}_{\text{g}}^{-2}{v}_{\text{K}},$$(40)

where k_mig is a constant that depends on the temperature and surface density gradients of the disk. Since the type-I migration rate of a protosatellite is directly proportional to its mass, the induced orbital decay is a threat to the growth of moons (see Johansen et al. 2019, for an investigation of this issue in the case of planets).

A rough estimate of the mass a protosatellite can reach before migration dominates its evolution can be inferred by equating its growth timescale, τ_gr ≡ M_s∕Ṁ_PA, with its type-I migration timescale, τ_I ≡ r∕ṙ (e.g., Canup & Ward 2002). In the 3D pebble accretion regime, this migration limited mass is $$\begin{array}{lll}{q}_{\text{mig}}=\hfill & \hfill & \frac{0.39}{k_{\text{mig}}}\frac{h_{\text{g}}^2}{\eta h_{\text{peb}}}\frac{{\dot{M}}_{\text{peb}}{\Omega }_{\text{K}}}{{\Sigma }_{\text{g}}{r}^2}\approx 3.85\times {10}^{-6}{e}^{-r^2/(2r_0^2)}{k}_{\text{mig}}^{-1}\hfill \\ \hfill & \hfill & \times {\left(\frac{r}{r_0}\right)}^{5/7}\left(\frac{{\dot{M}}_0}{3\times {10}^{-9}{M}_{\text{p}}{\text{yr}}^{-1}}\right){\left(\frac{M_{\text{CPD}}}{1.5\times {10}^{-3}{M}_{\text{p}}}\right)}^{-1}\hfill \\ \hfill & \hfill & \times {\left(\frac{h_{\text{peb}}/h_{\text{g}}}{0.15}\right)}^{-1}{\left(\frac{{\tau }_{\text{acc}}}{5\text{\hspace{0.17em}Myr}}\right)}^{1/7}{\left(\frac{M_{\text{p}}}{M_{\text{Jup}}}\right)}^{-3/14}{\left(\frac{r_{\text{out}}/r_0}{16}\right)}^{1/2},\hfill \end{array}$$(41)

where the numerical estimate is valid in the flaring part of the CPD (i.e., at r > r_tran). In the above expression, h_peb is the aspect ratio of the pebble disk. For a given pebble Stokes number and dimensionless vertical diffusion coefficient δ (see, e.g., Johansen et al. 2014, and references therein), the pebble aspect ratio is a fraction of the gaseous disk aspect ratio such that (Dubrulle et al. 1995) $$\frac{h_{\text{peb}}}{h_{\text{g}}}=\sqrt{\frac{\delta }{\text{St}+\delta }},$$(42)

and we have used St = 5 × 10⁻³ and δ = 10⁻⁴ to estimate h_peb in Eq. (41).

The value of q_mig depends only weakly on the mass of the planet and is therefore of the same order of magnitude in the Jovian and Saturnian CPDs. It however depends linearly on the pebble flux through the disk and the mass of the CPD. For our choice of parameters, the migration mass lies one to two orders of magnitude below the mass of the major satellites of Jupiter and Saturn (q_s ~ 10⁻⁴). Figure 9 presents numerical integrations of the growth tracks of protosatellites accreting pebbles while migrating inward through type-I migration. At first, growth largely dominates the evolution of the objects, but when approaching the migration mass, the growth curves of the protosatellites become nearly horizontal and they grow only little until they reach the inner edge of the CPD.

It is difficult to imagine that the mass of the CPDs were orders of magnitude lower, as this would render the capture and ablation of large planetesimals inefficient. Given that the migration timescale of full-grown satellites is on the order of 10³ –10⁴ yr only, the survival of the satellites against migration would require that both their formation and the dissipation of the CPDs happened on even shorter timescales. This conclusion is not very satisfying since, on the one hand, the dissipation mechanism of a circum-planetary disk and the timescale on which it operates are yet to be constrained, and on the other hand, such short accretion timescales for the moons would imply strong accretional heating, which would be difficult to reconcile with the inferred internal state of Callisto and Titan (e.g., Barr & Canup 2008)⁵.

Instead,it seems that the migration of the protosatellites must have been stopped to allow for the formation and survival of the moons. The most likely mechanism that could have prevented the protosatellites to migrate all the way down to the surface of their parent planet is the truncation of the CPD by an inner magnetic cavity (Sasaki et al. 2010; Ogihara & Ida 2012). The sharp drop of density that would occur at the truncation of the disk is indeed known to act as a type-I migration trap (e.g., Masset et al. 2006; Liu et al. 2017), thereby allowing planets or satellites to pile-up in resonant chains at the inner edge of the disk (e.g., Ogihara et al. 2010). Moreover, the magnetic coupling of the giant planets with their CPD and the opening of a magnetic cavity are the most promising mechanisms to explain their sub-critical rotation rate (Takata & Stevenson 1996; Batygin 2018). Such a scenario would also place the formation and early dynamical evolution of giant planets’ moons in the same context as that of compact systems around stars (e.g., Liu et al. 2017; Izidoro et al. 2017, 2019; Lambrechts et al. 2019). We discuss that possibility below.

Fig. 8 Pebble accretion efficiency as a function of distance from the planet and mass of the protosatellite (qs = Ms∕Mp is the ratio of the mass of the seed to that of the planet) in a CPD surrounding Jupiter (top panel) or Saturn (bottom panel). Here, pebbles with a constant Stokes number (St = 5 × 10−3) were assumed. In both panels, yellow dots correspond to the current mass and orbital distances of the major satellites of Jupiter (Galilean moons), or Saturn (Titan). The yellow dashed line shows the pebble isolation mass as a function of the distancefrom the planet, which is the maximum mass a satellite can reach through pebble accretion.

Fig. 9 Growth tracks of pebble accreting and migrating protosatellites in the circum-jovian CPD assuming Ṁ0 = 10−6 M⊕ yr−1, r0 = 10 RJup, and pebbles with a constant Stokes number St = 5 × 10−3. Here, clearly, protosatellites cannot grow to masses comparable to that of the Galilean moons (marked as yellow dots) before reaching the inner edge of the disk.

6.3 Growth of satellites in resonant chains

Protosatellites would rapidly migrate in the dense environment of the CPD and reach the inner edge of the disk before being able to grow to masses comparable to that of the observed satellite systems or the pebble isolation mass (Eq. (41), Fig. 9). If the disk inner edge is truncated by a magnetosphere, a first migrating protosatellite would be trapped there, whereas a subsequently migrating object would be caught in a mean motion resonance with the inner protosatellite, resulting in the build-up of a resonant train of moons over time.

The protosatellites caught in such resonant chains would continue to grow until reaching the pebble isolation mass or up to the point when the giant planet has cleaned up its feeding zone and the influx of pebbles to the CPD stops. In any case, most of the mass of the moons would be accreted while in the resonant chain, at a fixed distance from the planet.

Given these considerations, we show in Fig. 10 growth maps of protosatellites at fixed radial distances in the Jovian and Saturnian CPDs, where we have assumed that pebbles have fixed Stokes numbers, with St = 5 × 10⁻³, that the characteristic deposition radius of ablated material is r₀ = 20 R_Jup∕Sat, and that Ṁ₀ = 3 × 10⁻⁹ M_p yr⁻¹. In every panel, the solid white line shows the time at which the pebble isolation mass is reached.

The top panels of Fig. 10 are for cases in which the flux of pebbles through the CPD is simply given by Eq. (37). Due to the smaller aspect ratio in the inner regions of the CPDs, pebble accretion is more efficient closer to the planet (see also Fig. 8) and protosatellites grow faster there. At distances beyond r₀, growth becomes very slow due to the exponential cut-off in the pebble flux and the larger aspect ratio of the CPDs. The formation ofcompact satellite systems should thus be a quite natural outcome.

In the bottom panels of Fig. 10, we have considered the effect of the snowline on the growth of the protosatellites, as this latter could have played an important role (e.g., Ronnet et al. 2017). We estimated that in the CPDs, the snowline lies at a temperature of ~220–225 K⁶, corresponding to distances of ~8 R_Jup and ~4 R_Sat from Jupiter and Saturn, respectively. We have not modeled the detailed processes, such as sublimation, condensation and diffusive redistribution, occuring at the snowline, but simply assumed that the flux of pebbles was halved inside of the snowline due to water evaporating from pebbles. Including this reduction of the pebble flux, the region just outside the snowline becomes a sweet spot for the formation of moons, where protosatellites grow the fastest. Thiscould have some important implications since a protosatellite reaching the pebble isolation mass would substantially reduce the flux of pebbles inside of its orbit, thus hindering the growth of inner moons. This effect is probably less so relevant in the case of cold CPDs, such as Saturn’s, as the snowline stands very close to the planet.

Applying such a scenario to the formation of the Galilean moons, we constructed pebble accretion growth tracks of Galilean analogs trapped in a chain of 2:1 mean motion resonances, including the reduction of the pebble flux due to outer moons catching some fraction of the pebbles. These are presented in Fig. 11. Here we have assumed that Ganymede analog is growing at 10 R_jup and that the aspect ratio of the CPD at this location is h_g = 0.055 (corresponding to assuming that τ_acc = 10 Myr according toEq. (7)), such that the pebble isolation mass equals the mass of Ganymede, and the snowline is located in between the orbits of Ganymede’s analog and Europa’s analog. We note that the formation of the moons closer to Jupiter than where they are currently found is compatible with the subsequent tidal expansion of their orbits, at least in the case of the three inner ones that are part of the Laplace resonance system (Yoder & Peale 1981; Lainey et al. 2009). Another mechanism, such as a dynamical instability (see Sect. 7.3), is nevertheless necessary to explain that Callisto is currently not in resonance with the other moons. Finally, we have assumed that r₀ = 12.5 R_Jup, that pebbles had a constant Stokes number St = 5 × 10⁻³, and stopped the integration when Callisto’s analog (the outermost moon) had reached the current mass of Callisto.

Importantly, the flux of pebbles inside of a moon having reached pebble isolation is not entirely cut-off because of material being constantly deposited by the ablation of planetesimals. The flux of pebbles at a distance r₁ interior to a protosatellite that has reached pebble isolation at a distance r₂ is then $${\dot{M}}_{\text{peb}}={\dot{M}}_0\left(e^{-r_1^2/(2r_0^2)}-e^{-r_2^2/(2r_0^2)}\right).$$(43)

This still corresponds to an abrupt reduction of the pebble flux, as can be observed in the growth tracks of Io’s and Europa’s analogs in Fig. 11, whose growth rates are substantially reduced after Ganymede’s analog has reached the pebble isolation mass. At the end of the integration, the mass of Europa’s analog is ≈30% higher than that of the actual satellite, and the mass of Io’s analog is only ≈8% lower than the current mass of Io.

It is nonetheless important to note that the growth tracks presented here are merely illustrative and meant to demonstrate that a pebble accretion scenario, similar in many regards to the scenario recently proposed for the formation of close-in super-Earths (e.g., Lambrechts et al. 2019; Izidoro et al. 2019), could lead to the formation of a system similar to that of the Galilean moons. As discussed in Sect. 7.3 however, a wide diversity of possible systems should be expected and their dynamical evolution must be explored with numerical simulations accounting for the mutual gravitational interactions among the satellites.

Fig. 10 Growth maps of protosatellites through pebble accretion as a function of distance from their parent planet and time, in thecase of Jovian (left panels) and Saturnian (right panels) circum-planetary disks. A fixed Stokes number, St = 5 × 10−3, has been assumed, as well as a characteristic deposition radius of ablated material r0 = 20 RJup∕Sat and a mass accretion rate Ṁ0 = 3 × 10−9 Mp yr−1. Top panels show cases where the flux of pebbles through the disk is ${\dot{M}}_{\text{peb}}={\dot{M}}_0\exp [-r^2/(2r_0^2)]$. The bottom panels include the effect of a snowline, assuming that the flux of pebbles is halved inside of the snowline due to the sublimation of water ice. The white line in each panel indicates the time needed to reach the pebble isolation mass at a given radial distance, whereas the white dots in the top and bottom left panels mark the time at which the mass of the Galilean moons is reached at their present location (but see Fig. 11 for actual growth tracks of Galilean-like moons). The time needed to reach the mass of Titan at its present location (~20 RSat) exceeds 1 My for the assumed parameters, and we argue in Sect. 7.3 that Titan may have formed from the merging of lower mass moons following a dynamical instability.

7 Discussion

7.1 Additional effects that can influence the capture and ablation of planetesimals

7.1.1 Aspect ratio of the CPD

Here we have assumed an average aspect ratio of the CPD of 0.06, appropriate for a cold disk around a slowly accreting Jovian planet (τ_acc = 5 Myr). If CPDs are able to develop during earlier phases of the formation of a giant planet, when gas is accreted at a higher rate, they could be much thicker, with an average aspect ratio of 0.3 (D’Angelo & Podolak 2015; Szulágyi 2017), which would influence the efficiency of capture and ablation of planetesimals. Fujita et al. (2013) show that the energy of planetesimals captured on bound orbits is mostly dissipated at their closest approach from the planet. It is therefore possible to define a distance r_capt from the planet at which a sufficient amount of energy ΔE to allow capture is lost. The energy dissipation due to gas drag during closest approach is ΔE ≈ a_dragr_CA (Tanigawa & Ohtsuki 2010; Fujita et al. 2013), where r_CA is the distance at closest approach. Considering the case of neglible inclination (i.e., coplanar case), using Eqs. (16), (17), (20) and (1), and setting r_capt = r_CA, we find $$\Delta E\approx \frac{3}{8}\frac{C_{\text{D}}{v}_{\text{rel}}^2}{{\rho }_{\text{pl}}{R}_{\text{pl}}}\frac{{\Sigma }_{\text{out}}{r}_{\text{out}}^{3/2}}{2\pi h_{\text{g}}}{r}_{\text{capt}}^{-3/2}.$$(44)

It follows that $r_{\text{capt}}\propto h_{\text{g}}^{-2/3}$ if v_rel does not depend on the distance to the planet (the dispersion regime of Fujita et al. 2013). On the other hand, if planetesimals approach the planet at the escape velocity at distance r_CA (the shear regime of Fujita et al. 2013), then $v_{\text{rel}}\propto r_{\text{CA}}^{-1/2}$ and $r_{\text{capt}}\propto h_{\text{g}}^{-5/2}$. Although the dependency on the aspect ratio of the CPD is not very strong, it is noteworthy that gas drag assisted capture of planetesimals is more difficult in puffed-up disks. Thus, even if CPDs are able to form early on during the accretion history of gas giant planets, the capture and ablation of planetesimals might be inefficient until their CPDs cool down enough.

Fig. 11 Example of pebble accretion growth tracks of Galilean moons analogs trapped in their mutual 2:1 mean motion resonances in a Jovian CPD with an accretion timescale onto Jupiter τacc = 10 Myr (see text for details). The growth tracks are shown as thick colored lines, where the color gives an indication of time. For further readability, black dots on the growth tracks show time intervals of 0.1 Myr and the larger dots indicate the end state of the system after 0.5 Myr of accretion. The thin gray line connects the system at the time Ganymede reaches the pebble isolation mass, which is after ≈0.3 Myr for our adopted value of Ṁ0 = 3 × 10−9 Mp yr−1. Tidal migration after the dissipation of the CPD can bring the orbits of the three inner satellites to their current location (Yoder & Peale 1981). Another mechanism, such as a dynamical instability, is however required to account for the fact that the outermost moon, Callisto, is presently outside of the resonance system.

7.1.2 Inclination of planetesimals

We have considered here the evolution of planetesimals with initially small inclinations, hence mostly experiencing coplanar interactions with the CPD of the giant planet. In the extreme opposite case, planetesimals would cross the CPD almost vertically, as discussed in Mosqueira et al. (2010). Vertical disk crossers would interact with the gas in the CPD over a distance ≈ 2H_g (Mosqueira et al. 2010) and the energy dissipated would then become ΔE ≈ 2a_dragH_g. Thus, comparing the amount of energy dissipated through gas drag in a coplanar encounter (ΔE_∥) at distance r_CA, with that dissipated in the case of a vertical encounter (ΔE_⊥) at the same distance, we have crudely that ΔE_⊥∕ΔE_∥~ H_g∕r_CA = h_g ≪ 1, where we have assumed that v_rel is of the same order of magnitude in both the coplanar and vertical encounters. The capture and ablation of highly inclined planetesimals should therefore be much less efficient than that of coplanar objects, which is consistentwith the results of Fujita et al. (2013). However, it was found in Ronnet et al. (2018) that planetesimals could be delivered to the feeding zone of a giant planet with rather small inclinations (i ≲10⁻¹) and effectively captured within the CPD.

7.1.3 Break-up of planetesimals

The high dynamic (or ram) pressures experienced by planetesimals crossing the CPD at high relative speeds might cause them to break-up, an effect that has been ignored in the present study due to its complex treatment (see, e.g., Register et al. 2017). Fragmented planetesimals would have shorter ablation timescales which could change the distribution of material within the CPD. Assuming planetesimals are held together by their own gravity, for a given dynamic pressure $P_{\text{dy}}=(1/2){\rho }_{\text{g}}{v}_{\text{rel}}^2$ acting on them, they are subject to break-up if their radius is smaller than (Pollack et al. 1979) $$R_{\text{break-up}}=\sqrt{\frac{5}{4\pi }\frac{P_{\text{dy}}}{G{\rho }_{\text{pl}}^2}}.$$(45)

It is then possible to derive the distance r_break-up at which planetesimals of a given size might break-up. In the case of a coplanar encounter, with $v_{\text{rel}}=(\text{}\sqrt{2}\pm 1)v_{\text{K}}$, and our Jovian CPD with a constant aspect ratio h_g = 0.06, we find $$r_{\text{break-up}}\approx 4.3{\left(\text{}\sqrt{2}\pm 1\right)}^{4/7}{\left(\frac{R_{\text{pl}}}{100\text{\hspace{0.17em}km}}\right)}^{-4/7}{R}_{\text{Jup}}.$$(46)

Therefore, fracturing and break-up of planetesimals is likely to occur within the CPDs of Jovian planets, especially for the smaller planetesimals on retrograde orbits. However, comparing with the curves showing the distribution of ablated material in Fig. 5, it appears that for every planetesimal size, break-up would occur in regions of the CPD where we find ablation to be an efficient mechanism, such that it is unlikely to significantly alter the distribution of dust in the disk. It should however be noted that break-up could play a more important role for the capture and ablation of highly inclined planetesimals which could otherwise be inefficient as discussed above.

7.2 Required amount of pebbles to grow moons and implications

Irrespective of the proposed mechanism of solids delivery to a CPD, it seems that pebble accretion is the most likely channel of growth for the moons. In the scenario proposed by Canup & Ward (2002, 2006), it was assumed that the small dust grains brought with the inflow of gas onto the CPD could grow up to large (≳km) sizes and avoid a rapid inward radial drift, but this should in fact not be the case (e.g., Shibaike et al. 2017, and our Sect. 5). Even if a dust trap would exist in CPDs, as recently proposed by Drążkowska & Szulágyi (2018), this would only allow for the formation of the satellite seeds which should subsequently grow by pebble accretion. On the other hand, as investigated in Sect. 4 (see also Estrada et al. 2009; Mosqueira et al. 2010; Fujita et al. 2013), the capture of planetesimals should mainly result in the delivery of small dust grains due to the strong ablation of the planetesimals crossing a CPD. Finally, large icy building blocks would be incompatible with the observed composition gradient among the Galilean moons (Ronnet et al. 2017; Dwyer et al. 2013), pointing toward the fact that pebbles would be better building blocks.

In a pebble accretion scenario, however, the growing moons would only accrete a fraction of the drifting pebbles. A majority of pebbles would simply drift past the orbit of the moons, until one or more moons reach the pebble isolation mass and block the pebble flow outside their partial gas gap. As shown by Fig. 8, satellites comparable in mass to the Galilean moons or Titan catch ≈25% of the pebble flux they encounter, and, for the integration shown in Fig. 11, 1.5 × 10⁻³ M_Jup worth of pebbles were processed through the CPD to form a Galilean-like system whereas the actual mass of the system is about 2 × 10⁻⁴ M_Jup.

The total amount of pebbles necessary to grow a satellite can be expressed as $$f_{\text{peb}}\equiv \frac{M_{\text{peb}}}{M_{\text{p}}}={\displaystyle {\int }_{q_{\text{seed}}}^{q_{\text{sat}}}{\epsilon }_{\text{PA}}^{-1}}\text{d}q.$$(47)

Table 1 reports the amount of pebbles needed to grow Ganymede and Titan, the most massive moons of Jupiter and Saturn, respectively. Since smaller aspect ratios yield higher pebble accretion efficiencies, the minimum amount of pebbles required to grow a given satellite through pebble accretion is obtained when the aspect ratio of the disk is such that the pebble isolation mass corresponds to the mass of the satellite (a smaller aspect ratio would imply an isolation mass which is smaller than the mass of the satellite). Even in this case, the overall (integrated) pebble accretion efficiency is less than 10% for both Ganymede (≈8%) and Titan (≈7%), meaning that an amount of pebbles equivalent to more than 10 times the mass of the moons is required to allow for their formation (see Table 1). Although this depends on the initial mass of the satellite seeds, considering seeds that are initially an order of magnitude more massive yields an increase in the integrated pebble accretion efficiency of only about 2%.

It is therefore difficult to envision that many generations of satellites could form around giant planets such as proposed by Canup & Ward (2006; and more recently by Cilibrasi et al. 2018), even if their CPD develops well before their completion. If giant planets were able to efficiently cool down and contract while still accreting a lot of gas, allowing for a CPD to form early on, this latter would likely be very puffed-up, with h_g ~ 0.3 (D’Angelo & Podolak 2015; Szulágyi 2017). Accretion would be extremely inefficient in such an environment, thus preventing the growth of satellites.

The low accretion efficiency of Galilean-like moons strengthens the argument that the inflow of material accreted by a giant planet toward the end of its formation can not provide enough solids to grow the observed satellite systems (Ronnet et al. 2018). Dust grains larger than ≳0.1 mm are efficiently filtered by the gap opened by a giant planet in its native disk (Paardekooper & Mellema 2006; Paardekooper 2007; Zhu et al. 2012; Weber et al. 2018). For typical dust size distributions, this would mean that a fraction of only 1–10% of the dust mass could be efficiently entrained with the gas accreted by the giant planet (Zhu et al. 2012). The formation of the satellites through the inflow of material accreted by the giant planet requires that a total mass M_tot = f_peb∕ε_d is processed through the CPD, where ε_d is the dust-to-gas mass ratio of the material accreted by the giant planet at the epoch of formation of the satellites. If ε_d ~ 10⁻³ (that is 10% of the typical solar value), M_tot becomes comparable to or even larger than the mass of the planet.

7.3 The putative role of dynamical instabilities

In the present study, we have argued that the delivery of solids to a CPD through the capture and ablation of planetesimals yields to a formation scenario of moons which resembles closely that recently proposed for the formation of super-Earths (Izidoro et al. 2019; Lambrechts et al. 2019) and other compact systems around M-dwarf stars, such as the planets orbiting around Trappist-1 (Ormel et al. 2017). Drawing further the analogy with super-Earth systems (we note that in terms of mass ratio, the distribution of super-Earths peaks at a few 10⁻⁵, which is fairly similar to the massive moons of the Solar System; Wu 2019; Liu et al. 2019), the resonant chains formed through the pile-up of migrating protosatellites would be prone to dynamical instabilities (Izidoro et al. 2017, 2019; Lambrechts et al. 2019). These latter can result in mergers and breaking of orbital resonances.

The contrast between Saturn’s satellite system, whose mass is largely concentrated in the massive moon Titan, and the Galilean system, composed by four roughly equal mass satellites, might be indicative of a more violent dynamical history around Saturn, where an initially, somewhat Galilean-like system, could have undergone an instability to finally merge into a single Titan while spawning mid-sized moons (Asphaug & Reufer 2013; Sekine & Genda 2012). Even in the case of the Galilean system, it seems difficult to avoid that the outermost satellite Callisto was once part of the resonant chain of moons. A dynamical re-arrangement of the system after the dissipation of the CPD could have occurred as well.

It thus appears important to further explore the formation of moons around giant planets through numerical simulations that can account for the gravitational interactions among the protosatellites and follow the long-term evolution of the formed systems (that is, after gas dissipation). Here we have provided a parametrization that is relevant to this exploration and could be used to set up such numerical simulations.

8 Summary

In the past decade, new paradigms have emerged regarding several key processes in the theory of planet formation which challenge our current understanding of the formation of the giant planets’ moons. Firstly, the formation of large objects (i.e., planetesimals or satellitesimals) requires some special conditions to occur (e.g., Johansen et al. 2014), which are not met in most parts of the circum-planetary disks where the moons should form (Shibaike et al. 2017; Drążkowska & Szulágyi 2018). Secondly, as a giant planet carves a gap in its natal protoplanetary disk, it blocks the flux of inwardly drifting pebbles (Lambrechts et al. 2014), thereby isolating itself from the main source of solids needed to grow its satellites (Ronnet et al. 2018). Lastly, the α-disk model, commonly used to constrain the formation conditions of the moons (e.g., Canup & Ward 2002), likely provides a poor description of the structure and evolution of a circum-planetary disk (Fujii et al. 2014, 2017; Mori et al. 2019).

Here we have thus considered the case of a low viscosity circum-planetary disk, where the main source of energy at the midplane comes from irradiation by the central planet. As capture of planetesimals on initially heliocentric orbits should be the main solids delivery mechanism to the CPD in the late stages of formation of a giant planet (e.g., Estrada et al. 2009; Tanigawa et al. 2012, 2014; Ronnet et al. 2018), the fate of the captured objects was investigated considering their ablation through frictional heating. We find that ablation is an efficient mechanism, resulting in the capture of large planetesimals in the CPD delivering primarily solids in the form of small dust grains. A fraction of planetesimals can however survive ablation and remain as large objects (≳ 10 km) in the circum-planetary disk, from which protosatellites can start to form.

The dust grains provided by the continuous ablation of planetesimals grow rapidly into radially drifting pebbles, whose flux through the circum-planetary disk is regulated by the rate at which fresh material is supplied through ablation. Protosatellites growth then mainly proceeds through the accretion of pebbles of millimeter-centimeter sizes.

We find that the overall pebble accretion efficiency (~10%), and the fast migration timescale of full-grown moons in CPDs (10³ –10⁴ yr), argue against a scenario where many generations of satellites are able to form around giant planets (Canup & Ward 2006; Cilibrasi et al. 2018). It seems more reasonable instead that the migration of protosatellites is stopped at the inner edge of the CPD due to, for example, its truncation by a magnetic cavity (Takata & Stevenson 1996; Sasaki et al. 2010; Batygin 2018). A similar conclusion was reached by Shibaike et al. (2019).

It follows that giant planets’ moons could have mainly grown while being caught in chains of orbital resonances, and their formation process would overall be very similar to that of super-Earths and other compact systems such as the Trappist-1 planets (Izidoro et al. 2017, 2019; Ormel et al. 2017; Lambrechts et al. 2019). The protosatellites would be able to grow until reaching the pebble isolation mass, which is comparable to the current masses of the Galilean satellites and Saturn’s moon Titan (when the accretion timescale of the parent planet is ≳ 5 Myr), or up to the point when the giant planet has cleaned its feeding zone and the influx of pebbles vanishes.

Given that resonant chains of planets are prone to dynamical instabilities (Izidoro et al. 2017, 2019; Lambrechts et al. 2019), the formation and subsequent evolution of moon systems should be investigated using N-body simulations to quantify the importance of such dynamical events in setting a system’s final architecture and, perhaps, understand the differences among the Jovian and Saturnian systems.

Acknowledgements

The authors thank the anonymous referee for their careful reading of the manuscript and suggestions for improvement. T.R. is thankful to Michiel Lambrechts and Beibei Liu for helpful discussions. The authors also wish to warmly thank Alessandro Morbidelli for stimulating discussions that initiated this project. T.R. and A.J. were supported by the European Research Council (ERC Consolidator Grant 724687-PLANETESYS). A.J. was further supported by the Knut and Alice Wallenberg Foundation (grant number 2012.0150) and the Swedish Research Council (grant 2018-04867).

Appendix A Surface temperature of planetesimals

The expression given in Eq. (18) for the surface temperature of planetesimals corresponds to the temperature at equilibrium between the different heating and cooling mechanisms, that is, when d T_pl∕dt = 0, under given conditions (CPD temperature, density, and planetesimal’s relative velocity with respect to the gas). From D’Angelo & Podolak (2015), the evolution of the surface temperature of a planetesimal evolves according to $$\begin{array}{lll}\frac{4}{3}\pi \left[R_{\text{pl}}^3-{(R_{\text{pl}}-\delta )}^3\right]{\rho }_{\text{pl}}{C}_p\frac{d{T}_{\text{pl}}}{dt}=\hfill & \hfill & \frac{\pi }{8}{C}_D{\rho }_{\text{g}}{R}_{\text{pl}}^2{v}_{\text{rel}}^3\hfill \\ \hfill & \hfill & +4\pi R_{\text{pl}}^2{\sigma }_{\text{sb}}(T_{\text{d}}^4-T_{\text{pl}}^4)\hfill \\ \hfill & \hfill & +L_w{\dot{m}}_{\text{abl}}\hfill \end{array}$$(A.1)

where $C_{\text{p}}=1600\text{\hspace{0.17em}J\hspace{0.17em}}{\text{kg}}^{\text{-1}}{\text{K}}^{-1}$ is the heat capacity of the planetesimal, and $$\delta =0.3\frac{K_{\text{pl}}}{{\sigma }_{\text{sb}}{T}_{\text{pl}}^3}$$(A.2)

is the thickness of the isothermal layer that is directly affected by the heating and cooling sources (D’Angelo & Podolak 2015), with K_pl = 3 W m^–1 K^–1 is the thermal conductivity of the planetesimal.

Fig. A.1 Comparison between the distribution of material ablated from the surface of 100 km sized planetesimals, considering ablation of water ice (blue curve), or silicates (red curve). For the silicates, we adopted the following parameters, Ls = 5 × 106 J kg–1, and log Pv = 13.176 − 24 605.0∕Tpl (which isappropriate for enstatite; Cameron 1985). The silicates are slightly more difficult to ablate than water ice, hence depositing material a bit closer to the planet. We note however that differentiated planetesimals would most likely be composed of hydrated minerals which are weaker than the silicates we have considered here (e.g., Fegley 2000), which would bring the silicate curve closer to that of water ice. After a hundred planet’s orbit, in the case of silicates, we find that the equivalent of ~20% of the mass initially present in the feeding zone of the giant planet has been delivered in the CPD through ablation, compared to ~25% in the case of water ice.

In Fig. A.2, we show the results of the integration of Eq. (A.1) for a planetesimal with R_pl = 100 km, using a 4th order Runge–Kutta method, starting from T_pl = 100 K, with the parameters T_d = 190 K, ρ_g = 4 × 10⁻² kg m⁻³, and v_rel = 33 km s^–1, corresponding to a close approach at 10 R_Jup from Jupiter in the retrograde direction. We compare the results with the equilibrium solution given by Eq. (18) for the same parameters. An equilibrium temperature is reached in ~50 ms only, which is much shorter than the typical perihelion passage timescale, on the order of a few hours (Eq. (21)). The solution given by Eq. (18) agrees perfectly with the equilibrium reached by solving Eq. (A.1), thus validating its use in our numerical integrations.

Fig. A.2 Comparison between the full integration of the evolution of the surface temperature of a planetesimal with a fourth orderRunge-Kutta method (black dots), and the solution given by Eq. (18), shown as a red dashed line.

Appendix B Physical size of the dust particles in the CPD

In the Jovian CPD, dust particles would transition from the Epstein to the Stokes drag regime at a size corresponding to (e.g., Weidenschilling 1977) $$s_{\text{Stokes}}=\frac{9\lambda }{4}\approx 0.1\text{\hspace{0.17em}mm\hspace{0.17em}}{\left(\frac{r}{10R_{\text{Jup}}}\right)}^{3-\xi /2},$$(B.1)

where $\lambda \propto {\rho }_{\text{g}}^{-1}$ is the mean free path of the gas molecules, and ξ is the power-law index of the temperature in the CPD (T_d ∝ r^−ξ), which changes at the transition radius. In this case, the size of the particles is given by $$s_{•}^{(\text{Stokes})}=\sqrt{\frac{9\text{St}{h}_{\text{g}}r\lambda {\rho }_{\text{g}}}{4{\rho }_{•}}}.$$(B.2)

This is valid as long as the Reynolds number of the flow around the dust particles, Re = 4s_•v_rel∕c_gλ, lies below unity. The further transition to the nonlinear regime, at Re > 1, can only be approximately derived by analytical means, as v_rel is not known a priori. However, since the Stokes number of the dust particles should remain small in the CPDs, their relative speed with respect to the gas in the azimuthal direction should vanish (see, e.g., Eq. (10) of Johansen et al. 2014), and v_rel should then be equivalent to the radial drift velocity of the particles, hence v_rel ≈ 2Stηv_K. The transition from the Stokes to the nonlinear drag regime would then occur at a size corresponding to $$s_{\text{NL}}\approx \frac{c_{\text{g}}\lambda }{8\text{St}\eta v_{\text{K}}}\approx 5.6\text{\hspace{0.17em}mm\hspace{0.17em}}\left(\frac{{10}^{-2}}{\text{St}}\right){\left(\frac{r}{10R_{\text{Jup}}}\right)}^3.$$(B.3)

In this nonlinear regime, the size of the particles is given by $$s_{•}^{(\text{NL})}\approx \frac{9^{5/8}}{2^{1/2}}\frac{{\text{St}}^{7/8}{c}_{\text{g}}^{3/8}{(\eta v_{\text{K}})}^{1/4}}{{({\rho }_{•}{\Omega }_{\text{K}})}^{5/8}{\rho }_{\text{g}}^{1/8}}.$$(B.4)

Appendix C Pebble accretion efficiency

In the present study, we have used the pebble accretion efficiency as provided in Liu & Ormel (2018) and Ormel & Liu (2018). It is expressed as a combination of the 2D and 3D regimes of pebble accretion, $${\epsilon }_{\text{PA}}={\left({\epsilon }_{\text{2D}}^{-2}+{\epsilon }_{\text{3D}}^{-2}\right)}^{-1/2}.$$(C.1)

In the above expression, ε_2D is the pebble accretion efficiency in the planar approximation (which is valid when the accretion radius is larger than the scale height of the pebbles), and reads (Liu & Ormel 2018) $${\epsilon }_{\text{2D}}=0.32\sqrt{\frac{q_{\text{s}}}{{\eta }^2\text{St}}\left(\frac{\Delta v}{v_{\text{K}}}\right)},$$(C.2)

with q_s the ratio of the protosatellite’s mass to that of the central planet, and Δv the relative velocity between the accreting seed and the pebbles, which, following Ormel & Liu (2018), can be expressed as $$\Delta v={\left[1+5.7\left(\frac{q_{\text{s}}\text{St}}{{\eta }^3}\right)\right]}^{-1}\eta v_{\text{K}}+0.52{(q_{\text{s}}\text{St})}^{1/3}{v}_{\text{K}}.$$(C.3)

This expression ignores effects related to the eccentricity of the protosatellites, which is not considered here. It also ignores the contribution of the gas accretion flow within the disk to the velocity of pebbles, which can become important for pebbles with small Stokes numbers (Johansen et al. 2019). However, expressing the velocity of the inward accretion flow as v_acc= M_p∕(2πrΣ_gτ_acc), we find it to be about an order of magnitude slower than the radial drift speed of pebbles with St = 5 × 10⁻³.

In the 3D regime, the pebble accretion efficiency is expressed as (Ormel & Liu 2018), $${\epsilon }_{\text{3D}}=0.39\frac{q_{\text{s}}}{\eta h_{\text{peb}}},$$(C.4)

where h_peb is the scale height of the pebbles, reading (Dubrulle et al. 1995) $$\frac{h_{\text{peb}}}{h_{\text{g}}}=\sqrt{\frac{\delta }{\text{St}+\delta }},$$(C.5)

where δ = 10⁻⁴ is the vertical diffusion coefficient. The small value of δ is motivated by the fact that the turbulence level in the CPD should be low due to the suppression of the magneto-rotational instability (Fujii et al. 2017). Low turbulence levels are also expected in the case of protoplanetary disks (e.g., Bai 2017), which seems to be supported by observations measuring dust settling to levels of about 10% of the gas scale height (Pinte et al. 2016).

References

All Tables

All Figures

	Fig. 1 Surface density profile of the protoplanetary disk obtained by a 2-dimensional hydrodynamic simulation of a Jupiter mass planet in a disk with a constant aspect ratio HPPD ∕R = 0.05.
In the text

Fig. 2 Example of orbits of objects that are captured in the prograde (top panels) and retrograde (bottom panels) directions in the planetocentric frame. Each dot corresponds to the position of the planetesimal at a given timestep and the arrows mark the direction at which the objects approach the planet. The left row shows the “typical” orbital evolution during the gas drag assisted capture of planetesimals whereas the right row corresponds to the orbital evolution following lower probability events where planetesimals are captured on orbits with a large pericenter and their semimajor axis decays slowly (see Fig. 4; see also, Suetsugu & Ohtsuki 2017). Insets on the top and bottom left panels are close-up views of the trajectories near the pericenter of the orbits with X, Y ∈ [−10;10]. The color of the dots indicates the ablation rate from the surface of the planetesimal. Ablation occurs mainly around the point of closest approach to the planets where the density of gas is the highest. Due to their larger relative velocities with respect to the gas, and hence larger frictional heating, planetesimals on retrograde orbits (bottom) experience more sustained ablation with ablation timescales as short as τabl ≡ m∕ṁabl ≲hr. The planetesimals captured in the retrograde directions both reached our imposed cut-off radius of 10 m (at this small size, the objects couple with the gas and end up rotating in the prograde direction which can be more clearly seen in the bottom right panel). However, far from the pericenter of the orbit, ablation rapidly becomes negligible for planetesimals orbiting on both prograde and retrograde orbits (τabl ≳ 107–108 h $\gg {\Omega }_{\text{K}}^{-1}$).

In the text

Fig. 3 Top: evolution of the radial distance from the planet of a planetesimal captured in the prograde direction (panel a), and the retrograde direction (panel b), with an initial size Rpl = 100 km, as a function of time (normalized by the inverse of the orbital angular frequency of the planet ΩK). Bottom: corresponding evolution of the mass of the planetesimal due to ablation in the case of prograde (panel c) and retrograde (panel d) capture. Ablation essentially occurs at closest approach from the planet due to the strong frictional heating in the dense inner parts of the CPD.

In the text

Fig. 4 Final radius versus final semimajor axis of planetesimals in the circum-planetary disk (i.e., at the time of circularization of the orbit according to our criteria, ã ≤ 0.05 RH ~ 40 RJup and ẽ ≤ 0.1). Each planetesimal started with an initial radius Rpl = 100 km at the beginning of the simulation. On the right, the fraction of objects with a size Rpl > R0, defined as f> = N(Rpl > R0)∕Ntot, is reported for different threshold sizes R0 as indicated by the horizontal blue lines. The histogram at the top shows the corresponding distribution of the final semimajor axis of the planetesimals. The vast majority of planetesimals are captured at distances ≲ 10 RJup and are strongly, if not completely, ablated.

In the text

Fig. 5 Cumulative mass distribution of material ablated off of planetesimals with initial sizes of 100 km (solid black line), 10 km (dashed black line) and 1 km (dotted black line) as a function of the radial distance from the planet. Smaller planetesimals are captured at wider distances and deliver material outward compared to the larger objects. The light blue curves correspond to cumulative Rayleigh distributions, $F(r;r_0)=1-\exp \left[-r^2/(2r_0^2)\right]$, that we use as approximate fits to the distribution of ablated material. The number to the right of each blue curve corresponds to the value of the parameter r0 of the Rayleigh distribution.

In the text

Fig. 6 Top: dust-to-gas mass ratio in the circum-planetary disk of a Jupiter-mass planet for pebbles Stokes number of 10−3 (black solid lines) and 10−2 (black dashed lines), assuming a characteristic radius of deposition of ablated material r0 = 10 RJup and a global supply rate Ṁ0 = 1 M⊕ Myr−1, obtained through numerical integration of the advection-diffusion equation of the evolution of the dust surface density after 5 × 103 yr of integration. The red dotted curves show the corresponding analytical steady-state solutionsgiven by Eq. (32). The gray curves show the physical size of pebbles as a function of distance from the planet corresponding to Stokes number of 10−3 (solid) and 10−2 (dashed). Bottom: radial flux of pebbles through the CPD for the two considered Stokes number of the pebbles (overlapping black curves), along with the steady-state solution where ${\dot{M}}_{\text{peb}}={\dot{M}}_0{e}^{-r^2/(2r_0^2)}$ (red dotted curve).

In the text

Fig. 7 Comparison of the Stokes number of the dust particles in the CPDs of Jupiter (top) and Saturn (bottom) in the fragmentation limited regime (black solid lines), assuming vf = 1 m s−1 and α = 10−4, and the drift limited regime (black dashed lines), assuming Ṁ0 = 10−9 Mp yr−1 and r0 = 10 RJup∕Sat. The solid and dashed gray curves show the physical sizes of the dust grains in the fragmentation and drift limited regimes, respectively.

In the text

Fig. 8 Pebble accretion efficiency as a function of distance from the planet and mass of the protosatellite (qs = Ms∕Mp is the ratio of the mass of the seed to that of the planet) in a CPD surrounding Jupiter (top panel) or Saturn (bottom panel). Here, pebbles with a constant Stokes number (St = 5 × 10−3) were assumed. In both panels, yellow dots correspond to the current mass and orbital distances of the major satellites of Jupiter (Galilean moons), or Saturn (Titan). The yellow dashed line shows the pebble isolation mass as a function of the distancefrom the planet, which is the maximum mass a satellite can reach through pebble accretion.

In the text

	Fig. 9 Growth tracks of pebble accreting and migrating protosatellites in the circum-jovian CPD assuming Ṁ0 = 10−6 M⊕ yr−1, r0 = 10 RJup, and pebbles with a constant Stokes number St = 5 × 10−3. Here, clearly, protosatellites cannot grow to masses comparable to that of the Galilean moons (marked as yellow dots) before reaching the inner edge of the disk.
In the text

Fig. 10 Growth maps of protosatellites through pebble accretion as a function of distance from their parent planet and time, in thecase of Jovian (left panels) and Saturnian (right panels) circum-planetary disks. A fixed Stokes number, St = 5 × 10−3, has been assumed, as well as a characteristic deposition radius of ablated material r0 = 20 RJup∕Sat and a mass accretion rate Ṁ0 = 3 × 10−9 Mp yr−1. Top panels show cases where the flux of pebbles through the disk is ${\dot{M}}_{\text{peb}}={\dot{M}}_0\exp [-r^2/(2r_0^2)]$. The bottom panels include the effect of a snowline, assuming that the flux of pebbles is halved inside of the snowline due to the sublimation of water ice. The white line in each panel indicates the time needed to reach the pebble isolation mass at a given radial distance, whereas the white dots in the top and bottom left panels mark the time at which the mass of the Galilean moons is reached at their present location (but see Fig. 11 for actual growth tracks of Galilean-like moons). The time needed to reach the mass of Titan at its present location (~20 RSat) exceeds 1 My for the assumed parameters, and we argue in Sect. 7.3 that Titan may have formed from the merging of lower mass moons following a dynamical instability.

In the text

Fig. 11 Example of pebble accretion growth tracks of Galilean moons analogs trapped in their mutual 2:1 mean motion resonances in a Jovian CPD with an accretion timescale onto Jupiter τacc = 10 Myr (see text for details). The growth tracks are shown as thick colored lines, where the color gives an indication of time. For further readability, black dots on the growth tracks show time intervals of 0.1 Myr and the larger dots indicate the end state of the system after 0.5 Myr of accretion. The thin gray line connects the system at the time Ganymede reaches the pebble isolation mass, which is after ≈0.3 Myr for our adopted value of Ṁ0 = 3 × 10−9 Mp yr−1. Tidal migration after the dissipation of the CPD can bring the orbits of the three inner satellites to their current location (Yoder & Peale 1981). Another mechanism, such as a dynamical instability, is however required to account for the fact that the outermost moon, Callisto, is presently outside of the resonance system.

In the text

Fig. A.1 Comparison between the distribution of material ablated from the surface of 100 km sized planetesimals, considering ablation of water ice (blue curve), or silicates (red curve). For the silicates, we adopted the following parameters, Ls = 5 × 106 J kg–1, and log Pv = 13.176 − 24 605.0∕Tpl (which isappropriate for enstatite; Cameron 1985). The silicates are slightly more difficult to ablate than water ice, hence depositing material a bit closer to the planet. We note however that differentiated planetesimals would most likely be composed of hydrated minerals which are weaker than the silicates we have considered here (e.g., Fegley 2000), which would bring the silicate curve closer to that of water ice. After a hundred planet’s orbit, in the case of silicates, we find that the equivalent of ~20% of the mass initially present in the feeding zone of the giant planet has been delivered in the CPD through ablation, compared to ~25% in the case of water ice.

In the text

	Fig. A.2 Comparison between the full integration of the evolution of the surface temperature of a planetesimal with a fourth orderRunge-Kutta method (black dots), and the solution given by Eq. (18), shown as a red dashed line.
In the text