Statistical Study of Planetesimal-Driven Migration

Abstract

The abstract serves both as a general introduction to the topic and as a brief, non-technical summary of the main results and their implications. Authors are advised to check the author instructions for the journal they are submitting to for word limits and if structural elements like subheadings, citations, or equations are permitted.

Abstract

Recent exoplanet observations have revealed a diversity of exoplanetary systems, which suggests the ubiquity of radial planetary migration. One powerful known mechanism of planetary migration is planetesimal-driven migration (PDM), which makes planets undergo significant migration through gravitational scattering with planetesimals. Here, we present the results of our high-resolution self-consistent $N$-body simulations of PDM, in which gravitational interactions among planetesimals, the gas drag, and Type-I migration are all taken into account. Our results show that even small protoplanets can actively migrate through PDM. Moreover, a fair fraction of them migrate outward. This outward migration can give a solution for the “planet migration problem” caused by Type-I migration and explain the origin of Jovian planets

Current planetary formation theories are built upon the classical standard theories [1, 2, 3]. In these theories, it is assumed that the terrestrial planets and the cores of giant planets grow “in-situ”, from solid components of the protoplanetary disk around their current location. Many problems have been pointed out for this in-situ growth model. It has been difficult to explain the formation of Uranus and Neptune within the solar system’s lifetime [4, 5]. Protoplanets grown to the size of Mars experience Type-I migration and will drift toward the central star within a million years [6, 7]. Moreover, since the discovery of an exoplanet by [8], more than 5000 exoplanets have been found and a fair fraction of them are hot Jupiters or super-Earths (e.g., [9]), and their existence is not compatible with the classic theory of in-situ formation. These discoveries challenge the classical standard theories of planetary formation and necessitate consideration of planetary migration mechanisms.

In recent years, planetesimal-driven migration (PDM) has gained attention as a possible mechanism for migration [10, 11, 12, 13, 14, 15]. PDM is initiated by the emergence of asymmetry in the distribution of planetesimals around the planets [16]. The asymmetric distribution leads to an imbalance of the torques acting on the planets, resulting in the self-sustained planetary migration in one direction. The existence of plutinos suggests that Neptune has moved outward slowly during its formation process, and PDM is one possible mechanism for this outward migration [17, 18, 19, 20]. However, whether or not such an outward migration through PDM actually occurs has been an open question, as realistic simulations in which gravitational interactions among planetesimals, gas drag, and Type-I torque in planetary migration are taken into account have not been conducted [10, 11, 12, 13, 14, 15]. In this paper, we report the results of 570 such realistic self-consistent $N$-body simulations of PDM (all models are listed in Table 1 in SI1).

We first present the results of our standard model, followed by a discussion of the mass ratio effects between the protoplanets and planetesimals, and conclude with the implications of our findings.

The planetary migration through PDM

In this section, we present the result of our $N$-body simulations for the orbital evolution of one protoplanet through PDM. In all simulations, a single protoplanet is placed within a planetesimal disk. Both the planetesimal and gas disks are modeled as axisymmetric and radially smoothed structures, consistent with standard theories (refer to Methods). Due to computational resource constraints, the planetesimal disk is initially limited to the radial range of 2 au to 6 au.

We conducted 10 runs from the same initial protoplanet mass ($0.5~{}M_{\oplus}$) and varied only the initial particle distribution and another 10 runs changing the initial mass of the protoplanet. The first 10 runs are referred to as model 1 and the second 10 runs as model 2 hereafter. Our simulation results show that two runs exhibited outward PDM in both models 1 and 2. Furthermore, we found that the trend of outward PDM is independent of the planetary mass under the same initial disk conditions.

Fig. 1 shows the time evolution of semi-major axes of protoplanets in model 1. In this model, protoplanets exhibited outward PDM in two out of 10 runs. Note that we did not give any initial kick to the protoplanets. Each protoplanet is placed on a circular orbit. One can observe that once the outward migration started, it continued monotonically up to the outer cutoff of the planetary disk at 6 au. However, when the protoplanet reached the outer cutoff of the planetesimal disk, it reversed direction and began migrating inward. This change in the direction is clearly the effect of the outer cutoff of the disk. In the remaining eight runs, monotonic inward PDM took place and protoplanets reached to the inner cutoff of the planetesimal disk. Here, the change to outward migration was not observed, and protoplanets continued to drift inward through Type-I migration.

In SI2, we discuss how the protoplanet interacts with the surrounding planetesimals during its migration. Once the protoplanet starts to migrate, it scatters the cold planetesimals in front of it and exchanges the angular momentum. When the protoplanet scatters the planetesimals located outside, it generally gains angular momentum and moves outward. Thus, once the protoplanet starts to move outward (inward), it continues to scatter planetesimals located outside (inside), and gains (loses) angular momentum.

Figure 1: The time evolution of semi-major axes of protoplanets in model 1. Each solid line represents the result of 10 runs with different initial seeds. The red lines represent the results of runs shown in Fig. 7 in SI2. The greyed hatch region shows the initial planetesimal disk size.

