The Galilean Satellites Formed Slowly from Pebbles

As we show later in Section 2.4, the pebble isolation mass (PIM) in the CJD is close to the actual mass of Ganymede. Therefore, we investigate two models. In the first model, four planetesimals are captured, grow, and migrate (Model A). In the second model, three planetesimals are captured, the third satellite reaches pebble isolation, and the fourth satellite forms out of the pebbles trapped at the gas pressure maximum associated with reaching PIM (Model B). The methods used in the two models are the same except for the treatment of PIM and the way to calculate the growth of Callisto. In both models, the captured planetesimals slowly accrete the particles drifting toward Jupiter, here referred to as pebbles (∼10 cm). The interiors of the satellites are mainly heated by the radiogenic decay of ²⁶Al included in the accreted pebbles. Figure 1 represents the two models of our slow pebble accretion scenario.

Zoom In Zoom Out Reset image size

At first, we model the gas accretion rate as ${\dot{M}}_{{\rm{g}}}=0.2\exp (-(t-{t}_{\mathrm{gap}})/{t}_{\mathrm{dep}}))\,{M}_{{\rm{J}}}\,{\mathrm{Myr}}^{-1}$, where t, t_gap, t_dep, and M_J are the time after the formation of CAIs, gap opening time, gas depletion timescale of the CSD (we assume t_dep = 3 Myr), and current Jupiter mass M_J = 1.90 × 10³⁰ g, respectively. We consider a 1D viscous accretion CJD. We assume that the gas mass flux is uniform in the CJD and that it is equal to both the inflow mass flux to the disk and the gas accretion rate to Jupiter, and the mass of Jupiter grows by ${\dot{M}}_{{\rm{g}}}$ from 0.4 M_J at t = t_gap to 1.0 M_J at the end of the calculation. We fix the position of the edge of the cavity r_cav at the current position of Io. On the distribution of the disk temperature in our model, the position of the snow line r_snow is at Ganymede's orbit when the gap opens and moves to Europa's orbit at the end of the disk evolution.

We calculate the evolution of dust particles in the CJD using a 1D single-size analytical formula (Shibaike et al. 2017). They grow to pebbles and then drift to Jupiter because the gas disk rotates with sub-Kepler velocity, which is slower than the rotating velocity of the pebbles, so that they lose their angular momentum. We assume the dust inflow mass flux as $x{\dot{M}}_{{\rm{g}}}$, where x, the dust-to-gas accretion ratio, is a constant parameter with x = 0.0026. The pebble mass flux M_p is equal to $x{\dot{M}}_{{\rm{g}}}$ at the outer edge of the disk, but reduces inward because of filtering by satellites and the evaporation of ice. There are only rocky pebbles inside the snow line, and we assume that the pebble mass flux is halved. We also assume that the fragmentation of pebbles occurs when their collisional velocity becomes faster than 5 or 50 m s⁻¹ for inside or outside the snow line, respectively (Wada et al. 2009, 2013). We calculate the disk midplane temperature as ${T}_{{\rm{d}}}={(3{{GM}}_{\mathrm{cp}}{\dot{M}}_{{\rm{g}}}/(8\pi {\sigma }_{\mathrm{SB}}{r}^{3}))}^{1/4}g$, where g = (3/8τ + 1/(4.8τ))^1/4, and σ_SB and G are the Stefan–Boltzmann constant and the gravitational constant, respectively (Nakamoto & Nakagawa 1994). The Rosseland mean opacity is τ = κΣ_g, where the dust opacity is κ = 450r_gg for T_d ≥ 160 K and $\kappa =450{({T}_{{\rm{d}}}/160{\rm{K}})}^{2}{r}_{\mathrm{gg}}$ for T_d < 160 K, where r_gg is the ratio of the surface density of grains affecting the disk temperature to gas.

We calculate the mass of the growing seeds as

$$\begin{eqnarray}&&{M}_{{\rm{s}}}(t)={\int }_{{t}_{\mathrm{cap}}}^{t}{\dot{M}}_{{\rm{p}}}{P}_{\mathrm{eff}}\,{dt},\end{eqnarray} \tag{ 1 }$$

where t_cap and P_eff are the capture time of the seeds and their pebble accretion efficiencies. The pebble accretion efficiency depends on the mass of the seed and the Stokes number of the pebbles around it (Ormel & Liu 2018). For Callisto in Model B, which is located at a pressure bump, we need to modify Equation (1). See Section 2.4.

We also calculate the migration of the seeds by aerodynamic drag and Type I migration, which includes both inward and outward migration (Adachi et al. 1976; Weidenschilling 1977; Paardekooper et al. 2011; Ogihara et al. 2015). We finally consider the capture into 2:1 or 3:2 MMRs (Ogihara & Kobayashi 2013). The orbit of Callisto is fixed in Model B at the gas pressure maximum created by Ganymede.

We calculate the surface temperature, T_s(R_s), with the radii of the seeds, R_s, from the equilibrium of the emission and accretion heat of pebbles. We also calculate the cumulative heat of ²⁶Al decay, ΔT_fin(R), from when the pebbles (the ice mass fraction is 0.5) including ²⁶Al have been accreted on the seeds' surface to the end of the formation, where R is the distance from the center of the seeds. We can then estimate the final (i.e., maximum) internal temperature of the seeds at the point R by ${T}_{\mathrm{fin}}(R)={T}_{{\rm{s}}}(R)+{\rm{\Delta }}{T}_{\mathrm{fin}}(R)$ (Barr & Canup 2008). We do not include thermal diffusion, solid-state convection, and latent heat inside the seeds.

Below, we show the details of the methods.

After the gap formation, the gas accretion rate to Jupiter becomes much lower. In this work, we assume that the gas accretion rate is

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{g}}}={\dot{M}}_{{\rm{g}},\mathrm{gap}}\exp \left(-\displaystyle \frac{t-{t}_{\mathrm{gap}}}{{t}_{\mathrm{dep}}}\right)\ \ \ (t\gt {t}_{\mathrm{gap}}).\end{eqnarray} \tag{ 2 }$$

We also assume that the initial gas accretion rate and the gas depletion timescale are ${\dot{M}}_{{\rm{g}},\mathrm{gap}}=0.2\,{M}_{{\rm{J}}}\,{\mathrm{Myr}}^{-1}$ and t_dep = 3 × 10⁶ yr, where M_J is the current Jupiter mass. We set t = 0 as the time of the formation of CAIs. The mass of Jupiter grows by ${\dot{M}}_{{\rm{g}}}$ from 0.4 M_J at t = t_gap to 1.0 M_J at the end of the calculation. The accretion rate is equal to the inflow mass flux from the CSD to CJD because we assume that the flow of gas (and pebbles) in the whole region of the disk is semi-steady and all of the gas flow into the disk will be accreted by Jupiter eventually.

A cavity of the gas disk can open around Jupiter by the magnetic field of the planet. The time at which the inner cavity opens depends on the strength of the magnetic field. The position of the disk inner edge can be estimated from the balance of gas accretion and magnetic stress (Lovelace et al. 2011; Liu et al. 2017). If Jupiter has a dipole magnetic field, it is given by

$$\begin{eqnarray}\begin{array}{rcl}{r}_{\mathrm{cav}} & = & {\left(\displaystyle \frac{{B}_{\mathrm{cp}}^{4}{R}_{\mathrm{cp}}^{12}}{4{{GM}}_{\mathrm{cp}}{\dot{M}}_{{\rm{g}}}^{2}}\right)}^{1/7}\\ & = & \ 1.07{\left(\displaystyle \frac{{B}_{\mathrm{cp}}}{40\mathrm{Gauss}}\right)}^{4/7}{\left(\displaystyle \frac{{R}_{\mathrm{cp}}}{{R}_{{\rm{J}}}}\right)}^{12/7}\\ & & \times \ {\left(\displaystyle \frac{{M}_{\mathrm{cp}}}{0.4{M}_{{\rm{J}}}}\right)}^{-1/7}{\left(\displaystyle \frac{{\dot{M}}_{{\rm{g}}}^{2}}{0.2{M}_{{\rm{J}}}{\mathrm{Myr}}^{-1}}\right)}^{-2/7}\,[{R}_{{\rm{J}}}],\end{array}\end{eqnarray} \tag{ 3 }$$

