The Astrophysical Journal, 846:81 (10pp), 2017 September 1

https://doi.org/10.3847/1538-4357/aa8454© 2017. The American Astronomical Society. All rights reserved.

Satellitesimal Formation via Collisional Dust Growth in Steady Circumplanetary Disks
1

Yuhito Shibaike1, Satoshi Okuzumi1, Takanori Sasaki2, and Shigeru Ida3

Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Meguro-ku, Tokyo, 152-8551, Japan2
Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto, 606-8502, Japan
3
Earth-Life Science Institute, Tokyo Institute of Technology, Meguro-ku, Tokyo 152-8550, Japan
Received 2017 May 16; revised 2017 July 31; accepted 2017 August 2; published 2017 September 1

Abstract
The icy satellites around Jupiter are considered to have formed in a circumplanetary disk. While previous models
have focused on the formation of the satellites starting from satellitesimals, the question of how satellitesimals
themselves form from smaller dust particles has not yet been addressed. In this work, we study the possibility thatsatellitesimals form in situ in a circumplanetary disk. We calculate the radial distribution of the surface density andrepresentative size of icy dust particles that grow by colliding with each other and drift toward the central planet ina steady circumplanetary disk with a continuous supply of gas and dust from the parent protoplanetary disk. The
radial drift barrier is overcome if the ratio of the dust-to-gas accretion rates onto the circumplanetary disk, M˙ d M˙ g,
is high and the strength of turbulence, α, is not too low. The collision velocity is lower than the critical velocity offragmentation when α is low. Taken together, we ﬁnd that the conditions for satellitesimal formation via dust
coagulation are given by M˙ d M˙ g  1 and 10-4  a < 10-2 . The former condition is generally difﬁcult to achieve,suggesting that the in situ satellitesimal formation via particle sticking is viable only under extreme conditions. We
also show that neither satellitesimal formation via the collisional growth of porous aggregates nor via streaming
instability is viable as long as M˙ d M˙ g is low.
Key words: accretion, accretion disks – planets and satellites: formation – turbulence
Galilean satellites. Sasaki et al. (2010) also reproduced theSaturnian system.
One signiﬁcant problem for these satellite formation studies isthat how satellitesimals form from dust particles are unaddressed.
Canup & Ward (2006), Sasaki et al. (2010), and Ogihara & Ida(2012) assumed that the dust particles grow to satellitesimals assoon as they are delivered from the protoplanetary to thecircumplanetary disks. However, it is already known fromplanetesimal formation studies that dust growth to kilometersized bodies in protoplanetary disks can be hindered byphenomena such as the radial drift and the collisional
fragmentation of intermediate-sized particles (Whipple 1972;
Adachi et al. 1976; Weidenschilling 1977). Therefore, it is easy
to imagine that satellitesimal formation in circumplanetary disks
could suffer from similar difﬁculties.
In this work, we aim to answer the question of whether dust
particles can grow to satellitesimals by their direct collisional
growth in circumplanetary disks. We employ a simple onedimensional model in which we calculate the radial distributionof the surface density and typical size of dust particles in a
steadily accreting circumplanetary disk. We also consider only
icy dust particles and do not consider rocky particles. Although
we assume perfect sticking upon collision, fragmentation
occurs if the collision velocity is higher than a few m s-1
when the aggregates are mainly composed of silicate particles(e.g., Blum & Wurm 2008; Wada et al. 2009). The majority of
the Galilean satellites are indeed icy satellites: Europa is ~10%
and Ganymede and Callisto are ~50% ice by mass (Sohl et al.
2002). Our simple treatment allows us to explore a large
parameter space. The goal of this work is to derive theconditions under which satellitesimal formation via direct
coagulation of dust particles is viable.
We note that it has also been discussed that planetesimals
can be captured by gas drag from the circumplanetary disks(Fujita et al. 2013; Tanigawa et al. 2014; D’Angelo &1. Introduction
The gas giants in the solar system have large satellites whose
orbits are on the same plane. These regular satellites are
believed to have formed in gas disks around the planets, calledcircumplanetary disks, like protoplanetary disks that form
planets around stars (Lunine & Stevenson 1982). Numericalsimulations show that circumplanetary disks are by-products ofgas accretion from the parent protoplanetary disk to the gas
giants (e.g., Lubow et al. 1999; D’Angelo et al. 2002;
Tanigawa et al. 2012; Fung & Chiang 2016; Szulágyi et al.
2016).
Satellite formation in circumplanetary disks has been studied
extensively. Lunine & Stevenson (1982) proposed the minimummass disk model where the circumplanetary disk is static and its
solid mass is equal to the total mass of the satellites. Mosqueira &Estrada (2003a, 2003b), Estrada (2009), and Mosqueira et al.
(2010) showed that the Jovian and Saturnian satellites can formin the minimum-mass disks with a low-density outer portion,although the satellites would have to overcome rapid inward
migration induced by the interaction with the massive disk
(Miguel & Ida 2016). Canup & Ward (2002) and Ward & Canup(2010) proposed the gas-starved disk whose mass is regulated by
viscous accretion onto the central planet and by a gradual supplyof the gas from the parent protoplanetary disk. They showed thatthe temperature of the disk is sustained low enough to produce
icy satellites. Alibert et al. (2005) also developed a similaraccretion disk model that is consistent with a formation model for
Jupiter. Canup & Ward (2006) performed N-body simulations ofsatellite formation in the gas-starved disk, showing that the total
mass of the satellites formed in the disk is ~10-4 of the central
planet’s mass, as observed for the Jovian, Saturnian, and Uraniansatellite systems. Sasaki et al. (2010) and Ogihara & Ida (2012)
also studied the satellite formation in the gas-starved disk modelusing Monte Carlo and N-body simulations, respectively, andsuccessfully reproduced the number, masses, and orbits of the1

The Astrophysical Journal, 846:81 (10pp), 2017 September 1

Shibaike et al.

(23) and 25 of Fujii et al. 2014)M˙ g r 3 2 ⎛ 2 -1 2
⎞
2
⎜- r
Sg =
+ r b r -3 2 ⎟ ,⎝⎠2pr b n93Podolak 2015; Suetsugu et al. 2016; Suetsugu & Ohtsuki 2017).
Suetsugu et al. (2016) and Suetsugu & Ohtsuki (2017)
examined the captures and subsequent orbital evolutions ofplanetesimals. They showed that the capture hypothesis could
roughly reproduce the initial radial distribution of planetesimals(i.e., satellitesimals) assumed in the satellite formation modelby Canup & Ward (2006), Sasaki et al. (2010), and Ogihara &Ida (2012). However, if there is a gas gap around thecircumplanetary disk, the positive pressure gradient outsideof the gap removes low eccentricity planetesimals from thefeeding zone of the planet and the capture rate decreases(Kobayashi et al. 2012; Fujita et al. 2013; Tanigawa et al. 2014;
Suetsugu & Ohtsuki 2017).
The structure of this paper is as follows. In Section 2, we
describe our model for the circumplanetary gas disk andexplain how we treat the collisional growth and radial drift ofthe dust particles in the disk. In Section 3, we show the resultsof our calculations and derive the conditions for successful
satellitesimal formation via dust coagulation. In Section 4, we
discuss the feasibility of the conditions and consider otherprobabilities of satellitesimal formation. Finally, we conclude
this paper in Section 5.

