Long-term Hydrodynamic Simulations on the Planetesimals Trapped in the First-order Mean Motion Resonances

TL;DR

This study uses long-term hydrodynamic simulations to analyze how asymmetric disk structures influence the eccentricity and size thresholds of planetesimals trapped in first-order mean motion resonances, revealing significant deviations from previous symmetric assumptions.

Contribution

It provides revised analytical expressions for equilibrium eccentricity and minimum planetesimal size considering disk asymmetry effects, improving understanding of resonance trapping dynamics.

Findings

Eccentricity estimates are underestimated by 30-40% for j<5 resonances.

Asymmetric disks reduce equilibrium eccentricity and minimum size thresholds.

Eccentricity and size effects vary with planet mass and disk asymmetry levels.

Abstract

The resonant perturbations from planets are able to halt the drag-induced migration, and capture the inwardly drifting planetesimals into mean motion resonances. The equilibrium eccentricity of planetesimals in resonances, and the minimum size of planetesimal that can trigger resonance trapping, have been analyzed and formulated. However, the analytical works based on the assumption that the disk is axisymmetric, which is violated by the asymmetric structures developed by planets. We perform long-term 2D hydrodynamic simulations to study the dynamics of planetesimals in the $j:(j+1)$ first-order exterior resonances, and reexamine the theoretical expressions. We find the expression of equilibrium eccentricity underestimates the values for resonances with $j < 5$, in particular the 1:2 resonance that the underestimation can be $30 - 40\%$. Within the parameter space we explored, we find the equilibrium eccentricity and the minimum size are reduced in an asymmetric disk. The amount of discrepancy in eccentricity depends on the degree of asymmetric structures. For cases of Earth-sized planets, where the disk is less disturbed, the planetesimal's eccentricity can reach to the values predicted by our modified expression. For gaseous planets, however, the eccentricity can be $0.01 - 0.02$ smaller in value. We find the minimum size is 10 times smaller, and the factor seems to be independent of the planet's mass. The influences of asymmetric profiles on the eccentricity and the minimum size could affect the outcome of collisions between resonant and nonresonant planetesimals, and the amount of planetesimals migrated into the planet's feeding zone.

Tables

Table 1. Table 1: The minimum size of planetesimal s \text ​ m ​ i ​ n subscript 𝑠 \text 𝑚 𝑖 𝑛 s_{\text{min}} obtained from Eq. ( 21 ), and the employed planetesimal’s size s \text ​ p ​ l subscript 𝑠 \text 𝑝 𝑙 s_{\text{pl}} and the corresponding Stokes number for the simulation with an Earth-sized planet, μ = 10 − 5 𝜇 superscript 10 5 \mu=10^{-5}

Resonance	$s_{\text{text}min}$	$s_{\text{text}pl}$	Stokes Number
	(${10}^{-8}\text{text}au$/km)	(${10}^{-8}\text{text}au$/km)
1:2	60.0/89.8	100.0/149.6	$\sim {10}^6$
2:3	5.4/8.1	10.0/15.0	$\sim {10}^5$
3:4	2.6/3.9	5.0/7.5	$\sim {10}^5$
4:5	1.5/2.2	5.0/7.5	$\sim {10}^5$

Table 2. Table 2: The ratio of the numerical location and the nominal resonance location, Δ Δ \Delta , in different first-order resonances and planet-to-star mass ratios.

$\Delta$		$\mu ={10}^{-5}$	$\mu ={10}^{-4}$	$\mu ={10}^{-3}$
1:2	$1.0044$	$1.0042$	$1.0048$
2:3	$1.0031$	$1.0037$	-
3:4	$1.0027$	$1.0039$	-
4:5	$1.0025$	$1.0059$	-

Equations

$${\mathbf{F}}_D=-\frac{1}{2}{C}_D\pi s^2{\rho }_g\Delta V\mathbf{\Delta V},$$

$$\eta =-\frac{r}{2V_K^2{\rho }_g}\frac{dp}{dr}=\frac{1}{2}(1-\gamma )h^2$$