Fig. 2 shows the time evolution of semi-major axes of protoplanets in model 2. In this model, runs with initial planet mass $0.2~{}M_{\oplus}$ and $0.6~{}M_{\oplus}$ (represented by thin dashed and thick solid lines, respectively) showed outward PDM. As was the case in model 1, all planets exhibited monotonic PDM both inward and outward except for that in the run with $0.2~{}M_{\oplus}$. In this particular run, the direction of planet migration through PDM switched from inward to outward around $t=5\times 10^{3}$ years. Another notable feature of this figure is that the migrations of the least massive protoplanets ($0.1~{}M_{\oplus}$ and $0.2~{}M_{\oplus}$) seem slower than those of others.

Figure 2: The time evolution of semi-major axes of protoplanets in model 2. The thin blue solid, orange dashed, green dotted, red dashed-dotted, and purple double-dashed lines represent simulation results for initial planetary masses of 0.1, 0.2, 0.3, 0.4, and 0.5 $M_{\oplus}$, respectively. The thick brown solid, pink dashed, grey dotted, dark-yellow dashed-dotted, and cyan double-dashed lines represent simulation results for initial planetary masses of 0.6, 0.7, 0.8, 0.9, and 1 $M_{\oplus}$, respectively. The greyed hatch region shows the initial planetesimal disk size.

Figs. 1 and 2 indicate that planets can actively migrate within the protoplanetary disk through PDM both inward and outward. This result suggests that planets can migrate outward even in the presence of Type-I torque.

The effect of mass resolution

It has been assumed that whether PDM takes place or not depends on the mass ratio between planetesimals and the protoplanet [12, 14] (see also SI4). Here we investigate this dependency of PDM to the mass ratio using results of $N$-body simulations with the number of particles ranging from $10^{3}$ to $1.2\times 10^{5}$.

Fig. 3 shows the time evolution of the semi-major axes of protoplanets and their final distributions for our 350 $N$-body simulations. We have used six different values for the total number of planetesimals, $N=~{}10^{3},~{}2.5\times 10^{3},~{}5\times 10^{3},~{}10^{4},~{}3\times 10^{4}$, and $1.2\times 10^{5}$. They are referred to as models 3-1 to 3-5 and 4-1. For all models except 4-1, we tried 50 runs from different initial random number seeds, and for model 4-1 100 runs.

From Fig. 3 we can see that, as the number of planetesimals increases, planetary migration becomes more monotonic. In other words, for a smaller number of planetesimals, planetary migration is more stochastic. This trend highlights the distinct nature of planetary migration through PDM and can be quantitatively explained by examining the mean, variance, and flip frequency of the planetary semi-major axis across models 3-1 to 3-5 and 4-1.

Note that our result clearly indicates that PDM can drive the migration of protoplanets even when the mass ratio between the protoplanets and planetesimals is less than 10. Our result is different from the claims of a previous work [14] that the mass ratio should be more than 100 for PDM to be effective. We believe this difference is simply because we have performed a very large number of runs to quantitatively determine the distribution of results. Such a survey has not been conducted previously.

Figure 3: The time evolution of the semi-major axes of protoplanets and the final distribution of them in models 3-1 to 3-5 and 4-1. Each panel includes 50 runs with different initial seeds with the same number of particles except model 4-1 which includes 100 runs. Histograms are divided into bins of 0.01 au from 2.4 to 3.6 au. The two Gaussian functions within each panel are approximations of the histogram using a mixture Gaussian distribution with 2 clusters: $f_{\mathrm{pdf}}(x)=\sum_{i=1}^{2}\pi_{i}/(\sqrt{2\pi}\sigma_{i})\exp{\left[-(x-\mu_{i})^{2}/\sigma_{i}^{2}\right]}$. Here $x$, $\mu_{i}$, $\sigma_{i}$, and $\pi_{i}$ are the semi-major axis, the expected value, the standard deviation, and the weight of the Gaussian function. We note that when the index $i$ is 1, it represents the parameters of the Gaussian distribution for inward PDM, and when it is 2, it represents the parameters for outward PDM.

Figure 4: The standard deviation $\sigma_{1}$ and $\sigma_{2}$ of the final location of protoplanets plotted as functions of the number of planetesimals $N$.

Fig. 4 shows the standard deviation of the final location of protoplanets as a function of $N$. We can see that the deviation is roughly proportional to $1/\sqrt{N}$. This trend can be understood by looking at the frequency of the flipping of the direction of planetary migration $N_{\mathrm{flip}}$, since a higher flip frequency should imply that the PDM is more stochastic. Fig.5 shows $N_{\mathrm{flip}}$ as function of $N$. It indicates that $N_{\mathrm{flip}}$ is proportional to $N^{-1}$. Thus, in models with fewer particles, the increase in the flip frequency causes inhibition of the monotonous planetary migration through PDM. As a result, the width of the distribution of the final locations becomes larger. This $N^{-1}$ dependence can be understood as the effect of two-body relaxation between planetesimals and protoplanets [21] (see also SI5).

Figure 5: The average flip frequency of planetary migration direction per run for models 3-1 to 3-5. Given that a “flip” means a change in the planetary migration direction, we define a flip as the point in time when the time derivative of the planet’s semi-major axis becomes zero (see SI5). We use the Savitzky-Golay (SG) filter to smooth the time evolution of the planetary semi-major axis to investigate global changes in planetary migration direction [22].

Figure 6: The normalized average migration rate of models 3-1 to 3-5 and 4-1. Circles and squares represent inward and outward PDM, respectively.

Fig. 6 shows the average migration rates for models 3-1 to 3-5 and 4-1. The migration rates are normalized by the fiducial migration rate given by [10, 12],