where Ṁg is the mass accretion rate of the infall gas. The
simulation by Tanigawa et al. (2012) shows that r b » 20RJ for
the planet of Mcp = 0.4MJ (see also Fujii et al. 2014).
Assuming that r b scales with the Hill radius of the central
planet, we use r b = 27RJ for our 1MJ -mass planet.
We assume that the circumplanetary disk is viscously heated.
Then, the gas temperature at the midplane is given by
Nakamoto & Nakagawa (1994),
⎛ 9
⎞1 4
T=⎜
n Sg W2K ⎟ g (t ) ,⎝ 8sSB⎠⎛3
1 ⎞⎟1 4g (t ) = ⎜ t +⎝8
4.8t ⎠

We model the structure of the circumplanetary disk
following Fujii et al. (2014). Although some numericalsimulations suggested the possibility that gas near the midplane
spirals outward (Tanigawa et al. 2012; Fung & Chiang 2016;
Szulágyi et al. 2016), here we assume that the circumplanetary
disk is a viscous accretion disk with a continuous supply ofmaterial from the protoplanetary disk. The diffusion equation
for the gas surface density Sg of the circumplanetary disk isthen given by¶t⎤
1 ¶ ⎡ 1 2¶ 1 2(r n Sg) ⎥ + f ,
⎢⎣3r
⎦
¶r
r ¶r

(4 )

is a function of the Rosseland mean optical depth τ. In
principle, τ depends on the size distribution of the smallest dust
particles, which cannot be predicted with simple dust evolution
models as employed in this study. Lacking good knowledge of
τ, we opt to set g » 1. Since g  1 in general, the assumptiong » 1 yields a minimum estimate for the disk temperature. The
temperature can be up to three times higher than what isassumed here if the optical depth ranges between
10-2  t  10 2 . However, this uncertainty has little effecton the main results of this work because the dependence of ourresults on T is weak (see Section 3.2).
According to the simulation by Tanigawa et al. (2012), themass ﬂux Ṁg scales as Ṁg » 0.2SPPD rH2 WPPD , where SPPD and
WPPD are the gas surface density and orbital period of the parent
protoplanetary disk in the vicinity of the planet, respectively,
and rH is the planet’s Hill radius (see Figure 14 of Tanigawa
et al. 2012). At Jupiter’s orbit, the gas surface density is
143 gcm-2 according to the minimum-mass solar nebula modelof Hayashi (1981). For this value of SPPD , the accretion rate
onto a Jupiter-sized planet would be Ṁg » 200 MJ Myr-1.
However, a strong constraint on Ṁg can be obtained from thetemperature of the circumplanetary disk. In the left panel of
Figure 1, we plot the midplane temperature of our modeled
circumplanetary disk as a function of the distance from thecentral planets for three cases Ṁg = 2, 0.2, 0, 02, and
0.002 MJ Myr-1. For Ṁg = 2 and 0.2 MJ Myr-1, the temperatures are higher than the sublimation temperature of ice, which
is about 160 K. Such hot environments are unsuitable for theformation of icy regular satellites around Jupiter and Saturn.
Therefore, we only consider Ṁg  0.02 MJ Myr-1 in this work.
We plot the gas surface density of the disk for the two cases in
the right panel of Figure 1.

2.1. Circumplanetary Disk Model

=

(3 )

where sSB is the Stefan–Boltzmann constant and2. Methods

¶Sg

(2 )

(1 )

