An Analytic Criterion for Turbulent Disruption of Planetary Resonances
Konstantin Batygin Fred C. Adams

 
  
 

^1Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125 
^2Department of Physics, University of Michigan, Ann Arbor, MI 48109 
^3Department of Astronomy, University of Michigan, Ann Arbor, MI 48109 

kbatygin@gps.caltech.edu
 








 
Mean motion commensurabilities in multi-planet systems are an expected
outcome of protoplanetary disk-driven migration, and their relative
dearth in the observational data presents an important challenge to
current models of planet formation and dynamical evolution. One natural
mechanism that can lead to the dissolution of commensurabilities is
stochastic orbital forcing, induced by turbulent density fluctuations
within the nebula. While this process is qualitatively promising, the
conditions under which mean motion resonances can be broken are not
well understood. In this work, we derive a simple analytic criterion that
elucidates the relationship among the physical parameters of the
system, and find the conditions necessary to drive planets out of resonance. Subsequently, we confirm our findings with
numerical integrations carried out in the perturbative regime, as well
as direct N-body simulations. Our calculations suggest that turbulent
resonance disruption depends most sensitively on the planet-star mass
ratio. Specifically, for a disk with properties comparable to the
early solar nebula with α=10^-2, only planet pairs with cumulative mass ratios smaller than (m_1+m_2)/M≲10^-5∼3M_⊕/M_⊙ are
susceptible to breaking resonance at semi-major axis of order
a∼0.1AU. Although turbulence can sometimes compromise resonant
pairs, an additional mechanism (such as suppression of resonance
capture probability through disk eccentricity) is required to
adequately explain the largely non-resonant orbital architectures of
extrasolar planetary systems.
 

An Analytic Criterion for Turbulent Disruption of Planetary Resonances
    Konstantin Batygin^1 & Fred C. Adams^2,3
    December 30, 2023
======================================================================





§ INTRODUCTION
 

Despite remarkable advances in the observational characterization of
extrasolar planetary systems that have occurred over the last two
decades, planet formation remains imperfectly understood. With the
advent of data from large-scale radial velocity and photometric
surveys <cit.>, the origins of a newly
identified census of close-in Super-Earths (planets with orbital
periods that span days to months, and masses between those of the
Earth and Neptune) have emerged as an issue of particular
interest. Although analogs of such short-period objects are absent
from our solar system, statistical analyses have demonstrated that
Super-Earth type planets are extremely common within the Galaxy, and
likely represent the dominant outcome of planet formation <cit.>.

An elusive, yet fundamentally important aspect of the Super-Earth
conglomeration narrative is the role played by orbital transport.  
A key question is whether these planets experience accretion
in-situ <cit.>,
or if they migrate to their close-in orbits having formed at large
orbital radii, as a consequence of disk-planet interactions
<cit.>.
Although this question remains a subject of active research, a number
recent studies <cit.> have pointed
to a finite extent of migration as an apparent requirement for
successful formation of Super-Earths. Moreover, structural models 
<cit.> show that the majority of Super-Earths have
substantial gaseous envelopes, implying that they formed in gas-rich
environments, where they could have actively exchanged angular
momentum with their surrounding nebulae.

Establishment of mean motion resonances in multi-planet systems has
long been recognized as a signpost of the planetary migration
paradigm. Specifically, the notion that slow, convergent evolution of
orbits towards one another produces planetary pairs with orbital
periods whose ratio can be expressed as a fraction of (typically
consecutive) integers, dates back more than half a century
<cit.>.
While distinct examples of resonant planetary systems exist within the
known aggregate of planets[Archetypical examples of short-period resonant systems include GJ 876 <cit.>, Kepler-36 <cit.>, Kepler-79 <cit.>, and Kepler-227 <cit.>.], the
overall orbital distribution shows little preference for mean motion
commensurabilities (Figure <ref>). Therefore, taken at face
value, the paradigm of orbital migration predicts consequences for the
dynamical architectures of Super-Earths that are in conflict with the
majority of observations <cit.>. Accordingly,
the fact that mean motion commensurabilities are neither
  common nor entirely absent in the observational census of extrasolar
  planets presents an important challenge to the present understanding
  of planet formation theory.

Prior to the detection of thousands of planetary candidates by the
Kepler spacecraft, the expectations of largely resonant
architectures of close-in planets were firmly established by global
hydrodynamic, as well as N-body simulations
<cit.>. An important
distinction was drawn by the work of <cit.>, who pointed out that
resonances can be destabilized by random density fluctuations produced
by turbulence within the protoplanetary disk. Follow-up studies
demonstrated that a rich variety of outcomes can be attained as a
consequence of stochastic forcing within the disk
<cit.>, and that in specific cases, turbulence
can be conducive to the reproduction of dynamical architecture
<cit.>.





While the prediction of the infrequency of resonant systems made by
<cit.> was confirmed by the Kepler dataset, recent work
has shown that turbulent forcing is not the only mechanism through
which resonances can be disrupted. Specifically, the work of
<cit.> proposed that a particular relationship between
the rates of eccentricity damping and semi-major axis decay can render
resonances metastable, while <cit.> showed that
probability of resonance capture can be dramatically reduced in
slightly non-axisymmetric disks. In light of the ambiguity associated
with a multitude of theoretical models that seemingly accomplish the
same thing, it is of great interest to inquire which, if any, of the
proposed mechanisms plays the leading role in sculpting the
predominantly non-resonant architectures of known exoplanetary
systems.

Within the context of the aforementioned models of resonant
metastability and capture suppression, the necessary conditions for
passage through commensurability are relatively clear. Resonant
metastability requires the outer planet to be much more massive than
the inner planet <cit.>, while the capture suppression
mechanism requires disk eccentricities on the order of a few percent
to operate <cit.>. In contrast, the complex interplay
between planet-planet interactions, turbulent forcing, and dissipative
migration remains poorly quantified, making the turbulent disruption
mechanism difficult to decisively confirm or refute (see e.g., ). As a result, a key goal of this work is to
identify the regime of parameter space for which the stochastic
dissolution of mean motion resonances can successfully operate. In
doing so, we aim to gain insight into the evolutionary stages of young
planetary systems during which disk turbulence can prevent the
formation of resonant pairs of planets.

The paper is organized as follows. In Section <ref>, we present
the details of our model. In Section <ref>, we employ methods
from stochastic calculus to derive an analytic criterion for
turbulent disruption of mean motion resonances. In Section
<ref>, we confirm our results with both perturbative numerical
integrations and an ensemble of full N-body simulations. The paper
concludes in Section <ref> with a summary of our results and a
discussion of their implications.



§ ANALYTIC MODEL
 

The model we aim to construct effectively comprises three ingredients:
(1) first-order (k:k-1) resonant planet-planet interactions, (2) orbital migration and
damping, as well as (3) stochastic turbulent forcing. In this section,
we outline our treatment of each of these processes. A cartoon
depicting the geometric setup of the problem is shown in Figure
<ref>. Throughout much of the manuscript, we make the
so-called “compact" approximation, where we assume that the
semi-major axis ratio ξ≡ a_1/a_2→ 1. While
formally limiting, the agreement between results produced under this
approximation and those obtained within N-body integrations is
well-known to be satisfactory, particularly for k⩾3 (see, 
e.g., ), where the integer k 
specifies the resonance <cit.>. 