valid for r_cav < r_co, where r_co is the corotation radius (see Section 4.3). In the equation, B_cp, R_cp, G, and M_cp are the strength of the magnetic field of the central planet, radius of the central planet, gravitational constant, and mass of the central planet, respectively. Current Jupiter has a magnetic field, and its strength on the surface of the equational region of the planet is 4.2 G (Connerney 1993). Previous work, however, argued that the magnetic field was once stronger than the current one (Stevenson et al. 1983; Sánchez-Lavega 2004; Christensen et al. 2009). Therefore, we consider that the magnetic field is ≈40 G. In this case, the disk inner cavity and gap open at almost the same time (substituting the Jupiter radius R_J = 7.15 × 10⁹ cm at r_cav for Equation (3)). On the other hand, if the strength of the magnetic field is the same as the current one, the cavity only opens 14 Myr later than the gap opening. In this case, any satellite would be consumed by proto-Jupiter (t < 14 Myr), and too little material would remain to form the Galilean satellites after gap opening. We summarize the evolution of ${\dot{M}}_{{\rm{g}}}$ in Figure 2. In our model, we assume that r_cav is fixed at the current Io's orbit, 5.89 R_J, for simplicity.

Zoom In Zoom Out Reset image size

Based on the above gas accretion model, we calculate the evolution of the 1D CJD. The gas surface density of the viscous accretion disk is

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{g}}}=\displaystyle \frac{{\dot{M}}_{{\rm{g}}}{{\rm{\Omega }}}_{{\rm{K}}}}{3\pi \alpha {c}_{{\rm{s}}}^{2}}.\end{eqnarray} \tag{ 4 }$$

We assumed that the strength of turbulence is α = 10⁻⁴ because MRI is suppressed in the CJD (Fujii et al. 2014). The top left panel of Figure 3 represents the evolution of the gas surface density. The gas surface density becomes small as the gas accretion rate decreases.

Zoom In Zoom Out Reset image size

We assume that the CJD is viscously heated. The gas temperature in the midplane of the viscous accretion disk is given by (Nakamoto & Nakagawa 1994)

$$\begin{eqnarray}&&{T}_{{\rm{d}}}={\left(\displaystyle \frac{3{{GM}}_{\mathrm{cp}}{\dot{M}}_{{\rm{g}}}}{8\pi {\sigma }_{\mathrm{SB}}{r}^{3}}\right)}^{1/4}g,\end{eqnarray} \tag{ 5 }$$

where σ_SB is the Stefan–Boltzmann constant and

$$\begin{eqnarray}&&g={\left(\displaystyle \frac{3}{8}\tau +\displaystyle \frac{1}{4.8\tau }\right)}^{1/4}\end{eqnarray} \tag{ 6 }$$

is a function of the Rosseland mean optical depth τ = κΣ_g. In principle, the opacity κ depends on the size distribution of solid particles. However, the size distribution cannot be predicted from the simple dust evolution model as employed in this study. Therefore, we just assume r_gg, the ratio of the surface density of grains that affect the temperature to the gas surface density. Then, the opacity can be assumed as

$$\begin{eqnarray}\kappa =\left\{\begin{array}{ll}450{r}_{\mathrm{gg}} & {T}_{{\rm{d}}}\geqslant 160\,{\rm{K}}\\ 450{\left({T}_{{\rm{d}}}/160{\rm{K}}\right)}^{2}{r}_{\mathrm{gg}} & {T}_{{\rm{d}}}\lt 160\,{\rm{K}}.\end{array}\right.\end{eqnarray} \tag{ 7 }$$

The bottom left panel of Figure 3 represents the evolution of the temperature when r_gg = 1.7 × 10⁻⁷. We choose this value to get the thermal condition that the snow line is just inside current Ganymede's orbit. The gradients of curves change at the snow line, where the temperature is 160 K. The Rosseland mean optical depth τ becomes the smallest (∼1) slightly outside the snow line. The temperature decreases as the gas accretion rate reduces. The temperature depends on this ratio r_gg, and we discuss the dependence of the results on this parameter in detail in Section 4.

We calculate the distributions of the Stokes number and surface density of the drifting pebbles in the CJD. We only consider semi-steady conditions of pebbles because the evolution timescale of the pebbles is much shorter than those of the disk and satellites. In this case, the pebble mass flux in the CJD ${\dot{M}}_{{\rm{p}}}$ does not depend on the distance from Jupiter r. It is also equal to the mass flux of dust particles supplied from the parent CSD to the CJD,

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{p}}}=x{\dot{M}}_{{\rm{g}}},\end{eqnarray} \tag{ 8 }$$

where x is the ratio of the dust to gas accretion rates. We treat this ratio as a parameter and assume that the ratio does not depend on time, for simplicity. However, the mass flux of pebbles drifting inside the snow line is smaller than that outside because pebbles lose their H₂O ice inside the snow line. We assume that the rock mass fraction of pebbles outside the snow line is m_r = 0.5, so that ${\dot{M}}_{{\rm{p}}}$ inside the snow line becomes half of that of outside. We also consider the filtering effect by outer satellites.

Under these assumptions, we first calculate the Stokes number of the pebbles. When the Stokes number is determined by radial drift, it can be calculated by the following equation (modified version of Equation (15) of Shibaike et al. 2017):

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{St}}_{{\rm{p}}} & = & 0.23{\left(\displaystyle \frac{2}{3+2p+q}\right)}^{4/5}{\left(\displaystyle \frac{10}{18-39q}\right)}^{2/5}\\ & & \times \ {\left(\displaystyle \frac{{\dot{M}}_{{\rm{p}}}/{\dot{M}}_{{\rm{g}}}}{0.003}\right)}^{2/5}{\left(\displaystyle \frac{\alpha }{{10}^{-4}}\right)}^{1/5}\\ & & \times \ {\left(\displaystyle \frac{{T}_{{\rm{d}}}}{160{\rm{K}}}\right)}^{-2/5}{\left(\displaystyle \frac{{M}_{\mathrm{cp}}}{1{M}_{{\rm{J}}}}\right)}^{2/5}{\left(\displaystyle \frac{r}{10{R}_{{\rm{J}}}}\right)}^{-2/5},\end{array}\end{eqnarray} \tag{ 9 }$$

where p and q are the r exponents of the gas surface density and temperature (i.e., Σ_g ∝ r^−p and T_d ∝ r^−q).

However, fragmentation occurs if the collision velocity, in other words, the pebble-to-pebble relative velocity, is too fast. This relative velocity is

$$\begin{eqnarray}&&{v}_{\mathrm{pp}}=\sqrt{{\left({v}_{{\rm{r}}}/2\right)}^{2}+{v}_{{\rm{t}}}^{2}},\end{eqnarray} \tag{ 10 }$$

where v_r and v_t are the radial drift velocity of the pebbles and the relative velocity induced by turbulence, respectively. These two velocities are (Adachi et al. 1976; Weidenschilling 1977; Ormel & Cuzzi 2007)

$$\begin{eqnarray}&&{v}_{{\rm{r}}}=-2\displaystyle \frac{{\mathrm{St}}_{{\rm{p}}}}{{{\mathrm{St}}_{{\rm{p}}}}^{2}+1}\eta {v}_{{\rm{K}}},\end{eqnarray} \tag{ 11 }$$

where v_K = rΩ_K is the Kepler velocity and

$$\begin{eqnarray}&&{v}_{{\rm{t}}}=\sqrt{3\alpha }{c}_{{\rm{s}}}{{\mathrm{St}}_{{\rm{p}}}}^{1/2}.\end{eqnarray} \tag{ 12 }$$

The ratio of the pressure gradient force to the gravity of Jupiter is

$$\begin{eqnarray}&&\eta =-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{H}_{{\rm{g}}}}{r}\right)}^{2}\displaystyle \frac{\partial \mathrm{ln}{\rho }_{{\rm{g}}}{c}_{{\rm{s}}}^{2}}{\partial \mathrm{ln}r},\end{eqnarray} \tag{ 13 }$$

where ${\rho }_{{\rm{g}}}={{\rm{\Sigma }}}_{{\rm{g}}}/(\sqrt{2\pi }{H}_{{\rm{g}}})$ is the gas density at the midplane. If the Stokes number is limited by their fragmentation, it is

$$\begin{eqnarray}&&{\mathrm{St}}_{{\rm{p}}}=\displaystyle \frac{-3\alpha {c}_{s}^{2}+\sqrt{9{\alpha }^{2}{c}_{s}^{4}+4{\eta }^{2}{r}^{2}{{\rm{\Omega }}}_{{\rm{K}}}^{2}{v}_{\mathrm{cr}}^{2}}}{2{\eta }^{2}{r}^{2}{{\rm{\Omega }}}_{{\rm{K}}}^{2}},\end{eqnarray} \tag{ 14 }$$

