© ESO 2019

1 Introduction

The gaseous planets in the solar system, Jupiter and Saturn, have large regular satellites whose radii are larger than 1000 km, and their orbits are on the equatorial plane of the central planets. It is widely accepted that the gaseous planets were surrounded by a circumplanetary disk and that large satellites formed in the disks (e.g., Lunine & Stevenson 1982; Mosqueira & Estrada 2003; Canup & Ward 2006). Several hydrodynamic simulations revealed that circumplanetary disks can be formed as byproducts of the formation and growth of gaseous planets via mass accretion from the parent protoplanetary disk (e.g., Machida et al. 2010; Tanigawa et al. 2012).

Canup & Ward (2002) proposed the gas-starved disk model whose disk mass is regulated by the steady-state viscous accretion onto the central planet. Canup & Ward (2006) then performed N-body simulations of satellite formation in the gas-starved disk and revealed that the total mass of the satellites in circumplanetary disks is successfully reproduced. Moreover, Sasaki et al. (2010) found that the number, masses, and orbits of the Galilean satellites can be reproduced for the formation era when the disk had an inner cavity.

Jupiter has four large regular satellites called the Galilean satellites: Io, Europa, Ganymede, and Callisto. Saturn, in contrast, has only one large regular satellite, Titan, and many 100–1000 km-sized satellites. The inner three Galilean satellites orbit in a 4:2:1 mean motion resonance. One of the mechanisms for reproducing this orbital configuration is the stopping of the migration of Io near the bump in the surface density distribution and the subsequent resonant trapping of Europa and Ganymede (e.g., Sasaki et al. 2010; Fujii et al. 2017).

There have been two hypotheses to explain the formation of the bump in the surface density distribution of circumplanetary disks. Fujii et al. (2017) revealed that in the early stage of disk evolution, the disk temperature is high and the sublimation of silicate dust particles occurs within a certain orbital radius. This absence of dust particles changes the opacity (e.g., Bell & Lin 1994), causing the formation of the bump in the surface density distribution of the circumjovian disk. One problem in this scenario is whether the satellites can survive against the radial drift in the late stage of the circumjovian disk when the temperature of the disk edge is sufficiently lower than the sublimation temperature. Fujii et al. (2017) noted that the fate of satellites depends on the time evolution of the mass accretion rate onto Jupiter after the formation of planetary gas in the solar nebula (e.g., Dürmann & Kley 2015; Kanagawa et al. 2018). Another problem is that icy satellites, at a minimum, cannot survive longer than the disk dissipation timescale in such a high-temperature disk, because the temperature at the current location of icy satellites is about 1000 K at this stage (see Appendix A). We found that icy satellites at the ambient temperature of 500 K might be evaporated within 0.1 Myr, even if the satellite is as large as Ganymede.

On the other hand, Sasaki et al. (2010) proposed that the inner cavity of the circumjovian disk could be formed by the magnetic coupling between the disk and young Jupiter. This magnetic coupling can also explain why the spin rate of Jupiter is sufficiently slower than the break-up spin rate (e.g., Takata & Stevenson 1996; Batygin 2018). Nevertheless, one problem exists for this scenario: the spin period of young Jupiter should be longer than 40 h at the end of the Galilean satellite formation, as the innermost satellite should have formed at the location of the inner cavity (Ogihara & Ida 2012; Shibaike et al. 2019). This spin period is four times longer than the current spin period (approximately 10 h). Therefore, the planetary radius of young Jupiter may be at least twice as large as the current planetary radius (see Appendix B); although several evolutionary calculations have indicated that the planetary radius of young Jupiter might be approximately 1.5 times the current radius at 1 Myr after the formation of the planet (e.g., Fortney et al. 2011).

In this study, we propose another mechanism to create a bump in the surface density distribution of the circumjovian disk – the depletion of fine dust particles by photophoresis. Photophoresis is the momentum transfer mechanism whereby gas molecules collide with photo-irradiatied dust particles that exhibit a temperature gradient (e.g., Krauss & Wurm 2005; Loesche & Wurm 2012). The temperature gradient within dust particles is thought to form via direct irradiation from the central planet because the inner region of the circumjovian disk could be sufficiently depleted of dust particles. The inner region is therefore optically thin. The resulting dust-depleted inner region would have a higher ionization fraction, and would therefore admit increased accretion stress driven by magnetorotational instability (MRI; e.g., Balbus & Hawley 1991). This creates a bump in the gas surface density distribution and a pressure maximum at the boundary of the regions. We show that the photophoretic force in the circumjovian disk would overcome the gas drag, and the inward migration of dust particles would be stopped. We also reveal that this could form a pressure maximum and the bump in the surface density distribution of the circumjovian disk. In our model, the location of the bump is near the orbital radius of Io,and we do not need to consider the significant radius anomaly of young Jupiter. In addition, the temperatureof the disk at the current location of icy satellites could be lower than the sublimation temperature of H₂O ice, and the resonant-trapped icy satellites may survive in the disk. Consequently, photophoresis in the circumjovian disk has the potential to explain the formation of the Galilean satellites.

2 Models

2.1 Photophoresis and gas drag

In the circumplanetary disk with gas and dust, whether or not dust particles drift inward is determined by the balance of the residual gravity caused by the gas drag, F_D, and the photopho- retic force, F_Ph. In gaseous disks, other light-induced forces, such as the radiation pressure and the Poynting–Robertson effect, are negligible for μm- or mm-sized dust particles (Wurm & Krauss 2006, see also Appendix C).

