Delivery of Dust Particles from Protoplanetary Disks onto Circumplanetary Disks of Giant Planets

The principal regular satellites of gas giants are thought to be formed by the accumulation of solid materials in circumplanetary disks (CPDs). While there has been significant progress in the study of satellite formation in CPDs, details of the supply of satellite building blocks to CPDs remain unclear. We perform the orbital integration of solid particles in the protoplanetary disk (PPD) approaching a planet, considering the gas drag force by using the results of three-dimensional hydrodynamical simulations of a local region around the planet. We investigate the planetary mass dependence of the capture positions and the capture rates of dust particles accreting onto the CPD. We also examine the degree of dust retention in the accreting gas onto the CPD, which is important for determining the ratio of the dust-to-gas inflow rates, a key parameter in satellite formation. We find that the degree of dust retention increases with increasing planetary mass for a given dust scale height in the PPD. In the case of a small planet (M p = 0.2M Jup), most particles with insufficient initial altitudes in the PPD are isolated from the gas in the accreting region. On the other hand, in the case of a massive planet (M p = 1M Jup), dust particles can be coupled to the vertically accreting gas, even when the dust scale height is about 10%–30% of the gas scale height. The results of this study can be used for models of dust delivery and satellite formation in the CPDs of gas giants of various masses, including exoplanets.

Massive regular satellites, such as the Galilean satellites and Titan, are thought to be formed by accumulation of solids in circumplanetary disks (CPDs), gaseous disks formed around giant planets when the planet accretes gas of the protoplanetary disks (PPDs).In early models of satellite formation around Jupiter and Saturn, their CPDs were assumed to be isolated from the PPD.The total dust mass in each CPD is assumed to be consistent with the total current mass of the satellites in each system, which is a similar concept to the minimum mass solar nebula (MMSN) model (Lunine & Stevenson 1982).Mosqueira & Estrada (2003a,b) extended this model by considering CPDs composed of an optically thick inner region and an optically thin outer region.Satellite formation in such CPDs has been recently examined in detail by Miguel & Ida (2016).They found that the high gas density in the inner disk leads to very rapid migration of satellitesimals, and it is difficult to form the current Galilean satellites in such a CPD.
On the other hand, Canup & Ward (2002) proposed the gas-starved disk model for the CPD, in which gas and solids are continuously supplied from the PPD.Canup & Ward (2006) showed that the gas-starved disk can reproduce the common mass ratio of the satellite system to the host planet observed for giant planets in the solar system.In addition, the hydrodynamic simulations showed that a CPD is naturally formed around a sufficiently massive planet during its gas accretion phase (e.g., Tanigawa & Watanabe 2002;Machida et al. 2008;Tanigawa et al. 2012;Szulágyi et al. 2014Szulágyi et al. , 2016Szulágyi et al. , 2017)).Thus the picture of the gas-starved disk has been widely accepted, and various studies of satellite formation have been conducted with such CPDs (e.g., Canup & Ward 2002, 2006;Sasaki et al. 2010;Ogihara & Ida 2012).Recently, a model of satellite formation in a decretion disk, which has nearly vertical inflow in the inner part and outflow near the midplane, has also been proposed (Batygin & Morbidelli 2020).
In such models of satellite formation, the total amount and the size of solid particles supplied to the CPD are important, because they likely determine the timescale of satellite accretion and the composition of the formed satellites (e.g., Canup & Ward 2002;Ronnet et al. 2017).Fujita et al. (2013) studied capture of kilometer-size and larger planetesimals, and Suetsugu & Ohtsuki (2017) examined their capture and subsequent orbital evolution in the CPD.Such planetesimals can be captured in the inner part of the CPD where the gas density is sufficiently high and the planetesimals suffer strong gas drag.Suetsugu & Ohtsuki (2017) concluded that captured planetesimals would contribute to satellite formation more than the dust supplied with the gas if both the velocity dispersion of the planetesimals and the width of the planetesimal gap in the PPD are sufficiently small.If planetesimals can be captured near the outer edge of the gas gap due to the filtering effect (e.g., Rice et al. 2006), external gravitational perturbation from other planets or planetary embryos would be necessary to deliver them to the CPD (Ronnet et al. 2018).On the other hand, relatively small particles such as pebbles and dust would also be important building blocks of satellites.Shibaike et al. (2019) proposed the so-called slow-pebble-accretion scenario of Galilean satellite formation, where a small number of planetesimals captured by the CPD slowly accretes pebbles supplied into the CPD.The scenario can reproduce most characteristics of the current Galilean satellites if some reasonable parameters are assumed.While many of these studies assumed sufficient solids in the CPD, there are still many unknowns about the supply of solid particles into the CPD.
Since meter-sized and smaller particles are much more susceptible to the effects of gas drag than kilometer-size planetesimals, the gas field around the CPD needs to be taken into account for their delivery into the CPD.High-resolution hydrodynamic simulations showed that only the off-midplane gas in the PPD can accrete onto the CPD and the accretion occurs nearly vertically.In addition, part of the accreted gas flows radially outward through the planet's Hill sphere on the midplane unless the viscosity of the PPD is very small (e.g., Tanigawa et al. 2012;Szulágyi et al. 2014Szulágyi et al. , 2017)), which prevents small particles confined in the midplane of the PPD from accreting into the CPD (Tanigawa et al. 2014).Thus, assistance of vertical inflow is important for the delivery of small particles coupled with the gas to the CPD.Homma et al. (2020) investigated the supply to the CPD of solid particles with various sizes distributed vertically in the PPD, considering the gas field around the planet obtained from the local hydrodynamic simulation of Tanigawa et al. (2012).They found that small particles sufficiently stirred up by turbulence in the PPD can accrete onto the CPD with the vertically accreting gas.
However, Tanigawa et al. (2012Tanigawa et al. ( , 2014) ) and Homma et al. (2020) only examined the case where the Hill radius of the planet is the same as the scale height of the PPD, which corresponds to a planet with the mass of 0.4M Jup (M Jup is the current mass of Jupiter) at 5.2 au in the case of the MMSN, and the planetary-mass dependence was not investigated.Understanding of the planetary-mass dependence is important for the understanding of the satellite formation around planets with various masses, for example, Jupiter, Saturn, and exoplanets.Also, CPDs of exoplanets have been observed recently (e.g., Isella et al. 2019), and a better understanding of the planetary mass dependence is expected to contribute to modeling of solid delivery into the CPD for comparison with the observations.
In our recent paper (Maeda et al. 2022), we performed local three-dimensional high-resolution hydrodynamic simulations for various planetary masses (corresponding to 0.05 − 1M Jup at 5.2 au for the MMSN) with the same numerical code used in Tanigawa et al. (2012) (see also Machida et al. 2008), and revealed the planetary-mass dependence of the source region and mass accretion rate of gas accreting into the CPD.In the present work, we will examine the planetary-mass dependence of the supply of solid particles into the CPD using the gas fields obtained by Maeda et al. (2022).
The strucuture of this paper is as follows.In § 2, we will explain basic equations, gas drag calculation, and numerical settings.In § 3, we will show the results of orbital calculation, such as initial altitudes and radial locations in the PPD of dust particles captured in the CPD, and capture rates of dust particles.In § 4, we will derive the capture rate as a function of dust scale height, and also examine the degree of dust retention in accreting gas into the CPD.We will also derive the dust mass accretion rate onto the CPD and the ratio of dust-to-gas inflow rates by combining our results with the gap models for the surface densities of dust and gas obtained in previous works.We will also discuss model limitations of this study.We will summarize this work in § 5.

