Capture of Solids by Growing Proto-gas Giants: Effects of Gap Formation and Supply-limited Growth

TL;DR

This study examines how gap formation in circumstellar discs affects planetesimal capture during gas giant formation, revealing that gap effects both enhance and limit solid accretion, impacting heavy element enrichment.

Contribution

It introduces dynamical simulations of planetesimal capture considering gap formation, showing its dual role in increasing and decreasing accretion efficiency.

Findings

Gap formation reduces gas surface density and limits protoplanet growth.

Capture of planetesimals is at most ~3 Earth masses during gas accretion.

High initial solid surface density is required to explain Jupiter's heavy elements.

Abstract

Studies of internal structure of gas giant planets suggest that their envelopes are enriched with heavier elements than hydrogen and helium relative to their central stars. Such enrichment likely occurred by solid accretion during late formation stages of gas giant planets in which gas accretion dominates protoplanetary growth. Some previous studies performed orbital integration of planetesimals around a growing protoplanet with the assumption of an uniform circumstellar disc to investigate how efficiently the protoplanet captures planetesimals. However, not only planetesimals but also disc gas are gravitationally perturbed by the protoplanet in its late formation stages, resulting in gap opening in the circumstellar disc. In this study, we investigate the effects of such gap formation on the capture of planetesimals by performing dynamical simulations of planetesimals around a growing proto-gas giant planet. Gap formation reduces the surface density of disc gas, makes a steep pressure gradient and limits the growth rate of the protoplanet. We find that the first effect enhances the capture of planetesimals, while the others reduce it. Consequently the amount of planetesimals captured during the gas accretion is estimated to be at most ~3 Mearth. We conclude that the in-situ capture of planetesimals needs the initial solid surface density more than five times higher than that of the minimum mass solar nebula for explaining the inferred large amount of heavy element in Jupiter. For highly dense warm Jupiters, we would need additional processes enhancing the capture and/or supply of planetesimals.

Figures

Tables

Table 1. Table 1: Model settings.

Model	$f_{\mathrm{disc}}$	${\dot{M}}_p$	${\eta }_{\mathrm{gas}}$	$({\alpha }_{\nu },{\tau }_{\mathrm{disp}})$	Free Parameters	used in
$\mathrm{model}A$	1	constant	Eq. (9)	-	$R_{\mathrm{pl}},{\dot{M}}_p$	Sec. 3
$\mathrm{model}\mathrm{B1}$	Eq. (37)	constant	0	-	$R_{\mathrm{pl}},{\dot{M}}_p,{\alpha }_{\nu }$	Sec. 4
$\mathrm{model}\mathrm{B2}$	Eq. (37)	constant	Eq. (9)	-	$R_{\mathrm{pl}},{\dot{M}}_p,{\alpha }_{\nu }$	Sec. 4
$\mathrm{model}\mathrm{C0}$	Eq. (37)	Eq. (52)	Eq. (9)	$(\infty ,\infty )$	$R_{\mathrm{pl}}$	Sec. 5
$\mathrm{model}\mathrm{C1}$	Eq. (37)	Eq. (52)	Eq. (9)	$({10}^{-2},{10}^3)$	$R_{\mathrm{pl}}$	Sec. 5
$\mathrm{model}\mathrm{C2}$	Eq. (37)	Eq. (52)	Eq. (9)	$({10}^{-3},{10}^4)$	$R_{\mathrm{pl}}$	Sec. 5
$\mathrm{model}\mathrm{C3}$	Eq. (37)	Eq. (52)	Eq. (9)	$({10}^{-4},{10}^5)$	$R_{\mathrm{pl}}$	Sec. 5

Equations

$$\frac{{\mathrm{d}}^2{\mathbf{r}}_{\mathrm{pl}}}{\mathrm{d}{t}^2}=-G\frac{M_{\mathrm{s}}}{r_{\mathrm{pl},\mathrm{s}}^3}{\mathbf{r}}_{\mathrm{pl},\mathrm{s}}-G\frac{M_{\mathrm{p}}}{r_{\mathrm{pl},\mathrm{p}}^3}{\mathbf{r}}_{\mathrm{pl},\mathrm{p}}+{\mathbf{f}}_{\mathrm{gas}},$$

$${\mathbf{f}}_{\mathrm{gas}}=-\frac{1}{2m_{\mathrm{pl}}}{C}_{\mathrm{D}}\pi {R_{\mathrm{pl}}}^2{\rho }_{\mathrm{gas}}u\mathbf{u}.$$

$$C_{\mathrm{D}}\simeq {\left[{\left(\frac{24}{\mathcal{R}}+\frac{40}{10+\mathcal{R}}\right)}^{-1}+\frac{3\mathcal{M}}{8}\right]}^{-1}+\frac{(2-\omega )\mathcal{M}}{1+\mathcal{M}}+\omega ,$$

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

$$f_{\mathrm{disc}}(r_{\mathrm{pl}})={\Sigma }_{\mathrm{gas}}(r_{\mathrm{pl}})/{\Sigma }_0,$$

$$h_{\mathrm{g}}\equiv \frac{c_{\mathrm{s}}}{{\Omega }_{\mathrm{K}}},$$

$$v_{\mathrm{gas}}=v_{\mathrm{K}}\left(1-{\eta }_{\mathrm{gas}}\right)$$

$${\eta }_{\mathrm{gas}}$$

$${\alpha }_{\mathrm{disk}}\equiv -\frac{\mathrm{d}\ln {\Sigma }_{\mathrm{gas}}}{\mathrm{d}\ln r},{\beta }_{\mathrm{disk}}\equiv -\frac{\mathrm{d}\ln c_{\mathrm{s}}}{\mathrm{d}\ln r}.$$

$$\Delta t={\eta }_{\mathrm{ts}}\sqrt{\frac{|\mathbf{a}||{\mathbf{a}}^{(2)}|+|\dot{\mathbf{a}}{|}^2}{|\dot{\mathbf{a}}||{\mathbf{a}}^{(3)}|+|{\mathbf{a}}^{(2)}{|}^2}},$$

$$R_{\mathrm{p}}={\left(\frac{3M_{\mathrm{p}}}{4\pi {\rho }_{\mathrm{ave}}}\right)}^{1/3},$$

$$\frac{j}{N_{\mathrm{pl}}}=\frac{1}{M_{\mathrm{solid}}}{\int }_{a_{\mathrm{in}}}^{a_j}{\Sigma }_{\mathrm{pl}}\cdot 2\pi a\mathrm{d}a,$$

$$f_{\mathrm{e},\mathrm{i}}(e,i)=4\frac{{\Sigma }_{\mathrm{pl}}}{m_{\mathrm{pl}}}\frac{ei}{⟨e^2⟩⟨i^2⟩}\exp \left[-\frac{e^2}{⟨e^2⟩}-\frac{i^2}{⟨i^2⟩}\right],$$

$${\dot{M}}_{\mathrm{cap}}=m_{\mathrm{pl}}{n}_{\mathrm{pl}}{\Gamma }_{\mathrm{col}}{v}_{\mathrm{rel}},$$

$$R_{\mathrm{cap}}=R_{\mathrm{p}}{\left\{1+{\left(\frac{v_{\mathrm{esc}}}{v_{\mathrm{rel}}}\right)}^2\right\}}^{1/2},$$

$$\frac{M_{\mathrm{p}}}{M_{\mathrm{s}}}>\frac{1}{2}\frac{a_{\mathrm{p}}}{R_{\mathrm{p}}}{e}^2{i}^2,$$

$$\frac{M_{\mathrm{p}}}{M_{\mathrm{s}}}$$

$${\left(\frac{a_{\mathrm{p}}}{5.20\mathrm{AU}}\right)}^{3/4}{\left(\frac{e}{1.00\times 10^{-3}}\right)}^{3/2}{\left(\frac{i}{1.00\times 10^{-3}}\right)}^{3/2}.$$

$${\dot{M}}_{\mathrm{cap}}$$

$$\sim \sqrt{\frac{32\pi G{\rho }_{\mathrm{ave}}}{3}}{R_{\mathrm{p}}}^2{\Sigma }_{\mathrm{pl}}.$$

$$E_{\mathrm{J}}$$

$$U_{\mathrm{J}}$$

$$U_0=\frac{G\left(M_{\mathrm{s}}+M_{\mathrm{p}}\right)}{a_{\mathrm{p}}}\left[\frac{3}{2}+\frac{9}{2}{h}^2+O\left(h^3\right)\right].$$

$$\tilde{b}\equiv \frac{a-a_{\mathrm{p}}}{a_{\mathrm{p}}h},\tilde{e}\equiv \frac{e}{h},$$

$$h\equiv {\left(\frac{M_{\mathrm{p}}}{3M_{\mathrm{s}}}\right)}^{1/3}.$$

$$E_{\mathrm{J}}=-\frac{1}{2a}-\sqrt{a\left(1-e^2\right)}\cos i+\frac{3}{2}+\frac{9}{2}{h}^2+O\left(h^3\right)$$

$${\tilde{E}}_{\mathrm{J}}\equiv \frac{E_{\mathrm{J}}}{h^2}=\frac{1}{2}\left({\tilde{e}}^2+{\tilde{i}}^2\right)-\frac{3}{8}{\tilde{b}}^2+\frac{9}{2}+O\left(b,e,i,h\right),$$

$$\frac{\mathrm{d}{\tilde{E}}_{\mathrm{J}}}{\mathrm{d}t}=\frac{1}{h^2}\frac{\mathrm{d}{E}_{\mathrm{J}}}{\mathrm{d}t}-\frac{2}{3}\frac{E_{\mathrm{J}}}{h^2}\frac{1}{{\tau }_{\mathrm{grow}}},$$