For perfectly absorbing spherical particles whose radius is smaller than the mean free path of the gas molecules l_m, the photophoretic force F_Ph is given by (Beresnev et al. 1993; Wurm & Krauss 2006) $$F_{\text{Ph}}=\frac{0.5\pi a^{2}Ip}{3k_{\text{th}}T/a+12{\sigma }_{\text{SB}}{T}^4+p\sqrt{18k_{\text{B}}T/(\pi m_{\text{g}})}},$$(1)

where a is the radius of dust particles, I is the intensity of the radiative flux from the planet, p is the gas pressure, k_th is the thermal conductivity of dust particles, T is the temperature of the gas, σ_SB = 5.67 × 10⁻⁵ erg cm⁻² K⁻⁴ s⁻¹ is the Stefan–Boltzmann constant, k_B = 1.38 × 10⁻¹⁶ erg K⁻¹ is the Boltzmann constant, and m_g = 3.9 × 10⁻²⁴ g is the mean mass of gas molecules. We set a = 3 μm, which is consistent with the result of recent in-situ measurements of dust particles of the comet 67P/Churyumov-Gerasimenko (e.g., Bentley et al. 2016) and the estimation of the radius of constituent particles called monomers in the protoplanetary disk around the young protostar HL Tau (e.g., Okuzumi et al. 2016). The thermal conductivity of compact dust particles with no voids is about a few 10⁻¹   W m⁻¹ K⁻¹ for organic polymers (e.g., Choy 1977) and ~1 W m⁻¹ K⁻¹ for glass beads (e.g., Sakatani et al. 2017). Therefore we set k_th = 3 × 10⁴ erg s⁻¹ cm⁻¹ K⁻¹ (= 0.3 W m⁻¹ K⁻¹).

We checked the validity of the assumption that dust particles can be treated as perfectly absorbing particles for the radiation from the planet. The effective temperature of young Jupiter is approximately 1000 K (e.g., Fortney et al. 2011), and the characteristic wavelength of the radiation from the young Jupiter is λ ≃ 3 μm. Therefore, we can apply Eq. (1) for μm-sized dust particles whose size parameter x ≡ 2πa∕λ is larger than unity, that is, a ≳ 0.5 μm. We also confirmed the validity of the assumption of a < l_m (see Appendix D).

The rotation of dust particles is not considered in this study, although fast rotation results in a decrease in the photophoretic force. For μm-sized dust particles, the main source of random rotation is Brownian motion. We confirmed that the thermal rotation timescale is several orders of magnitude longer than the heat-conduction timescale and that the effect of rotation is negligible (see Appendix E).

The gas in a circumplanetary disk usually rotates slower than the Keplerian velocity owing to the negative pressure gradient in the accretion disk. In contrast, dust particles in the gas disk are not supported by the pressure gradient, and their stable orbital velocity is Keplerian; however, fine dust particles do rotate with the gas. Therefore, there is a gas drag on dust particles, and the residual gravity F_D is given by (Weidenschilling 1977; Wurm & Krauss 2006) $$F_{\text{D}}=-\frac{4\pi a^3}{3}\frac{{\rho }_{\text{p}}}{{\rho }_{\text{g}}}\frac{\text{d}p}{\text{d}r},$$(2)

where ρ_p = 2 g cm⁻³ is the density of dust particles, ρ_g is the gas density, and dp∕dr is the pressure gradient.

In the outer region of the circumplanetary disk, the residual gravity F_D is larger than the photophoretic force F_Ph, and the inward drift dominates, even in the presence of direct irradiation. However, there is a turnover distance where F_D = F_Ph, and the inward drift is prevented. Therefore, the inner region of the circumplanetary disk would be depleted of dust particles. Wurm & Krauss (2006) named this transition location the light barrier (see Fig. 1).

Fig. 1 Schematic of the circumplanetary disk with a dust-poor inner region formed by photophoresis. Both the pressure maximum and the bump in the surface density distribution are at the light barrier.

2.2 Circumplanetary disk model

We model the structure of the circumplanetary disk as Fujii et al. (2014) and Shibaike et al. (2017). In this study, we assume the circumplanetary disk is a steady-state viscous accretion disk with a supply from the protoplanetary disk. The gas surface density of the circumplanetary disk Σ is given by $$\Sigma =\frac{\dot{M}}{2\pi r_{\text{b}}}\frac{r^{3/2}}{\nu }\left(-\frac{2}{9}{r}^{-1/2}+\frac{2}{3}{r}_{\text{b}}{r}^{-3/2}\right),$$(3)

where Ṁ is the mass accretion rate, r_b is the radius of the region where the gas falls in, r is the distance from the planet, and ν is the kinematic viscosity. The numerical simulation by Tanigawa et al. (2012) found that r_b ≃ 20 R_J when the planet has a mass of 0.4M_J, where R_J = 7.15 × 10⁴ km and M_J = 1.9 × 10³⁰ g are the planetary radius and mass of Jupiter (see also Fujii et al. 2014), respectively. Assuming that r_b scales with the Hill radius of the planet, Shibaike et al. (2017) suggested that r_b = 27 R_J when the planet has a mass of 1M_J.

The gas temperature at the disk midplane T is given by (Nakamoto & Nakagawa 1994) $$T={\left[\frac{9\nu \Sigma {\Omega }_{\text{K}}{}^2}{8{\sigma }_{\text{SB}}}\left(\frac{3{\tau }_{\text{R}}}{8}+\frac{1}{2{\tau }_{\text{P}}}\right)\right]}^{1/4},$$(4)