where v_cr is the critical fragmentation speed (Okuzumi et al. 2016). This equation can be derived by substituting Equations (11), (12), and v_pp = v_cr for Equation (10).

Finally, the pebble surface density follows from the continuity equation,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{p}}}=\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{2\pi {{rv}}_{{\rm{r}}}}.\end{eqnarray} \tag{ 15 }$$

The right panels of Figure 3 represent the evolution of the Stokes number and the surface density of the dust particles (i.e., pebbles). The stair around 10 R_J is consistent with the position of the snow line. The snow line migrates inward because the temperature becomes lower as the gas accretion rate decreases. We also find that the Stokes number inside the snow line is smaller than that outside. Outside the snow line, the Stokes number is determined by drift (Equation (9)). Inside the snow line, the Stokes number is determined by fragmentation (Equation (14)) because rocky particles are more fragile than icy ones (Wada et al. 2009, 2013). We assume that the critical fragmentation speeds of rocky and icy pebbles are v_cr = 5 and 50 m s⁻¹, respectively. In the right panels of Figure 3 the minor stairs reflect the accretion of pebbles by the planets. As time goes on and the satellites grow larger, the pebble accretion efficiency increases, resulting in larger jumps (see Section 2.4). Although we do not consider the inner cavity in this calculation, the rocky pebbles should flow onto the planet with the gas because they are small enough to couple to the magnetospheric accretion flow of gas. Their dynamical timescale should be about the freefall timescale, and the stopping time of pebbles is much smaller than it because the top panel of Figure 3 shows that the Stokes number is about 0.02 around the cavity (5.89 R_J).

According to recent N-body simulations, the pebble accretion efficiency is well fitted by (Liu & Ormel 2018)

$$\begin{eqnarray}&&{P}_{\mathrm{eff}}={\left\{{\left(0.32\sqrt{\displaystyle \frac{{\mu }_{s}{\rm{\Delta }}v/{v}_{{\rm{K}}}}{{\mathrm{St}}_{{\rm{p}}}{\eta }^{2}}}\right)}^{-2}+{\left(0.39\displaystyle \frac{{\mu }_{{\rm{s}}}}{\eta {h}_{{\rm{p}}}}\right)}^{-2}\right\}}^{-1/2},\end{eqnarray} \tag{ 16 }$$

where μ_s = M_s/M_cp and h_p = H_p/r are the satellite-to-central planet mass ratio and the pebble aspect ratio, respectively. This equation combines two regimes of pebble accretion, 2D (the first term) and 3D (the second term) limits. If the pebble accretion radius is larger than the pebble scale height H_p, the first term is dominant. The pebble scale height can be derived analytically from the balance of their vertical sedimentation and diffusion (Youdin & Lithwick 2007),

$$\begin{eqnarray}&&{H}_{{\rm{p}}}={H}_{{\rm{g}}}{\left(1+\displaystyle \frac{{\mathrm{St}}_{{\rm{p}}}}{\alpha }\displaystyle \frac{1+2{\mathrm{St}}_{{\rm{p}}}}{1+{\mathrm{St}}_{{\rm{p}}}}\right)}^{-1/2},\end{eqnarray} \tag{ 17 }$$

where H_g = c_s/Ω_K is the gas scale height. The expression in the 2D limit depends on the approach velocity of the pebbles Δv, which is given by the Keplerian shear in the disk or else the disk headwind. The approach velocity is

$$\begin{eqnarray}&&{\rm{\Delta }}v{/v}_{{\rm{K}}}=0.52{\left({\mu }_{{\rm{s}}}{\mathrm{St}}_{{\rm{p}}}\right)}^{1/3}+\eta {\left\{1+5.7\left(\displaystyle \frac{{\mu }_{{\rm{s}}}}{{\eta }^{3}/{\mathrm{St}}_{{\rm{p}}}}\right)\right\}}^{-1},\end{eqnarray} \tag{ 18 }$$

where the first and second terms represent the shear and the headwind limits, respectively (Ormel & Liu 2018).

The pebble mass flux within an orbit of the seed is smaller than without because a fraction of drifted pebbles are captured by the embryo,

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{p}},\mathrm{in}}=(1-{P}_{\mathrm{eff}}){\dot{M}}_{{\rm{p}},\mathrm{out}},\end{eqnarray} \tag{ 19 }$$

where ${\dot{M}}_{{\rm{p}},\mathrm{in}}$ and ${\dot{M}}_{{\rm{p}},\mathrm{out}}$ are the fluxes inside and outside the embryo, respectively. The growth timescale of the seeds is t_grow = M_s/(dM_s/dt).

Efficient pebble accretion occurs when the satellite seed mass exceeds a critical mass (Ormel 2017),

$$\begin{eqnarray}&&{M}_{* }=\displaystyle \frac{{v}_{\mathrm{hw}}^{3}{t}_{\mathrm{stop}}}{8G}=\displaystyle \frac{1}{8}{\eta }^{3}{\mathrm{St}}_{{\rm{p}}}{M}_{\mathrm{cp}}\end{eqnarray} \tag{ 20 }$$

(settling regime). In this work, we assume that all seeds always accrete pebbles effectively. In our CJD model, we found that the gas aspect ratio h_g is about 0.1 in the whole disk regions and so η ∼ h_g² ∼ 10⁻². The Stokes number of the pebbles is St_p ∼ 0.1 (see Figure 3). The critical mass is then M_* ∼ 10⁻⁸M_cp, which is about 10 times smaller than the initial mass M_ini = 3 × 10²³g.

Pebble accretion stops if a gap structure forms around the seed (Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018; Johansen et al. 2019). The critical mass, known as the PIM, is M_iso ∼ h_g³M_cp, which is close to the largest Galilean satellite—Ganymede—in our CJD model. Because of this similarity, we consider two scenarios: (i) Ganymede did not reach the PIM (Model A), and (ii) Ganymede reached pebble isolation (Model B). In the latter case, we therefore define the PIM of Ganymede as its actual current mass, M_G = 1.48 × 10²⁶ g (see also Section 4.10). We also assume that if the mass of a satellite reaches the PIM, the growth of the satellite stops and the pebble mass flux inside the satellite (${\dot{M}}_{{\rm{p}},\mathrm{in}}$ in Equation (19)) becomes zero immediately.

At the gas pressure maximum, another satellite may form by pebble accretion because drifting pebbles pile up there. The growth rate of the satellite at the gas pressure maximum is given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{{\rm{s}}}}{{dt}}=\min ({R}_{\mathrm{col}}{{\rm{\Sigma }}}_{{\rm{p}}}{{\rm{\Omega }}}_{{\rm{K}}}{r}^{2},{\dot{M}}_{{\rm{p}}}),\end{eqnarray} \tag{ 21 }$$

where R_col = 2πηP_eff is the dimensionless accretion rate. Note that R_col does not contain an explicit η-dependence. The approach velocity of pebbles in P_eff, Δv/v_K, is then calculated by Equation (18) without the headwind term (the second term on the right-hand side). Kanagawa et al. (2018) argued that, in cases of planets in protoplanetary disks (PPDs), when a planet reaches its PIM, the dust-to-gas surface density ratio outside the gap formed by the planet becomes about unity and the dust-rich region (ring) widens outward. Therefore, we assume Σ_p = Σ_g here, but we find that the ${\dot{M}}_{{\rm{p}}}$ term in Equation (21) nonetheless limits the growth of the satellite: as soon as pebbles arrive at the pressure bump, they are accreted.

We calculate the internal temperature of satellites to estimate the levels of differentiation. The surface temperature of satellites with radius R_s, can be estimated by the following equation (Barr & Canup 2008). The temperature T_s(R_s) is

$$\begin{eqnarray}&&{\rho }_{{\rm{s}}}{C}_{{\rm{p}}}({T}_{{\rm{s}}}({R}_{{\rm{s}}})-{T}_{{\rm{d}}})\displaystyle \frac{{{dR}}_{{\rm{s}}}}{{dt}}=\displaystyle \frac{1}{2}\displaystyle \frac{{\dot{M}}_{{\rm{s}}}{u}_{{\rm{i}}}^{2}}{4\pi {R}_{{\rm{s}}}^{2}}-{\sigma }_{\mathrm{SB}}({T}_{{\rm{s}}}{\left({R}_{{\rm{s}}}\right)}^{4}-{T}_{{\rm{d}}}^{4}).\end{eqnarray} \tag{ 22 }$$