$$\frac{\mathrm{d}{E}_{\mathrm{J}}}{\mathrm{d}t}$$

$$+\frac{1}{e}\frac{\mathrm{d}e}{\mathrm{d}t}\frac{e^2}{\sqrt{1-e^2}}\cos i.$$

$$+\frac{1}{i}\frac{\mathrm{d}i}{\mathrm{d}t}\sqrt{1-e^2}i\sin i]$$

Full text

• planets and satellites: composition
• planets and satellites: formation
• planets and satellites: gaseous planets
• planet-disc interactions
• protoplanetary discs

††pubyear: 2019††pagerange: Capture of Solids by Growing Proto-gas Giants : Effects of Gap Formation and Supply-limited Growth–References

1 Introduction

For understanding of the formation of gas giant planets, planetary internal composition provides crucial information. Many studies of the internal structure of Jupiter and Saturn indicate that the gas giant planets contain much more heavy elements (or “metals”) in their interiors, relative to the solar composition (Guillot & Gautier, 2014, references therein). Indeed, using the updated gravitational moments observed by the Juno spacecraft, Wahl et al. (2017) recently estimated that the total mass of heavy elements contained inside Jupiter is as much as 24–46 $M_{\mathrm{\oplus}}$, while Jupiter would contain only several $M_{\mathrm{\oplus}}$ if its bulk composition were the solar. Saturn’s heavy-element content is also estimated to be as large as 13–28 $M_{\mathrm{\oplus}}$ (Saumon & Guillot, 2004). Likewise, for exoplanets, internal structure modeling of warm Jupiters based on their mass-radius relationships suggests that they contain greater amounts of heavy elements ($\sim 10$–100 $M_{\mathrm{\oplus}}$) than their host stars (Guillot et al., 2006; Miller & Fortney, 2011; Thorngren et al., 2016). Importantly, although the core-envelope partitioning of heavy elements is uncertain, not all the heavy elements are contained in the core, but the envelope is highly enriched with heavy elements in most of the internal structure models consistent with observations.

From a viewpoint of planetary accretion, it is also difficult to form cores massive enough to account for the inferred amounts of heavy elements ($\gtrsim 20M_{\mathrm{\oplus}}$), meaning that large fractions of the heavy elements are contained in the envelopes. In the core accretion model for planet formation, once a solid core grows to a critical mass, runaway gas accretion occurs to form a massive envelope (Mizuno, 1980; Bodenheimer & Pollack, 1986). In such a dense envelope, planetesimals are vaporized before reaching the core (Podolak et al., 1988; Pollack et al., 1996) and the heavy elements are distributed through the envelope. Although the critical core mass increases with core accretion rate, formation of such massive cores requires unrealistically high rates of core accretion (Ikoma et al., 2000; Ikoma et al., 2006). Conversely, if envelopes are highly enriched with heavy elements, the critical core mass is greatly reduced (Hori & Ikoma, 2011; Venturini et al., 2015).

Hence, the enrichment of heavy elements inside gas giant planets is considered to occur in gas accretion phases after the critical core mass is reached. A protoplanet grows in mass via gas accretion and, thus, expands its feeding zone. This results basically in further accretion of planetesimals, provided dynamical effects of the protoplanet on planetesimals are ignored (Pollack et al., 1996). However, in such gas accretion stages, since the protoplanet’s mass is usually comparable with the total mass of the surrounding planetesimals, the protoplanet has a non-negligible, dynamical effect on the planetesimals.

Dynamical capture process of planetesimals by a growing proto-gas giant planet was first studied by Zhou & Lin (2007). They performed numerical integration of orbital changes of planetesimals (massless test-particles) around a protoplanet growing at an assumed rate in a gaseous protoplanetary disc. They demonstrated that planetesimals are stirred gravitationally by the protoplanet and, then, change their orbits (i.e., migrate) because of disc-gas drag. Consequently, many of the planetesimals go out of the protoplanet’s feeding zone without being engulfed, ending up with a gap around the protoplanet. Also, some of them exterior to the outer boundary of the feeding zone are trapped at mean motion resonances to the protoplanet. Thus, the capture rate of planetesimals depends on the expansion of the protoplanet’s feeding zone and the excitation and damping of velocity dispersion of planetesimals by gravitational stirring and gas drag, respectively. Shiraishi & Ida (2008) conducted a systematic investigation of such capture processes by similar simulations to Zhou & Lin (2007) and gave detailed interpretation to the capture mechanisms. In addition, they derived semi-analytical formulae for the capture rate of planetesimals, which reproduced their simulation results well (also see Section 3.3).

Such a hypothesis of dynamical capture of planetesimals during the gas accretion phases may certainly explain the observed excess of heavy elements qualitatively, but not quantitatively. According to the above two studies, the total mass of planetesimals that the protoplanet captures till growing to the Jupiter mass is as small as $\sim$4–8 $M_{\mathrm{\oplus}}$ even in the case of a relatively massive circumstellar disc being about three times heavier than the minimum-mass solar nebula (MMSN, Hayashi, 1981). Such heavy element contents of the envelope are obviously lower than those inferred from recent internal structure modeling. Also, regarding the metallicities of warm gas giants, Hasegawa et al. (2018) recently showed that not their absolute values but the correlation with host stars’ metallicities are consistent with theoretical prediction derived by use of the semi-analytical formulae from Shiraishi & Ida (2008).

Those previous studies showed gas drag for damping the velocity dispersion of planetesimals has a crucial effect on the dynamical capture of planetesimals by a protoplanet. However, they simply assumed that the gaseous disc was uniform, namely, the gas surface density was spatially constant. In reality, protoplanets with mass of several tens of $M_{\mathrm{\oplus}}$ have gravitational effects not only on planetesimals, but also on disc gas. One of the important effects is the formation of a gap near the protoplanet’ orbit in the circumstellar gas disc (Lin & Papaloizou, 1979; Goldreich & Tremaine, 1980). Since gap formation brings about changes in the efficiency of damping of planetesimals’ random velocities by disc gas drag and in the rate of the protoplanet’s growth, it would affect the capture of planetesimals. In some cases, the gap becomes wider than the feeding zone and the surface density in the feeding zone becomes lower by several orders of magnitude than in the unperturbed disc (Kanagawa et al., 2016; Kanagawa et al., 2017). This is important partly because the thinner the disc gas is, the higher the capture rate of planetesimals (see Zhou & Lin, 2007; Shiraishi & Ida, 2008) and partly because the modified gas density profile affects the direction and velocity of radial drift of planetesimals (e.g. Adachi et al., 1976). In addition, once a gap is opened in the circumstellar disc, the growth of the protoplanet slows down (Tanigawa & Ikoma, 2007; Tanigawa & Tanaka, 2016), which, conversely, reduces the capture rate of planetesimals (Shiraishi & Ida, 2008). Thus, it is difficult to know in advance whether and how the capture rate of planetesimals increases or decreases by the effect of gap opening.

The main purpose of this study is to quantify the effects of the gap formation on the capture of planetesimals by a growing proto-gas giant planet. In Section 2, we describe the basic model and settings used in this study.

We investigate the capture of planetesimals step by step. In Section 3, we consider a protoplanet growing at a constant rate in an unperturbed circumstellar disc. We review the analytical interpretation for the capture of planetesimals given by Shiraishi & Ida (2008), because that is helpful for understanding our numerical results. Then, we compare results of our numerical integration with their analytical formulae for the capture rate of planetesimals.

After that, we consider the effects of the presence of a gap in the circumstellar disc in Section 4 and, furthermore, add the effects of time-variant rates of protoplanet growth regulated by disc gas supply in Section 5. In Section 6, we make a brief discussion of other effects that we have ignored. Finally we summarise and conclude this study in Section 7.

2 Basic Method and Model

2.1 Model Settings

In order to investigate the capture of planetesimals by a growing protoplanet, we perform orbital integration of a protoplanet of mass $M_{{\rm p}}$ and 10,000 planetesimals of mass $m_{{\rm pl}}$ around a central star of mass $M_{{\rm s}}$. The protoplanet is assumed to grow only via gas accretion with a rate $\dot{M}_{{\rm p}}$, which means that planetesimal accretion makes no contribution to protoplanetary growth. Planetesimals are treated as test particles, their mutual gravitational interaction being ignored, and are subject to gas drag in the circumstellar (or protoplanetary) disc. The protoplanet and central star orbit their common centre of mass; their orbits are assumed to be circular. Here we set that $M_{{\rm s}}$ = 1 $M_{\mathrm{\odot}}$ and the semi-major axis of the protoplanet $a_{{\rm p}}$ is 5.2 AU.

The equation of motion for a planetesimal is

[TABLE]

where ${\bf r_{{\rm pl}}}$, ${\bf r_{{\rm pl,s}}}$ and ${\bf r_{{\rm pl,p}}}$ are the position vectors of the planetesimal relative to the mass centre, relative to the star, and relative to the protoplanet, respectively, $r_{k}=|{\bf r}_{k}|$ ($k$ = "pl", "pl,s", and "pl,p"), $t$ is the time, $G$ is the gravitational constant, and $\bf f_{\rm gas}$ is the gas drag given by

[TABLE]

Here ${\bf u}$ is the planetesimal velocity relative to the ambient gas, $C_{{\rm D}}$ is the non-dimensional drag coefficient, $\rho_{{\rm gas}}$ is the gas density, $R_{{\rm pl}}$ is the planetesimal radius, and $u=|{\bf u}|$. In general, $C_{{\rm D}}$ is a function of the Reynolds number $\mathcal{R}$ (= $2\rho_{{\rm gas}}R_{{\rm pl}}u/\mu$; $\mu$ being dynamic viscosity) and the Mach number $\mathcal{M}$ (= $v_{{\rm pl}}/c_{{\rm s}}$; $c_{{\rm s}}$ being the sound speed). The dynamic viscosity is given by $\mu=(1/3)\rho_{{\rm gas}}c_{{\rm s}}l_{\rm p}$, where $l_{\rm p}$ is the mean free path of gas molecules for which we adopt the collision cross section of hydrogen molecules (= $2\times 10^{-15}{\mathrm{cm}}^{2}$). In this study, we use an approximated formulae for $C_{{\rm D}}$ written as (e.g. Tanigawa et al., 2014)