where Ω_K is the Keplerian angular velocity, and τ_R and τ_P are the Rosseland and Planck mean optical depths, respectively. The Keplerian angular velocity is given by ${\Omega }_{\text{K}}=\sqrt{GM/r^3}$, where G = 6.67 × 10⁻⁸ dyn cm² g⁻² is the gravitational constant, and M is the mass of the central planet. We set M = M_J in this study. The sound velocity c_s is $c_{\text{s}}=\sqrt{k_{\text{B}}T/m_{\text{g}}}$, the gas scale height of the disk h_g is h_g = c_s∕Ω_K, the gas density at the disk midplane is ${\rho }_{\text{g}}=\Sigma /(\sqrt{2\pi }{h}_{\text{g}})$, and the pressure at the disk midplane is p = ρ_gk_BT∕m_g. We employ the standard α-prescription for the kinematic viscosity ν as $\nu =\alpha c_{\text{s}}{}^2/{\Omega }_{\text{K}}$ (Shakura & Sunyaev 1973).

The Rosseland mean optical depth τ_R is given by $${\tau }_{\text{R}}=\frac{\Sigma }{2}\left(f_{\text{d}}{\kappa }_{\text{R},\text{d}}+{\kappa }_{\text{R},\text{g}}\right),$$(5)

where f_d is the scaling factor of the dust abundance, and κ_{R, d} and κ_{R, g} are the contributions of dust and gas to the Rosseland mean opacity of a disk with solar metallicity. Similarly, the Planck mean optical depth τ_P is given by $${\tau }_{\text{P}}=\frac{\Sigma }{2}\left(f_{\text{d}}{\kappa }_{\text{P},\text{d}}+{\kappa }_{\text{P},\text{g}}\right),$$(6)

where κ_{P, d} and κ_{P, g} are the contributions of dust and gas to the Planck mean opacity of a disk with solar metallicity. The dust and gas opacities are taken from Henning & Stognienko (1996) and Freedman et al. (2014), respectively (see Appendix F).

Because the dust particles in the solar nebula tend to settle down toward the midplane, the dust-to-gas mass ratio of the inflow gas is thought to be lower than the dust-to-gas mass ratio of the solar nebula (e.g., Tanigawa et al. 2012). In addition, the dust-to-gas mass ratio in the inner region and the outer region must differ by orders of magnitude by the light barrier. Thus, we assumed that f_d = 10⁻¹ for the outer region, and f_d = 10⁻⁶ for the inner region.

In this study, we assume that the abundance of nm- or sub-μm-sized dust particles is negligibly small, and κ_{R, d} and κ_{P, d} are determined by the abundance of μm-sized dust particles which affect the photophoretic force. Thus, we expect a drastic change of f_d at the light barrier. We also neglect the turbulent diffusion of dust particles and consider the boundary at the light barrier to be sharp for simplicity. We note that the diffusion of dust particles may be important if the diffused region becomes optically thick and the intensity of the radiative flux at the light barrier decreases. Hence, future studies on this point are essential, although a detailed discussion is beyond the scope of this study.

The depletion of dust also affects the turbulent viscosity of the disk. Fujii et al. (2011, 2014) found that the midplane of circumplanetary disks is in general MRI-inactive because the diffusion timescale of the magnetic field is short. In the MRI-inactive region, it is likely that hydrodynamic instabilities, for example, the convective overstability (e.g., Klahr & Hubbard 2014) and subcritical baroclinic instability (e.g., Lyra & Klahr 2011), could dominate the angular momentum transport with values of α ~ 10⁻³ (e.g., Lyra 2014)¹. In contrast, the disk might be (partly) MRI-active if it is significantly depleted of dust grains (e.g., Turner et al. 2014; Fujii et al. 2014), because the lower dust abundance leads to a higher ionization fraction of the disk gas. This is true when the dominant fraction of the total surface area of dust particles (where free charges can stick) is in μm-sized dust particles which will be stopped at the light barrier. In this case, the effective turbulent viscosity could be larger than the viscosity of the completely MRI-inactive disks. Thus, we assumed α = 10⁻³ for the outer region and α = 10⁻² for the inner region. We note that the choice of absolute value of α does not change the results drastically; the increase in α at the light barrier is the key to generating the bump in the surface density distribution (see Appendix G).

3 Results

3.1 Location of the light barrier

First, we calculate the location of the light barrier in the steady-state accretion disk. At the location of the light barrier, the intensity of the radiative flux from the planet I should be I = L∕(4πr²), where L is the luminosity of the planet. As the rapid gas accretion was turned off and the planetary mass reached approximately 1 M_J, the luminosity of young Jupiter rapidly collapsed to 10⁻⁵L_⊙ (e.g., Marley et al. 2007), where L_⊙ = 3.83 × 10³³ erg s⁻¹ is the solar luminosity. Therefore, the luminosity of young Jupiter was on the order of 10⁻⁵ L_⊙ to 10⁻⁶ L_⊙ within the first 10 Myr (e.g., Fortney et al. 2011), and we set L = 4 × 10⁻⁶L_⊙ in our calculation.

