An Analytic Criterion for Turbulent Disruption of Planetary Resonances

TL;DR

This paper derives an analytic criterion to determine when turbulent density fluctuations in protoplanetary disks can break mean motion resonances between planets, supported by numerical simulations, highlighting the importance of planet-star mass ratio.

Contribution

The authors present a new simple analytic criterion for turbulent resonance disruption, validated by numerical integrations and N-body simulations, linking physical parameters to resonance breaking conditions.

Findings

Resonance disruption depends strongly on planet-star mass ratio.

Only low-mass planet pairs (less than about 3 Earth masses) are susceptible to turbulence-induced resonance breaking.

Additional mechanisms are needed to explain the prevalence of non-resonant exoplanet systems.

Abstract

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 $\alpha=0.01$, only planet pairs with cumulative mass ratios smaller than $(m_1+m_2)/M\lesssim10^{-5}\sim3M_{\oplus}/M_{\odot}$ are susceptible to breaking resonance at semi-major axis of order $a\sim0.1\,$AU. 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.

Figures

Equations

$$\displaystyle\mathcal{H}=3\,\big{(}\varepsilon+1\big{)}\bigg{(}\frac{x^{2}+y^{2}}{2}\bigg{)}-\bigg{(}\frac{x^{2}+y^{2}}{2}\bigg{)}^{2}-2\,x,$$

$$\sigma =\sqrt{e_1^2+e_2^2-2e_1{e}_2\cos (\Delta \varpi )},$$

$$\displaystyle[A]=\frac{1}{2}\bigg{(}\frac{15}{4}\frac{k\,M}{m_{1}+m_{2}}\bigg{)}^{2/3}$$

$$\displaystyle[T]=\frac{1}{n}\bigg{(}\frac{5}{\sqrt{6}\,k^{2}}\frac{M}{m_{1}+m_{2}}\bigg{)}^{2/3},$$

$$x$$

$$y$$

$$\epsilon$$

$$\tilde{\omega }\equiv \arctan \left[\frac{e_2\sin {\varpi }_2-e_1\sin {\varpi }_1}{e_1\cos {\varpi }_1-e_2\cos {\varpi }_2}\right]$$

$$\displaystyle\tau_{\rm{wave}}=\frac{1}{n}\bigg{(}\frac{M}{m}\bigg{)}\bigg{(}\frac{M}{\Sigma\,a^{2}}\bigg{)}\bigg{(}\frac{h}{r}\bigg{)}^{4},$$

$$\displaystyle\frac{1}{a}\frac{d\,a}{d\,t}\equiv-\frac{1}{\tau_{\rm{mig}}}\simeq-\frac{4\,f}{\tau_{\rm{wave}}}\bigg{(}\frac{h}{r}\bigg{)}^{2}$$

$$\frac{1}{e}\frac{de}{dt}\equiv -\frac{1}{{\tau }_{\mathrm{dmp}}}\simeq -\frac{3}{4}\frac{1}{{\tau }_{\mathrm{wave}}}.$$

$${\mathcal{D}}_{\xi }=\frac{{\mathcal{D}}_a}{a^2}\sim 2{\mathcal{D}}_e\sim \frac{\alpha }{2}{\left(\frac{\Sigma a^2}{M}\right)}^{2}n,$$

$$\displaystyle\frac{d\chi}{dt}\simeq\frac{3}{2\,\langle a\rangle}\bigg{(}\frac{da_{1}}{dt}-\frac{da_{2}}{dt}\bigg{)},$$

$$d\chi =\frac{3}{2}\sqrt{2{\mathcal{D}}_{\xi }}dw-\frac{3}{2}\frac{\chi }{{\tau }_{\mathrm{mig}}}dt,$$

$$\delta \chi =\sqrt{\frac{3{\mathcal{D}}_{\xi }{\tau }_{\mathrm{mig}}}{2}}=\frac{1}{4}\frac{h}{r}\sqrt{\frac{3\alpha \Sigma ⟨a{⟩}^2}{f(m_1+m_2)}}.$$

$$\displaystyle\Delta\chi\simeq 5\bigg{[}\frac{\sqrt{k}\,(m_{1}+m_{2})}{M}\bigg{]}^{2/3}.$$