Basic Equations
We integrate orbits of particles taking account of the effects of the gas flow perturbed by the gravity of a planet, using the methods described in Homma et al. (2020) (see also Tanigawa et al. 2014).We consider a three-body problem for a central star, a planet, and a particle, assuming that the planet is on a circular orbit.We adopt the rotating local Cartesian coordinate system with origin at the planet, the x-axis pointing radially outward, the y-axis in the direction of the planet's orbital motion, and the z-axis in the vertical direction to the PPD midplane.We neglect the curvature of the protoplanetary disk, and define r and R as r ≡ (x 2 +y 2 +z 2 ) 1/2 and R ≡ (x 2 +y 2 ) 1/2 , respectively.When the masses of the planet and the particle are sufficiently smaller than that of the central star and the orbital eccentricity and inclination are small, the motion of the particle relative to the planet can be described by the Hill equations (Nakazawa et al. 1989;Ohtsuki 2012).The Hill equations can be normalized by the planet's Hill radius r H = (M p /(3M c )) 1/3 a and the inverse of orbital angular velocity of the planet Ω −1 K = (GM c /a 3 ) −1/2 , where M c , M p , and a are the masses of the central star and the planet, and the semi-major axis of the planet, respectively.The normalized Hill equations considering gas drag are described as (Fujita et al. 2013;Tanigawa et al. 2014;Homma et al. 2020): where tildes denote normalized quantities.Using the gas drag force F drag = −(1/2)C D πr 2 d ρ g ∆u∆u and the dust particle mass m d , the acceleration due to the gas drag force a drag = F drag /m d is expressed in the normalized notation as where C D , ρ d , ρ g , r d , ∆u are the gas drag coefficient, bulk density of the particle (we assume ρ d = 10 3 kg m −3 ), gas density, particle radius, and the velocity of the particle relative to the gas, respectively.

Hydrodynamic Simulations
We use the gas fields obtained by Maeda et al. (2022) with local three-dimensional hydrodynamical simulations to calculate the gas drag force for dust particles.The nested-grid method (e.g., Machida et al. 2005Machida et al. , 2006) ) is used for the hydrodynamical simulations to examine the gas flow in the vicinity of the CPD with high resolution while considering a wide area around the planet (x, y ∈ [−12h g , 12h g ], z ∈ [0, 6h g ], where h g is the gas scale height of the PPD).We denote the level of the nested-grid by l, where l = 1 corresponds to the whole simulation box.With each increment of the level of the nestedgrid, the size of the computational box is reduced by a factor of two, keeping the center of the box and the number of cells the same; the number of the cells for all grid levels is (n x , n y , n z ) = (64,64,16).We set the maximum grid level to be 11.
An important parameter in our simulation is the Hill radius of the planet normalized by the scale height of the PPD, denoted by Assuming that M c = M ⊙ and the temperature distribution of the PPD is given by the MMSN model (Hayashi 1981) as we obtain the relationship between the planetary mass M p and the parameter r H /h g as (5) Thus, r H /h g = 0.8, 1.0, and 1.36 correspond to the planetary masses of M p = 0.2, 0.4, and 1.0M Jup , respectively.In this work, we will examine these three cases.

Effects of Gas Drags
As for the gas drag coefficient C D , we use the following interpolation formula that reproduces known expressions in limiting cases (see Tanigawa et al. 2014;Homma et al. 2020): where R = 2r d ∆u/ν is the Reynolds number, and M = ∆u/c s is the Mach number.A correction factor w c is defined so that w c = 0.4 for R < 2 × 10 5 and w c = 0.2 for R > 2 × 10 5 .The kinetic viscosity is given by ν = 0.353 8/πc s l g (Chapman & Cowling 1970), where l g = m mol /(σ mol ρ g ) is the mean free path.We set the mass and collision cross section of a gas molecule to m mol = 3.9×10 −27 kg and σ mol = 2.0 × 10 −19 m 2 , respectively.
In the orbit calculation, we obtain the gas density ρ g and the relative velocity between gas and dust particles ∆u at the position of the particle from the results of hydrodynamic simulation, and calculate the gas drag force at each time step.We assume the gas surface density at Jupiter's orbit in the MMSN model (Hayashi 1981), Σ g,0 = 1434 kg m −2 , for the unperturbed region of the PPD.
While approaching the planet, sufficiently large particles undergo vertical oscillation across the midplane of the PPD due to the vertical component of the tidal force.In the case of the Epstein drag law, the critical radius r d,crit for this vertical oscillation can be expressed as r d,crit = v th ρ g /(2ρ d Ω K ), where v th = 8/πc s is the thermal velocity of the gas.If we assume that the planet is located at 5.2 au from the central star in the MMSN model and the depletion factor of the gas surface density is given by f H (f H = 1 corresponds to the MMSN model), v th = 1.1 × 10 3 m s −1 , and Ω K = 1.7 × 10 −8 rad s −1 , and the gas density at the midplane is ρ g = 1.5×10 −8 f H kg m −3 , thus we obtain r d,crit ≃ 0.5f H m. For example, when a gap is opened around the planetary orbit and f H = 0.1, we have r d,crit ≃ 5 cm.We examined cases with r d = 10 −4 − 10 m, but we will mostly focus on the cases of r d = 10 −4 − 10 −2 m (< r d,crit when f H ≳ 0.02), where the particles are easily coupled to the gas accreting into the CPD.

