The Growth of Protoplanets via the Accretion of Small Bodies in Disks Perturbed by the Planetary Gravity

We investigate the gas flow around a protoplanet with a circular orbit of radius $a$ in a protoplanetary disk around the host star with mass $M_{*}$ via hydrodynamic simulation. We assume a compressible, inviscid, non-isothermal, non-self-gravitating fluid. The basic equations of the hydrodynamic simulations are the equation of continuity, the Euler’s equation, and the energy conservation equation, given by

$$\frac{\partial\rho_{\mathrm{g}}}{\partial t}+\nabla\cdot(\rho_{\mathrm{g}}\mbox{\boldmath${v}$}_{\mathrm{g}})=0,$$

(1)

$$\left(\frac{\partial}{\partial t}+\mbox{\boldmath${v}$}_{\mathrm{g}}\cdot\nabla\right)\mbox{\boldmath${v}$}_{\mathrm{g}}=-\frac{\nabla p}{\rho_{\mathrm{g}}}+(\mbox{\boldmath${F}$}_{\mathrm{cor}}+\mbox{\boldmath${F}$}_{\mathrm{tid}}+\mbox{\boldmath${F}$}_{\mathrm{p}}),$$

(2)

	$\displaystyle\frac{\partial E}{\partial t}+\nabla\cdot[(E+p)\mbox{\boldmath${v}$}_{\mathrm{g}}]$	$\displaystyle=$	$\displaystyle\rho_{\mathrm{g}}\mbox{\boldmath${v}$}_{\mathrm{g}}\cdot(\mbox{\boldmath${F}$}_{\mathrm{cor}}+\mbox{\boldmath${F}$}_{\mathrm{tid}}+\mbox{\boldmath${F}$}_{\mathrm{p}})$		(3)
		$\displaystyle-$	$\displaystyle\frac{U(\rho_{\mathrm{g}},T)-U(\rho_{\mathrm{g}},T_{0})}{\beta/\Omega},$

where $\rho_{\mathrm{g}}$ is the gas density, $\mbox{\boldmath${v}$}_{\mathrm{g}}$ is the gas velocity, $p$ is the pressure, $\mbox{\boldmath${F}$}_{\mathrm{cor}}$ is the Coriolis force, $\mbox{\boldmath${F}$}_{\mathrm{tid}}$ is the tidal force, $\mbox{\boldmath${F}$}_{\mathrm{p}}$ is the gravitational force, $\Omega$ is the orbital frequency, $\beta$ is the dimensionless cooling timescale, $T$ and $T_{0}$ is the fluid and background temperatures, respectively, and $U$ and $E$ are, respectively, the internal and total energy densities, given by

$$U=\frac{p}{\gamma-1},$$

(4)

$$E=U+\frac{1}{2}\rho_{\mathrm{g}}v_{\mathrm{g}}^{2},$$

(5)

with the ratio of specific heat $\gamma=7/5$.

The forces in Eqs. (2) and (3) are given by $\mbox{\boldmath${F}$}_{\mathrm{cor}}=-2\Omega\mbox{\boldmath${e}$}_{z}\times\mbox{\boldmath${v}$}_{\mathrm{g}}$, $\mbox{\boldmath${F}$}_{\mathrm{tid}}=3x\Omega^{2}\mbox{\boldmath${e}$}_{x}-z\Omega^{2}\mbox{\boldmath${e}$}_{z}$, and

$$\mbox{\boldmath${F}$}_{\mathrm{p}}=\nabla\left(\frac{GM_{\mathrm{p}}}{\sqrt{r^{2}+r^{2}_{\mathrm{s}}}}\right)\left\{1-\mathrm{exp}\left[-\frac{1}{2}\left(\frac{t}{t_{\mathrm{inj}}}\right)^{2}\right]\right\},$$

(6)

where $\mbox{\boldmath${e}$}_{i}$ is the unit vector in the $i$-direction, $M_{\mathrm{p}}$ is the mass of the planet, $G$ is the gravitational constant, $r$ is the distance from the center of the planet, $r_{\mathrm{s}}$ is the softening length, and $t_{\mathrm{inj}}$ is the injection time of planetary gravity. To avoid the numerical effect, the gravity of the planet is gradually inserted using the injection time $t_{\mathrm{inj}}=0.5\Omega^{-1}$ (Ormel et al., 2015a). The planet have an atmosphere with radius $\sim R_{\rm B}$, where $R_{\rm B}\equiv GM_{\rm p}/c_{\rm s}^{2}$ is the Bondi radius and $c_{\rm s}$ is the isothermal sound speed. To resolve the atmospheric structure, we set $r_{\mathrm{s}}=0.1R_{\rm B}$ according to Kurokawa & Tanigawa (2018).

The last term in Eq. (3) is the cooling function according to the $\beta$ cooling model where the temperature $T$ relaxes toward $T_{0}$ in the timescale $\beta/\Omega$ (e.g., Gammie, 2001). In this study, we adopt the $\beta$ cooling model to obtain the density profile and flow around the planetary atmosphere according to Kurokawa & Tanigawa (2018), who showed the atmosphere is formed around the planet using the $\beta$ cooling model. We set $\beta=(R_{\mathrm{B}}\Omega/0.1c_{\rm s})^{2}$ according to the previous studies (Kurokawa & Tanigawa, 2018; Kuwahara & Kurokawa, 2020a, b).