$$\frac{\delta \chi }{\Delta \chi }$$

$$\displaystyle\times\Bigg{[}\,\frac{\Sigma\,\langle a\rangle^{2}}{k\,M}\sqrt{\frac{\Sigma\,\langle a\rangle^{2}}{m_{1}+m_{2}}}\,\Bigg{]}^{1/3}\gtrsim 1.$$

$$\displaystyle\mathcal{D}_{x}\simeq\alpha\,\bigg{(}\frac{\Sigma\,\langle a\rangle^{2}}{M}\bigg{)}^{2}\bigg{(}\frac{M}{\sqrt{k}\,(m_{1}+m_{2})}\bigg{)}^{4/3}.$$

$$\displaystyle\tau_{x}\simeq\,\bigg{(}\frac{128}{225}\frac{k\,M}{m_{1}+m_{2}}\bigg{)}^{1/3}\frac{k\,M}{\Sigma\,\langle a\rangle^{2}}\bigg{(}\frac{h}{r}\bigg{)}^{4},$$

$$dx=\sqrt{2{\mathcal{D}}_x}dw-\frac{x}{{\tau }_x}dt,$$

$$\frac{\delta x}{\Delta x}$$

$$\displaystyle\times\bigg{(}\frac{M}{m_{1}+m_{2}}\bigg{)}^{5/6}\gtrsim 1\,.$$

$$\displaystyle\frac{\delta x/\Delta x}{\delta\chi/\Delta\chi}\sim{5}{\sqrt{f}}\,\bigg{(}\frac{h}{r}\bigg{)}\,\bigg{[}\frac{k^{2}\,(m_{1}+m_{2})}{M}\bigg{]}^{1/3}\ll 1.$$

$$\displaystyle\frac{dx}{dt}=-3\,y\,\big{(}1+\varepsilon\big{)}+y\,\big{(}x^{2}+y^{2}\big{)}-\frac{x}{\tau_{\rm{dmp}}/[T]}$$

$$\displaystyle\frac{dy}{dt}=-2+3\,x\,\big{(}1+\varepsilon\big{)}-x\,\big{(}x^{2}+y^{2}\big{)}-\frac{y}{\tau_{\rm{dmp}}/[T]}$$

$$\displaystyle\frac{d\varepsilon}{dt}=\frac{2}{3}\bigg{(}\frac{[A]}{k\,\tau_{\rm{mig}}/[T]}-\frac{x^{2}+y^{2}}{\tau_{\rm{dmp}}/[T]}\bigg{)}+\mathcal{F}.$$

$$\frac{d\vec{v}}{dt}=-\frac{\vec{v}}{{\tau }_{\mathrm{mig}}}-\frac{2\vec{r}}{{\tau }_{\mathrm{dmp}}}\frac{\left(\vec{v}\cdot \vec{r}\right)}{\left(\vec{r}\cdot \vec{r}\right)},$$

Full text

1. 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 (Howard et al., 2012; Petigura et al., 2013; Batalha et al., 2013), 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 (Fressin et al., 2013; Foreman-Mackey et al., 2014; Mulders et al., 2015).

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 (Hansen & Murray, 2015; Chiang & Laughlin, 2013; Lee & Chiang, 2015, 2016), or if they migrate to their close-in orbits having formed at large orbital radii, as a consequence of disk-planet interactions (Goldreich & Tremaine, 1980; Tanaka et al., 2002; Crida et al., 2008; Kley & Nelson, 2012). Although this question remains a subject of active research, a number recent studies (Schlichting, 2014; Ogihara et al., 2015; D’Angelo & Bodenheimer, 2016) have pointed to a finite extent of migration as an apparent requirement for successful formation of Super-Earths. Moreover, structural models (Rogers, 2015) 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 (Goldreich, 1965; Allan, 1969, 1970; Sinclair, 1970, 1972). While distinct examples of resonant planetary systems exist within the known aggregate of planets111Archetypical examples of short-period resonant systems include GJ 876 (Rivera et al., 2010), Kepler-36 (Deck et al., 2012), Kepler-79 (Jontof-Hutter et al., 2014), and Kepler-227 (Mills et al., 2016)., the overall orbital distribution shows little preference for mean motion commensurabilities (Figure 1). 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 (Fabrycky et al., 2014). 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 (Lee & Peale, 2002; Quillen, 2006; Terquem & Papaloizou, 2007; Cresswell & Nelson, 2008). An important distinction was drawn by the work of Adams et al. (2008), 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 (Lecoanet et al., 2009; Ketchum et al., 2011; Horn et al., 2012), and that in specific cases, turbulence can be conducive to the reproduction of dynamical architecture (Rein & Papaloizou, 2009; Rein et al., 2010).