$$\left|\frac{\mathrm{d}a_{\mathrm{p}}}{\mathrm{d}t}\right|_{\mathrm{fid}}\approx\frac{4\pi\Sigma_{\mathrm{dust}}a_{\mathrm{p}}^{2}}{M_{\mathrm{Sun}}}\frac{a_{\mathrm{p}}}{T},$$ (1)

where $a_{\mathrm{p}},~{}\Sigma_{\mathrm{dust}}$ and $T$ are the planet’s semi-major axis, the surface density of the planetesimal disk, and the orbital period of the planet. Fig. 6 shows that the normalized migration rate is higher for large $N$ and reaches 0.5 for inward migration and 0.35 for outward migration. Our result is consistent with that of [12], which suggested that the normalized PDM rates were typically in the range of $0.2-0.7$ and that the inward migration rate is higher than the outward migration rate. However, we have clearly shown that PDM rate only gradually decreases when we reduce the number of planetesimals (and thus the mass ratio). The monotonous nature of PDM becomes less clear.

We conducted additional simulations both without Type-I torque and in Keplerian gas environments to examine the effects of gas drag and Type-I torque on planetary migration (see SI3 for more details). We confirmed that PDM is the primary driver of the migration, with other factors having a minor effect.

The “dynamic” planetary formation picture

We conducted 570 self-consistent $N$-body simulations of PDM, in which gravitational interactions among planetesimals, the gas drag, and Type-I migration are taken into account. In all simulations, we follow the evolution of a single protoplanet embedded in the radially smoothed planetesimal disk. We studied the effect of PDM on planetary formation using the results of $N$-body simulations. Our main findings are summarized as follows:

1. We found that planets can overcome the drift toward the central star caused by Type-I migration, with some planets migrating outward through PDM (Figs. 1 and 2).

2. Planets migration through PDM occurs for the mass ratio between protoplanet and planetesimals less than 10. This value is much smaller than those previously proposed [14, 15] (Figs. 3 and 6).

3. The degree of monotonicity in PDM is proportional to the number of particles $N$, and this can be explained by the two-body relaxation between planetesimals and planets (Fig. 5).

Our results indicate that even during the runaway growth [23, 24, 25], planetary embryos can actively migrate within the protoplanetary disk through PDM. Furthermore, the occurrence of the outward PDM in our simulations suggests that the core of a Jovian planet might have first formed in the inner region of the disk and migrated outward, and the radial distance of outward migration can be very large.

We need to reconsider the traditional view of “in-situ” planet formation scenarios and start to construct a new model, where PDM is taken into account as the primary mechanism of migration.

Methods

The disk model

We consider a gas disk around a Solar-type star. The gas disk is assumed to be axisymmetric and the gas distribution is similar to the classical model known as the minimum-mass solar nebula model (MMSN) [2]. Thus, the gas surface density, the gas density, and the gas disk temperature are given by:

$$\Sigma_{\mathrm{gas}} =2400f_{\mathrm{gas}}\left(\frac{r}{1\mathrm{au}}\right)^{p}\exp{\left(-\frac{t}{1\mathrm{Myr}}\right)}~{}\mathrm{g~{}cm}^{-2},$$ (2)

$$\rho_{\mathrm{gas}} =1.4\times 10^{-9}f_{\mathrm{gas}}\left(\frac{r}{1\mathrm{au}}\right)^{\alpha}\exp{\left(-\frac{t}{1\mathrm{Myr}}\right)}~{}\mathrm{g~{}cm}^{-3},$$ (3)

$$T =2.8\times 10^{2}\left(\frac{r}{1\mathrm{au}}\right)^{\beta}~{}\mathrm{K},$$ (4)

where $f_{\mathrm{gas}}$ and $r$ are the gas scaling factor and the radial distance from the central star. Following [2], we set values of $f_{\mathrm{gas}}$ to 0.71 and radial dependencies in equations (2) to (4), namely $p$, $\alpha$ and $\beta$ to -3/2, -11/4 and -1/2, respectively. We model the gas dissipation as an exponential decay with a timescale of 1 Myr.

Similar to the gas disk model, the dust distribution is given as:

$$\Sigma_{\mathrm{dust}}=\begin{cases}40f_{\mathrm{dust}}\left(\frac{r}{1\mathrm{au}}\right)^{p}~{}\mathrm{g~{}cm}^{-2}&(r<r_{\mathrm{ice}})\\ 42f_{\mathrm{ice}}\left(\frac{r}{1\mathrm{au}}\right)^{p}~{}\mathrm{g~{}cm}^{-2}&(r\geq r_{\mathrm{ice}}),\end{cases}$$ (5)

where $f_{\mathrm{dust}}$ and $r_{\mathrm{ice}}$ are the dust scaling factor and the snowline. The values of $f_{\mathrm{dust}}$ and $p$ are set to 0.71 and -3/2. The snowline in our model is assumed to be at 2.0 au, as it may have been closer to the Sun due to the viscous accretion of the gas disk and the Sun’s stellar evolution [26]. In addition, $f_{\mathrm{ice}}$ is the dust scaling factor beyond the snowline, set to $f_{\mathrm{ice}}=4f_{\mathrm{dust}}=2.84$. Here we set the value of $\Sigma_{\mathrm{dust}}$ to be four times that in the MMSN [14]. The scope of our simulations is restricted to the region beyond the snowline; thus, we use only the latter case of equation (5) to represent the dust distribution.