The time, velocity, and length are considered to be normalized by the reciprocal of the Keplerian orbital frequency $\Omega^{-1}$, the isothermal sound speed $c_{\mathrm{s}}$, and the gas scale hight $H\equiv c_{\mathrm{s}}/\Omega$, respectively. The basic equations are then characterized by a dimensionless number, given by

$$m=\frac{R_{\mathrm{B}}}{H}=\frac{GM_{\mathrm{p}}\Omega}{c_{\mathrm{s}}^{3}}.$$

(7)

The Hill radius is written as a function of $m$:

$$R_{\mathrm{H}}=\left(\frac{m}{3}\right)^{1/3}H.$$

(8)

We assume a solar mass host star and a disk temperature profile $T=270(a/1\ \mathrm{au})^{-1/2}$(i.e., the minimum mass solar nebula model, Weidenschilling, 1977b; Hayashi et al., 1985), so that $M_{\mathrm{p}}$ is given by (Kurokawa & Tanigawa, 2018)

$$M_{\mathrm{p}}\simeq 12m\left(\frac{a}{1\ \mathrm{au}}\right)^{3/4}M_{\oplus}.$$

(9)

We investigate the cases with $m=0.1$, $0.05$, and $0.03$, which correspond to $M_{\rm p}/M_{\oplus}=1.2$ $(4.0)$, $0.6$ $(2.0)$, and $0.36$ $(1.2)$ at $1$ $(5)$ au, respectively.

In this study, we use the hydrodynamic simulation code Athena++ (White et al., 2016; Stone et al., 2020). We choose the HLLC algorithm for a Riemann solver. Our simulations are performed in a spherical-polar coordinate $(r,\theta,\phi)$ centered on the planet. We focus on the embedded, non-gap-opening regime and thus the local simulations are likely to be valid. In order to resolve the flow structure around the planet in detail, we adopt a logarithmic grid for the radial dimension. In the polar angle direction, the cell size is proportional to $(3\psi^{2}+1)$, where $\psi$ is the angle from the midplane. The cell spacing near the midplane is smaller than near the pole (i.e., resolution near the midplane is higher than near the pole; Kurokawa & Tanigawa, 2018). The numerical resolution is set to be $128\times 64\times 128$ in the $r$, $\theta$, and $\phi$ directions, respectively. The computational domain ranges from $0$ to $2\pi$ for $\phi$, from $0$ to $\pi$ for $\theta$, and from $r_{\mathrm{inn}}$ to $r_{\mathrm{out}}$ for $r$, where $r_{\mathrm{inn}}$ and $r_{\mathrm{out}}$ are the radii at the inner and outer boundaries, respectively.

We assume the initial and outer boundary values in hydrodynamic simulations are those of a Kepler disk without pressure gradient term, therefore given by

$$\rho_{\mathrm{g},0}=\rho_{0}\ \mathrm{exp}\left[-\frac{1}{2}\left(\frac{z}{H}\right)^{2}\right],$$

(10)

$$\mbox{\boldmath${v}$}_{\mathrm{g},0}=-\frac{3}{2}\Omega x\mbox{\boldmath${e}$}_{y},$$

(11)

where $\rho_{0}$ is the density at the midplane and $z$ is the distance from the midplane. We focus on the shear regime of pebble accretion, where the shear velocity is more dominant than the headwind of the gas to determine the accretion rate. The condition satisfied to be in the shear regime is given by (Ormel, 2017)

	$\displaystyle M_{\mathrm{p}}$	$\displaystyle\gg$	$\displaystyle\frac{1}{8}\frac{v_{\mathrm{hw}}^{3}}{G\Omega St_{0}}$		(12)
		$\displaystyle\simeq$	$\displaystyle 2.0\times 10^{-4}M_{\oplus}\frac{1}{St_{0}}$
		$\displaystyle\times$	$\displaystyle\left(\frac{a}{1\ \mathrm{au}}\right)^{3/2}\left(\frac{v_{\mathrm{hw}}}{50\ \mathrm{m\ s^{-1}}}\right)^{3},$

where $v_{\rm hw}$ is the headwind of the gas and $St_{0}$ is the dimensionless stopping time of a particle (see Eq. 25). We perform simulations for $M_{\mathrm{p}}=0.36-1.2\ M_{\oplus}$ at $a=1\ {\rm au}$ ($M_{\rm p}=1.2-4.0\ M_{\oplus}$ at $a=5\ {\rm au}$) and $St_{0}=3.0\times 10^{-3}-1.0\times 10^{13}$, so that we ignore the headwind in our simulations.

We set the radius of the inner boundary according to that of the planet (Kuwahara & Kurokawa, 2020a). We assume the density of the planet $\rho_{\mathrm{pl}}=5\ \mathrm{g}/\mathrm{cm}^{3}$, so that the radius of the inner boundary is given by

	$\displaystyle r_{\mathrm{inn}}$	$\displaystyle=$	$\displaystyle\left(\frac{3M_{\mathrm{p}}}{4\pi\rho_{\mathrm{pl}}}\right)^{1/3}=\left(\frac{9M_{*}}{4\pi\rho_{\mathrm{pl}}}\right)^{1/3}\frac{R_{\rm H}}{a}$		(13)
		$\displaystyle\simeq$	$\displaystyle 3\times 10^{-3}m^{1/3}$
		$\displaystyle\times$	$\displaystyle\left(\frac{\rho_{\mathrm{pl}}}{5\ \mathrm{g\ cm^{-3}}}\right)^{-1/3}\left(\frac{M_{*}}{1M_{\odot}}\right)^{1/3}\left(\frac{a}{1\ \mathrm{au}}\right)^{-1}.$