[TABLE]

where $\omega$ is a correction factor, the value of which is 0.4 for $\mathcal{R}<2\times 10^{5}$ and 0.2 for $\mathcal{R}>2\times 10^{5}$.

The circumstellar disc is assumed to be vertically isothermal and, thus, the gas density $\rho_{{\rm gas}}$ is expressed with the surface density $\Sigma_{\rm gas}$ and the scale height of disc gas $h_{{\rm g}}$ as a function of the height from the mid-plane $z$ as

[TABLE]

For convenience we define the radial dependence of surface density $f_{\rm disc}(r_{{\rm pl}})$ as

[TABLE]

where $\Sigma_{\rm 0}$ is a reference surface density. The scale height is defined as

[TABLE]

where $\Omega_{\rm K}$ is the Keplerian angular velocity.

The gas in the circumstellar disc rotates with a velocity that differs from the Keplerian velocity $v_{{\rm K}}$, because of pressure gradient; namely

[TABLE]

with $\eta_{{\rm gas}}$ defined as

[TABLE]

where $r$ is the radial distance from the central star and $P_{{\rm gas}}$ is the gas pressure. For deriving the above equation, we have assumed $\eta_{{\rm gas}}\ll 1$ and used the ideal-gas relation $c_{{\rm s}}^{2}$ = $P_{{\rm gas}}/\rho_{{\rm gas}}$. Using this relation for $c_{{\rm s}}$ and Eq. (4) for $\Sigma_{\rm gas}$, $\eta_{{\rm gas}}$ can be written as

[TABLE]

where

[TABLE]

In all the simulations below, we set $\Sigma_{\rm 0}=1.43~{}\times~{}10^{2}~{}\mathrm{g}~{}{\mathrm{cm}}^{-3}$ and $c_{{\rm s}}=6.52~{}\times~{}10^{4}~{}\mathrm{cm}~{}{\rm s}^{-1}$, which correspond to the values at 5.2 ${\rm AU}$ in MMSN (Hayashi, 1981). Then, the aspect ratio of the circumstellar gas disc at the protoplanet orbit comes out to be $h_{{\rm g}}/r_{{\rm p}}\simeq 0.05$.

2.2 Numerical integration of planetesimal orbits

We integrate the equation of motion, Eq. (1), for 10 000 ($\equiv N_{\rm pl}$) planetesimals numerically, using the forth-order-Hermite integration scheme (Makino & Aarseth, 1992). For timesteps $\Delta t$, we adopt the method of Aarseth (1985), namely,

[TABLE]

where ${\bf a}$ is the acceleration of the planetesimal (= $\ddot{\bf r}_{\mathrm{pl}}$), ${\bf a}^{(k)}$ is the $k$th derivative of ${\bf a}$, and $\eta_{{\rm ts}}$ is an accuracy-controlling parameter. We have performed a benchmark test of the numerical code that we have newly developed in this study (see Appendix A). In the benchmark test, we have confirmed that $\eta_{{\rm ts}}=0.01$ yields sufficiently small energy error and, thus, use the value in all the numerical calculations presented below.

We judge that the planetesimal is trapped in the protoplanet, once (i) the planetesimal enters the protoplanet’s envelope or (ii) the Jacobi energy, $E_{{\rm J}}$, is negative in the Hill sphere (see Section LABEL:sec:Orb for the definition of $E_{{\rm J}}$), and, then, stop the orbital integration for the planetesimal. We set the physical radius of the protoplanet $R_{{\rm p}}$ as

[TABLE]

where $\rho_{{\rm ave}}$ is the protoplanet’s mean density. In our simulations, we set $\rho_{{\rm ave}}=0.125~{}\mathrm{g}\,\mathrm{cm}^{-3}$, which gives twice Jupiter’s radius for Jupiter’s mass. We will discuss the validity of this setting later.

2.3 Initial distribution of planetesimals

As described above, the protoplanet is assumed to grow only via gas accretion from its initial mass $M_{{\rm p,int}}$ to final mass $M_{{\rm p,fnl}}$. We set $M_{{\rm p,int}}=8.0\times 10^{-5}M_{\mathrm{\odot}}$ $(\lesssim 30~{}M_{\mathrm{\oplus}})$ and $M_{{\rm p,fnl}}=1.0\times 10^{-3}M_{\mathrm{\odot}}$ $(\gtrsim 1M_{\mathrm{J}})$. At the beginning of simulations, we assume that planetesimals that existed in the feeding zone of the protoplanet of mass $M_{{\rm p,int}}$

, inside which planetesimals have a possibility to be captured by the protoplanet (see Section 3.1 for the exact definition),

have already accreted in the protoplanet. Thus, we distribute planetesimals initially between the semi-major axes $a$ = $3.6~{}{\rm AU}$ ($\equiv a_{\rm in}$) and $6.8~{}{\rm AU}$ ($\equiv a_{\rm out}$), which covers the feeding zone of the protoplanet of mass $M_{{\rm p,fnl}}$, except for the zone between $4.7~{}{\rm AU}$ ($\equiv a_{\rm in,0}$) to $5.7~{}{\rm AU}$ ($\equiv a_{\rm out,0}$) (equivalent to the initial feeding zone). The initial distribution is visualized in Fig. 1$a$.

We set the initial values of the Keplerian orbital elements of planetesimals, semi-major axis $a$, eccentricity $e$, inclination $i$, longitude of ascending node $\Omega_{\rm loa}$, longitude of pericenter $\varpi$ and true anomaly $\psi$ in the following way. The initial total mass of planetesimals is equal to that of solids in the zones $a_{\rm in}\leq a\leq a_{\rm in,0}$ and $a_{\rm out,0}\leq a\leq a_{\rm out}$ in the minimum-mass solar nebula (MMSN; Hayashi, 1981), namely $6.7M_{\mathrm{\oplus}}$ ($\equiv M_{\rm solid}$). Note that a particle that we call a planetesimal in this study is, in fact, not a real planetesimal, but an imaginary super-particle that represents a group of many identical planetesimals. When simulating its orbital motion, we regard the particle as a real planetesimal of radius $R_{{\rm pl}}$ and mass $m_{{\rm pl}}$, while we regard it as a super-particle when estimating the total mass of captured planetesimals. The super-particles have the same mass $M_{\rm solid}/N_{\rm pl}$ ($N_{\rm pl}$ being the total number of the super-particles = 10000) and the same radius $R_{{\rm pl}}$. The semi-major axis of the particle labeled as $j$ ($j$ = 1, $\cdots$, $N_{\rm pl}$), $a_{j}$, is given in such a way that

[TABLE]

so that the mass distribution is consistent with the surface density of solids $\Sigma_{\rm pl}$, which is assumed to be spatially uniform in this study.

As for eccentricity and inclination, assuming that planetesimals are scattered by their mutual gravitational interaction, we adopt the Rayleigh distribution given by (Lissauer, 1993)

[TABLE]

as the initial $e$ and $i$ of planetesimals. We set ${\left\langle{e}^{2}\right\rangle}^{1/2}=0.001$ and ${\left\langle{i}^{2}\right\rangle}^{1/2}=0.0005$; however, choice of $\left\langle e^{2}\right\rangle$ and $\left\langle i^{2}\right\rangle$ scarcely affects our conclusions because viscous stirring by the protoplanet changes the distribution immediately. The orbital angles $\Omega_{\rm loa}$, $\varpi$ and $\psi$ are distributed uniformly. Finally, we transform these orbital elements into the Cartesian coordinate and then perform orbital integration.

2.4 Summary of case studies

The main focus of this study is on the effects of gap formation on capture of planetesimals by a growing protoplanet. We investigate such effects in a step-by-step fashion. The free parameters include $f_{\rm disc}$, $\dot{M}_{{\rm p}}$, and $R_{{\rm pl}}$. Firstly, we investigate the case where $f_{\rm disc}$ = 1 (i.e., no gap) and $\dot{M}_{{\rm p}}$ is constant with time. This model is regarded as the standard one and referred to as model A(${\rm Section~{}\ref{sec:Model1}}$). Secondly, we consider the effects of gap formation with constant $\dot{M}_{{\rm p}}$ (model B; ${\rm Section~{}\ref{sec:Model2}}$). Thirdly, we explore the case where the protoplanet embedded in a gapped disc grows with a mass-dependent rate (model C; ${\rm Section~{}\ref{sec:Model3}}$). For $R_{{\rm pl}}$ we consider four values, $1\times 10^{4},$ $10^{5}$, $10^{6}$, and $10^{7}\mathrm{cm}$ in each model. The cases and parameter values we consider are summarized in Table  1.

3 Capture of Planetesimals in a Uniform Circumstellar Disc

3.1 Basic Mechanism for Orbital Evolution of Planetesimals around a Growing Protoplanet

Before investigating the effects of a gap in the circumstellar disc on the capture of planetesimals by a growing protoplanet, we review the numerical results and analytical interpretation given by Shiraishi & Ida (2008), which is helpful for interpreting our numerical results. According to their analysis, the surface density of planetesimals in the so-called feeding zone is the most important factor that determines the capture rate of planetesimals during late stages of gas giant formation. Here, we show the relation between the capture rate and surface density of planetesimals, based on Shiraishi & Ida (2008). Then, we give an interpretation to the change in the surface density of planetesimals in terms of the Jacobi energy.