While the prediction of the infrequency of resonant systems made by Adams et al. (2008) 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 Goldreich & Schlichting (2014) proposed that a particular relationship between the rates of eccentricity damping and semi-major axis decay can render resonances metastable, while Batygin (2015) 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 (Deck & Batygin, 2015), while the capture suppression mechanism requires disk eccentricities on the order of a few percent to operate (Batygin, 2015). 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., Ketchum et al. 2011). 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 2, we present the details of our model. In Section 3, we employ methods from stochastic calculus to derive an analytic criterion for turbulent disruption of mean motion resonances. In Section 4, we confirm our results with both perturbative numerical integrations and an ensemble of full $N$-body simulations. The paper concludes in Section 5 with a summary of our results and a discussion of their implications.

2. 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 2. Throughout much of the manuscript, we make the so-called “compact” approximation, where we assume that the semi-major axis ratio $\xi\equiv a_{1}/a_{2}\rightarrow 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\geqslant 3$ (see, e.g., Deck et al. 2013; Deck & Batygin 2015), where the integer $k$ specifies the resonance (Murray & Dermott, 1999; Morbidelli, 2002).

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.

2.1. 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 (Sessin & Ferraz-Mello, 1984; Wisdom, 1986; Henrard, 1986). More recently, this formalism has been used to provide a geometric representation of resonant dynamics (Batygin & Morbidelli, 2013), study the onset of chaos (Deck et al., 2013), generalize the theory of resonant capture (Batygin, 2015), as well as to elucidate overstable librations (Deck & Batygin, 2015). 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 (Quillen, 2006). 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:

[TABLE]

where the variables $(\varepsilon,x,y)$ are defined below. A Hamiltonian of this form is typically referred to as the second fundamental model for resonance (Henrard & Lamaitre, 1983; Borders & Goldreich, 1984), and behaves as a forced harmonic oscillator at negative values of the proximity parameter, $\varepsilon$, while possessing a pendulum-like phase-space structure at large positive values of $\varepsilon$. 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., Peale 1976; Malhotra 1993; Deck et al. 2013). 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 $\varepsilon$ 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

[TABLE]

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

[TABLE]

where $m$ is planetary mass, $M$ is stellar mass, and $n=\sqrt{\mathcal{G}\,M/a^{3}}$ is the mean motion. Then, in the compact limit, the variables in the Hamiltonian become (Deck & Batygin, 2015):

[TABLE]

where $\xi=a_{1}/a_{2}$, and the quantity

[TABLE]

represents a generalized longitude of perihelion.

The specification of resonant dynamics is now complete. While application of Hamilton’s equations to equation (1) only yields the evolution of $\sigma$ and the corresponding resonant angle, the behavior of the individual eccentricities and apsidal lines can be obtained from the conserved222When the system is subjected to slow evolution of the proximity parameter $\varepsilon$, $\rho$ is no longer a strictly conserved quantity. Instead, $\rho$ becomes an adiabatic invariant that is nearly constant, except when the system encounters a homoclinic curve (Batygin & Morbidelli, 2013). quantity $\rho=m_{1}\,e_{1}^{2}+m_{2}\,e_{2}^{2}+m_{1}\,m_{2}\,e_{1}\,e_{2}\,\cos(\Delta\varpi)$. In addition, we note that the definitions of the variables (4) 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 $\xi\equiv a_{1}/a_{2}\rightarrow 1$, and is qualitatively equivalent to the Öpik approximation (Öpik, 1976).