where r is the distance from the central planet, f is the mass ﬂux
of the gas inﬂow from the protoplanetary to the circumplanetary disks, and ν is the kinematic viscosity. We employ thestandard α prescription (Shakura & Sunyaev 1973) and express
the viscosity as n = acs Hg , where cs is the isothermal soundspeed and Hg is the gas scale height. The sound speed is related
to the temperature as cs = kB T mg , with kB as the Boltzmann
constant and mg = 3.9 ´ 10-24 g as the mean molecular mass.
The gas scale height is given by Hg = cs WK , where
WK = GMcp r 3 is the Kepler frequency, and G and Mcp are
the gravitational constant and the central planet mass,
respectively. Unless otherwise noted, we assume Mcp to be
the Jupiter mass, MJ = 1.89 ´ 1030 g. Based on the results ofthe three-dimensional hydrodynamical simulation by Tanigawa
et al. (2012), Fujii et al. (2014) modeled f as f µ r -1 for r < r b
and f=0 for r > r b, where r b is the radius of the region where
the gas falls in. With this scaling for f, the steady-state solutionof Equation (1) can be analytically obtained as (see Equations

2.2. Dust Growth and Radial Drift
We treat the coagulation and radial drift of dust particles incircumplanetary disks. In particular, we focus on the steadystate where the radial distribution of the size and surface
density of the particles is constant over time. We approximate
the size distribution of the particles at each distance from thecentral planet by a narrow distribution peaked at mass md . With
2

The Astrophysical Journal, 846:81 (10pp), 2017 September 1

Shibaike et al.

Figure 1. Gas surface density and the midplane temperature of the circumplanetary disk. The red, blue, green, and purple curves in the left panel represent thetemperatures, where Ṁg = 2, 0.2, 0.02 , and 0.002 MJ Myr-1, respectively. The horizontal black line is the sublimation temperature of icy aggregates. The vertical
black dashed lines represent the current orbits of the Galilean satellites. The blue and purple curves in the left panel represent the gas surface density, where Ṁg = 0.02
and 0.002 MJ Myr-1, respectively. The strength of the turbulence is a = 10-5, 10-4, 10-3, and 10−2.

in terms of the Stokes number deﬁned by St = WK tstop . In
dense circumplanetary disks, one can safely assume thatparticles of Rd > 0.1 mm are much larger than the mean freepath of gas molecules, so the ﬂow around the particles can be
regarded as a continuous ﬂuid. Then, the Stokes number can be
expressed as

this approximation, the integro-differential equation governing
the evolution of particles in a disk can be rewritten into a
simple differential equation for md . We also approximate theradial distribution of dust inﬂow to the circumplanetary disk bya narrow peak lying at r = r b , the outer edge of the infall
region. With this approximation, the problem of obtaining theradial distribution of Sd and md in a steady-state reduces to asimple boundary-value problem. We discuss the validity ofthese assumptions in Section 4.4.
In addition to the above approximations, for the moment weassume that the collision velocities between the particles is so
low that their fragmentation is negligible. The issue of thefragmentation barrier will be separately treated in Section 3.3.
We also assume that the initial particle radius (i.e., the particle
radius at r = r b ) is 0.1 mm . This is the maximum size of theparticles that can diffuse into the gas gap against the outward
drift motion caused by the positive gas pressure gradient at the
gap’s outer edge (Zhu et al. 2012). We focus on the steady-state
where the radial distribution of the mass and surface density ofthe particles is independent of time.
In a steady-state, and in the absence of fragmentation, themass md of the radially drifting particles is determined as afunction of r by (Equation (5) of Sato et al. 2016)
vr

2
dm d
=
dr

p R d2 Dvdd
Hd

Sd ,

St =

(7 )

where Dvdg is the relative velocity between the dust particles
and the gas, and CD is a dimensionless coefﬁcient that depends
on the particle Reynolds number, Rep. According to Perets &Murray-Clay (2011), the coefﬁcient can be written asCD =24(1 + 0.27 Re p )0.43 + 0.47
Re p
´ (1 - exp ( - 0.04 Re0.38
p )).

(8 )

The particle Reynolds number is given by
Re p =(5 )4R d Dvdg
vth lmfp,(9 )where vth = 8 p cs is the thermal velocity andlmfp = mg (smol rg ) is the mean free pass of the gas withsmol = 2 ´ 10-15 cm2 the collisional cross section of the gas
molecules and rg = Sg ( 2p Hg ) is the gas density at the
midplane.
The scale height of the particles can be derived analyticallyfrom the balance of their vertical sedimentation and diffusion(Youdin & Lithwick 2007),

where Rd , Dvdd , vr , Hd , and Sd are the radius, collisionvelocity, drift speed, scale height, and surface mass densityof the particles, respectively. The particle mass is related to theparticle radius by md = (4p 3) Rd3 rint , where rint is the internal
density of the particles. We ﬁx rint = 1.4 g cm-3 but we
discuss the effects of changing it in Section 4.2. The inward
accretion rate of the particles isM˙d = - 2prvr Sd ,8 r int R d
WK ,3CD rg Dvdg⎛
St 1 + 2St ⎞-1 2⎟
Hd = Hg ⎜1 +.
⎝a 1 + St ⎠(6 )which is, in steady-state, constant over r, and is equal to theinfall rate set at r = rd as a boundary condition. We
numerically integrate Equation (5) with Equation (6) fromr = r b toward smaller r.
The radial drift velocity of a dust particle is determined by its
stopping time, tstop . In this study, we express the stopping time(10)

The radial drift velocity of the dust particles is (Whipple 1972;
Adachi et al. 1976; Weidenschilling 1977)
vr = - 2

3

St
St2 + 1

hv k ,(11)

The Astrophysical Journal, 846:81 (10pp), 2017 September 1

Shibaike et al.

the gap. This means that realistic values of Ṁg should be lessthan 2 MJ Myr-1, i.e., less than 1% of the accretion rate without
a gap. The infall rate that was adopted in the gas-starved disk
model of Canup & Ward (2002) is about 0.2 MJ Myr-1. The
value of Ṁg whose temperature is suitable for icy satellite
formation is lower than these estimated values. Although we do
not consider the decrease of the gas inﬂow in detail, the ﬁnal
phase of planetary formation should be suitable for satellite
formation (see also Section 4.1).
The ratio M˙ d M˙ g of the dust inﬂow rate to the gas inﬂow rateis chosen between 0.001–1. If we assume that the inﬂow has asolar composition and the dust particles are strongly coupledwith the gas, the ratio should be 0.01. However, the gap
structure of the gas around Jupiter dams the dust particles thatdrifted from the outer region of the protoplanetary disk toward
the Sun, and generates their similar gap structure. On the other
hand, the strong gas gradient may trigger a hydrodynamic
instability and disturb the gas inﬂow. This disturbance should
enhance the radial diffusion of small dust particles and makethem nearer to the central planet (Zhu et al. 2012; Tanigawa
et al. 2014).
The strength α of turbulence is varied from 10−5 to 10−2.
We note that the mechanisms that could drive turbulence in thecircumplanetary disk are highly uncertain. The magnetorotational instability, a viable mechanism driving turbulence inionized accretion disks, could operate on the surface of
circumplanetary disks (Turner et al. 2014), but might not
produce fully developed turbulence (Fujii et al. 2014).
Therefore, we cannot rule out that α of circumplanetary disks
falls below α. However, as we discuss in the following section,
turbulence that is too weak could make it difﬁcult forsatellitesimals to form within a realistic range of M˙ d M˙ g.

Table 1
Parameter Choice
Quantity

Description

Value

Ṁg
M˙ d M˙ g

Gas infall rate
Dust-to-gas infall rate ratio
Turbulence parameter

0.02, 0.002 MJ Myr-1
1, 0.1, 0.01, 0.001
10-5, 10-4, 10-3, 10-2αwhere v k = r Wk is the Kepler velocity and2
1 ⎛ Hg ⎞2 ¶ ln rg cs
h=- ⎜ ⎟
2 ⎝ r ⎠ ¶ ln r

(12)

is the ratio of the pressure gradient force to the gravity of thecentral planet.
The relative velocity between the dust particles (i.e.,collision velocity) is the root sum square
Dvdd =

DvB2 + Dvr2 + Dvf2 + Dvz2 + Dvt2 ,(13)

where Dv B, Dvr , Dvf, Dvz , and Dvt are the relative velocities
induced by Brownian motion, the radial drift, azimuthal drift,
vertical sedimentation, and turbulence (Okuzumi et al. 2012).
For collisions between equal-sized particles, the Brownianmotion-induced velocity can be written as Dv B =16kB T (pmd ) . The relative velocity induced by the radial
drift is Dvr = ∣vr (St1) - vr (St2)∣, where St1 and St2 are the
Stokes numbers of the two particles. We assume St2 = 0.5St1(see Section 2.4 in Sato et al. 2016) and vr is given by
Equation (11). The relative velocity induced by the azimuthal
drift is Dvf = ∣vf (St1) - vf (St2)∣, where vf = -hv K (1 + St2).
The relative velocity induced by the vertical motion is
Dvz = ∣vz (St1) - vz (St2)∣, where vz = -WK StHd (1 + St).
For the relative velocity induced by the turbulence, we use
the analytic formula derived from Ormel & Cuzzi (2007). The
formula has three limiting expressions:⎧ a c Re1 4∣St - St ∣ , St  Re-1 2,
s
1
2
1
t
t
⎪
1
2
1
2
⎪ 3a cs St1 , Re t
 St1  1,
Dvt = ⎨
1 2
⎛ 1
⎪
1 ⎞
+
⎟ , 1  St1.
⎪ a cs ⎜⎝ 1 + St1
1 + St 2 ⎠⎩3. Results
In this section, we present the results of our dust growth
calculations and explore the conditions under which satellitesimals can form through dust coagulation in circumplanetary disks.
3.1. Fiducial Calculations
Figure 2 shows the results of our ﬁducial calculations thatassume Ṁg = 0.02 MJ Myr-1 and a = 10-4 . We assume thatthe snowline is where the midplane temperature is 160 K and it
is at r = 10 RJ (see Figure 1). The top, middle, and bottom
panels represent the surface density, radius, and Stokes numberof the mass-dominating dust particles in the disk as a functionof the distance r from the central planet. Because the particles
grow and move inward at the same time, the particle size
increases with decreasing r. As already mentioned in theprevious section, we assume that the dust particles are 0.1 mm
in size when they are initially delivered from the protoplanetary
disk to r = r b (=27RJ ). However, the middle panel of 2
suggests that the assumption about the initial size of the dust
particles is not crucial because the particles immediately growat r » r b . The change of r b will not affect the steady-state
proﬁles either because the proﬁles should gradually approach
the approximated lines (see the dashed lines in the bottom
panel of Figure 2 and Equation (15)).
As the particles grow, their Stokes number St and inward
drift velocity ∣vr∣ increase in accordance with Equations (7) and(11). We ﬁnd that the radial drift becomes appreciable when

(14)

Here, Re t = n nmol is the turbulence Reynolds number, where
nmol = vth lmfp 2 is the molecular viscosity. We obtain therelative velocity between the solid materials and the gas, Dvdg ,by setting St1 = St and St2  0 in the above expressions forthe relative velocities.
2.3. Parameter Choice
Table 1 summarizes the parameter range explored in thisstudy. The gas infall rate Ṁg onto the circumplanetary disk istaken to be either 0.02 MJ Myr-1 or 0.002 MJ Myr-1. As
mentioned in Section 2.1, we do not consider a higher value ofṀg since the disk would become too hot for icy satellites toform. In reality, a giant planet carves a gap around its orbit. For
example, the hydrodynamical simulations by Kanagawa et al.
(2015) show that the gas surface density inside the gap isdepleted by a factor of more than 100 compared to that outside
4The Astrophysical Journal, 846:81 (10pp), 2017 September 1

Shibaike et al.

overcome this barrier of the St = 1 for growing to become
satellitesimals. After they achieve St = 1, the drift speed
becomes slower and they become jammed. The jam makes thecollisional rate higher, so the collisional growth speeds up. We
also found that this condition St > 1 is consistent withtgrow < tdrift when St ~ 1. Note that we did not consider the
possibility that the dust surface density near the snowline
increases because of sublimation or recondensation (e.g., Saito& Sirono 2011; Ros & Johansen 2013; Ida & Guillot 2016; Ida
et al. 2016; Schoonenberg & Ormel 2017).
3.2. Effects of the Dust and Gas Inﬂow Mass Fluxes
The amount of gas and dust that ﬂows to the circumplanetary
disk can be changed by the conditions of the central planet, theprotoplanetary disk, and the circumplanetary disk. We
investigate the effects of changing the gas and dust inﬂow
mass ﬂuxes. Figure 2 shows that the the dust surface density
increases with the dust-to-gas inﬂow mass ﬂux ratio M˙ d M˙ g.
The radius and Stokes number of the dust particles also havethe same features. The particles can grow to satellitesimals only
when M˙ d M˙ g = 1. This can be understood by using theapproximate analytical expression for the Stokes number(dotted lines in the ﬁgure). When a = 10-4 and Ṁg =0.02 MJ Myr-1, the gas surface density is so large that
Re p  1. For example, we found that the particle Reynolds
number Rep is about 103 when M˙ d M˙ g = 0.1 - 1 at r ~ 10 RJ .
In this case, the dimensionless coefﬁcient CD can beapproximated as a constant, CD » 0.5 (Newton’s friction law,see Equation (8)). The dust–dust and dust–gas relative
velocities can also be approximated as Dvdd » (1 2) vr andDvdg » vr because the turbulence is so week that Dvt is muchsmaller than Dvr (see Figure 5). The approximated Stokes
number can then be described as⎛ M˙ d M˙ g ⎞2 5 ⎛ a ⎞1 5⎟
St » 1.2 ⎜⎟ ⎜⎝ 1 ⎠ ⎝ 10-4 ⎠⎛ T ⎞-2 5 ⎛ Mcp ⎞2 5 ⎛ r ⎞-2 5
⎟
´⎜
,
⎜
⎟ ⎜
⎟
⎝ 160 K ⎠⎝ 1 MJ ⎠ ⎝ 10 RJ ⎠(15)

for small r. Equation (15) is derived by substituting Equation (6)into Equation (5) and integrating it. Here, the gas surface density
and the midplane temperature have been approximated as
Sg » M˙ g WK / (3pacs2 ) and T » (3GMcp M˙ g /(8psSB r 3))1 4 . The
scale height and radial drift velocity of the dust particles have also
been approximated as Hd » Hg (a St)1 2 and vr » -2Sthv K .
Equation (15) shows that the Stokes number is proportional to(M˙ d M˙ g )2 5 and reaches unity when M˙ d M˙ g = 1 and r = 10 RJ .
Once St exceeds unity, the radial drift velocity starts to decrease
with increasing particle size, and hence the particles growto satellitesimals. Neglecting the weak Ṁg dependence of
1 4
T (T µ M˙ g ), Sg is proportional to Ṁg . When Sg is high, thegas drag force that the dust particles receives is strong and St issmall. The collision rate of the dust particles becomes high whenṀd (i.e., Sd ) is large. The high collision rate promotes
satellitesimal formation. Even if the Newton’s friction law does
not apply, the Stokes number is proportional to (M˙ d M˙ g )2 3 or(M˙ d M˙ g )6 11 and the trend that the Stokes number is an increasing
function does not change (see Equations (16) and (17) in the
Appendix).

Figure 2. Steady-state proﬁles of the surface density Sd (top panel) and radius
Rd (middle panel) of dust particles that grow and drift in the circumplanetary
disk of Ṁg = 0.02 MJ Myr-1 and a = 10-4 . The bottom panel shows the
Stokes number St of the particles. The red, green, blue, and purple curves
correspond to M˙ d M˙ g = 1, 0.1, 0.01, and 0.001, respectively. The dashed lines
in the bottom panel show the prediction from the analytic estimate given by
Equation (15). Shaded gray is the region interior to the snowline, which liesat r = 10RJ .

their drift timescale tdrift = r ∣vr∣ becomes shorter than 30 times
the growth timescale tgrow = md (dmd dt ), in agreement withthe situation for dust evolution in protoplanetary disks(Okuzumi et al. 2012; Tsukamoto et al. 2017). For
M˙ d M˙ g = 1, we ﬁnd that the particles stop drifting and growto kilometer-sized satellitesimals at r ~ 10 RJ . The drift stalls
because the drift speed (normalized by hv K ) decreases withincreasing size as long as St > 1. Therefore, they have to
5The Astrophysical Journal, 846:81 (10pp), 2017 September 1

Shibaike et al.

Figure 4. Steady-state proﬁles of the Stokes number St of the particles whenṀg = 0.02 MJ Myr-1, and a = 10-5 (upper panel) and a = 10-2 (lower
panel). The dashed lines in the upper and lower panels show the predictions
from the analytic estimates given by Equations (15) and (17), respectively.

St for Ṁg = 0.002 MJ Myr-1, in the Appendix). Nevertheless,
we ﬁnd that dust particles grow beyond St = 1 only
when M˙ d M˙ g = 1. For ﬁxed M˙ d M˙ g, the Ṁg dependence of
St is indeed weak. When Ṁg = 0.02 MJ Myr-1, only T depends
-1 10
1 4
(Equation (15)).
on Ṁg (T µ M˙ g ), so St µ T -2 5 µ M˙ g
1
2
When Ṁg = 0.002 MJ Myr
or a = 10 , St µ T -1 ´
4 11
5 44
4
9
7
36
M˙ g µ M˙ g
µ M˙ gor St µ T -1 ´ M˙ g(Equations (16)and (17)).
Figure 3. Same as Figure 2, but for Ṁg = 0.002 MJ Myr-1 and a = 10-4 . The
dashed lines in the bottom panel show the prediction from the analytic estimate
given by Equation (16). Shaded gray is the region interior to the snowline,
which lies at r = 5RJ .

