Delivery of gas onto the circumplanetary disk of giant planets: Planetary-mass dependence of the source region of accreting gas and mass accretion rate

Gas accretion onto the circumplanetary disks and the source region of accreting gas are important to reveal dust accretion that leads to satellite formation around giant planets. We performed local three-dimensional high-resolution hydrodynamic simulations of isothermal and inviscid gas flow around a planet to investigate planetary-mass dependence of gas accretion band width and gas accretion rate onto circumplanetary disks. We examined cases with various planetary masses corresponding to M_p=0.05-1M_{Jup} at 5.2 au, where M_{Jup} is the current Jovian mass. We found that the radial width of the gas accretion band is proportional to M_p^{1/6} for the low-mass regime with M_p<0.2 M_{Jup} while it is proportional to M_p for the high-mass regime with M_p>0.2M_{Jup}. We found that the ratio of the mass accretion rate onto the circumplanetary disk to that into the Hill sphere is about 0.4 regardless of planetary mass for the cases we examined. Combining our results with the gap model obtained from global hydrodynamic simulations, we derive semi-analytical formulae of mass accretion rate onto circumplanetary disks. We found that the mass dependence of our three-dimensional accretion rates is the same as the previously-obtained two-dimensional case, although the qualitative behavior of accretion flow onto the CPD is quite different between the two cases.