The gas drag model

We employ a force formula of gas drag on planetesimals developed by [27]:

$$\boldsymbol{F}_{\mathrm{drag}}=-\frac{1}{2}C_{\mathrm{D}}\pi r_{\mathrm{p}}^{2}\rho_{\mathrm{gas}}|\Delta\boldsymbol{v}|\Delta\boldsymbol{v},$$ (6)

where $C_{\mathrm{D}}$, $r_{\mathrm{p}}$ and $\Delta\boldsymbol{v}$ are the gas drag coefficient, the planetesimal radius, and the relative velocity of the planetesimal to the gas disk. We assume that the gas disk is laminar, where the magnitude of the disk gas velocity is given by $v_{\mathrm{K}}(1-|\eta|)$ with the local Keplerian velocity $v_{\mathrm{K}}$. Here $\eta$ is a dimensionless quantity that characterizes the pressure gradient of the gas disk which is given by:

$$\eta=-\frac{1}{2}\frac{c_{\mathrm{s}}^{2}}{r^{2}\Omega^{2}}\left(\frac{\partial\log(\rho_{\mathrm{gas}}T)}{\partial\log r}\right),$$ (7)

where $c_{\mathrm{s}}$ and $\Omega$ are the speed of sound and the Keplerian orbital frequency. The functional forms of $c_{\mathrm{s}}$ and $\Omega$ are:

$$c_{\mathrm{s}} =\left(\frac{k_{\mathrm{B}}T}{\mu m_{\mathrm{H}}}\right)^{1/2}=1.0\times 10^{5}\left(\frac{T}{280\mathrm{~{}}{K}}\right)^{1/2}\mathrm{cm}~{}\mathrm{s}^{-1},$$ (8)

$$\Omega =\left(\frac{GM_{*}}{r^{3}}\right)^{1/2}=2.0\times 10^{-7}\left(\frac{r}{\mathrm{AU}}\right)^{-3/2}~{}~{}\mathrm{s}^{-1},$$ (9)

where $k_{\mathrm{B}}$, $\mu=2.34$, $m_{\mathrm{H}}$, $G$ and $M_{*}$ are the Boltzmann constant, the mean molecular weight of the gas, the mass of a hydrogen atom, the gravitational constant and the mass of the central star.

Type-I migration model

For Type-I migration, we follow the equation of motion proposed by [28]. The equation is given as:

$$\frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}t}=-\frac{v_{\mathrm{K}}}{2\tau_{a}}\boldsymbol{\mathrm{e}}_{\theta}-\frac{v_{r}}{\tau_{e}}\boldsymbol{\mathrm{e}}_{r}-\frac{v_{\theta}-v_{\mathrm{K}}}{\tau_{e}}\boldsymbol{\mathrm{e}}_{\theta}-\frac{v_{z}}{\tau_{i}}\boldsymbol{\mathrm{e}}_{z},$$ (10)

where $\boldsymbol{v}$, $\tau_{a}$, $\tau_{e}$ and $\tau_{i}$ are the planetary velocity, the evolution timescales for the semi-major axis, eccentricity, and inclination. These evolution timescales are given by Appendix D of [28]:

$$\tau_{a}^{-1} \simeq C_{\mathrm{T}}h^{2}\left[1+\frac{C_{\mathrm{T}}}{C_{\mathrm{M}}}(\hat{e}^{2}+\hat{i}^{2})^{1/2}\right]^{-1}t_{\mathrm{wave}}^{-1},$$ (11)

$$\tau_{e}^{-1} \simeq 0.780\left[1+\frac{1}{15}(\hat{e}^{2}+\hat{i}^{2})^{3/2}\right]^{-1}t_{\mathrm{wave}}^{-1},$$ (12)

$$\tau_{i}^{-1} \simeq 0.544\left[1+\frac{1}{21.5}(\hat{e}^{2}+\hat{i}^{2})^{3/2}\right]^{-1}t_{\mathrm{wave}}^{-1},$$ (13)

where $h=H/r$ is the gas disk’s aspect ratio, and the scale height of the disk $H$ is given as:

$$H=\frac{c_{\mathrm{s}}}{\Omega}=5.0\times 10^{11}\left(\frac{T_{\mathrm{m}}}{280\mathrm{~{}}{K}}\right)^{1/2}\left(\frac{r}{\mathrm{AU}}\right)^{3/2}~{}~{}\mathrm{cm}.$$ (14)

Additionally, $\hat{e}$ and $\hat{i}$ are the eccentricity and the inclination scaled with $h$. For equation (11), $C_{\mathrm{T}}=2.73-1.08p-0.87\beta$, $C_{\mathrm{M}}=-6(2p-\beta-2)$, where $p$ and $\beta$ are given in equations (2) and (4). The characteristic timescale $t_{\mathrm{wave}}$ is given by [7] as:

$$t_{\mathrm{wave}}^{-1}=\left(\frac{M_{\mathrm{p}}}{M_{*}}\right)\left(\frac{\Sigma r^{2}}{M_{*}}\right)h^{-4}\Omega_{\mathrm{K}},$$ (15)

where $M_{\mathrm{P}}$ is the planetary mass.

Numerical setup

We conducted a total of 570 self-consistent large-scale $N$-body simulations using GPLUM [29] on the Japanese supercomputer Fugaku. The following section details the numerical setup of all the simulations.