2.2. Planet-Disk Interactions

Dating back to early results on ring-satellite interactions (Goldreich & Tremaine, 1982), 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:

[TABLE]

where $\Sigma$ 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 (Mestel, 1963), the corresponding rates of eccentricity and semi-major axis decay are given by (Tanaka et al., 2002; Tanaka & Ward, 2004):

[TABLE]

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 $\mathcal{K}$. To this end, we note that the value of $\mathcal{K}\sim 10^{2}$ adopted by many previous authors (Lee & Peale, 2002; Ketchum et al., 2011) is in rough agreement with equation (7) which yields $\mathcal{K}\sim(h/r)^{-2}$.

While eccentricity damping observed in numerical simulations (e.g., Cresswell & Nelson 2008) is well matched by equation (7), 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 (Bitsch & Kley, 2011; Paardekooper, 2014). Although such corrections alter the migration histories on a detailed level, convergent migration followed by resonant locking remains an expected result in laminar disks (Coleman & Nelson, 2016). For simplicity, in this work, we account for this complication by introducing an adjustable parameter $f$ into equation (7).

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 (Adams et al., 2008). Such perturbations translate to effectively diffusive evolution of the eccentricity and semi-major axis (Laughlin et al., 2004; Nelson & Papaloizou, 2004). 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 Johnson et al. 2006; Adams et al. 2008; Okuzumi & Ormel 2013) to obtain the expressions

[TABLE]

where $\alpha$ is the Shakura-Sunayev viscosity parameter (Shakura & Sunyaev, 1973). Although non-ideal effects can modify the above expressions on the quantitative level (Okuzumi & Hirose, 2011), 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 $\alpha$ accordingly.

3. 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 $\chi\equiv n_{2}/n_{1}-k/(k-1)$ and $x$. Keep in mind that this latter quantity is directly proportional to the generalized eccentricity $\sigma$ (see equation [4]).

3.1. Diffusion of Semi-major Axes

In the compact limit $a_{1}\approx a_{2}$, the time evolution of the parameter $\chi$ can be written in the approximate form

[TABLE]

where $\langle a\rangle$ 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 $\mathcal{D}_{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 (Rein & Papaloizou, 2009). Additionally, for comparable-mass planets, we may adopt $\tau_{\rm{mig}}$ as a characteristic drift rate, replacing $m$ with $m_{1}+m_{2}$ in equation (7). 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

[TABLE]

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

[TABLE]

The approximate extent of stochastic evolution that the system can experience and still remain in resonance is given by the resonant bandwidth, $\Delta\chi$. At its inception333A resonance can only be formally defined when a homoclinic curve (i.e., a separatrix) exists in phase-space. For a Hamiltonian of the form (1), a separatrix appears at $\varepsilon=0$, along with an unstable (hyperbolic) fixed point, that bifurcates into two fixed points (one stable and one unstable) at $\varepsilon>0.$, the width of the resonance (Batygin, 2015) is given by

[TABLE]

Accordingly, a rough criterion for turbulent disruption of the resonance is

[TABLE]

Keep in mind that $\delta\chi$ is a measure of the width of the distribution in the variable $\chi$ due to stochastic evolution, whereas $\Delta\chi$ is the change in $\chi$ 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 3, where the expression (13) is shown as a function of the quantity $(m_{1}+m_{2})/M$, assuming system properties $\alpha=10^{-2}$, $h/r=0.05$, $\langle a\rangle=0.1\,$AU, $f=1$, and $k=3$. The disk profile is taken to have the form $\Sigma=\Sigma_{0}\,(r_{0}/r)$, with $\Sigma_{0}=1700\,\rm{g/cm}^{2}$ and $r_{0}=1\,$AU, such that the local surface density at $\langle a\rangle$ is $\langle\Sigma\rangle=17,000\,\rm{g/cm}^{2}$. Notice that the disruption criterion (13) 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,\dots 0.9,1\,\times\,\langle\Sigma\rangle$) are also shown, and color-coded accordingly.

While Figure 3 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 (Paardekooper, 2014). To this end, we note that setting $f=0$ in equation (13) yields $\infty>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 (Peale, 1976).