The total mass of planetesimals captured per unit time by the protoplanet (called simply the capture rate, hereafter) is given by

[TABLE]

where $n_{\mathrm{pl}}$ is the number density of planetesimals, $\Gamma_{\rm col}$ is the collision cross section of the protoplanet and $v_{{\rm rel}}$ is the planetesimal’s velocity relative to the protoplanet, which depends on its eccentricity and inclination. $\Gamma_{\rm col}$ depends on whether the capture radius $R_{\rm cap}$ is larger or smaller than the scale height of the planetesimal disc $h_{\rm pl}$. The former is defined as the threshold value of the impact parameter below which planetesimals can collide the protoplanet and is given with $v_{{\rm esc}}$ being the escape velocity from the protoplanet by

[TABLE]

provided the random velocity is high enough that the stellar gravity is negligible. Although not verified a priori, this condition is approximately achieved in the numerical simulations of Shiraishi & Ida (2008) and ours. For a massive protoplanet of interest in this study ($>1M_{\mathrm{\oplus}}$), $v_{{\rm esc}}>v_{{\rm rel}}$ and thus $R_{\rm cap}\sim R_{{\rm p}}v_{{\rm esc}}/v_{{\rm rel}}$ from Eq. (16). Also, $h_{\rm pl}\sim a_{{\rm p}}i$ and $v_{{\rm rel}}\sim ev_{{\rm K}}$. Thus, the condition $R_{\rm cap}>h_{\rm pl}$ is written in terms of protoplanetary mass as

[TABLE]

which is estimated to be

[TABLE]

This condition always holds true in this study, so that

the collision cross section is regarded approximately as a rectangle with height of $2h_{{\rm g}}$ and width of $2R_{\rm cap}$, namely

$\Gamma_{\rm col}\simeq 4R_{\rm cap}h_{\rm pl}$. Thus, Eq. (15) becomes

[TABLE]

This indicates that $\dot{M}_{{\rm cap}}$ depends only on $\Sigma_{\rm pl}$ for given $\rho_{{\rm ave}}$ and $R_{{\rm p}}$. Although $\rho_{{\rm ave}}$ and $R_{{\rm p}}$ also vary during protoplanetary growth, here we focus on the physical process that determines $\Sigma_{\rm pl}$.

According to Shiraishi & Ida (2008), using the surface density of planetesimals averaged over the feeding zone as $\Sigma_{\rm pl}$ in Eq. (3.1) (hereafter, called the effective surface density and denoted by $\Sigma_{{\rm eff}}$), one obtains an estimate for $\dot{M}_{{\rm cap}}$ that is consistent with their numerical results. The feeding zone is defined with the Jacobi energy $E_{{\rm J}}$ as (e.g. Hayashi et al., 1977),

[TABLE]

where $U_{{\rm J}}$ is the Jacobi potential defined as

[TABLE]

$v_{{\rm pl}}^{\prime}$ and ($x^{\prime},y^{\prime}$) are the velocity and position of the planetesimal, respectively, in the coordinate system co-rotating with the protoplanet, and $U_{0}$ is a constant determined according to the boundary condition. Here, we set $U_{0}$ such that $U_{{\rm J}}$ vanishes at the Lagrange L2 point and thus

[TABLE]

By definition, $\Sigma_{{\rm eff}}$ is controlled by the change in the Jacobi energy. Below, we give an interpretation to the orbital evolution of planetesimals around the growing protoplanet in terms of the Jacobi energy.

First, in Fig. 2, we show snapshots of the orbital evolution of planetesimals around the growing protoplanet at $M_{{\rm p}}$ = (a) 41 $M_{\mathrm{\oplus}}$ and (b) 149 $M_{\mathrm{\oplus}}$ for model A. Here, the horizontal and vertical axes are the normalized semi-major axis $\tilde{b}$ and eccentricity $\tilde{e}$, which are defined respectively as

[TABLE]

where $h$ is the reduced Hill radius defined as

[TABLE]

The growth rate $\dot{M}_{{\rm p}}$ is $1~{}\times~{}10^{{-2}}~{}M_{\mathrm{\oplus}}~{}{\mathrm{year}}^{-1}$ and the size of planetesimals $R_{{\rm pl}}$ is $1\times 10^{5}~{}\mathrm{cm}$. Planetesimals are color-coded according to their values of $E_{{\rm J}}$ and the red and blue symbols represent $E_{{\rm J}}>0$ and $<0$, respectively. The black dashed lines indicate the outer boundary of the feeding zone: Using the orbital elements and setting $M_{{\rm s}}$, $a_{{\rm p}}$ and $T_{{\rm K}}(a_{{\rm p}})/2\pi$ as the units of mass, length and time, respectively, Eq. (20) can be transformed as (e.g. Hayashi et al., 1977)

[TABLE]

or

[TABLE]

where $\tilde{i}$ is the normalized inclination ($\equiv i/h$). For $b,e,i,h\ll 1$, the shape of the feeding zone is fixed on the $\tilde{b}$–$\tilde{e}$ plane, as shown in Fig. 2.

The dynamic behaviour of planetesimals shown in Fig. 2 is outlined as follows. Outside the feeding zone, except for mean motion resonances to the protoplanet, the normalised eccentricities of planetesimals are much smaller than unity. This is because gravitational scattering by the protoplanet is so week far from the protoplanet that the gas drag keeps the planetesimals’ eccentricities small. At the mean motion resonances, the eccentricities are relatively high because of amplified gravitation scattering. Inside the feeding zone, the planetesimals experience close encounters with the protoplanet and are scattered into highly eccentric orbits. Since the gravitational scattering itself never changes the Jacobi energy, it merely distributes planetesimals along the constant Jacobi energy line and does not result in ejection from nor injection into the feeding zone.

The planetesimals go into the feeding zone through the protoplanet growth and get farther away from the feeding zone through the damping and/or orbital decay by disc gas drag. This is explained in terms of change in the Jacobi energy, which is expressed as

[TABLE]

where $\tau_{{\rm grow}}$ ($\equiv M_{{\rm p}}/\dot{M}_{{\rm p}}$) is the characteristic timescale of protoplanet growth and

[TABLE]

from Eq. (25). The characteristic timescales of change in $a$, $e$ and $i$ via gas drag are given for constant $C_{{\rm D}}$ as (Adachi et al., 1976; Kary et al., 1993)

[TABLE]

where

[TABLE]

Damping of eccentricity and inclination due to gas drag results in a decrease in the Jacobi energy. Change in the semi-major axis results in increase or decrease in the Jacobi energy depending on the disc gas flow ($\eta_{{\rm gas}}$) and orbital elements. When $C_{{\rm D}}$ is not constant, we have to integrate Eq. (27) for one orbital period directly to know the change of Jacobi energy; however, the qualitative characteristics concerning the changing rate of the Jacobi energy are almost the same.

Panel (c) of Fig. 2 shows the averaged values of ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t$ on the $\tilde{b}$–$\tilde{e}$ plane for $M_{{\rm p}}=100M_{\mathrm{\oplus}}$, which are obtained by integrating Eq. (27) for one orbital period with variable $C_{{\rm D}}$. The detailed method of the integration is shown in the Appendix B.

Planetesimals in the red area have their Jacobi energies enhanced (i.e., ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t>0$) and move toward the feeding zone, while their Jacobi energies decreased (i.e., ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t<0$) in the blue area. At the boundary (i.e., ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t=0$) indicated by the green solid line, the growth and damping timescales are similar to each other. The gray solid line shows the evolution path of a planetesimal. From this panel, one can predict how those planetesimals move around the growing protoplanet: Far from the protoplanet, since their equilibrium eccentricities are sufficiently small, the planetesimals stay in the red area and their Jacobi energies tend to increase. This means that those planetesimals are pushed toward the feeding zone(array 1). Once they enter the feeding zone ($E_{{\rm J}}>0$), their eccentricities become further higher due to strong gravitational scattering by the protoplanet(array 2). Many of the planetesimals move into the blue area and their Jacobi energies decrease. This means that the planetesimals are pushed out of the feeding zone(array 3).

The above analysis indicates that the capture of planetesimals depends on the effective surface density of planetesimals $\Sigma_{{\rm eff}}$ and also on the rate of change in the Jacobi energy, which depends on the rates of protoplanetary growth and damping of planetesimals’ velocity dispersion. The knowledge of the relation among the capture rate of planetesimals and these timescales is useful to understand the numerical results shown in succeeding sections.

3.2 Numerical Results

Here, we investigate the capture of planetesimals in a uniform circumstellar disc (model A; see Table 1). Figure 3 shows the temporal changes in cumulative mass of planetesimals captured by the protoplanet (simply the cumulative captured mass or $M_{{\rm cap}}$, hereafter) with protoplanetary growth for $\dot{M}_{{\rm p}}$ = $1\times 10^{-2}$ $M_{\mathrm{\oplus}}~{}{\mathrm{year}}^{-1}$ in panel (a) and for $\dot{M}_{{\rm p}}$ = $1\times 10^{-3}$ $M_{\mathrm{\oplus}}~{}{\mathrm{year}}^{-1}$ in panel (b). The black line represents the the gas-free case and the coloured lines represent the cases of different planetesimal size.