In every model, we employ an axisymmetric surface density distribution for the planetesimal disk as following equation (5) with inner and outer cutoffs, $r_{\mathrm{in}}$ and $r_{\mathrm{out}}$. We set $r_{\mathrm{in}}=2$ au and $r_{\mathrm{out}}=6$ au in all models. The physical density of the planetesimals is 2 $\mathrm{g~{}cm}^{-3}$. Initially, eccentricities and inclinations of planetesimals follow a Gaussian distribution with the dispersion $\langle e^{2}\rangle^{1/2}=2\langle i^{2}\rangle^{1/2}=2r_{\mathrm{Hill}}/a_{\mathrm{p}}$, where $r_{\mathrm{Hill}}$ is the Hill radius of the planetesimal with the semi-major axis $a_{\mathrm{p}}$ [30]. The Hill radius is given as:

$$r_{\mathrm{Hill}}=\left(\frac{M_{\mathrm{p}}}{3M_{*}}\right)^{1/3}a_{\mathrm{p}}.$$ (16)

To study the orbital evolution of planets through PDM, planetesimals totaling a mass ranging from 0.1 to 1.0$M_{\oplus}$ were extracted from the region starting at 3 au outward in the initial setup. These planetesimals were replaced by a single planet on a non-inclined circular orbit, positioned at the center of the gap created by the removed planetesimal ring. We note that in all our $N$-body simulations, perfect accretion is assumed in both planetesimal-planetesimal and planet-planetesimal collisions, with the effects of fragmentation disregarded.

The total of 570 simulations is categorized into three groups based on various factors: the initial number of particles ($N_{\mathrm{p,init}}$), the number of planetesimals converted into a planet ($M_{\mathrm{p}}/m_{\mathrm{p}}$), the initial planet mass ($M_{\mathrm{p}}$), computational time ($t$), the number of runs ($N_{\mathrm{run}}$) and the gas disk-planet interaction (for more details, see SI1 and Table 1).

Supplementary information

Supplementary Information I: Summary of simulation models

In this study, we conducted a total of 570 $N$-body simulations of PDM. These simulations are categorized into three groups. Below, we provide an overview of each group:

1. A group dedicated to investigating the long-term orbital evolution of planets due to PDM (models 1 and 2).

2. A group focused on examining the mass resolution dependency of PDM (models 3-1 to 3-5).

3. A group dedicated to studying the dependency of PDM on the interaction among the gas disk, planetesimals, and planets (models 4-1 to 4-3).

Table 1: List of models

Model	$N_{\mathrm{p,init}}$	$m_{\mathrm{p}}$ [$M_{\oplus}$]	$M_{\mathrm{p}}/m_{\mathrm{p}}$	$M_{\mathrm{p}}$ [$M_{\oplus}$]	Time [yr]	$N_{\mathrm{run}}$	Type-I	$F_{\mathrm{drag}}$
1	$1.2\times 10^{5}$	$4.9\times 10^{-4}$	1035	0.5	$10^{5}$	10	$\checkmark$	$\checkmark$
2	$1.2\times 10^{5}$	$4.9\times 10^{-4}$	$207\times$ (1–10)	$0.1\times$(1–10)	$10^{5}$	10	$\checkmark$	$\checkmark$
3-1	$1\times 10^{3}$	$5.88\times 10^{-2}$	9	0.5	$5\times 10^{3}$	50	$\checkmark$	$\checkmark$
3-2	$2.5\times 10^{3}$	$2.35\times 10^{-2}$	22	0.5	$5\times 10^{3}$	50	$\checkmark$	$\checkmark$
3-3	$5\times 10^{3}$	$1.17\times 10^{-2}$	43	0.5	$5\times 10^{3}$	50	$\checkmark$	$\checkmark$
3-4	$1\times 10^{4}$	$5.88\times 10^{-3}$	86	0.5	$5\times 10^{3}$	50	$\checkmark$	$\checkmark$
3-5	$3\times 10^{4}$	$1.96\times 10^{-3}$	259	0.5	$5\times 10^{3}$	50	$\checkmark$	$\checkmark$
4-1	$1.2\times 10^{5}$	$4.9\times 10^{-4}$	1035	0.5	$5\times 10^{3}$	100	$\checkmark$	$\checkmark$
4-2	$1.2\times 10^{5}$	$4.9\times 10^{-4}$	1035	0.5	$5\times 10^{3}$	100	N/A	$\checkmark$
4-3	$1.2\times 10^{5}$	$4.9\times 10^{-4}$	1035	0.5	$5\times 10^{3}$	100	N/A	$\checkmark$

The first group, which includes models 1 and 2 consists of 20 simulations. All these simulations were conducted with an initial number of particles $N=1.2\times 10^{5}$ and the simulation time up to $t=10^{5}$ years. In model 1, 1035 planetesimals, initially located outward from 3 au within the disk, were replaced by a single planet ($M_{\mathrm{P}}=0.5~{}M_{\oplus}$) with a non-inclined circular orbit. Model 2 involved replacing a swarm of planetesimals located outward from 3 au in the planetesimal disk with a single planet of varying mass on a non-inclined circular orbit. The mass of these planets ranged from 0.1 to 1 $M_{\oplus}$ in increments of 0.1 $M_{\oplus}$.