We introduce the reflective boundary condition for the inner boundary, but Kurokawa & Tanigawa (2018) reported this condition generates the unphysical energy flow in the results of simulations by Athena++. In order to prevent this unphysical affair from affecting the entire flow, we introduce the artificial cooling at the three inner cells ($\beta=10^{-5}$; Kurokawa & Tanigawa, 2018). The boundary condition in the azimuthal direction is set to the periodic boundary, which means $A(r,\ \theta,\ \phi)=A(r,\ \theta,\ \phi+2\pi)$ holds for an arbitrary physical quantity $A$. We calculate until the steady state flow is almost achieved at $t=t_{\rm end}$. We set $t_{\rm end}$ according to Kuwahara & Kurokawa (2020a). A summary of our simulation parameters is showed in Table 1, which are chosen from the values at $a=1\ {\rm au}$.

We calculate orbits of particles in the gas flow obtained from hydrodynamic simulations. The equation of motion of the particle is given by

	$\displaystyle\frac{d\mbox{\boldmath${v}$}}{dt}=\left(\begin{array}[]{ccc}2v_{y}\Omega+3x\Omega^{2}\\ -2v_{x}\Omega\\ -z\Omega^{2}\end{array}\right)-\frac{\mathrm{G}M_{\mathrm{p}}}{r^{3}}\left(\begin{array}[]{ccc}x\\ y\\ z\end{array}\right)+\frac{\mbox{\boldmath${F}$}_{\mathrm{drag}}}{m_{\mathrm{p}}},$		(20)
			(21)

where $\mbox{\boldmath${v}$}=(v_{x},v_{y},v_{z})$ is the velocity vector of the particle, $\mbox{\boldmath${r}$}=(x,y,z)$ is its position vector with respect to the center of the planet, and $m_{\mathrm{p}}$ is the mass of the particle. The first and second terms on the right-hand side of Eq. (21) are the Coriolis and tidal forces and the gravity force of the planet, respectively. The third term is the gas drag force given by (e.g., Adachi et al., 1976)

$$\mbox{\boldmath${F}$}_{\mathrm{drag}}=-\frac{C_{\mathrm{D}}}{2}\pi r_{\mathrm{p}}^{2}\rho_{\mathrm{g}}u\mbox{\boldmath${u}$},$$

(22)

where $r_{\mathrm{p}}$ is a particle radius, $\mbox{\boldmath${u}$}\equiv\mbox{\boldmath${v}$}-\mbox{\boldmath${v}$}_{\mathrm{g}},\ u=|\mbox{\boldmath${u}$}|$ is the velocity of a particle relative to the gas, and $C_{\mathrm{D}}$ is the gas drag coefficient. We consider centimeter to kilometer sized particles; $C_{\mathrm{D}}$ is given by the Stokes gas drag for small particles, while $C_{\mathrm{D}}$ is constant for large particles (Adachi et al., 1976). In addition, if particle velocities exceed the sound speed due to planetary gravity, the supersonic gas drag law should be applied. The gas drag coefficient, $C_{\mathrm{D}}$ is approximated to be (Adachi et al., 1976; Tanigawa et al., 2014).

$$C_{\mathrm{D}}=\frac{12\nu}{r_{\rm p}u}+\frac{(2-w)M}{1.6+M}+w,$$

(23)

where $\nu$ is the kinetic viscosity, $M\equiv u/c_{\mathrm{s}}$ is the Mach number, and $w=0.4$ is the correction factor. In Eq. (23), the first, second, and third terms indicate Stokes, quadratic, and supersonic drag coefficients, respectively, and we ignore Epstein drag for simplicity. In our simulation, we do not consider the evaporation or ablation. The kinetic viscosity is given by $\nu=l_{\rm mfp}c_{\rm s}/2$, where $l_{\rm mfp}=\mu m_{\rm H}/\rho_{\rm g}\sigma$ is the mean free path of the gas, $\mu=2.34$ is the mean molecular weight, $m_{\rm H}=1.67\times 10^{-24}\ {\rm g}$ is the mass of the proton, and $\sigma=2\times 10^{-15}\ {\rm cm}^{2}$ is the molecular collision cross section (Chapman & Cowling, 1970; Adachi et al., 1976). The stopping time of a particle is expressed by

$$t_{\mathrm{stop}}=\frac{m_{\mathrm{p}}u}{|\mbox{\boldmath${F}$}_{\mathrm{drag}}|}=\frac{4\rho_{\mathrm{p}}r_{\mathrm{p}}^{2}}{9\rho_{\mathrm{g}}c_{\mathrm{s}}l_{\mathrm{mfp}}}\left(1+\frac{2M+1.6w}{1.6+M}\cdot\frac{r_{\mathrm{p}}u}{6c_{\mathrm{s}}l_{\mathrm{mfp}}}\right)^{-1},$$