$$\displaystyle C_{D}=\left\{\begin{array}[]{ll}24Re^{-1}&\quad Re\leq 1\\ 24Re^{-0.6}&\quad 1<Re\leq 800\\ 0.44&\quad Re>800.\end{array}\right.$$

$${\mathbf{F}}_D=-\frac{4\pi }{3}{s}^2{\rho }_g{\bar{v}}_{th}\mathbf{\Delta V},$$

$$St=t_s{\Omega }_K,$$

$$t_s=-\frac{m_d\mathbf{\Delta V}}{{\mathbf{F}}_D},$$

$${\varphi }_i$$

$${\varphi }_o$$

$$\frac{a_i}{a_o}={\left(\frac{p}{p+q}\right)}^{2/3},$$

$$\frac{d{a}_o}{dt}\text{(res)}=-(j+1)\mu a_o{e}_o{n}_{o}C(\alpha )\sin \varphi$$

$$\frac{d{e}_o}{dt}\text{(res)}=-\frac{1}{2}\mu n_{o}C(\alpha )\sin \varphi ,$$

$$C(\alpha )=(2j+1+\alpha \frac{d}{d\alpha })b_{1/2}^{(j)}(\alpha )-\frac{1}{{\alpha }^2}{\delta }_{1j},$$

$$\frac{d{a}_o}{dt}\text{(drag)}=-2\frac{a_o}{{\tau }_0}[e_o\eta +(0.35-0.16\gamma )e_o^3]$$

$$\frac{d{e}_o}{dt}\text{(drag)}=-0.77\frac{e_o^2}{{\tau }_0},$$

$$\frac{da}{dt}\text{(res)}+\frac{da}{dt}\text{(drag)}=0$$

$$\frac{de}{dt}\text{(res)}+\frac{de}{dt}\text{(drag)}=0.$$

$$e_{\text{eq}}=1.14{\left(\frac{\Delta V/V_K}{j+1}\right)}^{1/2}.$$

$$s_{\text{min}}=\frac{{\rho }_g{a}_i(\Delta V/V_K{)}^{1/2}}{3{\rho }_d\mu C(\alpha )j^{3/2}},$$

$$\Delta t_{\text{RK5}}=\min (f_{ts}\times t_{s,\text{min}},0.2\times \Delta t_{\text{CFL}}).$$

$$v_{r,d}=\frac{S{t}^{-1}{v}_{r,g}-\eta V_K}{St+S{t}^{-1}},$$

$$\frac{|a-a_p|}{a_p}\le {\mu }^{2/7},$$

$$\lambda \sim \frac{4\times 10^9}{{\rho }_{g,\text{mid}}}\text{cm},$$

$$\frac{a_o}{a_i}={\left(\frac{j+1-{\dot{\varpi }}_o/n_o}{j}\right)}^{2/3}.$$

$$\Delta$$

$$\sim$$

$$\frac{\dot{\varpi }}{n}=-\left[\frac{F_D}{F_K}+\frac{1}{2}r\frac{d{F}_D/dr}{F_K}\right],$$

$$\frac{{\dot{\varpi }}_o}{n_o}=2\mu f_{s,1}(\alpha )+\frac{\mu }{2e_o}C(\alpha )\cos {\varphi }_o,$$

$$f_{s,1}(\alpha )=\frac{1}{8}\left[2\alpha \frac{d}{d\alpha }+{\alpha }^2\frac{d^2}{d{\alpha }^2}\right]b_{1/2}^{(0)}(\alpha ).$$

$$e_{\text{eq}}={\left[\frac{\Delta V/V_K}{0.77(j+1)-(0.35-0.16\gamma )}\right]}^{1/2},$$

Full text

hydrodynamics — methods: numerical — planet-disk interactions — planets and satellites: dynamical evolution and stability — protoplanetary disks

††journal: ApJ††software: IPython (Pérez & Granger, 2007), Matplotlib (Hunter, 2007), NumPy (van det Walt et al., 2011), SciPy (Jones et al., 2001), SymPy (Meurer et al., 2017)

1 Introduction

In protoplanetary disks, the dust particles follow the Keplerian orbits, but the gas rotates slightly less than the Keplerian velocity due to the support of the radial pressure gradient. The difference in velocities results in the “headwind” experienced by the dust particles that reduces their angular momenta and causes them to spiral toward the central star (Whipple, 1972; Adachi et al., 1976; Weidenschilling, 1977). When the drifting particles cross the orbit of a large object (e.g., protoplanets or planets), a fraction of the particles can be accreted onto the object due to the combined effects of gravitational attraction and gas drag. This process is known as pebble accretion, and provides a viable way for the rapid growth from protoplanets to planets (Ormel & Klahr, 2010; Lambrechts & Johansen, 2012).

On the other hand, the resonant perturbations from the planet are able to counteract the effects of gas drag, and trap the drifting particles into mean motion resonances (MMRs), before the orbital crossing. The resonance trapping mechanism was first investigated by Weidenschilling & Davis (1985), who derived the minimum size that the particles can be trapped in the first-order resonances. Patterson (1987) extended the work to high-order resonances. For particles smaller than the minimum size, the drag force is too large to be opposed by the resonant perturbation. The typical values of the minimum size ranges from a few kilometers to subkilometers.

For planetesimals trapped in resonances, their eccentricities are pumped by the resonant perturbations and damped by the gas drag, and eventually reach the equilibrium eccentricity with moderate values. The expression of the equilibrium eccentricity was derived by Weidenschilling & Davis (1985), and confirmed numerically via N-body simulations (e.g., Malhotra, 1993), and hydrodynamical simulations for the 1:2 resonance (Paardekooper, 2007; Ayliffe et al., 2012). The induced eccentricities can cause the overlap of orbits for planetesimals in different resonances, increase the collision rate between resonant planetesimals or between resonant and nonresonant planetesimals, and result in the destructive collisions (Weidenschilling & Davis, 1985). The fragments smaller than the minimum size can then pass through the resonances, and could be accreted by the planet.

The analytical studies on the resonance trapping mechanism assume an axisymmetric gaseous disk; however, the presence of a planet induces asymmetric structures in the disk, e.g., single/multiple spiral arms (e.g., Goldreich & Tremaine, 1979, 1980; Lee, 2016; Bae & Zhu, 2018; Miranda & Rafikov, 2019), a cavity/gap around the planet (Lin & Papaloizou, 1993; Bryden et al., 1999), and the eccentric modes (Kley & Dirksen, 2006; Teyssandier & Ogilvie, 2016; Lee et al., 2019). To study the influences of the asymmetric structures on the resonance trapping mechanism, hydrodynamics simulations or N-body simulations with a realistic gaseous profile are needed. On the other hand, the time the induced eccentricity takes to reach the equilibrium eccentricity can be hundreds to thousands of orbits, assuming no damping caused by gas drag (Patterson, 1987). This indicates that long-term simulations are required.

In this paper, we perform long-term hydrodynamic simulations to study the dynamics of planetesimals trapped in the first-order resonances, and reexamine the analytical expressions of the equilibrium eccentricity and the minimum size. The paper is structured as follows. In Section 2, we present the drag law adopted in this study, and introduce the mechanism of resonance trapping. The numerical method and the initial condition are described in Section 3. The simulation results are presented in Section 4. The conclusions and summary are in Section 5.

2 Theoretical model

2.1 Drag Law

The strength of drag force depends on the physical conditions of the dust particle and the surrounding gas. In cases of spherical particles, the drag force can be expressed in the form

[TABLE]

where $s$ is the particle radius, $\rho_{g}$ is the gas mass density, and $\Delta V$ is the relative velocity between the dust particle and the surrounding gas. The drag coefficient, $C_{D}$, is a function of three nondimensional quantities (see, e.g., Adachi et al., 1976):

The Mach number, $M=\Delta V/c_{s}$, where $c_{s}$ is the sound speed. 2. 2.

The Knudsen number, $Kn=\lambda/(2s)$, which is the ratio between the mean free path of gas molecules $\lambda$ and the particle diameter $2s$. 3. 3.

The Reynolds number, $Re=2s\Delta V/\nu$, where $\nu$ is the gas molecular kinematic viscosity.

In an isothermal disk, the relative velocity is $\Delta V\sim\eta V_{K}$, where $V_{K}$ is the Keplerian velocity,

[TABLE]

is the ratio of the gas pressure gradient to the stellar gravity in the radial direction, $\gamma$ is the slope of the mass density profile ($\rho_{g}\propto r^{\gamma}$), and $h$ is the aspect ratio of the gaseous disk. On the other hand, the sound speed is in the order of $h$, $c_{s}=hV_{K}$, and hence the Mach number $M\sim O(h)$. In a geometrically thin disk, $h\ll 1$, and we therefore neglect the effects of the Mach number on the drag coefficient.

We adopt the drag law in Weidenschilling (1977). For $Kn<2/9$, the drag coefficient is approximated by

[TABLE]

For $Kn>2/9$, we are in the Epstein regime. The drag force is given by

[TABLE]

where $\bar{v}_{th}=\sqrt{8/\pi}c_{s}$ is the mean thermal velocity. In order to make the drag force be continuous at the boundary ($Kn=2/9$), the gas molecular viscosity is set to $\nu=\bar{v}_{th}\lambda/2$.

One useful dimensionless quantity is the Stokes number (or dimensionless stopping time), for predicting the behavior of the dust particle influenced by the gas drag. The Stokes number is defined as

[TABLE]

where $\Omega_{K}$ is the Keplerian angular velocity, and $t_{s}$ is the stopping time given by

[TABLE]

and $m_{d}$ is the mass of the dust particle. The stopping time is the timescale the relative velocity decreases by a factor of $e$.

2.2 Resonance Trapping

For two bodies in the coplanar orbits, if at least one of the resonant arguments

[TABLE]

is in a libration and oscillates about either [math] or $\pi$, we say they are in the $p:(p+q)$ MMR (Murray & Dermott, 1999). The subscripts $i$ and $o$ denote the inner and outer bodies, respectively, and $\lambda$ is the mean longitude and $\varpi$ is the longitude of pericenter. The nominal resonance location, defined by

[TABLE]

provides an approximate location for the $p:(p+q)$ exterior resonance.

In this study, we focus on the $j:(j+1)$ first-order exterior resonances. The equations of motion for the outer body are (Weidenschilling & Davis, 1985; Murray & Dermott, 1999)

[TABLE]

[TABLE]

where $a$ is the semi-major axis, $e$ is the eccentricity, $\mu=m_{i}/m_{\star}$ is the planet-to-star mass ratio, $n$ is the mean motion, and $\alpha=a_{i}/a_{o}$ is the ratio of semi-major axes. The function $C$ is derived from the disturbing function:

[TABLE]

where $b^{(j)}_{1/2}$ is the Laplace coefficient, and $\delta_{1j}$ is the Kronecker delta function. For a body in the first-order exterior resonances, the stationary point is $\phi=\pi$, and hence its semi-major axis and eccentricity are pumped by the resonant perturbations.

On the other hand, the effects of gas drag on the dynamics of particles in Keplerian orbits were analyzed by Adachi et al. (1976). Given the specific force in the form of $f\propto\rho_{g}\Delta V^{2}$, which corresponds to the Stokes regime with $Re>800$, in the planar case their results give

[TABLE]

[TABLE]

where $\tau_{0}=t_{s}\eta$. Here, we assume $e_{o}>\eta$, and retain the formula corresponding the case of $e\gg\eta,i$ in Adachi et al. (1976).

The resonance trapping occurs if the effects of resonant perturbations and gas drag are equal and opposite:

[TABLE]

[TABLE]

The expression of equilibrium eccentricity can then be obtained by solving $de/da\text{(res)}=de/da\text{(drag)}$. To the leading term in $e$, this yields

[TABLE]

Then substituting Eq. (20) into Eq. (18), and letting $\sin\phi=-1$, the minimum size derived by Weidenschilling & Davis (1985) is

[TABLE]

where the outer body is assumed to be in the Stokes regime with $Re>800$, and $\rho_{d}$ is the bulk density. However, we note there is an error in this expression due to the incorrect relationship employed by Weidenschilling & Davis (1985), $n_{o}/n_{i}=(j/(j+1))^{3/2}$. With the correct relationship, $n_{o}/n_{i}=j/(j+1)$, the factor $j^{3/2}$ in the denominator should be $j^{2/3}(j+1)^{5/6}$ instead. This introduces an overestimation in the minimum size by a factor up to $1.8$ for $j=1$. In this study, we ignore this correction, and refer to the expression in Eq. (21) when mentioning the minimum size.

3 Numerical method

We perform the simulations using the hydrodynamical code FARGO3D version 1.3 (Benítez-Llambay & Masset, 2016), with the FARGO algorithm developed by Masset (2000). We implement the drag law outlined in Section 2.1 to include the drag force exerted by the gaseous disk onto the planets. The back reactions onto the gas phase are ignored in this study.

3.1 Algorithm for Gas Drag

The algorithms for the evolution of planets in FARGO3D are:

Update the velocity due to the disk potential with a time step constrained by the Courant–Friedrichs–Lewy condition, $\Delta t_{\text{CFL}}$. 2. 2.

Update the position and velocity, due to the potentials from the star and planets, via the fifth-order Runge-Kutta (RK5) integrator several times, with a smaller time step $\Delta t_{\text{RK5}}$. The value of $\Delta t_{\text{RK5}}$ is set to $0.2\times\Delta t_{\text{CFL}}$.

In our approach, the drag force is included in the source terms of velocity, and updated in real time, in the RK5 integration. We use the bilinear interpolation to approximate the density and velocity of the gas phase where the planet locates. In order to accurately simulate the dynamics of planets influenced by the drag force, the time step used in the RK5 integration must be smaller than the stopping time. Thus, we calculate the stopping time of all the embedded planets/planetesimals before step 2, and choose the minimum value (says $t_{s,\text{min}}$). The time step for each RK5 integration is then set to

[TABLE]

The factor $f_{ts}$ is set to $0.25$ to avoid the velocity difference decreasing by more than a factor of $e$ due to the drag force in one RK5 integration (Morbidelli & Nesvorny, 2012).

We test our code by comparing the drift velocities of dust particles obtained from the analytical expression and the numerical simulations. In each test run, we embed one particle in an axisymmetric disk. The particle is initially in a circular Keplerian orbit. The gaseous disk is assumed to be inviscid to avoid the disk evolution due to the shear force. We calculate the numerical drift velocity by averaging the simulation results over the last few orbits. On the other hand, the analytical expression of the drift velocity, $v_{r,d}$, is given by (Takeuchi & Lin, 2002)

[TABLE]

where the gas radial velocity is set to $v_{r,g}=0$. Figure 1 shows the inward drift velocity of dust particles with various sizes, obtained from Eq. (23) and simulations. Our simulation results agree with the theoretical expectation.

3.2 Disk Profile

We use the 2D fashion of the FARGO3D code. We assume a locally isothermal, razor-thin disk orbiting a $1~{}M_{\odot}$ star, where the aspect ratio is $h=0.05$. The disk profile is initially axisymmetric, and extends from $0.2$ to $4.0$ in the code unit. The unit of length is equal to $1$ au. The surface density is $1700$ g cm*-2* at $1$ au, corresponding to the Minimum Mass Solar Nebula model (Hayashi, 1981), and in proportion to $r^{-1}$ inferred from observations (Andrews et al., 2010). The disk’s viscosity is assumed to be uniform, and set to $\nu=10^{-5}$ in code units, which is equivalent to $\alpha=0.004$ at $1$ au in the alpha-disk prescription (Shakura & Sunyaev, 1973).

The disk is resolved by $384$ and $768$ grids in the radial and azimuthal directions, respectively. We use the damping boundary condition to reduce the unphysical reflecting waves near the boundary.

3.3 Planet/Planetesimal Profile

The planets and planetesimals are all embedded on the midplane of the gaseous disk. To explore the influences of the asymmetric structures in gaseous disks on the dynamics of planetesimals, we use three different planet-to-star mass ratios: $\mu=10^{-3}$, $10^{-4}$, and $10^{-5}$, which correspond to Jupiter, Neptune, and Earth, respectively. The planet is fixed in a circular Keplerian orbit, with the orbital radius $a_{p}=1$ au.

Duncan et al. (1989) showed that for a body within the region

[TABLE]

its behavior becomes chaotic. For $\mu=10^{-3}$, the region is $a\in[0.861,1.139]$ au, which is slightly smaller than the nominal resonance location of the 4:5 exterior resonance. Thus, we focus on the dynamics of planetesimals in $j:(j+1)$ resonances with $j<5$.

In each simulation, we place one planetesimal in circular orbit for each resonance. The sizes of planetesimals for simulations with $\mu=10^{-5}$ are listed in Table 3.3. The sizes are chosen to be slightly larger than the minimum size, based on the values in astronomical unit. For runs with $\mu=10^{-4}$ and $10^{-3}$, the sizes are $10$ and $100$ times smaller, respectively. We also show the corresponding Stokes number in Table 3.3. For a body with a large Stokes number, its dynamic is not strongly perturbed by the gas drag, and the drift velocity is small (see Figure 1). In order to reduce the simulation time required for the drag-induced migration, the initial orbital radii are slightly larger than the corresponding nominal resonance locations. We also adjust the locations such that the planetesimals are in near resonances with the embedded planet initially. In such a compact system, the planetesimals themselves are also in resonances. For simplicity, we ignore the interactions between planetesimals.

To calculate the drag force, we assume all the planetesimals are spherical, and the bulk density is set to $\rho_{d}=1.5$ g cm*-3*. The mean free path is adopted from Rice et al. (2006):

[TABLE]

where $\rho_{g,\text{mid}}$ is the gas mass density on the midplane of the disk.

Note that the force from disk potential is imparted onto the planet’s velocity directly in FARGO3D. For cases of subkilometer-sized planetesimals, the Hill sphere radius is about $10^{-7}$, which is too small to be used as the smoothing parameter. We find that the planetesimals are quite unstable in the simulations that employ the Hill radius as the smoothing length. Therefore, we use the thickness of the gaseous disk as the smoothing length, $\epsilon=0.4~{}H$, where $H$ is the pressure scale height.