The second group consists of 250 simulations. These simulations involved varying the initial number of particles, $N=1\times 10^{3},~{}2.5\times 10^{3},~{}5\times 10^{3},~{}1\times 10^{4},3\times 10^{4}$ and were conducted for each number of particles with 50 simulations each, over $t=5\times 10^{3}$ years. In all simulations, we made 50 different initial seeds to generate $N$-body realization of the planetesimal disk following the dust distribution given by equation (5) for each number of particles. In each simulation, the number of planetesimals initially located around 3 au from the central star within the disk, equivalent to $0.5~{}M_{\oplus}$, was replaced by a single planet ($M_{\mathrm{P}}=0.5~{}M_{\oplus}$) with a non-inclined circular orbit. We refer to it as model 3-1 to 3-5 depending on the number of particles.

The third group consists of 300 simulations. In these simulations, we conducted three types of simulations: (1) Type-I torque and the gas drag, (2) the gas drag, and (3) the Keplerian gas drag (given as equation (6), with the magnitude of the disk gas velocity by $v_{\mathrm{K}}$). Each type of simulation was run 100 times with different seeds until $t=5\times 10^{3}$ years. Same as the first and second groups, in all of these simulations, we followed the dust distribution given by equation (5). In all simulations, 1035 planetesimals initially located around 3 au from the central star within the disk were replaced by a single planet ($M_{\mathrm{P}}=0.5~{}M_{\oplus}$) with a non-inclined circular orbit. We refer to it as models 4-1 to 4-3 depending on the type of the gas disk-planet interactions used in the simulation.

Table.1 summarizes our simulation models. The surfaces of Table.1 show the names of the models, the initial number of particles $N_{\mathrm{p,init}}$, the initial mass of a planetesimal $m_{\mathrm{p}}$, the mass ratio $M_{\mathrm{p}}/m_{\mathrm{p}}$, the initial mass of a planet $M_{\mathrm{p}}$, the computational time $t$, the number of runs for each model, the presence of Type-I torque and the gas drag $F_{\mathrm{drag}}$. We note that the magnitude of the disk gas velocity for model 4-3 is set to $v_{\mathrm{K}}$.

Supplementary Information II: The detailed snapshots of model 1

The planetary migration through PDM is driven by the emergence of asymmetry in the density distribution within the planetesimal disk around the planets. Here we examine the interaction between a protoplanet and the planetesimal disk to understand the physical process of PDM.

Fig. 7 shows the snapshots of two simulations from model 1 in the eccentricity and semi-major axis plane. Three panels (a to c) on the left side of Fig. 7, the planet migrates inward while gravitationally scattering nearby planetesimals. Conversely, in the three panels (d to f) on the right side, the planet moves outward. In both inward and outward cases in Fig. 7, planetesimals located in the direction of planetary migration are more intensely scattered. This intense interaction between the planet and planetesimals in its migration path leads to an asymmetry in angular momentum transfer among planetesimals on both sides of the planet, thereby inducing monotonic inward/outward PDM. This result indicates that planets can potentially overcome inward Type-I migration and some migrate outward when planet-planetesimal interaction is taken into account.

Figure 7: Snapshots of two simulations from model 1. The vertical and horizontal axes represent eccentricity and semi-major axis. The left (a to c) and right (d to f) sides of the figure show the snapshots of inward PDM and outward PDM runs, arranged chronologically from the top to bottom panel. Red and black dots represent planets and planetesimals, respectively.

Supplementary Information III: The effect of disk-planet interaction

In the presence of a gas disk, planets undergo Type-I migration, and gas drag influences the dispersion of planetesimal velocities. This interaction plays a significant role in the process of PDM, highlighting the crucial influence of the gas disk on the overall dynamics of planetary migration[13, 15]. Here, we present the results of our simulations with and without Type-I torque and in Keplerian gas environments to examine the influence of the gas disk on PDM.

Fig. 8 shows the time evolution of semi-major axes of planets and histograms of the planet’ semi-major axes at the end of simulations in models 4-1 to 4-3. Horizontally across Fig. 8, the panels show the results accounting for (a) Type-I torque and the gas drag (model 4-1), (b) only the gas drag (model 4-2), and (c) the Keplerian gas drag (model 4-3), respectively. Each model includes the results of 100 runs with different seeds. Each histogram is divided into bins of 0.01 au from 2.5 to 3.5 au. Same as Fig. 3, the two Gaussian functions within each panel of Fig. 8 represent approximations of the histograms using a mixture Gaussian distribution with 2 clusters. The legends in each panel represent the parameter for mixture Gaussian distributions of models 4-1 to 4-3. These parameters are the expected values $\mu_{i}$, the variances $V_{i}$, and weights $\pi_{i}$ of the Gaussian distributions. When the index $i$ is 1, it represents the parameters of the Gaussian distribution for inward PDM, and when it is 2, it represents the parameters for outward PDM. It should be noted that the parameters $\pi_{1}$ and $\pi_{2}$ correspond to the fraction of samples in each respective cluster.

According to Fig. 8, a fair fraction of planets migrate outward through PDM, regardless of the presence of Type-I torque or gas drag. However, the values of $\pi_{i}$ in each model indicate that the proportion of inward migration is higher in model 4-1, which includes Type-I torque, compared to models 4-2 and 4-3, where Type-I torque is absent. This trend is attributed to the continuous inward kick provided by Type-I torque. Overall, our results indicate that although Type-I torque might favor inward PDM, outward PDM can occur universally, independent of the presence of a gas disk.