Settings of Orbital Integration for Solid Particles
We integrate Equation (1) using the eighth-order Runge-Kutta integrator as used in Homma et al. (2020).To focus on effects of the initial vertical distribution of particles, we assume that the particles are initially on circular orbits.This assumption is valid for the case of r d ≤ r d,crit , because eccentricities of such small particles are quickly damped by gas drag (see Figure 11 in Homma et al. 2020).Then, the parameters that determine the initial conditions of particles are the impact parameter b = a s − a and the initial altitude z 0 , where a s is the initial semi-major axis of the particle.We start orbital integration at the azimuthal boundary of the hydrodynamic simulation box; x > 0 and y 0 = L y /2 − L y /n y (Homma et al. 2020).We derive the initial x-coordinate from y 0 and the given impact parameter b as x 0 = b 2 − 8/y 0 (Ida & Nakazawa 1989).While large particles that are not well-coupled to the gas can be trapped at the pressure bump of the gas gap (Zhu et al. 2012;Weber et al. 2018), we consider those particles that enter the gap by some mechanisms and assume that the radial distribution of these particles within the gap is uniform.The actual surface density of particles within the gap should depend on their Stokes number; we will discuss dust accretion into the CPD considering the dust gap based on results of global simulation in § 4.3.
The conditions for the capture of particles and the termination of the orbital integration are the same as in Homma et al. (2020) (see also Tanigawa et al. 2014).Each particle is considered to be captured by the CPD when one of the following two conditions is met: (i) The orbit of an energetically captured particle becomes planet-centered and nearly circular in the CPD: Ẽ < 0, e p < 0.3, ãp < 0.5, |N w | ≥ 3.Here Ẽ is the energy integral given by In the above, e p and a p are the eccentricity and semi-major axis of the planet-centered orbit of the particle, respectively, and are given by (Murray & Dermott 1999): where v iner and j iner are the velocity and the specific angular momentum around the planet in the inertial frame, respectively.N w is the winding number around the planet (Iwasaki & Ohtsuki 2007).
(ii) Even if the particle is not in a circular orbit around the planet, it is considered to be captured if it is captured energetically within the planet's Hill sphere and has orbited around the planet a sufficient number of times: We terminate the orbital integration for each particle when any one of the above capture conditions (i) or (ii), or either of the following conditions is met: (iii) the particle collides with the planet: r < r p , where r p is the physical radius of the planet.We set r p = 10 −3 r H assuming the Jovian orbit (Ida & Nakazawa 1989;Ohtsuki 2012).(iv) The particle is far enough away from the planet: r > y 0 + b + (2e 0 + i 0 )a 0 , where e 0 , i 0 , and a 0 are the eccentricity, inclination, and semi-major axis of the particle's initial heliocentric orbit, respectively.

Examples of Particle Orbits
The orbital behavior of particles from the PPD to the CPD depends on their Stokes number St (Homma et al. 2020).Here we define St by the following equations using the planet-centered-Keplerian angular velocity Ω p near the planet (r < 0.5) and star-centered-Keplerian angular velocity Ω K elsewhere: where the stopping time t stop is described as *1 While we consider the capture of particles in both the prograde and retrograde directions, we find that all the particles captured in the prograde direction because the size of particles is sufficiently small and particles are dragged by the gas flow in the CPD.
Figure 1 shows examples of particle orbits and the change of the Stokes number along them in the case of z s = r H .In the case of r d = 10 m (red), St ≫ 1 and the effect of gas drag is sufficiently small, and the orbit is similar to that without gas drag (black).It oscillates in the vertical direction with the amplitude depending on the initial altitude and can be captured by the CPD when it happens to pass close to the planet.
In the cases of 0.1 ≤ r d ≤ 1 m (light blue and violet), St ∼ 1 is maintained while approaching the Hill sphere.In the case of r d = 0.1 m ≃ r d,crit , St = 0.5 and the particle gradually settles toward the mid-plane.When St < 0.5, particles cannot settle completely and approach the Hill sphere at an altitude slightly above the mid-plane and are captured from above by the CPD.On the other hand, when St ≳ 0.5 − 1, the vertical oscillation is not completely damped, and the particles approach the Hill sphere with a slight oscillation (see also Homma et al. 2020).The degree of oscillation damping depends on the density of the gas along the orbit, i.e., the depth of the gap formed around the planet's orbit.It should be noted that the gap depth obtained by our local simulation tends to be shallower than the one obtained by global simulation for the same planetary mass.Since the capture process for small particles with St ≲ 1 are sensitive to the depth of the gap, it is important to take such effects into account when examining planetary-mass dependence (Maeda et al. 2022, see also § 4.5).
The particles with r d ≤ 1 cm (orange, green, and blue) have St ≪ 1 and approach the Hill sphere while keeping a certain altitude, being more or less coupled to the gas flow.Then they move down nearly vertically with the accreting gas flow over the Hill sphere and become captured by the CPD.
Differences in the gas field around a planet due to differences in planetary mass affect particle orbits in the following ways: First, the impact parameter b of the accretion orbit is varied because of the structure of the accretion bands is modified (see § 3.2).Second, when the planetary mass is increased, the gas density around the planetary orbit tends to be reduced more due to gap opening, resulting in a larger Stokes number for the same particle size.This makes the particles with St ≃ 1 (r d = 10 −1 m) more easily settled toward the midplane while approaching the planet.