We calculated the force ratio of the photophoretic force to the gas drag, F_Ph ∕F_D, by using both the normal model (f_d = 10⁻¹ and α = 10⁻³) and the dust-poor model (f_d = 10⁻⁶ and α = 10⁻²). The force ratio calculated from the normal model is ${\left(F_{\text{Ph}}/F_{\text{D}}\right)}_{\text{n}}$, and the force ratio calculated from the dust-poor model is ${\left(F_{\text{Ph}}/F_{\text{D}}\right)}_{\text{p}}$. We then defined the location of the light barrier r_lb as the location where $\min \left\{{\left(F_{\text{Ph}}/F_{\text{D}}\right)}_{\text{n}},{\left(F_{\text{Ph}}/F_{\text{D}}\right)}_{\text{p}}\right\}=1$, and we applied the normal model for the outer region whose distance is larger than r_lb and applied the dust-poor model for the inner region whose distance is smaller than r_lb.

Figure 2 shows the force ratio F_Ph∕F_D under the assumption that the intensity of the radiative flux is I = L∕(4πr²). This estimate of I is true at the light barrier. We found that r_lb = 6.9 R_J, 6.1 R_J, and 4.2 R_J for the cases of Ṁ = 10⁻⁷M_J yr⁻¹, 10⁻⁸M_J yr⁻¹, and 10⁻⁹M_J yr⁻¹, respectively.Therefore, the photophoretic force stops the inward migration of dust particles, and the inner region of the circumjovian disk is depleted of fine dust particles. Moreover, the location of this light barrier is close to the current orbit of the innermost satellite, Io (r = 5.9 R_J).

We note that for the case of Ṁ = 10⁻⁹M_J yr⁻¹, the force ratio is F_Ph∕F_D = 1 not only at the light barrier (r = 4.2 R_J) but also at r = 10.9 R_J. We speculate that the true location of the light barrier is between 4.2 R_J and 10.9 R_J, although a detailed calculation of the dust and gas disk structure around the light barrier is needed to determine its location.

Behind the light barrier, the intensity of the radiative flux would decrease exponentially due to the increase in the optical depth from the light barrier. Therefore, the photophoretic force might only work around the light barrier, and the structure of the outer region might not be affected. However, the structure of the inner region must be affected by photophoresis, and our results do not change significantly. We also note that r_lb does not significantly depend on the radius of dust particles when the two conditions, x ≫ 1 and a ≪l_m, are satisfied (see Appendix H).

Fig. 2 Force ratio FPh∕FD under the assumption that I = L∕(4πr2). The solid lines correspond to the force ratio calculated from the normal model (fd = 10−1 and α = 10−3), and the dashed lines are the dust-poor model (fd = 10−6 and α = 10−2). The vertical lines represent the location of the light barrier. The line colors represent the mass accretion rate (Ṁ = 10−7MJ yr−1 for blue, Ṁ = 10−8MJ yr−1 for pink, andṀ = 10−9MJ yr−1 for green). As a reference, the current orbit of Io is r = 5.9 RJ.

Fig. 3 Pressure distribution of the disk. The solid lines correspond to the pressure p calculated from the normal model, and the dashed lines are the dust-poor model. The vertical lines represent the location of the light barrier, and a pressure maximum is formed at the light barrier.

3.2 Pressure and surface density distributions

Figure 3 shows the pressure distribution of the circumjovian disk. We found that a pressure maximum was formed at the light barrier. This structure originates from the change in the value of α at the light barrier in the steady-state accretion disk. We speculate that the radial drift of pebbles and satellitesimals might stop at the pressuremaximum (e.g., Haghighipour & Boss 2003) and help the inner region to maintain its dust-poor environment.

Figure 4 shows the surface density distribution. We found that the bump in the surface density distribution is formed at the light barrier when the mass accretion rate is lower than 10⁻⁷M_J yr⁻¹. The migration direction and speed depend on the mass of the satellite and the disk structure (e.g., Tanaka et al. 2002; Paardekooper et al. 2011; Romanova et al. 2019). Fujii et al. (2017) found that the innermost satellite is stopped by the convergent migration to the bump. From our calculation, a steep bump in the surface density distribution is formed at the light barrier. Thus, the innermost satellite could be trapped around the light barrier.

Moreover, we confirmed the possibility of the formation of the resonant-trapped Galilean satellites by the following migration of Europa and Ganymede (see Appendix I). We found that the migration timescale is sufficiently longer than the critical timescale for resonance capture, and the adjacent satellites (i.e., Io–Europa and Europa–Ganymede) can be captured in the 2:1 mean motion resonance in the circumjovian accretion disk.

Fig. 4 Surface density distribution of the disk. The solid lines correspond to the surface density Σ calculated from the normal model, and the dashed lines are the dust-poor model. The vertical lines represent the location of the light barrier, and a bump in the surface density distribution is formed at the light barrier.

Fig. 5 Temperature distribution of the disk. The solid lines correspond to the temperature T calculated from the normal model, and the dashed lines are the dust-poor model.

3.3 Disk temperature

Figure 5 shows the temperature distribution. We found that the temperature at the light barrier is lower than the sublimation temperature of silicate dust (≳1500 K; Lodders 2003) when Ṁ ≲ 10⁻⁷M_J yr⁻¹. Therefore, photophoresis works before the sublimation of fine dust particles. Moreover, the temperature at the current location of Ganymede (r = 15 R_J) becomes lower than the sublimation temperature of H₂O ice (≃200 K) when Ṁ ≲ 10⁻⁸M_J yr⁻¹. This result indicates that the Galilean satellites formed in the circumjovian disk with a mass accretion rate of ≲ 10⁻⁸M_J yr⁻¹.