Figure 8: The time evolution of the semi-major axes of protoplanets and the final distribution of them in models 4-1 to 4-3. The left to the right panels represent models with (a) Type-I torque and the gas drag, (b) the gas drag, and (c) the Keplerian gas drag, respectively. The two Gaussian functions within each panel are approximations of the histogram using a mixture Gaussian distribution with 2 clusters. The legends within each panel represent the parameters of the mixed Gaussian distribution.

Supplementary Information IV: The PDM criteria

Minton and Levison (2014) [14] identified five criteria that have to be satisfied simultaneously in order for planets to migrate due to PDM. Here we briefly explain these criteria and we check whether our results satisfy these criteria.

The five criteria are as follows:

1. Mass ration criterion:
   The amount of planetesimals within $5r_{\mathrm{Hill}}$ of the planet at $a_{\mathrm{p}}$ need to be larger than $3M_{\mathrm{p}}$.

2. Mass resolution criterion:
   The mass ratio between the planet and the planetesimal should be larger than 100.

3. Disk eccentricity criterion:
   The eccentricity of the planetesimal disk surrounding the planet has to be less than $5\times$ the reduced hill ($r_{\mathrm{Hill}}/a_{\mathrm{p}}$).

4. Crowded criterion:
   If there are two or more planets near each other in the disk, the larger and the smaller planets’ mass ratio needs to be greater than 10, for the larger planet to migrate by PDM.

5. Growth timescale criterion:
   The planet’s growth timescale, $\tau_{\mathrm{grow}}$, needs to be longer than the planet’s migration timescale, $\tau_{\mathrm{mig}}$. This is related to criterion 4, which implies that the bodies surrounding the planet should not become larger than the planet itself.

In our study, we examine the trajectory of a single planet in a uniform background. Thus, we check the first three criteria only, since the last two are related to the cases where there are multiple planets in the disk. Fig. 9 shows the time evolution of planetary orbits in models 3-1 to 3-5 and 4-1. These five trials were selected from one trial in each model where outward PDM occurred. Below, we examine the first three criteria for these five trials.

Figure 9: The time evolution of the semi-major axes of protoplanets in models 3-1 to 3-5 and 4-1. The thick solid line (black), thin solid line (purple), dashed line (red), dash-dotted line (green), dotted line (blue), and double-dotted dashed line (orange) represent the results from models 3-1, 3-2, 3-3, 3-4, 3-5, and 4-1, respectively.

The mass ratio criterion

Figure 10: The time evolution of the total planetesimal mass within 5 $r_{\mathrm{Hill}}$ outside the planet normalized by the mass of the planet in models 3-1 to 3-5 and 4-1. The thick solid line (black), thin solid line (purple), dashed line (red), dash-dotted line (green), dotted line (blue), and double-dotted dashed line (orange) represent the results from models 3-1, 3-2, 3-3, 3-4, 3-5, and 4-1, respectively. The thick-dashed black line represents $m_{\mathrm{tot}}/M_{\mathrm{p}}=1/3$, which is the critical value for criterion 1.

This criterion states that for monotonous PDM to occur, the planet’s mass $M_{\mathrm{p}}$ must be within a factor of 3 of the total mass of planetesimals $m_{\mathrm{tot}}$, which are within 5 $r_{\mathrm{Hill}}$ outside the planet. To check whether the mass ratio criterion is satisfied or not, we plot the total mass of planetesimals in front of the planet which exhibits outward PDM in Fig. 10. The figure shows the time evolution of the total planetesimal mass within 5 $r_{\mathrm{Hill}}$ outside the planet normalized by the mass of the planet. The thick-dashed black line represents $m_{\mathrm{tot}}/M_{\mathrm{p}}=1/3$, which is the critical value for this criterion. Fig. 10 indicates that this criterion is met in all models.
The mass resolution criterion

Figure 11: The time evolution of the planet mass normalized by the averaged planetesimal mass within 5 $r_{\mathrm{Hill}}$ outside the planet in models 3-1 to 3-5 and 4-1. The thick solid line (black), thin solid line (purple), dashed line (red), dash-dotted line (green), dotted line (blue), and double-dotted dashed line (orange) represent the results from models 3-1, 3-2, 3-3, 3-4, 3-5, and 4-1, respectively. The thick-dashed black line represents $M_{\mathrm{p}}/m=100$, which is the critical value for criterion 2.

This criterion tells us that the mass ratio between the planet and the planetesimal should be larger than 100 for the monotonous PDM to occur. To check this condition, we plot the time evolution of the planet mass normalized by the averaged planetesimal mass within 5 $r_{\mathrm{Hill}}$ outside the planet in Fig. 11. According to the figure, only models 3-5 and 4-1 satisfy this criterion. However, we can see that planets in other models also show general trends of migrating outward in Fig. 9. This result indicates that the outward migration of planets can occur at a lower mass ratio than was previously suggested by [14, 15].
The disk eccentricity criterion

Figure 12: The time evolution of the RMS eccentricity of planetesimals within 5 $r_{\mathrm{Hill}}$ outside the planet normalized by the reduced Hill in models 3-1 to 3-5 and 4-1. The thick solid line (black), thin solid line (purple), dashed line (red), dash-dotted line (green), dotted line (blue), and double-dotted dashed line (orange) represent the results from models 3-2, 3-3, 3-4, 3-5, and 4-1, respectively. The thick-dashed black line represents $\langle e^{2}\rangle^{1/2}\times\left(\frac{r_{\mathrm{Hill}}}{a_{\mathrm{p}}}\right)^{-1}=5$, which is the critical value for criterion 3.

