# The influence of cold Jupiters in the formation of close-in planets. I.
planetesimal transport

Sergio Best Instituto de Astrofísica, Pontificia Universidad Católica de
Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile  Millennium
Institute of Astrophysics MAS, Nuncio Monsenor Sotero Sanz 100, Of. 104,
Providencia, Santiago, Chile  Antranik A. Sefilian Alexander von Humboldt
Postdoctoral Fellow.  Astrophysikalisches Institut und Universitätssternwarte,
Friedrich-Schiller-Universität Jena, Schillergäßchen 2–3, D-07745 Jena,
Germany  Cristobal Petrovich Instituto de Astrofísica, Pontificia Universidad
Católica de Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile
Millennium Institute for Astrophysics, Chile

(Received XXX; Accepted YYY)

###### Abstract

The formation of a cold Jupiter (CJ) is expected to quench the influx of
pebbles and the migration of cores interior to its orbit, thus limiting the
efficiency of rocky planet formation either by pebble accretion and/or orbital
migration. Observations, however, show that the presence of outer CJs ($>1$ au
and $\gtrsim 0.3M_{\rm Jup}$) correlates with the presence of inner Super
Earths (at $<1$ au). This observation may simply be a result of an enhanced
initial reservoir of solids in the nebula required to form a CJ or a yet-to-
be-determined mechanism assisted by the presence of the CJ. In this work, we
focus on the latter alternative and study the orbital transport of
planetesimals interior to a CJ subject to the gravity and drag from a
viscously-evolving gaseous disk. We find that a secular resonance sweeping
inwards through the disk gradually transports rings of planetesimals when
their drag-assisted orbital decay is faster than the speed of the resonance
scanning. This snowplow-like process leads to large concentration (boosted by
a factor of $\sim 10-100$) of size-segregated planetesimal rings with aligned
apsidal lines, making their expected collisions less destructive due to their
reduced velocity dispersion. This process is efficient for a wide range of
$\alpha-$disk models and Jovian masses, peaking for $\sim 1-5M_{\rm Jup}$,
typical of observed CJs in radial velocity surveys. Overall, our work
highlights the major role that the disk’s gravity may have on the orbital
redistribution of planetesimals, depicting a novel avenue by which CJs may
enhance the formation of inner planetary systems, including super-Earths and
perhaps even warm and hot Jupiters.

planets and satellites: formation — protoplanetary disks — planet-disk
interactions

##  1 Introduction

Figure 1: Schematic cartoon describing the basic setup of our simulations. A
giant planet carves a gap in a depleting gaseous disk with a planetesimal disk
interior to its orbit. As the disk disperses, a secular resonance that excites
the planetesimals’ eccentricities moves inwards carrying them along due to gas
drag, accumulating material in the form of size-segregated planetesimal rings.
The final surface density of planetesimals inside $\sim 0.2$ au is then
largely enhanced.

###  1.1 Motivation

Recent exoplanet observations show that the presence of cold giant planets may
increase the probability of finding an inner super-Earth in the same system
(Zhu & Wu, 2018; Bryan et al., 2019). These works estimate that the
probability of having at least one super-Earth given a cold giant in the
exoplanetary system is $\sim 50-100\%$ (within 1$\sigma$; more recently
Rosenthal et al. 2022 also found such a correlation, albeit less significant).
Thus, understanding the role that cold Jovians have at promoting (or
inhibiting) the formation of inner super-Earths is a key step forward in
building a complete picture of planet formation.

According to current formation models that rely on the transport of rocky
material (e.g., pebbles or cores) from the outer regions to assemble super-
Earths in the inner au, the presence of a cold Jupiter is seen as an obstacle,
and thus generally predict an anti-correlation (e.g. Schlecker, M. et al.
2021). Indeed, in the pebble accretion model, planets form from planetesimal
cores that accrete pebbles until they become sufficiently large (Johansen &
Lacerda 2010; Ormel & Klahr 2010). However, if a giant planet is formed early
on in the outer regions of the protoplanetary disk, it would quench the flux
of pebbles to the inner system. This is because the typical core of the Jovian
reaches the pebble isolation mass (Lambrechts et al., 2014). This blockage
should halt the accretion of pebbles by the inner planetesimals, limiting the
formation of other planets interior to the initial gas giant as the remaining
pebbles would just drift towards the star too quickly. Also, the same pressure
bump forming the gas giant would quench the pebble flux (Izidoro et al.,
2021).

Instead of pebbles, the influx of solids may be dominated by cores undergoing
type-I migration that park into super-Earth orbits. However, the formation of
a giant planet would also prevent their migration, giving rise to an anti-
correlation between the two populations (Izidoro et al., 2015).

Finally, another possible explanation for the observed correlation is that the
presence of the cold Jupiter has no significant effect on the formation of a
super-Earth. Instead, the potential correlation simply emerges as a
consequence of birth disk conditions: disks with more solid material inside
$\sim 10$ au may promote the formation of both populations (see e.g. the
discussion in Zhu et al. 2018). This can happen either by having a more
massive disk or having a higher dust-to-gas ratio. Although there is still
debate about the exact distribution of material in the disk (Chiang &
Laughlin, 2013), the slope in the distribution of solids in multiple planet
systems seems to be significantly different from what models of gaseous disks
can explain. This points towards a rearrangement of the solid material
(Raymond & Cossou, 2014).

In this work, we explore yet another possibility where the same giant planet
blocking the pebbles or inwardly migrating planet cores could transport
copious amounts of planetesimals from $\gtrsim 1$ au to the inner regions.
This would imply that terrestrial planets in such a scenario grow by
planetesimal accretion, possibly assisted by the higher density of solids
required to form the Jovian.

###  1.2 Sweeping secular resonances

The mechanism that we invoke to transport planetesimals inwards is based on
the sweeping of secular resonances to drive eccentricity excitation and gas
drag to damp them, causing orbital decay. The sweeping arises from the
commensurability between the apsidal precession rates of the Jovian planet and
that of the planetesimals (Heppenheimer 1980; Ward 1981a), a process that has
been studied in various contexts, including the dynamics of asteroids in the
Solar System (Lemaitre & Dubru, 1991; Nagasawa et al., 2000; Zheng et al.,
2017a, b; Gong et al., 2019) and the excitation of planetary eccentricities
and/or inclinations (Hahn & Malhotra, 1999; Nagasawa et al., 2003; Petrovich
et al., 2019; Toliou et al., 2019).

This process can be understood as follows (see also Figure 1 for an
illustration). The gap-opening Jovian planet will undergo apsidal precession
with a positive frequency due to the disk’s gravity (Ward, 1981a). On the
other hand, non-gap-opening planets/planetesimals do so at a (generally)
negative rate due to the disk (Ward, 1981a) and positively due an external
perturber like in the case of the planetesimals with the Jovian(s). The
depletion of the disk can lead to instances when the frequencies match at a
given orbital radius, driving the excitation of the planetesimals’
eccentricities (see Figure 2 for an illustration). These eccentricities are
damped by the drag from the gas, pushing the planetesimals inwards, plowing
them through the inner system as the resonance moves inward as the disk
disperses. This can effectively transport a large percentage of solid material
available in the region.

A major goal of our work is to understand these dynamics for realistic disks
by accounting for gaps, inner cavities, and variable viscosities; expanding
previous treatments which typically use uniform power-law disks.

###  1.3 Structure

Our work is organized as follows. In Section 2 we describe the general setup
of our model, and present equations governing the dynamics of planetesimals.
In Section 3 we then present results showing the transport of material due to
our proposed mechanism, while the impact the alignment of planetesimals has on
their velocity dispersion is presented in in Section 4. Finally we discuss
some of the implications our findings may have for planet formation in Section
5 before giving our conclusions in Section 6.

##  2 Methods

In this section, we describe a semi-analytic model to analyze the long-term
dynamical evolution of planetesimals that are embedded within a viscously
evolving protoplanetary disk and perturbed by a single planet.

The planet, the disk, and the planetesimals within it are assumed to be
coplanar orbiting a central star of mass $M_{c}=1M_{\sun}$. We characterize
the orbit of a given planetesimal by its semimajor axis $a_{\rm p}$,
eccentricity $e_{\rm p}$, and longitude of pericenter $\varpi_{\rm p}$.
Similarly, orbital elements subscripted with J refer to the accompanying giant
planet.

In order to follow the evolution of the planetesimals we need to model their
interaction with the giant planet. Additionally, we need to include the
effects of the gaseous disk component. This component produces gas drag onto
the planetesimals, and drives apsidal precession of the planetesimals and the
giant. The following equations describe all the interactions considered:

$\displaystyle\frac{de_{\rm p}}{dt}$ $\displaystyle=$
$\displaystyle\left(\frac{de_{\rm p}}{dt}\right)_{\rm
drag}+\left(\frac{de_{\rm p}}{dt}\right)_{\rm CJ},$ (1)
$\displaystyle\frac{d\Delta\varpi}{dt}$ $\displaystyle=$
$\displaystyle\left(\frac{d\varpi_{\rm p}}{dt}\right)_{\rm
disk}+\left(\frac{d\varpi_{\rm p}}{dt}\right)_{\rm
CJ}-\left(\frac{d\varpi_{\rm J}}{dt}\right)_{\rm disk},$ (2)
$\displaystyle\frac{da_{\rm p}}{dt}$ $\displaystyle=$
$\displaystyle\left(\frac{da_{\rm p}}{dt}\right)_{\rm drag},$ (3)

where $\Delta\varpi=\varpi_{\rm p}-\varpi_{\rm J}$.

In subsequent subsections, we describe in detail each of these interactions,
as well as the evolution of the gaseous component of the disk and the initial
distribution of planetesimals.

###  2.1 Disk evolution

We consider a gaseous protoplanetary disk that evolves as an $\alpha$
accretion disk, without considering the influence of any other physical
effects (e.g. due the giant planet, photoevaporation winds, etc). In our most
general case, two different viscosities are considered (to the inner and outer
side of the giant planet), which is further described in Section 2.1.2. We
compute the time evolution of the gas surface density $\Sigma_{\rm gas}$ using
the classical 1D diffusion equation (Pringle, 1981):