3.3. Effects of the Strength of Turbulence
The strength of turbulence in the circumplanetary disk isalso a key parameter of the dust evolution. Figure 4 represents
the proﬁles of the dust particles in the case withṀg = 0.02 MJ Myr-1, and a = 10-5 (upper panel) anda = 10-2 (lower panel). In the case of a = 10-5, the Stokes
number of drifting particles is on average lower than those in
the ﬁducial case with a = 10-4 . Even if M˙ d M˙ g = 1, the dust
particles cannot grow to satellitesimals outside of the snowline
at r = 10 RJ . In the case of a = 10-2 , the Stokes number isslightly higher than that in the ﬁducial case. Equations (15)and (17) show that the Stokes number is actually proportional toa1 5 or a1 11 when Ṁg = 0.02 MJ Myr-1 (see also Equation (16)for Ṁg = 0.002 MJ Myr-1). However, the Stokes number is not

The key parameter of the dust evolution is not the pure gas
inﬂow mass ﬂux but the ratio of the dust and gas inﬂow mass
ﬂuxes. Figure 3 represents the distribution of the dust surface
density, dust radius, and Stokes number of the representative
dust particles for Ṁg = 0.002 MJ Myr-1 and a = 10-4 . In thiscase, the snowline lies at r = 5RJ (see Figure 1). The proﬁles ofthe surface density and radius of dust particles are lower thanthose for Ṁg = 0.02 MJ Myr-1 and a = 10-4 (Figure 2).
The radial proﬁles of St are steeper than those in the ﬁducial
case (we derive Equation (16), an analytic equation of
6

The Astrophysical Journal, 846:81 (10pp), 2017 September 1

Shibaike et al.

Figure 5. Dust–dust relative velocities with different turbulence strengths, a = 10-4 (left panel) and a = 10-2 (right panel). The other conditions are M˙ d M˙ g = 1 andṀg = 0.02 MJ Myr-1 in both panels. The black curves represent the dust–dust relative velocities (collision velocities). The red and blue curves represent the velocities
induced by only the turbulence and their radial drift, respectively. The black dashed lines are the critical velocity of fragmentation.

high enough for the particles to overcome the radial drift outside
the snowline unless M˙ d M˙ g  1.
Even if M˙ d M˙ g = 1, satellitesimals would not form via direct
dust growth because the relative velocity between the dust
particles would be too high to avoid collisional fragmentation.
Collision simulations by Wada et al. (2009) argued that icy dust
aggregates with monomers of 0.1 m m fragment upon collision
if the collision velocity is higher than 50 m s-1. Figure 5
represents the dust–dust relative velocities for different valuesof α. When a = 10-4 , the relative velocity is determined by
the radial drift speed and that induced by the turbulence is low.
When a = 10-2 , the dust–dust relative velocity is determinedby the turbulence because the relative velocity induced by the
turbulence is proportional to a (Equation (14)) and it
becomes 10 times higher than that for a = 10-4 . Figure 5
shows that the relative velocity exceeds 50 m s-1, indicatingthat collisional fragmentation would happen. Therefore,
satellitesimal formation via direct dust coagulation is unlikely
to occur in such strong turbulence. Note that experiments byGundlach & Blum (2015) showed that the fragmentation
occurs with a collision speed of ~10 m s-1 for icy aggregateswith monomers of 1 m m .

Figure 6. Condition for satellitesimal formation when Ṁg = 0.02 MJ Myr-1.
The green ticks indicate that dust particles grow to satellitesimals outside of thesnowline at r = 10 RJ . The red crosses indicate that the radial drift barrier
inhibits dust growth to satellitesimals. The blue triangles indicate that dust
particles grow to satellitesimals in the calculations, but the collision velocity(dust–dust relative velocity) is faster than the critical velocity of fragmentation,
50 m s-1. The dashed lines show the condition St = 1 at r = 10 RJ from
Equations (15) and (17). When Ṁg = 0.002 MJ Myr-1, the condition is thesame except that the drift barrier is overcome outside the snowline if
M˙ d M˙ g = 1 and a = 10-5 (asterisk).

3.4. Conditions for Satellitesimal Formation

4. Discussions

The results presented in the previous subsections can besummarized in Figure 6. This ﬁgure represents the condition
for satellitesimal formation when Ṁg = 0.02 MJ Myr-1. The
conditions are M˙ d M˙ g  1 and 10-4  a < 10-2 . The condition for breaking through the radial drift barrier is approximately given by M˙ d M˙ g > 6 ´ 10-3 ´ a-1 2 derived fromthe condition St > 1 at r = 10 RJ (see Equation (15)). Whenthe turbulence is strong (a  10-3), it is about M˙ d M˙ g >0.08 a-1 6 (see Equation (17)). The dashed lines in Figure 6
show the boundary of each condition. However, in the case of
a = 10-2 , the aggregate collision velocity is too high to avoid
collisional fragmentation. When Ṁg = 0.002 MJ Myr-1, thedrift barrier is overcome outside the snowline even if a = 10-5
because the line is at r = 5 RJ (asterisk in Figure 6).

4.1. Feasibility of the High Dust-to-gas Inﬂow Mass Flux Ratio
We found that the one of the conditions for satellitesimalformation is M˙ d M˙ g  1. However, this condition may be
difﬁcult to achieve. First, since the dust particles tend to settle
down toward the midplane, the inﬂow gas from the highaltitude is likely to be dust-poor gas (Tanigawa et al. 2012).
This effect must depend on the conditions of the turbulence andthe gas density of the region around the circumplanetary diskthat the accretion gas comes from (Equations (15) and (10)).
Second, the dust supply may not be enough to achieveM˙ d M˙ g  1. Dust particles are drifted from the outer region of
the protoplanetary disk. However, these particles have alreadygrown to the size of pebbles (centimeter-sized particles) until
7

The Astrophysical Journal, 846:81 (10pp), 2017 September 1

Shibaike et al.

Figure 7. Steady-state proﬁles of the Stokes number St of the highly porous (rint = 1.4 ´ 10-4 g cm-3) dust particles with a = 10-4 , and Ṁg = 0.02 (left panel) and
0.002 MJ Myr-1 (right panel). The dashed lines in both panels show the predictions from the analytic estimates given by Equation (15). We stopped the calculation
when the particle radius reached 100 km (the red curve in the left panel).

they reach around gas planets like Jupiter (e.g., Lambrechts &Johansen 2012; Okuzumi et al. 2012; Sato et al. 2016), so mostof them should be dammed at the outer edge of the gas gap by
the positive gas pressure gradient (e.g., Adachi et al. 1976; Zhu
et al. 2012). In this case, only a small part of the dust particles
can penetrate into the gas gap and ﬂow into the circumplanetarydisk, so the dust-to-gas mass inﬂow ﬂux ratio should be smaller
than unity.
One possibility for achieving the high ratio is by considering
satellitesimal formation in the ﬁnal phase of planetary formation.
Photoevaporation may increase the dust-to-gas mass ratio in
protoplanetary disks as time passes (e.g., Alexander et al.
2006a, 2006b). It is also considered that the gas ﬂux decreases
in the ﬁnal phase because the gas gap becomes wider and deeper(e.g., Kanagawa et al. 2015; Tanigawa & Tanaka 2016). Our
results actually suggested that a low gas inﬂow mass ﬂux issuitable for satellitesimal formation. The midplane temperature T
1 4
is almost proportional to Ṁg (Equation (3)). When the gas
inﬂow decreases, the disk becomes cooler and the snowline
moves inward (see Figure 1). This means that the area where icy
satellitesimals can form expands. Moreover, the collisional
velocity driven by turbulence weakly depends on the gas inﬂow
1 8rate, Dvt µ cs µ T1 2 µ M˙ g (Equation (14)). A low gas inﬂow
mass ﬂux may also contribute to overcoming the fragmentation
barrier.

We note that the growth timescale (i.e., St ) depends on theinternal density only in the Stokes regime (Okuzumi et al.
2012). Since the gas densities of protoplanetary disks are
generally much lower than those of circumplanetary disks, thegrowth timescale of highly porous dust aggregates in protoplanetary disks can be so small that they overcome the drift
barrier within the Stokes regime (Okuzumi et al. 2012).
4.3. Streaming Instability
Generating growing particle-density perturbations by streaming
instability is another planetesimal formation mechanism notapplicable to the collisional growth of the dust particles (Youdin& Goodman 2005). The difference between the velocities of thedust and gas drives the instability. The dust particles areconcentrated quickly in localized dense clumps, so they do not
drift to the Sun. This mechanism may also be applicable tosatellitesimal formation. Carrera et al. (2015) showed that thecondition that streaming instability is active depends on the dustto-gas surface density ratio and the Stokes number of the dust
particles. They found that the particle clumps can form if
Sd Sg > 0.02 in the most suitable Stokes number condition (seeFigure 8 in Carrera et al. 2015). However, our results showed thatthe ratio is much lower than the critical value (see Figure 8), so it
should be difﬁcult to form satellitesimals via the streaming
instability process. Note that water sublimation or recondensationnear the snowline can result in an enhancement in the dust surface
density and may be able to trigger streaming instability (Ida &Guillot 2016; Ida et al. 2016; Schoonenberg & Ormel 2017).

4.2. Effects of the Internal Density
We investigated the impact of changing the internal density of
dust particles on satellitesimal formation. Figure 7 represents the
Stokes number for rint = 1.4 ´ 10-4 g cm-3 and a = 10-4 . We
found that the conditions for satellitesimal formation did notchange from those with rint = 1.4 g cm-3 (see Figures 2 and 3).
This is because the Stokes number for small r can also beapproximated as Equation (15) in this case, (we found that
Re p  10 3 for r ~ 10 RJ and we have been able to assume
CD » 0.5) and this approximated Stokes number dose not dependon rint . Generally, the growth timescale takes a minimum value
within the Newton regime (Okuzumi et al. 2012). Therefore, theStokes number does not grow beyond the dashed lines in Figure 7
by changing rint unless it reaches unity. The impact of changing
rint is only that the ﬂuffy particles move earlier than the compact
particles from the Stokes regime (Re p  1) to the Newton regime.

4.4. Validity of the Single-size Approachin Circumplanetary Disks
We used a single-size approach to investigate the growth anddrift of dust particles in circumplanetary disks. As shown by
Sato et al. (2016), this approach is valid for the growth and driftof mass-dominating particles in protoplanetary disks, as long as
the collisional destruction of the particles is negligible (see also
Krijt et al. (2016), Okuzumi et al. (2016), and Tsukamoto et al.
(2017) for applications to dust growth in protoplanetary disks).
In principle, the single-size approximation breaks down when
there is more than one population of particles that dominates thetotal dust mass. Unlike our assumption that the dust inﬂow onto
the circumplanetary disk is concentrated at the outer edge of the
8The Astrophysical Journal, 846:81 (10pp), 2017 September 1

Shibaike et al.

to understand these effects because satellitesimals can form only
when the Stokes number achieves unity. From these parameterstudies, we have revealed the conditions for satellitesimalformation and discussed their feasibility. We have also discussed
the effects of changing the internal density of dust particles andthe possibility of satellitesimal formation by streaming instability. Our ﬁndings are summarized as follows.
1. The mass ﬂux of dust and gas ﬂowing into thecircumplanetary disk determines the evolution of dust
particles. Particularly, the dust-to-gas inﬂow mass ﬂux
ratio, M˙ d M˙ g , is important. A larger dust inﬂow provides
a higher dust density in the circumplanetary disk,
accelerating collisional dust growth. A small amount ofgas inﬂow makes the gas surface density of thecircumplanetary disk lower and the gas drag force that
dust particles receive weaker. This reduces the radial drift
speed of the particles and this is advantageous forcollisional growth. The approximated Stokes number isactually proportional to (M˙ d M˙ g )2 5, (M˙ d M˙ g )2 3, or(M˙ d M˙ g )6 11 (Equations (15)–(17)). As a result, satellitesimals can form only for M˙ d M˙ g  1. For ﬁxed M˙ d M˙ g,
the Stokes number depends on the gas inﬂow mass ﬂux
1 10
7 36
5 44
(Equations (15)–(17)), so
only: M˙ g , M˙ g , or Ṁg
the key parameter is not the pure gas (or pure dust) inﬂow
mass ﬂux but the ratio of the dust and gas mass inﬂow.
2. The strength of turbulence also affects the fate of dust
particles. When Ṁg = 0.02 MJ Myr-1, the approximated
Stokes number is proportional to a1 5 or a1 11(Equations (15) and (17)). As a result, satellitesimals
can form only for a  10-4 . Moreover, strong turbulence
increases the collisional velocity and it actually exceeds
the critical fragmentation velocity of 50 m s-1 (Wadaet al. 2009) when a = 10-2 (Figure 5).

Figure 8. Steady-state proﬁles of the dust-to-gas surface density ratio Sd Sg
with Ṁg = 0.02 and a = 10-4 . The dashed line shows the occurrence conditionof streaming instability in the most suitable case of St (Carrera et al. 2015).

gas inﬂow region, the dust inﬂow may be extended over a wide
area of the disk. In this case, the size distribution in the inner disk
regions may have two peaks, one for the drifting pebbles that
accreted in the outer disk region and one for the small dust grains
directly supplied to the inner disk regions. Actually, when the dust
particles strongly couple with the gas, the mass ﬂux of the dust
inﬂow should be proportional to that of the gas f µ r -1. The ratioof the dust mass ﬂux ﬂowing to the inside of r relative to the total
dust inﬂow mass ﬂux is then »r r b , which is still about 0.5 even
if r = 15 RJ (the orbit of Ganymede). Therefore, the actual size
distribution of the dust particles may be wide. We note that
our assumption of a concentrated dust inﬂow may have also
caused overestimation of the Stokes number, especially in the
2 5
2 3
6 11
in our model; seeouter disk regions (St µ M˙ d , M˙ d , or ṀdEquations (15)–(17)) and the satellitesimal formation could be
harder in reality.
On the other hand, even if the dust particles have a wide size
distribution, the smaller particles should grow rapidly and the
distribution will narrow. This is because the dust growthtimescale tgrow must be an increasing function of the dust mass
md in the Newton regime (tgrow µ md1 4 or md1 3 in our model;
see Equations (5)–(7), and the approximations in Section 3.2
and the Appendix), so the dust particles should grow as orderly
growth in most areas within the disk. However, the md
dependence of tgrow is weak and in the future we will have to
make sure of the validity of the single-size approach by using
full-size calculations of dust growth in circumplanetery disks.

In summary, the conditions for satellitesimal formation via
collisional dust growth are M˙ d M˙ g  1 and 10-4  a < 10-2
when Ṁg = 0.02 MJ Myr-1 (Figure 6). This result does notstrongly depend on the gas inﬂow mass ﬂux.
Our results also suggest that the porosity of dust particles does
not affect the condition for satellitesimal formation (Section 4.2).
The Stokes number of dust particles takes the largest value in the
Newton regime and the value does not depend on the internal
density (Equation (15)). Low internal density just makes theparticles go into the Newton regime quickly.
In reality, it would be difﬁcult to achieve the conditionM˙ d M˙ g  1 (Section 4.1). In protoplanetary disks, the majorityof dust in the vicinity of gas giants should settle down to the
midplane so that the inﬂow gas from the high altitude may be
dust-poor. The dust particles also have already grown topebble-sized particles. Such large particles may not accrete
onto circumplanetary disks, because they are easily trapped at
the edges of planet-carved gas gaps (Zhu et al. 2012).
Satellitesimal formation via streaming instability is also unfeasible(Section 4.3). The dust-to-gas surface density ratio Sd Sg is muchsmaller than the critical value 0.02 where the instability occurs.
However, the photoevaporation may contribute to satisfying
the condition for satellitesimal formation by increasing thedust-to-gas ratio of the inﬂow gas (Section 4.1). The ﬁnal phaseof planetary formation should also be suitable for satellitesimalformation because of the wide and deep gas gap. The smaller
the gas inﬂow mass ﬂux becomes, the lower temperature andthe slower collision velocity become.

5. Conclusions
We have investigated whether icy dust particles can form
satellitesimals by their pairwise collisional growth in circumplanetary disks with various conditions. We have calculated the
distribution of the surface density, radius, and Stokes number ofthe peak mass dust particles. Our model considered only steadyconditions and assumed that the temperature of the circumplanetary disk is almost like a minimum estimate. The gas inﬂow
and the circumplanetary disk model were based on the results of
Tanigawa et al. (2012) and Fujii et al. (2014). We have changed
the dust and gas inﬂow mass ﬂuxes and the strength of
turbulence in the disk so that we understood the effects of thesefactors on satellitesimal formation. We have also approximatedthe Stokes number (normalized stopping time) of dust particles
9The Astrophysical Journal, 846:81 (10pp), 2017 September 1

Shibaike et al.

The delivery of planetesimals from the protoplanetary disk
may be necessary for satellite formation around gas giants. As
we mentioned in Section 1, the planetesimal capture hypothesis
can roughly reproduce the initial radial distribution ofsatellitesimals assumed in the previous satellite formation
models (Suetsugu & Ohtsuki 2017). Moreover, the captured
planetesimals may efﬁciently accrete the drifting dust particles
that could not grow to satellitesimals in the circumplanetary
disk. Such a satellite formation scenario via pebble accretion
will be investigated in our future work.

ORCID iDs
Satoshi Okuzumi https://orcid.org/0000-0002-1886-0880
Takanori Sasaki https://orcid.org/0000-0003-1242-7290
References
Adachi, I., Hayashi, C., & Nakazawa, K. 1976, PThPh, 56, 1756
Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006a, MNRAS, 369, 216
Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006b, MNRAS, 369, 229
Alibert, Y., Mousis, O., & Benz, W. 2005, A&A, 439, 1205
Blum, J., & Wurm, G. 2008, ARA&A, 46, 21
Canup, R. M., & Ward, W. R. 2002, AJ, 124, 3404
Canup, R. M., & Ward, W. R. 2006, Natur, 441, 834
Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
D’Angelo, G., Henning, T., & Kley, W. 2002, A&A, 385, 647
D’Angelo, G., & Podolak, M. 2015, ApJ, 806, 203
Estrada, P. R., Mosqueira, I., Lissauer, J., D’Angelo, G., & Cruikshank, D.
2009, in Europa, ed. R. T. Pappalardo, W. B. McKinnon, & K. Khurana(Tucson, AZ: Univ. Arizona Press), 27
Fujii, Y. I., Okuzumi, S., Tanigawa, T., & ichiro Inutsuka, S. 2014, ApJ, 785,
101
Fujita, T., Ohtsuki, K., Tanigawa, T., & Suetsugu, R. 2013, AJ, 146, 140
Fung, J., & Chiang, E. 2016, ApJ, 832, 105
Gundlach, B., & Blum, J. 2015, ApJ, 798, 34
Hayashi, C. 1981, PThPS, 70, 35
Ida, S., & Guillot, T. 2016, A&A, 596, L3
Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72
Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2015, ApJL, 806, L15
Kobayashi, H., Ormel, C. W., & Ida, S. 2012, ApJ, 756, 70
Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. 2016, A&A, 586,
A20
Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32
Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001
Lunine, J. I., & Stevenson, D. J. 1982, Icar, 52, 14
Miguel, Y., & Ida, S. 2016, Icar, 266, 1
Mosqueira, I., Estrada, P., & Turrini, D. 2010, SSRv, 153, 431
Mosqueira, I., & Estrada, P. R. 2003a, Icar, 163, 198
Mosqueira, I., & Estrada, P. R. 2003b, Icar, 163, 232
Nakamoto, T., & Nakagawa, Y. 1994, ApJ, 421, 640
Ogihara, M., & Ida, S. 2012, ApJ, 753, 60
Okuzumi, S., Momose, M., iti Sirono, S., Kobayashi, H., & Tanaka, H. 2016,
ApJ, 821, 82
Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106
Ormel, C., & Cuzzi, J. 2007, A&A, 466, 413
Perets, H. B., & Murray-Clay, R. A. 2011, ApJ, 733, 56
Ros, K., & Johansen, A. 2013, A&A, 552, A137
Saito, E., & Sirono, S. 2011, ApJ, 728, 20
Sasaki, T., Stewart, G. R., & Ida, S. 2010, ApJ, 714, 1052
Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15
Schoonenberg, D., & Ormel, C. W. 2017, A&A, 602, A21
Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
Sohl, F., Spohn, T., Breuer, D., & Nagel, K. 2002, Icar, 157, 104
Suetsugu, R., & Ohtsuki, K. 2017, ApJ, 839, 66
Suetsugu, R., Ohtsuki, K., & Fujita, T. 2016, AJ, 151, 140
Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853
Tanigawa, T., Maruta, A., & Machida, M. N. 2014, ApJ, 784, 109
Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47
Tanigawa, T., & Tanaka, H. 2016, ApJ, 823, 48
Tsukamoto, Y., Okuzumi, S., & Kataoka, A. 2017, ApJ, 838, 151
Turner, N. J., Lee, M. H., & Sano, T. 2014, ApJ, 783, 14
Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ,
702, 1490
Ward, W. R., & Canup, R. M. 2010, AJ, 140, 1168
Weidenschilling, S. 1977, Icar, 44, 172
Whipple, F. 1972, in From Plasma to Planet, ed. A. Elvius (London:
Wiley), 211
Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
Youdin, A. N., & Lithwick, Y. 2007, Icar, 192, 588
Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, ApJ,
755, 6We thank the anonymous reviewer for the very usefulcomments. We also thank Yann Alibert, Sebastien Charnoz,
Joanna Drazkowska, Cornelis Dullemond, Paul Estrada, Yuri
Fujii, Ryuki Hyodo, Yamila Miguel, Chris Ormel, Judit
Szulágyi, and Takayuki Tanigawa for valuable discussions.
This work was supported by JSPS KAKENHI grant Numbers
JP15H02065, JP16K17661, JP16H04081, and JP16J09590.
Appendix
Approximation of the Stokes Number for
Low-gas-density Cases
Although we assume that CD is a constant in Section 3.2, thisassumption is not correct when the gas density is low via small
gas inﬂow mass ﬂux or strong turbulence. When Ṁg =0.002 MJ Myr-1 and a = 10-4 , or Ṁg = 0.002 MJ Myr-1 anda = 10-2 , the particle Reynolds number is within the range of
10-1  Re p  10 2 for r ~ 10 RJ . In this case, CD can beapproximated as CD » 12 Rep (Equation (8)). When Ṁg =0.002 MJ Myr-1 and a = 10-4 , the dust–dust and dust–gas
relative velocities are approximated as Dvdd » (1 2) vr andDvdg » vr . The Stokes number St can then be approximated as⎛ M˙ d M˙ g ⎞2 3 ⎛ a ⎞-1 9 ⎛
⎞4 9
M˙ g
⎜
⎟
St » 1.6 ⎜⎟⎜⎟⎝ 1 ⎠ ⎝ 10-4 ⎠⎝ 0.002 MJ Myr-1 ⎠⎞-2 9 ⎛ Mcp ⎞7 9 ⎛ r ⎞-13 9⎛ T ⎞-1⎛
r int
⎟ ⎜
,
´⎜
⎜
⎟ ⎜
⎟
⎟
⎝ 90 K ⎠ ⎝ 1.4 g cm-3 ⎠⎝ 1 MJ ⎠ ⎝ 10 RJ ⎠(16)

for small r. Unlike in Equations (15) and (17), St decreases withincreasing a, although the dependence is very weak. Strong
turbulence diffuses dust particles into the vertical direction andthereby reduces their collision rate (see Equations (5) and (10)).
However, this effect is canceled out by the particle collisionvelocity induced by turbulence, which increases with increasinga. When Ṁg = 0.02 MJ Myr-1 and a = 10-2 , the two relative
velocities are determined by the strength of turbulence (i.e.,
Dvdd » Dvdg » Dvt » 3a cs St11 2 , see Figure 5),⎛ M˙ d M˙ g ⎞6 11⎛ a ⎞1 11⎛
⎞4 11
M˙ g
⎜
⎟
St » 0.73 ⎜⎟⎜⎟⎝ 10-2 ⎠ ⎝ 0.02 MJ Myr-1 ⎠⎝ 0.1 ⎠⎞-2 11⎛ Mcp ⎞9 11⎛ r ⎞-15 11⎛ T ⎞-1⎛
r int
⎟ ⎜
,
´⎜
⎜
⎟
⎜
⎟
⎟
⎝ 160 K ⎠ ⎝ 1.4 g cm-3 ⎠⎝ 1 MJ ⎠⎝ 10 RJ ⎠(17)

for small r.

10