Being made up of analytic components, the model constructed here
cannot possibly capture all of the intricate details of the dynamical
evolution that planets are subjected to, within protoplanetary disks.
By sacrificing precision on a detailed level, however, we hope to
construct an approximate description of the relevant physical
processes that will illuminate underlying relationships. These
findings can then be used to constrain the overall regime over which
turbulent fluctuations can effect the dynamical evolution of nascent
planetary systems.



 §.§ Planet-Planet Interactions


In the late twentieth century, it was recognized that a perturbative
Hamiltonian that represents the motion of a massive pair of planets
residing on eccentric orbits, in the vicinity of a mean-motion
commensurability, can be cast into integrable form
<cit.>. More recently,
this formalism has been used to provide a geometric representation of
resonant dynamics <cit.>, study the onset of chaos
<cit.>, generalize the theory of resonant capture
<cit.>, as well as to elucidate overstable librations
<cit.>. A key advantage of this treatment is that it translates the full, unrestricted three-body problem into the same mathematical form as that employed for the well-studied circular restricted problem <cit.>. Here, we make use of this framework once again. Because detailed derivations of the aforementioned resonant normal
form are spelled out in the papers quoted above, we will not reproduce
it here, and instead restrict ourselves to employing the results. 

The Hamiltonian that describes planet-planet interactions in the
vicinity of a k:k-1 mean motion resonance can be written as follows: 

    ℋ = 3 (ε+1) ( x^2+y^2/2) - 
    ( x^2+y^2/2)^2 - 2 x,

where the variables (ε,x,y) are defined below. 
A Hamiltonian of this form is typically referred to as the second
fundamental model for resonance
<cit.>, and behaves as a
forced harmonic oscillator at negative values of the proximity
parameter, ε, while possessing a pendulum-like phase-space
structure at large positive values of ε. This integrable
model approximates the real N-body dynamics at low eccentricities and
inclinations, and formally assumes that the orbits do not cross, 
although this latter assumption is routinely violated without much
practical consequence (see, e.g.,
).  In the well-studied case
of the restricted circular three-body problem, the canonical variables
(x,y) are connected to the test particle's eccentricity and the
resonant angle, while ε is a measure of how close the
orbits are to exact resonance. Within the framework of the full
planetary resonance problem (where neither mass nor eccentricity of
either secondary body is assumed to be null), the variables take on
slightly more complex physical meanings.

In order to convert between Keplerian orbital elements and the
dimensionless canonical variables used here, we first define a
generalized composite eccentricity  

    σ= √(e_1^2 + e_2^2- 2 e_1 e_2 cos(Δϖ)),

where subscripts 1 and 2 refer to the inner and outer planets
respectively, e is eccentricity, and ϖ is the longitude of
periastron. Additionally, we define units of action and time 
according to 

    [A]= 1/2( 15/4k M/m_1+m_2)^2/3
       [T] = 1/n(5/√(6) k^2M/m_1+m_2)^2/3,

where m is planetary mass, M is stellar mass, and n=√( M/a^3) 
is the mean motion. Then, in the compact limit, the variables in the 
Hamiltonian become <cit.>:  

    x    = σ ( 15/4k M/m_1+m_2)^1/3cos(kλ_2-(k-1)λ_1-ω̃) 
      
    y    = σ ( 15/4k M/m_1+m_2)^1/3sin(kλ_2-(k-1)λ_1-ω̃) 
    ε   = 1/3( 15/4k M/m_1+m_2)^2/3(σ^2 - Δξ/k),

where ξ=a_1/a_2, and the quantity 

    ω̃≡arctan[
    e_2sin ϖ_2-e_1sin ϖ_1  e_1cos ϖ_1-e_2cos ϖ_2]
 
represents a generalized longitude of perihelion.





The specification of resonant dynamics is now complete. While
application of Hamilton's equations to equation (<ref>) only
yields the evolution of σ and the corresponding resonant angle,
the behavior of the individual eccentricities and apsidal lines can be
obtained from the conserved[When the system is subjected to
  slow evolution of the proximity parameter ε, ρ is
  no longer a strictly conserved quantity. Instead, ρ becomes an
  adiabatic invariant that is nearly constant, except when the system
  encounters a homoclinic curve <cit.>.]
quantity ρ=m_1 e_1^2+m_2 e_2^2+m_1 m_2 e_1 e_2 cos(Δϖ). In addition, we note that the definitions of the variables 
(<ref>) are independent of the individual planetary masses
m_1, m_2, and depend only on the cumulative planet-star mass ratio
(m_1+m_2)/M. This apparent simplification is a consequence of 
taking the limit ξ≡ a_1/a_2→ 1, and is qualitatively
equivalent to the Öpik approximation <cit.>.



 §.§ Planet-Disk Interactions