This equation represents the balance of the energy in the thin layer of the pebbles accreted during the unit time on the surface of the satellites. The left-hand side is the energy necessary for heating the thin layer. The terms of the right-hand side are the collisional energy of the accreted pebbles and the emission from the surface. The pebble–satellite collision velocity u_i can be estimated by

$$\begin{eqnarray}&&{u}_{{\rm{i}}}=\min (\sqrt{{v}_{\mathrm{esc}}^{2}+{\rm{\Delta }}{v}^{2}},{v}_{\mathrm{set}}),\end{eqnarray} \tag{ 23 }$$

where ${v}_{\mathrm{esc}}=\sqrt{2{{GM}}_{{\rm{s}}}/{R}_{{\rm{s}}}}$ and v_set = g_st_stop are the escape and settling velocities, respectively. The gravitational acceleration at the surface of the satellite is g_s = GM_s/R_s².

Heating by ²⁶Al is also very effective. The increase of the satellite internal temperature at R, the distance from the center of the satellite, during the formation is

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{T}_{\mathrm{fin}}(R) & = & \displaystyle \frac{1}{{C}_{{\rm{p}}}}{\displaystyle \int }_{{t}_{\mathrm{acc}}}^{{t}_{\mathrm{fin}}}{m}_{{\rm{r}}}{q}_{26}(t){dt}\\ & = & \ \displaystyle \frac{{m}_{{\rm{r}}}{q}_{\mathrm{26,0}}}{{C}_{{\rm{p}}}{\lambda }_{26}}\exp (-{\lambda }_{26}{t}_{\mathrm{acc}}),\end{array}\end{eqnarray} \tag{ 24 }$$

where t_acc and t_fin are the time when the satellite radius R_s was equal to R and the end of the formation, respectively (Barr & Canup 2008). We assume the initial heating rate by ²⁶Al, the specific heat, and the satellite density as q_26,0 = 1.82 × 10⁻⁷ W kg⁻¹, C_p = 1400 J kg⁻¹ K⁻¹, and ρ_s = 1.9 g cm⁻³, respectively, which are the values in the table (Barr & Canup 2008). Here we assume that all ²⁶Al heat has remained at the point until the end. Thermal diffusion can be ignored because it diffuses only ∼10 km in 10⁷ yr. Solid-state convection can also be ignored because it can start at t ∼ 10⁸ yr (Barr & Canup 2008). Latent heat is not included because the purpose of this estimation is to determine whether the satellites melt or not. The final (i.e., maximum) satellite internal temperature at R can then be estimated by the sum of the two heating sources, the accretion heating and the ²⁶Al heating,

$$\begin{eqnarray}&&{T}_{\mathrm{fin}}(R)={T}_{{\rm{s}}}(R)+{\rm{\Delta }}{T}_{\mathrm{fin}}(R).\end{eqnarray} \tag{ 25 }$$

Note that the ²⁶Al heating is dominant (see Section 4.8).

There are two main mechanisms that make satellites migrate in the disk, aerodynamic drag and Type I migration. We consider both mechanisms at the same time to calculate the migration of satellites. The aerodynamic drag migration velocity can be obtained by substituting the Stokes number of satellites, St_s, for St_p in Equation (11),

$$\begin{eqnarray}&&{v}_{\mathrm{ad}}=-2\displaystyle \frac{{\mathrm{St}}_{{\rm{s}}}}{{{\mathrm{St}}_{{\rm{s}}}}^{2}+1}\eta {v}_{{\rm{K}}}.\end{eqnarray} \tag{ 26 }$$

The Stokes number of satellites is

$$\begin{eqnarray}&&{\mathrm{St}}_{{\rm{s}}}=\displaystyle \frac{8}{3{C}_{{\rm{D}}}}\displaystyle \frac{{\rho }_{{\rm{s}}}{R}_{{\rm{s}}}}{{\rho }_{{\rm{g}}}\eta {v}_{{\rm{K}}}}{{\rm{\Omega }}}_{{\rm{K}}}.\end{eqnarray} \tag{ 27 }$$

The Type I migration velocity is

$$\begin{eqnarray}&&{v}_{{\rm{t}}1}={b}_{{\rm{t}}1}\left(\displaystyle \frac{{M}_{{\rm{s}}}}{{M}_{\mathrm{cp}}}\right)\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}{r}^{2}}{{M}_{\mathrm{cp}}}\right){\left(\displaystyle \frac{{v}_{{\rm{K}}}}{{c}_{{\rm{s}}}}\right)}^{2}{v}_{{\rm{K}}},\end{eqnarray} \tag{ 28 }$$

where b_t1 is the migration constant, which depends on the distribution of the temperature and gas surface density of the CJD (see Equation (10) of Ogihara et al. 2015). If b_t1 is negative or positive, the satellite migrates inward or outward, respectively. The migration timescale is then ${t}_{\mathrm{mig}}=r/| {v}_{\mathrm{ad}}+{v}_{\mathrm{mig}}|$.

According to recent N-body simulation, the critical migration timescale for capture into MMRs is (Ogihara & Kobayashi 2013)

$$\begin{eqnarray}&&{t}_{\mathrm{crit}}={C}_{\mathrm{MMR}}{\left(\displaystyle \frac{{M}_{\mathrm{in}}}{{M}_{\oplus }}\right)}^{-4/3}{\left(\displaystyle \frac{{M}_{\mathrm{cp}}}{{M}_{\odot }}\right)}^{4/3}{T}_{\mathrm{in}},\end{eqnarray} \tag{ 29 }$$

where M_in, M_⊕, and M_⊙ are the mass of the inner satellite, Earth, and Sun, respectively, and T_in is the orbital period of the inner satellite. If the migration timescale t_mig is longer than this critical timescale, the two bodies can be captured into the resonance. The capture coefficient C_MMR depends on the type of the resonance and the mass ratio of the two bodies. We summarize the values of C_MMR in Table 1 cited from Ogihara & Kobayashi (2013), which calculated numerical simulations (Ogihara & Kobayashi 2013). Note that our 1D model does not consider the eccentricity and inclination of the orbits.

We calculate the evolution of the mass and orbits of four planetesimals (Seeds A1 to A4) captured by the CJD one by one. We assume that the capture times of Seeds A1 to A4 are t_cap = 1.0, 1.25, 1.5, and 2.0 Myr, respectively. The disk condition is x = 0.0026, α = 10⁻⁴, and r_gg = 1.7 × 10⁻⁷. We also assume the initial mass and positions of the seeds as M_s,start = 3 × 10²³ g and r_s,start = 50 R_J, respectively. Figure 4 represents the evolution of the size and orbits of the seeds, and Table 2 lists the final mass and positions.

Zoom In Zoom Out Reset image size

Figure 4 and Table 2 show that the mass and sizes of the Galilean satellites can be reproduced very well. The dichotomy of the size between the inner and outer two satellites is created by the assumption that the pebble mass flux inside the snow line is half of that outside because icy pebbles evaporate inside the snow line.

The top panel of Figure 5 represents the pebble accretion efficiencies, P_eff, of Seeds A1 to A4. All pebble accretion efficiencies are smaller than ∼10%. The changes in the slope around M_s ≈ 5 × 10²³ g are caused by the stop of their migration due to the resonance traps. The mild changes in the gradients around M_s ≈ 2 × 10²⁴ g and M_s ≈ 10²⁵ g correspond to the transition of the pebble accretion regimes from the 3D to the 2D and from the headwind to the shear regimes, respectively. The Stokes number inside of the snow line is smaller than that outside by a factor of 2–3 (Figure 3), which makes the pebble accretion efficiencies of Io and Europa larger than those of Ganymede and Callisto in the 2D regime because P_eff is a decreasing function (see Equations (16) and (18)). The bottom panel represents the growth timescale of the seeds, t_grow. The timescale increases with mass and finally shoots up as the CJD disperses. It also shows that Ganymede grows fastest because the pebble mass flux becomes half after passing the snow line, which has more impact than the reduction in St in the 3D and the 2D shear regimes (see Equations (1), (16)–(18)). The small jump around M_s ≈ 5 × 10²⁵ g on Europa's curve is caused by the crossing of the snow line. The filtering effects are mildly important for the growth timescale. Note that their mass does not reach the PIM, M_iso (Equation (33)).

Zoom In Zoom Out Reset image size