We note, however, that the disk temperature depends on the assumed opacity (e.g., Canup & Ward 2002), and the opacity should be calculated from the size distribution and the true abundance of the dust particles. The disk temperature is also dependent on the vertical structure of the accretion disk (e.g., Mori et al. 2019), and the vertical distribution of dust particles might be altered by photophoresis itself (e.g., McNally & Hubbard 2015; McNally & McClure 2017). Thus, future studies in these areas are essential.

4 Conclusions

The inner three Galilean satellites are in a 4:2:1 mean motion resonance, and their orbital configuration may be originated from the stopping of the migration of Io around the bump in surface density distribution, although the formation mechanism of the bump is under debate. We found that photophoresis in the circumjovian disk stopped the inward migration offine dust particles. The resulting dust-depleted inner region would have a higher ionization degree, and thus admit increased turbulent viscosity. This also forms a bump in the surface density distribution, and the migration of the innermost large satellite terminated around the bump. The location of the bump is near the current orbit of Io. Therefore, photophoresis is a strong candidate for the origin of the resonant-trapped Galilean satellites.

Acknowledgements

We thank Taishi Nakamoto for meaningful comments. S.A. is supported by the Grant-in-Aid for JSPS Research Fellow (JP17J06861). Y.S. is supported by a JSPS KAKENHI Grant (JP15H02065). Parts of this work have 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.

Appendix A Evaporation of icy satellites

If the ambient temperature is much higher than the sublimation temperature of H₂ O ice, then the icysatellites might be evaporated within the dissipation timescale of the circumplanetary disk. Here, we estimate the evaporation timescale. The mass and radius of the largest satellite, Ganymede, are M_Ga = 1.48 × 10²⁶ g and R_Ga = 2.63 × 10⁸ cm, and the evaporation timescale τ_evap is evaluated as $$\begin{array}{ll}{\tau }_{\text{evap}}\hfill & \simeq \frac{L_{\text{evap}}{M}_{\text{sat}}}{4\pi R_{\text{sat}}{}^2{\sigma }_{\text{SB}}{T}^4}\hfill \\ \hfill & \simeq 4.1\times {10}^4\left(\frac{M_{\text{sat}}}{M_{\text{Ga}}}\right){\left(\frac{R_{\text{sat}}}{R_{\text{Ga}}}\right)}^{-2}{\left(\frac{T}{500\text{\hspace{0.17em}K}}\right)}^{-4}\text{\hspace{0.17em}yr},\hfill \end{array}$$(A.1)

where L_evap = 2.7 × 10¹⁰ erg g⁻¹ is the latent heat of sublimation of H₂O ice (Tanaka et al. 2013), M_sat is the mass of the satellite, and R_sat is the radius of the satellite. Therefore, the evaporation timescale is shorter than the dissipation timescale of the circumplanetary disk (a few Myr). Moreover, the evaporation timescale is proportional to the radius R_sat and inversely proportional to the fourth power of the ambient temperature T, such that not only large satellites but also icy pebbles and satellitesimals cannot survive in the high-temperature environment.

Appendix B Inner cavity and corotation radius

If the inner cavity is located around the orbit of Io (r = 5.9 R_J), the orbital configuration of the Galilean satellites can be reproduced (Ogihara & Ida 2012; Shibaike et al. 2019). If the planetary magnetic field can be coupled effectively to the inner region of a circumplanetary accretion disk, then the disk could be disrupted at the radius r_in (e.g., Königl 1991): $$\begin{array}{ll}{r}_{\text{in}}\hfill & ={\left(\frac{{\mu }^4}{2GM{\dot{M}}^2}\right)}^{1/7}\hfill \\ \hfill & \simeq 5.2R_{\text{J}}{\left(\frac{B}{150\text{\hspace{0.17em}G}}\right)}^{4/7}{\left(\frac{\dot{M}}{{10}^{-8}{M}_{\text{J}}{\text{yr}}^{-1}}\right)}^{-2/7},\hfill \end{array}$$(B.1)

where μ ~ BR³ is the dipole moment of the planet, B is the surface magnetic field of the planet, and R = R_J is the planetary radius. Although the current magnetic field of Jupiter is 7.8 G, according to the observations from the Juno spacecraft (Bolton et al. 2017; Connerney et al. 2017), young Jupiter might have a stronger magnetic field, on the order of 100 G (e.g., Sánchez-Lavega 2004; Reiners & Christensen 2010). Thus, we can expect that r_in ≃ 5.9 R_J could be achieved if the ionization degree of the circumjovian disk is high enough, and the planetary magnetic field is coupled to the disk.

However, there is another requirement to make an inner cavity in a steady-state disk. To maintain its steady-state accretion, the corotation radius r_co should not be smaller than the radius of the inner cavity r_in (e.g., Königl 1991). The corotation radius is given by $$r_{\text{co}}={\left(\frac{GM}{{\omega }_{\text{spin}}{}^2}\right)}^{1/3}\simeq 5.3R_{\text{J}}{\left(\frac{2\pi /{\omega }_{\text{spin}}}{36\text{\hspace{0.17em}h}}\right)}^{2/3},$$(B.2)