Initial Orbits of Captured Particles
In Figure 2, we show the distribution of the initial impact parameter b and initial altitude z s of the particles captured by the CPD.The different colors mean the different sizes of the particles.The results for r d = 0.1 cm (light blue) are shown as examples for the case of St ∼ 1 for comparison with other cases of St ≪ 1.In the case of the particles with St ≪ 1 (r d ≤ 10 −2 m), the initial orbits of captured particles distribute above the midplane continuously in the z-direction in a zonal pattern except for z s ≃ 0. Because particles of this size are greatly affected by the gas flow, the distribution of the dust accretion band is very similar to that of the gas accretion band (Tanigawa et al. 2012;Maeda et al. 2022).We confirm that sufficiently small particles initially near the midplane of the PPD cannot accrete into the CPD due to the radially outward gas flow near the midplane in the Hill sphere (Tanigawa et al. 2014;Homma et al. 2020;Maeda et al. 2022).
In the case of St ≪ 1 (r d ≤ 10 −2 m), the radial width of the dust accretion band measured in units of r H becomes wider with increasing planetary mass.For example, in the case of r d = 0.1 mm (blue), while the width of the dust accretion band is ∆b cap ≲ 0.15r H for M p = 0.4M Jup , ∆b cap ≳ 0.3r H at maximum for M p = 1M Jup .This expansion of ∆b cap is pronounced in the case of smaller particles.This is because the smaller particles are more affected by the gas flow and the width of the gas accretion band is wider for larger planetary masses (Maeda et al. 2022).In the case of M p = 1M Jup , ∆b cap is larger at z s ≳ r H , which is attributed to the fact that the gas accretion band becomes slightly wider at z ≳ r H (see Fig. 7 in Maeda et al. 2022).Also, ∆b cap is almost independent of z s when r d = 0.1 m because the particles settle toward the midplane before they approach the planet.Thus, the range of b that leads to accretion is similar to that of z s ∼ 0, regardless of initial altitudes.The colors show the size of the particle; light blue, orange, green, blue corresponds to r d = 0.1, 10 −2 , 10 −3 , and 10 −4 m, respectively.The vertical axis is plotted at every z s = 0.1r H , but slightly shifted up and down to avoid overlapping of the lines.