$\displaystyle\frac{\partial\Sigma_{\rm gas}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{1}{r}\frac{\partial}{\partial
r}\left[3r^{1/2}\frac{\partial}{\partial r}(\nu\Sigma_{\rm
gas}r^{1/2})\right],$ (4) $\displaystyle\nu$ $\displaystyle=$
$\displaystyle\alpha c_{s}H,$ (5)

and account for a gap profile produced by the giant planet a posteriori; see
Section 2.1.3 for details. Here, $\alpha$ is the viscosity parameter, $c_{s}$
is the sound speed, and H is the scale height, which in our case of a locally
isothermal disk, correspond to $c_{s}\propto r^{-1/4}$ and $H\propto r^{5/4}$
(with $H\approx 0.04$ au at $1$ au).

####  2.1.1 Initial profile

We consider a gas surface density which initially has a power-law profile with
a cutoff at some outer radius:

$\Sigma_{\rm gas}(r)=\Sigma_{\rm gas,0}\left(\frac{r}{1\rm
au}\right)^{-\gamma}\exp\left[-(r/r_{\rm cut})^{2-\gamma}\right],$ (6)

where $\Sigma_{\rm gas,0}$ is chosen so that the total mass of gas is 1% of
the stellar mass. Motivated by observations, we adopt a $\gamma$ index of 1
here and $r_{\rm cut}$ to be 25 au. Here, we note that variations in these
parameters do not affect the evolution of planetesimals significantly as the
shape of the profile reaches a steady state irrespective of the initial value
of $\gamma$ (see Appendix A). Additionally, the planetesimals/planet’s orbital
precession depends more on the total mass of the outer disk and not so much on
its outer edge (see Appendix C).

####  2.1.2 Viscosity

The evolution of the disk’s surface density depends on the gas viscosity; see
e.g. Eq. (4). This directly impacts the rate of depletion of a region; namely,
the higher the viscosity, the quicker the region depletes its gas. At an
iceline or the boundary of a dead zone, it is generally thought that the gas
viscosity changes111In Fig. 11 we see a surface density transition at this
location., serving as a dust trap where a giant planet can form early on
(Guilera et al., 2020). Following this line of reasoning, we use the same
position for both the viscosity transition and the position of the giant
planet (namely, $a_{\rm J}=3$ au, unless stated otherwise). In this work we
will refer to these viscosities as $\alpha_{\rm in}$ and $\alpha_{\rm out}$
(inside and outside the boundary, respectively). This gives us the freedom of
separating the depletion of the inner and outer disks which, as we shall see,
has important consequences for the evolution of the secular resonance location
(Appendix C).

####  2.1.3 Density gap opened by the Jovian

As the planet accretes the surrounding gas, it will open a gap in the gas disk
if it is sufficiently massive (see the red region of panels B1 to B3 in Fig.
3). We need to keep track of this effect, as the drag force and precession
rate of the planetesimals depend on the gas surface density. In principle,
this can be done self-consistently by adding a torque at the position of the
giant planet to the viscous evolution equation (Eq. (4)). For simplicity,
however, we model the gap using the profile and code provided in Duffell
(2020) for each time step independently. For reference, an approximate
expression for the gap width $\Delta_{\rm gap}$ is provided by Equation (4) of
Kanagawa et al. (2016):

$\frac{\Delta_{\rm gap}}{a_{\rm J}}=0.41\left(\frac{M_{\rm
J}}{M_{c}}\right)^{1/2}h^{-3/4}\alpha^{-1/4},$ (7)

where $M_{\rm J}$ is the mass of the giant, $h$ is the disk’s aspect ratio at
the position of the Jovian, and the viscosity of the inner disk ($\alpha_{\rm
in}$) is considered for the gap. We use this simplified expression for our
model in Appendix C.

###  2.2 Drag force

The planetesimals embedded in the disk feel a drag force when the gas in their
vicinity moves at a non-Keplerian speed. We can quantify this difference in
velocity with the factor $\eta(r)$. From hydrostatic equilibrium, a disk with
volumetric gas density $\rho\propto r^{-\Gamma}$ and temperature $T\propto
r^{-\beta}$ has an azimuthal velocity given by (e.g. Adachi et al., 1976):

$v_{\rm gas}(r)=v_{\rm K}(r)\left[1-2\eta(r)\right]^{1/2}\simeq v_{\rm
K}(r)\left[1-\eta(r)\right],$ (8)

with

$\eta(r)=\frac{\pi(\Gamma+\beta)}{16}\frac{c_{m}^{2}}{v_{K}^{2}}=\frac{(\Gamma+\beta)}{2}h^{2},$
(9)

where $v_{K}$ is the Keplerian velocity, and we have used the fact that
$c_{s}/v_{K}$ is equal to the aspect ratio of the disk $h$, and the mean
thermal speed is $c_{m}^{2}=\frac{8}{\pi}c^{2}_{s}$.

In our models, we use $\beta=1/2$ and we can obtain $\Gamma$(r)
locally222Computing $\Gamma$(r) locally is important because the surface
density quickly deviates from a power law due to the disk evolution,
especially near the edges of the disk, the gap of the planet, and the
viscosity transition. by calculating:

$\Gamma(r)=-\frac{r}{\rho}\frac{d\rho}{dr}.$ (10)

By averaging the drag force over one orbit, Adachi et al. (1976) find the
following approximate expressions for the migration of the planetesimals:

$\displaystyle\left(\frac{da_{\rm p}}{dt}\right)_{\rm drag}$
$\displaystyle=-\dfrac{2a_{\rm p}}{\tau_{0}}\sqrt{\dfrac{5}{8}e_{\rm
p}^{2}+\eta^{2}}\left[\eta+\left(\dfrac{\Gamma}{4}+\dfrac{5}{16}\right)e_{\rm
p}^{2}\right],$ (11) $\displaystyle\left(\frac{de_{\rm p}}{dt}\right)_{\rm
drag}$ $\displaystyle=-\dfrac{e_{\rm p}}{\tau_{0}}\sqrt{\dfrac{5}{8}e_{\rm
p}^{2}},$ (12)

where

$\displaystyle\tau_{0}$ $\displaystyle=\left(\dfrac{\pi C_{\rm
D}}{2m}s^{2}\rho v_{K}\right)^{-1}.$

where $s$ is the radius of the planetesimal, $C_{D}$ the drag coefficient, and
$\rho$ the density of planetesimals which we take to be 1 g/cm3. Note that we
are not considering the inclinations of the planetesimals in our work, as
these should be suppressed due to the Jovian’s inclination being close to
zero.

###  2.3 Apsidal precession due to the disk gravity

Next, we consider the gravitational effects of the protoplanetary disk on both
the planetesimals embedded within and the planet. Given that the disk is taken
to be axisymmetric, it simply causes the longitudes of periapsis of the
planetary and planetesimal orbits to precess.

Since these objects are surrounded by the disk, the gas can be infinitely
close, leading to a numerical (not physical) problem which causes a divergence
in the gravitational force (see Sefilian & Rafikov, 2019, for a detailed
discussion). To overcome this problem, we soften the gravitational potential
using the softening prescription of Hahn (2003) so that (Sefilian & Rafikov,
2019):

$\dfrac{1}{|\textbf{r}-\textbf{r'}|}=\left[r^{2}+r^{\prime
2}-2rr^{\prime}\cos(\theta)+h^{2}(r^{2}+r^{\prime 2})\right]^{-1/2},$ (13)

Accordingly, the usual Laplace coefficients,

$\displaystyle b_{\rm s}^{\rm(m)}(\alpha)$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{\rm
0}^{\pi}\frac{\cos(m\theta)d\theta}{\big{(}1+\alpha^{2}-2\alpha\cos\theta\big{)}^{s}},$
(14)

are modified which now read as follows:

$\displaystyle\mathcal{B}_{\rm s}^{\rm(m)}(\alpha)$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{\rm
0}^{\pi}\frac{\cos(m\theta)d\theta}{\big{(}1+\alpha^{2}-2\alpha\cos\theta+(1+\alpha^{2})h^{2}\big{)}^{s}}.$

Here, we note that the aspect ratio is taken to be a distance-independent
constant in Hahn (2003). In our model, however, we assume that $h$ varies with
semimajor axis. As a consequence, we use its value at the position of the
object. This approximation should be valid since the smoothening only affects
the potential in the neighbourhood of the object.

Then, we can calculate the precession due to an axisymmetric disk using the
following equation from Sefilian & Rafikov (2019):

$\displaystyle\left(\frac{d\varpi_{i}}{dt}\right)_{\rm
disk}=\frac{2G}{n_{i}a_{i}^{3}}$ $\displaystyle\bigg{[}\int_{\rm r_{\rm
in}}^{\rm a_{i}}\mu(r)\phi\left(\frac{r}{a_{\rm i}}\right)dr$
$\displaystyle+\int_{\rm a_{i}}^{\rm r_{\rm outbox}}\frac{a_{\rm
i}}{r}\mu(r)\phi\left(\frac{a_{\rm i}}{r}\right)dr\bigg{]},$ (16)

where

$\displaystyle\phi(x)=\frac{1}{8}x\left[\mathcal{B}_{\rm
3/2}^{(1)}(x)-3xh^{2}(2+h^{2})\mathcal{B}_{\rm 5/2}^{(0)}(x)\right],$ (17)

where i = {p, J}, $n_{i}=\sqrt{G(M_{c}+m_{i})/a_{i}^{3}}$ its mean motion,
$r_{\rm outbox}$ the outer limit of the simulation box (1000 au) and $\mu(r)$
is the gas mass density per unit semimajor axes, $\mu(r)=2\pi r\Sigma_{\rm
gas}$. Equation (16) generally gives a negative value for the planetesimals
(subscript p) because they are embedded in the gas, but a positive value for
the giant planet (subscript J) because of the gap the latter opens (see e.g.
Ward, 1981b; Sefilian & Rafikov, 2019).

It is worth noting that the singularity discussed above does not arise when
considering objects orbiting outside the disk, i.e., where $\Sigma_{\rm
gas}(r)=0$, such as the planet in our model. We also point out that there are
analytic frameworks that allow for the computation of the disk-induced
precession without introducing any softening parameter (see e.g. Davydenkova &
Rafikov (2018)). Despite this, however, we opted to adopt the softening
prescription of Hahn (2003), as the resulting softening parameter $h$ is the
same as the disk’s aspect ratio. This has the added effect of accounting for
the gaseous disk not being razor-thin which is a more realistic
representation.

###  2.4 Planetesimal-Jupiter interaction

We now consider the orbital evolution of the planetesimals due to the Jovian.
Making use of the classical disturbing function $\mathcal{R}$ expanded to
second order in eccentricities (e.g. Eq. (7.8) in Murray & Dermott (2000)),
one can write:

$\displaystyle\left(\frac{de_{\rm p}}{dt}\right)_{\rm CJ}$
$\displaystyle=\dfrac{1}{n_{\rm p}a_{\rm p}^{2}e_{\rm
p}}\dfrac{\partial\mathcal{R}}{\partial\varpi_{\rm p}}=A_{\rm pJ}e_{\rm
J}\sin(\varpi_{\rm p}-\varpi_{\rm J}),$ (18)
$\displaystyle\left(\frac{d\Delta\varpi}{dt}\right)_{\rm CJ}$
$\displaystyle=\dfrac{1}{n_{\rm p}a_{\rm p}^{2}e_{\rm
p}}\dfrac{\partial\mathcal{R}}{\partial e_{\rm p}}=A_{\rm pp}+A_{\rm
pJ}\dfrac{e_{\rm J}}{e_{\rm p}}\cos(\varpi_{\rm p}-\varpi_{\rm J}),$ (19)

where $n_{\rm p}$ is the planetesimal’s mean motion. In Equations (18) and
(19), the coefficients $A_{\rm pp}$ and $A_{\rm pJ}$ are given by:

$\displaystyle A_{\rm pp}$ $\displaystyle=$ $\displaystyle\dfrac{n_{\rm
p}}{4}\dfrac{M_{\rm J}}{M_{\rm c}+m_{\rm p}}\left(\frac{a_{\rm p}}{a_{\rm
J}}\right)^{2}b_{\rm 3/2}^{(1)}(a_{\rm p}/a_{\rm J}),$ (20) $\displaystyle
A_{\rm pJ}$ $\displaystyle=$ $\displaystyle-\dfrac{n_{\rm p}}{4}\dfrac{M_{\rm
J}}{M_{\rm c}+m_{\rm p}}\left(\frac{a_{\rm p}}{a_{\rm J}}\right)^{2}b_{\rm
3/2}^{(2)}(a_{\rm p}/a_{\rm J}),$ (21)

see e.g. Murray & Dermott (2000).

###  2.5 Secular Resonance

Now that all the sources for the apsidal precession of the Jovian and
planetesimals are defined, we proceed to describe when a secular resonance
takes place. Loosely speaking, this kind of resonance occurs when we have a
matching of the apsidal precession rates (in our case between the planet and
some planetesimals). This would mean that Eq. (2) is close to zero in the
absence of the gas drag which aligns (or anti-aligns) the orbits. A
planetesimal in this secular resonance will increase its eccentricity up to a
maximum value (see Eq. (B2) for an estimate). Once the planetesimal acquires
enough eccentricity, its migration will be significant (Eq. (11)). For
illustrative purposes, in Fig. 2 we plot the precession rate of the
planetesimals in the fiducial case evolving in time as a function of semimajor
axis. Also, overlapping, is the precession rate of the Jovian (in green) which
decreases in time as the disk depletes. We show the matching of the
frequencies as a circle at each time. As one can see, the resonance starts
close to $1.5$ au and moves inwards towards the inner edge at $0.1$ au. In the
next section (Section 3), we shall demonstrate how this sweeping of the
resonance could transport material.

###  2.6 Planetesimal distribution

Finally, with the basic elements of the model in place, we describe how the
planetesimals themselves are distributed throughout the disk. To allocate the
number of planetesimals at the beginning of the simulation, we consider 200
logarithmically spaced rings spanning from 0.1 au to 1.5 au. Each ring
contains planetesimals such that their combined mass is 1% of the
corresponding gas mass for the ring. In turn, the size distribution for
planetesimals follows a power law:

$\frac{dN}{dM}\propto M^{\rm-p},$ (22)

where N is the number of planetesimals in each of the 50 logarithmically
spaced mass bins M (from 1 km to 100 km with a density of 1 g/cm3). Looking at
Eq. (22), it can be easily seen that with $p=2$, all sizes of planetesimals
contribute to the mass budget with the same total mass. Values below 2 are top
heavy (meaning, the bigger planetesimals contribute more to the mass budget),
while values above 2 are weighted towards the smaller planetesimals. In what
follows we use $p=2$ as to not prioritize any size of planetesimals. The exact
distribution of planetesimals will not significantly affect our results as we
will discuss the evolution of each size of planetesimals in what follows.

Our choice for the distribution in semimajor axis of the planetesimals is 0.1
au to 1.5 au. The lower limit is simply the inner edge of the disk and the
upper limit is introduced because planetesimals close to the planet would be
evactuated by the mean motion resonances (MMRs). We discuss this further in
Section 5.3.3.

Figure 2: Precession rate of the planetesimals due to both the disk and planet
gravity as a function of semimajor axis. The different lines are for different
times as described by the colorbar. The green lines represent the precession
rate of the giant planet which decreases in time as the disk depletes. For
reference, the resonance position as a function of time is marked with a
circle symbol corresponding to the same colorbar. The resonance position moves
from approximately 1.5au to 0.2au at the end of the 10 Myrs. A second
resonance can also be seen at the inner edge, but this one spans a narrow
range from 0.1 au to 0.15 au. Figure 3: Evolution of the planetesimals and the
gaseous disk. Row A shows the eccentricity of planetesimals as a function of
semimajor axis for a sub-sample of sizes ($s=1$ km, $10$ km, and $100$ km.).
For reference, the corresponding forced eccentricity solutions obtained by
solving Equations (1) and (2)) are also shown in red, green, and blue dashed
lines, respectively. Row B shows the gas surface density in blue with the gap
carved by the planet in red. The solid black line represents the surface
density of planetesimals, and for reference the initial profile is shown with
a dashed black line. Row C shows the cosine of the relative apsidal angle of
planetesimals with respect to the precessing giant planet. For reference, the
forced apsidal angles are also shown using the same color scheme as Row A.
Each column corresponds to a different time: from left to right, $t=0.01$,
$0.5$, and $10$ Myrs. See Section 3.1 for further details. An animated version
of this figure is available online.