First, we investigate the dependence of the cumulative captured mass $M_{{\rm cap}}$ on gas drag or, equivalently, planetesimal size in the higher-$\dot{M}_{{\rm p}}$ case (Fig. 3a). Planetesimals are captured most efficiently in the gas-free case. In the case of $R_{{\rm pl}}$ = $1~{}\times~{}10^{7}~{}\mathrm{cm}$ for which gas drag is less effective, the cumulative captured mass increases quite similarly until $M_{{\rm p}}=200~{}M_{\mathrm{\oplus}}$, but more slowly for $M_{{\rm p}}>200~{}M_{\mathrm{\oplus}}$ than in the gas-free case. As planetesimal size decreases, the coloured line starts to deviate from the black line at a smaller mass, that is, $M_{{\rm p}}=130~{}M_{\mathrm{\oplus}}$ for $R_{{\rm pl}}$ = $1~{}\times~{}10^{6}~{}\mathrm{cm}$, while $M_{{\rm p}}=60~{}M_{\mathrm{\oplus}}$ for $R_{{\rm pl}}$ = $1~{}\times~{}10^{5}~{}\mathrm{cm}$. In the case of $R_{{\rm pl}}$ = $1~{}\times~{}10^{4}~{}\mathrm{cm}$, the cumulative captured mass always increases more slowly than that in the gas-free case. Thus, the total captured mass ($M_{\rm cap}^{\rm total}$: $M_{{\rm cap}}$ at $M_{{\rm p}}=M_{{\rm p,fnl}}$) decreases with planetesimal size.

In the lower-$\dot{M}_{{\rm p}}$ case (Fig. 3b), the cumulative captured mass is smaller than in the higher-$\dot{M}_{{\rm p}}$ case, except for the gas-free case where the total captured mass is independent of $\dot{M}_{{\rm p}}$. In the case of $R_{{\rm pl}}$ = $1~{}\times~{}10^{4}~{}\mathrm{cm}$, no planetesimals are captured by the protoplanet.

Such dependence can be interpreted, in terms of the Jacobi energy, as follows. In the gas-free case, since the damping timescale is infinity, one obtains

[TABLE]

from Eqs. (27) and (3.1). The above equation indicates that the Jacobi energy keeps increasing up to $\tilde{E}_{{\rm J}}=9/2$ and final Jacobi energy which a planetesimal has at $M_{{\rm p}}=M_{{\rm p,fnl}}$ is determined only by the initial and final conditions of the integration. The growth rate $\dot{M}_{{\rm p}}$ never affects the value, which means how many planetesimals enter the feeding zone depends on $M_{{\rm p,fnl}}$, but not on $\dot{M}_{{\rm p}}$. That is the reason why the total captured mass takes the same value regardless of $\dot{M}_{{\rm p}}$.

When gas drag is effective, some of the planetesimals inside the feeding zone are pushed out before being captured by the protoplanet. When the protoplanet is small in mass, its growth timescale is so short relative to the damping timescale that the planetesimals keep staying in the feeding zone and $\Sigma_{{\rm eff}}$ is almost the same as that in the gas-free case. Then, the protoplanet grows and thereby the growth timescale becomes long enough that the gas drag is effective, reducing $\Sigma_{{\rm eff}}$ and thereby making the cumulative captured mass smaller relative to the gas-drag case. How much $\Sigma_{{\rm eff}}$ is reduced depends on those two timescales, because ${\rm d}{\tilde{E}}_{\rm J}/{\rm d}t$ also depends on them. In summary, the cumulative captured mass becomes smaller for smaller planetesimals or slower protoplanet growth.

3.3 Comparison with Analytical Formulae from Shiraishi &

Ida (2008)

As described in Introduction, the previous study, Shiraishi & Ida (2008), investigated the capture of planetesimals in the unperturbed disc and derived semi-analytical formulae for the capture rate of planetesimals that explain their numerical results. Here, we compare our results with the analytical formulae from Shiraishi & Ida (2008).

They used a simplified gas drag formula,

[TABLE]

where the damping timescale $\tau_{{\rm damp}}$ was assumed to be unchanged through their simulations. Also they assumed that the disc gas rotates around the central star with the Kepler velocity, namely $\eta_{{\rm gas}}=0$. They found that the capture process can be divided into two different regimes. The first regime is such that no clear gap is formed in the swam of planetesimals that surrounds the protoplanet (called the planetesimal disc, hereafter). This regime is achieved when $\tau_{{\rm grow}}<\tau_{{\rm damp}}$, which means the protoplanet’s growth dictates the capture of planetesimals. Thus, $\dot{M}_{{\rm cap}}$ is approximately determined only by $\tau_{{\rm grow}}$. The second regime in which $\tau_{{\rm grow}}>\tau_{{\rm damp}}$ is such that a clear gap is formed in the planetesimal disc. In this regime, $\dot{M}_{{\rm cap}}$ depends not only on $\tau_{{\rm grow}}$ but also on $\tau_{{\rm damp}}$. The second regime is, therefore, affected by assumptions regarding the structure of the surrounding gas disc.

Figure 4 shows the change in $\dot{M}_{{\rm cap}}$ with the protoplanetary growth. The semi-analytical formula for the first regime is shown by the gray dashed line, and those for the second regime with different values of $\tau_{{\rm damp}}$ are shown by the gray dotted lines. Firstly, in the cases of no drag and the largest planetesimals ($R_{{\rm p}}=1\times 10^{7}$ cm), the black and red lines are almost parallel to the gray dashed one, since $\tau_{{\rm grow}}<\tau_{{\rm damp}}$ (i.e., the first regime). Secondly, the green and blue lines ($R_{{\rm p}}=1\times 10^{6}$ and $1\times 10^{5}$ cm, respectively) are initially similar to the gray dashed one. As the protoplanet grows, its growth timescale ($\propto M_{{\rm p}}$) becomes longer and the regime switches to the other one. Indeed, the green and blue lines depart from the dashed line at $M_{{\rm p}}\sim 130M_{\mathrm{\oplus}}$ and $60M_{\mathrm{\oplus}}$, respectively, and shift toward the lines for the second regime. Finally, as for the smallest planetesimals ($R_{{\rm p}}=1\times 10^{4}$ cm; magenta line), the capture process is in the second regime from the beginning on. It turns out that the numerical results are considerably different from the semi-analytical formulae in the second regime: the calculated $\dot{M}_{{\rm cap}}$ decreases with $M_{{\rm p}}$ more rapidly than that given by the semi-analytical formulae and show a jagged profile. The difference between the numerical and semi-analytical $\dot{M}_{{\rm cap}}$ is more striking for smaller planetesimals.

Obviously, such inconsistency between the numerical and semi-analytical results comes from the different treatments for gas drag. Our simulations demonstrate that both the variable damping timescale of planetesimals’ eccentricities and the sub-Kepler rotation of the disc gas (i.e., radial drift of planetesimals), which Shiraishi & Ida (2008) ignored, have great impacts on the capture rate of planetesimals. As shown in the next sections, if a gap is formed in the gas disc, both the timescales of eccentricity damping and radial drift are significantly changed and, thus, the condition of $\tau_{{\rm grow}}>\tau_{{\rm damp}}$ is readily achieved. It will turn out that the capture rate of planetesimals changes more complexly with protoplanet growth and is significantly different from prediction by the semi-analytical formulae derived in Shiraishi & Ida (2008).

4 Effect of Change in Damping Timescale of Planetesimals

Gap formation affects the surface density profile $f_{\rm disc}$ and, thereby, the rotation velocity $v_{{\rm K}}(1-\eta_{{\rm gas}})$ in a circumstellar disc. The change in surface density leads to changing the gas drag strength and the damping timescales, $\tau_{{\rm damp,}a}$, $\tau_{{\rm damp,}e}$ and $\tau_{{\rm damp,}i}$ almost equally. The change in $\eta_{{\rm gas}}$ alters the gas drag direction and speed, which results in changing $\tau_{{\rm damp,}a}$ more largely than the other timescales.

We investigate these effects on the capture of planetesimals step by step in this section. First, we describe the model for the surface density profile in the gap in Section 4.1. In Section 4.2, we investigate the capture of planetesimals using the gap model, but assuming that $\eta_{{\rm gas}}$ is zero. Then, we take into account the effect of $\eta_{{\rm gas}}$ according to the pressure gradient of the disc model in Section 4.3.

4.1 Gap Model

To take into account the presence of a gap around the protoplanet’s orbit in the circumstellar disc, we use the empirical formula for the radial density profile in the gap that Kanagawa et al. (2017) derived from their two-dimensional hydrodynamic simulations.

The width of the gap, $\Delta_{{\rm gap}}$, inside which the surface density is smaller than a threshold value, $\Sigma_{\rm th}$, is given approximately as

[TABLE]

where $r_{{\rm p}}$ is the semi-major axis of the planetary orbit, and $K^{{}^{\prime}}$ is a non-dimensional quantity given by

[TABLE]

$\alpha_{{\rm\nu}}$ is the viscosity factor introduced by Shakura & Sunyaev (1973). The radial distribution of the surface density in the circumstellar disc with the gap is given as

[TABLE]

where

[TABLE]

and $\Sigma_{{\rm min}}$ is the surface density at the gap bottom, which is given by

[TABLE]

with

[TABLE]

For the gapped-disc models, the hydrodynamic stability of the surface density profile must be checked by the well-known Rayleigh criterion (e.g. Chandrasekhar, 1961; Tanigawa & Ikoma, 2007),

[TABLE]

where $x=|r-{a_{{\rm p}}}|$ and $\Omega_{\rm p}$ is the angular velocity of the protoplanet. In the gap model given by Eq. (37), the Rayleigh criterion is violated in the vicinity of the protoplanet. The unrealistically steep gradient yields such a large difference in velocity between gas and planetesimals that the gas drag makes a considerable effect on the dynamics of the planetesimals. Therefore, we modify the gap structure so that the surface density profile is smoothly connected with what fulfills the Rayleigh condition marginally. We denote the point of connection by $x_{\rm m}$ and the gas velocity at that point by $v_{\rm gas,m}$. Thus, at $x<x_{\rm m}$, the gas velocity is given by