where ω_spin is the angular velocity of the planetary spin. The rotation period of current Jupiter is 10 h; thus, the angular velocity of young Jupiter should have been about four times lower than that of the current Jupiter if the inner cavity were created by the magnetic coupling between young Jupiter and the circumjovian disk. In contrast, if r_co > r_in, a steady-state accretion disk cannot form, but a disk with no accretion might form. Instead, the spin rate of the central planet would decrease by propeller-driven outflows (e.g., Romanova et al. 2005). Then, the spin-down of the planet increases the corotation radius, reaching a stable state with r_co ≃ r_in, as Sasaki et al. (2010) predicted.

If we model the interior structure of the young Jupiter as a polytropic body with a polytropic index of 1.5, the spin angular momentum L_spin is given by (Batygin 2018) $$L_{\text{spin}}\simeq 0.21M{R}^2{\omega }_{\text{spin}}.$$(B.3)

Consequentially, the planetary radius of young Jupiter would be as large as twice that of the current Jupiter if we assumed conservationof spin angular momentum. Evolutionary calculations of the planetary interior demonstrated, however, that the planetary radius of young Jupiter should have been about 1.5 times the current radius at 1 Myr after the formation of Jupiter (e.g., Fortney et al. 2011). Therefore, we need to consider the significant radius anomaly of young Jupiter in this scenario.

We note that Sasaki et al. (2010) argued that Io could form at its current position, even if the inner cavity was at the current corotation radius (r = 2.25 R_J); however, subsequent N-body simulations by Ogihara & Ida (2012) revealed that the positions of the cavity and the innermost satellite must have been the same.

Appendix C Radiation pressure

The radiation pressure is caused by the direct momentum transfer of photons. The radiation force F_rad is given by $$F_{\text{rad}}=\frac{I}{c}\pi a^2,$$(C.1)

where c = 3.00 × 10¹⁰ cm s⁻¹ is the speed of light. Here, we compare the photophoretic force F_Ph and the radiation force F_rad. The force ratio of the radiation force to the photophoretic force is $$\begin{array}{ll}\frac{F_{\text{rad}}}{F_{\text{Ph}}}\hfill & \simeq \frac{6k_{\text{th}}T}{\text{cap}}\hfill \\ \hfill & \simeq 6\times {10}^{-2}\left(\frac{T}{300\text{\hspace{0.17em}K}}\right){\left(\frac{a}{3\mu \text{m}}\right)}^{-1}{\left(\frac{p}{{10}^2\text{\hspace{0.17em}dyn\hspace{0.17em}}{\text{cm}}^{-2}}\right)}^{-1}.\hfill \end{array}$$(C.2)

We then found that the radiation force is negligibly weaker than the photophoretic force in the gaseous circumplanetary disks (see also Loesche & Wurm 2012).

Appendix D Mean free path of gas molecules

The mean free path of the gas molecules l_m is given by $$l_{\text{m}}=\frac{m_{\text{g}}}{{\sigma }_{\text{mol}}{\rho }_{\text{g}}},$$(D.1)

where σ_mol = 2 × 10⁻¹⁵ cm² is the collisional cross section of gas molecules (Adachi et al. 1976). In the circumjovian disk, l_m is larger than 10⁻² cm at the location of the light barrier (see Fig. D.1), and the photophoretic force for small dust particles of a = 3 μm is given by Eq. (1), as we assumed.

For spherical particles whose radius is smaller than the mean free path of the gas molecules, a < l_m, the photophoretic force F_Ph is given by (Beresnev et al. 1993; Wurm & Krauss 2006) $$F_{\text{Ph}}\simeq \frac{\pi a^{3}Ip}{6k_{\text{th}}T},$$(D.2)

and therefore F_Ph does not depend on l_m. On the other hand, when the dust radius is larger than the mean free path of the gas molecules, a > l_m, the photophoretic force is approximately given by (Takeuchi & Krauss 2008) $$F_{\text{Ph}}\simeq \frac{\pi a^{3}Ip}{6k_{\text{th}}T}\times 4.5{\left(\frac{l_{\text{m}}}{a}\right)}^2,$$(D.3)

and F_Ph decreases with increasing of a. In contrast to F_Ph, the residual gravity caused by gas drag F_D is independent of a∕l_m (Weidenschilling 1977): $$F_{\text{D}}=-\frac{4\pi a^3}{3}\frac{{\rho }_{\text{p}}}{{\rho }_{\text{g}}}\frac{\text{d}p}{\text{d}r}.$$(D.4)

Fig. D.1 Mean free path of gas molecules in the disk. The solid lines correspond to the mean free path lm calculated from the normal model, and the dashed lines are the dust-poor model.

Therefore, F_Ph∕F_D is independent of a for the case of a < l_m, while F_Ph∕F_D is inversely proportional to the square of the particle radius for the case of a > l_m.

This means that the radial drift of cm-sized compact pebbles could not be prevented by photophoresis itself, although the radial drift would be stopped because the pressure maximum formed at the light barrier. We note, however, that the thermal conductivityof dust aggregates is exceedingly lower than that of compact dust particles (e.g., Arakawa et al. 2017, 2019), and the radial drift of aggregated pebbles might be stopped by photophoresis, even if their radius is larger than the mean free path of the gas molecules.

Appendix E Rotationof dust particles

If particles rotate faster than the heat conduction timescale, the temperature gradient would be reduced, and the photophoretic force might be suppressed. For μm-sized dust particles, the main source of particle rotation is the Brownian motion. The thermal rotation timescale of a spherical particle, τ_rot, is (Krauss & Wurm 2005; Loesche & Wurm 2012) $${\tau }_{\text{rot}}=\sqrt{\frac{8{\pi }^2{a}^5{\rho }_{\text{p}}}{45k_{\text{B}}T}}\simeq 1.4\times {10}^{-2}{\left(\frac{a}{3\mu \text{m}}\right)}^{5/2}{\left(\frac{T}{300\text{\hspace{0.17em}K}}\right)}^{-1/2}\text{\hspace{0.17em}s}.$$(E.1)