##  3 Results: Transport

###  3.1 Fiducial case

Table 1: Initial conditions of the disk and planet for our simulations.  Parameter | Fiducial value | Other values
---|---|---
$M_{\rm J}[M_{\rm Jup}]$ | 3 | [0.1 – 7]
$a_{\rm J}[au]$ | 3 | 5
$e_{\rm J}$ | 0.05 | 0.01
$\alpha_{\rm in}$ | $10^{-3}$ | [$10^{-4}$ – $10^{-2}$], [$3\times 10^{-3}$]
$\alpha_{\rm out}$ | $10^{-4}$ | [$10^{-5}$ – $10^{-3}$],[$3\times 10^{-3}$]

We simulate the evolution of planetesimals ranging from 1 km to 100 km,
keeping track of their semi major axes, eccentricities, and apsidal angles
(with respect to the Jovian) for $10$ Myrs; see Equations (1)–(3). All
parameters considered for the fiducial case are listed in Table 1 and in text
form here. In the fiducial case, we invoke a giant planet of 3 Jupiter masses
(a typical cold giant planet) at 3 au (the same position as the snowline), and
$e_{\rm J}=0.05$ (which is reasonable as the planet’s eccentricity can grow up
to the order of the local aspect ratio (Papaloizou & Larwood, 2000; Duffell &
Chiang, 2015)). The inner and outer viscosity parameters for the disk are
$\alpha_{\rm in}=10^{-3}$ and $\alpha_{\rm out}=10^{-4}$. We start the
simulation with a gas disk mass of 1% of a solar mass and a solid surface
density of 1% of the gas disk mass. As planetesimals do not interact among
themselves, the amount of available solid material can be easily scaled up or
down later and we shall express our results in fraction of transported mass.

The results of our fiducial simulation are summarized in Figure 3, where we
plot the radial profiles of the planetesimal eccentricities (row A), surface
densities of the gas and planetesimals (row B), and planetesimal apsidal
angles (relative to the precessing planet; row C). Results are shown at three
different times: from left to right, $t=0.01$, $t=0.5$, and $t=10$ Myr,
respectively. For clarity purposes, only a sub-sample of three planetesimal
sizes are shown: namely, 1, 10, and 100 km. Starting from panel A1 we see the
planetesimals oscillating around the forced eccentricity profiles (dashed
lines). These forced eccentricities are calculated by simultaneously solving
Eq. (1) and (2) looking for a steady state (which is different for each
planetesimal size). The different terms for these equations are described
throughout Section 2. We consider each time step to be independent, using the
previously calculated gas surface density as described in Section 2.1. As the
resonance location moves towards the star, it carries the planetesimals with
it (panels A1 to A3), thus sweeping up the solid material. Whether the
planetesimals move with the resonance or not is explored in Section 3.2. Note
the corresponding increase in solid surface density (black line) in panels B1
to B3. We will quantify this increase in surface density in Section 3.3.

Some things to further note. First, the largest planetesimals of 100 km depend
on the resonance position to migrate, while the smaller ones of 1 km migrate
almost irrespective of the resonance position, as the planet’s perturbation
already excites eccentricities that are large enough to drive their migration
(see Section 3.2).

Second, in panels C1 to C3, we see planetesimals close to the planet become
aligned with the orbit of the planet (i.e., $\Delta\varpi=0$), while
planetesimals further away become anti aligned (i.e., $\Delta\varpi=\pi$).
This means that there is a transition between these two alignments which
coincides with the resonance position. As a consequence, planetesimals near
the resonance point are not well aligned with each other, especially for
smaller sizes, thus increasing their velocity dispersion (as we will further
discuss in Section 4).

Finally, we also note the occurrence of a second secular resonance near the
inner edge of the disk. This is due to the sharp positive slope of the surface
density at this edge which contributes a positive component to the precession
rate of the planetesimals; see Fig. 2. This additional resonance is not
important for our purposes as it transports very little material (see the
black line at the inner edge of panels B1 and B2).

The fiducial case was chosen because the timescale for the migration of its
planetesimals remains shorter than the timescale for the movement of its
resonance (for the first 10 Myrs). As such, the planetesimals can ride the
wave of the secular resonance which results in a very efficient transport of
material. Despite this, the fiducial case is not the case that transports the
most material in our simulations as some of the planetesimals are left behind
the wave (see Fig. 3 A3). To better understand why this happens, next we will
explore the condition for transporting planetesimals in more detail.

###  3.2 Condition for transport

We can obtain an approximate condition for the effective transport of
material. For a planetesimal to remain ahead of the wave, its migration
timescale $\tau_{\rm mig}$, at the peak of the eccentricity ($e_{\rm peak}$),
must be shorter than the timescale for the movement of the resonance position
$\tau_{\rm res}$. As we show in Appendix B, this translates into the following
criterion:

$\displaystyle\tau_{\rm res}\gtrsim\tau_{\rm
mig}\approx\frac{1}{2}\tau_{0}^{-1/2}\left(A_{\rm pJ}e_{\rm J}\right)^{-3/2}$
$\displaystyle\approx$ $\displaystyle 40{\rm
Myr}\left(\frac{M_{c}}{1M_{\rm\sun}}\right)\left(\frac{3M_{\rm Jup}}{M_{\rm
J}}\right)^{3/2}\left(\frac{\Sigma_{\rm gas}}{10\rm
g/cm^{2}}\right)^{1/2}\left(\frac{0.1}{h}\right)^{1/2}$ (23)
$\displaystyle\left(\frac{s}{100\rm km}\right)^{-1/2}\left(\frac{a_{\rm
J}}{3{\rm au}}\right)^{6}\left(\frac{a_{\rm p}}{0.25{\rm
au}}\right)^{-9/2}\left(\frac{e_{\rm J}}{0.05}\right)^{-3/2}.$

where s is the radius of the planetesimal. In order to illustrate the use of
this expression, in Fig. 4 we observe the transport of material and pay
attention to the time at which the transport stops working. Panel 4A shows
that the large planetesimals (10 - 100 km) need the resonance to move inwards,
while stopping their migration if at any moment the resonance runs past them
(i.e. they no longer fulfill Eq. (3.2)), leaving the planetesimals behind as
can be seen for the 15 Myrs snapshot. At 10 Myrs, the migration rate of the
100 km planetesimals is around 40 Myrs according to Eq. (3.2), which is above
the timescale for the resonance movement ($\sim\Delta t\times a/\Delta
a\sim\rm 5Myrs\times 0.2au/(0.05au)\sim 20Myrs$).

On the other hand, the small planetesimals (1 - 5 km) all migrate regardless
of the position of the resonance because a forced eccentricity below the peak
value suffices to drive efficient migration (see Fig. 4B). Because of our
choice of p in Eq. (22), all sizes contribute the same amount of mass to the
surface density, and so we still see migration of material even when the
resonance moves too quickly (albeit, less efficiently), as can be seen in Fig.
4C.

The difference in the migration timescales for different sizes translates into
a segregation of planetesimals. This segregation helps to keep the
planetesimals better aligned, as when comparing two planetesimals identical in
size, their relative velocity would be zero.

Figure 4: Radial profiles of the surface density distribution of
planetesimals. In each panel, results are shown relative to the initial
profile for different times: 0.1, 0.5, 3, 10, and 15 Myrs. The vertical dashed
lines correspond to the positions of the secular resonance at each considered
timestep. Panel A shows only the big planetesimals from 10 to 100 km. These
planetesimals need to be near the resonance to migrate and when the resonance
moves too quickly, the wave cannot catch up on time. Panel B shows only the
smaller planetesimals from 1 to 5 km. Note that their migration happens more
or less irrespective of the resonance position. Panel C shows all considered
sizes, i.e., from 1 to 100 km.

###  3.3 Amount of transported material

We now quantify the amount of material transported by the sweeping secular
resonance. We do this by simply integrating the surface density of solid
material in places where there is an excess compared to the starting density,
i.e.:

$M_{\rm t}=\int_{\rm r_{\rm in}}^{\rm r_{\rm trans}}\left(\Sigma_{\rm
p}-\Sigma_{0}\right)dA=\int_{\rm r_{\rm trans}}^{\rm r_{\rm
out}}\left(\Sigma_{0}-\Sigma_{\rm p}\right)dA,$ (24)

with $r_{\rm in}$ and $r_{\rm out}$ the limits for the planetesimal
distribution, and $r_{\rm trans}$ (the transition point at $\Sigma_{\rm
p}=\Sigma_{0}$) where we go from having an excess of planetesimals to a
deficit, compared to the initial distribution (see Fig. 4). An approximation
to this equation can be made, where we assume $\Sigma_{\rm p}=0$ for the
rightmost expression of Eq. (24). In other words, we assume all the region
depleted of material (swept by the resonance) has no solid material left and
set $r_{\rm trans}=a_{\rm res}$. By doing this, we find the following
approximation for the fraction of transported mass:

$\displaystyle f\equiv\frac{M_{\rm t}}{M_{0}}$ $\displaystyle\approx$
$\displaystyle\int_{\rm\max(0.1,a_{\rm f})}^{\min(1.5,a_{\rm
i})}\Sigma_{0}dA\bigg{/}\int_{\rm 0.1au}^{\rm 1.5au}\Sigma_{0}dA$ (25)
$\displaystyle\approx$ $\displaystyle\frac{\min(1.5\rm au,a_{\rm
i})-\max(0.1\rm au,a_{\rm f})}{\rm 1.5au-0.1au},$

where $a_{i}$ and $a_{f}$ represent the initial and final positions of the
resonance, respectively, and $M_{0}$ is the initial solid mass in the disk. In
the fiducial case, this fraction is about 0.73 after 10 Myrs.

Next we use Eq. (24) to compare the amount of transported material for
different simulations with parameters different from the fiducial.

###  3.4 Dependence on different parameters

For consistency, we initiate all disks with the same gas/solid surface density
and planetesimal distribution so we can report the efficiency of each
simulation as simply the fraction of transported mass. To begin with, we note
that the amount of transported mass depends on the movement of the resonance
location (Equation 25). In Appendix C we derive a simple equation for the
position of the resonance, finding that to a good approximation one can write:

$\displaystyle\left(\frac{a_{\rm res}}{a_{\rm
J}}\right)^{3/2}=\frac{0.02M_{\rm out}\Delta^{3}+3M_{\rm
in}\Delta^{-3/2}}{M_{\rm J}},$ (26)

where $M_{\rm in}$ and $M_{\rm out}$ are the masses of the inner and outer
disks, respectively, and $\Delta$ is the gap width as defined in Eq. C8.
Looking at Equation (26), we can see the resonance position depends on the
Jovian’s mass, as well as the inner and outer disk masses. In what follows we
will explore how these parameters affect the transported mass. The inner and
outer disk mass evolution will be parameterized by the $\alpha$ viscosity
parameters. We will also briefly discuss the dependence of the transported
mass on the Jovian’s position and its orbital eccentricity.

In Fig. 5 we vary 2 of these parameters (panel A varying disk viscosities and
panel B varying the mass of the Jovian). We will discuss each panel in its
corresponding subsections. Also, note that much more material is transported
compared to the case in which disk gravity is not accounted for (in Fig. 5,
for both panels, compare the orange and black lines, with and without disk
gravity respectively, for the fiducial parameters). Next, we will explore the
dependence of the fraction of transported mass on each parameter.

Figure 5: Fraction of transported material as a function of time. The total
initial mass of planetesimals is 1.8 earth masses, equivalent to 1% of the gas
mass in the inner disk. Different coloured curves in panel A correspond to
different values of $\alpha_{\rm in}$ as highlighted therein; namely, 1$\times
10^{-4}$, 1$\times 10^{-3}$, 2.2$\times 10^{-3}$, and 4.6$\times 10^{-3}$.
Results in panel B correspond to different planetary masses, namely, 0.7, 1,
3, 5, and 7 $M_{\rm Jup}$, and are shown using purple, red, orange, green, and
blue curves, respectively. In both panels, the black curves represent the case
where the gravitational effects of the disk is not taken into account. The
colored circle symbols correspond to rtrans as defined in Eq. (24); see the
colorbar. All other planet–disk parameters are set equal to their fiducial
values (see Table 1).  Figure 6: Transported mass fraction at the end of the
10 Myr simulation for different combinations of $\alpha_{\rm in}$ and
$\alpha_{\rm out}$ in simulations which are otherwise identical to the
fiducial case (Table 1). The red lines mark the contours of constant
$\alpha_{\rm out}/\alpha_{\rm in}$; with a ratio of 1, 0.1 and 0.01 (top to
bottom). The diamond symbols correspond to the same simulations that are
depicted in Fig. 5(A) with the same color coding: i.e., with $\alpha_{\rm in}$
= 1$\times 10^{-4}$, 1$\times 10^{-3}$, 2.2$\times 10^{-3}$, and 4.6$\times
10^{-3}$ (shown in red, orange, green, and blue, respectively).

####  3.4.1 Dependence on disk viscosities

In Fig. 5A we can see that higher than fiducial $\alpha_{\rm in}$ values
result in less efficient transport as the resonance moves too quickly (see Eq.
(3.2)) compared with the fiducial case. On the other hand, for the lower than
fiducial $\alpha_{\rm in}$, we would expect the transported mass to catch up
if we let the disk evolve for the same number of viscous times as the
fiducial. In fact, our results depend on the time we let the disk evolve (10
Myrs in all simulations used for Fig. 5), so the values stated here would
actually be a lower limit. If the disk lives longer, the resonance would also
move slower, thus transporting material more efficiently.

A more exhaustive search of the parameter space for viscosities (following the
quoted values for $\alpha$ viscosity parameters listed in Table 1) only
showing the final fraction of transported mass can be found on Fig. 6 (some
simulations from Fig. 5A have been marked on this plot). The peak of
transported material in Fig. 6 depends on the Jovian mass and in the case of 3
$M_{\rm Jup}$ the simulations transporting the most material are the ones with
a ratio $\alpha_{\rm out}/\alpha_{\rm in}\sim$ 0.1 - 1. Consistently, these
models with $\alpha_{\rm out}<\alpha_{\rm in}$ thought to be more
realistic333Because of the higher irradiation of the star in the inner disk,
more ionized material couples with the magnetic field of the star promoting
the Magnetorotational Instability and hence higher turbulence which translates
into a higher viscosity (Gammie, 1996; Armitage, 2022).. As we had mentioned
previously, if $\alpha_{\rm in}$ is too high, then very little material is
transported as the resonance location moves too quickly (which also explains
the lower values of $f$ in the right part of the parameter space of Figure 6).
In other cases, similar ratios of viscosities give similar fractions of
transported material because the inner and outer masses remain proportional.

####  3.4.2 Dependence on mass of the giant planet

We now turn our attention to the dependence of transported mass on planetary
mass. To this end, we consider the evolution of planetesimals in otherwise
identical disk–planet systems but differing in the value of the planetary mass
within the range $[0.1-7]M_{\rm Jup}$ ; see Table 1. The results are
summarized in Figure 5(B). We generally find that simulations adopting the
fiducial planetary mass of 3 $M_{\rm Jup}$ transport the highest fraction of
solid material. This can be understood as follows. Increasing the planet’s
mass results in the resonance condition being established closer to the star,
while making the planet smaller makes the resonance location move closer to
the planet where there is no solid material in our simulations (see Section
5.3.3). Both of these effects reduce the transported mass by reducing the
range of solid material swept by the resonance.

