It is generally accepted that the four major Galilean satellites formed out of the gas disk that accompanied Jupiters formation. However understanding the specifics of the formation process is challenging as both small particles pebbles and the satellites are subject to fast migration processes. Here we hypothesize a new scenario for the origin of the Galilean system based on the capture of several planetesimal seeds and subsequent slow accretion of pebbles. To halt migration we invoke an inner disk truncation radius and other parameters are tuned for the model to match physical dynamical compositional and structural constraints. In our scenario it is natural that Ganymedes mass is determined by pebble isolation. Our slow pebble accretion scenario then reproduces the following characteristics 1 the mass of all the Galilean satellites 2 the orbits of Io Europa and Ganymede captured in mutual 21 mean motion resonances 3 the ice mass fractions of all the Galilean satellites and 4 the unique icerock partially differentiated Callisto and the complete differentiation of the other satellites. Our scenario is unique to simultaneously reproduce these disparate properties.
['accretion', 'accretion disks', 'planet–disk interactions', 'planets and satellites: dynamical evolution and stability', 'planets and satellites: formation', 'planets and satellites: individual: Jupiter', 'Galilean satellites', 'planets and satellites: interiors', 'Astrophysics - Earth and Planetary Astrophysics']
The Galilean Satellites Formed Slowly from Pebbles
{'THE GALILEAN SATELLITES FORMED SLOWLY FROM PEBBLES': 'Y uhito S hibaike 1,2 , C hris W. O rmel 3,4 , S higeru I da 1 , S atoshi O kuzumi 5 , and T akanori S asaki 6  \n- 1 Earth-Life Science Institute, Tokyo Institute of Technology, Meguro-ku, Tokyo 152-8550, Japan\n- 2 Physikalisches Institut & NCCR PlanetS, Universitaet Bern, CH-3012 Bern, Switzerland\n- 3 Anton Pannekoek Institute, University of Amsterdam, Science Park 904,1090GE Amsterdam, The Netherlands\n- 4 Department of Astronomy, Tsinghua University, Haidian DS 100084, Beijing, China\n- 5 Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Meguro-ku, Tokyo, 152-8551, Japan\n- 6 Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto, 606-8502, Japan', 'ABSTRACT': "It is generally accepted that the four major (Galilean) satellites formed out of the gas disk that accompanied Jupiter's formation. However, understanding the specifics of the formation process is challenging as both small particles (pebbles) as well as the satellites are subject to fast migration processes. Here, we hypothesize a new scenario for the origin of the Galilean system, based on the capture of several planetesimal seeds and subsequent slow accretion of pebbles. To halt migration, we invoke an inner disk truncation radius, and other parameters are tuned for the model to match physical, dynamical, compositional, and structural constraints. In our scenario it is natural that Ganymede's mass is determined by pebble isolation. Our slow-pebble-accretion scenario then reproduces the following characteristics: (1) the mass of all the Galilean satellites; (2) the orbits of Io, Europa, and Ganymede captured in mutual 2:1 mean motion resonances; (3) the ice mass fractions of all the Galilean satellites; (4) the unique ice-rock partially di GLYPH<11> erentiated Callisto and the complete di GLYPH<11> erentiation of the other satellites. Our scenario is unique to simultaneously reproduce these disparate properties.  \nKeywords: accretion, accretion disks - planets and satellites: dynamical evolution and stability - planets and satellites: formation - planets and satellites: individual (Jupiter, Galilean satellites) - planets and satellites: interiors - planet-disk interactions", '1. INTRODUCTION': "The four large satellites around Jupiter were discovered by Galileo Galilei over 400 years ago. The physical, dynamical, compositional, and structural properties of the satellite system are well known. The inner three satellites, Io, Europa, and Ganymede, are captured in mutual 2:1 mean motion resonances. The satellites' masses are similar at about 10 GLYPH<0> 4 M J (Jupiter mass). The satellites' ice mass fractions increase with their distance from Jupiter: Io is dry, Europa consist of 6 GLYPH<0> 9 wt% ice, while Ganymede and Callisto have ice mass fractions of about 50 wt% (Kuskov & Kronrod 2005). Uniquely, Callisto features an undi GLYPH<11> erentiated internal structure (Schubert et al. 2004).  \nAny model describing the formation of the Galilean satellites, must obey these basic observational constraints. However, previous models have only explained parts of these characteristics and the models are inconsistent with each other. Given that the satellites reside in the same plane, it is natural to assume that they formed in a gas disk surrounding Jupiterthe circum-Jovian disk (CJD)analogous to the formation of planets in circum-stellar disks (CSDs) (Canup & Ward 2002; Mosqueira & Estrada 2003). Previous works showed  \nthat if enough satellitesimals (km-sized bodies) exist in the disk, satellites with the current Galilean satellites mass can form by two-body collisions (Canup & Ward 2006; Sasaki et al. 2010). A problem with this scenario, however, is the formation of satellitesimals. Unless the CJD features pressure reversals or the gas flows outward on the midplane, the progenitor dust grains will not have had the time to conglomerate, because of the strong radial drift (Shibaike et al. 2017; Dra˙zkowska & Szul'agyi 2018). Alternatively, planetesimals expelled from the CSD can be captured by the CJD by virtue of the large density of the latter disk. A problem with this planetesimal-capture scenario is, however, that the growing Jupiter and the opened gas gap push the planetesimals out from Jupiter's feeding zone, rendering the capture rate very low (Hayashi et al. 1977; Fujita et al. 2013).  \nAnother problem formation models face is that of strong radial migration of planetesimals and satellites by aerodynamic and tidal forces. Previous works argued that today's Galilean satellites are the final survivors of a history in which an earlier generation of satellites repeatedly formed, but were lost because they migrated into Jupiter (Canup & Ward 2006; Sasaki et al. 2010; Ogihara & Ida 2012; Cilibrasi et al. 2018).  \nThis model is therefore very ine GLYPH<14> cient in its use of (already  \nscarce) solid material supplied to the CJD. An alternative idea, which we invoke in our model, is to stop the inward migration by virtue of a cavity of the gas disk around Jupiter opened by a strong magnetic field of the planet. The resonance chain of the inner three satellites' orbits is actually consistent with a scenario where the cavity halted the migration of Io and then Europa and Ganymede were captured into the resonances one by one (Sasaki et al. 2010; Ogihara & Ida 2012).  \nThe variation in the ice fraction also constitutes a formidable modeling challenge (Heller & Pudritz 2015a,b). In particular, the ice mass fraction of Europa (6 GLYPH<0> 9%)rather small compared to Ganymede and Callistois hard to explain by the previous satellitesimal-accretion scenario because the rocky and icy satellitesimals must be radially mixed beyond the snowline (Dwyer et al. 2013). Recently, a scenario where Europa accretes small icy particles, which dehydrated interior to the snowline, has been suggested to explain the small ice fraction (Ronnet et al. 2017). Another idea is to reproduce the fraction by accreting small particles by invoking the inward movement of the snowline (the place in the CJD corresponding to T = 160 K) in the final growth phases of the satellites (Canup & Ward 2009). A formation scenario for the TRAPPIST-1 system which has various ice mass fractions of planets may be able to apply to the Galilean satellites, where the seeds of the planets form around the snowline and accrete pebbles during their migration toward the star (Ormel et al. 2017).  \nFinally, the dichotomy between the di GLYPH<11> erentiated Ganymede and the partially di GLYPH<11> erentiated Callisto requires tuned conditions. Given that Ganymede and Callisto have similar mass and compositions, it is natural to likewise expect a similar thermal history for these neighboring satellites (Barr & Canup 2008). Thermal evolution after their formation by the release of the gravitational energy during the di GLYPH<11> erentiation and / or the impacts during the Late Heavy Bombardment might be able to make the dichotomy but still tuned conditions were needed (Friedson & Stevenson 1983; Barr & Canup 2010).  \nHere we show that we have succeeded to construct a new scenario, 'Slow-Pebble-Accretion scenario', which reproduces the physical, dynamical, compositional, and structural properties of the satellite system simultaneously and consistently. Naturally, in order to match the above constraints, our model is characterized by a number of parameters. Therefore, we do not argue that the Galilean satellite system is an inevitable result. It is a possible result, and with other choice of parameters, di GLYPH<11> erent satellite systems are reproduced. Nevertheless, a unique and strong point of our scenario is that the small amount of ice in Europa and the dichotomy between the internal structures of Ganymede and Callisto can be reproduced simultaneously by a single scenario with plausible assumptions.  \nThe structure of this paper is as follows. First, in Section  \n2, we explain the methods used in this work. We then show the results of our calculations in Section 3 and discuss the validity of our assumptions In Section 4. These two sections are the key parts of this paper. We also add some discussion in Section 5 and conclude our work in Section 6.", '2. METHODS': 'In order to create a scenario for the origin of the Galilean satellites consistent with almost all of the available constraints, we rely on combining a variety of modeling tools. We first provide an executive summary of the overall model (Section 2.1) and then discuss its elements in more detail.', '2.1. Model Summary': "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 the models, the captured planetesimals slowly accrete the particles drifting toward Jupiter, here referred to as pebbles ( GLYPH<24> 10 cm). The interiors of the satellites are mainly heated by the radiogenic decay of 26 Al included in the accreted pebbles. Figure 1 represents the two models of our slow-pebble-accretion scenario.  \nAt first, we model the gas accretion rate as ˙ M g = 0 : 2 exp( GLYPH<0> ( t GLYPH<0> t gap) = t dep)) M J Myr GLYPH<0> 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 the current Jupiter mass M J = 1 : 90 GLYPH<2> 10 30 g, respectively. We consider a 1-D viscous accretion CJD. We assume that the gas mass flux is uniform in the CJD and it is equal to both the inflow mass flux to the disk and gas accretion rate to Jupiter, and the mass of Jupiter grows by ˙ M 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 snowline r snow is at Ganymede's orbit when the gap opens and moves to Europa's orbit at the end of the disk evolution.  \nWe calculate the evolution of dust particles in the CJD using a 1-D 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 ˙ M 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 ˙ M g at the outer edge of the  \nFigure 1 . Two models of the new formation scenario for the Galilean satellites. Model A) 1 Myr) Jupiter grows to GLYPH<25> 0 : 4 M J. The gas accretion rate decreases due to a gap around the CJD and then an inner cavity around Jupiter opens. Pebbles drifting from the outer region of the CSD pile up at the pressure maximum of the gap. Only small dust particles coupled with gas are supplied to the CJD. Some planetesimals form from the pebbles at the pressure maximum. 1 : 5 Myr) Three planetesimals are captured by the CJD and migrate toward Jupiter by aerodynamic drag. The innermost one stops at the edge of the inner cavity and the two other planetesimals are captured into 2:1 mean motion resonances one by one. The position of the snow line is just inside the third satellite. 2 Myr) The fourth planetesimal is captured by the disk. It migrates inward quickly and is captured into a 2:1 resonance with the third one. The di GLYPH<11> erence in capture time between the third and fourth planetesimals creates the dichotomy of their interior ice-rock di GLYPH<11> erentiation. 30 Myr) The four planetesimals grow to the same sizes with the current Galilean satellites. The gas accretion rate decreases much because of the depletion of the parent CSD. The snowline then moves to just inside the second satellite (Europa) and small quantities of icy pebbles are accreted onto its outermost shell. GLYPH<24> 100 Myr) The CJD has disappeared and the fourth satellite (Callisto) escapes from the resonance. The rock-metal di GLYPH<11> erentiation occurs in Io by tidal heating and in Europa and Ganymede by long-lived radiogenic heating. Model B) Until the third planetesimal is captured by the CJD, Model B is the same with Model A. 10 Myr) Ganymede reaches its PIM; a gas pressure maximum forms and drifting pebbles are trapped at the maximum, resulting in the termination of the growth of Io, Europa, and Ganymede. Europa accretes a small munber of icy pebbles before the termination of growth due to inward migration of the snowline. The seed of Callisto forms at the pressure maximum out of the trapped material. 30 Myr) Callisto grows to its current size, accreting the trapped material. GLYPH<24> 100 Myr) The gas disk disappears and Callisto is scattered because its orbit is close to that of Ganymede.  \n<!-- image -->  \ndisk, but reduces inwards because of filtering by satellites and evaporation of ice. There are only rocky pebbles inside the snowline and we assume that the pebble mass flux is halved. We also assume that the fragmentation of pebbles occur when their collisional velocity becomes faster than 5 or 50 m s GLYPH<0> 1 for inside or outside the snowline, respectively (Wada et al. 2009, 2013). We calculate the disk midplane temperature as T d = (3 GM cp ˙ M g = (8 GLYPH<25>GLYPH<27> SB r 3 )) 1 = 4 g , where g = (3 = 8 GLYPH<28> + 1 = (4 : 8 GLYPH<28> )) 1 = 4 , GLYPH<27> SB and G are the StefanBoltzmann constant and the Gravitational constant, respectively (Nakamoto & Nakagawa 1994). The Rosseland mean opacity is GLYPH<28> = GLYPH<20> GLYPH<6> g, where the dust opacity is GLYPH<20> = 450 r gg for T d GLYPH<21> 160 K and GLYPH<20> = 450( T d = 160 K) 2 r gg for T d < 160 K, where r gg is the ratio of the surface density of grains a GLYPH<11> ecting the disk temperature to gas.  \nWe calculate the mass of the growing seeds as  \nM s( t ) = Z t t cap ˙ M p P e GLYPH<11> dt ; (1)  \nwhere t cap and P e GLYPH<11> are the capture time of the seeds and their pebble accretion e GLYPH<14> ciencies. The pebble accretion e GLYPH<14> ciency 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 Eq. (1). See Section 2.4.  \nWe 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 mean motion resonances (Ogihara & Kobayashi 2013). The orbit of Callisto is fixed in Model B at the gas pressure maximum created by Ganymede.  \nWecalculate 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 26 Al decay, GLYPH<1> T fin( R ), from when the pebbles (the ice mass fraction is 0.5) including 26 Al have been accreted on the  \nseeds' 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 fin( R ) = T s( R ) + GLYPH<1> T fin( R ) (Barr & Canup 2008). We do not include thermal di GLYPH<11> usion, solid-state convection, and latent heat inside the seeds.  \nBelow, we show the details of the methods.", '2.2. Circum-Jovian Disk model': "After the gap formation, the gas accretion rate to Jupiter becomes much lower. In this work, we assume that the gas accretion rate is  \n˙ M g = ˙ M g ; gap exp GLYPH<0> t GLYPH<0> t gap t dep ! ( t > t gap) : (2)  \nWe also assume that the initial gas accretion rate and the gas depletion timescale are ˙ M g ; gap = 0 : 2 M J Myr GLYPH<0> 1 and t dep = 3 GLYPH<2> 10 6 yr, where M J is the current Jupiter mass. We set t = 0 as the timing of the formation of CAIs. The mass of Jupiter grows by ˙ M 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 circum-stellar disk (CSD) to circum-Jovian disk (CJD) because we assume that the flow of gas (and pebbles) in whole region of the disk is semi-steady and all of the gas flow into the disk will be accreted by Jupiter eventually.  \nA cavity of the gas disk can open around Jupiter by the magnetic field of the planet. The timing that 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  \nr cav = 0 B B B B B @ B 4 cp R 12 cp 4 GM cp ˙ M 2 g 1 C C C C C A 1 = 7 = 1 : 07 B cp 40 Gauss ! 4 = 7 R cp R J ! 12 = 7 GLYPH<2> M cp 0 : 4 M J ! GLYPH<0> 1 = 7 0 B B B B B @ ˙ M 2 g 0 : 2 M J Myr GLYPH<0> 1 1 C C C C C A GLYPH<0> 2 = 7 [ R J] ; (3)  \nvalid 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 Gauss (Connerney 1993). Previous work, however, argued that the magnetic field was once stronger than the current one (Stevenson et al. 1983; S'anchez-Lavega 2004; Christensen et al. 2009). Therefore, we consider that the magnetic field is GLYPH<25> 40 Gauss. In this case, the disk inner cavity and gap open at almost same time (substituting the Jupiter radius R J = 7 : 15 GLYPH<2> 10 9 cm at r cav for Eq. (3)). On the other hand, if the strength of the magnetic field is the same  \nFigure 2 . Model of the gas accretion rate. The black solid curve represents the evolution of the gas accretion rate in our model. After the gap opens at t = t gap, the accretion rate decreases exponentially. The red and blue regions represent the gas accretion rates required for the gap opening if the strength of Jupiter's magnetic field is B cp GLYPH<25> 40, or 4 Gauss (i.e., the current strength), respectively. The red and blue dashed lines represent the time (after the formation of CAIs) that the inner cavity opens if B cp GLYPH<25> 40 or 4 Gauss, respectively.  \n<!-- image -->  \nwith 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 < 14Myr) and too little material would remain to form the Galilean satellites after gap opening. We summarize the evolution of ˙ M 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.  \nBased on the above gas accretion model, we calculate the evolution of the 1-D CJD. The gas surface density of the viscous accretion disk is,  \nGLYPH<6> g = ˙ M g GLYPH<10> K 3 GLYPH<25>GLYPH<11> c 2 s : (4)  \nWe assumed that the strength of turbulence is GLYPH<11> = 10 GLYPH<0> 4 because MRI is suppressed in the CJD (Fujii et al. 2014). The upper 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.  \nWe assume that the CJD is viscously heated. The gas temperature in the midplane of the viscous accretion disk is given by (Nakamoto & Nakagawa 1994),  \nT d = 3 GM cp ˙ M g 8 GLYPH<25>GLYPH<27> SB r 3 ! 1 = 4 g ; (5)  \nwhere GLYPH<27> SB is the Stefan-Boltzmann constant and,  \ng = 3 8 GLYPH<28> + 1 4 : 8 GLYPH<28> ! 1 = 4 (6)  \nis a function of the Rosseland mean optical depth GLYPH<28> = GLYPH<20> GLYPH<6> g. In principle, the opacity GLYPH<20> depends on the size distribution of solid particles. However, the size distribution cannot be  \nFigure 3 . Evolution of the circum-Jovian disk and the pebbles. The left upper and lower panels represent the evolution of the gas surface density and the temperature of the midplane, respectively. The black horizontal line is the sublimation temperature of water ice, 160 K. The right upper and lower panels represent the evolution of the Stokes number of the drifting pebbles and the pebble surface density in Model A, respectively. The color variations of the both panels represent the time after the gap opens ( t GLYPH<0> t gap).  \n<!-- image -->  \npredicted 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 a GLYPH<11> ect the temperature to the gas surface density. Then, the opacity can be assumed as  \nGLYPH<20> = 8 > > > < > > > : 450 r gg T d GLYPH<21> 160 K 450 ( T d = 160 K) 2 r gg T d < 160 K : (7)  \nThe lower left panel of Figure 3 represents the evolution of the temperature when r gg = 1 : 7 GLYPH<2> 10 GLYPH<0> 7 . We choose this value to get the thermal condition that the snowline is just inside current Ganymede's orbit. The gradients of curves change at the snowline where the temperature is 160 K. The Rosseland mean optical depth GLYPH<28> becomes the smallest ( GLYPH<24> 1) slightly outside the snowline. 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.", '2.3. Pebble Growth and Radial Drift': '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 ˙ M p does not depend on the distance from Jupiter r . It is also equal to the mass flux of dust particles supplied from  \nthe parent CSD to the CJD,  \n˙ M p = x ˙ M g ; (8)  \nwhere 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 snowline is smaller than that of outside because pebbles lose their H2O ice inside the snowline. We assume that the rock mass fraction of pebbles outside the snowline is m r = 0 : 5 so that ˙ M p inside the snowline becomes half of that of outside. We also consider the filtering e GLYPH<11> ect by outer satellites.  \nUnder 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 Eq. (15) of Shibaike et al. (2017)),  \nStp = 0 : 23 2 3 + 2 p + q ! 4 = 5 10 18 GLYPH<0> 39 q ! 2 = 5 GLYPH<2> ˙ M p = ˙ M g 0 : 003 ! 2 = 5 GLYPH<18> GLYPH<11> 10 GLYPH<0> 4 GLYPH<19> 1 = 5 GLYPH<2> GLYPH<18> T d 160 K GLYPH<19> GLYPH<0> 2 = 5 M cp 1 M J ! 2 = 5 r 10 R J ! GLYPH<0> 2 = 5 ; (9)  \nwhere p and q are the r exponents of the gas surface density and temperature (i.e., GLYPH<6> g / r GLYPH<0> p and T d / r GLYPH<0> q ).  \nHowever, fragmentation occurs if the collision velocity, in  \nother words, the pebble-to-pebble relative velocity is too fast. This relative velocity is  \nv pp = q ( v r = 2) 2 + v 2 t ; (10)  \nwhere 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),  \nv r = GLYPH<0> 2 Stp Stp 2 + 1 GLYPH<17> v K ; (11)  \nwhere v K = r GLYPH<10> K is the Kepler velocity, and,  \nv t = p 3 GLYPH<11> c sStp 1 = 2 : (12)  \nThe ratio of the pressure gradient force to the gravity of Jupiter is,  \nGLYPH<17> = GLYPH<0> 1 2 H g r ! 2 @ ln GLYPH<26> g c 2 s @ ln r ; (13)  \nwhere GLYPH<26> g = GLYPH<6> g = ( p 2 GLYPH<25> H g) is the gas density at the midplane. If the Stokes number is limited by their fragmentation, it is  \nStp = GLYPH<0> 3 GLYPH<11> c 2 s + q 9 GLYPH<11> 2 c 4 s + 4 GLYPH<17> 2 r 2 GLYPH<10> 2 K v 2 cr 2 GLYPH<17> 2 r 2 GLYPH<10> 2 K ; (14)  \nwhere v cr is the critical fragmentation speed (Okuzumi et al. 2016). This equation can be derived by substituting Eq. (11), (12), and v pp = v cr for Eq. (10).  \nFinally, the pebble surface density follows from the continuity equation,  \nGLYPH<6> p = ˙ M p 2 GLYPH<25> rv r : (15)  \nThe 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 snowline. The snowline migrates inward because the temperature becomes lower as the gas accretion rate decreases. We also find that the Stokes number inside the snowline is smaller than that outside. Outside the snowline, the Stokes number is determined by drift (Eq. (9)). While, inside the snowline, the Stokes number is determined by fragmentation (Eq. (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 GLYPH<0> 1 , 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 e GLYPH<14> ciency 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 free-fall timescale, and the stopping time of pebbles is much smaller than it because  \nthe upper panel of Figure 3 shows that the Stokes number is about 0 : 02 around the cavity (5 : 89 R J).', '2.4. Pebble Accretion E GLYPH<14> ciency': "According to recent N-body simulations, the pebble accretion e GLYPH<14> ciency is well-fitted by (Liu & Ormel 2018)  \nP e GLYPH<11> = 8 > > > < > > > : 0 B B B B B B @ 0 : 32 s GLYPH<22> s GLYPH<1> v = v K Stp GLYPH<17> 2 1 C C C C C C A GLYPH<0> 2 + 0 : 39 GLYPH<22> s GLYPH<17> h p ! GLYPH<0> 2 9 > > > = > > > ; GLYPH<0> 1 = 2 ; (16)  \nwhere GLYPH<22> 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, 2-D (the first term) and 3-D (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 di GLYPH<11> usion (Youdin & Lithwick 2007),  \nH p = H g 1 + Stp GLYPH<11> 1 + 2Stp 1 + Stp ! GLYPH<0> 1 = 2 ; (17)  \nwhere H g = c s = GLYPH<10> K is the gas scale height. The expression in the 2-D limit depends on the approach velocity of the pebbles GLYPH<1> v , which is given by the Keplerian shear in the disk or else the disk head wind. The approach velocity is  \nGLYPH<1> v = v K = 0 : 52( GLYPH<22> sStp) 1 = 3 + GLYPH<17> ( 1 + 5 : 7 GLYPH<22> s GLYPH<17> 3 = Stp !) GLYPH<0> 1 ; (18)  \nwhere the first and second terns represent the shear and the head wind limits, respectively (Ormel & Liu 2018).  \nThe pebble mass flux within an orbit of the seed is smaller than without because a fraction of drifted pebbles are captured by the embryo,  \n˙ M p ; in = (1 GLYPH<0> P e GLYPH<11> ) ˙ M p ; out ; (19)  \nwhere ˙ M p ; in and ˙ M p ; out are the fluxes inside and outside the embryo, respectively. The growth timescale of the seeds is t grow = M s = ( dM s = dt ).  \nE GLYPH<14> cient pebble accretion occurs when the satellite seed mass exceeds a critical mass (Ormel 2017),  \nM GLYPH<3> = v 3 hw t stop 8 G = 1 8 GLYPH<17> 3 Stp M cp (20)  \n(settling regime). In this work, we assume that all seeds always accrete pebbles e GLYPH<11> ectively. In our CJD model, we found that the gas aspect ratio h g is about 0.1 in the whole disk regions and so GLYPH<17> GLYPH<24> h 2 g GLYPH<24> 10 GLYPH<0> 2 . The Stokes number of the pebbles is Stp GLYPH<24> 0 : 1 (See Figure 3). The critical mass is then M GLYPH<3> GLYPH<24> 10 GLYPH<0> 8 M cp, which is about 10 times smaller than the initial mass M ini = 3 GLYPH<2> 10 23 g.  \nPebble 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 pebble isolation mass (PIM), is M iso GLYPH<24> h 3 g M cp, which is close  \nto 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's 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 GLYPH<2> 10 26 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 ( ˙ M p ; in in Eq.(19)) becomes zero immediately.  \nAt 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,  \ndM s dt = min( R col GLYPH<6> p GLYPH<10> K r 2 ; ˙ M p) ; (21)  \nwhere R col = 2 GLYPH<25>GLYPH<17> P e GLYPH<11> is the dimensionless accretion rate. Note that R col does not contain an explicit GLYPH<17> -dependence. The approach velocity of pebbles in P e GLYPH<11> , GLYPH<1> v = v K, is then calculated by Eq. (18) without the headwind term (the second term of the right-hand side). Kanagawa et al. (2018) argued that, in cases of planets in 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 GLYPH<6> p = GLYPH<6> g here but we find that the ˙ M p term in Eq.(21) nonetheless limits the growth of the satellite: as soon as pebbles arrive at the pressure bump, they are accreted.", '2.5. Satellite Internal Temperature': 'We calculate the internal temperature of satellites to estimate the levels of di GLYPH<11> erentiation. 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  \nGLYPH<26> s C p( T s( R s) GLYPH<0> T d) dR s dt = 1 2 ˙ M s u 2 i 4 GLYPH<25> R 2 s GLYPH<0> GLYPH<27> SB( T s( R s) 4 GLYPH<0> T 4 d ) : (22)  \nThis 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 form the surface. The pebble-satellite collision velocity u i can be estimated by  \nu i = min( q v 2 esc + GLYPH<1> v 2 ; v set) ; (23)  \nwhere v esc = p 2 GM s = R s and v set = g s t stop are the escape and settling velocities, respectively. The gravitational acceleration at the surface of the satellite is g s = GM s = R 2 s .  \nHeating by 26 Al is also very e GLYPH<11> ective. The increase of the satellite internal temperature at R , distance from the center of  \nthe satellite, during the formation is  \nGLYPH<1> T fin( R ) = 1 C p Z t fin t acc m r q 26( t ) dt = m r q 26 ; 0 C p GLYPH<21> 26 exp( GLYPH<0> GLYPH<21> 26 t acc) ; (24)  \nwhere 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 26 Al, the specific heat, and the satellite density as q 26 ; 0 = 1 : 82 GLYPH<2> 10 GLYPH<0> 7 Wkg GLYPH<0> 1 , C p = 1400 J kg GLYPH<0> 1 K GLYPH<0> 1 , and GLYPH<26> s = 1 : 9 g cm GLYPH<0> 3 , which are the values in the table (Barr & Canup 2008). Here, we assume that all 26 Al heat has remained at the point until the end. Thermal di GLYPH<11> usion can be ignored because it di GLYPH<11> uses only GLYPH<24> 10 km in 10 7 years. Solid-state convection can also be ignored because it can start at t GLYPH<24> 10 8 year (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 be then estimated by the sum of the two heating sources, the accretion heating and the 26 Al heating,  \nT fin( R ) = T s( R ) + GLYPH<1> T fin( R ) : (25)  \nNote that the 26 Al heating is dominant (see Section 4.8).', '2.6. Satellite Migration': '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, Sts, for Stp in Eq. (11),  \nv ad = GLYPH<0> 2 Sts Sts 2 + 1 GLYPH<17> v K : (26)  \nThe Stokes number of satellites is,  \nSts = 8 3 C D GLYPH<26> s R s GLYPH<26> g GLYPH<17> v K GLYPH<10> K : (27)  \nThe Type I migration velocity is,  \nv t1 = b t1 M s M cp ! GLYPH<6> g r 2 M cp ! v K c s ! 2 v K ; (28)  \nwhere b t1 is the migration constant, which depends on the distribution of the temperature and gas surface density of the CJD (see Eq. (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 mig = r = j v ad + v mig j .', '2.7. Capture into Mean Motion Resonances': 'According to recent N-body simulation, the critical migration timescale for capture into mean motion resonances is (Ogihara & Kobayashi 2013),  \nt crit = C MMR M in M GLYPH<8> ! GLYPH<0> 4 = 3 M cp M GLYPH<12> ! 4 = 3 T in ; (29)  \nTable 1 . Capture coe GLYPH<14> cient C MMR for 2:1 or 3:2 MMRs adopted from Ogihara & Kobayashi (2013)  \nwhere M in, M GLYPH<8> , and M GLYPH<12> 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 coe GLYPH<14> cient 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 1-D model does not consider the eccentricity and inclination of the orbits.', '3. RESULTS': 'We calculate the evolution of the mass and orbits of the satellites in the two models, Models A and B. We first show the results of Model A (Section 3.1). We then also show the results of Model B and compare them with those of Model A (Section 3.2).', '3.1. Model A': "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 ; 2 : 0 Myr, respectively. The disk condition is x = 0 : 0026, GLYPH<11> = 10 GLYPH<0> 4 , and r gg = 1 : 7 GLYPH<2> 10 GLYPH<0> 7 . We also assume the initial mass and positions of the seeds as M s ; start = 3 GLYPH<2> 10 23 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.  \nFigure 4 and Table 2 show that the mass (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 snowline is half of that of outside because icy pebbles evaporate inside the snowline.  \nThe upper panel of Figure 5 represents the pebble accretion e GLYPH<14> ciencies, P e GLYPH<11> , of Seeds A1 to A4. All pebble accretion e GLYPH<14> ciencies are smaller than GLYPH<24> 10%. The changes in the slope around M s GLYPH<25> 5 GLYPH<2> 10 23 g are caused by the stop of their migration due to the resonance traps. The mild changes in the gradients around M s GLYPH<25> 2 GLYPH<2> 10 24 g and M s GLYPH<25> 10 25 g correspond to the transition of the pebble accretion regimes from the 3-D to the 2-D and from the head wind to the shear regimes, respectively. The Stokes number inside of the snowline is smaller than that of outside by a factor of 2-3 (Figure 3), which makes the pebble accretion e GLYPH<14> ciencies of Io and Europa larger than that of Ganymede and Callisto in the 2-D  \nregime because P e GLYPH<11> is a decreasing function (see Eqs. (16), (18)). The lower 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 they passing the snowline, which has more impact than the reduction in St in the 3-D and the 2-D-shear regimes (see Eqs. (1), (16), (17), and (18)). The small jump around M s GLYPH<25> 5 GLYPH<2> 10 25 g on Europa's curve is caused by the crossing of the snowline. The filtering e GLYPH<11> ects are mildly important for the growth timescale. Note that their mass does not reach the pebble isolation mass, M iso (Eq. (33)).  \nFigure 4 also shows that all seeds migrate quickly ( < 3 GLYPH<2> 10 5 year) 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, di GLYPH<11> ers 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).  \nEach Galilean satellite has the di GLYPH<11> erent 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 GLYPH<0> 9%) can be reproduced by the migration of the snowline 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 which is suitable for reproducing the ice mass fractions of all the Galilean satellites. In our model ice sublimation occurs instantaneously at the snowline, 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.  \nWe find that, in order to avoid di GLYPH<11> erentiation of Callisto by 26 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 26 Al heat. The dichotomy of their internal ice-rock differentiation can be created by the di GLYPH<11> erence in their capture time, 0 : 5 Myr, because the half-life of 26 Al is 0 : 717 Myr. The long growth timescale ( GLYPH<24> 10 7 yr) is the reason why such di GLYPH<11> erent capture time is allowed. If the growth timescale is shorter and the di GLYPH<11> erence in the capture time is the same, the final mass of Ganymede and Callisto would end up too large. Indeed, it is di GLYPH<14> cult to make the dichotomy between the internal structures of the two satellites by the classical  \nCanup-Ward formation scenario where the growth timescale is < 10 6 yr (Barr & Canup 2008). Note that once rock-metal cores form in Europa and Ganymede, they can also di GLYPH<11> erentiate metallic cores by long-lived radiogenic heating (Spohn & Breuer 1998), and Io can be di GLYPH<11> erentiated 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.", '3.2. Model B': "We then calculate the evolution of the satellites in a situation that 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.  \nFigure 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 e GLYPH<11> ect of Callisto. In Model A, Seed A4 starts to decrease the pebble mass flux inside its orbit at t GLYPH<0> t gap = 1 Myr. Conversely, Seed B4 appears only at t GLYPH<0> 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.  \nFigure 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 = (3 M J)) 1 = 3 r 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 p 3 r H ; 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).  \nWe 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  \nthat in Model A but not much smaller than the actual value, 6 GLYPH<0> 9% (Kuskov & Kronrod 2005).  \nThe 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 26 Al radiogenic decay does not heat the satellite since Seed B4 starts to grow at later time (9 : 54 Myr), which is naturally achieved in Model B. Moreover, the di GLYPH<11> erence 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.", '4. ASSESSMENT': 'Although our new slow-pebble-accretion scenario reproduces most of the important properties of the current Galilean satellites - mass, orbits, ice mass fractions, and internal structures - many parameters needed to be tuned in order to meet the constraints. We summarize the key assumptions in Table 3 and discuss each of them in the following paragraphs. The first column lists the key assumptions of our model and the second and third columns describe the motivations for them. Before discussing the sensitivity of the results to these parameters in detail, we would like to remark that the many constraints available for the Jovian system necessarily compels us to adopt specific parameter choices. Nonetheless, most of assumptions are supported by (or not inconsistent with) previous predictions or observations, and (as we will see below) our model would cope with a modest level of parameter variation. Our scenario, for the first time, reproduces most of the characteristics of the Galilean satellites simultaneously and consistently. Here, we show the detailed assessment of each of the assumptions.', '4.1. Early Formation of Jupiter': "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 26 Al heat. This is consistent with a resent 'early Jupiter formation' scenario showing that solid materials in the solar system were spatially separated by that time (Kruijer et al. 2017).", '4.2. Gas Inflow to the CJD': '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, ˙ M g, decreases as an exponential decay with the  \n0.6  \non  \nti  \nfrac  \ns  \nas  \nm  \ne  \nc  \nI  \nFigure 4 . Evolution of the four satellites (Model A). The solid curves represent the positions of the evolving seeds (A1 to A4) at the time after the gap opens (i.e., t GLYPH<0> t gap). The sizes of the circles represent the radii of the seeds and the current Galilean satellites. The color scales of the curves range from gray to dark blue for the increasing ice mass fractions of the seeds. The current ice mass fraction of the satellites (the mean values of the estimates by Kuskov & Kronrod (2005)) are also shown as the color scale of the circles and the white lines in the column. The blue dashed curve represents the position of the snowline. The position of the edge of the inner cavity is fixed at the current orbit of Io (the vertical dotted line). The horizontal dotted line represents the gas depletion timescale of the CSD. The values of the final mass, positions, and ice mass fractions are shown in Table 2.  \n<!-- image -->  \n0  \nFigure 5 . Pebble accretion e GLYPH<14> ciency and growth timescale (Model A). The purple, green, light blue, and orange curves in the upper panel represent the pebble accretion e GLYPH<14> ciencies of Seeds A1 to A4, respectively. The lower panel is the growth timescale of the seeds.  \n<!-- image -->  \ntimescale 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 that, however, 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,b). In this case, the migration of snowline 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.', '4.3. Inner Cavity': "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'anchezLavega 2004; Christensen et al. 2009). The position of the  \nedge of the cavity can be estimated from the balance of the gas accretion rate and magnetic stress by Jupiter's magnetic field (Eq. (3)). Since the gas accretion keeps decreasing after the gap opens, the edge moves outward and the innermost seed (Io) may also migrates 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 further than the current, r co = 2 : 25 R J.  \nFirst, we check if 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,  \nGLYPH<0> 2s = 1 2 M s v K v t1 : (30)  \nOnthe 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 onesided corotation and Lindblad torque is,  \nGLYPH<0> 1s ; co = 2 : 46 GLYPH<6> g r 2 M cp ! 0 B B B B @ GLYPH<22> s h 3 g 1 C C C C A 1 = 2 M s( r GLYPH<10> K) 2 ; (31)  \nand  \nGLYPH<0> 1s ; Lin = GLYPH<0> 0 : 65 GLYPH<6> g r 2 M cp ! GLYPH<22> s h 3 g M s( r GLYPH<10> K) 2 ; (32)  \nrespectively (Liu et al. 2017). Figure 8 shows that the onesided corotation torque Io (Seed A1) receives is much larger than the negative one-sided Lindblad torque and the negative two-sided torque of 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 : 5 t 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 that, however, the timescale of Io's migration by the one-sided corotation torque is GLYPH<24> 10 3 year and this is much shorter than the migration timescale of the edge, GLYPH<24> 10 6 year. Therefore, Io may have moved much further than its current orbit with the disk inner edge moves outward.  \nHowever, this outward migration of the inner edge should stop at the corotation radius r co where the Keplerian frequency of the disk equals to 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 GLYPH<25> 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  \non  \nti  \nfrac  \ns  \nas  \nm  \ne  \nc  \nI  \nFigure 6 . Evolution of the four satellites (Model B). Same as Figure 4 but the evolution of Seeds B1 to B4 are shown. The time that Seed B3 reaches its PIM is t iso ; G = 9 : 54 Myr after the gap opening. The values of the final mass, positions, and ice mass fractions are shown in Table 2.  \n<!-- image -->  \nTable 2 . Final properties of the satellites  \n- 1 by definition\n- 2 Estimated current ice mass fractions by Kuskov & Kronrod (2005).  \nFigure 7 . Internal temperature of Ganymede and Callisto. The light blue and orange curves represent the final ( t GLYPH<0> t gap = 30 Myr) internal temperature of Ganymede (Seeds A3 and B3) and Callisto (Seed A4 and B4), respectively. The solid and dashed curves are the results of Model A and B, respectively. The black curve is the melting point of Callisto cited by Figure 5 in Barr & Canup (2008).  \n<!-- image -->  \nFigure 8 . Torque which the satellites receive (Model A). The solid and dashed purple curves are the Lindblad and corotation one-sided torque that Io (Seed A1) receives, respectively. The green, light blue, and orange solid curves are the two-sided torque of Europa (Seed A2), Ganymede (Seed A3), and Callisto (Seed A4), respectively. The black curve is the sum of the negative torque; the Lindblad one-sided torque of Io and the two-sided torque of the other three satellites.  \n<!-- image -->  \nlarger 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 GLYPH<25> 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 GLYPH<25> 4 : 7 R J. We can then consider two scenarios of Io formation. In the first one, Io formed around r GLYPH<25> 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, espe-  \ncially 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 snowline 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 that the corotation radius when the radius of Jupiter is GLYPH<25> 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).  \nThe ionization degree of the disk inner region, GLYPH<31> e GLYPH<17> n e = n n, where n e and n n are respectively number density of electron and neutral gas, is also important. Without enough ionization, the disk gas can not 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 GLYPH<17> VH g =GLYPH<21> > 1, where V is the relative velocity between the rotating magnetic field and the disk gas, and GLYPH<21> is the magnetic di GLYPH<11> usivity (Takata & Stevenson 1996). We note that this magnetic Reynolds number is GLYPH<24> 10 3 times larger than the Elsasser number, the critical number for the MRI activation, GLYPH<3> = v 2 Az = ( GLYPH<21> GLYPH<10> K), where v Az is the z component of the Alfv'en velocity. If Jupiter and its magnetic field spin rigidly at the current speed, V GLYPH<25> 3 GLYPH<2> 10 4 ms GLYPH<0> 1 , H g GLYPH<25> 1 : 1 GLYPH<2> 10 7 m, and GLYPH<21> GLYPH<25> 0 : 74 =GLYPH<31> e at r = 1 R J with the disk temperature is T = 1000 K (Blaes & Balbus 1994). The condition for angular momentum transport is then GLYPH<31> e & 10 GLYPH<0> 12 GLYPH<0> 10 GLYPH<0> 11 . Takata & Stevenson (1996) estimated the ionization degree of the circum-Jovian disk including the e GLYPH<11> ects of galactic cosmic rays and 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 GLYPH<31> e . 10 GLYPH<0> 16 at the midplane, GLYPH<31> e GLYPH<24> 10 GLYPH<0> 15 at the altitude of one scale height, and GLYPH<31> e GLYPH<24> 10 GLYPH<0> 13 GLYPH<0> 10 12 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 GLYPH<31> e . 10 GLYPH<0> 14 at the midplane, GLYPH<31> e GLYPH<24> 10 GLYPH<0> 13 at the altitude of one scale height, and GLYPH<31> e GLYPH<24> 10 GLYPH<0> 11 GLYPH<0> 10 GLYPH<0> 10 at three scale heights. Therefore, the condition for angular momentum transport can be achieved at the 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  \nTable 3 . Summary of the key assumptions and their validity  \n- 1 Only in Model A  \ndrift 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 future.  \nMoreover, 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.", '4.4. Strength of Turbulence': 'We have also assumed that the strength of turbulence of the viscous accretion disk is GLYPH<11> = 10 GLYPH<0> 4 , 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 a GLYPH<11> ects many properties of the system, for example, the disk temperature, the pebble accretion rates, and the migration speeds.', '4.5. Disk Temperature': "To get the disk thermal condition suitable for reproducing the ice mass fractions of the satellites, we needed that Rosseland mean optical depth GLYPH<28> is around unity at Ganymede's orbit. Figure 9 represents the r gg dependences of the snowline, satellite mass, and ice fraction of satellites in the case of the fixed orbits. The top panel shows that when r gg = 1 : 7 GLYPH<2> 10 GLYPH<0> 7 , the snowline 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 Eq. (6)). If r gg is larger than this, the first term is dominant until the almost 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 GLYPH<2> 10 GLYPH<0> 7 . However, if the gas accretion rate ˙ M g is smaller than this assumption and the dust-to-gas accretion rate ratio x is correspondingly larger, the snowline 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 snowline 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 can not get enough ice. These results can be understood by the top panel.  \nHowever, 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 GLYPH<20> p GLYPH<24> 1 = ( GLYPH<26> int ; p a p), where GLYPH<26> int ; p and a p are the pebble internal density and the pebble radius, respectively. If the pebbles are flu GLYPH<11> y, the opacity of pebbles could be GLYPH<20> p GLYPH<24> 1 = (10 GLYPH<0> 3 GLYPH<2> 10) = 10 2 cm 2 g GLYPH<0> 1 (Kataoka et al. 2014;  \nShibaike et al. 2017). In this case, the mean optical depth of pebbles is GLYPH<28> p = GLYPH<20> p GLYPH<6> p GLYPH<24> 10 2 GLYPH<2> 10 GLYPH<0> 2 = 1 (see Figure 3). This is consistent with the Rosseland mean optical depth in our models, GLYPH<28> = 450 r gg GLYPH<6> g GLYPH<24> 500 GLYPH<2> (2 GLYPH<2> 10 GLYPH<0> 7 ) GLYPH<2> 10 4 GLYPH<24> 1 (see Eq. (7)). Note that the gas (molecular) opacity ( GLYPH<24> 10 GLYPH<0> 5 GLYPH<0> 10 GLYPH<0> 4 cm 2 g GLYPH<0> 1 , consistent with an optical depth of GLYPH<24> 0 : 1 GLYPH<0> 1 with GLYPH<6> g GLYPH<24> 10 4 g 1 cm GLYPH<0> 2 ) may also a GLYPH<11> ect the disk temperature but the e GLYPH<11> ect should be limited (Mizuno 1980). Second, in reality Eqs.(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.  \nNote that, we simply assumed the snowline as the position where the disk temperature is 160 K after previous models of PPDs, the gas pressure of the CJD is much higher than that of PPDs, the partial water vapor pressure and the condensation temperature should also be higher. In this case, the evolution track of snowline 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 chose a larger ˙ M g but it means a shorter formation timescale and may make the di GLYPH<11> erence between the internal structures of Ganymede and Callisto smaller.", '4.6. A Little Supply of Solid Material to the CJD': '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 ( GLYPH<24> 0 : 1 mm) can form as fragments of the collisions of the solids trapped at the gas pressure maximum (Kobayashi et al. 2012).  \nFurthermore, 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 which have high eccentricities in the CSD (Fujita et al. 2013); planetesimals can form at the  \n2  \nFigure 9 . E GLYPH<11> ects of the thermal condition of the disk. The top, middle, and bottom panels represent the evolution of the snowline, the final satellite mass, and the final ice mass fractions of satellites, respectively. The color variations represent the di GLYPH<11> erence in the grain-to gas surface density ratio. The black lines and diamonds are the current properties of the Galilean satellites (Kuskov & Kronrod 2005). We put seeds of satellites on 1 to 100 R J with the initial mass of M s ; start = 3 GLYPH<2> 10 23 g and fix their positions in this calculation. The dust-to-gas accretion ratio is x = 0 : 0021.  \n<!-- image -->  \ngas 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).', '4.7. Properties of the Captured Planetesimals': 'We consider the cases that the planetesimals which 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 that planetesimals are captured by multiple approaches to Jupiter and become circu-  \nlarized around GLYPH<25> 50 R J (long-lived prograde captured orbits) (Suetsugu et al. 2016).  \nWe also fix the mass of the captured planetesimals as 3 GLYPH<2> 10 23 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 e GLYPH<14> cient pebble accretion (i.e., in the settling regime) for the adopted disk properties (see Eq. (20)). Smaller planetesimals, on the other hand, will stay small and may be scattered or accreted by large ones. In the settling regime, any di GLYPH<11> erence in the initial mass does not change the final mass because the growth timescale increases with mass (see the lower 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 Eqs. (26), (29), and Table 1). Therefore, the initial seed mass must be similar or larger than our standard value.', '4.8. Capture-time of the Planetesimals': "Each capture time of the planetesimals t cap are assumed to make the dichotomy of the internal structures of the Galilean satellites; Callisto is minimally or only modestly di GLYPH<11> erentiated but the other satellites are fully di GLYPH<11> erentiated. It is difficult to make the dichotomy between the internal structures of Ganymede and Callisto by the di GLYPH<11> erence in the heat of accretion or 26 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 di GLYPH<11> erentiated by 26 Al heat provided the seed of Callisto was captured by the disk at a su GLYPH<14> ciently late time.  \nFigure 10 represents the internal temperature of Callisto where its seed is captured at di GLYPH<11> erent timing. We fix the position of Callisto at the current orbit in this calculation. Without the 26 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 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).  \nOnly in the cases where 26 Al radiogenic heating is considered, the internal temperature exceeds the melting point of water ice. Figure 10 shows that this condition is achieved if most of the 26 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 time of the seeds of Ganymede and Callisto were 1 : 5 and 2 : 0 Myr after the CAI formation, respectively. In Model B, Callisto's  \nFigure 10 . Internal temperature of satellites. The color variations represent the di GLYPH<11> erence in the gap opening time, in other words, the time that the pebble accretion starts. The black solid curve represents the internal temperature without 26 Al decay heating (only accretion heating). The dotted black curve is the melting point of Callisto (Barr & Canup 2008). We put seeds of satellites on 26 : 3 R J (the current orbit of Callisto) at t = 0 ; 1 ; 2 ; 3 Myr with the initial mass of M s ; start = 3 GLYPH<2> 10 23 g and fix their positions in this calculation. The dust-to-gas accretion ratio is x = 0 : 0021.  \n<!-- image -->  \nseed formed at the gas pressure bump made by Ganymede at 10 : 5 Myr, when Ganymede reaches its PIM.  \nWe also note that if there is less 26 Al in pebbles (Larsen et al. 2016), the internal temperature of satellites must be lower. The curve of the internal temperature in R sT s space should approach the black curve (w / o 26 Al) in Figure 10. In that case, not only Callisto but also Ganymede should be undi GLYPH<11> erentiated during their formation.", '4.9. Fragmentation of Pebbles': '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 of 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 GLYPH<0> 1 , respectively. As a result, the size (i.e., Stokes number) of pebbles is determined by fragmentation inside the snowline although drift determines it outside the snowline (see Section 2.3). In this work, we have adopted the critical fragmentation speeds from numerical studies of collisions of pebbles with 0 : 1 GLYPH<22> 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.  \nWe note that there is still uncertainty about the fragmentation threshold. A recent experimental work argued that the critical speeds are slower than our assumptions, GLYPH<24> 1 m s GLYPH<0> 1 for both the icy and rocky pebbles (Musiolik & Wurm 2019).  \nOn the other hand, Kimura et al. (2015) and Steinpilz et al. (2019) claimed that the critical fragmentation speed is faster than our fiducial value, GLYPH<24> 10 m s GLYPH<0> 1 for silicate particles. The upper 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 snowline but also outside. On the other hand, if v cr = 10 m s GLYPH<0> 1 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 lower 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 GLYPH<0> 1 (dashed purple) is similar to the fiducial case when GLYPH<11> = 10 GLYPH<0> 5 , r gg = 1 : 7 GLYPH<2> 10 GLYPH<0> 8 , and x = 0 : 0011. The final mass with the higher critical speeds (yellow) is also similar to the fiducial case. Therefore, the e GLYPH<11> ect of the di GLYPH<11> erence of the critical fragmentation speed on our formation scenario is limited.', '4.10. Pebble Isolation Mass': 'In Model B, we define that the PIM of Ganymede is its actual current mass, 1 : 48 GLYPH<2> 10 26 g. Here, we show this assumption is plausible. An expression of the PIM estimated by 2-D simulations Ataiee et al. (2018) is,  \nM iso = h 3 g p 37 : 3 GLYPH<11> + 0 : 01 GLYPH<2> 8 > > > < > > > : 1 + 0 : 2 0 B B B B B B @ p GLYPH<11> h g s 1 Stp 2 + 4 1 C C C C C C A 0 : 7 9 > > > = > > > ; M cp : (33)  \nWe find that this mass is 7 : 65 GLYPH<2> 10 25 g with the parameters of Model B, which is about 1 : 9 times smaller than the mass of Ganymede. The PIM derived from 2-D simulations can be about 1.5-2 times smaller than the expression by 3-D simulations and the di GLYPH<11> erence can become larger with weak turbulence at least in protoplanetary disks (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.', "5.1. Subsequent Evolution of Ganymede's Internal Structure": 'We conclude that Ganymede\'s interior di GLYPH<11> erentiates into the rocky core and the icy mantle because of the 26 Al heat in the accreted pebbles. However, strictly speaking, it can melt only the region of 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 (undi GLYPH<11> erentiated) icy-rocky crust and the pure icy mantle by the Rayleigh-Taylor (R-T) instability. Here, we consider a  \nFigure 11 . E GLYPH<11> ects of the fragmentation of pebbles. The upper and lower panels represent the Stokes number of pebbles ( t = t gap) and the final mass of satellites, respectively. The purple, light blue, and yelow curves are those with v cr = 1 m s GLYPH<0> 1 (rocky and icy pebbles), v cr = 5 m s GLYPH<0> 1 (rocky) and 50 m s GLYPH<0> 1 (icy), and v cr = 10 m s GLYPH<0> 1 (rocky) and 50 m s GLYPH<0> 1 (icy), respectively. We put seeds of satellites on 1 to 100 R J with the initial mass of M s ; start = 3 GLYPH<2> 10 23 g and fix their positions in this calculation. The dust-to-gas accretion ratio is x = 0 : 0021. The dashed purple curves are the results with v cr = 1 m s GLYPH<0> 1 (rocky and icy), GLYPH<11> = 10 GLYPH<0> 5 , r gg = 1 : 7 GLYPH<2> 10 GLYPH<0> 8 , and x = 0 : 0011.  \n<!-- image -->  \nthree-layered interior structure: (1) the pristine (undi GLYPH<11> erentiated) icy-rocky crust (1000 GLYPH<0> 2600 km from the center), (2) the pure icy mantle (600 GLYPH<0> 1000 km), and (3) the rocky core (0 GLYPH<0> 600 km), as the initial condition. We estimate the depth of the layers by assuming the density of the pure ice, rocky core, and pristine crust is GLYPH<26> ice = 1 : 4 g cm GLYPH<0> 3 , GLYPH<26> rock = 3 : 0 g cm GLYPH<0> 3 and GLYPH<26> s = 1 : 9 g cm GLYPH<0> 3 , respectively.  \nFirst, we check whether the Rayleigh-Taylor 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),  \nGLYPH<17> crit = " ( n GLYPH<0> 1) 1 = n C L GLYPH<1> GLYPH<26> 2 n # GLYPH<18> Z 0 L GLYPH<19> ( n GLYPH<0> 1) = n GLYPH<1> GLYPH<26> g s( R s) Lt RT ; (34)  \nwhere n , C L GLYPH<1> GLYPH<26> , and Z 0 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 GLYPH<1> GLYPH<26> = 0 : 76, and  \nZ 0 = 1 km (see Rubin et al. (2014)). The di GLYPH<11> erence of density is GLYPH<1> GLYPH<26> = GLYPH<26> s GLYPH<0> GLYPH<26> ice = 500 kg m GLYPH<0> 3 and g s( R s) = 0 : 54 m s GLYPH<0> 1 is the gravitational field of satellite at R s( = 1000 km), distance from the center. The lengthscale over which the viscosity varies significantly is  \nL = nRT 0 E a T 0 j dT up = dz j ; (35)  \nwhere R , T 0, and E a are the gas constant, the temperature of at the interface between the two layers, and the activation energy, respectively. We assume them as T 0 = 370 K and E a = 49 kJ mol GLYPH<0> 1 (see Figure 7 and Rubin et al. (2014)). According to Figure 7, the temperature gradient across the upper crust is T 0 = j dT up = dz j = 0 : 14 K km GLYPH<0> 1 , where z is the vertical distance from the interface. We can then calculate the critical viscosity, GLYPH<17> crit = 7 : 6 GLYPH<2> 10 13 ( 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  \nGLYPH<17> up = GLYPH<17> ref exp " A T ref T up GLYPH<0> 1 !# ; (36)  \nwhere GLYPH<17> 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 coe GLYPH<14> cient A is about 20 GLYPH<0> 25 and the reference viscosity of pure ice I and pure ice V and VI, are GLYPH<17> ref GLYPH<24> 10 14 [Pa s] and GLYPH<24> 10 16 GLYPH<0> 10 17 [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 Eq. (36) and equate GLYPH<17> up with GLYPH<17> crit, the timescale of the R-T instability is 10 9 GLYPH<0> 10 10 year for A = 20 and 10 11 GLYPH<0> 10 12 year for A = 25. However, if we substitute T up = 170 K, only 20 K larger, the timescale becomes 10 GLYPH<0> 100 times shorter. Therefore, if the upper crust was heated up 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 Gaymede within the age of the satellite. One interesting observational fact of Ganymede is that there are two regions on the surface of the satellite: very primitive region ( GLYPH<24> 4 Ga) and relatively newer region ( GLYPH<24> 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.', '5.2. In-situ Formation of the Seeds': 'Finally, we highlight the possibility that the seeds of the satellites are not supplied from the CSD but form in-situ around the snowline of the CJD by streaming instability  \n(Schoonenberg & Ormel 2017). Ormel et al. (2017) and Schoonenberg et al. (2019) have used this idea towards understanding the properties of the TRAPPIST-1 system. In their model planets consecutively form at the snowline, then migrate inwards. 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 planetesimalcapture 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 snowline migrate inward and are trapped into 2:1 MMRs one by one, the orbits of the third seed must be around the snowline and no more seeds can form there.  \nHowever the main problem with this scenario is that in-situ satellitesimal formation requires dust-to-gas ratios GLYPH<24> 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 35 at the snowline Schoonenberg & Ormel (2017), would not be enough by far.', '6. CONCLUSIONS': "We developed a new formation scenario for the Galilean satellites using the pebbles drifting toward Jupiter in the circum-Jovian disk. As the seeds of the satellites, we assumed that only several planetesimals had been captured at given timing. Such conditions should be easier to occur than those assumed in previous satellitesimal-accretion scenarios because it is considered that dust particles in a circumplanetary disk drift toward the central planet before they grow to satellitesimals and a gap structure of planetesimals prevents them from being captured by the disk (Fujita et al. 2013; Shibaike et al. 2017).  \nWe first calculated the evolution of the CJD, a gas accretion disk, which was determined by the given gas inflow mass flux reducing exponentially, and the fixed strength of turbulent viscosity. We then calculated 1-D (radial distribution) representative-size evolution of the pebbles including their collisional growth, aerodynamic drag, fragmentation and the distinction of icy / rocky pebbles by the position of snowline. Finally, the growth by accretion of the pebbles and 1-D radial orbital evolution of the seeds including aerodynamic drag, Type I migration, and simple checks of resonance capture were calculated. We simultaneously calculated their internal thermal evolution by pebble accretion and 26 Al decay heating, and the pebble filtering e GLYPH<11> ects by outer satellite seeds.  \nIn contrast to the previous scenarios which only explain parts of the characteristics of the Galilean system and the sce-  \nnarios are inconsistent with each other, we found that our new scenario can reproduce the following characteristics simultaneously and consistently with possible assumptions. First, it can reproduce the mass distribution of all the Galilean satellites even in the case of a very small amount of material supply to the CJD. Second, Io, Europa, and Ganymede are captured into 2:1 resonances one by one because the inner cavity opens by the strong magnetic field of Jupiter and halts the migration of Io at the edge of the cavity. Third, Europa accretes small quantities of icy particles in the final phase of its formation because the snowline moves inward as the gas accretion rate onto the CJD decreases. Therefore, Europa's rocky core and icy mantle are explained naturally. The ice mass fractions of the three other satellites are also reproduced because the orbits are inside, respectively outside, the snowline. Finally, our model explains why only Callisto stays undi GLYPH<11> erentiated and why the other satellites melt by radiogenic heating of 26 Al decay. The di GLYPH<11> erence in the capture time of the planetesimals a GLYPH<11> ects the total 26 Al heat they got and their internal structures but not their final mass because of their slow growth.  \nWe also considered a model that Ganymede reaches its pebble isolation mass (PIM) and Callisto forms by the solid material trapped at the gas pressure maximum created by Ganymede. In this case, the evolution of the inner three satellites are almost same with the first model and the mass of Callisto is also reproduced. The orbits of Ganymede and Callisto are so close that Callisto could be scattered outward after the CJD disappears, resulting in Callisto's unique orbit not trapped in any resonance. The unmelted Callisto is also reproduced and the dichotomy between the internal temperatures of the two satellites is larger than that in the first model. This is because the interior of Callisto is only heated up by (weak) accretion heating and 26 Al heat is negligible due to the late start of growth.  \nWe thank Yasuhito Sekine for valuable discussion. This work is supported by JSPS KAKENHI JP16J09590, 15H02065 and 16K17661, and MEXT KAKENHI 18H05438. C.W.O is supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). Part of this research was initiated during his stay as a visiting Professor at the Tokyo Institute of Technology. This work has been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). Y.S. acknowledges the financial support of the SNSF.", 'A. GAS PRESSURE MAXIMUM FORMED BY GANYMEDE': 'In Model B, we considered that Callisto formed at the gas pressure maximum created by Ganymede. The local gas surface density around Ganymede can be given by (Kobayashi et al. 2012)  \nGLYPH<6> g ; local = GLYPH<6> g exp[ GLYPH<0> ( r GLYPH<0> r G = l ) GLYPH<0> 3 ] ; (A1)  \nwhere r G = 14 : 8 R J is the orbital radius of Ganymede (in Model B) and  \nl = 2 6 6 6 6 6 4 8 81 GLYPH<25> r 2 G GLYPH<10> K( r G) GLYPH<23> 0 B B B B @ M G M 2 cp 1 C C C C A 2 3 7 7 7 7 7 5 1 = 3 r G ; (A2)  \nwhere M G = 1 : 48 GLYPH<2> 10 26 g is the mass of Ganymede (in Model B). Then the local ratio of the pressure gradient and gravity is  \nGLYPH<17> local = c 2 s 2 r 2 GLYPH<10> 2 K " p + q 2 + 3 2 GLYPH<0> 3 l 3 r ( r GLYPH<0> r G) 4 # ; (A3)  \nwhere p and q are the r exponents of the gas surface density and temperature (i.e., GLYPH<6> g / r GLYPH<0> p and T d / r GLYPH<0> q ). Note that in Kobayashi et al. (2012), q is defined as the r exponents of the sound speed (see the paper for further explanations about the above equations). At the pressure maximum, GLYPH<17> local = 0, so we solve this equation and gain the position of the pressure maximum. When T = 150 K and M cp = M J, and the other parameters are consistent with those in Model B, we find that the pressure maximum is at r = 17 : 0 R J and it is almost independent on these parameters. Therefore, in Model B, we fix the orbital radius of Seed B4 as r = 17 : 0 R J.', 'B. HEAT RELEASED BY R-T INSTABILITY': 'One possible heat source to melt the whole interior of the satellite is the potential energy released during the overturn of the crusts by the R-T instability. We roughly estimate the heat released by the instability and how much the internal temperature rises. If the pure ice mantle (consistent with the half mass of the di GLYPH<11> erentiated region, M ice( = M rock) = 4 GLYPH<2> 10 21 kg) is lifted up to the surface of the satellite, the thickness of the pure ice crust will be h ice = 120 km. The mass of the pure ice (and rock) is much smaller than the whole mass of Ganymede ( M G), the increase in the temperature of the undi GLYPH<11> erentiated pristine crust uniformly heated up by the released potential energy ( GLYPH<1> Q RT) can be roughly estimated by the following equations,  \nGLYPH<1> T RT = GLYPH<1> Q RT C p M G GLYPH<25> GM ice C p( R G GLYPH<0> h ice) GLYPH<26> s GLYPH<26> ice GLYPH<0> 1 ! 0 B B B B @ 1 GLYPH<0> R 2 rock GLYPH<26> rock ( R G GLYPH<0> h ice) 2 GLYPH<26> s 1 C C C C A GLYPH<25> 20 [K] ; (B4)  \nwhere R G is the radius of Ganymede. This temperature increase is, unfortunately, not enough for the whole di GLYPH<11> erentiation of the interior of the satellite. To melt the whole region, about 100 K increase of the temperature would be needed (see Figure 7). However, we neglect the two kinds of heat provided by (1) the released potential energy during the di GLYPH<11> erentiation of the region of 1000 km from the center and (2) the central captured planetesimal which has been heated up by 26 Al decay (Wakita & Sekiya 2011). Such heat may be transported with the lifting up pure icy mantle and so be able to contribute to di GLYPH<11> erentiate the region remaining unmelted.', 'B.1. Parameters, Constants, and Variables': 'We summarize the parameters and constants used in this work in Table B1. We also summarize the valuables in Table B2 and B3.  \nTable B1 . Parameters and constants  \nTable B2 . Variables  \nTable B3 . Variables', 'REFERENCES': '- Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, The Astrophysical Journal, 747, 47\n- Tanigawa, T., & Tanaka, H. 2016, The Astrophysical Journal, 823, 48 Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62\n- Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ,\n- 702, 1490  \nWakita, S., & Sekiya, M. 2011, Earth, planets and space, 63, 1193 Weidenschilling, S. 1977, Icarus, 44, 172  \n- Yang, C. C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80  \nYoder, C. F., & Peale, S. J. 1981, Icarus, 47, 1 Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588  \n- Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, The Astrophysical Journal, 755, 6'}