Distribution of captured planetesimals in circumplanetary gas disks and implications for accretion of regular satellites

TL;DR

This study investigates how captured planetesimals distribute in circumplanetary disks and their role in forming regular satellites, highlighting the importance of gas drag, orbital dynamics, and planetesimal properties.

Contribution

It provides the first detailed orbital integration analysis of captured planetesimal distribution and their potential contribution to satellite formation.

Findings

Captured planetesimals reach an equilibrium distribution.

Prograde orbits dominate over retrograde ones.

Captured planetesimals can significantly enhance solid material density.

Abstract

Regular satellites of giant planets are formed by accretion of solid bodies in circumplanetary disks. Planetesimals that are moving on heliocentric orbits and are sufficiently large to be decoupled from the flow of the protoplanetary gas disk can be captured by gas drag from the circumplanetary disk. In the present work, we examine the distribution of captured planetesimals in circumplanetary disks using orbital integrations. We find that the number of captured planetesimals reaches an equilibrium state as a balance between continuous capture and orbital decay into the planet. The number of planetesimals captured into retrograde orbits is much smaller than those on prograde orbits, because the former ones experience strong headwind and spiral into the planet rapidly. We find that the surface number density of planetesimals at the current radial location of regular satellites can be significantly enhanced by gas drag capture, depending on the velocity dispersions of planetesimals and the width of the gap in the protoplanetary disk. Using a simple model, we also examine the ratio of the surface densities of dust and captured planetesimals in the circumplanetary disk, and find that solid material at the current location of regular satellites can be dominated by captured planetesimals when the velocity dispersion of planetesimals is rather small and a wide gap is not formed in the protoplanetary disk. In this case, captured planetesimals in such a region can grow by mutual collision before spiraling into the planet, and would contribute to the growth of regular satellites.

Equations

$$\ddot{x}$$

$$\ddot{y}$$

$$\ddot{z}$$

$$\displaystyle\mbox{\boldmath$F$}_{\rm drag}=-\frac{1}{2}C_{\rm D}\pi r^{2}_{\rm s}\rho_{\rm gas}u\mbox{\boldmath$u$},$$

$$\ddot{\tilde{x}}$$

$$\ddot{\tilde{y}}$$

$$\ddot{\tilde{z}}$$

$$\displaystyle\mbox{\boldmath$\tilde{a}$}_{\rm drag}=\mbox{\boldmath$a$}_{\rm drag}/(R_{\rm H}\Omega^{2})=-\frac{3}{8}C_{\rm D}\frac{\rho_{\rm gas}}{\tilde{r}_{\rm s}\rho_{\rm s}}\tilde{u}\mbox{\boldmath$\tilde{u}$},$$

$$E=\frac{1}{2}\left({\dot{\tilde{x}}}^2+{\dot{\tilde{y}}}^2+{\dot{\tilde{z}}}^2\right)+U(\tilde{x},\tilde{y},\tilde{z}),$$

$$U(\tilde{x},\tilde{y},\tilde{z})=-\frac{1}{2}\left(3{\tilde{x}}^2-{\tilde{z}}^2\right)-\frac{3}{\tilde{R}}+\frac{9}{2}.$$

$${\rho }_{\mathrm{gas}}=\frac{\Sigma }{\sqrt{2\pi }h}\exp \left(-\frac{z^2}{2h^2}\right),$$

$$\Sigma ={\Sigma }_{\mathrm{rd}}{\left(\frac{r}{r_{\mathrm{d}}}\right)}^{-p},c_{\mathrm{s}}=c_{\mathrm{rd}}{\left(\frac{r}{r_{\mathrm{d}}}\right)}^{-q/2}$$

$$v_{\mathrm{gas}}=(1-\eta )v_{\mathrm{K}},$$

$$\eta =\frac{1}{2}\frac{h^2}{r^2}\left(p+\frac{q+3}{2}+\frac{q}{2}\frac{z^2}{h^2}\right).$$

$$\displaystyle\mbox{\boldmath$\tilde{a}$}_{\rm drag}=-\zeta\tilde{r}^{-\gamma}{\rm exp}\left(-\frac{\tilde{z}^{2}}{2\tilde{h}^{2}}\right)\tilde{u}\mbox{\boldmath$\tilde{u}$},$$