Given the strong dependence of the fraction of transported mass on the
planet’s mass, we provide a more detailed picture in Figure 7(B). It is
evident that the transported mass peaks at about 2 or 3 Jupiter masses. To
explain this we need to look at the range that the resonance sweeps over.

We approximate the position of the resonance from the simulations (blue and
green solid curves for initial and final positions, respectively) using Eq.
(26) (blue and green dashed lines) and the transported material with Eq. (25)
(purple line on panel B). The result is shown as the red dashed line in panel
B and provides a rough description of the the transported mass. Three regions
can be distinguished in both panels. (i) $M_{\rm J}\lesssim 0.5M_{\rm Jup}$:
the final resonance location is outside the region with material ($0.1-1.5$
au) and an insignificant amount of material is transported; (ii) $0.5M_{\rm
Jup}\lesssim M_{\rm J}\lesssim 2M_{\rm Jup}$: the fraction increases as the
resonance sweeping region gets in until it peaks when the starting resonance
position enters the region with material, (iii) $M_{\rm J}\gtrsim 2M_{\rm
Jup}$: as the initial resonance location approaches the inner edge of the
disk, the material transported goes down again. This simple model we describe
in detail in Appendix C together with Eq. (25) gives us the same qualitative
behavior as full simulations. This includes the upward and downward slope and
the peak in the Jovian’s mass at around 3 $M_{\rm Jup}$ as can be seen in Fig.
7B.

Figure 7: The relationship between the radial range over which the resonance
sweeps and the fraction of transported mass of planetesimals. Panel A shows
the range covered by the resonance as a function of $M_{\rm J}$, starting at
the blue solid line and ending at the green solid one. The approximate initial
and final resonance positions as obtained using Eq. (26) are shown with blue
and green dashed curves, respectively. Panel B shows the amount of transported
material, both as measured from the simulations (purple curve) and as
estimated using Eq. (25) and (26) (shown using red dashed curve). See the text
(Section 3.4.2) for details.

Furthermore, we can see in Fig. 8 that this shape, including the peak at
around $\sim 1-3$ Jupiter masses, is very similar for many combinations of
inner and outer viscosities. This is appealing, especially given that the
typical measured masses of cold Jupiters in exoplanet systems fall within this
range (e.g. Fulton et al. (2021)). Finally, to explain the slight difference
in shape for the green and red curves shown in Fig. 8 at low Jupiter masses,
we note that the final position of the resonance for these cases are already
inside the region with solid material (because the outer region gets depleted
faster in these cases), so we do not see a first change in slope as with the
blue and orange curves.

Figure 8: Fraction of transported mass of planetesimals at the end of the 10
Myr simulation as a function of Jupiter masses for different combinations of
inner and outer viscosities as shown in the legend. The blue curve corresponds
to the fiducial viscosities (Table 1).

####  3.4.3 Dependence on other parameters

As mentioned before, other parameters can affect the transport of material.
The most important ones, besides the parameters mentioned in the previous
subsections, are the Jupiter’s eccentricity and position.

According to our simple model described in Appendix C, the position of the
resonance scales with the Jovian’s position (Equation 26). In this sense, we
would expect the system to just scale with the position of the Jupiter.
However, the available solid material reservoir would be different. Also, the
resonance position movement timescale will be affected. Namely, if the
resonance starts with the Jovian further away, it will travel faster through
the inner disk.

In the case of the eccentricity, from Eq. (B2), one can see that the
eccentricity peak of the resonance scales with the Jovian’s eccentricity. This
will directly impact the migration timescale as can be seen from Eq. (3.2). A
lower Jovian’s eccentricity makes the migration timescale longer, hindering
the transport. We have tested this by lowering $e_{J}$ to 0.01 and observed
little transport, and only of the smallest planetesimals. However, we expect
that gap-opening Jovians may retain some eccentricity due to disk-planet
interactions with values around 0.07 to 0.09 according to the simulations by
Duffell & Chiang (2015).

We further discuss the effect of other parameters not mentioned here such as
the disk evolution time in Section 5.

##  4 Results: Velocity Dispersion

In the previous section, we demonstrated that a significant amount of solid
material can be transported from regions close to the giant planet to the
inner regions of the disk. In order to asses whether this may result in planet
formation we need to calculate the velocity dispersion ($\sigma_{v}$), which
is a function of eccentricity and $\Delta\varpi$, for the planetesimals in
each individual ring. This will help us estimate the collisional velocity of
this planetesimals once we include collisions in our simulations.

To this end, we first introduce the planetesimal eccentricity vector as
follows:

$\mathbf{e}=(e\cos\varpi,e\sin\varpi),$ (27)

and an analog to a center of mass velocity

$\textbf{e}_{\rm CM}=\frac{\Sigma_{i}\textbf{e}_{i}m_{i}}{\Sigma_{i}m_{i}},$
(28)

where i = 1, …, N for all the planetesimals that fall into a specific ring.
Then we compare each individual eccentricity vector of the planetesimals with
this center of mass eccentricity using the following relation:

$\langle e^{2}\rangle=\frac{\Sigma_{\rm i}(\textbf{e}_{i}-\textbf{e}_{\rm
CM})^{2}m_{i}}{\Sigma_{i}m_{i}}.$ (29)

Finally, the dispersion in velocity would be proportional to the dispersion of
these eccentricity vectors, such that,

$\sigma_{v}=v_{K}\sigma_{e}=v_{K}\sqrt{\langle e^{2}\rangle}.$ (30)

In Fig. 9A we show the velocity dispersion for our fiducial case as a function
of time and semimajor axis. We observe that in regions near the resonance,
$\sigma_{v}$ increases as the eccentricity becomes larger when they approach
$e_{\rm peak}$. Additionally, there is a region of misalignment surrounding
this position (see Fig. 3c). Inner to this region, the velocities can get as
low as $\sim 10-100$ m/s (which is lower than if the effect of the disk
gravity was not included; see e.g. Fig. 10). This lower velocity dispersion is
also due to the segregation in sizes of the planetesimals (as seen in Fig. 4),
because planetesimals of similar size tend to be more aligned. It is worth
noting that even if these higher velocity dispersion near the resonance result
in destructive collisions, the smaller fragmented planetesimals would migrate
faster and accumulate in the region with a lower velocity dispersion. We also
include the surface density of solids in Fig. 9B as a function of time and
semimajor axis compared to the initial distribution of solids. Here we can see
that most of the material concentrates around the resonance, the region swept
by the resonance gets depleted and accumulated ahead.

Figure 9: Panel A shows the velocity dispersion of planetesimals computed
according to Eq. (30) as a function of time and semimajor axis, using results
of the fiducial simulation (Table 1). The red line marks the position of the
resonance. Note that the velocity dispersion increases near the resonance due
to eccentricities being higher there and also the apsidal angle alignment
varies quickly near the peak as can be seen in Fig 3c. Panel B shows the
corresponding surface density of solids compared to the initial solid
distribution. Note that most of the material is near the resonance and the
surface density interior to the resonance increases with time. Figure 10: Same
as Fig. 9 but comparing the fiducial case with the case where no disk gravity
(NDG) was considered. Both in log scale.

##  5 Discussion

We have shown that, by including the effect of disk gravity for an evolving
gaseous disk, a large fraction of the solid material can be transported. This
transporting effect seems to be very robust, as it works with a wide range of
typical Jovian masses and disk viscosities. This effect is also robust against
different choices for the planetesimal distribution.

In the following, we first discuss the relevance of including the disk gravity
in an evolving disk in Section 5.1. Then, what implications our mechanism has
on inner planet formation in Section 5.2. We mention some caveats in Section
5.3 and finally, some future work we hope to include for the next paper in
Section 5.4.

###  5.1 The importance of disk gravity in an evolving disk

One crucial ingredient in our model is the inclusion of the disk gravity. To
highlight its importance, we now compare our fiducial results with the case
where disk gravity is not included which, from hereon, we refer to as the NDG
case. In this case there is no secular resonance, and so, only the smaller
planetesimals would migrate. In Fig. 10A we compare the velocity dispersion of
planetesimals (according to Eq. (30)) in the fiducial case to that in the NDG
case. We can see the peak in eccentricity promotes higher dispersion
velocities near the resonance, but away from it, lower velocity dispersions of
planetesimals can be achieved compared to the NDG case. This is appealing,
because previous works which did not include disk gravity in their
calculations (e.g. Kortenkamp et al. (2001), Guo & Kokubo (2021)) do already
report on the planet-formation-friendly environment that the Jovian, together
with the gas drag, could produce. In this context, the inclusion of disk
gravity, which further reduces the velocity dispersion, would further enhance
the possibility of planet formation.

The effects of this disk-induced secular resonance have also been studied in
the shake-up model for the Solar System (Nagasawa et al., 2005), as well as in
the context of planet formation in and around binaries (Rafikov & Silsbee,
2015a, b; Silsbee & Rafikov, 2021). It is worth noting that these studies
concentrate more on the coupled eccentricity–apsidal angle evolution of
planetesimal orbits, but not on their transport. Additionally, these works
either do not include the evolution of the disk or use a simplistic model
(e.g., a uniformly- decaying power-law disk). Instead, our work shows that the
evolution of the disk is crucial to predict the evolution of the position of
the secular resonance and the subsequent transport of material.

###  5.2 Relevance to inner planet formation

Coming back to what we set off to explore, namely an explanation for the
correlation observed between cold giants and super-Earths. We have proven
that, thanks to the effect of planetesimal transport by the sweeping of the
secular resonance, a significant amount of material with small velocity
dispersion is accumulated in the inner region of the disk ( $\lesssim$ 0.3au).
Additionally, the transported solid mass is segregated in rings by size. All
of this occurs by the time the surface density of gas is largely depleted
($\leq 10{\rm g/cm^{2}}$, Fig. 3) which could promote the formation of planets
as seen in Morbidelli (2020). This environment would be ideal for
planetesimals to assemble into super-Earths (Lee & Chiang, 2016). Another
possibility would be to form warm or hot Jupiters, if the transported solid
material is sufficient and there is still enough gas to accrete (Batygin et
al., 2016; Boley et al., 2016). This may be the case as the population of hot
Jupiters generally have external companions (Bryan et al., 2016), similar to
the close-in Super-Earths.

It is also possible that the material transported results in a cascade of
collisions (especially for planetesimals at and around the resonance where the
eccentricites are relatively large). This could have the effect of grinding
the material down which then would migrate faster producing a ring of dust or
a recycled flux of pebbles that would affect subsequent planet formation by
forming planetesimals via streaming instability (Drążkowska et al., 2016).