Figure 4 also shows that all seeds migrate quickly (<3 × 10⁵ yr) by aerodynamic drag (not by Type I migration) and are captured into 2:1 resonances one by one from the inner ones. After the seeds are captured into the resonances, they grow by pebble accretion without migration and keep their orbits on the current ones. The position of Seed A4, on the other hand, differs from the real orbit of Callisto. Seed A4 is also captured into a 2:1 resonance with Seed A3. However, Callisto's orbit may be able to expand rapidly, escaping from the resonance by the resonant dynamical tide works (Fuller et al. 2016).

Each Galilean satellite has a different ice mass fraction, and in particular, the low ice mass fraction of Europa is very unique (Table 2). We find that this low ice mass fraction of Europa (6%–9%) can be reproduced by the migration of the snow line at the final phase of the formation (Table 2). Figure 4 shows that Seed 2 (Europa) accretes icy pebbles after 10 Myr. Although the ice mass fraction strongly depends on the disk temperature profile, there is a disk condition that is suitable for reproducing the ice mass fractions of all the Galilean satellites. In our model ice sublimation occurs instantaneously at the snow line, in contrast to Ronnet et al. (2017), where pebbles only gradually lose their ice. Europa also naturally acquired an icy surface on top of a rocky interior, because the satellite accretes dry pebbles before accreting ice-rich pebbles.

We find that in order to avoid differentiation of Callisto by ²⁶Al heat, its seed must be captured by the disk late enough. The solid curves of Figure 7 represent the internal temperature of Ganymede (Seed A3; light blue) and Callisto (Seed A4; orange). The first one is higher than the melting point of Callisto (black), and the second one is lower than it. This means that Callisto does not melt but Ganymede may melt by ²⁶Al heat. The dichotomy of their internal ice-rock differentiation can be created by the difference in their capture time, 0.5 Myr, because the half-life of ²⁶Al is 0.717 Myr. The long growth timescale (∼10⁷ yr) is the reason why such a different capture time is allowed. If the growth timescale is shorter and the difference in the capture time is the same, the final mass of Ganymede and Callisto would end up too large. Indeed, it is difficult to make the dichotomy between the internal structures of the two satellites by the classical Canup-Ward formation scenario, where the growth timescale is <10⁶ yr (Barr & Canup 2008). Note that once rock-metal cores form in Europa and Ganymede, they can also differentiate metallic cores by long-lived radiogenic heating (Spohn & Breuer 1998), and Io can be differentiated completely by tidal heating after its formation with its current orbit (Peale et al. 1979). See also Section 5.1 discussing the subsequent evolution of Ganymede's internal structure.

We then calculate the evolution of the satellites in a situation in which three planetesimals (Seeds B1 to B3) are captured at first and the fourth seed (B4) forms and grows at the gas pressure bump of Seed B3. We assume that Seed B4 is put at its orbit and it starts to grow by pebble accretion when the mass of Seed B3 reaches its PIM, assumed to be equal to Ganymede's mass. We fix the orbit of Seed 4 at the gas pressure maximum formed by Seed B3, r = 17.0 R_J (see Appendix A). The other initial conditions are identical to those in Model A. Figure 6 represents the evolution of the size and orbits of the seeds, and Table 2 lists the final mass and positions.

Zoom In Zoom Out Reset image size

Figure 6 and Table 2 show that the mass of all the satellites is reproduced by Model B as well. We find that the evolution of the inner three satellites is identical to Model A, except for the filtering effect of Callisto. In Model A, Seed A4 starts to decrease the pebble mass flux inside its orbit at t−t_gap = 2 Myr. Conversely, Seed B4 appears only at t−t_gap = 9.54 Myr when Seed B3 reaches its PIM. Therefore, the inner satellites in Model B become slightly more massive than those in Model A. The growth of Callisto stops by the disappearance of the disk. Although we do not include the PIM of Callisto, the satellite should not reach pebble isolation under our assumption that the mass of the largest satellite, Ganymede, is defined as its PIM and PIM barely depends on the distance from Jupiter.

Figure 6 shows that the inner three satellites are captured into 2:1 MMRs also in Model B. On the other hand, in Model B we assume that the orbit of Callisto is at the gas pressure maximum made by Ganymede, r = 17.0 R_J. The distance between the two orbits, 2.2 R_J, is consistent with 5.0 r_H,G, where r_H,G = (M_G/(3M_J))^1/3r_G is the Hill radius of Ganymede and r_G = 14.8 R_J is the orbital radius of Ganymede (in Model B). This distance is slightly wider than $2\sqrt{3}{r}_{{\rm{H}},{\rm{G}}}$ but much narrower than 10 r_H,G. Therefore, Callisto is so close to Ganymede that the satellite should be scattered after the gas disk disappears and may reproduce its orbit, which is not captured in any resonance (Chambers et al. 1996).

We find that it is also possible in Model B to reproduce the small amount of ice on Europa (Table 2). Figure 6 shows that Seed B2 accretes icy pebbles only at the end of the growth because Seed B3 reaches its PIM and traps drifting pebbles at its orbit just after Seed B2 starts to accrete icy pebbles. As a result, Europa's ice mass fraction is only 2.9%, smaller than that in Model A but not much smaller than the actual value, 6%–9% (Kuskov & Kronrod 2005).

The dashed curves in Figure 7 represent the internal temperature of Seeds B3 (light blue) and B4 (orange). The internal temperature of Seed B3 is similar to Seed A3 and higher than the melting point of Callisto (black). The internal temperature of Seed B4 is lower than the melting point in the whole interior region and is almost uniform. This is because ²⁶Al radiogenic decay does not heat the satellite since Seed B4 starts to grow at a later time (9.54 Myr), which is naturally achieved in Model B. Moreover, the difference between the internal temperatures of Ganymede and Callisto is larger than that of Model A. We also find that the temperature is dominated by accretion heating and depends on the dust-to-gas accretion rate ratio.

Zoom In Zoom Out Reset image size

In our scenario, Jupiter has to grow so large that the gap structure opens and the pebbles can be trapped at the gas pressure maximum around Jupiter at 1 Myr after the formation of the calcium–aluminum-rich inclusions (CAIs) to make the dichotomy of the internal structures of Ganymede and Callisto by ²⁶Al heat. This is consistent with a recent "early Jupiter formation" scenario showing that solid materials in the solar system were spatially separated by that time (Kruijer et al. 2017).

We have assumed that the mass of young Jupiter is 0.4 M_J when the gap forms and it becomes 1 M_J at the end of the calculation. In our models, the mass flux of the gas inflow to the CJD, ${\dot{M}}_{{\rm{g}}}$, decreases as an exponential decay with the timescale of t_dep = 3 Myr, which is consistent with the observations of the lifetime of the CSDs (Haisch et al. 2001). It is generally believed that even after a gap forms, the gas accretion onto the planet (i.e., the gas inflow to the CPD) continues and its flux is determined by the gas surface density of the CPD inside the gap (Tanigawa & Tanaka 2016). Therefore, our assumption about the gas inflow is plausible. We note, however, that if the viscous gas accretion rate onto the Sun is lower than the photoevaporation rate, the gas surface density of the CSD and the gas accretion rate onto the CJD quickly decrease (Alexander et al. 2006a, 2006b). In this case, the migration of the snow line must be quick and the supply mechanism of ice to Europa may not work well, but the accretion of partially dehydrated pebbles argued in Ronnet et al. (2017) may be an alternative way.

In our model, we have fixed the position of the disk inner edge as the current orbit of Io, 5.89 R_J, for simplicity, but this assumption is plausible. It is generally accepted that Jupiter's magnetic field was stronger than the current one, resulting in magnetospheric accretion and opening of the cavity in the latter phase of Jupiter's formation, i.e., after the gap formed around the CJD (Stevenson et al. 1983; Sánchez-Lavega 2004; Christensen et al. 2009). The position of the edge of the cavity can be estimated from the balance of the gas accretion rate and magnetic stress by Jupiter's magnetic field (Equation (3)). Since the gas accretion keeps decreasing after the gap opens, the edge moves outward and the innermost seed (Io) may also migrate together (Liu et al. 2017). Then, the edge should stop at the corotation radius. It is generally believed that Jupiter was larger than today (Burrows et al. 1997; Fortney et al. 2011). Conceivably, its spin frequency was lower and the corotation radius was located farther than the current, r_co = 2.25 R_J.

First, we check whether Io can stop its inward migration at the edge and move outward together or not. The two-sided (i.e., normal Type I migration) torque that a satellites receives is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{2{\rm{s}}}=\displaystyle \frac{1}{2}{M}_{{\rm{s}}}{v}_{{\rm{K}}}{v}_{{\rm{t}}1}.\end{eqnarray} \tag{ 30 }$$