$$\begin{array}{rl}\zeta & \equiv \frac{3}{8\sqrt{2\pi }}\frac{C_{\mathrm{D}}}{r_{\mathrm{s}}{\rho }_{\mathrm{s}}}\frac{{\Sigma }_{\mathrm{rd}}}{{\tilde{h}}_{\mathrm{rd}}}{d}^{\gamma }\\ & =3\times 10^{-7}{C}_{\mathrm{D}}{\left(\frac{r_{\mathrm{s}}}{1\mathrm{km}}\right)}^{-1}{\left(\frac{{\rho }_{\mathrm{s}}}{1\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}\\ & \times \left(\frac{{\Sigma }_{\mathrm{rd}}}{1\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}\right){\left(\frac{{\tilde{h}}_{\mathrm{rd}}}{0.06}\right)}^{-1}{\left(\frac{d}{0.2}\right)}^{\gamma }.\end{array}$$

$$r_{\mathrm{s}}\simeq 0.03\left(\frac{{\Sigma }_{\mathrm{rd}}}{1{\mathrm{g}\mathrm{cm}}^{-2}}\right){\left(\frac{{\rho }_{\mathrm{s}}}{1{\mathrm{g}\mathrm{cm}}^{-3}}\right)}^{-1}{\zeta }^{-1}\mathrm{cm}.$$

$${\Sigma }_{\mathrm{dust}}$$

$${\delta }_{\mathrm{pla}}={\Sigma }_{\mathrm{pla}}/{\Sigma }_{\mathrm{dust}},$$

$$T_{\mathrm{damp}}=m_{\mathrm{s}}u/F_{\mathrm{drag}}=\frac{2m_{\mathrm{s}}h}{C_{\mathrm{D}}\pi r_{\mathrm{s}}^2\Sigma u},$$

$$T_{\mathrm{vs}}=\frac{1}{n_{\mathrm{pla}}\pi r_{\mathrm{g}}^{2}v},$$

$$T_{\mathrm{vs}}=\frac{v^4}{4G^2{m}_{\mathrm{s}}^2\pi n_{\mathrm{s},\mathrm{cap}}{\Omega }_{\mathrm{p}}\ln \Lambda },$$

$$v_{\mathrm{eq}}$$

$$T_{\mathrm{col}}$$

$$T_{\mathrm{col}}/T_K=\frac{1}{2n_{\mathrm{s},\mathrm{cap}}{\pi }^2{r}_{\mathrm{s}}^2\left(1+v_{\mathrm{esc}}^2/v_{\mathrm{eq}}^2\right)}{\left(\frac{{\Omega }_{\mathrm{p}}}{\Omega }\right)}^{-1}.$$

$$T_{\mathrm{decay}}/T_{\mathrm{K}}=\frac{1}{2\pi \eta }\frac{1+g^2}{2g}{\left(\frac{{\Omega }_{\mathrm{p}}}{\Omega }\right)}^{-1},$$

$$g=a_{\mathrm{drag}}/(u{\Omega }_{\mathrm{p}})=\zeta \eta {\tilde{r}}^{-\gamma +1}.$$

Full text

††slugcomment: Draft version March 23, 2017

1. INTRODUCTION

Regular satellites of giant planets have nearly circular and coplanar prograde orbits, and are thought to have formed in circumplanetary disks (Canup & Ward, 2009; Estrada et al., 2009). Gas and solids are supplied into the circumplanetary disk from the vicinity of the giant planet’s orbit in the protoplanetary disk. Canup & Ward (2002) proposed the so-called gas-starved disk model for the formation of regular satellites of giant planets, where the satellites are formed in a waning circumplanetary disk at the very end stage of giant planet formation. Canup & Ward (2006) performed N-body simulation of satellite formation based on the above model, and showed that mass fraction of satellite system relative to the host planet is regulated to $\sim 10^{-4}$, which is consistent with observations. Later studies revised this model and tried to explain the difference between the satellite systems of Jupiter and Saturn (Sasaki et al., 2010; Ogihara & Ida, 2012). On the other hand, the so-called solid enhanced minimum mass disk model proposed by Mosqueira & Estrada (2003) considers a circumplanetary disk consisting of two different parts of surface density in order to explain the rock-ice ratio of the Galilean satellites. Recently, satellite formation based on this model was examined semi-analytically (Miguel & Ida, 2016). Although growth and orbital evolution of satellites in the circumplanetary disk have been investigated in detail by these works, the distribution of solid materials that are building blocks of satellites in circumplanetary disks has been poorly understood.

There are two kinds of processes for the supply of solid materials into circumplanetary disks. If incoming bodies are sufficiently small to be coupled with the gas flow, they are brought to the circumplanetary disk with the gas inflow from the protoplanetary disk (Canup & Ward, 2002). Tanigawa et al. (2014) examined capture process of solid materials with various sizes, using results of hydrodynamic simulations of gas flow around a growing giant planet (Machida et al., 2008; Tanigawa et al., 2012) and three-body orbital integration. They found that accretion efficiency of dust-size particles initially on circular orbits is rather low when the gas surface density of the protoplanetary disk is that of the minimum mass nebula model. They also found that such small solid bodies can be supplied to the circumplanetary disk when the gas surface density becomes smaller as a results of disk dispersal. Another process for the supply of solid bodies is capture of planetesimals by gas drag from the circumplanetary disk (Fujita et al., 2013; Tanigawa et al., 2014; D’Angelo & Podolak, 2015). In the case of large planetesimals that are decoupled from the gas flow, effects of gas drag on their orbits become significant only when they pass through the dense part of the circumplanetary disk in the vicinity of the planet. Fujita et al. (2013) performed three-body orbital integrations and examined capture of planetesimals from their heliocentric orbits into the circumplanetary disk. They found that planetesimals approaching the circumplanetary disk in the retrograde direction (i.e., in the direction opposite to the motion of the gas in the disk) are more easily captured by gas drag owing to the larger velocity relative to the gas. They also found that capture rates decrease with increasing size of planetesimals because strong gas drag is needed to capture large ones (see also Tanigawa et al., 2014). D’Angelo & Podolak (2015) examined orbital evolution of planetesimals in the vicinity of Jupiter’s orbit, using global orbital integration including gas drag, disruption via ram pressure, and mass loss through ablation. They found that planetesimals are captured from both interior and exterior orbits. Although they also show the distribution of orbital elements of captured planetesimals in the circumjovian disk, they did not examine their orbital evolution in the disk.

Recently, we examined the capture and subsequent orbital evolution of planetesimals that are decoupled from the gas flow in circumplanetary gas disks using three-body orbital integration. In Suetsugu et al. (2016), we investigated the orbital evolution and the distribution of orbital elements of planetesimals captured by relatively strong gas drag, where planetesimals lose sufficient amount of energy to become captured in a single encounter with the planet. On the other hand, in Suetsugu & Ohtsuki (2016), we examined the capture and subsequent orbital evolution of irregular satellites in waning circumplanetary gas disks. We found that irregular satellites are easily captured and are likely to survive when the gas drag is strong and the disk dispersal takes place in a short timescale. However, in these works, radial distribution of solid materials in the circumplanetary disk was not studied in detail. Furthermore, the abundance of captured planetesimal-size bodies relative to dust-size particles in the circumplanetary disk is important to understand which solid materials mainly contribute to the growth of regular satellites. The distribution and size of solid materials in the circumplanetary disk would also influence the timescale of satellite accretion and their chemical composition such as the ice-rock ratio (Sekine & Genda, 2012; Dwyer et al., 2013).

In the present work, as an extension of our recent works (Suetsugu et al., 2016; Suetsugu & Ohtsuki, 2016), we examine radial distribution of planetesimals captured in circumplanetary disks. We perform orbital integration for planetesimals that are decoupled from the gas inflow but are affected by gas drag from the circumplanetary gas disk. We examine capture of planetesimals from their heliocentric orbits by gas drag from the circumplanetary gas disk and orbital evolution of captured planetesimals in the disk. Using our numerical results, we obtain distribution of the surface number density of captured planetesimals in the circumplanetary disk. In Section 2, we describe basic equations, disk model, and numerical methods used in the present work. We show dynamical evolution and distribution of planetesimals captured from heliocentric circular orbits in Section 3. In Section 4, we investigate the effect of the dynamical states of planetesimals on the surface number density. In Section 5, we discuss implications of our numerical results for accretion of regular satellites. Section 6 summarizes our results.

2. THE MODEL

2.1. Numerical Procedures

We consider a local coordinate system ($x,y,z$) centered on a planet (Figure 1). We assume that the planet ($M$) is on a circular orbit with semi-major axis $a_{0}$, and is embedded in a disk of planetesimals. Also, the planet is assumed to have a circumplanetary gas disk, whose mid-plane coincides with the planet’s orbital plane. In the present work, we assume a disk of equal-mass planetesimals ($m_{\rm s}\ll M$), and neglect gravitational interaction among them. At azimuthal locations far from the planet, planetesimals are assumed to have uniform radial distribution. Planetesimals are supplied through the azimuthal boundaries at $y=\pm 50R_{\rm H}$ ($R_{\rm H}$ is the planet’s Hill radius; $R_{\rm H}=a_{0}h_{\rm H}=a_{0}\{\left(M+m_{\rm s}\right)/(3M_{\odot})\}^{1/3}$), which is sufficiently far from the planet to neglect its gravitational effect. Supplied planetesimals are assumed to be distributed radially uniformly in the protoplanetary disk; we examine the effect of non-uniform radial distribution in Section 4. We integrate the orbits of planetesimals by solving Hill’s equation with the effect of gas drag from the circumplanetary disk. The motion of the planetesimal is expressed as (Ohtsuki, 2012; Suetsugu & Ohtsuki, 2016),

[TABLE]

where $\Omega$ is the planet’s orbital angular frequency, and $R=\sqrt{x^{2}+y^{2}+z^{2}}$ is the distance between the centers of the planet and the planetesimal. \mbox{\boldmatha}_{\rm drag}\equiv\mbox{\boldmathF}_{\rm drag}/m_{\rm s} is the acceleration due to the gas drag force \mbox{\boldmathF}_{\rm drag} given by

[TABLE]

where $C_{\rm D}$ is the drag coefficient (we assume $C_{\rm D}=1$; Appendix A), $r_{\rm s}$ is the radius of the planetesimal, $\rho_{\rm gas}$ is the gas density, and u=|\mbox{\boldmathu}| is the velocity of the planetesimal relative to the gas. We scale the distance by mutual Hill radius $R_{\rm H}$ and time by $\Omega^{-1}$. Then, we can express the above equation of motion in a non-dimensional form as

[TABLE]

where tildes denote non-dimensional quantities. The non-dimensional acceleration due to gas drag can be written as

[TABLE]

where $\rho_{\rm s}$ is the internal density of planetesimals. When the gas drag can be neglected, we can confirm from Equation (2.1) the conservation of energy given as

[TABLE]

with

[TABLE]

We remove planetesimals when their distance from the planet becomes sufficiently large or when collision with the planet is detected. We assume that the physical size of the planet ($r_{\rm p}$) relative to its Hill radius, $r_{\rm p}/R_{\rm H}$, is $10^{-3}$, which corresponds to the physical size of a planet at Jupiter’s orbit.

The number of planetesimals in the circumplanetary disk increases due to capture by gas drag and decreases due to collision with the planet after orbital decay. However, as we will show below, the number reaches an equilibrium state by a balance between continuous capture of planetesimals and loss to the planet due to orbital decay. Since the timescale needed to reach a quasi-steady state depends on the dynamical properties of planetesimals and the strength of gas drag, we perform orbital integration for a sufficiently long period of time ($200-500T_{\rm K}$, where $T_{\rm K}$ is the orbital period of the planet).

Unperturbed surface number density of the planetesimal disk $n_{\rm s}$ is one of the important parameters in order to determine the distribution of captured planetesimals in the circumplanetary disk. Thus, we briefly estimate the value of $n_{\rm s}$. We can obtain $n_{\rm s}$ when we set the solid surface density in the protoplanetary disk $\Sigma_{\rm ppd,solid}$ and the mass of the planetesimal $m_{\rm s}$ (recall that we assume equal-sized planetesimals). If we adopt the minimum mass solar nebula model (Hayashi, 1981), the surface density of solids in the protoplanetary disk is $\Sigma_{\rm ppd,solid}\simeq 3.3$ gcm${}^{-2}\simeq 9\times 10^{25}$ g$R_{\rm H,J}^{-2}$ at the orbit of Jupiter ($R_{\rm H,J}$ is the Hill radius of Jupiter). In the present work, we will examine the distribution of captured planetesimals in circumplanetary disks for cases with a wide range of planetesimal sizes ($r_{\rm s}\sim 1-10^{4}$ m). If we assume that the size of the planetesimal is $r_{\rm s}=300$ m, the surface number density scaled by the Hill radius of Jupiter can be written by $\tilde{n}_{\rm s}=8.2\times 10^{11}\left(\rho_{\rm s}/1{\rm gcm^{-3}}\right)^{-1}\left(r_{\rm s}/300{\rm m}\right)^{-3}$ $R_{\rm H,J}^{-2}$. However, if we adopt such a realistic surface number density, the number of orbits to integrate per one timestep is $\sim 10^{14}$ (Even if we assume a large size of planetesimals to reduce the surface number density to be $\tilde{n}_{\rm s}=8.2\times 10^{5}$, we have to integrate $\sim 10^{8}$ orbits). Here, we should recall that we neglect mutual gravity between planetesimals, thus their orbits are not influenced by the choice of $\tilde{n}_{\rm s}$. Therefore, we adopt a value of $\tilde{n}_{\rm s}$ feasible for our simulation (typically, $\tilde{n}_{\rm s}=400$), and discuss the enhancement of the surface number density in the circumplanetary disk relative to the background value 111In this case, the number of planetesimals in the simulation cell is $\sim 10^{4}$ for the case of planetesimals initially on circular orbits, and is $\sim 10^{5}$ orbits for the case of planetesimals initially on eccentric and inclined orbits. We confirmed that the surface number density of captured planetesimals in the circumplanetary disk is unchanged when $\tilde{n}_{\rm s}$ is increased to 1000. Therefore, the above value of $\tilde{n}_{\rm s}$ is sufficiently large to accurately obtain the distribution of captured planetesimals in the circumplanetary disk..

2.2. Gas Drag

For those planetesimals that are large enough to be decoupled from the inflowing gas, gas drag from the circumplanetary disk becomes important only when they pass through the dense part of the disk in the vicinity of the planet, where the gas distribution can be approximated to be axisymmetric (Tanigawa et al., 2012; Fujita et al., 2013). Thus, in the present work, we assume that the radial distribution of the gas density of the circumplanetary disk is axisymmetric, its radial dependence is given by a power law, and its vertical structure is isothermal. Under these assumptions, the gas density can be written by

[TABLE]

where $h=c_{\rm s}/\Omega_{\rm p}$ is the scale height of the circumplanetary disk ($\Omega_{\rm p}$ is the Keplerian orbital frequency around the planet), and

[TABLE]

are the gas surface density and sound velocity, respectively, with $r=\sqrt{x^{2}+y^{2}}$ being the horizontal distance from the planet in the mid-plane. In the above, $r_{\rm d}=dR_{\rm H}$ is a typical length scale roughly corresponding to the effective size of the circumplanetary disk, and $\Sigma_{\rm rd}$ and $c_{\rm rd}$ are the surface density and sound velocity at that radial location. In our calculations, we set $d=0.2$ and $p=3/2$ based on results of hydrodynamic simulations (Machida et al., 2008; Tanigawa et al., 2012), and also assume $q=1/2$ as a simple model (Fujita et al., 2013; Suetsugu & Ohtsuki, 2016; Suetsugu et al., 2016). In our orbital calculations we turn on gas drag when planetesimals enter the planet’s Hill sphere, in order to avoid effects of artificial cutoff at $r=r_{\rm d}$. Because the gas density decreases rapidly with increasing distance from the planet, this assumption on gas drag does not affect results of our calculations222For those small planetesimals whose motion is strongly affected by the gas flow, the effect of gas drag is important even outside of the planet’s Hill sphere (Tanigawa et al., 2014). However, calculations of long-term evolution of captured bodies including such an effect are time-consuming, and in the present work, we focus on the orbital evolution under gas drag within the Hill sphere, by making the above assumption.. Gas elements in the disk are assumed to rotate in circular orbits around the planet with a velocity slightly lower than the Keplerian velocity due to radial pressure gradient, i.e.,

[TABLE]

where $v_{\rm K}$ is the Keplerian velocity around the planet at the radial location considered. Using Equations (7) and (8), $\eta$ can be written as (Tanaka et al., 2002)

[TABLE]

When the gas density is given by Equation (7), Equation (4) can be rewritten as (Fujita et al., 2013)

[TABLE]

where $h=h_{\rm rd}\left(r/r_{\rm d}\right)^{(3-q)/2}$ with $h_{\rm rd}$ being the scale height at $r=r_{\rm d}$, and $\gamma\equiv p+(3-q)/2$, and $\zeta$ is the non-dimensional parameter representing the strength of gas drag defined by

[TABLE]

We set $\tilde{h}_{\rm rd}\equiv h_{\rm rd}/R_{\rm H}=0.06$ in the present work (Tanigawa et al., 2012; Fujita et al., 2013) 333In the following, we assume icy planetesimals with $\rho_{\rm s}=1$gcm*-3* when converting $\zeta$ into $r_{\rm s}$. The value of $\zeta$ for rocky planetesimals becomes smaller than that for icy planetesimals with the same size. This leads to smaller capture rates, while the lifetime of captured bodies in the circumplanetary disk becomes longer. As a result, the influence of the material density on the radial distribution of captured bodies is expected to be small.. If we assume the gas surface density based on the gas-starved disk model (Canup & Ward, 2002), we have $\Sigma_{\rm rd}\simeq 1$gcm*-2*. For the above fiducial values of $C_{\rm D}$, $d$, and $\tilde{h}_{\rm rd}$, the relation between $\zeta$ and the planetesimal size $r_{\rm s}$ can be written as (see Fujita et al., 2013; Suetsugu & Ohtsuki, 2016)

[TABLE]

In the present work, we will examine cases with $\zeta=10^{-8}-10^{-4}$, which roughly corresponds to planetesimal sizes on the order of $r_{\rm s}\sim 1-10^{4}$ m.

3. CASE OF INITIALLY HELIOCENTRIC CIRCULAR ORBITS

3.1. Supply of Planetesimals into Circumplanetary Disk

Figure 2 shows an example of snapshots of the spatial distribution of planetesimals supplied from heliocentric circular orbits in the case of $\zeta=10^{-6}$ ($r_{\rm s}\simeq 300$ m). Figure 2(a) shows the distribution of planetesimals in the vicinity of the planet. Most of planetesimals are scattered by the gravity of the planet or pass by the planet’s Hill sphere, while some of them enter the Hill sphere. Planetesimals that enter the planet’s Hill sphere are affected by gas drag from the circumplanetary gas disk, and those planetesimals that lose a sufficient amount of energy become permanently captured (Figure 2(b)). In Figure 2, blue and green marks show planetesimals permanently captured by gas drag ($E<0$); blue and green represent those captured into prograde and retrograde orbits, respectively. On the other hand, gray dots represent free planetesimals just passing by the planet. Figure 2(c) and (d) show blow-ups of the distribution in the vicinity of the planet. Typically, those captured into retrograde orbits tend to have larger velocity relative to the gas than those on prograde orbits. Thus, planetesimals in retrograde orbits spiral into the planet more quickly, and the innermost part of the circumplanetary disk is dominated by planetesimals captured into prograde orbits (Figure 2(d)). Also, free planetesimals rarely exist in the vicinity of the planet, because they become permanently captured by strong gas drag in such a dense part of the disk.

Figure 3 shows evolution of the number of planetesimals permanently captured within the planet’s Hill sphere for the case shown Figure 2. The red curve represents the total number of captured planetesimals. The number of captured planetesimals $N_{\rm cap}$ can be determined by a balance between the increase by the capture of planetesimals and the decrease due to collision with the planet after orbital decay. We find that the total number first increases gradually until $t\simeq 35T_{\rm K}$ due to capture of planetesimals. After $t\simeq 35T_{\rm K}$, the system reaches an equilibrium state with $N_{\rm cap}\simeq 2000$. In an unperturbed state without the planet and the circumplanetary disk, the number of planetesimals passing within the planet’s Hill sphere is $\tilde{n}_{\rm s}R_{\rm H}^{2}\pi\simeq 1200$ when $\tilde{n}_{\rm s}=400$. Also, the number of free planetesimals within the Hill sphere is about 1200 in this case. Thus, the number density of planetesimals within the Hill sphere is significantly enhanced due to the effects of the planet and the disk. The blue and green lines represent the number of planetesimals permanently captured in the prograde direction and the retrograde direction, respectively. We find that the number of planetesimals captured in retrograde orbits is significantly smaller because of their rapid orbital decay into the planet.

On close inspection of Figure 3, we find that the number of bodies captured in prograde orbits increases very rapidly at the very beginning ($t\lesssim 5T_{\rm K}$), and then its growth rate slows down and continues gradual increase at $5T_{\rm K}\lesssim t\lesssim 35T_{\rm K}$. The rapid increase in the first stage is due to capture of bodies in a single encounter with the planet (Figure 4; red and green lines). On the other hand, the gradual increase in the second stage comes from a different type of capture orbits, which maintain larger pericenter distance and survive longer in the disk (Figure 4, blue line). Dynamical behavior of these orbits is described in more detail in Appendix B.

3.2. Radial Distribution of Captured Planetesimals

Figure 5 shows the surface number density of captured planetesimals $n_{\rm s,cap}$ scaled by the background surface number density $n_{\rm s}$, as a function of their radial distance from the planet. These plots are created by averaging the distributions at different times after the system reached the equilibrium state. The distribution of planetesimals captured into prograde orbits (blue line) show different slopes depending on the radial distance. In the outer region with $r/R_{\rm H}\gtrsim 0.04$, the distribution can be roughly approximated by a power law $n_{\rm s,cap}\propto r^{-\alpha}$ with $\alpha\simeq 2$. Then, with decreasing distance, the slope becomes steeper ($\alpha\simeq 4.5$) at $0.01\lesssim r/R_{\rm H}\lesssim 0.04$, and then the distribution levels off at $r/R_{\rm H}\simeq 0.01$. With further decreasing distance, it increases again and, finally, levels off again at $r/R_{\rm H}\lesssim 0.003$. We find that the above behavior of the surface number density distribution reflects dynamical evolution of captured planetesimals in the circumplanetary disk. Those captured by a single encounter with the dense part of the gas disk accumulate in the vicinity of the planet and are short-lived because of rapid orbital decay, while those captured in the outer part are typically long-lived, with their pericenter avoiding penetration into the dense part of the disk. Captured bodies in the outer part ($r/R_{\rm H}\gtrsim 0.01$) are dominated by the latter long-lived ones. Immediately after being captured, their orbits decay rather rapidly due to large orbital eccentricities; this corresponds to the region with $\alpha\simeq 2$ ($r/R_{\rm H}\gtrsim 0.04$). Then, their eccentricities are damped by gas drag and orbits become nearly circular. Then the rate of orbital decay decreases (Adachi et al., 1976; Kobayashi, 2015); the region with the steeper slope with $\alpha\simeq 4.5$ ($0.01\lesssim r/R_{\rm H}\lesssim 0.04$) corresponds to such slowly migrating bodies on the relatively long-lived orbits. With further decreasing radial distance, the rate of orbital decay increases because of higher gas density, which leads to the nearly flat distribution at $r/R_{\rm H}\lesssim 0.01$.

On the other hand, the increase of the surface number density in the inner region ($0.003\lesssim r/R_{\rm H}\lesssim 0.007$) is caused by planetesimals captured into short-lived orbits (the red and green lines in Figure 4). In fact, the outer boundary of this region roughly corresponds to the prograde capture radius where planetesimals can be captured by a single passage through the gas disk (Fujita et al., 2013). Since such captured planetesimals spiral into the planet quickly due to strong gas drag (Appendix B), the surface number density distribution becomes nearly flat at $r/R_{\rm H}\lesssim 0.003$.

The distribution of planetesimals captured into retrograde orbits can be divided into two regions: the outer region at $r/R_{\rm H}\gtrsim 0.01$, where the distribution can be approximated by a single power-law with $\alpha\simeq 1.2$, and the inner region with nearly flat distribution. In the outer region, bodies are captured into relatively long-lived orbits but migrate inward more rapidly than the prograde ones in the same region due to strong headwind. Owing to rapid orbital decay, the distribution levels off inside the critical radius for capture into retrograde orbits by a single encounter (Fujita et al. 2013; green vertical dashed line).

We also plot the distribution of free planetesimals for comparison. Although its surface number density is slightly enhanced due to the effect of the planet’s gravitational focusing, it is comparable to the surface number density of the background unperturbed planetesimal disk ($n_{\rm s}$). The surface number density of free planetesimals vanishes inside the prograde capture radius ($r\lesssim 0.005R_{\rm H}$), because planetesimals become permanently captured in that region due to gas drag.

3.3. Dependence of Distribution on Planetesimal Size

Next, we examine the dependence of the distribution of captured planetesimals on the size of planetesimals, by changing the values of the gas drag parameter ($\zeta=10^{-8}-10^{-4}$). Figures 6 show the evolution of the number of captured planetesimals for $\zeta=10^{-5}$ (corresponding to $r_{\rm s}\simeq 30$ m for the gas-starved disk model), $10^{-7}$ ($r_{\rm s}\simeq 3$ km), and $10^{-8}$ $(r_{\rm s}\simeq 30$ km), respectively. General behavior is similar to the case with $\zeta=10^{-6}$ (Figure 3), but, the timescale to reach the equilibrium state increases with deceasing strength of gas drag because the rates of capture and orbital decay are reduced. We find that the number of planetesimals captured in prograde orbits in the case of $\zeta=10^{-8}$ ($r_{\rm s}\simeq 30$ km) is significantly smaller than the other cases, because the capture radius for prograde orbits in this case is smaller than the physical size of the planet. On the other hand, the number of planetesimals captured in retrograde orbits in this case is somewhat larger, because the lifetime of planetesimals in retrograde orbits becomes longer under such weak gas drag. Figure 7 shows the number of captured planetesimals in the disk in the equilibrium state as a function of $\zeta$ (or $r_{\rm s}$). We find that the number hardly depends on the planetesimal size for the intermediate values of $\zeta$, while it decreases for small and large values of $\zeta$ for the reasons described above. Figure 8 shows the surface number density distribution for various values of $\zeta$. The general features of the distribution are similar to the case of $\zeta=10^{-6}$ shown in Figure 5, which suggests that these features are common for a wide range of gas drag parameters. In the case of $\zeta=10^{-5}$, planetesimals can be captured by gas drag even in the outer parts of the disk where the gas density is rather low. Thus, the surface number density of captured planetesimals shifts radially outward. On the other hand, when $\zeta=10^{-7}$, planetesimals need to pass through the denser part of the disk in order to be captured, thus, the distribution shifts radially inward. In the case of $\zeta=10^{-8}$ shown in Figure 8(c), the inner region corresponding to single-encounter capture into prograde orbits disappears because the prograde capture radius is smaller than the planet’s physical radius.

In some cases, captured planetesimals change their direction of orbital motion from retrograde to prograde due to strong headwind after they become permanently captured (e.g., Tanigawa et al., 2014; Suetsugu et al., 2016). In the case of a planet at Jupiter’s orbit, the change of the orbital direction takes place before planetesimals spiral into the planet when $\zeta\gtrsim 10^{-5}$ (Suetsugu et al., 2016). The above range of $\zeta$ corresponds to $r_{\rm s}\lesssim 10$ m if we assume the gas density based on the gas-starved disk model (Canup & Ward, 2002). Figure 9 shows the surface number density of planetesimals captured into retrograde orbit as a function of radial distance for the cases with $\zeta=10^{-4}$ and $10^{-5}$. The number density vanishes in the vicinity of the planet because planetesimals captured in retrograde orbits change the direction of orbital motion. We confirmed that the surface number density of captured planetesimals on prograde orbits increases due to this effect in the innermost region of the disk when $\zeta=10^{-4}$.

4. DEPENDENCE OF THE DISTRIBUTION OF CAPTURED PLANETESIMALS ON DYNAMICAL STATES OF INCOMING PLANETESIMALS

So far, we have examined capture of planetesimals that are initially on heliocentric circular orbits with uniform radial distribution in the protoplanetary disk. Here, we examine effects of dynamical states of incoming planetesimals, i.e., eccentricities and inclinations of their initial heliocentric orbits, and non-uniform radial distribution in the protoplanetary disk. We assume that eccentricities and inclinations follow a Rayleigh distribution with given rms values (we assume $\langle e_{\rm H}^{2}\rangle^{1/2}=2\langle i_{\rm H}^{2}\rangle^{1/2}$, where $e_{\rm H}$ and $i_{\rm H}$ are the eccentricity and inclination scaled by $R_{\rm H}/a_{0}$), and that orbital phase angles are randomly distributed in [0-2$\pi$]. As for the radial distribution of planetesimals in the protoplanetary disk, a sufficiently massive planet opens a gap both in the gas and solid component of the disk, and interactions between solid bodies and the gas near the edge of the gap is important for the supply of solids into the circumplanetary disk (Paardekooper & Mellema, 2006; Rice et al., 2006; Ayliffe et al., 2012; D’Angelo & Podolak, 2015). Exact treatment of such interactions requires global simulation of such a system, which is beyond the scope of the present work. Instead, in the present work, we adopt a simple model for the non-uniform radial distribution of planetesimals in the protoplanetary disk (Fujita et al., 2013), by assuming that planetesimals in the vicinity of the planet’s orbit have been removed and a gap centered on the planet’s orbit with a half width $W_{\rm gap}$ (scaled by $R_{\rm H}$) is formed.

In Figure 10, the solid lines represent results for the cases with uniform radial distribution of planetesimals with different velocity dispersions; $\zeta=10^{-6}$ ($r_{\rm s}\simeq 300$ m) in Panels (a), (b) and (c), while $\zeta=10^{-5}$ ($r_{\rm s}\simeq 30$ m) in (d), (e), and (f). Figures 10(a) and (d) show the distribution of planetesimals captured in prograde orbits. We find that the surface number density in the inner region that corresponds to capture by a single encounter ($r\lesssim 0.02R_{\rm H}$) decreases significantly with increasing velocity dispersion (Fujita et al., 2013). The bump that we find in the case of initially circular orbits and corresponds to the capture radius is smoothed out in the case with non-zero velocity dispersions, because the capture radius depends on the planetesimals’ velocity dispersion and those with various eccentricities and inclinations are captured in the cases with non-zero velocity dispersions. We find that the dependence of the distribution in the outer region ($r\gtrsim 0.02R_{\rm H}$) on the velocity dispersion is rather weak. As we mentioned above, the distribution in this region is determined by a balance between capture into relatively long-lived orbits whose types depend on the velocity dispersion (Suetsugu et al., 2011; Suetsugu & Ohtsuki, 2013; Suetsugu et al., 2016), and rather rapid orbital decay due to large eccentricities of planet-centered orbits. The resultant enhancement of the surface number density in this region compared to the background value is not large ($n_{\rm s,cap}/n_{\rm s}\lesssim 10$).

Figure 10(b) and  10(e) show similar plots for planetesimals captured into retrograde orbits. In this case, the surface number density for the case with small but non-zero velocity dispersion is larger than that for the initially circular orbits, because of contribution from long-lived capture orbits that only appear in cases of non-zero velocity dispersions (Suetsugu et al., 2016; Suetsugu & Ohtsuki, 2016). With further increase of velocity dispersion, the surface number density decreases for the same reasons in the case of prograde orbits. Figure 10 (c) and (f) show that free planetesimals exist even in the vicinity of the planet in the case with non-zero velocity dispersions, because limited energy dissipation allows them to avoid permanent capture.

In Figure 10, we also plotted with dashed and dotted lines results for the cases with a gap in the planetesimal disk. In this case, only planetesimals initially on orbits with $|b_{\rm H}|\geq W_{\rm gap}$ can approach the planet and the circumplanetary disk, where $b_{\rm H}=(a-a_{0})/R_{\rm H}$ is the difference in the semi-major axes of the planetesimal and the planet scaled by the planet’s Hill radius (Figure 1). Figure 10(a) shows that the surface number density of captured bodies on prograde orbits in the inner region decreases due to the effect of the gap, while that in the outer region is hardly affected. Figure 10(b) and  10(e) show that the surface number densities of bodies on retrograde orbits as well as free planetesimals decrease significantly; thus the dominance of the bodies on prograde orbits become still notable in the presence of the gap.

Most of regular satellites of giant planets have orbits with semi-major axes of $\sim 0.005R_{\rm H}-0.03R_{\rm H}$. Interestingly, our numerical results show that the surface number density of captured planetesimals is significantly enhanced in such regions compared to the background value. We will further discuss implications of our results for satellite accretion in Section 5. Using numerical results presented above, here we calculate averaged surface number density of captured planetesimals $n_{\rm sat}$ for the radial region between $r=0.003R_{\rm H}$ and $r=0.03R_{\rm H}$. Since planetesimal on prograde orbits are dominant in such a region, we calculate $n_{\rm sat}$ using only planetesimals captured into prograde orbits. Figure 11 shows that $n_{\rm sat}$ decreases with increasing random velocity of planetesimals, regardless of the presence of the gap in the protoplanetary disk. Extrapolating on numerical results, we find that enhancement of the surface number density of planetesimals captured by gas drag would be negligible in the radial location corresponding to the current orbits of the regular satellites when $\langle e_{\rm H}^{2}\rangle^{1/2}\gtrsim 8$. On the other hand, in the case of non-uniform radial distribution of planetesimals ($W_{\rm gap}>3$), capture of planetesimals hardly takes place in the shear-dominated regime $\langle e^{2}_{\rm H}\rangle^{1/2}\lesssim 1$ because planetesimals approaching the planet are removed by the gap (Fujita et al., 2013). Therefore, moderate values of the velocity dispersion seem to be preferable for the supply of planetesimals as building blocks of regular satellites in the presence of a gap.

5. IMPLICATIONS FOR SATELLITE ACCRETION

5.1. Planetesimal-to-Dust Ratio in Circumplanetary Disks

In the previous sections, we have shown that planetesimals captured from heliocentric orbits are supplied into the vicinity of the current orbits of regular satellites. Here, we compare the surface density of captured planetesimals with that of dust, in order to examine their relative contribution to satellite accretion.

We assume that the surface density of the gas in the protoplanetary disk at the time of satellite accretion considered here is given by $\Sigma_{\rm ppd,gas}=f_{\rm dep}\times\Sigma_{\rm MMSN}$, where $\Sigma_{\rm MMSN}$ is the surface density of the minimum mass solar nebula (Hayashi, 1981) and $f_{\rm dep}$ is the depletion factor due to disk dispersal. We consider a simple model for a circumplanetary disk at the last stage of giant planet formation as follows. On the basis of hydrodynamic simulations of gas accretion onto a forming giant planet and the gas-starved model for the circumplanetary disk, we set $f_{\rm dep}\simeq 10^{-3}$ (Appendix C). As for the circumplanetary disk, we assume here that the power-law distribution of the gas surface density ($\Sigma$) given by Equation (8) can be extended to $r=R_{\rm H}$, where it equals $\Sigma_{\rm ppd,gas}$, i.e., $\Sigma=\Sigma_{\rm ppd,gas}(r/R_{\rm H})^{-p}$. In this case, if the dust-to-gas ratio in the circumplanetary disk is given by $\chi$, the dust surface density in the circumplanetary disk is given by

[TABLE]

On the other hand, we assume that the ratio of the surface density of planetesimals in the protoplanetary disk ($\Sigma_{\rm ppd,pla}$) to that of the gas is given by $\chi_{\rm ppd}$, i.e., $\Sigma_{\rm ppd,pla}=\chi_{\rm ppd}\Sigma_{\rm ppd,gas}$. As for the surface density distribution of planetesimals in the circumplanetary disk, $\Sigma_{\rm pla}(r)$, we use our numerical results obtained in the previous sections, assuming that it equals $\Sigma_{\rm ppd,pla}$ at $r=R_{\rm H}$, i.e., $\Sigma_{\rm pla}(R_{\rm H})$ = $\chi_{\rm ppd}\Sigma_{\rm ppd,gas}$. Note that we use the total surface number of planetesimals captured into both prograde and retrograde orbits. From $\Sigma_{\rm dust}$ and $\Sigma_{\rm pla}$ obtained above, we can calculate the ratio of $\Sigma_{\rm pla}$ to $\Sigma_{\rm dust}$ in the circumplanetary disk as

[TABLE]

as a function of the radial distance $r$ for a given value of $\chi_{\rm ppd}/\chi$.

Figure 12 shows the surface density of the gas (solid straight line), dust (dashed straight line), and captured planetesimals (curves) for the case of initially on heliocentric circular orbits with $\chi_{\rm ppd}=\chi=0.01$. We find that solid material in the disk is dominated by captured planetesimals in the inner region, while dust particles are dominant in the outer region. The range of the region dominated by captured planetesimals shift inward with decreasing gas drag parameter (i.e., increasing size of planetesimals). Figure 12(b) shows the planetesimal-to-dust ratio as a function of the radial distance, where we can confirm the above tendency. In this figure, we also plot the case with $\chi_{\rm ppd}/\chi=10$ ($\chi_{\rm ppd}=5\times 10^{-2}$ and $\chi=5\times 10^{-3}$), because the dust abundance in the circumplanetary disk would be reduced by, for example, the filtering effect due to accumulation at the gas density bump in the protoplanetary disk (Ayliffe et al., 2012). In this case, the planetesimal-dominated region becomes extended almost to the outer edge of the circumplanetary disk ($r\lesssim 0.2R_{\rm H}$). However, the planetesimal-dominated region shrinks when incoming planetesimals have non-zero velocity dispersions in the protoplanetary disk (Figure 13). The region further shrinks if there is a wide gap in the radial distribution of planetesimals in the protoplanetary disk (Figure 10). From these results, we conclude that the contribution of captured planetesimals would become comparable to or even larger than that of dust particles as building blocks of regular satellites when the velocity dispersion of planetesimals and the width of their gap in the protoplanetary disk are rather small, while their contribution becomes small in the case of larger velocity dispersions and/or a wide gap.

5.2. Timescale of Accretion of Captured Planetesimals

Using the surface density distribution of captured planetesimals that we obtained in the previous sections, we can evaluate the timescale of collision among them and examine whether they grow by accretion before spiraling into the planet. For this purpose, first, we estimate the relative velocity between captured planetesimals, assuming that it is determined by damping due to gas drag and enhancement due to viscus stirring among the planetesimals. The timescale of velocity damping by gas drag is given by

[TABLE]

where $u\simeq(e_{\rm p}^{2}+i_{\rm p}^{2}+\eta^{2})^{1/2}v_{\rm K}$ is the relative velocity between the gas and planetesimals ($i_{\rm p}$ is an inclination of the planet-centered orbit), and we used $\rho_{\rm gas}\simeq\Sigma/h$. On the other hand, the timescale of viscus stirring among captured planetesimals is given by

[TABLE]

where $n_{\rm pla}$ is number density of captured planetesimals in the circumplanetary disk, $r_{\rm g}=2Gm_{\rm s}/v^{2}$ is the impact parameter for the 90∘ deflection, and $v=(e_{\rm p}^{2}+i_{\rm p}^{2})^{1/2}v_{\rm K}$. Since $n_{\rm pla}$ can be approximately given by $n_{\rm pla}\simeq n_{\rm s,cap}\Omega_{\rm p}/v$, $T_{\rm vs}$ is rewritten as

[TABLE]

where ${\rm ln}\Lambda$ comes from the effect of distant encounters, and ${\rm ln}\Lambda\sim 1$ in the dispersion-dominated regime (Ohtsuki et al., 2002). Setting $T_{\rm damp}=T_{\rm vs}$, the equilibrium velocity dispersion $v_{\rm eq}$ is obtained as (Kokubo & Ida, 2000; Kobayashi et al., 2010)

[TABLE]

where $v_{\rm esc}=\sqrt{2Gm_{\rm s}/r_{\rm s}}$ is the escape velocity of planetesimals, and we have neglected the relative velocity due to the pressure gradient.

Using this equilibrium velocity dispersion, the timescale of collision between planetesimals is written as

[TABLE]

Substituting $\sigma_{\rm col}=\pi r_{\rm s}^{2}\left(1+v_{\rm esc}^{2}/v_{\rm eq}^{2}\right)$ for the collision cross section, we have

[TABLE]

On the other hand, the timescale of the orbital decay due to gas drag can be written as (Adachi et al., 1976)

[TABLE]

where

[TABLE]

Figure 14(a) shows $T_{\rm col}$ and $T_{\rm decay}$ as a function of the radial distance from the planet in the case of $\zeta=10^{-5}$ ($r_{\rm s}\simeq 30$m). In the outer region of the circumplanetary disk with $r\gtrsim 0.04R_{\rm H}$, $T_{\rm decay}<T_{\rm col}$, and captured planetesimals migrate inward before they grow by accretion. On the other hand, in the inner region where the surface density rapidly increases ($r\lesssim 0.03R_{\rm H}$; Figure 10), $T_{\rm col}<T_{\rm decay}$ and their mutual collision becomes important.

Figure 14(b) shows the impact velocity $(v_{\rm imp}=\left(v_{\rm eq}^{2}+v_{\rm esc}^{2}\right)^{1/2})$ between captured planetesimals as a function of the radial distance in the above case. Although the impact velocity is higher than the escape velocity, it is rather low ($\sim 0.01$ ms*-1*). Recent impact simulations with the effect of material strength show that bodies can accrete and grow in such low-velocity collisions (Jutzi & Asphaug, 2015). Thus, we can expect collisional growth of captured planetesimals in the inner region of the circumplanetary disk. On the other hand, the surface density of captured planetesimals largely depends on the random velocity of the planetesimals disk. As a result, $T_{\rm col}$ increases with increasing velocity dispersion, and is comparable to or larger than $T_{\rm decay}$ when $\langle e_{\rm H}^{2}\rangle^{1/2}=2-4$, even in the inner part of the disk. In this case, captured planetesimals would spiral into the planet before they grow significantly by accretion.

In the case of $\zeta=10^{-6}$ ($r_{\rm s}\simeq 300$ m; Figure 14(c)), the collision timescale in the inner region of the disk becomes somewhat longer than the case with $\zeta=10^{-5}$ ($r_{\rm s}\simeq 30$ m). The velocity dispersion becomes lower also due to the lower surface density of captured planetesimals (Figure 14(d)). On the other hand, the orbital decay timescale becomes longer, because of weaker gas drag. As a result, $T_{\rm col}$ becomes shorter than $T_{\rm decay}$ in the vicinity of the planet in the inner region of the disk, and growth of planetesimals by accretion can be expected there. As we mentioned in Section 4, the surface number density of captured planetesimals is significantly enhanced in the regions corresponding to the current radial location of regular satellites of giant planets. Our results suggest that once a significant amount of planetesimals are captured in such regions, they would contribute to the growth of regular satellites of giant planets. Further studies including mutual gravitational interactions between captured planetesimals are desirable to quantitatively examine their effect on the growth of satellites.

6. SUMMARY

In the present work, we examined the distribution of planetesimals captured in circumplanetary disks using orbital integration. We found that the number of captured planetesimals in the circumplanetary disk reaches an equilibrium state as a balance between continuous capture from their heliocentric orbits by gas drag and orbital decay into the planet. Typically, the number of planetesimals captured into retrograde orbits is much smaller than that on prograde orbits, because those captured on retrograde orbits experience strong headwind and spiral into the planet rapidly. Using our numerical results, we obtained surface number density distribution of permanently captured planetesimals in the circumplanetary disk in the equilibrium state. We found that the surface number density of captured planetesimals is significantly enhanced at the current location of regular satellites, and that it depends on the dynamical states of planetesimals, i.e., their initial velocity dispersion in the protoplanetary disk and the presence of a gap in their radial distribution; larger velocity dispersions and a wider gap both reduce the efficiency of planetesimal capture, thus decreases their surface number density in the disk.

We also examined the surface density ratio of dust and captured planetesimals in the circumplanetary disk. We found that solid material at the current location of regular satellites is dominated by captured planetesimals when the random velocity of planetesimals is rather small and a wide gap is not formed in the protoplanetary disk. In this case, captured planetesimals can grow by mutual collision before they spiral into the planet due to gas drag. Therefore, captured planetesimals would significantly contribute to the growth of regular satellites. Further studies including mutual gravitational interactions between captured planetesimals are necessary to quantitatively examine their effect on satellite accretion. Also, their importance on satellite accretion relative to dust particles depends on dust-to-planetesimal ratio in the protoplanetary disk and in the circumplanetary disk. Although, the influence of the material density on the radial distribution of captured bodies is expected to be small, mass of captured planetesimals would be reduced by ablation during the passage through the circumplanetary disk (Fujita et al., 2013; D’Angelo & Podolak, 2015; Suetsugu & Ohtsuki, 2016). For example, in Suetsugu & Ohtsuki (2016), we examined the effect of ablation during capture of planetesimals by a waning circumplanetary gas disk. We found that the effect of ablation in the case of planetesimals captured into planet-centered orbits with large semi-major axes ($\gtrsim 0.1R_{\rm H}$) is insignificant, because the mass loss due to ablation mostly occurs during a short period of time when the bodies pass through the dense part of the circumplanetary disk near their pericenter. On the other hand, when their semi-major axes become sufficiently small due to orbital decay by gas drag, the captured planetesimals constantly pass through the dense part of the disk, and their mass loss would be significant. Further studies of the influence of ablation of captured planetesimals on satellite accretion would also be desirable.

We thank Takayuki Tanigawa and Toshihiko Kadono for discussion, and the anonymous reviewer for useful comments on the manuscript. This work was supported by JSPS Grants-in-Aid for JSPS Fellows (12J01826) and Scientific Research B (22340125 and 15H03716). Part of numerical calculations were performed using computer systems at the National Astronomical Observatory of Japan.

APPENDIX A.

Gas Drag Coefficient for Small Planetesimals

In the present work, we examine cases with $r_{\rm s}=1-10^{4}$ m. The assumption of $C_{\rm D}=1$ is reasonable for large planetesimals, but the dependence of $C_{\rm D}$ on the Mach number and the Reynolds number would be important for small bodies. In order to check the validity of the above assumption, we examined values of the drag coefficient of planetesimals orbiting in the circumplanetary gas disk based on the gas-starved disk model (Canup & Ward, 2002), and taking account of its dependence on the Mach number and the Reynolds number following Tanigawa et al. (2014). For planetesimals captured into retrograde orbits, we found that $C_{\rm D}\simeq 2$ regardless of their sizes, because of their large velocity relative to the gas. For those captured into prograde orbits, $C_{\rm D}$ can become as large as $5-7$ in the outer part of the disk with $r/R_{\rm H}\gtrsim 0.1$, but the values become smaller and close to unity in the inner part of the disk; even in the case of $r_{\rm s}=1$m, we found $1\lesssim C_{\rm D}\lesssim 3$ at $r/R_{\rm H}\lesssim 0.02$, which corresponds to the current radial locations of the Galilean satellites. Therefore, we think that the above assumption is reasonable for the size range considered in the present work, although a more exact treatment is required for still smaller bodies.

APPENDIX B.

Orbital Evolution After Capture

In the case of capture of planetesimals initially on heliocentric circular orbits into planet-centered prograde orbits shown in Figure 3, the number of captured planetesimals rapidly increases for $0\lesssim t/T_{\rm K}\lesssim 5$, and its growth slows down for $5\lesssim t/T_{\rm K}\lesssim 35$. This can be explained by the difference in the manner of orbital decay after capture. Figures 4 show orbital behavior of planetesimals captured from heliocentric circular orbits when $\zeta=10^{-6}$, for three cases with slightly different initial conditions. The red and green lines show the case of rapid orbital decay after permanent capture, while the blue line shows the case of slow decay. The time evolution of several quantities for these orbits are shown in Figure 15. In the case of long-lived orbits shown by the blue line (Figures 4, 15(a) and 15(b)), first, semi-major axis of the planet-centered orbit $a_{\rm p}$ decreases rapidly, then the body experiences a phase of rather slow orbital decay before spiraling into the planet. This is because the body maintains a rather large pericenter distance after permanently captured, in contrast to the case of the other short-lived orbits (Figure 15(c)). Figure 15(d) shows the evolution of the Stokes number $(St\equiv u/(a_{\rm drag}\Omega_{\rm p}^{-1}))$. The evolution of the Stokes number strongly depends on the eccentricity of the planet-centered orbit, $e_{\rm p}$. When $e_{\rm p}$ is large, the radial region swept by the planetesimal is very wide, and the relative velocity between the planetesimal and the gas significantly changes depending on the distance from the planet. This causes the large variation of the Stokes number in Figure 15(d). After $e_{\rm p}$ becomes sufficiently small, the oscillation in the Stokes number ceases because the radial distance from the planet becomes nearly constant.

APPENDIX C.

A Simple Model for the Circumplanetary Disk at the Last Stage of Giant Planet Formation

Hydrodynamic simulations of gas accretion onto a forming giant planet show the formation of circumplanetary disk, whose gas density distribution is nearly axisymmetric in the vicinity of the planet. Although detailed structure of the disk depends on the parameters adopted in simulations, here we consider a simple model on the basis of recent high resolution simulations (Machida et al., 2008; Tanigawa et al., 2012). These simulations show that the size of the circumplanetary disk is given by $r_{\rm d}\simeq 0.2R_{\rm H}$ with $R_{\rm H}$ being the planet’s Hill radius, and the gas surface density at this disk edge is approximately 10 times enhanced compared to the unperturbed surface density of the protoplanetary disk, $\Sigma_{\rm ppd,gas}$. We assume that $\Sigma_{\rm ppd,gas}$ at this last stage of giant planet formation is depleted by a factor $f_{\rm dep}$ compared to the gas surface density of the minimum mass solar nebula model, $\Sigma_{\rm MMSN}$; $\Sigma_{\rm ppd,gas}=f_{\rm dep}\Sigma_{\rm MMSN}$. On the other hand, the typical value of the gas surface density at the outer edge of the circumplanetary disk in the gas-starved disks model is about 1 gcm*-2* (Canup & Ward, 2002), while the gas surface density at Jupiter’s orbit is on the order of 102 gcm*-2* in the minimum mass solar nebula model. Therefore, if $f_{\rm dep}\simeq 10^{-3}$, the gas surface density of the protoplanetary disk at Jupiter’s orbit at the time of satellite accretion is $\Sigma_{\rm ppd,gas}\simeq 10^{-1}$ gcm*-2*, which gives the surface density of the circumplanetary disk roughly consistent with the gas-starved disk model.