3.2. 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 $\sigma$ diffuses with the coefficient $\sqrt{2}\,\mathcal{D}_{e}$. Accounting for conversion factors between conventional quantities and the dimensionless coordinates (given by equation [3]), we obtain

[TABLE]

Similarly, the damping timescale takes the form

[TABLE]

where, as before, we adopted the total planetary mass as an an approximation for $m$ in the expression (7).

In direct analogy with equation (10), we obtain the stochastic equation for the time evolution of $x$,

[TABLE]

so that the distribution of $x$ is characterized by the standard deviation $\delta x=\sqrt{\mathcal{D}_{x}\,\tau_{x}}$. At the same time, we take the half-width of the resonant separatrix to be given by $\Delta x=2$ (e.g., Batygin & Morbidelli 2013; Deck et al. 2012). Combining these two results, we obtain a second criterion for resonance disruption, i.e.,

[TABLE]

3.3. 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 (13) and (17) 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

[TABLE]

The fact that this expression evaluates to a number substantially smaller than unity means that diffusion in semi-major axes (equation [13]) dominates over diffusion in eccentricities (equation [17]) as a mechanism for disruption of mean-motion commensurabilities. Although the relative importance of $\mathcal{D}_{a}$ compared to $\mathcal{D}_{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.

4. 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 (13) through numerical integrations. We first present a perturbative approach (Section 4.1) and then carry out a series of full $N$-body simulations (Section 4.2).

4.1. Perturbation Theory

The dynamical system considered here is described by three equations of motion, corresponding to the variations in $x$, and $y$, and $\varepsilon$. Although the resonant dynamics itself is governed by Hamiltonian (1), 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 $\tau_{\rm{dmp}}$ as the decay timescale for the generalized eccentricity, $\sigma$, and take $\tau_{\rm{mig}}$ as the characteristic orbital convergence time. The full equations of motion are then given by:

[TABLE]

In the above expression, $\mathcal{F}$ 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\,\zeta\sin(\pi\,t/\Delta t)/\Delta t$, where $\zeta$ is a Gaussian random variable. The pulse time interval was taken to be $\Delta t=0.1$, and the standard deviation of $\zeta$ was chosen such that the resulting diffusion coefficient $\mathcal{D}_{\zeta}=\sigma_{\zeta}^{2}/\Delta t$ matched that given by equation (8).

Note that here, we have opted to only implement stochastic perturbations into the equation that governs the variation of $\varepsilon$. 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 (18).

Turbulent fluctuations aside, the equation of motion for the parameter $\varepsilon$ indicates that there exists an equilibrium value of the generalized eccentricity $\sigma_{\rm{eq}}=\sqrt{\tau_{\rm{dmp}}/(2\,k\,\tau_{\rm{mig}})}$ that corresponds to stable capture into resonance. Analogously, the equilibrium value of $x_{\rm{eq}}=\sigma_{\rm{eq}}\sqrt{2[A]}$ parallels the strictly real fixed point of Hamiltonian (1). As a result, if we neglect the small dissipative contributions and set $dx/dt=0,dy/dt=0,x=x_{\rm{eq}},y=0$ in the first and second equations in expression (19), we find an equilibrium value of the proximity parameter, $\varepsilon_{\rm{eq}}$, that coincides with resonant locking. An ensuing crucial point is that if resonance is broken, the system will attain values of $\varepsilon$ substantially above the equilibrium value $\varepsilon_{\rm{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 3. In particular, we adopted $\alpha=10^{-2}$, $h/r=0.05$, $\langle a\rangle=0.1\,$AU, $f=1$, and $k=3$. Additionally, we again chose a surface density profile with $\Sigma_{0}=1700\,$g/cm2 at $r_{0}=1\,$AU, that scales inversely with the orbital radius, such that the nominal surface density at $r=\langle a\rangle$ is $\langle\Sigma\rangle=17,000\,$g/cm2. We also performed a series of simulations that span a lower range of surface densities ($0.1,0.2,\dots,0.9,1\,\times\,\langle\Sigma\rangle$). All of the integrations were carried out over a time span of $\tau_{\rm{mig}}=100\,\tau_{\rm{wave}}$, with the system initialized at zero eccentricity ($x_{0}=y_{0}=0$), on orbits exterior to exact commensurability ($\epsilon_{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 3, the qualitative expectations for the outcomes of these simulations (as dictated by equation [13]) 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 4 depicts numerically computed evolution of $\varepsilon$ for the full range of local surface densities under consideration (color-coded in the same way as in Figure 3) as a function of time. These numerical results are in excellent agreement with our theoretical expectations from Section 3. The proximity parameter always approaches its expected equilibrium value $\varepsilon_{\rm{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 $\Sigma=0.7\,\langle\Sigma\rangle$ and $\Sigma=0.9\,\langle\Sigma\rangle$ show the system breaking out of resonance within a single orbital convergence time, $\tau_{\rm{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 5 shows the phase-space counterpart of the evolution depicted in Figure 4. Specifically, the $x$-$y$ projections of the system dynamics are shown for cases with surface densities $\Sigma=0.3\,\langle\Sigma\rangle$ (blue) and $\Sigma=\langle\Sigma\rangle$ (black), where the background depicts the topology of the Hamiltonian (1). In each panel, the black curve designates the separatrix of $\mathcal{H}$, given the equilibrium value of the proximity parameter $\varepsilon=\varepsilon_{\rm{eq}}$. The background color scale is a measure of the value of $\mathcal{H}$. The thee equilibrium points of the Hamiltonian are also shown, as transparent green dots.

As in Figure 4, 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 $\Sigma=\langle\Sigma\rangle$. 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.

4.2. $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 (13), we have carried out a series of direct $N$-body simulations. The integrations utilized a Burlisch-Stoer integration scheme (e.g., Press et al. 1992) 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 Papaloizou & Larwood (2000), we incorporated the orbital decay and eccentricity damping using accelerations of the form:

[TABLE]

where $\vec{v}$ and $\vec{r}$ denote the orbital velocity and radius respectively.444Note 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.1\,$AU at every time step555Qualitatively, this procedure is equivalent to changing the unit of time at every time step (Deck & Batygin, 2015).. The characteristic timescales $\tau_{\rm{mig}}$ and $\tau_{\rm{dmp}}$ were kept constant, given by equation (7), adopting identical physical parameters of the disk to those employed above. Finally, following previous treatments (Adams et al., 2008; Rein & Papaloizou, 2009; Lecoanet et al., 2009), 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 (8). For completeness, we have also included the leading order corrections due to general relativity (Nobili & Roxburgh, 1986).

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 $\tau_{\rm{mig}}\simeq 1.5\times 10^{3},10^{4}$, and $10^{5}$  years respectively) over a time span of $0.1\,$Myr. The numerical results are shown in Figure 6, and show excellent agreement with the analytic criterion from equation (13). 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 $\phi$[3:2] = $3\lambda_{2}-2\lambda_{1}-\varpi_{1}$ and $\psi$[3:2] = $3\lambda_{2}-2\lambda_{1}-\varpi_{2}$, shown in red and blue in the bottom panel of Figure 6. 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 (13) yields $\delta\chi/\Delta\chi\sim 1$, perfectly exemplifies the transitory regime. As shown in the top panel of Figure 6, 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 $\sim 15,000$ years. Although evolution within this resonance is relatively long-lived, the bottom panel of Figure 6 shows that the corresponding resonant angles $\phi$[7:6] = $7\lambda_{2}-6\lambda_{1}-\varpi_{1}$ (green) and $\psi$[7:6] = $3\lambda_{2}-2\lambda_{1}-\varpi_{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 6, we have confirmed that the presented solutions are indeed representative of the evolutionary outcomes. As a result, we conclude that the analytic expression (13) represents an adequate description of the requirement for resonance disruption, consistent with the numerical experiments.

5. Conclusion

While resonant locking is an expected outcome of migration theory (Cresswell & Nelson, 2008; Ogihara et al., 2015), the current sample of exoplanets shows only a mild tendency for systems to be near mean motion commensurabilities (Winn & Fabrycky, 2015). 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 5.1), followed by a conceptual interpretation of the calculations (Section 5.2), and finally a discussion of the implications (Section 5.3).

5.1. 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 [13]). 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 3).

The analytic estimate (13) 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 $\varepsilon$ (Section 4.1). The results confirm the analytical prediction that given nominal disk parameters, systems with mass ratios smaller than $(m_{1}+m_{2})/M\sim 10^{-5}\sim 3M_{\oplus}/M_{\odot}$ are forced out of resonance by turbulence, whereas systems with larger mass ratios survive (Figure 4). We have also performed full $N$-body simulations of the problem (Section 4.2). 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 6). 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 (Hayashi, 1981), 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\sim 10^{-5}-10^{-4}$ (Figure 1), 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 [13] 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 (Adams et al., 2008; Rein & Papaloizou, 2009; Lecoanet et al., 2009; Ketchum et al., 2011), 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 (Deck et al., 2012; Batygin et al., 2015).

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 (Crida et al., 2006; Duffell & MacFadyen, 2013). However, one complication regarding this issue is that the damping rate of eccentricity is also reduced due to the gap (e.g., Duffell & Chiang 2015). 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.

5.2. 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 (Adams et al., 2008). 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 (Rein & Papaloizou, 2009). 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 (Kley & Nelson, 2012).

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 $\phi$ that obeys d’Alembert rules (e.g., see Chapter 8 of Murray & Dermott 1999). 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 [1]; Wisdom 1986).

In the parameter range relevant to the problem at hand, the behavior of Hamiltonian (1) is well-approximated by that of a simple pendulum (Henrard & Lamaitre, 1983). Specifically, the equilibrium value of $\varepsilon$ dictates the value of the pendulum’s action, $\Phi$, at which zero-amplitude libration of the angle $\phi$ can occur, as well as the location of the separatrix. Correspondingly, oscillation of the angle $\phi$ translates to variations of the action $\Phi$, which is in turn connected to the eccentricities (equations [4]) as well as the semi-major axes, through conservation of the generalized Tisserand parameter (Batygin & Morbidelli, 2013).

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 $\Phi$-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 $\varepsilon$ (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 (Tanaka & Ward, 2004), 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 $\varepsilon$ — ends up being more important for purposes of moving planets out of mean motion resonance (see equation [18]).

5.3. 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 (13) 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 Horn et al. (2012).

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 (Paardekooper, 2014; Bitsch et al., 2015). 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 (Ormel & Okuzumi, 2013), which would weaken the degree of stochastic forcing. Both of these effects can be incorporated into the criterion of equation (13) by lowering the migration factor $f$ and the value of $\alpha$ 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 (Balbus & Hawley, 1991), 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 (Gammie 1996; see also Bai & Stone 2013), and ionization by cosmic rays can be suppressed in the outer disk (Cleeves et al., 2013). Indeed, suppressed levels of ionization are now inferred from ALMA observations of young star/disk systems (Cleeves et al., 2015), 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 (Nelson et al., 2013; Lin & Youdin, 2015).

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 (13), we find that increasing the local surface density by an order of magnitude ($\Sigma=10\langle\Sigma\rangle=170,000\,$g/cm2) while lowering the orbital convergence rate a hundred-fold ($f=0.01$) and retaining $h/r=0.05$, $\langle a\rangle=0.1\,$AU, $\alpha=0.01$ yields $(m_{1}+m_{2})/M\simeq 2\times 10^{-4}\sim 60M_{\oplus}/M_{\odot}$ as the critical mass ratio, thus explaining the full range of values shown in Figure 1. Correspondingly, rough agreement between observations and the stochastic migration scenario is reproduced in the work of Rein (2012), 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 $\sim$few Myr lifetime of the nebula, as long as $\Sigma\gtrsim 0.1\langle\Sigma\rangle=170\,$g/cm2. 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 (Goldreich & Schlichting, 2014; Xu & Lai, 2016) and capture probability suppression (Batygin, 2015). The first mechanism requires that the outer planet is more massive then the inner planet to compromise resonance (Deck & Batygin, 2015). 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_{\rm{disk}}\sim 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 (Mittal & Chiang, 2015).

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.