INTRODUCTION
Recent developments of numerical hydrodynamic simulations have made progress in understanding gas flow around protoplanets. When a protoplanet accretes gas from the protoplanetary disk (PPD), a circumplanetary disk (CPD) naturally forms around a sufficiently massive protoplanet while a spherically symmetric envelope forms around a small protoplanet (e.g., Tanigawa & Watanabe 2002;Machida et al. 2008;Machida 2009;Ayliffe & Bate 2009;Tanigawa et al. 2012;Gressel et al. 2013;Szulágyi et al. 2014Szulágyi et al. , 2016Szulágyi et al. , 2017Fung et al. 2019;Lambrechts et al. 2019;Schulik et al. 2019Schulik et al. , 2020. Whether envelopes or CPDs form depends not only on the planetary mass but also on the surface temperature of the planet and the equation of state (Szulágyi et al. 2016(Szulágyi et al. , 2017Fung et al. 2019). High-resolution three-dimensional hydrodynamic simulations demonstrate a meridional circulation of the gas near a gap opened by a giant planet in the PPD: gas flows into the gap at a high altitude over the CPD, and the planet pushes the gas back into the PPD while accreting part of the gas (Tanigawa et al. 2012;Szulágyi et al. 2014;Morbidelli et al. 2014;Fung & Chiang 2016). Similar circulating flows have also been reported in simulations around less massive planets corresponding to super-Earths and sub-Neptunes (Ormel et al. 2015;Fung et al. 2015;Cimerman et al. 2017;Kurokawa & Tanigawa 2018;Kuwahara et al. 2019).
Growing giant planets accumulate gas and solid particles via their CPD, and principal regular satellites, such as the Galilean satellites and Titan, are thought to be formed in the CPD. Thus, accretion of gas and solid particles onto the CPD is important for both planet and satellite formation. Kilometer-sized or larger planetesimals can be captured by gas drag from the CPD (Fujita et al. 2013). Some of the captured planetesimals can be ablated and become particles with the size of pebbles or dust, and the remaining part can become satellitesimals (Ronnet & Johansen 2020). However, capture efficiency of planetesimals significantly decreases when a gap in the planetesimal disk forms at the planet's orbit and/or the velocity dispersion of the planetesimals is too large (Fujita et al. 2013;Suetsugu et al. 2016;Suetsugu & Ohtsuki 2017), unless the planetesimal reservoir exterior to the gap is perturbed by other planetary embryos (Ronnet et al. 2018). On the other hand, dust particles accreted with gas are also potential building blocks of satellites and would become a major contributor in circumstances where planetesimal capture is inefficient (Canup & Ward 2002Suetsugu & Ohtsuki 2017;Shibaike et al. 2019;Batygin & Morbidelli 2020).
The motion of sufficiently small particles delivered from the PPD into the CPD is strongly influenced by the gas flow around the planet, such as the meridional circulation mentioned above. For example, using high-resolution local three-dimensional hydrodynamic simulations assuming isothermal and inviscid gas around a growing giant planet, Tanigawa et al. (2012) found that the gas vertically accretes onto the CPD but flows radially outward in the midplane, while the direction of the radial flow in the midplane of the CPD would depend on the viscosity of the gas ). Using orbital integration of particles that takes account of the gas flow obtained by Tanigawa et al. (2012), Tanigawa et al. (2014) examined accretion of particles initially confined within the midplane of the PPD. They found that mm-sized or smaller particles (i.e., those with the Stokes number significantly smaller than unity) cannot accrete into the CPD because of the strong influence of the outflowing gas in the midplane. On the other hand, Homma et al. (2020) performed similar calculations but took account of the vertical motions of the gas and particles. They found that sufficiently small particles coupled to the vertically accreting gas can accrete onto the CPD when they are vertically stirred up in the PPD. However, they also found that the dust-to-gas ratio in the gas accreting onto the CPD can become significantly smaller than that in the PPD if the turbulence in the PPD is weak (corresponding to the turbulence parameter α ∼ 10 −4 , Shakura & Sunyaev 1973). More recently, Bi et al. (2021) and Szulágyi et al. (2022) carried out global three-dimensional hydrodynamic simulations with dust particles. They found that the dust particles are stirred up vertically well above the midplane due to the meridional circulation, where the spiral wakes created by the planet bring up the dust from the midplane. Thus, it is crucial for understanding satellite formation around gas giants to study the delivery of gas and dust particles onto the CPD from the PPD using three-dimensional high-resolution hydrodynamic simulations, such as the one shown in Tanigawa et al. (2012). However, Tanigawa et al. (2012) (and also Tanigawa et al. 2014 andHomma et al. 2020) examined only the case of r H /h g = 1, where r H is the planet's Hill radius and h g is the scale height of the gas in the PPD, corresponding to a planetary-mass of 0.4M Jup at the current Jovian orbit for the minimum-mass solar nebula model (MMSN;Hayashi 1981). Since formation of principal regular satellites is likely to occur during planet formation, it is necessary to understand how the flow pattern and amount of gas and dust particles delivered onto the CPD depend on planetary mass. The dependence of gas and dust supply on planetary mass will also help us understand satellite formation in different planet systems, such as Jupiter and Saturn, or exoplanet systems.
In a series of papers, we will investigate planetary-mass dependence of delivery of gas and dust particles onto CPDs. As the first step of understanding planetary-mass dependence of dust supply onto CPDs, we focus on gas accreting onto CPDs in this paper. Accretion of dust particles that are not completely coupled to the gas will be investigated in our subsequent paper. Details of the structure of the CPD depend not only on the planetary mass but also on various other factors such as the equation of state, viscosity, and surface temperature of the planet (Machida 2009;Machida et al. 2010;Szulágyi et al. 2014Szulágyi et al. , 2016Szulágyi et al. , 2017Gressel et al. 2013;Fung et al. 2019). In contrast, the picture that the gas around a sufficiently massive planet vertically accretes onto the CPD is rather insensitive to numerical settings (Machida et al. 2010;Szulágyi et al. 2014Szulágyi et al. , 2017Gressel et al. 2013). Thus, we consider the simplest case of an isothermal and inviscid fluid, following Tanigawa et al. (2012). We perform three-dimensional high-resolution hydrodynamic simulations with various planetary masses and investigate the dependence of gas accretion on planetary mass.
The structure of this paper is as follows. In § 2, we will explain basic equations and numerical settings. In § 3, we will show the results of hydrodynamic simulations. In § 4, we will derive the semi-analytical formulae of mass accretion rate of gas by combining our results with the gap model obtained from global hydrodynamic simulations . We will also discuss comparisons with previous works in § 4. We will summarize this work in § 5.

Basic equations
The model used in this work is similar to those described in Machida et al. (2008) and Tanigawa et al. (2012). We consider a giant planet on a circular orbit embedded in a PPD. We erect a 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 define r and R as r ≡ x 2 + y 2 + z 2 and R ≡ x 2 + y 2 , respectively.
As in Tanigawa et al. (2012), we assume that the gas is inviscid and isothermal, and denote its density, velocity, and pressure as ρ, v, and P , respectively. We normalize time by the inverse of the Keplerian angular velocity Ω −1 K = (GM c /a 3 ) −1/2 , length by the scale height of the PPD h g , and density by Σ 0 /h g , where G, M c , a, and Σ 0 are the gravitational constant, mass of the central star, semi-major axis of the planet, and unperturbed surface density of the PPD, respectively. In this case, velocity is normalized by sound speed c s = h g Ω K . In the following, hats denote normalized quantities. Then, the basic equations can be written as (for details, see Machida et al. 2008) where e z = (0, 0, 1) is the unit vector in the z-direction, and In the above,Φ tidal is the tidal potential of the central star,Φ p is the gravitational potential of the planet, and a constant (9/2)r 2 H is added in Equation (4) so thatΦ = 0 at the two Lagrangian points L 1 and L 2 of the planet;r H is the normalized Hill radius of the planet denoted bŷ where M p is the mass of the planet. This parameterr H is the only one physical parameter in our formalization. Note thatr H = 1 corresponds to the planetary mass M p equal to the so-called thermal mass, 3(h g /a) 3 M c (e.g., Ginzburg & Chiang 2019) To determine the quantitative relationship betweenr H and M p , we need to give a semimajor axis dependence of the disk scale height. If we consider an optically thin disk and assume that the temperature distribution is T = 280 K (a/1 au) −1/2 (Hayashi 1981), the scale height can be derived as h g = c s Ω −1 K = 4.93 × 10 9 a 1 au 5/4 m.
Using Equations (7) and (8), we obtain On the other hand, recent MHD simulations show that the temperature in magnetically accreting disks is lower than that in viscous accretion disks (Mori et al. 2019). If we consider a colder disk and assume that the temperature distribution is T = 100 K (a/1 au) −3/7 (Kusaka et al. 1970), the scale height can be derived as and from Equations (7) and (10) (11) Figure 1 shows the relationship betweenr H and M p (i.e., Equations (9) and (11)) in the case of a = 5.2 au and a = 9.6 au, respectively, for the two disk models described above. To investigate the dependence of gas accretion on planetary mass, we ran 11 simulations varying the value ofr H and analyzed quasi-steady state for each run. We ran the simulations with 11 values ofr H in the range of r H = 0.5 − 1.36, corresponding to M p = 0.05M Jup − 1M Jup for the Hayashi disk (Equation (9)) and M p = 0.01M Jup − 0.2M Jup for the cold disk (Equation (11)) at the current Jovian orbit, respectively.
In the following, when we refer to M p assumed in each simulation, it is obtained from Equation (9) assuming a = 5.2 au (the current Jovian orbit) and M c = M (the solar mass), unless explicitly stated otherwise.   (9), which corresponds to the Hayashi model (T = 280 K (a/1 au) −1/2 ). The blue lines are drawn from Equation (11), which corresponds to the cold disk (Kusaka) model (T = 100 K (a/1 au) −3/7 ). The solid and dashed lines show the cases for a = 5.2 au (the Jovian orbit) and 9.6 au (the Saturnian orbit), respectively. The symbols correspond to the values ofr H we adopted in our simulations.

Numerical settings
We used a three-dimensional nested-grid hydrodynamic simulation code (e.g., Machida et al. 2005Machida et al. , 2006. Our numerical settings are basically the same as in Tanigawa et al. (2012). Assuming symmetry about theẑ = 0 plane, we set the sizes of the computational domain tox ∈ [−12, 12],ŷ ∈ [−12, 12], andẑ ∈ [0, 6], i.e., the length of each side of the simulation box is (L x ,L y ,L z ) = (24, 24, 6). We set the nested-grid composed of 11 levels at maximum (l max = 11), where l is the level of the nested-grid. As the level increases, the length of each side is halved keeping the planet at the origin, while the numbers of cells in each grid are kept constant as (n x , n y , n z ) = (64, 64, 16). The finest cell width isL x /64/2 11−1 = 3.66 × 10 −4 , which roughly corresponds to one-fourth of the physical size of a planet with a mean density of 1 g cm −3 .
Initial and boundary conditions are also the same as in Tanigawa et al. (2012). Initially, the gas density distribution is in hydrostatic equilibrium in the z-direction, and the Keplerian velocity field is assumed. We set the unperturbed boundary in the x-direction, the mixed condition of unperturbed and periodic boundaries in the y-direction, and the mirror condition in the z-direction, respectively (for details, see Tanigawa et al. 2012). We set the smoothing length r sm = 2.44 × 10 −4 h g to avoid singularity of the planet gravity, but no sink cells are set in the present work.
We set l = 1 at the start of each simulation and add deeper levels of grids (l > 1) gradually with time. We empirically determined to create a new level with l = 2 att = 20, l = 3 att = 50, and l = 4 att = 150, respectively, and deeper levels with l ≥ 5 are created when the gas flow in each level is sufficiently relaxed (see also Tanigawa et al. 2012). We observed that non-negligible oscillations appeared in the density and velocity fields for the case ofr H > 1.0 (i.e., M p > 0.4M Jup ), while quasistationary flow patterns can be obtained forr H ≤ 1.0. In order to analyze the field data, while we use snapshots at the end of each simulation for the case ofr H ≤ 1.0, we obtain time-averaged values of the density and velocity forr H > 1.0 as described below. For l ≤ 4, we take a time average for t ≥ 150 because numerical oscillations mainly occur in the region of l ≤ 4. We can reduce the effect of numerical oscillation by taking a time average of the snapshots over a sufficiently long duration compared to the period of the oscillation. For l ≥ 5, we use snapshots at the end of each simulation (t ∼ 170). The spatial region corresponding to l ≥ 5 is roughly within the Hill sphere and the gas distribution is nearly axisymmetric, thus discontinuity in the gas distribution and flow patterns between this region and the other region with l ≤ 4 is negligible. Thus, we combine time-averaging data for l ≤ 4 and snapshots for l ≥ 5. We confirmed that the discontinuities of combined gas fields do not significantly affect our qualitative results, while this method might affect quantitative results to some degree, such as gas accretion rates (Fung et al. 2019). , where a CPD with high density and an axisymmetric structure can be confirmed within the Hill sphere. The vertical inflow from high altitude to the CPD along the z-axis and a weak outflow at the midplane can be observed in the lower panels. . The upper and lower panels correspond to the planes where the cell closest to z = 0 and y = 0, respectively. The lemon-shaped region enclosed by the gray line in each panel is the planet's Hill sphere (its surface is defined byΦ = 0 in Equation (4)). Sizes of velocity vectors are appropriately chosen in each panel.  When the planetary mass is high enough, radial structure of the gas flow within the Hill radius is characterized by the Hill radius. Figure 4 shows radial distributions of radial velocity v R (top) and the z-component of specific angular momentum j z (bottom) around the planet at the midplane for the high mass cases (r H = 1.0 − 1.36; M p = 0.4M Jup − 1M Jup ). *1 Both radial velocity and specific angular momentum are azimuthally averaged around the planet and normalized by the values for the Keplerian rotation around the planet (v Kep,p and j Kep,p , respectively). Note that the specific angular momentum is measured in the rotational frame. We found that the radial distributions of the scaled radial velocity and the scaled specific angular momentum do not depend on planetary mass appreciably. We also found that gas elements at R < 0.2r H (left side of the vertical dashed line) rotate with a velocity over 80% of the Keplerian velocity around the planet. Defining this region as a CPD, the radius of the CPD is 0.2r H , i.e., proportional to the Hill radius of the planet forr H ≥ 1.  For the low mass cases, on the other hand, radial structure of the gas flow within the Hill radius is characterized by the Bondi radius r B = GM p /c 2 s . Figure 5 shows radial distributions of v R and j z for the low mass cases (r H = 0.5 − 0.8; M p = 0.05M Jup − 0.2M Jup ). We found that v R and j z shown as a function of the horizontal radial distance scaled by the Bondi radius do not depend on planetary mass appreciably. Gas elements at R < 0.07r B (left side of the vertical dashed line) rotate with a velocity over 80% of the Keplerian velocity around the planet. Thus, the radius of the CPD in this case is proportional to the Bondi radius forr H ≤ 0.8 (M p ≤ 0.2M Jup ).

Structure of gas flow within the Hill sphere
Our results that the radii of CPDs are proportional to the Bondi/Hill radius in the low-/high-mass regimes are consistent with the results for the isothermal and inviscid cases in Fung et al. (2019). The transition planetary mass between the two regimes we found is roughly consistent with the expectation of Fung et al. (2019) (r H 0.7; q thermal 4 in their notation).

Gas accretion onto circumplanetary disks
Accretion of small particles onto a CPD is largely influenced by the flow pattern of the gas around the planet. Previous studies show that the gas coming from a limited range of radial and vertical region in the PPD can accrete onto the CPD (e.g., Tanigawa et al. 2012, we call this region "accretion band", hereafter). The initial radial location of dust particles accreting onto a CPD is related to gas accretion band (Homma et al. 2020). In this section, we focus on the radial and vertical structure of the gas accretion band, which is important for accretion of gas and dust particles into the CPD. Note that the structure of the gas accretion band is hardly influenced by the gap depth, which we cannot precisely evaluate by our local simulation.
In Figures 6 and 7, starting radial and vertical positions of streamlines (x 0 and z 0 , respectively) at y = 12 are classified by color according to their final destination (Tanigawa et al. 2012;Homma et al. 2020). Those gas elements with large x 0 pass by the planet (passing; shown in red), while those with small x 0 move along horseshoe orbits (horseshoe; shown in green). Some of the streamlines between passing and horseshoe regions reach within r = 0.2r H (accreting; shown in dark blue) *2 , and others enter the planet's Hill sphere (whose surface is defined by the equipotential surface witĥ Φ = 0 in Equation (4)) and then escape from it (recycling; shown in light blue). We confirmed the three-dimensional structure of gas accretion, which was reported by previous works that examined cases for M p ∼ M Jup (Machida et al. 2008(Machida et al. , 2010Tanigawa et al. 2012;Gressel et al. 2013;Szulágyi *2 For the low-mass cases, we found that the gas streamlines which reach within r = 0.2r H also reaches within r = 0.07r B (< 0.2r H ) and the same accretion bands are obtained for both criteria of accretion. This means that the gas which reaches within r = 0.2r H loses its angular momentum at the shock surface on the CPD and flows toward the inner region of the CPD (see also Figure 15 in Tanigawa Figure 6 shows the streamlines with (a)ẑ 0 = 0.5 and (b)ẑ 0 = 0, respectively, forr H = 1.36 (M p = 1M Jup ). The off-midplane gas elements encounter the planetary shock and change their direction slightly upward. Then, they accrete onto the CPD (Figure 6a). On the other hand, the gas elements which originate from the midplane cannot enter the Hill sphere, because radially outward flow dominates at the midplane (Figure 6b; see also Figure 2 of Tanigawa et al. 2012). It should be noted that the path of the accreting streamline likely depends on temperature structure. According to global radiative-hydrodynamic simulations (Schulik et al. 2020), the spiral arm tilt, which depends on the vertical structure of temperature, controls the direction of accreting gas flow. The direction change observed in Figure 6a is less clear than the one found in Schulik et al. (2020). This is probably due to the difference in the temperature distribution. In our isothermal simulation, the vertical distribution of temperature does not have a gradient, which likely affects the tilt of the spiral arms. Figure 7 shows radial and vertical distributions of each type of streamline at the starting points (x 0 , z 0 ) for planetary masses ofr H = 0.8, 1.0, and 1.36 (M p = 0.2, 0.4, and 1M Jup ), respectively. In all these cases, we found that the accretion bands (the dark blue regions) do not exist near the midplane (z 0 0) but are distributed above the midplane continuously in the z-direction in a zonal pattern. We also found that the radial (x-direction) width of the accretion band increases as the planetary mass increases. This indicates that a more massive planet can acquire gas and dust particles from a radially wider region in the PPD. horseshoe (moving along a horseshoe orbit), dark-blue: accreting (reaching within r = 0.2r H ), light-blue: recycling (streamlines which enter the planet's Hill sphere and then escape from it). The lemon-shaped region enclosed by the dashed line in each panel is the planet's Hill sphere (its surface is defined byΦ = 0 in Equation (4)).
Next, we evaluate the radial width of the accretion band and its dependence on planetary mass. Let the non-dimensional radial width of accretion band at a scaled starting altitudeẑ 0 beŵ(ẑ 0 ). Assuming the hydrostatic equilibrium in the z-direction, the vertical density distribution of gas is given byρ We take an average ofŵ(ẑ 0 ) weighted by the above vertical density distribution aŝ where we assumedẑ 0,max = 3 because contribution of gas accretion fromẑ 0 > 3 is negligible. Figure 8 shows the dependence of the averaged width of the accretion band as a function ofr H , which depends on the planetary mass as shown in Equation (7). The crosses are obtained from our numerical results. The dashed line is the fitting function for the low-mass regime: and the solid line shows the fitting function for the high-mass regime: Thus, planetary-mass dependences of the averaged accretion band width are different between the low-mass and high-mass regimes;w ∝ r , respectively. Note that the planetary-mass dependence ofw for the low-mass regime has larger uncertainties due to the small number of cases we examined.  (14) and (15), respectively.

Mass accretion rate onto circumplanetary disks 3.4.1. Vertical structure of accreting gas
As shown in the previous section, the gas accreting into the CPD comes from high altitudes, and the vertical structure of the accretion band depends on the planetary mass. On the other hand, the gas density in the unperturbed PPD decreases with increasing altitude. In order to understand which vertical part of the PPD contributes most to the accretion into the CPD, we examine the mass accretion rate as a function of the initial altitudeẑ 0 of the accreting gas. Figure 9 shows the mass accretion rate per unit length in the z-direction of the gas originating from an altitudeẑ 0 (> 0) into the CPD defined aŝ wherex 0,min andx 0,max are the values ofx 0 at the inner and outer boundaries of the accretion band atẑ 0 , respectively. Equation (16) includes only the gas accreting from x 0 > 0 and z 0 > 0. We show the plots off CPD (ẑ 0 ) for three cases of different planetary masses corresponding tor H = 0.8, 1.0, and 1.36. We found that a major contribution to the mass accretion comes from lower altitudes for higher planetary masses; the peak of the mass accretion rate is at z 0 /h g 1.0, 0.6, and 0.3 forr H = 0.8, 1.0, and 1.36, respectively. This reflects the vertical structure of the accretion band (see Figure 7), which becomes wider at lower z 0 and is distributed down to lower altitudes with increasing planetary mass. We also show in Figure 9 the cumulative distribution of the mass accretion rate normalized by the total mass accretion rate with the dashed lines. We found that about 80% of gas accreting into the CPD comes from z 0 1.5h g , and the contribution from lower altitudes increases with increasing planetary mass, as we described above. For example, the contribution from z ≤ 0.3h g is 1.6%, 10%, and 20% forr H = 0.8, 1.0, and 1.36, respectively.

Mass accretion rate evaluated from our simulation
As shown in Figure 7, only part of the gas entering the planet's Hill sphere comes sufficiently close to the planet to accrete into the CPD, and the rest is expelled back to the PPD as part of the meridional circulation (e.g., Tanigawa et al. 2012;Morbidelli et al. 2014;Fung & Chiang 2016). Thus, the mass accretion rate of the gas into the CPD (Ṁ acc,CPD ) is generally smaller than that into the Hill sphere (Ṁ acc,Hill ), and the ratioṀ acc,CPD /Ṁ acc,Hill is an important quantity when considering the delivery of small dust particles coupled to the gas into the CPD as building blocks of satellites. Here we calculateṀ acc,CPD andṀ acc,Hill from the above analysis of the streamlines and using the following equations:Ṁ acc,CPD = 4 andṀ acc,Hill = 4 wheref Hill (ẑ 0 ) is the mass accretion rate into the Hill sphere per unit length in the z-direction, which is obtained by using both recycling and accreting streamlines in a similar calculation to Equation (16) (see Figures 6 and 7). The factor of four is multiplied in Equation (17) and Equation (18) to account for the symmetry about the x = 0 and z = 0 planes. We did not adopt the sink cell. Thus, in the quasi-steady state of the gas flow in our simulations, the accreted gas accumulates in the CPD while part of it is expelled back toward the protoplanetary disk as an outflow near the midplane of the CPD (the top panel of Figures 4 and 5). In order to precisely derive the amount of the gas that is accreted in the CPD but then expelled back toward the protoplanetary disk as an outflow from the midplane region of the CPD, we will need to investigate the long-term viscous evolution of the CPD. However, it is difficult to investigate it in our present simulation where the gas is assumed to be inviscid and is beyond the scope of the present work. Thus, in the present work, we derived the accretion rate onto the CPD from the analysis of the streamlines that reach within r = 0.2r H (i.e., the dark blue region in Figure 7). In Figure 10, we show the planetary mass dependence ofṀ acc,CPD andṀ acc,Hill (top panel), and the ratioṀ acc,CPD /Ṁ acc,Hill (bottom panel). We found thatṀ acc,Hill /Σ 0 h 2 g Ω K 0.4 whenr H = 1, which is consistent with Tanigawa et al. (2012, their Figure 14), who found thatṀ acc,Hill /Σ 0 h 2 g Ω K 0.2 for the mass accretion rate onto the z > 0 surface of the CPD forr H = 1. We also found thaṫ M acc,CPD /Ṁ acc,Hill 0.4 regardless of the planetary mass. This means that gas accretion rate onto the CPD can be estimated from that into the Hill sphere in global simulation at least in the planetary mass range corresponding tor H = 0.5−1.36, as long as the resolution of the simulation is sufficient to resolve the gas flow on the Hill-radius scale. Also, the fact that the ratioṀ acc,CPD /Ṁ acc,Hill is nearly independent of the planetary mass indicates that this ratio is insensitive to whether the gap opens or not, probably because it is determined by the distributions of the accretion band and recycling region (Figure 7), which is likely to be insensitive to the gap depth. It should be noted that absolute values of the mass accretion rates shown in Figure 10 depend on the gap depth at the planet's orbit, which depends on numerical settings such as simulation time and resolution (e.g., Kanagawa et al. 2017). Thus, in § 4.1, we will discuss the corrected values of mass accretion rate using the gap model obtained from global hydrodynamic simulations .

Semi-analytical formulae of mass accretion rate
Here, we consider a physical interpretation of the planetary-mass dependence of mass accretion rates onto the CPD and derive semi-analytical formulae by combining our results with the gap model obtained from global hydrodynamic simulations , because our local simulation cannot reproduce the gap depth accurately as we mentioned before.
Let us consider the gas flowing through the area of the accretion band at (x, y) = (x ab , L y /2), where x ab is a typical value of the x-coordinate of the accretion band. In the coordinate system rotating with the planet, the velocity of the gas passing the area of the accretion band is (3/2)x ab Ω K assuming the Keplerian shear flow at y = L y /2. Then, we can express the mass accretion rate onto the CPD asṀ where Σ g,ab is the surface density of the gas at (x, y) = (x ab , L y /2), and D is the "accretion area" of the protoplanetary disk per unit time (Tanigawa & Tanaka 2016). Tanigawa & Tanaka (2016) used the expression for D for the gas accretion onto a growing planet based on the isothermal and inviscid two-dimensional hydrodynamic simulation of Tanigawa & Watanabe (2002). *3 Here we assume that D for the mass accretion rate on the CPD can be expressed in terms of the averaged accretion band width obtained from our three-dimensional hydrodynamic simulation ( § 3.3) as The factor of two is multiplied in the r.h.s. of Equation (20) to account for the symmetry about x = 0 plane. For simplicity, we assume x ab = 2.2r H , regardless of the planetary mass (Figure 7). Substituting x ab = 2.2r H and the expression forw using Equation (14) or (15) In Figure 11, we plot the values of D calculated by Equations (21) and (22), and compare them with the results of our hydrodynamic simulation, where D is derived byṀ acc,CPD /Σ g,ab,sim with Σ g,ab,sim being the surface density at the accretion band obtained by our hydrodynamic simulation (Figure 13). For comparison, we also plot the values of D obtained from the two-dimentional simulation by Tanigawa & Watanabe (2002), which was used in Tanigawa & Tanaka (2016). We found that the results of our hydrodynamic simulation can be well reproduced by Equation (22) in the high-mass regime, while they somewhat deviate from those calculated by Equation (21) in the low-mass regime. This indicates that our model given by Equations (19) and (22) reproduces our numerical results well in the high-mass regime, but it is not satisfactory for the low-mass regime. The assumption that x ab is independent of M p may be too simple in the low-mass regime, and effects of different flow patterns for cases with r H < h g would also be important. In any case, it is difficult to construct an accurate model in the low-mass regime from our results, because in the present work we mainly focus on large-mass planets with a CPD in relation to satellite formation, and we ran only a few simulations for the low-mass regime. It should be noted that the planetary-mass dependence given by Equation (22) is consistent with the one obtained by the two-dimensional simulation (Tanigawa & Watanabe 2002). Although the qualitative behavior of accretion flow onto the CPD such as the strength of shocks is quite different between the two-and three-dimensional cases, the mass dependence of the accretion rates remains the same. Planetary-mass dependence of D obtained by our hydrodynamic simulation (= M acc,CPD /Σ g,ab,sim ; red crosses) and the one calculated from Equations (21) (dot-dashed line) and (22) (dashed line). The values of D obtained by two-dimensional simulation of Tanigawa & Watanabe (2002) are also plotted (dotted line). Following Tanigawa & Tanaka (2016), we assume a = 5.2 au.
In the following, we assume that our simulation results can be approximated by Equations (19) and (22) for the whole range of planetary-mass we examined, and derive mass accretion rate using the empirical formula for the gap surface density obtained from global simulations  given as Σ min where Σ min is the surface density at the bottom of the gap. We assume Σ g,ab = Σ min because x ab = 2.2r H is located near the bottom of the gap.
In Figure 12, we plotted semi-analytical formula obtained from Equation (19) and Equation (22) (red line), and our simulation results. Our simulation results derived by Equation (17) are shown with the blue open circles (raw values), and those corrected to reproduce the true gap depth by multiplying by Σ min /Σ g,ab,sim (corrected values) are shown with the green filled circles. We find that the semianalytical formula for the mass accretion rate onto the CPD we obtained can reproduce the numerical results well in the high-mass regime, although they somewhat underestimate the numerical results in the low-mass regime. In the high-mass regime, the planetary-mass dependence of our formula is consistent with the accretion rate onto the planet derived in Tanigawa & Tanaka (2016), but the absolute values of the accretion rates derived from our three-dimensional simulation are about 60% of those obtained in Tanigawa & Tanaka (2016). Ginzburg & Chiang (2019) argued that, after the onset of runaway gas accretion but when the planet's Hill radius is still smaller than the disk scale height (i.e.,r H < 1 in our simulations), the gas accretion rate is expected to match the Bondi accretion rate (∝ M 2 p ) rather than the empirical formula obtained from the results of hydrodynamic simulations (∝ M 4/3 p ; Tanigawa & Tanaka 2016) in the limit of low planetary mass (r B r H ), because the envelope boundary should be approximately spherically symmetric and the Keplerian shear velocity at r B is negligible in such a case. It should be noted that our results derived on the basis of the analysis of the accretion bands obtained from hydrodynamic simulations cannot be applied to such low-mass cases. This work (Eq.19,22) simulation, raw simulation, corrected Figure 12. Comparison of the mass accretion rates obtained by our hydrodynamic simulation (marks) and semi-analytical formulae (lines). The red line represents our semi-analytic formula (Equations (19) and (22)), and the dashed line is that derived by Tanigawa & Tanaka (2016) using the results of two-dimensional hydrodynamic simulation. The raw values of our simulation results derived by Equation (17)

Comparisons with previous works on gas accretion band
Some previous works using high resolution local simulation also reported the increase of radial width of the gas accretion band as the planet mass increases. Tanigawa & Watanabe (2002) performed two-dimensional hydrodynamic simulations and found that the radial width of the accretion band increases as the planet becomes more massive. Their detailed analysis of streamlines showed that energy dissipation at the spiral shock around the planet determines the radial width of the accretion band. The energy dissipation increases as the planet becomes more massive. However, our threedimensional simulations did not show such significant energy dissipation at the spiral shock around the planet. For the three-dimensional accretion, energy dissipation at the shock surface on the CPD is much stronger than that at the spiral arm of the planet (see Figures 10 and 11 in Tanigawa et al. 2012). The structure of the accreting flow in the three-dimensional simulation is qualitatively different from the two-dimensional one. Further studies are needed to explain the planetary-mass dependence of accretion band width in the three-dimensional case with a simple model.
Although the planetary-mass dependence of accretion band width in the three-dimensional case is also examined by Machida et al. (2010), their definition of the accretion band width is different from ours. They defined the width by the difference between x 0,min and x 0,max defined for the whole accretion area in the x 0 − z 0 plane, without measuring the radial width of the accretion band at each altitude. The z-dependence of the accretion band width is important when discussing dust accretion onto the CPD (e.g., Homma et al. 2020) because vertical distribution of dust particles varies with turbulence strength in the PPD and deviates from gas distribution unless turbulence in the PPD is sufficiently strong.

Effects of thermodynamics
In this study, we use the isothermal equation of state (EOS). In the late stage of planet formation, the gas surface density is expected to be low. Fluid in low-density regions mainly controls the accreting gas flow, and isothermal assumption should be valid in this region. We note that some regions near the midplane could be optically thick and not well-modeled by isothermal gas because the gap depth is moderately deep in our mass range (Figure 13 in Appendix). Fung et al. (2019) showed that the flow patterns are similar between isothermal and adiabatic cases when planetary mass is large enough (r H ≥ 1.0). This is because CPD dynamics are dictated by gravity, and the effect of EOS becomes small. Forr H < 1.0, on the other hand, the gravity of the planet is not so large, and EOS may affect the flow pattern around the planet. Machida et al. (2010) investigated thermal effects on mass accretion rate and found that mass accretion rates with the isothermal and adiabatic EOSs are almost the same when planetary mass is small (r H ≤ 1.0). However, its planetary-mass dependence is slightly different forr H ≥ 1.0, and the isothermal case showed stronger dependence on planetary mass. This trend seems different from Fung et al. (2019).
As we mentioned in § 3.3, global radiation hydrodynamic simulations of Schulik et al. (2020) showed that the vertical temperature gradient causes the spiral arm tilt, which changes the direction of accreting flow. This indicates that the vertical temperature gradient likely affects the gas accretion onto CPD. Thus, analysis of the structure of the gas accretion band onto CPD with more realistic conditions is desirable (e.g., Szulágyi et al. 2016Szulágyi et al. , 2017Fung et al. 2019).

Effects of Other Numerical Settings
Settings of a numerical domain can affect the circumplanetary structure and gas accretion. Using global simulation, Gressel et al. (2013) reported that the CPD can have a non-axisymmetry structure, the side facing the central star being thicker. In contrast, the CPD of our local simulation is roughly axisymmetry around the planet. The non-axisymmetry in the global case may lead to different structures of the gas accretion bands in regions interior and exterior to the planet's orbit. However, the basic characteristics of the structure of gas flow near the CPD (i.e., vertical accreting flow and outflow in the midplane) should be similar in both regions.
The radius of the CPD can depend on the effective viscosity (sum of the physical and the numerical viscosities). Szulágyi et al. (2014) showed that the larger effective viscosity tends to produce a larger CPD. In our inviscid simulation, we found that the CPD radius is about 0.2r H for the massive planet cases and 0.08r B for the small planet cases, but these values may become larger for a viscous fluid. However, the qualitative results such as the CPD radius being proportional to the planet's Hill radius and the Bondi radius for massive and small planets, respectively, are likely not to change.

CONCLUSIONS
We performed three-dimensional high-resolution hydrodynamic simulations to investigate the gas flow in the local region around the planet. The main purpose of the present work was to understand the planetary mass dependence of the accreting gas flow onto the CPD, because it is important for delivery of satellite building blocks to the CPD around planets with various masses and was not reported in the previous work (Tanigawa et al. 2012). We performed simulations forr H = r H /h g = 0.5 − 1.36, corresponding to the planetary masses of 0.05M Jup − 1M Jup at the current Jovian orbit for the MMSN model (Hayashi 1981). The main results are as follows: CPD such as the strength of shocks is quite different between the two cases. We also found that the mass dependence of the accretion rate is different between the two-and three-dimensional cases for the low-mass regime (M p 0.2M Jup ), but our semi-analytical model is too simple to reproduce the numerical results for the low-mass regime and further studies are needed to construct a more accurate model.
When dust particles are small and strongly coupled with the accreting gas, the mass accretion rate of dust particles should be strongly influenced by the gas. Since actual motion of dust particles is expected to decouple from the gas (e.g., Homma et al. 2020), orbital integration of particles is needed to study supply of dust onto CPD. We will report results of orbital calculation of dust particles in such cases in our subsequent paper.
In this work, we focused on planets whose mass is below the Jovian mass. However, more massive planets have been observed beyond the Solar system (e.g., Liu & Ji 2020). Understanding satellite formation around exoplanet systems is important to construct a more general theory of satellite formation. Recently, some CPD candidates have been found (Ginski et al. 2018;Isella et al. 2019;Tsukagoshi et al. 2019;Benisty et al. 2021) as well as exomoon candidates (e.g., Teachey & Kipping 2018;Kipping et al. 2022). Further studies of gas and dust accretion onto the CPD of planets with various masses will help us understand satellite formation in our Solar System as well as in exoplanet systems.

ACKNOWLEDGMENTS
We thank Hidekazu Tanaka for helpful discussion. We also thank the anonymous reviewer for the careful report that helped us improve the manuscript. T.T. sincerely appreciate Willy Kley for all of his energetic research activities, which had been motivating T.T. for decades. This work was supported by JSPS KAKENHI Grant numbers JP22J10202, JP18K11334, JP21H00043, JP22H01286, JP19K14787, JP15H02065, and JP20K04051. N.M. gives thanks for the support by the Kobe University Doctoral Student Fellowship and JSPS Research Fellowship for Young Scientists. This research used the computational resources of the 2020 and 2021 Koubo Kadai on Earth Simulator (NEC SX-ACE) at JAMSTEC, supercomputing resources at Cyberscience Center, Tohoku University, and the general-purpose PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.

A. GAS SURFACE DENSITY AT THE ACCRETION BAND
The planetary-mass dependence of the gas surface density at the accretion band, Σ g,ab,sim , is shown in Figure 13 with red crosses. Σ g,ab,sim /Σ 0 has values near unity in the low-mass limit, i.e., gap opening is little progressed in this regime. On the other hand, it has a planetary-mass dependence of Σ g,ab,sim ∝ M −1 p in the high-mass regime (r H 1.0). The deviation between our results and those obtained from the empirical formulae obtained by global simulation, Equations (23) and (24), (dotdashed curve) increases with increasing planetary mass. This deviation is caused by the use of local simulation in the present work.  Figure 13. Dependence of surface density at the accretion band Σ g,ab , i.e., surface density of gas at (x, y) = (2.2r H , L y /2), onr H = r H /h g . The dashed line shows the function proportional to M −1 p . The dotted-dashed curve shows surface density at the gap bottom derived by Kanagawa et al. (2015) for the case of α = 4 × 10 −3 .