We defer a detailed anaylisis of these possibilities to a future work and next
discuss the implications of our mechanism for exoplanet systems and the Solar
System.

####  5.2.1 Implications to exo-planetary systems

Once we implement collisional evolution of planetesimals, we could start
making more specific predictions about the correlation between inner planets
and outer giants. For example, the California Legacy Survey (Rosenthal et al.,
2022) found that this correlation breaks down if they consider outer giants
more massive than about 0.38 MJ and closer to the star (0.3–3 au). Similar
studies, including the Kepler Giant Planet Survey (Weiss et al., 2023) will
shed light on the various trends (if any) between the gas giant properties and
thoose of inner sub-Neptunes that will allow us to compare with our models.

Eventually, we hope to have a more complete panorama with observations from
TESS and GAIA ($\sim 100$ systems with super-Earths coexisting with outer
Jovians with absolute masses in the range $1-5$ au (Espinoza, Zhu & Petrovich
2023, in prep). Furthermore, with future observations from JWST on atmospheric
compositions, we could better determine where the material that forms the
inner planets came from (as the planetesimal composition varies with its
position relative to the icelines), providing us again with more constraints
for our model.

####  5.2.2 Implications to the Solar System

This mechanism is also relevant for the Solar System as it also has cold giant
planets. Our results indicate that, during the gas disk depletion, Jupiter
would have transported copious amounts of solid material as it falls within
the transported mass peak observed in Fig. 8 (albeit, somewhat to the lower
end of the peak). However, we have further checked that accounting for the
effects of a second giant planet (i.e. Saturn), puts a constraint on the
minimum value the precession rate of Jupiter’s orbit can have (see Fig. 2),
reducing the range of sweeping for the secular resonance, which reduces the
transported mass.

We plan on expanding our analysis by including the effects of a second giant
planet in a future paper. This kind of simulation has already been performed
by Bromley & Kenyon (2017) to explain the low mass of Mars and the asteroid
belt by focusing on the depletion of material that is swept by the secular
resonance. This work, however, makes few remarks on the material that is swept
to the inner regions. In our work we also see this effect of depletion left
behind after the secular resonance has swept. Nonetheless, we care more about
the accumulation of material at and interior to the resonance position.

There is an alternative besides the possibility of this transported solid
material forming planets. The planetesimals migrating inwards away from
Jupiter into the inner region of the disk could have plowed through earlier
formed planets in our Solar System, which could be an explanation of why we do
not have a multi-planet compact system of super-Earths as many of the other
exoplanetary systems do. This plowing of material moved by Jupiter is very
similar to the one described in the Grand Tack model (Batygin & Laughlin,
2015), but without fine tuned positions and migration histories of Jupiter and
Saturn.

###  5.3 Caveats

We now discuss some of our model assumptions and caveats and how we expect
these to impact our results.

####  5.3.1 Jovian’s growth

A more realistic simulation would have the Jovian begin as a core and undergo
runaway accretion. The timescale for the formation of the core is a few Myrs
(Pollack et al., 1996; Guilera et al., 2020). During this time, the resonance
will remain close to the planet so no sweeping will take place. On the other
hand, the runaway accretion is too quick ($\sim 10^{5}\rm yrs$) for the system
to react, so the planetesimals just align themselves with the Jovian once it
stops growing. We have checked that including the Jovian’s growth in our
simulations is equivalent to starting with a partially depleted disk.

####  5.3.2 Outside the giant planet

In our study, we focused on the inner region of the disk, ignoring the outer
region, because of the correlation seen between giant cold planets and inner
super-Earths. However, there could also be a traveling secular resonance and
an accumulation of material in the outer region that could create a second
giant planet external to the first one as explored in Guo & Kokubo (2021).
This could work together with the dust traps that are known to be created at
the edges of giant planets (Eriksson et al., 2022) to form planetesimals and
subsequently a second planet outside.

####  5.3.3 Beyond the Lagrange-Laplace approximation: MMRs and higher-order
terms

We have neglected the effect of higher-order terms in eccentricity
($\mathcal{O}[e^{4}]$), and all terms which included inclinations. In general,
the eccentricity rarely increases above $\simeq 0.3$ and it only does so at
early times when the resonance is closest to the planet (e.g., panel A1 of
Fig. 3). Furthermore, this eccentricity is already large enough to promote
migration of the larger planetesimals. Thus, the Laplace-Lagrange (LL)
approximation is often a reasonably good approximation. Adding small
inclinations (up to second order) would effectively replicate the effect from
the eccentricities as it enters in the gas drag with a similar dependence (Eq.
(11) - (12)) (Adachi et al., 1976), while remaining decoupled from the
eccentricities in the LL approximation.

We have to keep in mind, that as the eccentricity of the planetesimals
increase substantially near the secular resonance, it may be different at this
position, which would affect our condition for transport. As seen in Guo &
Kokubo (2021), planetesimals in regions close to an MMR do not align well in
eccentricity nor apsidal angle. So the other resonances could actually grind
down material closer to the Jovian, which would then align better and migrate
faster (because of its smaller size), so it could help with the transport of
material, or the material could get stuck in an MMR resonance. For the
fiducial distance chosen of 3 au, we avoid the major MMRs when considering
planetesimals interior to 1.5 au.

####  5.3.4 Other disk dispersal effects

We presented a very simple model for the evolution of the gas disk. At the end
of our simulations, we even still have a massive outer disk with $M_{\rm
out}\sim$ 5MJ, which would need to be dissipated by other means like photo-
evaporation or winds once an inner cavity has been formed. In this context,
other dissipative effects which mainly affect the inner or outer disk
separately (e.g. external irradiation), might have a significant effect on the
evolution of the position of the secular resonance. Indeed, as shown in
Appendix C the position and movement of the secular resonance depend
significantly on the masses of the inner and outer gas disks.

###  5.4 Future work

In this work, we have assumed that planetesimals do not interact with each
other (neither gravitationally nor collisionally), even within regions of the
disk with high surface densities of solids. If we want to actually form the
inner planets, it will be necessary to implement the collisions between
planetesimals and follow on the collisional outcome aided by the relative
velocity of the encounter, focusing mainly on their evolution in size as they
accumulate in the inner disk. This will help to separate the cases where we
could get super-Earth formation inner to the Jovian’s orbit.

Also, throughout this paper we considered all planetesimals to be in the same
plane along with the giant planet and the gas disk. A similar term to Eq. (11)
can be added to account for inclinations in their evolution and we would need
to generalize the velocity dispersion to include inclinations.

On another note, as was already mentioned, the first Jovian might help with
the formation of a second giant planet exterior to its orbit. In this sense,
including a second Jovian might be interesting, not least for applicability to
the Solar System.

##  6 Conclusion

In this work, we have studied the secular orbital dynamics of planetesimals
interior to a cold Jupiter (semi-major axes $\gtrsim 3$ au) affected by the
gravitational and drag forces from a viscously-evolving gaseous disk.

We find that the planetesimals are efficiently transported from $\sim 1$ au to
the inner regions $\sim 0.1-0.3$ au owing to the sweeping of secular
resonances and drag for a wide range of disk and cold Jupiter parameters. As a
result, most planetesimals interior to the cold Jupiter end up in a high
surface density ring (density boosted by $\sim 1-2$ orders of magnitude) in
these regions. In contrast, ignoring the disk gravity yields only a modest
fraction of transported mass (limited to the smaller size planetesimals). Our
main findings can be summarized as follows:

  * •

The inward motion of the planetesimals is most prominent for those which can
remain ahead of the inwardly-moving resonance. This favors smaller
planetesimals and sets a limit for the migration given by comparing the
timescales from the motion of the resonance and the migration of planetesimals
(condition given by Eq. [3.2]).

  * •

The transport efficiency peaks at cold Jupiter masses of $1-5M_{\rm Jup}$ for
a wide range of disk viscosities. Below this range, the transport efficiency
decreases as the resonance only sweeps near the planet, while very massive
Jovians start the sweeping already close to the star so almost no material is
transported.

  * •

The planetesimals move inwards in size-segregated rings with aligned
eccentricity vectors, thus lowering their, otherwise catastrophically large,
velocity dispersion. Although the velocity dispersion remains large near the
peak in eccentricity, their values interior to the resonance are lower
compared to the simulations that ignore the disk’s gravity.

These results point to the possibility of a region of the disk well suited for
planet formation. The amount of transported material is so enhanced in some
cases that it may even lead to the formation of massive cores before the disk
has dissipated, triggering runaway accretion (possible formation path for hot
and warm Jupiters). Overall, our work establishes the major role of the disk’s
gravity in the redistribution of planetesimals. Future work including the
collisional evolution of these stirred and size-segregated rings is needed to
quantify whether the formation of cores is indeed promoted.

S.B. acknowledges additional support from the National Agency for Research and
Development (ANID) / Scholarship Program / Doctorado Nacional grant 2021 -
21211921 and CASSACA grant CCJRF2105. A.A.S. acknowledges support by the
Alexander von Humboldt Foundation through a Humboldt Research Fellowship for
postdoctoral researchers. A.A.S. also thanks Jihad Touma and the America
University of Beirut, Lebanon, for hospitality during the early stages of this
work. C.P. acknowledges support from ANID Millennium Science Initiative-
ICN12_009, CATA-Basal AFB-170002, ANID BASAL project FB210003, FONDECYT
Regular grant 1210425, CASSACA grant CCJRF2105, and ANID+REC Convocatoria
Nacional subvencion a la instalacion en la Academia convocatoria 2020
PAI77200076.

##  Appendix A Disk evolution

Throughout our work, following the work of Guilera & Sándor (2017), we used a
full implicit Crank-Nicholson method to evolve the density profile, solving
Eq. (4) considering zero torques at the boundaries. As an example, we show in
Figure 11 the evolution of the surface density for our fiducial two-alpha
model (Table 1). We note the development of a peak in the density at the
location where the viscosity transition occurs (between $\alpha_{\rm
in}=10^{-3}$ to $\alpha_{\rm out}=10^{-4}$). This shape can be more pronounced
depending on the sharpness of the transition between the
viscosities.444Because no planetesimals go through the transition and the
precession rate of the giant planet depends much more on the amount of
material outside its orbit (because at its position, gas is severely
depleted), this transition is not relevant to our simulations. We have also
checked that, regardless of the initial slope in the inner disk, the slope of
the distribution tends to converge to the same value of about $\gamma=0.5$
within the first 1 Myr. This observation will become important for our
calculations in Appendix C where we estimate the precession rate of the
planetesimals.

Figure 11: Viscous evolution of the gaseous surface density for 20 Myr for our
fiducial case (see Section 3.1). Main features are the inner region depleting
much faster than the outer disk due to viscosity ($\alpha_{\rm
in}=10^{-3},\alpha_{\rm out}=10^{-4}$). We also see a transient behaviour at
the beginning of the simulation located at the boundary of the two regions
which depends on the width of the viscosity transition.

##  Appendix B Migration timescale

To derive a lower limit for the migration timescale of planetesimals, we need
to consider planetesimals at the peak of the forced eccentricity, right at the
secular resonance. To begin with, we note that planetesimal orbits exactly at
the resonance are perpendicularly aligned with the giant (Fig. 3). We thus set
$\cos\Delta\varpi=0$ in Equation (19) which accounts for perturbations due to
the gas drag and the Jovian, finding that the peak eccentricity $e_{\rm peak}$
is given by:

$\frac{de_{\rm p}}{dt}=A_{\rm pJ}e_{\rm J}-\frac{e_{\rm
p}^{2}}{\tau_{0}}\sqrt{\frac{5}{8}}=0$ (B1)

$e_{\rm peak}=\sqrt{A_{\rm pJ}e_{\rm J}\tau_{0}\sqrt{\frac{8}{5}}}.$ (B2)

Replacing Eq. (B2) into Eq. (11) gives us the migration timescale:

$\frac{1}{\tau_{\rm mig}}=\left\lvert\frac{1}{a_{\rm p}}\frac{da_{\rm
p}}{dt}\right\rvert=\frac{2}{\tau_{0}}\sqrt{\frac{5}{8}e_{\rm
peak}^{2}+\eta^{2}}\left[\eta+\left(\frac{\Gamma}{4}+\frac{5}{16}\right)e_{\rm
peak}^{2}\right]\approx\frac{2}{\tau_{0}}\sqrt{\frac{5}{8}e_{\rm
peak}^{2}}\left(\eta+\frac{1}{2}e_{\rm peak}^{2}\right).$ (B3)

As $\eta\sim\frac{1}{2}h^{2}\sim 10^{-4}$ it is safe to assume that $e_{\rm
p}^{2}$ is bigger than $\eta^{2}$, thus:

$\frac{1}{\tau_{\rm mig}}=\sqrt{\frac{A_{\rm pJ}e_{\rm J}}{\tau_{0}}}(A_{\rm
pJ}e_{\rm J}\tau_{0}+h^{2}).$ (B4)

Assuming a planetesimal density of 1 g/cm3 and $h^{2}\sim e_{\rm peak}$, we
get the following equation considering the parameters of the fiducial case:

$\displaystyle\tau_{\rm res}>\tau_{\rm
mig}\approx\frac{1}{2}\tau_{0}^{-1/2}\left(A_{\rm pJ}e_{\rm J}\right)^{-3/2},$
$\displaystyle\approx$ $\displaystyle 40{\rm
Myr}\left(\frac{M_{c}}{1M_{\sun}}\right)\left(\frac{M_{\rm J}}{3M_{\rm
Jup}}\right)^{-3/2}\left(\frac{\Sigma_{g}}{10\rm
g/cm^{2}}\right)^{1/2}\left(\frac{h}{0.1}\right)^{-1/2}\left(\frac{s}{100\rm
km}\right)^{-1/2}\left(\frac{a_{\rm J}}{3\rm au}\right)^{6}\left(\frac{a_{\rm
p}}{0.25\rm au}\right)^{-9/2}\left(\frac{e_{\rm J}}{0.05}\right)^{-3/2}.$ (B5)

where, we remind, $s$ is the radius of the planetesimal.

##  Appendix C Resonance location: approximate treatment

To understand how the resonance position $a_{\rm res}$ depends on the
parameters of our simulations, we will approximate the two-$\alpha$ disk model
by two profiles with total masses of $M_{\rm in}$ and $M_{\rm out}$,
respectively and slope of $\gamma$ for both disks. We also denote by $r_{\rm
in}$ and $r_{\rm out}$ the inner and outer edges of the gap carved by the
planet, respectively. We assume that $r_{\rm cut}$, the cutoff radius for the
disk, is sufficiently far away from the gap ($r_{\rm cut}/r_{\rm out}\gg 1$).

Let us first consider the giant planet. Since it opens a gap around its orbit,
its apsidal precession rate $\dot{\varpi}_{J}$ can be computed by accounting
for the gravity of the inner and outer disks. Since the latter is more
massive, i..e, $M_{\rm out}\gg M_{\rm in}$, it is expected to dominate
$\dot{\varpi}_{J}$. Accordingly, we can employ the disturbing function due to
the outer disk, with the latter being the dominant (simply because ).
Accordingly, we approximate $\dot{\varpi}_{J}$ using the disturbing function:

$\displaystyle\mathcal{R_{\rm out}}$ $\displaystyle=$
$\displaystyle\left[\int_{\rm r_{\rm out}}^{\rm r_{\rm cut}}\frac{1}{4}a_{\rm
J}^{2}n_{\rm J}^{2}\frac{(2-\gamma)M_{\rm out}r^{1-\gamma}}{r_{\rm
out}^{2}M_{c}}\left(\frac{a_{\rm J}}{r}\right)^{2}b^{(1)}_{\rm 3/2}(a_{\rm
J}/r)du\right]\left[\frac{1}{2}e_{\rm J}^{2}+\mbox{cst.}\right],$ (C1)

to arrive at (see e.g. Petrovich et al., 2019; Sefilian et al., 2021):

$\displaystyle\dot{\varpi}_{\rm J}$ $\displaystyle=$
$\displaystyle\frac{3n_{\rm J}}{4}\frac{M_{\rm out}}{M_{c}}\left(\frac{a_{\rm
J}}{r_{\rm out}}\right)^{3}\mathcal{B}(\gamma,r_{\rm out},r_{\rm cut}).$ (C2)

Here, $\mathcal{B}$ is a correction term that accounts for approximating the
Laplace coefficient $b_{3/2}^{(1)}(\alpha)$ as $\approx 3\alpha$; it is a weak
function of the ratio $r_{\rm cut}/r_{\rm out}$ and the surface density slope
$\gamma$. For our fiducial parameters, we find that $\mathcal{B}\approx 0.02$.

Let us now consider a planetesimal embedded in the inner disk. Its orbit will
undergo apsidal precession at a rate $\dot{\varpi}_{\rm p}$ determined by the
inner and outer disks, as well as the planet. We find that by and large, the
contribution the planet dominates over that of the outer disk, and so we
neglect the effects of the latter in the following. Accordingly, we can write

$\displaystyle\dot{\varpi}_{\rm p}$ $\displaystyle=$ $\displaystyle A_{\rm
pp}+A_{\rm pJ}\frac{e_{\rm J}}{e_{\rm p}}\cos(\varpi_{\rm p}-\varpi_{\rm
J})-\frac{M_{\rm
in}}{M_{c}}\bigg{(}1-\frac{\gamma}{2}\bigg{)}\left(\frac{r_{\rm in}}{a_{\rm
p}}\right)^{\gamma-2}n_{\rm p}\mathcal{F}(\gamma,r_{\rm in}),$ (C3)

which can be further simplified by noting that, at the resonance, the cosine
of the apsidal angles is zero (Row C of Fig. 3). Here, the factor
$\mathcal{F}$ is proportional to $\xi(1+\eta_{2})$, where $\xi$ is a
normalization factor and $\eta_{2}$ is a correction for a finite disk; see
appendix B of Tamayo et al. (2015). For our fiducial parameters, we have
$\mathcal{F}\approx 3$. Finally, taking the limit of $a_{p}/a_{\rm J}\ll 1$,
we arrive at:

$\displaystyle\dot{\varpi}_{\rm p}$ $\displaystyle=$
$\displaystyle\frac{3n_{\rm p}}{4}\frac{M_{\rm J}}{M_{c}}\left(\frac{a_{\rm
p}}{a_{\rm J}}\right)^{3}-\frac{M_{\rm
in}}{M_{c}}\bigg{(}1-\frac{\gamma}{2}\bigg{)}\left(\frac{r_{\rm in}}{a_{\rm
p}}\right)^{\gamma-2}n_{\rm p}\mathcal{F}$ (C4)

A secular resonance occurs at the radial location where the precession rates
of the planet and the planetesimals match, i.e.,
$\dot{\varpi}_{p}=\dot{\varpi}_{J}$, which yields the following relation:

$\displaystyle\frac{3}{4}\frac{M_{\rm out}}{M_{c}}\left(\frac{a_{\rm
J}}{r_{\rm out}}\right)^{3}n_{\rm J}B=\frac{3}{4}\frac{M_{\rm
J}}{M_{c}}\left(\frac{a_{\rm res}}{a_{\rm J}}\right)^{3}n_{\rm p}-\frac{M_{\rm
in}}{M_{c}}\bigg{(}1-\frac{\gamma}{2}\bigg{)}\left(\frac{r_{\rm in}}{a_{\rm
res}}\right)^{\gamma-2}n_{\rm p}F$ (C5)

Assuming that the gap around the planet is of similar width on either side of
the orbit such that $\Delta=a_{\rm J}/r_{\rm out}\approx r_{\rm in}/a_{\rm
J}$, and noting that $n_{\rm p}/n_{\rm J}=(a_{\rm res}/a_{\rm J})^{-3/2}$,
Equation (C5) can be simplified to read:

$\displaystyle M_{\rm out}\Delta^{3}\mathcal{B}=\left[M_{\rm
J}\left(\frac{a_{\rm res}}{a_{\rm
J}}\right)^{3}-(4-2\gamma)\Delta^{\gamma-2}\left(\frac{a_{\rm res}}{a_{\rm
J}}\right)^{2-\gamma}M_{\rm in}\right]\left(\frac{a_{\rm res}}{a_{\rm
J}}\right)^{-3/2},$ (C6)

In Fig. 11, we can see that the slope of the surface density becomes shallower
as time goes on. When the slope approaches $\gamma=0.5$ at late times,
Equation (C6) is easy to solve, which yields:

$\displaystyle\left(\frac{a_{\rm res}}{a_{\rm
J}}\right)^{3/2}=\frac{0.02M_{\rm out}\Delta^{3}+3M_{\rm
in}\Delta^{-3/2}}{M_{\rm J}}.$ (C7)

We note that a better approximation of the position of the resonance can be
obtained if we account for the exact surface density of the inner disk (rather
than approximating it as a simple power-law). Indeed, at around $r\approx
1-1.5$ au, the index $\gamma$ increases due to the giant planet’s gap, and in
the range $\approx 0.1-0.2$ au, $\gamma$ becomes negative as the slope inverts
(see Fig. 3). Additionally, we can use the approximation for the gap width in
Eq. (7) from Kanagawa et al. (2016) to define a more accurate value for
$\Delta$, such that:

$\displaystyle\Delta=1-\frac{\Delta_{\rm gap}}{2a_{\rm J}}\approx 0.85.$ (C8)

This completes our derivation of an approximate expression for $a_{\rm res}$.

## References

  * Adachi et al. (1976) Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756, doi: 10.1143/PTP.56.1756
  * Armitage (2022) Armitage, P. J. 2022, arXiv e-prints, arXiv:2201.07262, doi: 10.48550/arXiv.2201.07262
  * Batygin et al. (2016) Batygin, K., Bodenheimer, P. H., & Laughlin, G. P. 2016, ApJ, 829, 114, doi: 10.3847/0004-637X/829/2/114
  * Batygin & Laughlin (2015) Batygin, K., & Laughlin, G. 2015, Proceedings of the National Academy of Science, 112, 4214, doi: 10.1073/pnas.1423252112
  * Boley et al. (2016) Boley, A. C., Granados Contreras, A. P., & Gladman, B. 2016, ApJ, 817, L17, doi: 10.3847/2041-8205/817/2/L17
  * Bromley & Kenyon (2017) Bromley, B. C., & Kenyon, S. J. 2017, AJ, 153, 216, doi: 10.3847/1538-3881/aa6aaa
  * Bryan et al. (2019) Bryan, M. L., Knutson, H. A., Lee, E. J., et al. 2019, The Astronomical Journal, 157, 52, doi: 10.3847/1538-3881/aaf57f
  * Bryan et al. (2016) Bryan, M. L., Knutson, H. A., Howard, A. W., et al. 2016, ApJ, 821, 89, doi: 10.3847/0004-637X/821/2/89
  * Chiang & Laughlin (2013) Chiang, E., & Laughlin, G. 2013, Monthly Notices of the Royal Astronomical Society, 431, 3444–3455, doi: 10.1093/mnras/stt424
  * Davydenkova & Rafikov (2018) Davydenkova, I., & Rafikov, R. R. 2018, ApJ, 864, 74, doi: 10.3847/1538-4357/aad3ba
  * Drążkowska et al. (2016) Drążkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105, doi: 10.1051/0004-6361/201628983
  * Duffell (2020) Duffell, P. C. 2020, The Astrophysical Journal, 889, 16, doi: 10.3847/1538-4357/ab5b0f
  * Duffell & Chiang (2015) Duffell, P. C., & Chiang, E. 2015, The Astrophysical Journal, 812, 94, doi: 10.1088/0004-637x/812/2/94
  * Eriksson et al. (2022) Eriksson, L. E. J., Ronnet, T., Johansen, A., et al. 2022, A&A, 661, A73, doi: 10.1051/0004-6361/202142391
  * Fulton et al. (2021) Fulton, B. J., Rosenthal, L. J., Hirsch, L. A., et al. 2021, The Astrophysical Journal Supplement Series, 255, 14, doi: 10.3847/1538-4365/abfcc1
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355, doi: 10.1086/176735
  * Gong et al. (2019) Gong, M., Zheng, X., Lin, D. N. C., et al. 2019, ApJ, 883, 164, doi: 10.3847/1538-4357/ab3e70
  * Guilera & Sándor (2017) Guilera, O. M., & Sándor, Z. 2017, A&A, 604, A10, doi: 10.1051/0004-6361/201629843
  * Guilera et al. (2020) Guilera, O. M., Sándor, Z., Ronco, M. P., Venturini, J., & Miller Bertolami, M. M. 2020, A&A, 642, A140, doi: 10.1051/0004-6361/202038458
  * Guo & Kokubo (2021) Guo, K., & Kokubo, E. 2021, The Astronomical Journal, 162, 115, doi: 10.3847/1538-3881/ac0895
  * Hahn (2003) Hahn, J. M. 2003, The Astrophysical Journal, 595, 531–549, doi: 10.1086/377195
  * Hahn & Malhotra (1999) Hahn, J. M., & Malhotra, R. 1999, The Astronomical Journal, 117, 3041, doi: 10.1086/300891
  * Heppenheimer (1980) Heppenheimer, T. A. 1980, Icarus, 41, 76, doi: 10.1016/0019-1035(80)90160-8
  * Izidoro et al. (2021) Izidoro, A., Bitsch, B., & Dasgupta, R. 2021, ApJ, 915, 62, doi: 10.3847/1538-4357/abfe0b
  * Izidoro et al. (2015) Izidoro, A., Raymond, S. N., Morbidelli, A., Hersant, F., & Pierens, A. 2015, The Astrophysical Journal, 800, L22, doi: 10.1088/2041-8205/800/2/l22
  * Johansen & Lacerda (2010) Johansen, A., & Lacerda, P. 2010, Monthly Notices of the Royal Astronomical Society, doi: 10.1111/j.1365-2966.2010.16309.x
  * Kanagawa et al. (2016) Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, PASJ, 68, 43, doi: 10.1093/pasj/psw037
  * Kortenkamp et al. (2001) Kortenkamp, S. J., Wetherill, G. W., & Inaba, S. 2001, Science, 293, 1127, doi: 10.1126/science.1062391
  * Lambrechts et al. (2014) Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35, doi: 10.1051/0004-6361/201423814
  * Lee & Chiang (2016) Lee, E. J., & Chiang, E. 2016, ApJ, 817, 90, doi: 10.3847/0004-637X/817/2/90
  * Lemaitre & Dubru (1991) Lemaitre, A., & Dubru, P. 1991, Celestial Mechanics and Dynamical Astronomy, 52, 57, doi: 10.1007/BF00048587
  * Morbidelli (2020) Morbidelli, A. 2020, A&A, 638, A1, doi: 10.1051/0004-6361/202037983
  * Murray & Dermott (2000) Murray, C. D., & Dermott, S. F. 2000, Solar System Dynamics (Cambridge University Press)
  * Nagasawa et al. (2003) Nagasawa, M., Lin, D. N. C., & Ida, S. 2003, ApJ, 586, 1374, doi: 10.1086/367884
  * Nagasawa et al. (2005) Nagasawa, M., Lin, D. N. C., & Thommes, E. 2005, The Astrophysical Journal, 635, 578, doi: 10.1086/497386
  * Nagasawa et al. (2000) Nagasawa, M., Tanaka, H., & Ida, S. 2000, AJ, 119, 1480, doi: 10.1086/301246
  * Ormel & Klahr (2010) Ormel, C. W., & Klahr, H. H. 2010, Astronomy and Astrophysics, 520, A43, doi: 10.1051/0004-6361/201014903
  * Papaloizou & Larwood (2000) Papaloizou, J. C. B., & Larwood, J. D. 2000, Monthly Notices of the Royal Astronomical Society, 315, 823, doi: 10.1046/j.1365-8711.2000.03466.x
  * Petrovich et al. (2019) Petrovich, C., Wu, Y., & Ali-Dib, M. 2019, AJ, 157, 5, doi: 10.3847/1538-3881/aaeed9
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62, doi: 10.1006/icar.1996.0190
  * Pringle (1981) Pringle, J. E. 1981, ARA&A, 19, 137, doi: 10.1146/annurev.aa.19.090181.001033
  * Rafikov & Silsbee (2015a) Rafikov, R. R., & Silsbee, K. 2015a, ApJ, 798, 69, doi: 10.1088/0004-637X/798/2/69
  * Rafikov & Silsbee (2015b) —. 2015b, ApJ, 798, 70, doi: 10.1088/0004-637X/798/2/70
  * Raymond & Cossou (2014) Raymond, S. N., & Cossou, C. 2014, Monthly Notices of the Royal Astronomical Society: Letters, 440, L11–L15, doi: 10.1093/mnrasl/slu011
  * Rosenthal et al. (2022) Rosenthal, L. J., Knutson, H. A., Chachan, Y., et al. 2022, ApJS, 262, 1, doi: 10.3847/1538-4365/ac7230
  * Schlecker, M. et al. (2021) Schlecker, M., Mordasini, C., Emsenhuber, A., et al. 2021, A&A, 656, A71, doi: 10.1051/0004-6361/202038554
  * Sefilian & Rafikov (2019) Sefilian, A. A., & Rafikov, R. R. 2019, MNRAS, 489, 4176, doi: 10.1093/mnras/stz2412
  * Sefilian et al. (2021) Sefilian, A. A., Rafikov, R. R., & Wyatt, M. C. 2021, ApJ, 910, 13, doi: 10.3847/1538-4357/abda46
  * Silsbee & Rafikov (2021) Silsbee, K., & Rafikov, R. R. 2021, A&A, 652, A104, doi: 10.1051/0004-6361/202141139
  * Tamayo et al. (2015) Tamayo, D., Triaud, A. H. M. J., Menou, K., & Rein, H. 2015, ApJ, 805, 100, doi: 10.1088/0004-637X/805/2/100
  * Toliou et al. (2019) Toliou, A., Tsiganis, K., & Tsirvoulis, G. 2019, Celestial Mechanics and Dynamical Astronomy, 132, 1, doi: 10.1007/s10569-019-9942-0
  * Ward (1981a) Ward, W. R. 1981a, Icarus, 47, 234, doi: 10.1016/0019-1035(81)90169-X
  * Ward (1981b) —. 1981b, Icarus, 47, 234, doi: 10.1016/0019-1035(81)90169-X
  * Weiss et al. (2023) Weiss, L. M., Isaacson, H., Marcy, G. W., et al. 2023, arXiv e-prints, arXiv:2304.00071, doi: 10.48550/arXiv.2304.00071
  * Zheng et al. (2017a) Zheng, X., Lin, D. N. C., & Kouwenhoven, M. B. N. 2017a, ApJ, 836, 207, doi: 10.3847/1538-4357/836/2/207
  * Zheng et al. (2017b) Zheng, X., Lin, D. N. C., Kouwenhoven, M. B. N., Mao, S., & Zhang, X. 2017b, ApJ, 849, 98, doi: 10.3847/1538-4357/aa8ef3
  * Zhu et al. (2018) Zhu, W., Petrovich, C., Wu, Y., Dong, S., & Xie, J. 2018, ApJ, 860, 101, doi: 10.3847/1538-4357/aac6d5
  * Zhu & Wu (2018) Zhu, W., & Wu, Y. 2018, The Astronomical Journal, 156, 92, doi: 10.3847/1538-3881/aad22a