[TABLE]

The surface gas density is given by integrating the equation of motion,

[TABLE]

inward from $x$ = $x_{\rm m}$. When the surface density thus obtained is smaller than that of the gap bottom obtained by Eq. (41), we assume $f_{\rm disc}=\Sigma_{{\rm min}}/\Sigma_{\rm 0}$. Examples of the surface density profile are shown in Fig. 5. In this gap model, the gap width and depth depend on the protoplanet mass and disc gas viscosity.

We regard the growth rate of the protoplanet $\dot{M}_{{\rm p}}$, the size of planetesimals $R_{{\rm pl}}$ and the viscosity parameter $\alpha_{{\rm\nu}}$ as free parameters. Firstly, we focus on the strength of gas drag. The surface density profile is given by Eq. (37), but disc gas rotates around the central star with the Kepler velocity (namely, $\eta_{{\rm gas}}=0$). We call this model the zero-$\eta_{{\rm gas}}$ model (or ${\rm model~{}B1}$). Secondly, we focus on the direction of gas drag. The surface density profile is given by Eq. (37) and $\eta_{{\rm gas}}$ is calculated from Eq. (9). We call this model the nonzero-$\eta_{{\rm gas}}$ model (or ${\rm model~{}B2}$). Parameters and models are summarised in Table  1.

4.2 Effect of gas drag strength

Figure 6 shows the cumulative captured mass as a function of the protoplanet mass for $\dot{M}_{{\rm p}}$ = $1\times 10^{-2}$ $M_{\mathrm{\oplus}}{\mathrm{year}}^{-1}$ obtained in the case of the gapped disc with Kepler rotation (i.e., $\eta_{{\rm gas}}=0$; ${\rm model~{}B1}$). In panels (a) and (b), $R_{{\rm pl}}$ = $1\times 10^{5}$ and $1\times 10^{4}$ $\mathrm{cm}$, respectively. The lines are colour-coded according to the value of $\alpha_{{\rm\nu}}$ except for the black one representing the no-gap case.

In both panels, because of reduced gas drag due to gap formation, the cumulative captured mass becomes larger than in the unperturbed (no-gap) disc and larger with decreasing $\alpha_{{\rm\nu}}$. Such an effect turns out to be noticeable for $\alpha_{{\rm\nu}}\lesssim 10^{-3}$. This is because the gap becomes wider than the feeding zone for $\alpha_{{\rm\nu}}\lesssim 10^{-3}$, as found in Fig. 5. In the small-$R_{{\rm pl}}$ case shown in panel (b), where the damping timescale is shorter than in the high-$R_{{\rm pl}}$ case, a larger number of planetesimals are eliminated from the feeding zone, reducing the cumulative captured mass as a whole.

For a clear understanding of the effect of gap formation on the capture of planetesimals, we introduce an enhancement factor defined as

[TABLE]

where $M_{\rm A}^{\rm total}$ and $M_{\rm B1}^{\rm total}$ are the total captured mass in ${\rm model~{}A}$ and ${\rm model~{}B1}$, respectively. Figure 7 shows the calculated relationship between $f_{\rm enh,B1}$ and $\alpha_{{\rm\nu}}$ for two different values of $\dot{M}_{{\rm p}}$ and four different values of $R_{{\rm pl}}$. It is found that $f_{\rm enh,B1}$ ranges between 1.0 and 3.6, depending on $\alpha_{{\rm\nu}}$, $R_{{\rm pl}}$ and $\dot{M}_{{\rm p}}$. Note that for $\dot{M}_{{\rm p}}$ = $1\times 10^{-3}~{}M_{\mathrm{\oplus}}~{}{\mathrm{year}}^{-1}$ and $R_{{\rm pl}}$ = $1\times 10^{4}~{}\mathrm{cm}$, $f_{\rm enh,B1}$ is not defined because $M_{\rm A}^{\rm total}=0$.

The enhancement factor $f_{\rm enh,B1}$ increases with decreasing $\alpha_{{\rm\nu}}$ for given $R_{{\rm pl}}$ and $\dot{M}_{{\rm p}}$. Again, this is because the gap becomes wider and deeper with decreasing $\alpha_{{\rm\nu}}$, so that eccentricity damping and elimination of planetesimals become less effective. Also, $f_{\rm enh,B1}$ becomes larger with decreasing $R_{{\rm pl}}$ except for the cases of ${\alpha_{{\rm\nu}}~{}=}~{}1~{}\times 10^{-2}$. When $R_{{\rm pl}}$ is as large as $1\times 10^{7}\mathrm{cm}$, the timescale of eccentricity damping is much longer than the protoplanet’s growth timescale, so that the elimination of planetesimals from the feeding zone is ineffective even in the unperturbed disc. Thus, the effect of the gap is negligibly small and therefore $f_{\rm enh,B1}=1$. However, for $R_{{\rm pl}}\leq 1\times 10^{6}\mathrm{cm}$, eccentricity damping occurs efficiently enough that a significant fraction of planetesimals are pushed away from the feeding zone. The effect of the gap reduces such an elimination effect and thus $f_{\rm enh,B1}$ is larger for smaller $R_{{\rm pl}}$.

4.3 Effect of gas flow velocity

Figure 8 is the same as Fig. 6, but for the nonzero-$\eta_{{\rm gas}}$ model (${\rm model~{}B2}$). In this case, unlike in Fig. 6, the cumulative captured mass is smaller in the gapped disc (coloured lines) than in the unperturbed disc (black line), except for the case of $\alpha_{{\rm\nu}}$ = $1\times 10^{-5}$ (magenta line). On the gap slopes, $|\eta_{{\rm gas}}|$ takes large values, so that planetesimals drift rapidly in the directions away from the protoplanet, resulting in reduction in the number of planetesimals that enter the feeding zone and thereby in the capture rate of planetesimals. Such an effect appears basically regardless of $\alpha_{{\rm\nu}}$, because the gap extends beyond the feeding zone (see Fig. 5) and $|\eta_{{\rm gas}}|$ is sufficiently large even near the edges of the gap slopes in the nonzero-$\eta_{{\rm gas}}$ models. The reason why the cumulative captured mass for $\alpha_{{\rm\nu}}=1\times 10^{-5}$ and $\dot{M}_{{\rm p}}=1\times 10^{-2}M_{\mathrm{\oplus}}\,\mathrm{year}^{-1}$ (magenta line) is larger than that for the unperturbed disc case (black line) is that the effect of reduction in gas drag strength, which was shown in section 4.2, dominates over that of change in pressure gradient. Thus, compared to the unperturbed disc model, the total captured mass is larger (smaller) when the effect of gas drag strength is more (less) effective than that of $\eta_{{\rm gas}}$.

For a clear understanding of the effect of $\eta_{{\rm gas}}$ on the capture of planetesimals, we introduce an enhancement factor for ${\rm model~{}B2}$ defined as

[TABLE]

where $M_{\rm B2}^{\rm total}$ is the total mass of captured planetesimals in ${\rm model~{}B2}$. Figure 9 shows $f_{\rm enh,B2}$ versus $\alpha_{{\rm\nu}}$ for two different values of $\dot{M}_{{\rm p}}$ and four different values of $R_{{\rm pl}}$. The results are divided into three groups: In the first group including the cases of ($R_{{\rm pl}}$ [$\mathrm{cm}$], $\dot{M}_{{\rm p}}$ [$M_{\mathrm{\oplus}}{\mathrm{year}}^{-1}$]) = ($1\times 10^{7}$, $1\times 10^{-2}$), ($1\times 10^{6}$, $1\times 10^{-2}$) and ($1\times 10^{7}$, $1\times 10^{-3}$), $f_{\rm enh,B2}$ is almost unity. This is because the timescale of change in semi-major axis $\tau_{{\rm damp,}a}$ is much longer than the growth timescale of the protoplanet $\tau_{{\rm grow}}$, meaning that the radial drift on the gap slopes is too slow to affect the cumulative captured mass $M_{{\rm cap}}$. In the second group including the cases of $(R_{{\rm pl}},\dot{M}_{{\rm p}})=(1\times 10^{4},1\times 10^{-2})$ and $(1\times 10^{5},1\times 10^{-3})$, $f_{\rm enh,B2}$ is significantly smaller than unity because $\tau_{{\rm damp,}a}\ll\tau_{{\rm grow}}$ on the gap slopes. In the last group including the cases of $(R_{{\rm pl}},\dot{M}_{{\rm p}})=(1\times 10^{5},1\times 10^{-2})$ and $(1\times 10^{6},1\times 10^{-3})$, $\tau_{{\rm damp,}a}\sim\tau_{{\rm grow}}$ on the gap slopes and thus $f_{\rm enh,B2}$ takes intermediate values between zero and unity.

4.4 Change of Jacobi Energy

As indicated in Eq. (3.1), the capture rate of planetesimals $\dot{M}_{{\rm cap}}$ is proportional to the effective surface density of planetesimals inside the feeding zone $\Sigma_{{\rm eff}}$ (see Section 3.1 for its definition), which depends on the Jacobi energy $E_{{\rm J}}$. Here, we give an interpretation to numerical results obtained in Sections 4.2 and 4.3, in terms of the Jacobi energy explained in Section 3.1.

First, we focus on the effect of decrease in gas drag strength. Figure 10 shows the lines of ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t=0$, below which ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t>0$ (also see Fig. 2c), for different choices of $\alpha_{{\rm\nu}}$.

Values of ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t$ are calculated by numerically integrating Eq. (27) with $i=0$ and variable $C_{D}$ for an orbital period.