The heat conduction timescale within a dust particle, τ_cond, is (Krauss & Wurm 2005; Loesche & Wurm 2012) $${\tau }_{\text{cond}}=\frac{a^2}{k_{\text{th}}}{\rho }_{\text{p}}{c}_{\text{p}}\simeq 6.0\times {10}^{-5}{\left(\frac{a}{3\mu \text{m}}\right)}^2\text{\hspace{0.17em}s},$$(E.2)

where c_p = 10⁷ erg g⁻¹ K⁻¹ is the specific heat (e.g., Ciesla et al. 2004). We found that τ_rot ≫ τ_cond, and the effect of particle rotation is negligible for μm-sized dust particles.

In addition to Brownian motion, dust-dust collisions spin up dust particles. However, this effect is negligibly small for a dust-to-gas mass ratio lower than unity (see McNally & McClure 2017).

Appendix F Dust and gas opacities

We calculated the opacity of the disk using standard temperature-dependent opacity models (Henning & Stognienko 1996; Freedman et al. 2014).

Fig. F.1 Planck mean opacity of gas κP, g calculated by Freedman et al. (2014) (gray solid lines) and the analytical fit used in this study (blue dashed line).

The contribution of dust to the Rosseland mean opacity of a disk with solar metallicity, κ_R,d, is given by Henning & Stognienko (1996), and we use the iron-poor model of the study. The contribution of dust to the Rosseland mean opacity of a disk with solar metallicity, κ_{P, d}, is given by κ_{P, d} = 2.4κ_{R, d} (Nakamoto & Nakagawa 1994).

The contribution of gas to the Rosseland mean opacity of a disk with solar metallicity, κ_R,g, is given by Freedman et al. (2014), and we use their analytical fit of κ_R,g for the solar metallicity gas. We also evaluated κ_{P, g}, using Table 3 of Freedman et al. (2014) as $${\kappa }_{\text{P},\text{g}}={\left[15{\left(T/100\text{\hspace{0.17em}K}\right)}^{-8}+0.15{\left(T/100\text{\hspace{0.17em}K}\right)}^{0.8}\right]}^{-1}{\text{cm}}^2{\text{g}}^{-1}.$$(F.1)

Within the temperature range of 100–2000 K, the pressure dependence of κ_{P, g} is small. Therefore, we assume κ_{P, g} is only the function of the temperature (see Fig. F.1).

We note that the dominant contributions to the Rosseland mean opacity of gas are in opacity windows, while the dominant contributions to the Planck mean opacity of gas are strong bands (see Fig. 2 of Freedman et al. 2014). Therefore, κ_{P, g} is several orders of magnitude larger than κ_{R, g} in general.

Appendix G Dependence on the turbulent viscosity

In the main text, we assumed α = 10⁻³ for the outerregion and α = 10⁻² for the inner region. However, the choice of absolute values of α does not change the results significantly,but the increase in α at the light barrier is the key to generate the bump in the surface density distribution.

Figure G.1 shows the force ratio F_Ph∕F_D for various viscosity parameters. Here, we assume that the value of α for the inner region is ten times higher than that for the outer region, and the mass accretion rate is Ṁ = 10⁻⁸M_J yr⁻¹. The line colors represent the value of α for the outerregion (α = 10⁻³ for pink, α = 5 × 10⁻⁴ for orange, and α =2.5 × 10⁻⁴ for blue). We found that the location of the light barrier is r_lb = 6.1 R_J, 7.5 R_J, and 9.0 R_J for the cases of α =10⁻³, 5 × 10⁻⁴, and 2.5 × 10⁻⁴, respectively. The location of the light barrier is not strongly dependent on α and we can adjust r_lb = 5.9 R_J by choosing appropriate values of the luminosity L and mass accretion rate Ṁ.

Fig. G.1 Force ratio FPh∕FD. The solid lines correspond to the force ratio calculated from the normal model, and the dashed lines are the dust-poor model. The vertical lines represent the location of the light barrier. The line colors represent the value of α for the outer region (α = 10−3 for pink, α = 5 × 10−4 for orange, and α = 2.5 × 10−4 for blue). Weassume that the value of α for the inner region is ten times higher than that for the outer region, and the mass accretion rate is Ṁ = 10−8MJ yr−1.

Figure G.2 shows the surface density distribution of the disk for various viscosity parameters. We confirmed that a bump in the surface density distribution forms at the light barrier in all cases. Therefore, the general view of the disk structure hardly depends on the choice of absolute values of α, provided that the turbulence of the inner region is stronger than that of the outer region.

Fig. G.2 Surface density distribution of the disk for various viscosity parameters. The solid lines correspond to the surface density Σ calculated from the normal model, and the dashed lines are the dust-poor model. The vertical lines represent the location of the light barrier, and a bump in the surface density distribution is formed at the light barrier.

Appendix H Dependence on the particle radius

For μm-sized dust particles in the circumjovian disk, we can simplify Eq. (1) as follows: $$F_{\text{Ph}}\simeq \frac{\pi a^{3}Ip}{6k_{\text{th}}T}.$$(H.1)