On the other hand, a satellite at the disk edge receives a strong one-sided positive corotation torque that pushes the satellite outward (Liu et al. 2017; Romanova et al. 2019). The one-sided corotation and Lindblad torque are

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{1{\rm{s}},\mathrm{co}}=2.46\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}{r}^{2}}{{M}_{\mathrm{cp}}}\right){\left(\displaystyle \frac{{\mu }_{{\rm{s}}}}{{h}_{{\rm{g}}}^{3}}\right)}^{1/2}{M}_{{\rm{s}}}{\left(r{{\rm{\Omega }}}_{{\rm{K}}}\right)}^{2}\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{1{\rm{s}},\mathrm{Lin}}=-0.65\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}{r}^{2}}{{M}_{\mathrm{cp}}}\right)\displaystyle \frac{{\mu }_{{\rm{s}}}}{{h}_{{\rm{g}}}^{3}}{M}_{{\rm{s}}}{\left(r{{\rm{\Omega }}}_{{\rm{K}}}\right)}^{2},\end{eqnarray} \tag{ 32 }$$

respectively (Liu et al. 2017). Figure 8 shows that the one-sided corotation torque that Io (Seed A1) receives is much larger than the negative one-sided Lindblad torque and the negative two-sided torque that the other satellites (Seeds B2, B3, and B4) receive, respectively. Since the sum of the negative torque is smaller than the positive torque, even if the other satellites push Io through the chain of the resonance, Io keeps its position on the edge of the cavity. It is also argued that if the migration timescale of the edge (i.e., 3.5t_dep) is shorter than that of the migration of the satellite by the one-sided corotation torque, the disk edge leaves the satellite there (Liu et al. 2017). We found, however, that the timescale of Io's migration by the one-sided corotation torque is ∼10³ yr, and this is much shorter than the migration timescale of the edge, ∼10⁶ yr. Therefore, Io may have moved much farther than its current orbit when the disk inner edge moves outward.

Zoom In Zoom Out Reset image size

However, this outward migration of the inner edge should stop at the corotation radius r_co, where the Keplerian frequency of the disk equals the spin frequency of Jupiter. When r_cav > r_co, there will be two possibilities: the angular momentum will be transferred from Jupiter to the disk and then the gas accretion will stop, or otherwise the corotation radius and the disk edge will move outward together and then the accretion will continue (Takata & Stevenson 1996; Liu et al. 2017). Although the current corotation radius is r_co ≈ 2.25 R_J, Jupiter at its time of formation was much larger than it is today (Burrows et al. 1997; Fortney et al. 2011), and this means that the corotation radius was also larger than the current one if the conservation of the angular momentum of Jupiter is assumed. Considering the transport of the angular momentum from Jupiter to the disk, the angular momentum should have been conserved since the disk disappeared. According to a formation model of Jupiter, the radius of the planet was ≈1.75 R_J after its rapid gas accretion and it decreased little by little (Lissauer et al. 2009). When the radius of Jupiter is 1.75 R_J, the corotation radius should be r_co ≈ 4.7 R_J. We can then consider two scenarios of Io formation. In the first one, Io formed around r ≈ 4.7 R_J, slightly interior to the r = 5.89 R_J of our fiducial model, and then moved outward after the disk dissipated. The satellites, especially the inner ones, could move outward by the tidal force from Jupiter (Yoder & Peale 1981). The outer ones would be pushed by the inner ones and move outward with them because of the resonance. In this case, the position of the snow line should have been more inside than the fiducial case in this work, but this thermal condition could easily be reproduced by another parameter set. The second possibility is that Io was not the innermost satellite. If a body was present at r = 3.7 R_J, Io would have been situated at r = 5.7 R_J if they were trapped in a 2:1 resonance. This orbit is consistent with the corotation radius when the radius of Jupiter is ≈1.5 R_J, and this radius can be achieved during the contraction of Jupiter. The innermost body may have been broken by the tidal force of Jupiter when it has entered inside the Roche limit. Current Io, trapped in the Laplace resonance, actually moves inward little by little because of the tidal dissipation, and the innermost body may have also experienced such inward migration (Lainey et al. 2009).

The ionization degree of the disk inner region, χ_e ≡ n_e/n_n, where n_e and n_n are the number densities of electron and neutral gas, respectively, is also important. Without enough ionization, the disk gas cannot couple with the rotating magnetic field of Jupiter, and the inner cavity does not open. The angular momentum transfer can occur if the magnetic Reynolds number R_m ≡ VH_g/λ > 1, where V is the relative velocity between the rotating magnetic field and the disk gas and λ is the magnetic diffusivity (Takata & Stevenson 1996). We note that this magnetic Reynolds number is ∼10³ times larger than the Elsässer number, the critical number for the MRI activation, Λ = v_Az²/(λΩ_K), where v_Az is the z component of the Alfvén velocity. If Jupiter and its magnetic field spin rigidly at the current speed, V ≈ 3 × 10⁴ m s⁻¹, H_g ≈ 1.1 × 10⁷ m, and λ ≈ 0.74/χ_e at r = 1 R_J with disk temperature T = 1000 K (Blaes & Balbus 1994). The condition for angular momentum transport is then χ_e ≳ 10⁻¹²–10⁻¹¹. Takata & Stevenson (1996) estimated the ionization degree of the CJD including the effects of galactic cosmic rays, radioactive isotope decay, and so on. They assumed the minimum mass disk model by Lunine & Stevenson (1982), whose gas surface density is 100 times larger than that of our model, and argued that the ionization degree is χ_e ≲ 10⁻¹⁶ at the midplane, χ_e ∼ 10⁻¹⁵ at the altitude of one scale height, and χ_e ∼ 10⁻¹³–10¹² at three scale heights. If the ionization degree is inversely proportional to the gas surface density, although it is a very rough estimate, the ionization degree in our disk model should be χ_e ≲ 10⁻¹⁴ at the midplane, χ_e ∼ 10⁻¹³ at the altitude of one scale height, and χ_e ∼ 10⁻¹¹–10⁻¹⁰ at three scale heights. Therefore, the condition for angular momentum transport can be achieved at an altitude higher than the scale height. A possible mechanism to curve the inner cavity of the gas disk is that the transport of angular momentum (i.e., the radial gas flow) is dominated in the upper region and the quick vertical relaxation to hydrostatic conditions provides an accompanying vertical upward drift of gas (see Takata & Stevenson 1996). However, this discussion is based on a very rough estimate; more detailed investigation (e.g., MHD simulations) of the disk inner region should be carried out in the future.

Moreover, a recent work by Arakawa & Shibaike (2019) found that photophoresis in the CJD stops the inward drift of dust particles near the orbit of Io and carves an inner cavity in a broad range of initial conditions. This result also supports our assumption of the inner cavity strongly.

We have also assumed that the strength of turbulence of the viscous accretion disk is α = 10⁻⁴, consistent with the inability of the MRI to operate in the CJD (Fujii et al. 2014). From the assumptions of the strength of turbulence and the gas accretion rate, we calculate the gas surface density of the disk. The surface density affects many properties of the system, for example, the disk temperature, the pebble accretion rates, and the migration speeds.

To get the disk thermal condition suitable for reproducing the ice mass fractions of the satellites, we needed to know that the Rosseland mean optical depth τ is around unity at Ganymede's orbit. Figure 9 represents the r_gg dependences of the snow line, satellite mass, and ice fraction of satellites in the case of the fixed orbits. The top panel shows that when r_gg = 1.7 × 10⁻⁷, the snow line migrates inward from just inside Ganymede's orbit and stops just inside Europa's. This is because the first term is dominant in the opacity factor g in the first half of its evolution and then the second term becomes dominant as the gas surface density decreases (see Equation (6)). If r_gg is larger than this, the first term is dominant until almost the end, and vice versa. The middle panel shows that the mass distribution of the Galilean satellites can be reproduced only when r_gg = 1.7 × 10⁻⁷. However, if the gas accretion rate ${\dot{M}}_{{\rm{g}}}$ is smaller than this assumption and the dust-to-gas accretion rate ratio x is correspondingly larger, the snow line will again be near Ganymede and the mass distribution of the satellites will be reproduced. The bottom panel shows that the r_gg dependence of the ice mass fraction is large. When r_gg is larger than the most suitable case, the slope inside the first position of the snow line is gentler and the ice mass fraction of Europa is larger. On the other hand, when r_gg is smaller, the slope is steeper and Europa cannot get enough ice. These results can be understood by the top panel.