Panel (a) represents the zero-$\eta_{{\rm gas}}$ disc models (${\rm model~{}B1}$) with $M_{{\rm p}}=100M_{\mathrm{\oplus}}$, $\dot{M}_{{\rm p}}$ = $1\times 10^{-2}~{}M_{\mathrm{\oplus}}~{}{\mathrm{year}}^{-1}$ and $R_{{\rm pl}}=1\times 10^{4}~{}\mathrm{cm}$; the black solid line represents the unperturbed disc model. In every case, planetesimals with low $\tilde{e}$ ($\lesssim 1$) have positive values of ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t$ outside the feeding zone, which means that they take a similar evolution path shown in Fig. 2c until they enter the feeding zone. After entering the feeding zone, however, planetesimals take different paths depending on the gap size. As illustrated in Fig. 10a, the area of ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t>0$ expands with decreasing $\alpha_{{\rm\nu}}$ inside the feeding zone. Because of the expansion of the area of ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t>0$ and $E_{{\rm J}}>0$, the planetesimals can stay inside the feeding zone for a longer time and $\Sigma_{{\rm eff}}$ becomes higher in a deeper and wider gap.

Second, we focus on the effect of the change in gas flow direction and velocity. Panel (b) of Fig. 10 is the same as panel (a), but for the nonzero-$\eta_{{\rm gas}}$ disc model (${\rm model~{}B2}$). The values of $\dot{M}_{{\rm p}}$ and $R_{{\rm pl}}$ adopted here are different from those in panel (a) and are $1\times 10^{-3}~{}M_{\mathrm{\oplus}}~{}{\mathrm{year}}^{-1}$ and $1\times 10^{5}~{}{\mathrm{cm}}$, respectively. The change of gas flow inside the gap turns out to make a great change to the line of ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t=0$. In the three cases of $\alpha_{{\rm\nu}}\leq 1\times 10^{-3}$ (green, blue and magenta lines), planetesimals are unable to reach the feeding zone, because the evolution path is truncated by the line of ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t=0$ (solid line) before it reaches the line of $E_{{\rm J}}=0$ (dashed line). Even when the domain of ${\rm d}\tilde{E}_{{\rm J}}/{\rm d}t>0$ expands in the feeding zone, $\Sigma_{{\rm eff}}$ never increases due to the prohibited supply of planetesimals. In the case of $\alpha_{{\rm\nu}}=1\times 10^{-2}$, planetesimals barely reach the feeding zone when $M_{{\rm p}}\lesssim 100M_{\mathrm{\oplus}}$, but are never captured when $M_{{\rm p}}\gtrsim 100M_{\mathrm{\oplus}}$.

5 Effect of Change in Growth Timescale of Protoplanet

In the previous sections, we have assumed a constant growth rate of the protoplanet. In reality, however, the protoplanet grows at a time-dependent rate, which depends on the protoplanet’s mass and surrounding gas density, namely the gap structure. Here, adopting a simple model of protoplanetary growth via gas accretion, we demonstrate that the capture of planetesimals is affected by growth modes of the protoplanet.

5.1 Model of Protoplanet’s Growth via Gas Accretion

We consider protoplanet’s growth via gas accretion after the onset of the runaway gas accretion. In the first phase, the protoplanet’s growth is controlled by the Kelvin-Helmholtz contraction of the gaseous envelope (Bodenheimer & Pollack, 1986; Tajima & Nakagawa, 1997; Ikoma et al., 2000). We adopt the fitting formula of the growth rate numerically derived by Tajima & Nakagawa (1997),

[TABLE]

which is valid for $M_{{\rm p}}\gtrsim$ 30 $M_{\mathrm{\oplus}}$. Since the growth rate increases rapidly with increasing mass of the protoplanet, it exceeds the maximum supply rate of disc gas at a certain mass (Tanigawa & Watanabe, 2002). Thus, the second phase is characterised by supply-limited gas accretion and called the supply-limited phase.

In the second phase, we adopt the model of the supply-limited gas accretion formulated by Tanigawa & Ikoma (2007). Based on the numerical results of Tanigawa & Watanabe (2002), which demonstrated that disc gas flows into the protoplanet’s Roche lobe through the two bands around $x_{{\rm gap}}\equiv{|r-a_{{\rm p}}|}\sim 2R_{{\rm H}}$, Tanigawa & Ikoma (2007) derived the maximum supply of disc gas as

[TABLE]

where

[TABLE]

and $\Sigma_{\mathrm{band}}$ is the surface density of disc gas at $x_{{\rm gap}}=2R_{{\rm H}}$. The surface density $\Sigma_{\mathrm{band}}$ is affected not only by gap formation, but also global disc evolution, in particular, disc dispersal. Since the detailed treatment of global disc evolution is beyond the scope of this study, we assume that the surface gas density outside the gap $\Sigma_{\rm 0}$ is given simply as

[TABLE]

where $\Sigma_{{\rm 0,int}}$ is the initial, unperturbed disc surface density, $\tau_{{\rm disp}}$ is the timescale of disc dissipation, and $t_{{\rm c}}$ is a parameter that determines the onset of disc dissipation. The surface density in the accretion band $\Sigma_{\mathrm{band}}$ is given by combination of Eqs. (38) and (51), which is applied to Eq. (49) to obtain the maximum supply rate. In summary, the protoplanet’s growth rate is given as the smaller of $\dot{M}_{{\rm p,R}}$ and $\dot{M}_{{\rm p,S}}$, namely

[TABLE]

Figure 11 shows (a) the protoplanet’s growth rate and (b) the protoplanet’s mass thus calculated for three different choices of $\alpha_{{\rm\nu}}$ and $\tau_{{\rm disp}}$ as indicated in each panel. The value of $t_{{\rm c}}$ is chosen so that the protoplanet mass is 1 $M_{\mathrm{J}}$ at the end of integration. The black solid lines represent the case without gap formation nor disc dispersal (i.e., $\alpha_{{\rm\nu}}\rightarrow\infty$ and $\tau_{{\rm disp}}\rightarrow\infty$), which Shiraishi & Ida (2008) assumed. The two parameters control the protoplanet mass at which the transition from the contraction-limited and supply-limited growth occurs. The smaller the viscosity $\alpha_{{\rm\nu}}$, the smaller the transition mass is, because the gap is wider, reducing the gas surface density in the accretion band. In the early supply-limited growth phase, the growth rate decreases gradually, but drops rapidly in the late phase, because of disc dispersal. In every case, the growth rate in the late phase is smaller by several orders of magnitude than that in the case without gap formation nor disc dispersal. Such a low growth rate is expected to reduce the capture rate of planetesimals, as seen above. On the other hand, the phase of low growth rate lasts long, as shown in the Fig. 11b. The growth models are labelled in the order of growth speed as model C0-C3, namely, model C0 (black) is the fastest model and model C3 (blue) is the slowest model. The parameters used in those models are listed in Table  1. Below we investigate such competing effects on the total captured mass of planetesimals.

5.2 Calculation of capture of planetesimals

Here we perform orbital integration of planetesimals around the protoplanet growing at a rate modelled above. The integration is continued until the surface density becomes low enough that the protoplanet mass increases no more, namely, $\tau_{{\rm grow}}>1\times 10^{7}$ year ($\gg\tau_{{\rm disp}}$). The initial protoplanet mass $M_{{\rm p,int}}$ is assumed to be the same value adopted in the previous sections. The calculated cumulative mass of captured planetesimals is shown in Fig. 12a for four different sets of $\alpha_{{\rm\nu}}$ and $\tau_{{\rm disp}}$ (including $\alpha_{{\rm\nu}}\rightarrow\infty$ and $\tau_{{\rm disp}}\rightarrow\infty$; i.e., no gap and disc dispersal) and two different choices of $R_{{\rm pl}}$. In Fig. 12b, results are summarised in the total captured mass $M_{\rm cap}^{\rm total}$ with the size of the planetesimals $R_{{\rm pl}}$.

Firstly, we focus on $M_{\rm cap}^{\rm total}$. For the cases of $R_{{\rm pl}}=10^{7}\mathrm{cm}$, the fastest growth mode (C1; red) yields the largest total captured mass ($M_{\rm cap}^{\rm total}\sim 2.3M_{\mathrm{\oplus}}$). As shown in Section 4, the change in the damping timescale due to the gap formation hardly affects the dynamics of such large-sized planetesimals. $M_{\rm cap}^{\rm total}$ takes a similar value in models C0 (black) and C2 (green), but smaller in model C3 (blue). Thus, it follows that this difference in $M_{\rm cap}^{\rm total}$ comes from the change in the growth timescale, namely the difference of the growth model. For the cases of $R_{{\rm pl}}=10^{4}\mathrm{cm}$, the changes in both timescales result in smaller $M_{\rm cap}^{\rm total}$. In models C2 an C3, no small planetesimals are captured by the protoplanet.

Secondly, the slower the protoplanet growth is, the more rapidly the cumulative captured mass increases with the protoplanet mass in the early stage of planetary growth. This is because a given increase in protoplanet mass takes a longer time in the case of slower growth, so that more conjunctions and thus collisions between the planetesimals and protoplanet occur at a given protoplanet mass. For example, at $M_{{\rm p}}=100~{}M_{\mathrm{\oplus}}$, the slowest growth model (C3; blue) yields the largest value of $M_{{\rm cap}}$ among other models. Such a difference in timing of planetesimal bombardment may affect the distribution of heavy elements inside gas giant planets.

6 Other Effects

6.1 Capture of Planetesimals before Gap Formation

In this study, we have started the calculations at $M_{{\rm p}}\sim 30M_{\mathrm{\oplus}}$ to focus on the effect of gap formation, ignoring the amount of planetesimals that the protoplanet captures until $M_{{\rm p}}\sim 30M_{\mathrm{\oplus}}$. Here we make a simple estimate on how much planetesimals should be captured until $M_{{\rm p}}\sim 30M_{\mathrm{\oplus}}$.