This criterion states that the RMS eccentricity of planetesimals has to be less than $5\times$ the reduced Hill to maintain the constant migration rate. To see this condition, we plot the time evolution of the RMS eccentricity of planetesimals within 5 $r_{\mathrm{Hill}}$ outside the planet normalized by the reduced Hill in models 3-1 to 3-5 and 4-1 in Fig. 12. We can see that the planets in all models meet this criterion.

In models 3-5 and 4-1, the planets satisfy all three criteria. On the other hand, in models 3-1 to 3-4, the planets do not satisfy criterion 2. However, even in models 3-1 to 3-4 with mass ratios less than 100, there are cases where the planets exhibit a relatively monotonic outward-PDM trend. This trend indicates that the potential for planets to move actively by PDM within the protoplanetary disk even during the early protoplanet stage when the mass ratio between planetesimals and planets is not significantly high.

Supplementary Information V: The effect of perturbations from planetesimals on PDM

A planet embedded in the planetesimal disk migrates by gravitationally scattering surrounding planetesimals. Consequently, the migrating planet is influenced by the two-body relaxation arising from interactions with these planetesimals. The impact of the two-body relaxation appears as a random walk in planetary migration. Here we examine the influence of the two-body relaxation from individual planetesimals on a planet’s PDM.

Figure 13: The time evolution of the planetary semi-major axis (upper panel) and its time derivative (lower panel) for one run each from models 3-1 to 3-5 and 4-1 (from upper left to lower right). The red lines represent smoothed planetary semi-major axes and their time derivatives using the Savitzky-Golay (SG) method [22]. The green and the purple dashed lines represent the semi-major axis $\pm$0.05 au from the initial planetary semi-major axis, and we define these lines as the upper (or lower) boundary where PDM begins.

Fig. 13 shows the time evolution of the planetary semi-major axis (upper panel) and its time derivative (lower panel) for one run each from models 3-1 to 3-5 and 4-1. In the figure, the red lines represent smoothed planetary semi-major axes and their time derivatives using the Savitzky-Golay (SG) method [22]. To smooth the data, we used an SG method with a second-degree polynomial for data points spanning $4\times 10^{2}$ years before and after each point. A timescale of $4\times 10^{2}$ years is selected due to the synodic periods associated with closed encounters. Specifically, this pertains to encounters between a $0.5$ Earth-mass planet in a circular orbit near 3 au and planetesimals situated at a distance of and planetesimals located at a distance of $1~{}r_{\mathrm{Hill}}$ from the planet, which approximate duration of $4\times 10^{2}$ years. The blue lines in the lower panels represent the boundary where the time derivative of the planetary semi-major axis becomes zero. Thus, the points where the blue and red lines intersect mark the moments when the direction of the planet’s migration flips. Comparing the results of each run, it is evident that runs with a larger number of particles result in fewer flips in planetary migration. In addition, the results from our total of 350 runs also indicate a proportional decrease in the flip count with $1/N$, as shown in Fig. 5. This $1/N$ dependence of planetary migration flips can be explained by two-body relaxation between a planet and planetesimals. Specifically, the variance in velocity changes due to two-body relaxation between a planet and planetesimals is proportional to $1/N$ in the limit where the velocity of a planet becomes zero [21].

Indeed, the increase in velocity dispersion of planets caused by two-body relaxation between planets and planetesimals slows down the monotonic planetary migration driven by PDM (Fig. 6). However, the overall effect of PDM is stronger than two-body relaxation, resulting in planets migrating relatively monotonously in one direction.

Acknowledgments This work was supported by MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (Structure and Evolution of the Universe Unraveled by Fusion of Simulation and AI; Grant Number JPMXP1020230406) and used computational sources of supercomputer Fugaku provided by the RIKEN Center for Computational Science (Project ID: hp230204). Test simulations in this paper were also carried out on a Cray XC50 system at the Centre for Computational Astrophysics (CfCA) of the National Astronomical Observatory of Japan (NAOJ).

Declarations

• Funding
  This work was supported by MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (Structure and Evolution of the Universe Unraveled by Fusion of Simulation and AI; Grant Number JPMXP1020230406)

• Conflict of interest/Competing interests (check journal-specific guidelines for which heading to use)
  The authors declare no competing interests.

• Ethics approval
  Not applicable

• Consent to participate
  Not applicable

• Consent for publication
  Not applicable

• Availability of data and materials
  The data that support the findings of this study are available from the the corresponding author on request.

• Code availability
  The codes used to generate these results are available from the corresponding author on request. The source code for GPLUM is hosted on GitHub and is available for download from https://github.com/YotaIshigaki/GPLUM under the MIT license.

• Authors’ contributions
  Tenri Jinno: Conceptualization, Methodology, Software, Formal analysis, Investigation, Data Curation, Writing - Original Draft, and Supervision
  Takayuki R. Saitoh: Methodology, Software, and Writing - Review & Editing
  Junichiro Makino: Methodology, Software, Writing - Review & Editing, Project administration, and Funding acquisition
  All authors have agreed to the published version of the manuscript.

If any of the sections are not relevant to your manuscript, please include the heading and write ‘Not applicable’ for that section.