Distribution of the Locations of Particle Capture in the CPD
Figure 3 shows the distribution of the coordinates (R, z) of particles when their capture conditions are met in the range 0 < z s ≤ r H for the cases of two different sizes.We show the results for only two cases of r d = 10 −2 m and 10 −4 m.The results for r d = 10 −3 m are similar to the case of r d = 10 −4 m and are not shown here.In agreement with Homma et al. (2020), particles accreting with the gas are widely supplied onto the surface of the CPD.We found that the particles with r d = 10 −4 m tend to be captured at more inner region of the CPD than those with 10 −2 m.This suggests that small particles accreting from high altitudes can be delivered into more inner regions of the CPD with the gas (Homma et al. 2020).On the other hand, if particles in the PPD undergo significant settling in approaching the planet, they are captured in the midplane near the outer edge of the CPD.Although the size of the CPD is commonly r CPD ≃ 0.2r H when r H /h g ≳ 1 (e.g., Figure 4 in Maeda et al. 2022), we found that the Hill-scaled radial positions of particle captured on the CPD are more widely distributed in the case of larger planets.This may reflect the fact that the size of the CPD measured in units of r H and/or the gas density near the outer edge of the CPD is somewhat larger for a larger planetary mass (Figures 3 and 4

Capture Rate of Dust Particles onto the CPD
Next, we calculate the capture rates of particles onto the CPD.The capture rates in this work are obtained by counting the number of captured particles per unit time and unit surface number density of particles in the PPD, and P cap is normalized by r 2 H Ω K .For r d < r d,crit , the capture rate of particles which have initial altitude z s is described by (Fujita et al. 2013;Tanigawa et al. 2014;Homma et al. 2020): where φ cap is defined such that φ cap = 1 for capture orbits and φ cap = 0 otherwise.We adopted d b = 5 × 10 −4 and performed calculations by varying z s in increments of 0.1.
Figure 4 shows the capture rates for r d = 10 −4 m.Although P cap is non-dimensionalized by the Hill-scaling, the overall values of P cap are larger for larger planetary masses, reflecting the increase in the width of the gas accretion band scaled by r H (Figure 2).The increase of P cap at z s ≳ r H for M p = 1M Jup is attributed to the fact that the gas accretion band becomes slightly wider in this altitude range, as mentioned in § 3.2.We also found that P cap = 0 at z s = 0 for all the three cases, owing to the outward gas flow near the midplane in the Hill sphere ( § 3.2).The increase of the capture rate of dust with small z s (≲ 0.5r H ) for M p = 1M Jup can be explained by the extension of the gas accretion band to lower altitudes (Maeda et al. 2022).
Results for larger particle sizes are shown in Figure 5, with those for r d = 10 −4 m with the dashed lines for comparison.We find that the results are quite similar for all the three particle sizes, because in almost all these cases dust particles have sufficiently small Stokes numbers and are coupled well with the gas.As an exception, a significant deviation from the other cases can be seen in the case of r d = 10 −2 m at z s ≳ 1.5h g when M p = 1M Jup .This is because the gas density is very low at such high altitudes (z s ≳ 2h g ), so the Stokes number of the particles is large.Integrating the capture rates per unit surface number density of particles in the PPD (i.e., P cap (r d , z s ) given by Equation ( 12)) over z s with a given vertical distribution of particles, we can calculate the total capture rates of particles into the CPD.The degree of vertical stirring of dust particles in the PPD can be expressed by the dust scale height h d .Homma et al. (2020) considered dust stirring due to turbulence and used the dust scale height h d expressed in terms of the turbulent viscosity parameter α (Shakura & Sunyaev 1973) and the Stokes number of particles St (Youdin & Lithwick 2007).Then they obtained total capture rates and mass accretion rates of dust particles for several combinations of α and St.However, turbulence is not the only mechanism that stirs up dust particles in the PPD.For example, recent studies show that the meridional circulation caused by an embedded planet can stir up dust particles (Bi et al. 2021;Szulágyi et al. 2022).
Therefore, rather than focusing on any specific mechanism for dust particle stirring, in the present work we will obtain the total capture rate P cap,total , treating the dust scale height h d as a parameter.We assume that the vertical distribution of particles is given by the Gaussian distribution with a scale height h d ; Then the total capture rate of particles P cap,total with this vertical distribution is given by Figure 6 shows the plots of P cap,total (h d ).The behavior of P cap,total (h d ) can be roughly explained by the z s -dependence of the width of the dust accretion band (Figure 2), i.e., the total capture rate is larger for a larger planetary mass and a larger dust scale height.Bi et al. (2021) investigated effects of planet-driven meridional circulation on vertical stirring of dust using three-dimensional hydrodynamic simulation considering the back-reaction of dust on the gas.They showed that even in the case of a planet with a Saturnian mass (M p ≃ 0.3M Jup ), dust particles are easily stirred up to h d ≃ 0.6h g .In the case of h d = 0.6h g in Figure 6, we find that P cap,total ≃ 0.4 for M p = 0.2M Jup , ≃ 0.6 for M p = 0.4M Jup , and ≃ 1.2 for M p = 1M Jup , respectively.That is, although P cap,total is normalized by r 2 H Ω K as mentioned above, its value for M p = 1M Jup is about three times larger than for M p = 0.2M Jup as compared to the simple Hill-scaling.As for the dependence on particle size in the size range considered here (0.1 mm -1 cm), the total capture rate is roughly determined by the dust scale height and the planetary mass, suggesting that the particle size dependence is not significant.However, since the dust scale height likely depends on the particle size (or the Stokes number), the total capture rate likely depends on the dust size via the dust scale height.

Degree of Dust Retention in Accreting Gas
Recent theoretical models of satellite formation in the CPD suggest that one of the key parameters for the formation of satellitesimals and satellites in the CPD is the ratio of dust-to-gas inflow rates (Shibaike et al. 2017(Shibaike et al. , 2019)).Using P cap,total we obtained above, we can derive the mass accretion rate of dust into the CPD from (Homma et al. 2020) where Σ d is the dust surface density at the position of the dust accretion band.Now, following Maeda et al. (2022), we consider the gas flowing through the area of the accretion band at (x, y) = (x g,ab , L y /2), where x g,ab is a typical value of the x-coordinate of the accretion band.The velocity of the gas passing through the area of the accretion band is (3/2)x g,ab Ω K in the coordinate system rorating with the planet.Then, the mass accretion rate of the gas into the CPD is expressed as: where D g is the accretion area of the gas (Tanigawa & Tanaka 2016).We obtain wg by averaging the width of the accretion band w g (z) over z with the gas density distribution (Maeda et al. 2022): Similarly, the dust accretion rate into the CPD can be obtained using the accretion area of the dust particles D d and the averaged width of the dust accretion band wd as Then the ratio of dust-to-gas inflow rates χ is given by In the limit of a very small dust size, dust particles are completely coupled to the gas and we have h d = h g and wd = wg , thus D d = D g .In this case, we obtain χ = Σ d /Σ g , i.e., the ratio of dust-to-gas inflow rates is determined by the ratio of dust-to-gas surface density at the accretion band (feeding zone).On the other hand, if the dust-gas coupling is not perfect, h d < h g due to dust settling, and dust is depleted in the gas accreting onto the CPD, resulting in D d /D g < 1 (Tanigawa et al. 2012).
Here, we call D d /D g as degree of dust retention of accreting gas into the CPD.
The surface densities of dust and gas at the accretion band depend on the depth and width of the gap in the PPD, which depends on various parameters, such as the viscosity of the PPD.In order to avoid uncertainty related to such uncertain parameters, we will focus on the degree of dust retention D d /D g to examine the planetary mass dependence of the ratio of dust-to-gas inflow rates in the accretion band.For the range of planetary masses considered in this study corresponding to 0.8 ≤ r H /h g ≤ 1.36, the width of the gas accretion band is approximately given by (Maeda et al. 2022) Substituting Equation ( 23) into Equation ( 17) and setting x g,ab = 2.2r H , we obtain As for the accretion of dust, from Equations ( 15) and ( 19), we obtain From Equations ( 24) and ( 25), we can obtain the degree of dust retention of accreting gas into the CPD as where P cap,total can be obtained from our numerical results.
Figure 7 shows the plots of D d /D g obtained from Equation ( 26).We confirm that D d /D g ≃ 1 when dust particles are perfectly coupled to the gas and h d /h g ≃ 1, while D d /D g decreases with decreasing dust scale height due to insufficient dust stirring in the PPD.The decrease in D d /D g at small dust scale heights is more pronounced for smaller planetary masses, because the accretion band of dust and gas becomes narrower at low altitudes for small planetary masses (Figure 2).On the other hand, for r H /h g = 1.36 (M p = 1M Jup ), D d /D g ≳ 0.1 even when h d /h g = 0.1.When the surface densities of dust and gas at the accretion band are given, we can determine the ratio of dust-to-gas inflow rates into the CPD combining them with the degree of dust retention D d /D g obtained above.The ratio of surface densities of dust and gas at the accretion band depends on the depth and width of gap opened by the planet and the radial position of the accretion band relative to the gap.As an example, in § 4.3, we derive dust mass accretion rates onto the CPD and the ratio of dust-to-gas inflow rates in the case where the surface densities of dust and gas at the accretion band are given using simple models based on results of global hydrodynamic simulation.
The surface density of dust at the accretion band depends on dust size.Particles with St ∼ 1 are easily trapped at the pressure bump of the gap edge and the surface density becomes high (Zhu et al. 2012;Weber et al. 2018), and even planetesimals may form there (Eriksson et al. 2020;Shibaike & Alibert 2020, 2023).Then collision and fragmentation of the planetesimals at the pressure bump can provide a significant amount of dust, which can diffuse into the gap (Kobayashi et al. 2012;Stammler et al. 2023).While recent studies using hydrodynamic simulations including gas and dust components obtain empirical formulae for the depth and width of the dust gap (Zhang et al. 2018), detailed studies of dynamical behavior of dust and gas around the gap is still in progress.Further investigation of the surface density distribution of dust and gas around the planetary orbit will be important to obtain the ratio of dust-to-gas inflow rates into the CPD and derive constraints on satellite formation.

Applications 4.3.1. Dust Mass Accretion Rates
In order to obtain dust mass accretion rates onto the CPD using Equation (15) with our numerical results for P cap,total , we need to assume the dust surface density in the accretion band.Here, we derive dust mass accretion rates using the surface density at the bottom of the gap.Zhang et al. (2018) performed two-dimensional two-fluid hydrodynamic simulations for dust and gas, and obtained empirical formulae for the width and depth of the dust gap.The formulae obtained for the depth of the dust gap is where Σ d,gap and Σ d,peak are the bottom and peak values of the dust surface density distribution, respectively.The values of the fitting parameters C and D for various values of the Stokes number of the dust particles are obtained from the hydrodynamic simulations (see Table 2 in Zhang et al. 2018).We assume α PPD , the efficiency of the viscous angular momentum transport outside the gap, is independent of the dust vertical diffusion.
Figure 8 shows the relationship between the depth of the dust gap and r H /h g .We find that the dust gap becomes deeper with increasing planetary mass and increasing dust size.Now, we obtain the dust mass accretion rate onto the CPD from Equation ( 15).We assume that the dust surface density at the accretion band is equal to the value at the bottom of the dust gap.Then, where Σ d,0 is the dust surface density at the unperturbed region of the PPD.Here, we assume Σ d,0 = Σ d,MMSN f dust , where Σ d,MMSN and f dust are the dust surface density of the MMSN model at 5.2 au (Hayashi 1981) and the depletion factor of the dust surface density (f dust = 1 corresponds to the MMSN model), respectively.We further assume f dust = 1 for simplicity.The ratio between the peak and bottom values of the dust surface density Σ d,gap /Σ d,peak is obtained using Equation ( 27).
Since the dust particles treated in this study are sufficiently small to avoid significant accumulation at the gas pressure bump, we assume f peak = 1.
Figure 9 shows the plots of the dust mass accretion rates as a function of h d /h g for various values of planetary mass, particle size, and the viscosity parameter α PPD of the PPD.It should be noted that we assume f dust = 1 for simplicity in the calculation of the mass accretion rates presented here, the amount of dust in the PPD should be less than the MMSN value (i.e., f dust ≪ 1) in the late stage of planet formation.Therefore, it is unlikely that the core mass will become unnaturally large due to excessive dust supply, as would be obtained by simply multiplying the mass accretion rate shown here by the growth timescale of the planet.
We find that the mass accretion rate and its dependence on dust scale height are nearly independent of the planetary mass for r d = 10 −2 m.As we mentioned in § 4.1, P total,cap in the case of strong vertical stirring of particles (i.e., h d /h g ≃ 1) becomes about three times larger when the planetary mass increases from 0.2M Jup to 1M Jup (Figure 6).In this case, the planet's Hill radius increases by about a factor of 1.7, from r H /h g = 0.8 to 1.36, thus r 2 H in Equation ( 15) increases by a factor of (1.7) 2 ∼ 3. On the other hand, from the upper panel of Figure 8, Σ d,gap in the case of r d = 10 −2 m becomes about one order of magnitude smaller when the planetary mass increases from 0.2M Jup (r H /h g = 0.8) to 1M Jup (r H /h g = 1.36).As a result, the effects of increased planetary mass on the above three quantities that appear on the r.h.s. of Equation ( 15) nearly cancel out, resulting in the mass accretion rate that is nearly independent of the planetary mass in the case of r d = 10 −2 m.On the other hand, in the case of r d = 10 −4 m, the change of Σ d,gap due to the increase in the planetary mass is small, so the effects of increased P total,cap and r H with increased planetary mass are dominant.Therefore, the mass accretion rate significantly increases with increasing planetary mass in this case.In all the cases shown in Figure 9, mass accretion rates are larger for smaller particles due to the shallower gap opening for smaller particles (i.e., Σ d,gap becomes larger, as shown in Figure 8).

Ratio of Dust-to-Gas Inflow Rates
Using the results of this work together with the empirical formulae of the surface densities of gas and dust at the gap bottom (Kanagawa et al. 2015;Zhang et al. 2018), we obtain the ratio of dustto-gas inflow rates onto the CPD.Using the total capture rate P cap,total obtained in the present work and the gas mass accretion rate Mg obtained by the hydrodynamic simulation of Maeda et al. (2022), the ratio of dust-to-gas inflow rates is obtained by the following equation (Homma et al. 2020): where we substituted the dust surface density at the dust gap obtained by Equation ( 27) into Σ d .In the above, M ′ g is the corrected gas surface density at the gap using the empirical formulae obtained by the global hydrodynamic simulation (Kanagawa et al. 2015; see also Maeda et al. 2022), which is obtained as where Σ g,glo and Σ g,loc are the gas surface densities at the gap bottom obtained by the global simulation of Kanagawa et al. (2015) and the local simulation of Maeda et al. (2022), respectively.
Figure 10 shows the ratio of dust-to-gas inflow rates as a function of M p /M Jup .In the case of small dust scale height (h d /h g ≲ 0.1), the ratio of dust-to-gas inflow rates increases with increasing planetary mass.For example, in the case of h d /h g = 0.1 and α PPD = 10 −2 , the ratio of dust-togas inflow rates increases by about one order of magnitude when the planetary mass increases from 0.2M Jup (r/h g = 0.8) to 1M Jup (r/h g = 1.36).This reflects the fact that the dust accretion band extends down to lower altitudes for a large planetary mass ( § 4.2), which significantly influences the dust accretion rates when the vertical stirring of dust particles is limited.On the other hand, for the larger dust scale height (h d /h g ≳ 0.3), the ratio of dust-to-gas inflow rates is nearly independent of planetary mass or even decreases with increasing planetary mass.In this case, the effect of the deeper gap (thus lower dust surface density) for a larger planetary mass is dominant.Thus, the planetary mass dependence of the ratio of dust-to-gas inflow rates into the CPD varies depending on the degree of vertical stirring of dust particles in the PPD.4.3.3.Implication for Satellite Formation Models Shibaike et al. (2017) investigated satellitesimal formation by calculating collisional growth of dust in a one-dimensional steady CPD.They found that satellitesimals can be formed if the ratio of dustto-gas inflow rates is χ ≃ 1.On the other hand, in the case of smaller χ, they found that particles that have grown from small sizes with St ≪ 1 to larger sizes with St ≃ 1 quickly fall to the central planet and cannot form satellitesimals. Figure 10 shows that χ ≲ 10 −2 even when the PPD viscosity parameter is as high as α PPD = 10 −2 and the dust scale height is sufficiently large, indicating that the conditions for the formation of satellitesimals in Shibaike et al. (2017) cannot be satisfied.Recently, Shibaike & Mori (2023) investigated satellitesimal formation in the CPD with magnetic wind-driven accretion and found that such CPDs can relax the condition for satellitesimal formation (χ ≥ 0.02).However, even this relaxed condition seems difficult to satisfy in our results shown in Figure 10.
On the other hand, Shibaike et al. (2019) proposed "slow-pebble-accretion scenario", which can reproduce major characteristics of the current Galilean satellite systems, such as their orbits, masses, and compositions.In this scenario, several planetesimals captured by the CPD grow by accretion of pebbles that are continuously supplied to the CPD, while they slowly move inward in the CPD to form satellites.They assumed that the ratio of dust-to-gas inflow rates is χ = 0.0026, independent of the stages of the growth of the host planet.We plotted the value of χ (= 0.0026) that corresponds to the one assumed in Shibaike et al. (2019) with the dashed lines in Figure 10. Figure 10 shows that the value of χ assumed in Shibaike et al. (2019) can be achieved in cases where α PPD is large and/or the dust scale height is large, although the values of χ that we obtained vary somewhat depending on adopted dust size and planetary mass.It should be noted that the ratio of the dust-to-gas inflow rates would become lower if the depletion of solids in the PPD (i.e., f dust < 1; see Section 4.3.1) in the late stage of giant planet formation is taken into account.
According to the polarization observation of HD 163296 by ALMA, the dust scale height in the gap is estimated to be h d = 0.3 − 0.7h g and the dust size is r d ≲ 10 −4 m (Ohashi & Kataoka 2019).From our results, the condition for satellitesimal formation obtained by Shibaike et al. (2019) can be achieved in cases where α PPD ∼ 10 −2 when f dust = 1 is assumed.Cilibrasi et al. (2018) performed a population synthesis simulation of satellite formation around a Jupiter-mass planet using the results of radiative hydrodynamical simulation obtained by Szulágyi et al. (2017).In their model, dust particles inflow onto the CPD with the same profile as the gas with the ratio of dust-to-gas inflow rates (0.03 − 0.5), and the satellitesimals were assumed to grow via streaming instability.Figure 10 shows that such high ratio of dust-to-gas inflow rates (∼ 10 −2 − 10 −1 ) can only be achieved under limited conditions of α PPD ≳ 10 −2 , h d /h g ≃ 1, and M p ≃ 0.2M Jup when f dust = 1 is assumed.However, it has recently been shown that planet-driven meridional circulation could increase the ratio of dust-to-gas inflow rates to ∼ 10 −2 − 10 −1 (Szulágyi et al. 2022).These conditions of high ratio of dust-to-gas inflow rates could be achieved by considering global flows in the PPD.Bae et al. (2022) estimated the dust-to-gas mass ratio of the CPD candidate of AS 209b as ≲ 2 × 10 −4 − 9 × 10 −4 , suggesting that the dust in the CPD is more depleted than the typical ISM environment.They inferred that the depletion is presumably due to limited supply and/or rapid radial drift of dust within the CPD.Our results support the former, although we cannot rule out the possibility of the latter.It should be noted that the results presented here depend on the assumed dust/gas surface densities, which are described in § 4.3.1.
In the present work, we assume the isothermal gas for simplisity and did not consider sublimation of icy dust or changes in dust composotion.For example, in the case of Jupiter's orbit, substituting a = 5.2 au in Equation (4) yields 123 K.However, the temperature within the CPD may increase due to viscous heating (e.g., Canup & Ward 2002), and effects of shock heating during accretion onto the CPD may also be important.In this case, if icy dust particle sublimate, the amount of solids supplied to the CPD could be smaller than that shown here.In our model, the composition of the dust only affects its assumed bulk density, and the results will not change significantly as long as the densities are similar.Further studies that take account of the changes in temperature and composition during accretion is important for better understanding of dust supply to CPDs (e.g., Ronnet & Johansen 2020).

Model Limitations
In the present work, we investigated delivery of dust particles from the PPD to the CPD by performing orbital integration of dust particles under the influence of the gas flows in the vicinity of a planet obtained by high-resolution local hydrodynamic simulation, but some simplifications were introduced.First, the local hydrodynamic simulation allows us to investigate the gas flow near the planet with high resolution, but as mentioned earlier, it is difficult to accurately represent the gap structures.We used the assumption of the surface densities of gas and dust based on the gap model derived from global simulation (Kanagawa et al. 2015;Zhang et al. 2018) to estimate the mass accretion rates of gas (Maeda et al. 2022) and dust ( § 4.3.1),and the ratio of dust-to-gas inflow rates ( § 4.3.2).However, direct use of global simulation is desirable to obtain realistic gap structures for both gas and dust simultaneously.Second, in our simulation, the gas is assumed to be isothermal and inviscid for simplicity, and radiative transport is not included.While we have focused on dust accretion into the CPD from a dynamical point of view assuming the locally isothermal PPD and CPD (Equation ( 4)) in this work, the temperature distribution in the CPD and its evolution is important for studies of compositional evolution of dust particles and forming satellites in the CPD (e.g., D'Angelo & Podolak 2015; Ronnet & Johansen 2020).Third, despite the assumption that the gas flow is laminar in the unperturbed state, we assume that dust particles are vertically stirred by some mechanism in the PPD.We treated the degree of vertical stirring as a parameter using the dust scale height, which is assumed to be independent of the efficiency of the angular momentum transport of the PPD.Furthermore, although the influence of the gas flow on the orbital motion of dust particles is taken into account, the back reaction that dust exerts on the gas is not considered.
Recently, treating these effects more realistically, global simulations of PPDs with gas-dust twocomponent fluids have been performed (e.g., Binkert et al. 2021;Szulágyi et al. 2022).The mechanisms that can stir up dust particles have also been investigated in detail, such as MRI (Balbus & Hawley 1991) or meridional circulation (Fung & Chiang 2016).Szulágyi et al. (2022) performed three-dimentioanl global hydrodynamic simulation that took into account radiative transport, gas viscosity, and back reactions from dust to gas, and investigated the dust supply to the CPD.Their results showed that vertical stirring of dust particles by spiral wakes excited by the planet's gravity can play an important role in supplying solids to the CPD.Cilibrasi et al. (2023) performed three-dimentioanl global radiative magnetohydrodynamics (MHD) simulations.They found that a meridional circulation exists in the radiative MHD case, but high turbulent viscosity prevents gap opening, leading to a higher accretion rate than the hydrodynamic case.These simulations require a large amount of computer resource, but simulations with sufficiently high resolution will be important to examine, for example, distribution of captured solids in the CPD.It will also be important to clarify the degree of dust retention and its dependence on various parameters using global simulation with sufficiently high resolution.Li et al. (2023) performed three-dimensioanl global simulations of PPDs over gap-opening timescales and explore the planetary mass dependence of mass accretion rates onto the planet.They found that the scaling law is different from the local simulation of Maeda et al. (2022), probably due to the shearing box approximation not considering the global accretion of the PPD.This effect may also affect dust accretion.
Our calculations do not take into account the size distribution of dust particles.Recently, Karlin et al. (2023) performed multiple dust grain size simulations and showed that the dust behavior itself does not change with or without consideration of size distribution.This indicates that the superposition of the results of single-size simulations can be used to obtain the dust accretion rate onto the CPD in the presence of dust size distribution.

CONCLUSIONS
In this work, we performed orbital integration of solid particles approaching a planet, taking account of drag from the gas flow perturbed by the planet obtained by Maeda et al. (2022) using three-dimensional hydrodynamical simulations of a local region around a planet.We investigated the planetary-mass dependence of the orbits, capture positions, and capture rates of dust particles accreting onto the CPD.In addition, we examined the planetary-mass dependence of the degree of dust retention in accreting gas onto the CPD, which determines the ratio of dust-to-gas inflow rates, a key parameter in models of satellite formation in the CPD.The main results are as follows: 1. Similarly to the gas accretion band examined in our previous work (Maeda et al. 2022), the initial orbits in the PPD of particles captured in the CPD distribute above the midplane continuously in the z-direction in a zonal pattern except near the midplane.For dust particles with St ≪ 1, the radial width of the dust accretion band measured in units of the planet's Hill radius becomes wider with increasing planetary mass.
2. Dust particles accreting with the gas are supplied into a radially extended region of the surface of the CPD.Smaller particles tend to be captured in more inner regions of the CPD.Also, dust particles are captured in radially wider regions of the CPD in the case of a larger planetary mass.
3. The degree of dust retention in accreting gas onto the CPD increases with increasing planetary mass for a given dust scale height in the PPD.For the case of a small planet (M p = 0.2M Jup ), particles with insufficient altitude in the PPD cannot get assisted by the accreting gas flow and cannot accrete onto the CPD.On the other hand, for the case of a massive planet (M p = 1M Jup ), even in the case with a small dust scale height of h d = 0.1h g , the ratio of dust-to-gas inflow rates in the case of particles with r d = 10 −4 to 10 −2 m can be as large as 10 − 30% of the dust-to-gas mass ratio in the PPD (Figure 7).
The results of this study can be used to develop models of satellite formation around gas giants of various masses, e.g., Saturn, Jupiter, and exoplanets.Further investigation on the delivery of dust particles onto the CPD using improved numerical models, such as global simulation of dust and gas components taking into account of their back reactions will be desirable for better understanding of dust supply onto the CPD ( § 4.5).Also, the amount of dust supplied to the CPD is highly dependent on the dust scale height in the PPD, the structures of the gas and dust gaps, and the size distribution of dust particles in the PPD.Further observational studies on these quantities as well as their planetary mass dependence will be desirable to derive more information.In particular, given that giant planets are expected to cause meridional circulation (Teague et al. 2019;Bi et al. 2021;Szulágyi et al. 2022), it is highly important to observe PPDs in the late stages of giant planet formation.Recently, CPD candidates around exoplanets have been detected, such as PDS70c and AS 209b (Isella et al. 2019;Bae et al. 2022), and a comparison of model calculations of dust evolution in the CPD with observations of dust thermal emissions is also being attempted (Shibaike & Mordasini, submitted).These theoretical and observational studies will contribute to the development of more realistic satellite formation models.

Figure 1 .
Figure1.Examples of particle orbits and change of the Stokes number along the orbits in the case of ( b, zs ) = (2.3,1.0) and (2.4,1.0) for M p = 0.4M Jup (top panel) and M p = 1M Jup (bottom panel), respectively.The colors show the size of particles; red, violet, light blue, orange, green, and blue correspond to r d = 10, 1, 0.1, 10 −2 , 10 −3 , and 10 −4 m, respectively.The black lines show the case without gas drag.The dashed lines in the bottom panels show the Hill sphere of the planet.

Figure 2 .
Figure 2. Distribution of the initial orbital parameters (b, z s ) of captured particles (dust accretion band).The colors show the size of the particle; light blue, orange, green, blue corresponds to r d = 0.1, 10 −2 , 10 −3 , and 10 −4 m, respectively.The vertical axis is plotted at every z s = 0.1r H , but slightly shifted up and down to avoid overlapping of the lines.

Figure 3 .
Figure 3. Radial and vertical positions (R, z) where the particle capture conditions are met for the case of 0 < z s ≤ r H . Orange and blue points represent results for r d = 10 −2 m and 10 −4 m, respectively.The two black vertical dashed lines show the current radial positions of Io and Callisto.The results for r d = 10 −3 m are not shown here, but are similar to the case of r d = 10 −4 m.

Figure 4 .Figure 5 .
Figure 4. Capture rate P cap as a function of the initial altitude z s for r d = 10 −4 m.The colors show the difference of the planetary mass; green, black, red correspond to M p = 0.2M Jup , 0.4M Jup , and 1M Jup , respectively.

Figure 6 .
Figure6.Total capture rate P cap,total as a function of h d /h g .Top, middle, and bottom panels show the results for M p = 0.2, 0.4, and 1M Jup , respectively.Colors show the dust size, and orange, green, and blue correspond to r d = 10 −2 , 10 −3 , and 10 −4 m, respectively.

Figure 7 .
Figure7.Degree of dust retention of accreting gas into the CPD as a function of r H /h g .The top, middle, bottom panels show the case of r d = 10 −2 , 10 −3 , and 10 −4 m, respectively.Color means the dust scale height in the PPD; violet, green, blue, and red correspond to h d /h g = 0.1, 0.3, 0.6, and 1.0, respectively.

Figure 8 .
Figure8.Relationship between the depth of dust gap and r H /h g obtained by the empirical formulae ofZhang et al. (2018).Colors represent results for different values of the viscous parameter of the PPD α PPD : red, blue, and green correspond to α PPD = 10 −2 , 10 −3 , and 10 −4 , respectively.Top, middle, and bottom panels show the cases of r d = 10 −2 , 10 −3 , and 10 −4 m, respectively.

Figure 9 .
Figure9.Dust mass accretion rate onto the CPD as a function of the dust scale height h d , where the effects of dust gap is taken into account.We assumed f dust = 1 for simplicity, neglecting dust depletion in the PPD in the late stage of giant planet formation (see text for details).Colors represent the results for different values of the viscous parameter of the PPD α PPD : red, blue, and green correspond to α PPD = 10 −2 , 10 −3 , and 10 −4 , respectively.Top, middle, and bottom panels show the cases of r d = 10 −2 , 10 −3 , and 10 −4 m, respectively.