Zoom In Zoom Out Reset image size

However, this constraint of marginal optical thickness could be plausible. First, the opacity of the drifting pebbles is consistent with our assumption. In the geometric limit (pebbles larger than the wavelength), the opacity of pebbles is given by κ_p ∼ 1/(ρ_int,pa_p), where ρ_int,p and a_p are the pebble internal density and the pebble radius, respectively. If the pebbles are fluffy, the opacity of pebbles could be κ_p ∼ 1/(10⁻³ × 10) = 10² cm² g⁻¹ (Kataoka et al. 2014; Shibaike et al. 2017). In this case, the mean optical depth of pebbles is τ_p = κ_pΣ_p ∼ 10² × 10⁻² = 1 (see Figure 3). This is consistent with the Rosseland mean optical depth in our models, τ = 450r_ggΣ_g ∼ 500 × (2 × 10⁻⁷) × 10⁴ ∼ 1 (see Equation (7)). Note that the gas (molecular) opacity (∼10⁻⁵–10⁻⁴ cm² g⁻¹, consistent with an optical depth of ∼0.1–1 with Σ_g ∼ 10⁴ g¹ cm⁻²) may also affect the disk temperature, but the effect should be limited (Mizuno 1980). Second, in reality Equations (5) and (6) may not be applicable because the dissipation of accretion energy occurs in the disk upper regions if the disk is laminar and wind-driven accreting. This implies that optically thick disks can still have cold midplanes (Hirose et al. 2009; Mori et al. 2019). Under these conditions, our model will also work with higher r_gg.

Note that we simply assumed the snow line as the position where the disk temperature is 160 K after previous models of PPDs, that the gas pressure of the CJD is much higher than that of PPDs, and that the partial water vapor pressure and the condensation temperature should also be higher. In this case, the evolution track of the snow line should be shifted closer to Jupiter, causing the increase of the ice mass fraction of Europa. Our model will be able to reproduce the ice-depleted Europa again if we choose a larger ${\dot{M}}_{{\rm{g}}}$, but it means a shorter formation timescale and may make the difference between the internal structures of Ganymede and Callisto smaller.

One of the strong points compared with previous scenarios is that our slow pebble accretion scenario only needs a small amount of solid material (i.e., dust particles and planetesimals). The total amount of dust needed to drive the growth in the slow pebble accretion scenario is also modest and is smaller than that for the classical satellitesimal accretion scenarios (Canup & Ward 2006; Shibaike et al. 2017). Our scenario can reproduce the mass of the current Galilean satellites even with the dust-to-gas accretion ratio is as low as x = 0.0026, which is smaller than the solar composition by a factor of four. This is consistent with the fact that the gas flowing into the CJD is depleted in solids (Canup & Ward 2002; Tanigawa et al. 2012). Because of the gap formation in the CJD caused by Jupiter, only small particles (≲0.1 mm) can overcome the ensuing gas pressure gradient to end up in the CJD (Zhu et al. 2012), In addition, the accreted gas should be supplied from high altitude, where gravitational settling of larger particles limits the amounts of dust (Tanigawa et al. 2012). On the other hand, small amounts of particles (∼0.1 mm) can form as fragments of the collisions of the solids trapped at the gas pressure maximum (Kobayashi et al. 2012).

Furthermore, only four planetesimals need to be captured by the CJD. This assumption is consistent with the following facts: after the gap opens in the CSD, such captures only appear for the planetesimals that have high eccentricities in the CSD (Fujita et al. 2013), planetesimals can form at the gas pressure maximum of the gap, and those that are scattered by a large body will have high eccentricities (Kobayashi et al. 2012; Ronnet et al. 2018).

We consider the cases in which the planetesimals that have formed in the CJD are captured by the CJD and become the seeds of the satellites. The assumptions of a starting location of r_s,start = 50 R_J with zero eccentricity are plausible. According to a previous work, there are orbits in which planetesimals are captured by multiple approaches to Jupiter and become circularized around ≈50 R_J (long-lived prograde captured orbits; Suetsugu et al. 2016).

We also fix the mass of the captured planetesimals as 3 × 10²³ g (about 300 km), which is consistent for the properties of the captured orbits (Suetsugu et al. 2016). These planetesimals are large enough to start growing by efficient pebble accretion (i.e., in the settling regime) for the adopted disk properties (see Equation (20)). Smaller planetesimals, on the other hand, will stay small and may be scattered or accreted by large ones. In the settling regime, any difference in the initial mass does not change the final mass because the growth timescale increases with mass (see the bottom panel of Figure 5). In addition, planetesimals smaller than the assumed value are likely to end up in higher-order resonances (i.e., 3:2 instead 2:1), by virtue of their faster migration (see Equations (26) and (29) and Table 1). Therefore, the initial seed mass must be similar to or larger than our standard value.

Each capture time of the planetesimals t_cap is assumed to make the dichotomy of the internal structures of the Galilean satellites; Callisto is minimally or only modestly differentiated, but the other satellites are fully differentiated. It is difficult to make the dichotomy between the internal structures of Ganymede and Callisto by the difference in the heat of accretion or ²⁶Al decay in the previous Canup-Ward formation scenario (Canup & Ward 2006; Barr & Canup 2008). In our slow pebble accretion scenario, we found that only Ganymede was differentiated by ²⁶Al heat, provided the seed of Callisto was captured by the disk at a sufficiently late time.

Figure 10 represents the internal temperature of Callisto where its seed is captured at different times. We fix the position of Callisto at the current orbit in this calculation. Without the ²⁶Al heat, the surface of Callisto would remain at the same temperature as the disk. Heating by pebble accretion can be neglected, despite the fact that the ice-rich (big) pebbles impact the surface at the escape velocity, i.e., similar to planetesimal heating (see Section 2.5). However, whereas in the classical planetesimal-driven accretion scenario accretion proceeds quickly, accretion timescales in our pebble-driven scenario are quite long, 10 Myr (except for Seed B4). This timescale is much longer than the melting critical accretion timescale, 0.6 Myr (Barr & Canup 2008).

Zoom In Zoom Out Reset image size

Only in the cases where ²⁶Al radiogenic heating is considered does the internal temperature exceed the melting point of water ice. Figure 10 shows that this condition is achieved if most of the ²⁶Al in the material of Callisto has decayed before its formation starts, in other words, t_gap(=t_cap) > 2 Myr. Also, if Ganymede starts to form by 1.5 Myr after the formation of CAIs and Callisto starts later than 2.0 Myr, the dichotomy between the two satellites should be explained. Actually, in Model A we assumed that the capture times of the seeds of Ganymede and Callisto were 1.5 and 2.0 Myr after the CAI formation, respectively. In Model B, Callisto's seed formed at the gas pressure bump made by Ganymede at 10.5 Myr, when Ganymede reaches its PIM.

We also note that if there is less ²⁶Al in pebbles (Larsen et al. 2016), the internal temperature of satellites must be lower. The curve of the internal temperature in R_s–T_s space should approach the black curve (w/o ²⁶Al) in Figure 10. In that case, not only Callisto but also Ganymede should be undifferentiated during their formation.

If the collision speed of pebbles is too fast, they cannot merge and grow larger but fragment. The collision speeds depend on the Stokes number of pebbles, so there is a limit to the size of pebbles. We have assumed that the critical fragmentation speeds of rocky and icy pebbles are v_cr = 5 and 50 m s⁻¹, respectively. As a result, the size (i.e., Stokes number) of pebbles is determined by fragmentation inside the snow line, although drift determines it outside the snow line (see Section 2.3). In this work, we have adopted the critical fragmentation speeds from numerical studies of collisions of pebbles with 0.1 μm sized monomers (Wada et al. 2009, 2013). This monomer size is consistent with the typical size of the grains constituting interplanetary dust particles of presumably cometary origin (Rietmeijer 1993). Moreover, the results of previous experimental work are consistent with those of the numerical simulations (Blum & Wurm 2000; Poppe et al. 2000; Gundlach & Blum 2014). Therefore, our assumption of the critical fragmentation speed is plausible.