Thus the force ratio F_Ph∕F_D is given by $$\frac{F_{\text{Ph}}}{F_{\text{D}}}\simeq \frac{{\rho }_{\text{g}}Ip}{8{\rho }_{\text{p}}{k}_{\text{th}}T}{\left(\frac{\text{d}p}{\text{d}r}\right)}^{-1},$$(H.2)

and the location of the light barrier r_lb is independent of the particle radius a, because both F_Ph and F_D are proportional to the cubic of a.

Appendix I Condition for resonance capture

The Galilean satellites are believed to have migrated in the circumjovian disk by exchanging their angular momentum with the disk. If the migration timescale, τ_mig, is longer than the critical migration timescale for resonance capture, τ_crit, the satellites are captured in the 2:1 mean-motion resonance. Ogihara & Kobayashi (2013) performed N-body simulations and revealed that the critical migration timescale for resonance capture τ_crit is given by $${\tau }_{\text{crit}}\simeq 2\times {10}^5{\left(\frac{M_{\text{inner}}}{M_{\text{Io}}}\right)}^{-4/3}{\Omega }_{\text{K}}{}^{-1},$$(I.1)

where M_Io = 8.9 × 10²⁵ g is the mass of Io and M_inner is the mass of the inner satellite. The migration timescale τ_mig is given by (e.g., Paardekooper et al. 2011; Fujii et al. 2017) $$\begin{array}{ll}{\tau }_{\text{mig}}\hfill & \simeq \frac{1}{5}{\left(\frac{M_{\text{sat}}}{M}\right)}^{-1}{\left(\frac{r^2\Sigma }{M}\right)}^{-1}{\left(\frac{c_{\text{s}}}{r{\Omega }_{\text{K}}}\right)}^2{\Omega }_{\text{K}}{}^{-1}\hfill \\ \hfill & \simeq 7\times {10}^7\left(\frac{T}{300\text{\hspace{0.17em}K}}\right){\left(\frac{\Sigma }{300\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-2}}\right)}^{-1}{\left(\frac{r}{9R_{\text{J}}}\right)}^{-1}{\left(\frac{M_{\text{sat}}}{M_{\text{Eu}}}\right)}^{-1}{\Omega }_{\text{K}}{}^{-1},\hfill \end{array}$$(I.2)

where M_Eu = 4.8 × 10²⁵ g is the mass of Europa. We found that τ_mig ≫ τ_crit and, therefore, the adjacent satellites (i.e., Io–Europa and Europa–Ganymede) can be captured in the 2:1 mean motion resonance.

The migration timescale τ_mig is on the order of a few 10⁵ years for the Galilean satellites, and they can migrate in the circumjovian disk. We note, however, that τ_mig might be longer than the dissipation timescale of the disk when the mass of the satellite is far smaller than 10⁻²M_Eu. This fact indicates that Europa and Ganymede should be larger than 10⁻² times the current mass at the time of resonance capture.

References

All Figures

	Fig. 1 Schematic of the circumplanetary disk with a dust-poor inner region formed by photophoresis. Both the pressure maximum and the bump in the surface density distribution are at the light barrier.
In the text

Fig. 2 Force ratio FPh∕FD under the assumption that I = L∕(4πr2). The solid lines correspond to the force ratio calculated from the normal model (fd = 10−1 and α = 10−3), and the dashed lines are the dust-poor model (fd = 10−6 and α = 10−2). The vertical lines represent the location of the light barrier. The line colors represent the mass accretion rate (Ṁ = 10−7MJ yr−1 for blue, Ṁ = 10−8MJ yr−1 for pink, andṀ = 10−9MJ yr−1 for green). As a reference, the current orbit of Io is r = 5.9 RJ.

In the text

	Fig. 3 Pressure distribution of the disk. The solid lines correspond to the pressure p calculated from the normal model, and the dashed lines are the dust-poor model. The vertical lines represent the location of the light barrier, and a pressure maximum is formed at the light barrier.
In the text

	Fig. 4 Surface density distribution of the disk. The solid lines correspond to the surface density Σ calculated from the normal model, and the dashed lines are the dust-poor model. The vertical lines represent the location of the light barrier, and a bump in the surface density distribution is formed at the light barrier.
In the text

	Fig. 5 Temperature distribution of the disk. The solid lines correspond to the temperature T calculated from the normal model, and the dashed lines are the dust-poor model.
In the text

	Fig. D.1 Mean free path of gas molecules in the disk. The solid lines correspond to the mean free path lm calculated from the normal model, and the dashed lines are the dust-poor model.
In the text

	Fig. F.1 Planck mean opacity of gas κP, g calculated by Freedman et al. (2014) (gray solid lines) and the analytical fit used in this study (blue dashed line).
In the text

Fig. G.1 Force ratio FPh∕FD. The solid lines correspond to the force ratio calculated from the normal model, and the dashed lines are the dust-poor model. The vertical lines represent the location of the light barrier. The line colors represent the value of α for the outer region (α = 10−3 for pink, α = 5 × 10−4 for orange, and α = 2.5 × 10−4 for blue). Weassume that the value of α for the inner region is ten times higher than that for the outer region, and the mass accretion rate is Ṁ = 10−8MJ yr−1.

In the text

	Fig. G.2 Surface density distribution of the disk for various viscosity parameters. The solid lines correspond to the surface density Σ calculated from the normal model, and the dashed lines are the dust-poor model. The vertical lines represent the location of the light barrier, and a bump in the surface density distribution is formed at the light barrier.
In the text