The total mass of planetesimals initially existing inside the feeding zone of a protoplanet of $30~{}M_{\mathrm{\oplus}}$ at $5.2~{}{\rm AU}$ in MMSN is $\sim~{}3.3~{}M_{\mathrm{\oplus}}$ (Hayashi, 1981). According to numerical simulations of the runaway gas accretion via the Kelvin-Helmholtz contraction of the gaseous envelope (e.g., Ikoma et al., 2000), the protoplanet growth timescale is longer than $10^{5}~{}\mathrm{year}$ just after the onset of the runaway gas accretion. Such a growth timescale is longer than that in the model A with $\dot{M}_{{\rm p}}=10^{-3}~{}M_{\mathrm{\oplus}}~{}{\mathrm{year}}^{-1}$. According to our numerical results shown in Fig. 3, even in the case of large planetesimals such as $R_{{\rm pl}}=1\times 10^{7}~{}\mathrm{cm}$, the fraction of captured planetesimals is smaller than $0.20$ and, thus, the cumulative captured mass of planetesimals at $M_{{\rm p}}=30M_{\mathrm{\oplus}}$ is $\sim 0.66$ $M_{\mathrm{\oplus}}$ at most. Also in Shiraishi & Ida (2008), the numerical results show that the cumulative captured mass at $M_{{\rm p}}=30~{}M_{\mathrm{\oplus}}$ is smaller than $1~{}M_{\mathrm{\oplus}}$. Hence, a protoplanet is estimated to capture a greater amount of planetesimals after gap formation than it does until that. Thus, we conclude that the capture of planetesimals after gap formation is important to investigate the origin of the heavy elements observed inside gas giant planets.

6.2 Effective Capture Radius of the Growing Protoplanet

In all the simulations above, we assumed a constant mean density of the protoplanet $\rho_{{\rm ave}}=0.125~{}\mathrm{g}~{}/~{}{\mathrm{cm}}^{3}$, with which we calculated the protoplanet’s radius $R_{{\rm p}}$. During the formation stage, however, the protoplanet’s mean density is thought to be much lower and variable, because continuous gas accretion brings a large amount of entropy, making the protoplanet’s envelope greatly inflated. In addition, in the supply-limited growth stage, a circum-planetary disc is formed around the accreting protoplanet (e.g., Tanigawa & Watanabe, 2002). Both the expanding envelope and circumplanetary disc reduce the kinetic energy of incoming planetesimals and thereby enhance the capture probability of planetesimals. Here, we make a brief discussion on such an effect.

The interaction between planetesimals and protoplanet’s atmosphere before the onset of runaway gas accretion is well studied (e.g., Podolak et al., 1988; Inaba & Ikoma, 2003; Valletta & Helled, 2019). According to these studies, the gas drag from the atmosphere enhances the capture probability greatly. A capture radius of a protoplanet, inside which planetesimals lose parts of their kinetic energy enough to stay in the envelope, is ten times or more larger than the radius of its solid core. During the runaway gas accretion stage, the capture radius would be more enhanced by the inflated envelope. Although the capture radius enhanced by the circumplanetary disc is smaller than that of the inflated envelope in runaway gas accretion phase, it is several times larger than that without circumplanetary disc (e.g., Fujita et al., 2013; Tanigawa et al., 2014). In both cases, however, it is easy to say that the capture radius is much larger than the protoplanet radius we used in this study. Detailed investigation of the effects of such enhanced capture radius is one of our future studies.

7 Summary and Conclusion

Gas giant planets including Jupiter, Saturn, and some hot/warm Jupiters are known to contain much larger amounts of heavy elements relative to their host stars. Such heavy elements were likely brought by planetesimals in late stages of the growth of the proto-gas giants. In this study, performing numerical simulations of the dynamics of planetesimals around a protoplanet growing in mass via gas accretion, we have investigated the efficiency for planetesimals to be captured by the protoplanet.

Our main focus has been on the effects of the presence of a gap near the protoplanet’s orbit in the circumstellar gaseous disc, which previous studies ignored. Gap formation leads to reducing the disc surface density, enhancing the velocity shear between the planetesimals and disc gas by changing the gradients of surface density (or pressure), and slowing the protoplanet’s growth, the effects of which we have investigated step by step in detail.

The reduction in the disc surface density makes damping of planetesimals’ random velocities less efficient, which leads to enhancing the capture rate of planetesimals; 2. 2.

The enhancement of the velocity shear between the planetesimals and disc gas accelerates radial drift of planetesimals, which leads to reducing the capture rate; 3. 3.

The slowed protoplanet’s growth suppresses the rate of expansion of the protoplanet’s feeding zone, which leads to reducing the capture rate, whereas it increases the duration of protoplanet growth and thereby increases the number of conjunctions of planetesimals with the protoplanet, which leads to early capture of planetesimals.

Our main findings are as follows:

• •

The second effect (accelerated radial drift) always dominates over the first one (inefficient damping) in the parameter ranges considered in this study, thereby reducing the total mass of captured planetesimals relative to the case with no gap (see Section 4).

• •

The third effect (suppressed disc gas supply) makes no change or reduction in the total mass of captured planetesimals and, on the other hand, makes major planetesimal accretion occur earlier (see Section 5).

In conclusion, gap formation leads to a reduction in the total mass of captured planetesimals.

From our numerical simulations, we estimate that a Jupiter-mass planet captures, at most, 3 $M_{\mathrm{\oplus}}$ after the onset of runaway gas accretion until disc dispersal at 5.2 AU in the MMSN (see Section 6.1). For such large amounts of heavy elements inferred for several gas giants to accrete (e.g., Wahl et al., 2017), the solid surface density has to be more than five times higher than the MMSN value. Also, although detailed investigation being needed, the early accretion of planetesimals that we have found in Section 5 may be important for explaining the fact that heavy elements are depleted in the outer envelope and instead concentrated in the inner envelope and central core of Jupiter (Wahl et al., 2017).

For explaining dense warm Jupiters that are inferred to contain $>100~{}M_{\mathrm{\oplus}}$ heavy elements (Thorngren et al., 2016), we would need additional capture and/or supply processes such as planetary radial migration, during which the protoplanet would encounter planetesimals outside the feeding zone (Alibert et al., 2005). Such an investigation will be done in our forthcoming paper.

Appendix A Benchmark Test of Our Orbital Integration Code

For a benchmark test, we consider a Jupiter-mass planet in a circular Keplerian orbit and integrate the motion of test particles that orbit a central star of solar mass. We distribute the planetesimals using the same method of our simulation and check energy errors during numerical integration; the energy error is defined as

[TABLE]

where $E_{{\rm Jacobi,0}}$ is the initial Jacobi energy of the particle. The energy error comes from two components, the gravitational force of the central star and that of the protoplanet. The relative distance from the central star is almost constant and the energy error accumulates linearly with time. On the other hand, the relative distance from the protoplanet changes largely and the energy error depends on the distance of the closest encounter. We check the each component considering the cases of with and without close encounters with the protoplanet. Figure 13 shows the some of the results of numerical integration. The upper panels show the change in the energy error during numerical integration. The lower panels show the change in the relative distance between the particle and the planet. The left column shows the results for the case where the particle experiences no close encounter with the planet, whereas the right column shows the results for the case where the particle does so. In the case with no close encounter, the energy error accumulates linearly with time and it is expected that the energy error would be suppressed smaller than $\sim 10^{-7}$ after the $10^{5}\mathrm{years}$ integration. On a close encounter with the planet, the energy error increases sharply and its magnitude depends on the distance. It turns out that increase in energy error due to close encounters is suppressed to as small as $\sim 10^{-7}$ for $\eta_{{\rm ts}}$ = 0.01, if $r_{{\rm pl,p}}>10^{-3}{\rm AU}$. In our numerical study, all particles getting inside the radius of the planet are removed from the calculation and the protoplanet growth is stopped at Jupiter mass. Thus, setting the accuracy-controlling parameter as $\eta_{{\rm ts}}$ = 0.01 and considering the expanding envelope which has a mean density of $0.125\mathrm{g}/{\mathrm{cm}}^{3}$, the maximum energy error in the end of our simulation is as small as $\sim 10^{-7}$.

Appendix B Change of the Jacobi Energy

To derive the variations of the orbital elements $a$, $e$ and $i$, we consider an orbital system of cylindrical coordinates $(R,\phi,\zeta)$. In this coordinate, the velocity of planetesimals ${\bf v_{{\rm pl}}}$ is given as

[TABLE]

where, $\psi$ is a true anomaly. The velocity of disk gas ${\bf v_{{\rm gas}}}$ is given as

[TABLE]

where $l$ is the distance from the central axis of the circumstellar disc, and

[TABLE]

where $\varpi$ is a longitude of pericenter. Using these equations, the relative velocity between the planetesimals and ambient disc gas is given as ${\bf u}={\bf v_{{\rm gas}}}-{\bf v_{{\rm pl}}}$. The variations of the orbital elements are given by

[TABLE]

where ${\bf F}$ is the force acting on the planetesimals given by

[TABLE]

In our simulation, $F=f_{\rm gas}$. The averaged change of the Jacobi energy for one orbital period is given as

[TABLE]

Using above equations and Eq. (27), we know the dependence of the change of Jacobi energy on the growth and damping timescales.

Acknowledgements

We thank Tristan Guillot, Julia Venturini, Yann Alibert and Yuhiko Aoyama for valuable discussions and suggestions. We are also grateful to the reviewer Keiji Ohtsuki for his careful reading and constructive comments. This work is supported by JSPS KAKENHI Grant Numbers 17H01153 and 18H05439 and JSPS Core-to-Core Program “International Network of Planetary Sciences”. Numerical computations were carried out on the Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.