We note that there is still uncertainty about the fragmentation threshold. A recent experimental work argued that the critical speeds are slower than our assumptions, ∼1 m s⁻¹ for both the icy and rocky pebbles (Musiolik & Wurm 2019). On the other hand, Kimura et al. (2015) and Steinpilz et al. (2019) claimed that the critical fragmentation speed is faster than our fiducial value, ∼10 m s⁻¹ for silicate particles. The top panel of Figure 11 shows that the Stokes number of pebbles with the slower critical speeds (purple) is smaller than the fiducial case (light blue). The size of pebbles is determined by fragmentation not only inside the snow line but also outside. On the other hand, if v_cr = 10 m s⁻¹ for rocky pebbles (yellow), the Stokes number is larger than the fiducial case and is determined by the drift except for the region inside several Jupiter radii. However, the bottom panel shows that although the final mass of the satellites with the slow critical speeds (solid purple) is smaller than the fiducial case (light blue), it can be compensated by changing the other parameters. The distribution of the final mass with v_cr = 1 m s⁻¹ (dashed purple) is similar to the fiducial case when α = 10⁻⁵, r_gg = 1.7 × 10⁻⁸, and x = 0.0011. The final mass with the higher critical speeds (yellow) is also similar to the fiducial case. Therefore, the effect of the difference of the critical fragmentation speed on our formation scenario is limited.

Zoom In Zoom Out Reset image size

In Model B, we define that the PIM of Ganymede is its actual current mass, 1.48 × 10²⁶ g. Here we show that this assumption is plausible. An expression of the PIM estimated by 2D simulations (Ataiee et al. 2018) is

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{iso}} & = & {h}_{{\rm{g}}}^{3}\sqrt{37.3\alpha +0.01}\\ & & \times \ \left\{1+0.2{\left(\displaystyle \frac{\sqrt{\alpha }}{{h}_{{\rm{g}}}}\sqrt{\displaystyle \frac{1}{{{\mathrm{St}}_{{\rm{p}}}}^{2}}+4}\right)}^{0.7}\right\}{M}_{\mathrm{cp}}.\end{array}\end{eqnarray} \tag{ 33 }$$

We find that this mass is 7.65 × 10²⁵ g with the parameters of Model B, which is about 1.9 times smaller than the mass of Ganymede. The PIM derived from 2D simulations can be about 1.5–2 times smaller than the expression by 3D simulations, and the difference can become larger with weak turbulence at least in PPDs (Ataiee et al. 2018; Bitsch et al. 2018; Johansen et al. 2019). Therefore, the PIM of Ganymede must be around the actual mass of the satellite. We also note that Ganymede is the largest satellite, so only that satellite should reach its PIM.

We conclude that Ganymede's interior differentiates into the rocky core and the icy mantle because of the ²⁶Al heat in the accreted pebbles. However, strictly speaking, it can melt only the region 1000 km from the center of Ganymede (see Figure 7). Therefore, we should consider the subsequent evolution of the interior, an overturn of the pristine (undifferentiated) icy-rocky crust and the pure icy mantle by the Rayleigh–Taylor (R-T) instability. Here we consider a three-layered interior structure, (1) the pristine (undifferentiated) icy-rocky crust (1000–2600 km from the center), (2) the pure icy mantle (600–1000 km), and (3) the rocky core (0–600 km), as the initial condition. We estimate the depth of the layers by assuming that the densities of the pure ice, rocky core, and pristine crust are ρ_ice = 1.4 g cm⁻³, ρ_rock = 3.0 g cm⁻³, and ρ_s = 1.9 g cm⁻³, respectively.

First, we check whether the R-T instability occurs or not. The condition depends on the viscosity of the upper layer, which is the pristine crust in this situation. If the viscosity is lower than the following critical value, the instability occurs with the timescale t_RT (Rubin et al. 2014),

$$\begin{eqnarray}&&{\eta }_{\mathrm{crit}}=\left[{\left(n-1\right)}^{1/n}\displaystyle \frac{{C}_{{\rm{L}}{\rm{\Delta }}\rho }}{2n}\right]{\left(\displaystyle \frac{{Z}_{0}}{L}\right)}^{(n-1)/n}{\rm{\Delta }}\rho {g}_{{\rm{s}}}({R}_{{\rm{s}}}){{Lt}}_{\mathrm{RT}},\end{eqnarray} \tag{ 34 }$$

where n, C_LΔρ, and Z₀ are the index of the stress related to the strain rate, a dimensionless quantity determined by the geometry and rheology, and the initial perturbation amplitude, respectively. We assume n = 1.8, C_LΔρ = 0.76, and Z₀ = 1 km (see Rubin et al. 2014). The difference of density is Δρ = ρ_s − ρ_ice = 500 kg m⁻³, and g_s(R_s) = 0.54 m s⁻¹ is the gravitational field of the satellite at R_s(=1000 km), the distance from the center. The length scale over which the viscosity varies significantly is

$$\begin{eqnarray}&&L=\displaystyle \frac{{{nRT}}_{0}}{{E}_{{\rm{a}}}}\displaystyle \frac{{T}_{0}}{| {{dT}}_{\mathrm{up}}/{dz}| },\end{eqnarray} \tag{ 35 }$$

where R, T₀, and E_a are the gas constant, the temperature of the interface between the two layers, and the activation energy, respectively. We assume T₀ = 370 K and E_a = 49 kJ mol⁻¹ (see Figure 7 and Rubin et al. 2014). According to Figure 7, the temperature gradient across the upper crust is ${T}_{0}/| {{dT}}_{\mathrm{up}}/{dz}| \,=0.14\,{\rm{K}}\,{\mathrm{km}}^{-1}$, where z is the vertical distance from the interface. We can then calculate the critical viscosity, η_crit = 7.6 × 10¹³(t_RT/year) [Pa s]. The viscosity of the upper layer should strongly depend on the temperature. According to the Arrhenius functions, the viscosity can be estimated by

$$\begin{eqnarray}&&{\eta }_{\mathrm{up}}={\eta }_{\mathrm{ref}}\exp \left[A\left(\displaystyle \frac{{T}_{\mathrm{ref}}}{{T}_{\mathrm{up}}}-1\right)\right],\end{eqnarray} \tag{ 36 }$$

where η_ref, T_up, and T_ref = 273 K are the reference viscosity, the temperature of the upper crust, and the reference temperature (melting point of pure ice), respectively. A constant coefficient A is about 20–25, and the reference viscosities of pure ice I and pure ice V and VI are η_ref ∼ 10¹⁴ [Pa s] and ∼10¹⁶–10¹⁷ [Pa s], respectively (Rubin et al. 2014; Shoji & Kurita 2014). If the volume rate of rock is smaller than 75%, the viscosity of rock-mixed ice is almost the same as that of pure ice (Durham et al. 2010). The minimum temperature of the upper crust is about 150 K (see Figure 7). If we substitute this temperature for Equation (36) and equate η_up with η_crit, the timescale of the R-T instability is 10⁹–10¹⁰ yr for A = 20 and 10¹¹–10¹² yr for A = 25. However, if we substitute T_up = 170 K, only 20 K larger, the timescale becomes 10–100 times shorter. Therefore, if the upper crust was heated up to 20 K by long-lived radiogenic heating or the released gravitational energy of halfway R-T instability itself (see Appendix B), R-T instability has been able to occur in the upper crust of Ganymede within the age of the satellite. One interesting observational fact of Ganymede is that there are two regions on the surface of the satellite: a very primitive region (∼4 Ga) and a relatively newer region (∼2 Ga). This characteristic is consistent with our estimate of the evolution of Ganymede's interior; the newer surface had been the inner pure icy mantle and was transported to the surface by the R-T instability around 2 Ga. On the other hand, Callisto's interior did not melt, and the surface has been kept as a primitive crust.

Finally, we highlight the possibility that the seeds of the satellites are not supplied from the CSD but form in situ around the snow line of the CJD by streaming instability (Schoonenberg & Ormel 2017). Ormel et al. (2017) and Schoonenberg et al. (2019) have used this idea toward understanding the properties of the TRAPPIST-1 system. In their model planets consecutively form at the snow line and then migrate inward. Such in situ formation of the seeds is in particular attractive to Model B, as it would no longer rely on capturing planetesimals. The properties of the satellites in this case will not change from those in the planetesimal-capture case because the seeds must migrate quickly and grow at each current orbital position of the satellites in both of the cases. The number of the Galilean satellites is also naturally explained in Model B. If the seeds forming at the snow line migrate inward and are trapped into 2:1 MMRs one by one, the orbits of the third seed must be around the snow line, and no more seeds can form there.

However, the main problem with this scenario is that in situ satellitesimal formation requires dust-to-gas ratios of ∼1 (Carrera et al. 2015; Yang et al. 2017), which is not reached by far because of the very fast pebble drift in the CJD (Shibaike et al. 2017). The increase in ice surface density by a factor of 3–5 at the snow line (Schoonenberg & Ormel 2017) would not be enough by far.