(24)

where $\rho_{\mathrm{p}}$ is the density of a particle. We investigate the orbits of particles with different radii. We introduce the Stokes parameter, $St$ defined as the stopping time multiplied by $\Omega$. The initial particle has $u\ll c_{\rm s}$ so that the gas drag is mainly given in the Stokes gas drag regime. This initial Stokes parameter is approximated to be

$$St_{0}=4.5\times 10^{-4}\left(\frac{\rho_{\mathrm{p}}}{1\ \mathrm{g}\ \mathrm{cm}^{-3}}\right)\left(\frac{c_{\mathrm{s}}}{1\ \mathrm{km\ s^{-1}}}\right)^{-1}\left(\frac{r_{\mathrm{p}}}{1\ \mathrm{cm}}\right)^{2},$$

(25)

where the value of $c_{\mathrm{s}}$ is chosen as that approximately at $1\ \mathrm{au}$ in the minimum mass solar nebula model. Note that $St$ continuously changes during the calculation. The value of $St$ can be less than $St_{0}$ due to the gas drag for $u\gg c_{\rm s}$. We use $St_{0}$ given in Eq. (25) as a parameter instead of the radii of particles.

A particle is launched from a starting point ($x_{\mathrm{s}},y_{\mathrm{s}},z_{\mathrm{s}}$), where $y_{\mathrm{s}}$ is fixed at $y_{\mathrm{s}}=40R_{\mathrm{H}}$, and $x_{\mathrm{s}}$ and $z_{\mathrm{s}}$ are varied (Ida & Nakazawa, 1989; Ormel & Klahr, 2010). The particle is initially well coupled to the gas for $St_{0}\ll 1$, while the motion is determined by Kepler’s laws for $St_{0}\gg 1$. We thus set the initial condition according to the gas motion for $St_{0}<1$ or the orbital elements for $St_{0}>1$. If the orbital eccentricities of particles are much smaller than the Hill radius of the planet divided by $a$, the collision rate is independent of eccentricities (Ida & Nakazawa, 1989). We set initially circular orbits. The initial conditions are given as follows.

	$\displaystyle z_{\rm s}$	$\displaystyle=$	$\displaystyle\left\{\begin{array}[]{cc}z_{0}\ \ \ (St_{0}<1),\\ ia_{\rm s}\ \mathrm{sin}\ \omega\equiv z_{0}\ \mathrm{sin}\ \omega\ \ \ (St_{0}>1),\\ \end{array}\right.$
	$\displaystyle\mbox{\boldmath${v}$}_{0}$	$\displaystyle=$	$\displaystyle\left(0,-\frac{3}{2}\Omega x_{\mathrm{s}},v_{z,0}\right),$		(28)
	$\displaystyle v_{z,0}$	$\displaystyle=$	$\displaystyle\left\{\begin{array}[]{cc}0\ \ \ (St_{0}<1),\\ z_{0}\ \mathrm{cos}\ \omega\ \ \ (St_{0}>1),\\ \end{array}\right.$

where $i$ is the initial inclination of the particle, $a_{\rm s}$ is its initial semimajor axis and $\omega$ is its longitude of ascending node. We assume $\omega$ is distributed uniformly in the range between 0 and $\pi$, and take the average for the calculation of the collision rate. For $St_{0}<1$, we ignore the $z$-component of the tidal force in Eq. (21) assuming the balance with turbulent stirring according to Kuwahara & Kurokawa (2020a, b).

The lower limits of $x_{\mathrm{s}}$ and $z_{\mathrm{0}}$ are set to $x_{\mathrm{s}}=0.0001H$ and $z_{\mathrm{0}}=0$ and the spacial intervals of $x_{\mathrm{s}}$ and $z_{\mathrm{0}}$ are $0.0001H$ and $0.02H$. The upper limits of $x_{\mathrm{s}}$ and $z_{\mathrm{0}}$ are set to a few disk scale hights. We calculate orbits only coming from positive $x$ and $y$ because of the symmetry.

We use $\rho_{\mathrm{g}}$, $\mbox{\boldmath${v}$}_{\mathrm{g}}$, and $c_{\mathrm{s}}$ given from the hydrodynamic simulation at $t=t_{\mathrm{end}}$, at which the fluid is in a quasi-steady state. The starting point of the particle is out of the computational domain of our hydrodynamic simulation ($r_{\mathrm{out}}$), so that we set the velocity and density of gas in $r>r_{\mathrm{out}}$ are the same as the outer boundary conditions given in Eqs. (10) and (11). In our hydrodynamic simulations for $m=0.03$ and $m=0.05$, unexpected vortices are appeared in the horseshoe structure. These vortices are found by Kuwahara & Kurokawa (2020a, b), which are caused by the low resolution of the distant place from the planet. In these cases, we only use the results of the hydrodynamic simulation in $r<0.3H\ (m=0.03)$ and in $r<0.55H\ (m=0.05)$, respectively. All physical quantities obtained via hydrodynamic simulations are the discrete data, so that we interpolate physical quantities using the linear interpolation method (Kuwahara & Kurokawa, 2020a).

Orbital integration is terminated if any one of the following conditions is satisfied. (i) A particle goes far away; $|y|>40R_{\mathrm{H}}$. (ii) A particle collides with the planet; $r<r_{\mathrm{inn}}$, where $r_{\mathrm{inn}}$ is the inner boundary for hydrodynamic simulations and the physical radius of the planet in Eq. (13).

We numerically integrate Eq. (21) using the Runge-Kutta-Fehlberg scheme (Fehlberg, 1969; Eshagh, 2005; Ormel & Klahr, 2010). This integration method controls each time step comparing a fourth-order solution with a fifth-order solution. We set the relative error tolerance of $10^{-8}$, which ensures numerical convergence (Ormel & Klahr, 2010; Visser & Ormel, 2016).

We show the sketch of the orbital calculation of particles in Fig. 1.

First of all, we perform two-dimensional orbital calculations in which particles have orbits in the midplane ($z=0$). Figure 5 shows the trajectories of the particles with different $St_{0}$.

For the large Stokes parameter, $St_{0}=1.0\times 10^{13}$ (Fig. 5a), the trajectories are almost the same as those in the gas-free case because particles are too large to be affected by the gas drag (Ida & Nakazawa, 1989). Particles collide with the planet for the initial positions ($x_{\rm s}$) in three discrete bands (Fig. 5a).

For $St_{0}=1.0\times 10^{2}$ (Fig. 5b), the orbits outside the Hill sphere is almost the same as those for $St_{0}=1.0\times 10^{13}$. However, the particles feel strong gas drag in the atmosphere. All the particles entering the Bondi sphere collide with the planet.

For $St_{0}\gg 1$, the orbits of particles are almost independent of $St_{0}$ and similar to the gas-free ones unless $r<R_{\mathrm{H}}$. In the Hill sphere, the velocities of particles can be larger than the sound velocity and the gas density increases (see Fig. 4). Therefore, gas drag effectively damps the kinetic energies of particles, which induces the capture of passing particles.

For $St_{0}=1.0$ (Fig. 5c), the orbits of particles are similar to those for $St_{0}=1.0\times 10^{2}$. The orbits of particles in $y<0$ are closer to the streamlines than those for $St_{0}=10^{2}$ are. In addition, all particles entering the Hill sphere accrete onto the planet.

For $St_{0}=1.0\times 10^{-2}$ (Fig. 5d), the orbits are the almost same as the steady streamlines in Fig. 2. The horseshoe width is much smaller than that for $St_{0}\gtrsim 1$. Particles are prevented from accreting onto the planet by the horseshoe flow and Keplerian shear flow, and the particles coming from the narrow band between the horseshoe and Keplerian shear flows are allowed to collide with the planet. This result is consistent with previous studies (Popovas et al., 2018; Kuwahara & Kurokawa, 2020a, b; Homma et al., 2020).

We calculate the two-dimensional specific collision rate, $P_{\mathrm{col}}=P_{\mathrm{col,2D}}$, based on the orbital calculations. The definition of $P_{\mathrm{col,2D}}$ is given by

$$P_{\mathrm{col},2D}=2\int_{0}^{\infty}\Phi(x_{\mathrm{s}})|v_{0,y}|dx_{\mathrm{s}}=3\Omega\int_{0}^{\infty}\Phi(x_{\mathrm{s}})x_{\mathrm{s}}dx_{\mathrm{s}},$$

(31)

where $\Phi(x_{\mathrm{s}})=1$ if the particle collides with the planet and zero otherwise, and $v_{0,y}$ is $y$-component of the initial velocity in Eq. (28). In order to account for the collision from both positive $x_{\rm s}$ and negative $x_{\rm s}$, we multiply the integral over positive $x_{\rm s}$ by two.

Figure 6 shows the collision rate for $m=0.1$ as a function of $St_{0}$. For $St_{0}\gtrsim 10^{12}$, the collision rate converges to a constant value, which is the same as the gas-free limit (Ida & Nakazawa, 1989; Inaba et al., 2001).

For $St_{0}\gtrsim 1$, the collision rate increases with decreasing $St_{0}$. For large $St_{0}$, particles are scattered away due to close encounters with the planet, so that $P_{\mathrm{col,2D}}$ is small. However for small $St_{0}$, particles feel strong gas drag during close encounters because of high atmospheric density. Our hydrodynamic simulations indicate that density enhancement occurs even at $r\gtrsim R_{\rm B}$ if $r<R_{\rm H}$ (see Fig. 4). For $1\lesssim St_{0}\lesssim 10^{2}$, particles can be captured at the outer edge of the atmosphere, $r\sim R_{\rm H}$. In addition, a close encounter with the planet accelerates the velocity, which can exceed the sound velocity. The super-sonic gas drag effectively reduces the kinetic energy prior to scattering. Particles are then captured by atmosphere, which enhances $P_{\mathrm{col,2D}}$ (Inaba & Ikoma, 2003).

For $St_{0}\lesssim 1$, $P_{\mathrm{col,2D}}$ decreases with decreasing $St_{0}$. Small particles are well coupled to the gas. The flow pattern of the gas due to horseshoe and Keplerian shear limits the accretion of particles onto the planet (see Fig. 5d). Therefore, the collision rate sharply decreases with decreasing $St_{0}$.

First, we show the case of large particles (for $St_{0}>1$). For $ia_{\rm s}$ or $z_{0}$ $\ll R_{\rm H}$, the 3D orbits projected on the $x$-$y$ plane are almost the same as the 2D orbits. Figure 7a shows the orbits of particles projected on the $y$-$z$ plane for $St_{0}=1.0\times 10^{3}$. For $i\gg R_{\rm H}/a_{\rm s}$, particles pass through the planet. Particles passing near the planet can be accreted onto the planet.

Figure 7b shows the orbits of the smaller particles projected on the $y$-$z$ plane for $St_{0}=1.0\times 10^{-2}$. Particles can be accreted if they pass near the planet as well as the case of $St_{0}\gtrsim 1$.

We calculate the three-dimensional specific collision rate. We then need the vertical distribution of particles for the calculation.

For $St_{0}>1$, we set the Rayliegh type distribution for inclinations (Ida & Makino, 1992). Assuming the random orbital phases, we obtain the collision rate as,

	$\displaystyle P_{\rm col,3D}(i^{*})=\frac{6\Omega}{\pi i^{*2}}\int_{0}^{\infty}\int_{0}^{\pi}\int_{0}^{z_{0,\mathrm{max}}/a_{\rm s}}\Phi(x_{\rm s},i,\omega)x_{\rm s}i$
	$\displaystyle\times\mathrm{exp}\left(-\frac{i^{2}}{i^{*2}}\right)dx_{\rm s}d\omega di$
	,		(32)

where $i^{*}$ is the dispersion of $i$ and $z_{0,\mathrm{max}}$ is the upper limit of $z_{0}$. Note that the vertical distribution of particles under this assumption corresponds to a Gaussian distribution for $z_{0}$ with the scale hight of $a_{\rm s}i^{*}/\sqrt{2}$ (see Appendix B).

For $St_{0}<1$, we assume the Gaussian distribution function of $z_{0}$. The collision rate is given by

	$\displaystyle P_{\rm col,3D}(H_{\rm d})=\frac{3\Omega}{\sqrt{2\pi}H_{\rm d}}\int_{0}^{\infty}\int_{0}^{\infty}\Phi(x_{\rm s},z_{0})$
	$\displaystyle\times x_{\rm s}\mathrm{exp}\left(-\frac{z_{0}^{2}}{2H_{\rm d}^{2}}\right)dx_{\rm s}dz_{0},$		(33)

where $H_{\rm d}$ is the scale hight of particles.

Figure 8 shows $P_{\rm col,3D}(i^{*})$ for $St_{0}>1.0$ as a function of $H_{\mathrm{d}}=a_{\rm s}i^{*}/\sqrt{2}$, using the relation between $H_{\rm d}$ and $i^{*}$ given in Eq. (B6). For $H_{\rm d}\ll R_{\mathrm{H}}$, $P_{\rm col,3D}$ is independent of $H_{\rm d}$ and almost the same as $P_{\rm col,2D}$. For $H_{\rm d}\gg R_{\mathrm{H}}$, $P_{\rm col,3D}$ decreases with $H_{\rm d}$. This is caused by the passing by without collisions for high inclination particles (see Fig. 7a).

Figure 9 shows $P_{\rm col,3D}(H_{\rm d})$ for $St_{0}<1.0$. The dependence of $P_{\rm col,3D}$ on $H_{\rm d}$ is similar to that for $St_{0}>1$. For $St_{0}=0.003$, $P_{\rm col,3D}(H_{\rm d})$ increases around the $H_{\rm d}=R_{\rm H}$. This increase is caused by the vertical gas flow that enters from high latitudes as seen in Fig. 3 (see also Kuwahara & Kurokawa, 2020a).

First, we confirm the analytic formula in the gas free limit that means particles and $St_{0}$ are extremely large. The two-dimensional collision rate for the gas-free limit was derived in the previous studies, given by (Ida & Nakazawa, 1989; Inaba et al., 2001)

$$\frac{P_{\mathrm{col},\mathrm{free}0}(R_{\rm pl})}{R_{\mathrm{H}}^{2}\Omega}=11.3\sqrt{R_{\mathrm{pl}}/R_{\mathrm{H}}},$$

(34)

where $R_{\mathrm{pl}}$ is the radius of the planet. This formula is in agreement with our simulations for the huge Stokes number. Note that Eq. (34) is valid for $R_{\rm pl}\ll R_{\rm H}$. We improve the formula for $R_{\rm pl}\sim R_{\rm H}$ in §4.1.2.

Second, we derive the analytic formula for the atmospheric capture regime. Particles entering planetary atmosphere have high velocities due to planetary gravity. We estimate the terminal velocity of particles at $r=R_{\mathrm{B}}$ as

	$\displaystyle u$	$\displaystyle\sim$	$\displaystyle\frac{GM_{\mathrm{p}}}{R_{\rm B}^{2}}\frac{St_{0}}{\Omega}$		(35)
		$\displaystyle=$	$\displaystyle\frac{St_{0}}{m}c_{\mathrm{s}}.$

For $St_{0}\gtrsim m$, the velocity can be greater than the sound velocity, so that the super-sonic gas drag given by the second term in Eq. (23) is effective in the atmosphere.

The kinetic energies of particles are damped by gas drag in the atmosphere. Once the particles’ energies are smaller than the gravitational energy at the Hill radius, the orbits of the particles are bound by the planet. The kinetic energy at infinity is negligible. The condition for capture of particles is then given by

$$\Delta E<-\frac{GM_{\mathrm{p}}m_{\mathrm{p}}}{R_{\mathrm{H}}},$$

(36)

where $\Delta E$ is the energy loss due to gas drag in the atmosphere. We estimate $\Delta E$ from the integral of the energy loss along the particle orbit in the atmosphere with radius $R_{\rm atm}$, given by

	$\displaystyle\Delta E(q)$	$\displaystyle\approx$	$\displaystyle-\zeta\int\pi r_{\mathrm{p}}^{2}\rho_{\mathrm{atm}}(r)u^{3}\ dt$
		$\displaystyle=$	$\displaystyle-4\zeta\pi GM_{\mathrm{p}}r_{\mathrm{p}}^{2}\int^{R_{\mathrm{atm}}}_{q}\frac{\rho_{\mathrm{atm}}(r)}{\sqrt{r(r-q)}}\ dr,$

where $\zeta$ is the number of close encounters between the planet and the particle and we assume a parabolic orbit with pericenter distance $q$ for the particle and then we use the relations $u^{2}=2GM_{\rm p}/r$ and $dr/dt=\sqrt{2GM_{\rm p}(r-q)}/r$ for the derivation. We multiply two because the integral range is the half of a parabola. The number of close encounters and the atmospheric radius are determined according to our simulation and we adopt $\zeta=2$ and $R_{\rm atm}=0.5\ R_{\rm H}$. The following condition gives the captured radius $R_{\mathrm{cap}}$

$$\Delta E(R_{\mathrm{cap}})=-\frac{GM_{\mathrm{p}}m_{\mathrm{p}}}{R_{\mathrm{H}}}.$$

(38)

We derive $R_{\mathrm{cap}}$ satisfing Eq. (38).

According to Inaba & Ikoma (2003), we obtain $P_{\rm col}$ from $P_{\rm col,free0}$ using $R_{\rm cap}$ instead of $R_{\rm pl}$. However, $P_{\rm col,free0}$ is valid for $R_{\rm pl}\ll R_{\rm H}$. Figure 10 shows the collision rate in the gas-free case as a function of the planetary radius. For $R_{\rm pl}>0.05R_{\rm H}$, we obviously need to modify the formula. Therefore, we give the collision rate

$$P_{\mathrm{col},\mathrm{free}}=\mathrm{MIN}\left[\mathrm{MAX}\left(P_{\rm col,free0},\ P_{\rm col,free1}\right),\ P_{\rm col,free2}\right],$$

(39)

where

	$\displaystyle\frac{P_{\rm col,free1}(R_{\rm pl})}{R_{\rm H}^{2}\Omega}$	$\displaystyle=$	$\displaystyle 1.06\times 10^{3}(R_{\rm pl}/R_{\rm H})^{2},$		(40)
	$\displaystyle\frac{P_{\rm col,free2}(R_{\rm pl})}{R_{\rm H}^{2}\Omega}$	$\displaystyle=$	$\displaystyle 4.66(R_{\rm pl}/R_{\rm H})^{0.15}.$		(41)

For the atmospheric collision rate, we use $R_{\rm cap}$ instead of $R_{\rm pl}$ and then obtain

$$P_{\rm col,atm}=P_{\rm col,free}(R_{\rm cap}).$$

(42)

Third, we derive the analytic formula for the settling regime. There are two regimes for settling; the cases of the supersonic gas drag and the Stokes gas drag.

For $St_{0}\lesssim 1$, Ormel & Klahr (2010) derived $P_{\rm col}$ for constant $St$ in the Keplerian shear flow. They obtained the analytic formula by comparing the encounter time with the settling time. For a particle with an impact parameter $x_{\rm s}$, the encounter time is estimated to be $t_{\mathrm{enc}}=x_{\rm s}/|v_{0,y}|$, where $|v_{0,y}|=3\Omega x_{\rm s}/2$ is the encounter velocity. The settling time needed for particles to settle down to the planet is given by $t_{\mathrm{set}}=x_{\rm s}/u_{\mathrm{term}}$, where $u_{\mathrm{term}}$ is the terminal velocity determined by the force balance between the gravitational force and gas drag.

As explained above, for $St_{0}\gtrsim m$, the velocity of the particle approaching the planet exceeds the sound velocity, so that the supersonic gas drag, the second term in Eq. (23), is effective. The equation of force balance between the gravity and gas drag is given by

$$\frac{GM_{\mathrm{p}}m_{\mathrm{p}}}{r^{2}}=\pi r_{\mathrm{p}}^{2}\rho_{\mathrm{g}}u_{\mathrm{term}}^{2},$$

(43)

so that the terminal velocity is expressed by

$$u_{\mathrm{term}}=\sqrt{\frac{4GM_{\mathrm{p}}\rho_{\mathrm{p}}}{3\rho_{\mathrm{g}}x_{\rm s}^{2}}}\left(\frac{9\rho_{\mathrm{g}}c_{\mathrm{s}}l_{\mathrm{mfp}}}{4\rho_{\mathrm{p}}}\frac{St_{0}}{\Omega}\right)^{1/4}.$$

(44)

For $t_{\rm enc}\gtrsim t_{\rm set}$, particles are accreted onto the planets. The accretion condition is given as

$$C_{1}\frac{x_{\rm s}}{|v_{0,y}|}>\frac{x_{\rm s}}{u_{\mathrm{term}}},$$

(45)

where $C_{1}$ is the constant value on the order of unity. The impact parameter $x_{\rm s}$ required for the settling is given by

$$x_{\rm s}<x_{\rm ss}=\left(\frac{8}{3}\right)^{1/4}\sqrt{C_{1}}\left(\frac{\rho_{\mathrm{p}}}{\rho_{\mathrm{g}}}\frac{c_{\mathrm{s}}}{R_{\mathrm{H}}\Omega}\frac{l_{\mathrm{mfp}}}{R_{\mathrm{H}}}St_{0}\right)^{1/8}R_{\rm H}.$$

(46)

If $x_{\rm ss}$ is much larger than the half width of the horseshoe orbit, $P_{\rm col}$ is given by

$$P_{\rm col}=3\Omega x_{\rm ss}^{2}.$$

(47)

Using Eqs. (46) and (47), we obtain

$$\frac{P_{\mathrm{col},\mathrm{ss}}}{R_{\mathrm{H}}^{2}\Omega}=2\sqrt{6}C_{1}\left(\frac{\rho_{\mathrm{p}}}{\rho_{\mathrm{g}}}\frac{c_{\mathrm{s}}}{R_{\mathrm{H}}\Omega}\frac{l_{\mathrm{mfp}}}{R_{\mathrm{H}}}St_{0}\right)^{1/4}.$$

(48)

On the other hand, if the terminal velocity of particles is smaller than the sound speed, we consider $St\approx St_{0}$ even in the encounter. In that case, the Stokes gas drag, the first term in Eq. (23), is effective. The terminal velocity is given by

$$u_{\mathrm{term}}=\frac{GM_{\mathrm{p}}}{r^{2}}\frac{St_{0}}{\Omega}.$$

(49)

Same as above, Eq. (49) gives

$$x_{\rm ss}=(2C_{1}St_{0})^{1/3}R_{\mathrm{H}}.$$

(50)

Using Eqs. (47) and (50), we obtain the collision rate,

$$\frac{P_{\mathrm{col},\mathrm{set}}}{R_{\mathrm{H}}^{2}\Omega}=3\left(2C_{1}St_{0}\right)^{2/3}.$$

(51)

According to the our simulations, $C_{1}=1.5$.

Finally, we derive the analytic formula considering the effects of the horseshoe and outflow around the Bondi radius. In our simulation, the outflow with a few percent of the sound speed is found around the Bondi radius, as shown in the previous study (Kuwahara et al., 2019). For smaller $St_{0}$, the horseshoe and outflow around the Bondi radius is effective. The collision rate is then given by $P_{\mathrm{col}}=2\Delta x_{\rm s}|v_{0,y}|$, where $\Delta x_{\rm s}=|x_{\rm ss}-r_{\mathrm{HS}}|$ is the difference between the impact parameter and the half horseshoe width $r_{\mathrm{HS}}$. The particles passing around $r\sim R_{\rm B}$ are accreted onto the planet if they can enter the Bondi radius (see Fig. 5d). We consider the passing particles are distributed from $r=R_{\rm B}$ to $r=R_{\rm B}+\Delta r$. The mass conservation gives $\Delta x_{\rm s}|v_{0,y}|=\Delta rv_{\theta}$. If the particle at $r=R_{\rm B}+\Delta r$ enters the Bondi radius with the radial drift of $\Delta r$ during the encounter, the particle then accrete onto the planet. Therefore, the accretion condition is given by the passing time comparable to or longer than the drift time; $R_{\rm B}/v_{\theta}\gtrsim\Delta r/v_{r}$. We thus estimate $P_{\rm col}$ as

$$P_{\rm col}=2R_{\rm B}v_{r}.$$

(52)

Assuming the outflow velocity of $\xi c_{\rm s}$ with $\xi\ll 1$, $v_{r}$ is expressed by the terminal velocity $u_{\rm term}$ and the outflow velocity $\xi c_{\rm s}$, given by

$$v_{r}=\frac{GM_{\rm P}}{r^{2}}\frac{St_{0}}{\Omega}-\xi c_{\rm s}.$$

(53)

Using Eqs. (52) and (53), we obtain

$$\frac{P_{\mathrm{col},\mathrm{ho}}}{R_{\mathrm{H}}^{2}\Omega}=2\frac{R_{\mathrm{B}}}{R_{\mathrm{H}}}\left[\frac{3St_{0}}{\left(R_{\mathrm{B}}/R_{\mathrm{H}}\right)^{2}}-\xi\sqrt{\frac{3}{\left(R_{\mathrm{B}}/R_{\mathrm{H}}\right)}}\right].$$

(54)

We adopt $\xi=0.1m$, based on the results of hydrodynamic simulations for $m=0.03-0.1$. Note that the radial outflow velocity determined by $\xi$ is much smaller than the flow velocity around $r\sim R_{\rm B}$ derived by Kuwahara et al. (2019). We use the $m$ dependence of $\xi$ given by their formula.