Dating back to early results on ring-satellite interactions
<cit.>, it has been evident that planets can
exchange orbital energy and angular momentum with their natal
disks. For planets that are not sufficiently massive to open gaps
within their nebulae, this exchange occurs through local excitation 
of spiral density waves (i.e., the so-called “type-I" regime), and
proceeds on the characteristic timescale: 

    τ_wave= 1/n(M/m) 
    (M/Σ a^2) (h/r)^4,

where Σ is the local surface density, and h/r is the aspect
ratio of the disk. For an isothermal equation of state and a surface
density profile that scales inversely with the orbital radius 
<cit.>, the corresponding rates of eccentricity and semi-major
axis decay are given by <cit.>: 

    1/ad a/d t≡ - 1/τ_mig≃ - 4 f/τ_wave(h/r)^2
       1/ed e/d t≡ - 1/τ_dmp≃ - 
    3/41/τ_wave.

A different, routinely employed approach to modeling disk-driven
semi-major axis evolution is to assume that it occurs on a timescale
that exceeds the eccentricity decay time by a numerical factor
𝒦. To this end, we note that the value of 𝒦∼10^2 
adopted by many previous authors <cit.> is in rough 
agreement with equation (<ref>) which yields 
𝒦∼(h/r)^-2.

While eccentricity damping observed in numerical simulations 
(e.g., ) is well matched by equation
(<ref>), state-of-the-art disk models show that both the rate
and direction of semi-major axis evolution can be significantly
affected by entropy gradients within the nebula <cit.>.
Although such corrections alter the migration histories on a detailed
level, convergent migration followed by resonant locking remains an
expected result in laminar disks <cit.>.  For
simplicity, in this work, we account for this complication by
introducing an adjustable parameter f into equation (<ref>).

In addition to acting as sources of dissipation, protoplanetary disks
can also drive stochastic evolution. In particular, density
fluctuations within a turbulent disk generate a random gravitational
field, which in turn perturbs the embedded planets <cit.>.  Such
perturbations translate to effectively diffusive evolution of the
eccentricity and semi-major axis <cit.>.  In the ideal
limit of MRI-driven turbulence, the corresponding eccentricity and
semi-major axis diffusion coefficients can be constructed from
analytic arguments (e.g., see ) 
to obtain the expressions 

    _ξ = _a/a^2∼ 2 _e∼α/2( Σ a^2/M)^2 n,

where α is the Shakura-Sunayev viscosity parameter 
<cit.>. Although non-ideal effects can modify the
above expressions on the quantitative level <cit.>, for the
purposes of our simple model we neglect these explicit corrections. We
note, however, that such details can be trivially incorporated into
the final answer by adjusting the value of α accordingly.



§ CRITERION FOR RESONANCE DISRUPTION
 

With all components of the model specified, we now evaluate the
stability of mean motion resonances against stochastic
perturbations. In order to obtain a rough estimate of the interplay
between turbulent forcing, orbital damping, and resonant coupling, we
can evaluate the diffusive progress in semi-major axis and
eccentricity against the width of the resonance. Specifically, the
quantities, whose properties we wish to examine are χ≡
n_2/n_1 - k/(k-1) and x. Keep in mind that this latter quantity 
is directly proportional to the generalized eccentricity σ 
(see equation [<ref>]).



 §.§ Diffusion of Semi-major Axes
 

In the compact limit a_1≈ a_2, the time evolution of the 
parameter χ can be written in the approximate form 

    dχ/dt≃3/2 (d a_1/dt - d a_2/dt),

where  is a representative average semi-major axis. For the
purposes of our simple model, we treat a_1 and a_2 as uncorrelated
Gaussian random variables with diffusion coefficients _a; we note however, that in reality significant correlations may exist
between these quantities and such correlations could potentially alter the nature of the random walk <cit.>. Additionally, for comparable-mass planets, we
may adopt  as a characteristic drift rate, replacing m with
m_1+m_2 in equation (<ref>). Note that this assumption leads
to the maximum possible rate of orbital convergence.

With these constituents, we obtain a stochastic differential equation of the form

    dχ = 3/2√(2 _ξ)  dw - 3/2χ/ dt,

where w represents a Weiner process (i.e., a continuous-time random
walk; ). The variable χ will thus take on a distribution of values
as its evolution proceeds.  Adopting the t→∞
standard deviation of the resulting distribution function as a
characteristic measure of progress in χ, we have: 

    δχ = √(3 _ξ /2) = 
    1/4h/r√(3 α Σ ^2/f (m_1+m_2)).


The approximate extent of stochastic evolution that the system can
experience and still remain in resonance is given by the resonant
bandwidth, Δχ. At its inception[A resonance can
  only be formally defined when a homoclinic curve (i.e., a separatrix)
  exists in phase-space. For a Hamiltonian of the form (<ref>),
  a separatrix appears at ε=0, along with an unstable
  (hyperbolic) fixed point, that bifurcates into two fixed points (one
  stable and one unstable) at ε>0.], the width of the
resonance <cit.> is given by 

    Δχ≃ 5 [ √(k) (m_1+m_2)/M]^2/3.

Accordingly, a rough criterion for turbulent disruption of the resonance is
[box=]align
δχ/Δχ   ∼1/20h/r
M/m_1+m_2 √(3 α/f)

   ×[   Σ ^2/k M 
√(Σ ^2/m_1+m_2)   ]^1/3 ≳1.


Keep in mind that δχ is a measure of the width of the 
distribution in the variable χ due to stochastic evolution, 
whereas Δχ is the change in χ necessary to 
compromise the resonance. 

The above expression for resonance disruption depends sensitively on
the planet-star mass ratio. This relationship is illustrated in Figure
<ref>, where the expression (<ref>) is
shown as a function of the quantity (m_1+m_2)/M, assuming system
properties α=10^-2, h/r=0.05, ⟨ a⟩ = 0.1AU,
f=1, and k=3. The disk profile is taken to have the form
Σ=Σ_0 (r_0/r), with Σ_0=1700 g/cm^2 and
r_0=1AU, such that the local surface density at ⟨ a⟩
is ⟨Σ⟩ = 17,000 g/cm^2. Notice that the 
disruption criterion (<ref>) also depends on (the square
root of) the surface density of the disk.  A family of curves
corresponding to lower values of the surface density (i.e.,
0.1,0.2,…0.9,1 × ⟨Σ⟩) are also shown,
and color-coded accordingly.


While Figure <ref> effectively assumes a maximal rate of orbital convergence, we reiterate that hydrodynamical simulations suggest that both the speed and sense of type-I migration can have a wide range of possible values <cit.>. To this end, we note that setting f=0 in equation (<ref>) yields ∞>1, meaning that in the case of no net migration, an arbitrarily small turbulent viscosity is sufficient to eventually bring the resonant angles into circulation. Furthermore, a negative value of f, which corresponds to divergent migration, renders our criterion meaningless, since resonance capture cannot occur in this instance <cit.>.






 §.§ Diffusion of Eccentricities
 

An essentially identical calculation can carried out for stochastic
evolution of x (or y). To accomplish this, we assume that the
generalized eccentricity σ diffuses with the coefficient
√(2) _e. Accounting for conversion factors between
conventional quantities and the dimensionless coordinates (given by
equation [<ref>]), we obtain 

    _x≃α (Σ ^2/M)^2 
    (M/√(k) (m_1+m_2))^4/3.

Similarly, the damping timescale takes the form 

    τ_x≃ (128/225k  M/m_1+m_2)^1/3k M/Σ ^2(h/r)^4,

where, as before, we adopted the total planetary mass as an an
approximation for m in the expression (<ref>).

In direct analogy with equation (<ref>), we obtain 
the stochastic equation for the time evolution of x, 

    dx = √(2 _x)  dw - x/τ_x dt,

so that the distribution of x is characterized by the standard 
deviation δ x = √(_x τ_x). At the same time, we
take the half-width of the resonant separatrix to be given by 
Δ x=2 (e.g., ). Combining 
these two results, we obtain a second criterion for resonance 
disruption, i.e., 

    δ x/Δ x   ∼(h/r)^2  (√(2) k/15)^1/3√(α Σ ^2/M)
       ×(M/m_1+m_2)^5/6≳ 1  .




 §.§ Semi-major Axis vs Eccentricity
 

In order to construct the simplest possible model that still captures
the dynamical evolution adequately, it is of interest to evaluate the
relative importance of stochastic evolution in the degrees of freedom
related to the semi-major axis and eccentricity. Expressions
(<ref>) and (<ref>) both represent conditions
under which resonant dynamics of a planetary pair will be short-lived,
even if capture occurs. To gauge which of the two criteria is more
stringent, we can examine the ratio 

    δ x/Δ x/δχ/Δχ∼5√(f) ( h/r)  [ k^2 (m_1+ m_2)/M]^1/3≪ 1.

The fact that this expression evaluates to a number substantially
smaller than unity means that diffusion in semi-major axes (equation
[<ref>]) dominates over diffusion in eccentricities
(equation [<ref>]) as a mechanism for disruption of
mean-motion commensurabilities. Although the relative importance of
_a compared to _e is not obvious a priori, it likely
stems in large part from the fact that the orbital convergence
timescale generally exceeds the eccentricity damping timescales by a
large margin.






§ NUMERICAL INTEGRATIONS
 

In order to derive a purely analytic criterion for turbulent
disruption of mean motion resonances, we were forced to make a series
of crude approximations in the previous section. To assess the
validity of these approximations, in this section we test the
criterion (<ref>) through numerical integrations.  We
first present a perturbative approach (Section <ref>) and then
carry out a series of full N-body simulations (Section <ref>).



 §.§ Perturbation Theory

 

The dynamical system considered here is described by three equations
of motion, corresponding to the variations in x, and y, and
ε. Although the resonant dynamics itself is governed by
Hamiltonian (<ref>), to account for the stochastic and
dissipative evolution, we must augment Hamilton's equations with terms
that describe disk-driven evolution. As before, we adopt
τ_dmp as the decay timescale for the generalized
eccentricity, σ, and take τ_mig as the
characteristic orbital convergence time. The full equations of motion
are then given by:

    dx/dt = -3 y (1+ε)+y (x^2+y^2) - x/τ_dmp/[T]
       dy/dt = -2+3 x (1+ε) - x (x^2 + y^2 ) - y/τ_dmp/[T]
       dε/dt = 2/3([A]/k τ_mig/[T] - x^2+y^2/τ_dmp/[T])+ℱ.


In the above expression, ℱ represents a source of
stochastic perturbations. For computational convenience, we
implemented this noise term as a continuous sequence of analytic
pulses, which had the form 2 ζsin(π t/Δ t )/Δ t, 
where ζ is a Gaussian random variable. The pulse time interval
was taken to be Δ t=0.1, and the standard deviation of ζ
was chosen such that the resulting diffusion coefficient
𝒟_ζ = σ_ζ^2/Δ t matched that given
by equation (<ref>).




Note that here, we have opted to only implement stochastic
perturbations into the equation that governs the variation of
ε. Qualitatively, this is equivalent to only retaining
semi-major axis diffusion and neglecting eccentricity diffusion. To
this end, we have confirmed that including (appropriately scaled)
turbulent diffusion into equations of motion for x and y does not
alter the dynamical evolution in a meaningful way, in agreement with
the discussion surrounding equation (<ref>).

Turbulent fluctuations aside, the equation of motion for the parameter
ε indicates that there exists an equilibrium value of the
generalized eccentricity
σ_eq=√(τ_dmp/(2 k τ_mig)) that
corresponds to stable capture into resonance. Analogously, the
equilibrium value of x_eq=σ_eq√(2[A])
parallels the strictly real fixed point of Hamiltonian (<ref>).
As a result, if we neglect the small dissipative contributions and set
dx/dt=0,dy/dt=0,x=x_eq,y=0 in the first and second equations
in expression (<ref>), we find an equilibrium value of the
proximity parameter, ε_eq, that coincides with
resonant locking. An ensuing crucial point is that if resonance is
broken, the system will attain values of ε substantially
above the equilibrium value ε_eq.

In order to maintain a close relationship with the results presented
in the preceding section, we retained the same physical parameters for
the simulations as those depicted in Figure <ref>. In
particular, we adopted α=10^-2, h/r=0.05, ⟨ a
⟩=0.1AU, f=1, and k=3. Additionally, we again chose a
surface density profile with Σ_0=1700g/cm^2 at
r_0=1AU, that scales inversely with the orbital radius, such that
the nominal surface density at r=⟨ a ⟩ is ⟨Σ⟩=17,000g/cm^2. We also performed a series of
simulations that span a lower range of surface densities
(0.1,0.2,…,0.9,1 × ⟨Σ⟩). 
All of the integrations were carried out over a time span of
τ_mig=100 τ_wave, with the system initialized
at zero eccentricity (x_0=y_0=0), on orbits exterior to exact
commensurability (ϵ_0=-1).

We computed three sets of evolutionary sequences, corresponding to
planet-star mass ratios (m_1+m_2)/M=10^-6,10^-5, and 10^-4.
As can be deduced from Figure <ref>, the qualitative
expectations for the outcomes of these simulations (as dictated by
equation [<ref>]) are unequivocally clear. Resonances
should be long-term stable for (m_1+m_2)/M=10^-4 and long-term
unstable for (m_1+m_2)/M=10^-6. Meanwhile, temporary resonance
locking, followed by turbulent disruption of the commensurability
should occur for (m_1+m_2)/M=10^-5.

Figure <ref> depicts numerically computed evolution of
ε for the full range of local surface densities under
consideration (color-coded in the same way as in Figure <ref>)
as a function of time. These numerical results are in excellent
agreement with our theoretical expectations from Section <ref>.
The proximity parameter always approaches its expected equilibrium
value ε_eq for large mass ratios
(m_1+m_2)/M=10^-4 (top panel), but never experiences long-term
capture for small mass ratios (m_1+m_2)/M=10^-6 (bottom
panel). Resonance locking does occur for the intermediate case
(m_1+m_2)/M=10^-5 (middle panel). However, two evolutionary
sequences corresponding to Σ=0.7 ⟨Σ⟩ and
Σ=0.9 ⟨Σ⟩ show the system breaking out of
resonance within a single orbital convergence time, τ_mig.
It is sensible to assume that other evolutionary sequences within this
set would also break away from resonance if integrations were
extended over a longer time period.

Figure <ref> shows the phase-space counterpart of the
evolution depicted in Figure <ref>. Specifically, the x-y
projections of the system dynamics are shown for cases with surface
densities Σ=0.3 ⟨Σ⟩ (blue) and
Σ=⟨Σ⟩ (black), where the background depicts
the topology of the Hamiltonian (<ref>). In each panel, the
black curve designates the separatrix of ℋ, given the
equilibrium value of the proximity parameter
ε=ε_eq. The background color scale is a
measure of the value of ℋ. The thee equilibrium points of
the Hamiltonian are also shown, as transparent green dots.

As in Figure <ref>, three representative ratios of
planet mass to stellar mass are shown.  In the right panel (for mass ratio
  (m_1+m_2)/M=10^-4), turbulent diffusion plays an essentially
negligible role and the system approaches a null libration amplitude 
under the effect of dissipation. In the middle panel (for mass ratio
  (m_1+m_2)/M=10^-5), resonant capture is shown, but the libration
amplitude attained by the orbit is large, particularly in the case of
Σ=⟨Σ⟩. In the left panel (for mass ratio
  (m_1+m_2)/M=10^-6), the trajectory is initially advected to high
values of the action, but inevitably breaks out of resonance and
decays towards the fixed point at the center of the internal
circulation region of the portrait.



 §.§ N-body Simulations

 

In order to fully evaluate the approximations inherent to the
perturbative treatment of the dynamics employed thus far, and to
provide a conclusive test of the analytic criterion
(<ref>), we have carried out a series of direct N-body
simulations. The integrations utilized a Burlisch-Stoer
integration scheme (e.g., ) that included the full
set of 18 phase space variables for the 3-body problem consisting of
two migrating planets orbiting a central star.  For the sake of
definiteness, the physical setup of the numerical experiments was
chosen to closely mirror the systems used in the above discussion.
Specifically, two equal-mass planets were placed on initially circular
orbits slightly outside of the 2:1 mean motion resonance, so that the
initial period ratio was 0.45. The planets were then allowed to
evolve under the influence of mutual gravity, as well as disk-driven
convergent migration, orbital damping, and turbulent perturbations.




Following <cit.>, we incorporated the
orbital decay and eccentricity damping using accelerations 
of the form: 

    dv⃗/dt = -v⃗/τ_mig - 
    2 r⃗/τ_dmp(v⃗·r⃗)/(r⃗·r⃗),

where v⃗ and r⃗ denote the orbital velocity and radius
respectively.[Note that we have neglected disk-induced damping 
of the orbital inclination, because of the planar setup of the problem.]  
While both planets were subjected to eccentricity
damping, inward (convergent) migration was only experienced by the
outer planet. Simultaneously, for computational convenience, the
semi-major axis of the inner planet was re-normalized to a_1=0.1AU
at every time step[Qualitatively, this procedure is equivalent
to changing the unit of time at every time step <cit.>.].  
The characteristic timescales τ_mig and τ_dmp
were kept constant, given by equation (<ref>), adopting
identical physical parameters of the disk to those employed above. 
Finally, following previous treatments <cit.>, 
turbulent fluctuations were introduced into the equations of motion
through random velocity kicks, whose amplitude was tuned such that the
properties of the diffusive evolution of an undamped isolated orbit
matched the coefficients from equation (<ref>). For completeness, we have also included the leading order corrections due to general relativity <cit.>.

As in the previous sub-section, we computed the orbital evolution of
three representative cases with mass ratios
(m_1+m_2)/M=10^-4,10^-5, and 10^-6 (corresponding to migration timescales of τ_mig≃1.5×10^3,10^4, and 10^5  years respectively) over a time span of
0.1Myr. The numerical results are shown in Figure <ref>,
and show excellent agreement with the analytic criterion from equation
(<ref>). In particular, the system with mass ratio 
(m_1+m_2)/M=10^-4 exhibits long-term stable capture into a 3:2
mean motion resonance, as exemplified by the ensuing low-amplitude
libration of the resonant angles ϕ[3:2] = 3λ_2 - 2
λ_1-ϖ_1 and ψ[3:2] = 3λ_2 - 2 λ_1-ϖ_2,
shown in red and blue in the bottom panel of Figure <ref>. 
Correspondingly, both the period ratio (top panel) and the
eccentricities (middle panel) rapidly attain their resonant
equilibrium values, and remain essentially constant throughout the
simulation.

The case with mass ratio (m_1+m_2)/M=10^-5, for which equation
(<ref>) yields δχ / Δχ∼1,
perfectly exemplifies the transitory regime. As shown in the top
panel of Figure <ref>, where this experiment is represented
in gray, the system exhibits temporary capture into the 3:2 as well
as the 4:3 commensurabilities, and subsequently locks into a
meta-stable 7:6 resonance at time ∼15,000 years. Although evolution
within this resonance is relatively long-lived, the bottom panel of
Figure <ref> shows that the corresponding resonant angles
ϕ[7:6] = 7λ_2 - 6λ_1-ϖ_1 (green) and
ψ[7:6] = 3λ_2 - 2 λ_1-ϖ_2 (gray) maintain 
large amplitudes of libration, due to the nearly perfect balance
between orbital damping and turbulent excitation. As a result, the
system eventually breaks out of its resonant state. After a period of
chaotic scattering, the orbits switch their order, and the period
ratio increases.

Finally, the case with mass ratio (m_1+m_2)/M=10^-6 represents a
system that never experiences resonant locking. As the period ratio
evolves towards unity (purple curve in the top panel), encounters with
mean motion commensurabilities only manifest themselves as impulsive
excitations of the orbital eccentricities (purple/orange curves in the
middle panel) of the planets. As such, the planets eventually
experience a brief phase of close encounters, and subsequently
re-enter an essentially decoupled regime, after the orbits reverse. 

We note that because turbulence introduces a fundamentally stochastic
component into the equations of motion, each realization of the N-body
simulations is quantitatively unique. However, having carried out tens
of integrations for each set of parameters considered in Figure
<ref>, we have confirmed that the presented solutions are
indeed representative of the evolutionary outcomes. As a result, we 
conclude that the analytic expression (<ref>) represents
an adequate description of the requirement for resonance disruption, 
consistent with the numerical experiments.



§ CONCLUSION
 

While resonant locking is an expected outcome of migration theory <cit.>, the current sample of exoplanets shows only a mild tendency for systems to be
near mean motion commensurabilities <cit.>. Motivated by
this observational finding, this paper derives an analytic criterion
for turbulent disruption of planetary resonances and demonstrates its
viability through numerical integrations. Our specific results are
outlined below (Section <ref>), followed by a conceptual
interpretation of the calculations (Section <ref>), and finally a
discussion of the implications (Section <ref>).



 §.§ Summary of Results
 
 

The main result of this paper is the derivation of the constraint
necessary for turbulent fluctuations to compromise mean motion
resonance (given by equation [<ref>]).  This criterion
exhibits a strong dependence on the ratio of planetary mass to stellar
mass, but also has significant dependence on the local surface
density. That is, turbulence can successfully disrupt mean motion
resonances only for systems with sufficiently small mass ratios and/or large
surface densities (see Figure <ref>).

The analytic estimate (<ref>) for the conditions required
for turbulence to remove planet pairs from resonance has been verified by
numerical integrations. To this end, we have constructed a model of disk-driven resonant dynamics
in the perturbative regime, and have calculated the time evolution of
the resonance promiximity parameter ε (Section
<ref>).  The results confirm the analytical prediction that given nominal disk parameters, systems with mass ratios
smaller than (m_1+m_2)/M∼10^-5∼ 3M_⊕/M_⊙ are forced out of resonance by turbulence, whereas
systems with larger mass ratios survive (Figure <ref>). We have
also performed full N-body simulations of the problem (Section
<ref>). These calculations further indicate that planetary systems with
small mass ratios are readily moved out of resonance by turbulent
fluctuations, whereas systems with larger mass ratios are not (Figure
<ref>). Accordingly, the purely analytic treatment, simulations performed within the framework of perturbation
theory, and the full N-body experiments all yield consistent results. 

For circumstellar disks with properties comparable to the minimum mass
solar nebula <cit.>, the results of this paper suggest that compact Kepler-type planetary systems are relatively close to the border-line for stochastic disruption
of primordial mean motion commensurabilities. Nonetheless, with a cumulative mass ratio that typically lies in the range of (m_1+m_2)/M∼10^-5 - 10^-4 (Figure <ref>), the majority of these planets are
sufficiently massive that their resonances can survive in the face of
turbulent disruption, provided that the perturbations operate at the expected amplitudes
(this result also assumes that the stochastic fluctuations act over a
time scale that is comparable to the migration time).

Given critical combinations of parameters (for which equation [<ref>] evaluates to a value of order unity), resonant systems can ensue, but
they routinely come out of the disk evolution phase with large libration
amplitudes. This effect has already been pointed out in previous work
<cit.>, which focused primarily on numerical
simulations with limited analytical characterization. Importantly, this notion suggests that the stochastic forcing mechanism may be critical to setting up
extrasolar planetary systems like GJ 876 and Kepler-36 that exhibit
rapid dynamical chaos <cit.>.

Although this work has mainly focused on the evolution of sub-Jovian planets, we can reasonably speculate that turbulent fluctuations are unlikely to strongly affect mean motion
resonances among giant planets. In addition to having mass ratios well
above the critical limit, the influence that the disk exerts on large
planets is further diminished because of gap-opening <cit.>.  However, one
complication regarding this issue is that the damping rate of
eccentricity is also reduced due to the gap (e.g., ). Since both the excitation
and damping mechanisms are less effective in the gap-opening regime,
a minority of systems could in principle allow for excitation to dominate.  



 §.§ Conceptual Considerations
 
 










The analysis presented herein yields a practical measure that informs the outcome of dynamical evolution of multi-planetary systems embedded in turbulent protoplanetary disks. While numerical experiments confirm that the analytic theory indeed provides an acceptable representation of perturbed N-body dynamics, the phenomenological richness inherent to the problem calls for an additional, essentially qualitative account of the results. This is the purpose of the following discussion.

Within the framework of our most realistic description of the relevant physics (i.e., the N-body treatment), the effect of turbulent fluctuations is
to provide impulsive changes to the planet velocities. The turbulence
has a coherence time of order one orbital period, so that the
fluctuations provide a new realization of the random gravitational
field on this time scale <cit.>. With these impulses, the orbital elements of
the planets, specifically the semi-major axis a and eccentricity
e, execute a random walk. In other words, as the elements vary, the changes in a and e accumulate
in a diffusive manner <cit.>. Simultaneously, the interactions between planets and the spiral density waves they induce in the nebula lead to smooth changes in the orbital periods, as well as damping of the planetary eccentricities <cit.>. 





In contrast with aforementioned disk-driven effects, the bandwidth of a planetary resonance is typically described in terms of maximal libration amplitude of a critical angle ϕ that obeys d'Alembert rules (e.g., see Chapter 8 of ). Thus, the conceptual difficulty lies in connecting how the extrinsic forcing of orbital elements translates to the evolution of this angle. Within the framework of our theoretical model, this link is enabled by the Hamiltonian model of mean motion resonance (equation [<ref>]; ).

In the parameter range relevant to the problem at hand, the behavior of Hamiltonian (<ref>) is well-approximated by that of a simple pendulum <cit.>. Specifically, the equilibrium value of ε dictates the value of the pendulum's action, Φ, at which zero-amplitude libration of the angle ϕ can occur, as well as the location of the separatrix. Correspondingly, oscillation of the angle ϕ translates to variations of the action Φ, which is in turn connected to the eccentricities (equations [<ref>]) as well as the semi-major axes, through conservation of the generalized Tisserand parameter <cit.>.

In this picture, there are two ways to drive an initially stationary pendulum to circulation: one is to perturb the ball of the pendulum directly (thereby changing the energy-level of the trajectory), and the other is to laterally rock the base (thus modulating the separatrix along the Φ-axis). These processes are directly equivalent to the two types of diffusion considered in our calculations. That is, [1] diffusion in the dynamic variables x and y themselves (explicitly connected to eccentricities and resonant angle) is analogous to direct perturbations to the ball of the pendulum, while [2] diffusion in the proximity parameter ε (explicitly connected to the semi-major axes) corresponds to shaking the base of the pendulum back and forth.

Meanwhile, consequences of eccentricity damping and convergent migration are equivalent to friction that acts to return the ball of the pendulum back to its undisturbed state, and restore the separatrix to its equilibrium position, respectively. In the type-I migration regime however, eccentricity damping by the disk is far more efficient than orbital decay <cit.>, meaning that the ball of the pendulum is effectively submerged in water, while the base of the pendulum is only subject to air-resistance (in this analogy). As a result, the latter process — diffusion in proximity parameter ε — ends up being more important for purposes of moving planets out of mean motion resonance (see equation [<ref>]).























































 §.§ Discussion

 




The work presented herein suggests that turbulent forcing is unlikely to be the single dominant effect that sculpts the final orbital distribution of exoplanets. At the same time, the functional form of expression (<ref>) yields important insight into the evolutionary aspects of the planet formation process. Particularly, because the resonance disruption criterion
depends on the disk mass, it implies a certain time-dependence of the
mechanism itself (as the nebula dissipates, the critical mass ratio below which the mechanism operates
decreases from a value substantially above the Earth-Sun mass
ratio, to one below). This means that even though the turbulent disruption mechanism becomes ineffective
in a weaning nebula, it may be key to facilitating growth in the early stages of evolution of
planetary systems, by allowing pairs of proto-planets to skip over mean-motion commensurabilities and merge, instead of forming resonant chains. In essence, this type of dynamical behavior is seen in the large-scale numerical experiments of <cit.>.

For much of this work, the system parameters that we use effectively
assume a maximum rate of orbital convergence. Because the quantitative
nature of migration can change substantially in the inner nebula, the
actual rate of orbital convergence may be somewhat lower
<cit.>.  This change would make planetary resonances more
susceptible to stochastic disruption. At the same time, we have not
taken into account the inhibition of the random gravitational field
through non-ideal magnetohydrodynamic effects <cit.>,
which would weaken the degree of stochastic forcing. Both of these
effects can be incorporated into the criterion of equation
(<ref>) by lowering the migration factor f and the
value of α accordingly.  However, because both of these quantities
appear under a square root in the expression, the sensitivity
of our results to these corrections is not expected to be extreme.

This work assumes that turbulence operates in circumstellar
disks at the expected levels. The presence of turbulence is most commonly attributed to the magneto-rotational-instability <cit.>,
which in turn requires the disk to be sufficiently ionized. Although
the innermost regions of the disk are expected to be ionized by
thermal processes, dead zones could exist in intermediate part of the
disk (; see also ), and ionization by cosmic rays can be suppressed
in the outer disk <cit.>. Indeed, suppressed levels of
ionization are now inferred from ALMA observations of young star/disk
systems <cit.>, implying that the assumption of sufficient
ionization — and hence active MRI turbulence — is not guarenteed. At the same time, our model is agnostic towards the origins of turbulent fluctuations themselves, and can be employed equally well if a purely hydrodynamic source of turbulence were responsible for angular momentum transport within the nebula <cit.>.


In light of the aforementioned uncertainties inherent to the problem at hand, it is of considerable interest to explore if simply adjusting the parameters can, in principle, yield consistency between the model and the observations. That is, can reasonable changes to the migration rate, etc., generate agreement between the turbulent resonance disruption hypothesis and the data? Using equation (<ref>), we find that increasing the local surface density by an order of magnitude (Σ=10⟨Σ⟩=170,000g/cm^2) while lowering the orbital convergence rate a hundred-fold (f=0.01) and retaining h/r=0.05, ⟨ a⟩=0.1AU, α=0.01 yields (m_1+m_2)/M≃2× 10^-4∼60M_⊕/M_⊙ as the critical mass ratio, thus explaining the full range of values shown in Figure <ref>. Correspondingly, rough agreement between observations and the stochastic migration scenario is reproduced in the work of <cit.>, where the amplitude of turbulent forcing was tuned to give consistency with data.

Although this line of reasoning may appear promising, it is important to note that as the disk accretes onto the star, the local surface density will diminish, causing the critical mass ratio to decrease as well. Meanwhile, even with a reduction factor of f=0.01, the type-I migration timescale remains shorter than the ∼few Myr lifetime of the nebula, as long as Σ≳0.1⟨Σ⟩=170g/cm^2. As a result, we argue that any realistic distribution of the assumed parameters is unlikely to allow turbulence to provide enough resonance disruption to explain the entire set of observations.

If disk turbulence does not play the defining role in generating an
observational census of extrasolar planets that is neither dominated
by, nor devoid of, mean motion resonances, than what additional
processes are responsible for the extant data set? As already
mentioned in the introduction, there are two other ways in which
planets can avoid resonant locking – resonant metastability <cit.> and capture
probability suppression <cit.>.  The first mechanism requires that the outer
planet is more massive then the inner planet to compromise
resonance <cit.>. As a result, observed resonant systems would almost always
have a more massive inner planet, but this ordering is not reflected
in the data. On the other hand, the (second) capture suppression
mechanism requires disk eccentricities of order
e_disk∼0.02 to explain the data. Importantly, disk eccentricities of
this magnitude (and greater) are not only an expected result of
theoretical calculations, they are invoked to explain observations of
asymmetric glow of dust <cit.>.

In conclusion, turbulent fluctuations probably do not explain the
entire ensemble of observed planetary systems, which exhibit only a weak preference for mean motion commensurability. In addition to
turbulent forcing, many other physical processes are likely at work,
where perhaps the most promising mechanism is capture suppression due to nonzero
disk eccentricities. Nonetheless, a subset of exotic planetary systems
that exhibit large-amplitude resonant librations likely require a
turbulent origin. The relative duty cycle of this mechanism, and
others, poses an interesting problem for further exploration.



Acknowledgments: We would like to thank Juliette Becker, Tony Bloch, Wlad Lyra and Chris Spalding for useful discussions, as well as the referee, Hanno Rein, whose insightful report led to a considerable improvement of the manuscript. K.B. acknowledges support from the NSF AAG program AST1517936, and from Caltech. F.C.A. acknowledges support from the NASA Exoplanets Research Program NNX16AB47G, and from the University of Michigan. 



 

[Adams et al.(2008)]alb
Adams, F. C., Laughlin, G., & Bloch, A. M. 2008, ApJ, 683, 1117

[Allan(1969)]Allan1969 
Allan, R. R. 1969, AJ, 74, 497

[Allan(1970)]Allan1970
Allan, R. R. 1970, Celest. Mech., 2, 121

[Balbus & Hawley(1991)]mri 
Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

[Bai & Stone(2013)]2013ApJ...769...76B Bai, X.-N., & Stone, J. M. 2013, , 769, 76 

[Batalha et al.(2013)]Batalha 
Batalha, N. M. et al. 2013, ApJS, 204, 24

[Batygin(2015)]Batygin2015 
Batygin, K. 2015, MNRAS, 451, 2589

[Batygin et al.(2015)]BatDeckHol 
Batygin, K., Deck, K. M., & Holman, M. J. 2015, AJ, 149, 167

[Batygin et al. (2011)]batynmorbid
Batygin, K., & Morbidelli, A. 2011, Celest. Mech., 111, 219

[Batygin & Morbidelli(2013)]BatMorby2013b
Batygin, K., & Morbidelli, A. 2013, A&A, 556, A28

[Bitsch & Kley(2011)]Bitsch
Bitsch, B., & Kley, W. 2011, A&A, 536, A77

[Bitsch et al.(2015)]2015A   A...575A..28B Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, , 575, A28 

[Borders & Goldreich(1984)]BordersGoldreich1984
Borderies, N., & Goldreich, P. 1984, Celest. Mech., 32, 127

[Chiang & Laughlin(2013)]ChiangLaughlin
Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444 

[Cleeves et al.(2013)]cleevesone 
Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5 

[Cleeves et al.(2015)]cleeves 
Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I.
2015, ApJ, 799, 204 

[Coleman & Nelson(2016)]ColemanNelson
Coleman, G.A.L., & Nelson, R. P. 2016, MNRAS, 457, 2480

[Cresswell & Nelson(2008)]CresswellNelson2008 
Cresswell, P., & Nelson, R. P. 2008, A&A, 482, 677

[Crida et al.(2006)]2006Icar..181..587C Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587 

[Crida et al.(2008)]Crida2008
Crida, A., Sándor, A., & Kley, W., 2008, A&A, 483, 325 

[D'Angelo & Bodenheimer(2016)]DAngeloBodenheimer2016 D'Angelo, G., & Bodenheimer, P. 2016, , 828, 33 

[Deck et al.(2012)]Decketal 
Deck, K. M. et al. 2012, ApJ, 755, 21  

[Deck et al.(2013)]Deck2013
Deck, K. M., Payne, M., & Holman, M. J. 2013, ApJ, 774, 129 

[Deck & Batygin(2015)]DeckBatygin 
Deck, K. M., & Batygin, K. 2015, ApJ, 810, 119 

[Duffell & MacFadyen(2013)]2013ApJ...769...41D Duffell, P. C., & MacFadyen, A. I. 2013, , 769, 41 


[Duffell & Chiang(2015)]DuffellChiang2015 Duffell, P. C., & Chiang, E. 2015, , 812, 94 


[Fabrycky et al.(2014)]Fabrycky2014
Fabrycky, D. C. et al. 2014, ApJ, 790, 146

[Foreman-Mackey et al.(2014)]Forman-Mackey 
Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ,
795, 64 

[Fressin et al.(2013)]Fressin2013 Fressin, F., Torres, G., Charbonneau, D., et al. 2013, , 766, 81 

[Gammie(1996)]gammie 
Gammie, C. F. 1996, ApJ, 457, 355

[Goldreich(1965)]Goldreich1965 
Goldreich, P. 1965, MNRAS, 130, 159

[Goldreich & Schlichting(2014)]GoldShicht2014 
Goldreich, P., & Schlichting, H. E. 2014, AJ, 147, 32

[Goldreich & Tremaine(1980)]GoldreichTremaine1980 
Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425

[Goldreich & Tremaine(1982)]GoldreichTremaineRing  
Goldreich, P., & Tremaine, S. 1982, ARA&A, 20, 249  

[Hansen & Murray(2015)]MurrayHansen
Hansen, B.M.S., & Murray, N. 2015, MNRAS, 448, 1044 

[Hayashi(1981)]1981PThPS..70...35H Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35 

[Henrard(1986)]Henrard1986 
Henrard, J., & Lemaitre, A. 1986, Celest. Mech., 39, 213 

[Henrard & Lamaitre(1983)]HenrardLemaitre1983 Henrard, J., & Lamaitre, A. 1983, Celestial Mechanics, 30, 197 

[Horn et al.(2012)]lyra 
Horn, B., Lyra, W., Mac Low, M-M., & Sándor, Z. 2012, ApJ, 750, 34  

[Howard et al.(2012)]Howard
Howard, A. W. et al. 2012, ApJS, 201, 15

[Johnson et al.(2006)]johnson 
Johnson, E. T., Goodman, J., & Menou, K. 2006, ApJ, 647, 1413

[Jontof-Hutter et al.(2014)]2014ApJ...785...15J Jontof-Hutter, D., Lissauer, J. J., Rowe, J. F., & Fabrycky, D. C. 2014, , 785, 15 


[Ketchum et al.(2011)]ketchum 
Ketchum, J. A., Adams, F. C., & Bloch, A. M. 2011, ApJ, 726, 53   

[Kley & Nelson(2012)]KleyNelson2012
Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211

[Laughlin et al.(2004)]lsa
Laughlin, G., Steinacker, A., & Adams, F. C. 2004, ApJ, 608, 489

[Lecoanet et al.(2009)]lecoanet
Lecoanet, D., Adams, F. C., & Bloch, A. M. 2009, ApJ, 692, 659  

[Lee & Peale(2002)]leepeale
Lee, M.-H., & Peale, S. J. 2002, ApJ, 567, 596

[Lee & Chiang(2015)]Lee2015 
Lee, E. J., & Chiang, E. 2015, ApJ, 811, 41 

[Lee & Chiang(2016)]Lee2016 
Lee, E. J., & Chiang, E. 2016, ApJ, 817, 90

[Lin & Youdin(2015)]2015ApJ...811...17L Lin, M.-K., & Youdin, A. N. 2015, , 811, 17 

[Malhotra(1993)]Malhotra1993
Malhotra, R. 1993, Nature, 365, 819

[Mestel(1963)]Mestel63
Mestel, L. 1963, MNRAS, 126, 553

[Mills et al.(2016)]2016Natur.533..509M Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, , 533, 509 

[Mittal & Chiang(2015)]2015ApJ...798L..25M Mittal, T., & Chiang, E. 2015, , 798, L25 

[Morbidelli(2002)]Morby 
Morbidelli, A. 2002, Modern Celestial Mechanics: Aspects of Solar 
System Dynamics (London: Taylor & Francis) 

[Mulders et al.(2015)]Mulders2015 Mulders, G. D., Pascucci, I., & Apai, D. 2015, , 798, 112

[Murray & Dermott(1999)]md99
Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics 
(Cambridge: Cambridge Univ. Press) 

[Nelson & Papaloizou(2004)]nelson
Nelson R. P., & Papaloizou J.C.B. 2004, MNRAS, 350, 849

[Nelson et al.(2013)]2013MNRAS.435.2610N Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, , 435, 2610 

[Nesvorný & Morbidelli(2008)]morbyone  
Nesvorný, D., & Morbidelli, A. 2008, Icarus, 688, 636 

[Nobili & Roxburgh(1986)]genrel Nobili, A., & Roxburgh, I. W. 1986, IAUS, 114, 105

[Ogihara et al.(2015)]Ogihara2015insitu Ogihara, M., Morbidelli, A., & Guillot, T. 2015, , 578, A36 

[Ogihara et al.(2015)]Ogihara2015mig Ogihara, M., Morbidelli, A., & Guillot, T. 2015, , 584, L1 


[Okuzumi & Hirose(2011)]okuzumi 
Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65 

[Okuzumi & Ormel(2013)]OkuzumiOrmel2013 Okuzumi, S., & Ormel, C. W. 2013, , 771, 43 


[Öpik(1976)]Opik1976
Opik, E. J. 1976, Interplanetary Encounters: 
Close-range gravitational interactions (New York: Elsevier) 

[Ormel & Okuzumi(2013)]ormelpaper
Ormel, C. W., & Okuzumi, S. 2013, ApJ, 771, 44

[Paardekooper(2014)]Paardekooper Paardekooper, S.-J. 2014, , 444, 2031 

[Papaloizou & Larwood(2000)]PapaloizouLarwood2000
Papaloizou, J.C.B., & Larwood, J. D. 2000, MNRAS, 315, 823

[Peale(1976)]Peale1976
Peale, S. J. 1976, ARA&A, 14, 215

[Petigura et al.(2013)]Petigura 
Perigura, E. A., Howard, A. W., & Marcy, G. W. 2013, PNAS, 110, 19273

[Press et al.(1992)]Press1992 
Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery,
B. P. 1992, Numerical Recipes in FORTRAN: The Art of Scientific
Computing (Cambridge: Cambridge Univ. Press)

[Quillen(2006)]quill
Quillen, A. C. 2006, MNRAS, 365, 1367 

[Rein & Papaloizou(2009)]rein 
Rein, H., & Papaloizou, J.C.P. 2009, A&A, 497, 595  

[Rein et al.(2010)]2010A   A...510A...4R Rein, H., Papaloizou, J. C. B., & Kley, W. 2010, , 510, A4 

[Rein(2012)]Rein2012 Rein, H. 2012, , 427, L21 

[Rivera et al.(2010)]2010ApJ...719..890R Rivera, E. J., Laughlin, G., Butler, R. P., et al. 2010, , 719, 890 




[Rogers(2015)]Rogers2015 
Rogers, L. A. 2015, ApJ, 801, 41 

[Schlichting(2014)]Schlichting
Schlighting, H. 2014, ApJ, 795, 15 

[Sessin & Ferraz-Mello(1984)]SessinFerraz-Mello1984
Sessin, W., & Ferraz-Mello, S. 1984, Celest. Mech., 32, 307

[Sinclair(1970)]Sinclair1970
Sinclair, A. T. 1970, MNRAS, 148, 325

[Sinclair(1972)]Sinclair1972
Sinclair, A. T. 1972, MNRAS, 160, 169

[Shakura & Sunyaev(1973)]ShakuraSunayev1973
Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
 
[Tanaka et al.(2002)]Tanaka2002
Tanaka, H., Takeuchi, T., Ward, W. R. 2002, ApJ, 565, 1257

[Tanaka & Ward(2004)]Tanaka2004 
Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388 

[Terquem & Papaloizou(2007)]TerquemPap2007
Terquem, C., Papaloizou, J.C.B. 2007, ApJ, 654, 1110

[Van Kampen(2001)]vankampen
Van Kampen, N. G. 2001, Stochastic Processes in Physics 
and Chemistry (Amsterdam: North Holland) 

[Weiss & Marcy(2014)]WeissMarcy2014 
Weiss, L. M., & Marcy, G. W. 2014, ApJ, 783, 6 

[Winn & Fabrycky(2015)]WinnFabrycky2015 Winn, J. N., & Fabrycky, D. C. 2015, , 53, 409 

[Xu & Lai(2016)]2016arXiv161106463X Xu, W., & Lai, D. 2016, arXiv:1611.06463 

[Wisdom(1986)]Wisdom1986
Wisdom, J. 1986, Celest. Mech., 38, 175

 

