MEAN MOTION RESONANCES IN EXTRASOLAR PLANETARY SYSTEMS WITH TURBULENCE, INTERACTIONS, AND DAMPING

As a starting point, we use a version of the formalism presented in MD99. To start, this treatment is limited to the circular restricted three-body problem and hence considers the libration angle ϕ for two planets in resonance to have the form

$$\begin{equation} \phi = j_1 \lambda _2 + j_2 \lambda + j_4 \varpi, \end{equation} \tag{ 1 }$$

where λ and λ₂ are the mean longitudes of the two planets and ϖ is the longitude of periapse of the inner planet. In this case, the equation of motion for the resonance angle ϕ reduces to that of a pendulum, i.e.,

$$\begin{equation} {d^2 \phi \over dt^2} + \omega _0^2 \sin \phi = 0 \,. \end{equation} \tag{ 2 }$$

The natural oscillation frequency ω₀ of the pendulum—the libration frequency of the resonance angle—is given by

$$\begin{equation} \omega_0^2 = - 3 j_2^2 {\cal C}_r \Omega e^{|j_4|} \,, \qquad {\rm where} \qquad {\cal C}_r = \mu \Omega \, \alpha f_d(\alpha) \,. \end{equation} \tag{ 3 }$$

Here, Ω and e are the mean motion and eccentricity of the inner planet, and (j₂, j₄) are integers that depend on the type of resonance being considered. The parameters ${\cal C}_r$, e, and Ω are assumed to be constant for purposes of determining the frequency. The mass ratio μ = m₂/M_*, where m₂ is the mass of the outer planet and M_* is the stellar mass. The quantity αf_d(α) results from the expansion of the disturbing function (MD99), and the parameter α is the ratio of the semimajor axes of the two planets, i.e., α ≡ a₁/a₂. Note that this approximation scheme (adapted from MD99) neglects the terms of order ${\cal O}(\mu)$; we consider higher-order terms in Section 4.

Many of the observed (candidate) resonant systems are in or near the 2:1 mean motion resonance; for the sake of definiteness, we focus on that case. This resonance is also generally the strongest. For the 2:1 mean motion resonance, the integer j₂ = −1 and αf_d(α) ≈ −3/4 (MD99). In this case, the natural oscillation frequency ω₀ of the libration angle is given by

$$\begin{equation} \omega _0^2 \approx {9 \over 4} \mu e^{|j_4|} \, \Omega ^2 \,. \end{equation} \tag{ 4 }$$

Typical planet masses in observed extrasolar planetary systems are of order a Jovian mass so that μ ∼ 10⁻³. The eccentricity can vary over a wide range, with e ∼ 0.1 in order of magnitude, but the median eccentricity is closer to e ∼ 0.3. The relevant values of |j₄| = (1,2), and we generally use |j₄| = 2 (see the discussion of Section 2.4). As a result, the ratio of frequencies ω₀/Ω ∼ 10⁻² and hence the period of the libration angle is expected to be ∼100 orbits. With this frequency specified, we define a dimensionless time variable

$$\begin{equation} \tau = \omega _0 t \,, \end{equation} \tag{ 5 }$$

and write the pendulum equation in dimensionless form

$$\begin{equation} {d^2 \phi \over d\tau ^2} + \sin \phi = 0. \end{equation} \tag{ 6 }$$

Note that the energy scale associated with the potential well of the resonance, modeled here as a pendulum, is much smaller than the binding energy of the planets in their orbits. Neglecting dimensionless numbers of order unity, we can write the specific energy E_p of the pendulum as

$$\begin{equation} E_p \sim \omega _0^2 a^2 \sim \mu e^{|j_4|} \, {G M_\ast \over a} \sim \mu e^{|j_4|} \, E_{\rm orb} \sim 10^{-4} \, E_{\rm orb} \,, \end{equation} \tag{ 7 }$$

where E_orb is the energy (per unit mass) of the planet in the gravitational potential well of the star. As result, we estimate that $E_p / E_{\rm orb} = {\cal O} (\mu e^{|j_4|})$, so that the potential well of the resonance is about 10⁴ times shallower than the potential well of the star. As a result, planets can be removed from resonance much more easily than they can be removed from orbit (by turbulent fluctuations—note that planets are often ejected by scattering after orbit crossing).

The net effect of turbulence is to provide stochastic forcing perturbations on the mean motion resonance, which is modeled here as a pendulum. To include this stochastic process, the equation of motion (6) can thus be modified to take the form

$$\begin{equation} {d^2\phi \over d \tau ^2} + \left[ 1 + \eta _k\delta ([\tau ]-\Delta \tau) \right] \sin \phi = 0 \,, \end{equation} \tag{ 8 }$$

where η_k sets the amplitude of the turbulent forcing and Δτ sets the time interval. These quantities are derived in Paper I (including a discussion of where to introduce the turbulent term) and are discussed further below.

Briefly, the timescale Δτ is the time (in dimensionless units) required for the turbulence in the disk to produce an independent realization of the fluctuations. This timescale is typically one or two orbit times at the disk location in question (LSA04; Nelson 2005). Since the relevant part of the disk extends outside the orbit of the outer planet, and since the outer planet has twice the period of the inner planet, the timescale Δτ ≈ 4ω₀(2π/Ω) = 8πω₀/Ω.

Next we need to specify the forcing strength over the time interval Δτ. In the original derivation of the pendulum equation (MD99), one finds

$$\begin{equation} {d^2 \phi \over d t^2} \sim {d \Omega \over dt} \sim - {\Omega \over 3} \left({1 \over J} {dJ \over dt} \right) \,, \end{equation} \tag{ 9 }$$

where the second approximate equality follows from the relationship between the mean motion and the orbital angular momentum. Converting to dimensionless form and integrating over one cycle to produce a discrete quantity, one finds

$$\begin{equation} \eta _k \approx - {1 \over 3} {\Omega \over \omega _0} \left({\Delta J \over J} \right)_k \,. \end{equation} \tag{ 10 }$$

The amplitudes [(ΔJ)/J]_k have been calculated for a variety of cases using MHD simulations (e.g., Nelson & Papaloizou 2003, 2004; LSA04; Nelson 2005). In basic terms, the torque exerted on a planet by the disk will be a fraction of the benchmark scale T_D = 2πGΣrm_P, where Σ is the disk surface density (Johnson et al. 2006). The amplitude for angular momentum variations is thus given by ΔJ = f_TΓ_RT_D(8π/Ω), where the factor 8π/Ω is the time over which one realization of the turbulence acts, f_T ∼ 0.05 is the fraction of the benchmark scale T_D that applies for a disk without a planet, and Γ_R ∼ 0.1 is the reduction factor due to the production of a gap in the disk (and hence loss of disk material). Including all of these factors (see Paper I), the relative fluctuation amplitude is given by [(ΔJ)/J]_rms ∼ 10⁻⁴ under typical conditions. With this typical amplitude for [(ΔJ)/J]_k and Ω/ω₀ ∼ 100 (see above), we expect that η_rms ∼ 0.005.

In general, the fluctuation amplitude will vary from case to case over a range of at least an order of magnitude due to different levels of turbulence and varying surface densities, or equivalently, varying disk masses. The amplitude quoted above is applicable for planets of Jovian mass. Lighter planets produce smaller gaps and have larger values of the factor Γ_R and hence larger forcing amplitudes [(ΔJ)/J]_k. Note that even planets that are too small to clear a gap will still have an effective value of Γ_R < 1 because the wakes produced by the planet compromise the turbulent fluctuations in the vicinity (Nelson & Papaloizou 2004). In addition, different disk systems will retain their disk mass for a range of lifetimes, and this effect leads to another source of system to system variation. Finally, even for a given initial disk mass and disk lifetime, the amount of turbulence that resonant systems experience depends on how late (in the life of the disk) the planets are formed and/or captured into resonance (e.g., Thommes et al. 2008). These variations thus allow for a wide range of possible outcomes.

For completeness, we note that turbulence is not the only possible source of stochastic fluctuations acting on planets in mean motion resonance. For example, a residual disk of planetesimals can provide a granular background gravitational potential and act in a qualitatively similar manner (Murray-Clay & Chiang 2006). The formalism developed here can be used with any source of stochastic fluctuations with a given amplitude η_k and forcing time interval Δτ = ω₀(Δt).

An important technique used to analyze the behavior of stochastic differential equations—including those describing mean motion resonance—is to work in terms of phase-space variables. For the stochastic pendulum considered here, we rewrite the equation of motion in the form

$$\begin{equation} {d \phi \over d \tau } = V \, \qquad {\rm and} \qquad {d V \over d \tau } = - \left[ 1 + \eta _k\delta ([\tau ]-\Delta \tau) \right] \sin \phi \,. \end{equation} \tag{ 11 }$$

The probability distribution function P(τ, ϕ, V) for these phase-space variables obeys the Fokker–Planck equation (e.g., Mallick & Marcq 2004, hereafter MM04) which takes the form

$$\begin{equation} {\partial P \over \partial \tau } + V {\partial P \over \partial \phi } - \sin \phi {\partial P \over \partial V} = {1 \over 2} D \sin ^2 \phi {\partial ^2 P \over \partial V^2 } \,. \end{equation} \tag{ 12 }$$

As outlined in Paper I, the phase-space diffusion constant D is specified by the amplitude η_k of the turbulent fluctuations and the time interval Δτ required for them to attain an independent realization, i.e.,

$$\begin{equation} D = {\big\langle \eta _k^2 \big\rangle \over \Delta \tau } \,. \end{equation} \tag{ 13 }$$

In most cases of interest, the angle ϕ varies more rapidly than the velocity V (see Paper I and MM04). As a result, we can average the Fokker–Planck Equation (12) over the angle ϕ to obtain the corresponding time evolution equation for the ϕ-averaged probability distribution function P(τ, V). Throughout this paper, we use the same symbol (here, P) for the distribution function both before and after an averaging procedure; this choice simplifies the notation, but could result in some ambiguity, although the relevant version of the distribution function should be clear from the context (this issue also arises in Section 4). After averaging, the resulting Fokker–Planck Equation (12) becomes a basic diffusion equation

$$\begin{equation} {\partial P \over \partial \tau } = {D \over 4} {\partial ^2 P \over \partial V^2} \,, \end{equation} \tag{ 14 }$$

which can be solved to obtain the result

$$\begin{equation} P(\tau, V) = {1 \over (\pi D \tau)^{1/2} } \exp \left[ - {V^2 \over D \tau } \right] \,. \end{equation} \tag{ 15 }$$

This solution shows that the energy of the pendulum can grow large at long times. In this limit, the kinetic energy dominates the potential energy term, and we approximate the energy using E ≈ V²/2. The resulting probability distribution function for the energy then can be written in the form

$$\begin{equation} P (E, \tau) = \left({2 \over \pi D \tau } \right)^{1/2} E^{-1/2} \exp \left[ - {2E \over D \tau } \right] \,. \end{equation} \tag{ 16 }$$

For this solution for the distribution function, the expectation value of the energy grows linearly with time in the long time limit, i.e., 〈E〉 = Dτ/4. As the mean energy increases, the probability of the system remaining bound in a low-energy (resonant) state decreases. Here we use E > 1 as the requirement for the libration angle to circulate and hence for the resonance to be compromised (see Paper I for further discussion). The probability P_b of remaining bound is thus given by

$$\begin{eqnarray} { P_{\, {\rm b}} }(\tau) &=& \left({2 \over \pi D \tau } \right)^{1/2} \int _0^1 {dE \over \sqrt{E} } \exp \left[ - {2E \over D \tau } \right]\nonumber\\ &=& {2 \over \sqrt{\pi }} \int _0^{z_0} {\rm e}^{-z^2} dz \, = { {\rm Erf} }(z_0) \,, \end{eqnarray} \tag{ 17 }$$

where z₀ = (2/Dτ)^1/2 and where Erf(z) is the error function (e.g., Abramowitz & Stegun 1970). In the limit of late times, when z ≪ 1, the error function has the asymptotic form ${ {\rm Erf} }(z) \sim 2z/\sqrt{\pi }$, and the probability P_b that the planetary system remains in resonance is given by

$$\begin{equation} { P_{\, {\rm b}} }(\tau) \approx \left({8 \over \pi D \tau } \right)^{1/2} = \left({2 \over \pi \langle E \rangle } \right)^{1/2} \,, \end{equation} \tag{ 18 }$$

where this expression is valid for sufficiently late times when Dτ ≫ 1.

The third and final technique that we use to study this problem is direct numerical integration of the planetary systems. In other words, we integrate the full set of 18 phase-space variables of the corresponding three-body problem (using a B–S integration scheme, e.g., Press et al. 1992). Turbulence is included by applying discrete velocity perturbations at regular time intervals; for the sake of definiteness, the forcing intervals are chosen to be twice the orbital period of the outer planet (four times the period of the inner planet). Both components of velocity in the plane of the orbit are perturbed randomly, but the vertical component of velocity is not changed. The amplitude of the velocity perturbations is then allowed to vary.

For the numerical experiments carried out for this paper, we start a large ensemble of N_ens planetary systems in the 2:1 mean motion resonance, and then integrate forward in time using the velocity perturbations (as described above) to model turbulent fluctuations. Throughout this work, we use the resonance angle defined by Equation (42) for interactive systems such as GJ876 (see Section 4), and its analog in the circular, restricted three-body problem given by Equation (1) with |j₄| = 2. Our numerical experiments show that this choice of resonance angle allows the systems to remain in bound (resonant) states for longer times than for the standard case where |j₄| = 1 (see Paper I and Section 4). Since we are primarily interested in how turbulent fluctuations can compromise mean motion resonance, we want to consider the cases that are most robust (longest lived). Finally, we note that the analysis of this paper can be performed for a host of other forms for the resonance angle.

Over the course of time, systems move both in and out of resonance. To monitor this behavior, we determine the maximum libration amplitude of the resonance angle over a fixed time interval. This timescale is chosen to be somewhat longer than the libration time, which is much shorter than the evolution time, i.e., the timescale over which a substantial fraction of the systems leave resonance. Note that the libration time is typically ∼100 orbits (see Equation (4)), whereas the evolution time is of order 10⁶ orbits. For the sake of definiteness, we consider a planetary system to be in resonance when its maximum libration angle is less than 90°; for larger maximum libration angles, we consider the resonance to be compromised. Note that the boundary at 90° is somewhat arbitrary; fortunately, however, the statistics for the fractions of bound systems are insensitive to this value (see Paper I).

Although the observed set of extrasolar planetary systems shows an astonishing degree of diversity, here we consider only two particular examples that bracket the possibilities. For the first case, we consider analogs of the GJ876 planetary system, which has two planets that are observed to be deep in a 2:1 mean motion resonance (Marcy et al. 2001). This system has a relatively small central star (with mass of only M_* = 0.32 M_☉) and relatively large planets (with masses m = 0.79 and m₂ = 2.52 m_J). The eccentricities are substantial, with e = 0.26 for the inner planet and e₂ = 0.034 for the outer planet. The periods are approximately 30 and 60 days, respectively (see Rivera et al. 2005 for further details). With these parameters, the GJ876 system is the most highly interactive system observed to date. Although comparable solar systems with larger planetary masses and/or higher orbital eccentricities would be even more interactive, they would also most likely be unstable (Laughlin & Chambers 2001). Since the pendulum model for mean motion resonance represents the simplest possible model, and in particular does not allow for interactions between the planets, the GJ876 system (its analogs considered here) is about as far as possible from the idealized one-variable pendulum model. Because of its highly interactive nature, the GJ876 system analogs take a long time to integrate and we consider ensembles of N_ens = 10³.

At the other end of parameter space, we consider a model solar system with relatively little interaction between the two planets. In this case, we take the outer planet to be a Jupiter mass planet (m₂ = 1 m_J) in a 150 day orbit around a 1.0 M_☉ star. The second, inner planet is taken to be a "Super-Earth" with mass m = 0.01m_J. The period of the inner planet is taken to be 75 days, and the system is started in the 2:1 mean motion resonance. The orbital eccentricities are e = 0.1 for the inner planet and e₂ = 0.01 for the outer planet. This latter system is thus well approximated by the circular restricted three-body problem, which is used in the derivation of the pendulum model for mean motion resonance (DM99). We thus expect ensembles of this class of solar system to follow the predictions of the stochastic pendulum model. For these systems, the integrations can be run faster and we consider ensembles of N_ens = 10⁴.

The effects of turbulence on mean motion resonance are illustrated in Figures 1–3. The first two of these figures show the fraction of resonant states as a function of time for an ensemble of solar systems that started out with the configuration of GJ876, and Figure 3 describes results for the less interactive Super-Earth systems.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 shows the various ways to account for the fraction of bound states. The solid curve shows the ratio of the number of planets in resonance to the number of planets in the original sample. However, some planets are lost with time: systems that leave resonance often experience orbit crossings, which can lead to planetary ejection. The dotted curve shows the fraction of systems that have not ejected a planet. Finally, the dashed curve shows the ratio of the number of planets in resonance to the number of planets that remain in orbit. For the remainder of this paper, we use this latter ratio (number of systems in resonance over the number of surviving systems) as the relevant fraction of bound states (this ratio is the most directly observable). Note that a substantial fraction of the systems that leave resonance lose a planet. This ejection takes place because of orbit crossing and subsequent planet–planet scattering. Systems are protected from such behavior while they remain in resonance, but interactive systems (e.g., GJ876) readily lose planets when resonance is compromised.

Figure 2 shows how the time evolution of the bound fraction depends on the level of turbulence. For an ensemble of GJ876 systems, the figure shows the fraction of bound states as a function of time for two levels of turbulent fluctuations, where (Δv)/v = 0.0002 and 0.0004. Note that in both cases, the bound fraction is nearly a straight line on this log–linear plot, so that the time evolution is nearly exponential with a well defined decay rate, or, equivalently, half-life (see Section 4 for further discussion). Note also that the diffusion constant D sets the timescale for systems to leave resonance, and that D depends on the square of the fluctuation amplitude, and hence the decay timescales differ by a factor of 4. The dashed lines shown in Figure 2 show a purely exponential decay, with timescales that differ by exactly a factor of four. The agreement shown in Figure 2 thus confirms our theoretical expectations. The error bars plotted on the dashed lines show the expected amplitudes of root-N fluctuations, which are somewhat smaller than the actual variations found in the numerical results; the presence of planet–planet interactions (see also Section 4) increases the level of these variations.

Figure 3 shows the time evolution of the fraction of systems remaining in resonance for an ensemble of Super-Earth systems. In this set of experiments, two different levels of turbulence are used, so that the amplitudes of the velocity perturbations are (Δv)/v = 0.0001 and 0.0002. This class of systems is close to the idealized, circular restricted three-body problem, and hence the behavior of the resonance is expected to follow the predictions of the stochastic pendulum model. As shown in Figure 3 (by dashed lines), the fraction of bound states decreases as a power law in time, with P_b ∼ t^−1/2 at late times, roughly as predicted in Section 2.3. Given that the velocity perturbations differ by a factor of 2, one expects the coefficients of the power laws to also differ by a factor of 2. However, the numerical results show a somewhat wider gap, namely a factor of ∼2.6. Although the size of the leading coefficient is correct in order of magnitude, the following three ambiguities arise in determining its predicted value from the analytic formulation: the first issue is the mismatch between the simplified treatment of the circular, restricted three-body problem and full (18 variable) integrations; as a result, the libration frequency of the numerical system is not exactly the same as that given by the approximation of Equation (4). The second difference is due to the presence of a partial barrier; after systems leave resonance, they have a somewhat lower probability ρ_ret of re-entering resonance from an excited state because the physical system (with 18 phase-space variables) is more complicated than the one-variable pendulum (see Paper I for further discussion). This effect should not be overly severe for these Super-Earth systems, as they are close to the regime of validity of the pendulum approximations (MD99), but the partial barrier is still present. Yet another relevant process is that planetary ejection can take place. In three-body numerical integrations of Earth-like planets with Jovian companions, when the periastron is the same as in these systems, the expectation value of the ejection time is about 3 × 10⁵ orbits, or 6 × 10⁴ yr (David et al. 2003). In spite of these complications, however, the basic trends shown by the numerical results are in qualitative (and rough quantitative) agreement with the expectations of the simplest pendulum model. Finally, we note that the error bars shown in Figure 3 show the amplitude of root-N fluctuations, which are roughly the same amplitude as the variations seen in the numerical results.

Now we can compare Figure 2 with Figure 3: for the interactive GJ876 systems of Figure 2, the fraction of bound states P_b is presented as a log–linear plot, so that straight lines (as found in the simulations) correspond to exponential decay. For the Super-Earth systems of Figure 3, however, the fraction P_b is presented as a log–log plot, where straight lines (as found in these simulations) correspond to power-law decay. Since exponential decay is much stronger than power-law decay, we expect that highly interactive systems will leave mean motion resonance much more readily than less interactive systems. This effect is accounted for in the analytic model developed in Section 4. In particular, this model predicts power-law decay of the fraction of bound states in the regime of minimal interactions between planets, and exponential decay in the regime of strong interactions.

In this section, we derive the equations of motion to describe the evolution of one particular resonance angle and the corresponding orbital elements. Note that we cannot use the same resonance angle as in Section 2. In that case, the simple pendulum equation resulted from the circular restricted three-body problem (see MD99). In this context, we want to allow for interactions between the planets (so we must move beyond the restricted three-body problem) and hence the outer planet can attain nonzero and time-varying eccentricity (so the circular approximation is no longer applicable). We thus consider the following (nonstandard) resonance angle

$$\begin{equation} \phi = 2\lambda _2 - \lambda + 2\varpi _2 - 2\varpi \,, \end{equation} \tag{ 42 }$$

where λ is the mean longitude and ϖ is the argument of periastron. Note that this resonance angle is the analog (beyond the circular, restricted three-body approximation) of that from Equation (1) for the case |j₄| = 2. Throughout this derivation, the variables associated with the outer planet are denoted with the subscript "2,', whereas the variables associated with the inner planet are left with no subscript.

For our chosen resonance angle, Equation (42), the lowest-order part of the disturbing function for the inner planet is given by

$$\begin{equation} \langle \mathcal {R} \rangle =\frac{Gm_2}{a_2} \big[ \mathcal {R}^{{\rm (sec)}}_D+e^2e_2^2[ f_d(\alpha)+f_e(\alpha) ]\cos \phi \big] \end{equation} \tag{ 43 }$$

and that for the outer planet is

$$\begin{equation} \langle \mathcal {R}_2 \rangle = \frac{Gm}{a} \big[ \alpha \mathcal {R}^{{\rm (sec)}}_D+e^2e_2^2 [ \alpha f_d(\alpha)+f_i(\alpha) ]\cos \phi \big] \end{equation} \tag{ 44 }$$

where e is the eccentricity, m is the mass, (a, a₂) are the semimajor axes, and G is the gravitational constant. The functions f_d, f_e, f_i are defined and discussed in MD99. To lowest order, the direct secular contribution is then given by

$$\begin{equation} \mathcal {R}^{{\rm (sec)}}_D=\big(e^2+e_2^2\big)f_{s,1}(\alpha) + e e_2 f_{s,2}(\alpha)\cos (\varpi _2-\varpi). \end{equation} \tag{ 45 }$$

The corresponding equations of motion for the orbital elements of the inner planet then become

$$\begin{equation} \dot{\Omega } =-\frac{3}{a^2}\frac{\partial \mathcal {R}}{\partial \lambda }= -\frac{3Gm_2}{a^2a_2}e^2e_2^2\left[ f_d(\alpha)+f_e(\alpha)\right]\sin \phi, \end{equation} \tag{ 46 }$$

$$\begin{eqnarray} \dot{e} &=&-\frac{1}{\Omega a^2e}\frac{\partial \mathcal {R}}{\partial \varpi } =-\frac{Gm_2}{\Omega a^2a_2e} [ee_2f_{s,2}(\alpha)\sin (\varpi _2-\varpi) \nonumber\\ &&+\, 2e^2e_2^2\left[f_d(\alpha)+f_e(\alpha)\right]\sin \phi ], \end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} \dot{\varpi } &=&\frac{1}{\Omega a^2e}\frac{\partial \mathcal {R}}{\partial e}= \frac{Gm_2}{\Omega a^2a_2e}[2ef_{s,1}(\alpha) + e_2f_{s,2}(\alpha)\cos (\varpi _2-\varpi)\nonumber\\ &+&2ee_2^2\left[f_d(\alpha)+ f_e(\alpha)\right]\cos \phi] \,. \end{eqnarray} \tag{ 48 }$$

To complete the set, note that $\dot{\epsilon } = e^2 \dot{\varpi }/2$, where is the mean longitude at epoch. However, because $\dot{\epsilon }$ is e² times smaller than $\dot{\varpi }$, we can ignore time derivatives of in favor of time derivatives of ϖ. This assumption is valid within the range of validity for these equations, which are the lowest-order terms in an expansion in the eccentricities (e, e₂). The equations of motion for the orbital elements of the outer planet are similar, i.e.,

$$\begin{equation} \dot{\Omega }_2 = -\frac{3}{a_2^2} \frac{\partial \mathcal {R}_2}{\partial \lambda _2}= \frac{6Gm}{a_2^2a}e^2e_2^2 \left[ \alpha f_d(\alpha)+f_i(\alpha)\right]\sin \phi, \end{equation} \tag{ 49 }$$

$$\begin{eqnarray} \dot{e}_2 &=& \frac{Gm}{\Omega _2 a_2^2ae_2} \big[ee_2 \alpha f_{s,2}(\alpha)\sin (\varpi _2-\varpi)\nonumber\\ &&+ 2e^2e_2^2 \left[\alpha f_d(\alpha)+f_i(\alpha)\right]\sin \phi\big], \end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray} \dot{\varpi }_2 &=& \frac{Gm}{\Omega _2a_2^2ae_2}[2e_2 \alpha f_{s,1}(\alpha) + e \alpha f_{s,2}(\alpha)\cos (\varpi _2-\varpi)\nonumber\\ &&+2e^2e_2 \left[\alpha f_d(\alpha)+f_i(\alpha)\right]\cos \phi ]. \end{eqnarray} \tag{ 51 }$$

Next we rewrite λ as λ = Ωt + , and assume that $\dot{\Omega }t \ll \Omega$, and that $\ddot{\epsilon }\ll \ddot{\varpi }$. The time evolution of the resonance angle is then given by

$$\begin{equation} \ddot{\phi }\approx 2\dot{\Omega _2} - \dot{\Omega }+2\ddot{\varpi _2}-2\ddot{\varpi } \,. \end{equation} \tag{ 52 }$$

Also, one can write GM_* = Ω²a³ = Ω²₂a³₂, where M_* is the mass of the central star. Using this latter result, we define

$$\begin{equation} \mathcal {C} _r \equiv \frac{Gm_2}{\Omega a^2 a_2} \left[f_d(\alpha)+f_e(\alpha)\right]= \frac{m_2}{M_\ast } \Omega \alpha \left[f_d(\alpha)+f_e(\alpha)\right], \end{equation} \tag{ 53 }$$

$$\begin{equation} \mathcal {C} _{s1} \equiv \frac{Gm_2}{\Omega a^2a_2}f_{s,1}(\alpha)= \frac{m_2}{M_\ast } \Omega \alpha f_{s,1}(\alpha), \end{equation} \tag{ 54 }$$

$$\begin{equation} \mathcal {C} _{s2} \equiv \frac{Gm_2}{\Omega a^2a_2} f_{s,2}(\alpha)=\frac{m_2}{M_\ast } \Omega \alpha f_{s,2}(\alpha), \end{equation} \tag{ 55 }$$

and define $\mathcal {C}^{(2)} _r, \mathcal {C}^{(2)} _{s1}, \mathcal {C}^{(2)} _{s2}$ in an analogous way (where the superscript now refers to the second planet). With this choice of notation, the equations of motion take the form

$$\begin{equation} \dot{\Omega } =-3 \mathcal {C} _r \Omega e^2e_2^2\sin \phi, \end{equation} \tag{ 56 }$$

$$\begin{equation} \dot{e} =- \mathcal {C} _{s2}e_2\sin \theta -2 \mathcal {C} _r ee_2^2\sin \phi, \end{equation} \tag{ 57 }$$

$$\begin{equation} \dot{\varpi } =2 \mathcal {C} _{s1}+ \mathcal {C} _{s2}\frac{e_2}{e} \cos \theta + \mathcal {C} _r e_2^2\cos \phi, \end{equation} \tag{ 58 }$$

$$\begin{equation} \dot{\Omega }_2 =6 \mathcal {C}^{(2)} _r \Omega _2 e^2e_2^2\sin \phi, \end{equation} \tag{ 59 }$$

$$\begin{equation} \dot{e}_2 = \mathcal {C}^{(2)} _{s2}e\sin \theta +2 \mathcal {C}^{(2)} _r e^2e_2\sin \phi, \end{equation} \tag{ 60 }$$

$$\begin{equation} \dot{\varpi }_2 =2 \mathcal {C}^{(2)} _{s1}+ \mathcal {C}^{(2)} _{s2} \frac{e}{e_2}\cos \theta + \mathcal {C}^{(2)} _r e^2\cos \phi, \end{equation} \tag{ 61 }$$

where θ ≡ ϖ₂ − ϖ.

After taking a second time derivative of ϖ and ϖ₂, we can derive the equation of motion, which gives $\ddot{\phi }$ in terms of the variables (ϕ, e, e₂, Ω, Ω₂, θ). For simplicity, we also leave the variables ($\dot{e},\dot{e}_2,\dot{\theta }$) in the equation, but these can all be rewritten in terms of the variables listed above. The resulting equation of motion then takes the form

$$\begin{eqnarray} \ddot{\phi } &=& 3e^2e_2^2 \left(4 \mathcal {C}^{(2)} _r \Omega _2 + \mathcal {C} _r \Omega \right) \sin \phi\nonumber\\ &&+ 2\left[ \mathcal {C}^{(2)} _{s2}\frac{\dot{e}}{e_2}\left(1 + \frac{ \mathcal {C} _{s2}e_2^2}{ \mathcal {C}^{(2)} _{s2}e^2}\right) + \mathcal {C} _{s2}\frac{\dot{e}_2}{e} \left(1 + \frac{ \mathcal {C}^{(2)} _{s2}e^2}{ \mathcal {C} _{s2}e_2^2}\right)\right]\cos \theta\nonumber\\ &&+ 2 \left(\mathcal {C}^{(2)} _{s2}\frac{e}{e_2}- \mathcal {C} _{s2}\frac{e_2}{e}\right) \dot{\theta } \sin \theta + 2\left(\mathcal {C}^{(2)} _re\dot{e}- \mathcal {C} _re_2\dot{e}_2\right) \cos \phi\nonumber\\ &&+ 2\big(\mathcal {C}^{(2)} _re^2- \mathcal {C} _re_2^2\big)\dot{\phi }\sin {\phi }. \end{eqnarray} \tag{ 62 }$$

Consider the ordering of the terms in Equation (62). All of the time derivative terms (i.e., those including $\dot{e}$ or $\dot{\phi }$), have an extra factor of $\mathcal {C}$, which introduces another factor of μ = m/M_*. All of the terms in the equation are thus second order in m/M_*, except those coming from the $\dot{\Omega }$ and $\dot{\Omega }_2$ terms, which are only first order. As a result, to leading order, we recover the usual pendulum equation of motion to describe a mean motion resonance, i.e.,

$$\begin{equation} \ddot{\phi }= 3 e^2 e_2^2 \left(4 \mathcal {C}^{(2)} _r \Omega _2 + \mathcal {C} _r \Omega \right) \sin \phi, \end{equation} \tag{ 63 }$$

where the effective frequency of the system is a linear combination of the frequencies of the "oscillators" for each planet separately. Note that the frequency of this pendulum contains extra factors of eccentricity compared with the simplified case considered in Section 2; this difference arises because the two pendulum equations, Equations (2) and (63), correspond to different resonance angles.

Now that we have derived the "classic" equation of motion for ϕ, without fluctuations, we must include the effects of turbulence. In this treatment, we incorporate turbulent forcing as a set of discrete impulses in the equation of motion for the eccentricity, i.e.,

$$\begin{eqnarray} \dot{e}&=& - \mathcal {C} _{s2}e_2\sin \theta -2 \mathcal {C} _r ee_2^2\sin \phi + \xi\quad {\rm where}\nonumber\\ \xi &=& \eta _k \delta \left([t] - \Delta t \right). \end{eqnarray} \tag{ 64 }$$

As shown above, the equations of motion for all of the variables include eccentricity, and hence this ansatz for including turbulence is convenient. These impulses acting on the eccentricity then produce corresponding impulses in the equation of motion for the resonance angle ϕ, and in the equations for ϖ and ϖ₂. Since the orbital elements of the external planet should not be directly perturbed by an impulse acting on the inner planet, the validity of including impulses in the latter equation is somewhat subtle. However, the effect shows up in the second derivative, which acts like a force, and which in turn changes when the inner planet experiences a perturbation (here due to turbulence). Finally, we note that the distribution of impulses produced by Equation (64) is more complicated than the uniform distribution of velocity perturbations used in our numerical simulations. Because of this complication, the amplitude of the eccentricity perturbations will be different from the effective amplitude of the perturbations acting on the other variables (e.g., the velocity $V = {\dot{\phi }}$).

We can now write out the Fokker–Planck equation, which includes six variables at this stage of derivation. In order to obtain tractable results, we average over most of the oscillating variables. To start, we note that the equation simplifies considerably after we average over the angle ϕ itself. Although this approximation is often used (e.g., MM04), its implementation in this case is somewhat complicated. Here, the frequency of the largest oscillations of ϕ (those lowest order in $\mathcal {C}$) is proportional to $\sqrt{ \mathcal {C} }$, whereas the primary frequencies of the other quantities are proportional to $\mathcal {C}$; note that $\mathcal {C} \ll 1$, so that ϕ oscillates with a higher frequency. Further, the terms for which this frequency difference does not hold are small and can be neglected. As a result, ϕ has effectively a higher frequency and can be averaged out (see Paper I and MM04 for further justification). After averaging, the Fokker–Planck equation becomes

$$$ \frac{\partial P}{\partial t} = \mathcal {C} _{s2}e_2\sin \theta \frac{\partial P}{\partial e}- \mathcal {C}^{(2)} _{s2}e\sin \theta \frac{\partial P}{\partial e_2} - 2 \left(\mathcal {C}^{(2)} _{s1} - \mathcal {C} _{s1} \right) {\partial P \over \partial \theta } $$$

$$$ +\left(\mathcal {C}^{(2)} _{s2}\frac{e}{e_2}P- \mathcal {C} _{s2} \frac{e_2}{e}P\right)\sin \theta - \left(\mathcal {C}^{(2)} _{s2} \frac{e}{e_2} - \mathcal {C} _{s2}\frac{e_2}{e} \right) {\partial P \over \partial \theta } \cos \theta $$$

$$\begin{eqnarray} &-& 2\left(\mathcal {C}^{(2)} _{s2}\frac{e}{e_2}- \mathcal {C} _{s2} \frac{e_2}{e}\right)\dot{\theta }\sin (\varpi _2-\varpi) \frac{\partial P}{\partial V} + \frac{{ \mathcal {D} }_e}{2}\frac{\partial ^2 P}{\partial e^2} \nonumber\\ &+& 2 { \mathcal {D} }_{ev} \mathcal {C}^{(2)} _{s2}\frac{1}{e_2}\left(1+\frac{ \mathcal {C} _{s2}e_2^2}{ \mathcal {C}^{(2)} _{s2}e^2}\right) \cos \theta \frac{\partial ^2 P}{\partial V \partial e}\nonumber \end{eqnarray}$$

$$\begin{equation} + 2 { \mathcal {D} }_v \left\lbrace \left[ \mathcal {C}^{(2)} _{s2}\frac{1}{e_2} \left(1+\frac{ \mathcal {C} _{s2}e_2^2}{ \mathcal {C}^{(2)} _{s2}e^2}\right) \cos \theta \right]^2 + \frac{1}{2} \left(\mathcal {C}^{(2)} _r e \right)^2\right\rbrace \frac{\partial ^2 P}{\partial V^2} \,. \end{equation}$$

Because of the difference in perturbation amplitudes for the eccentricity e and for the velocity $\dot{\phi }$ = V, the diffusion constants ${ \mathcal {D} }_e$, ${ \mathcal {D} }_{ev}$, and ${ \mathcal {D} }_v$ are not equal (in general).

Next we consider how the angle θ evolves with time. We start with the equation of motion for θ in the form

$$\begin{eqnarray} \ddot{\theta }&=& \left[\frac{\dot{e}}{e_2} \left(\mathcal {C}^{(2)} _{s2}+ \mathcal {C} _{s2} \frac{e_2^2}{e^2}\right)-\frac{\dot{e}_2}{e} \left(\mathcal {C} _{s2}+ \mathcal {C}^{(2)} _{s2} \frac{e^2}{e_2^2}\right)\right]\cos \theta\nonumber\\ &&+ \left(\mathcal {C}^{(2)} _{s2}\frac{e}{e_2}- \mathcal {C} _{s2}\frac{e_2}{e}\right) \dot{\theta }\sin \theta\nonumber\\ &&+ 2\big(\mathcal {C}^{(2)} _re\dot{e}+ \mathcal {C} _re_2\dot{e}_2\big) \cos \phi\nonumber\\ && +2\big(\mathcal {C}^{(2)} _re^2- \mathcal {C} _re_2^2\big)\dot{\phi }\sin \phi. \end{eqnarray}$$

As discussed above, in this approximation we average over the (rapidly oscillating) variable ϕ. This averaging reduces the equation of motion for θ to the form

$$\begin{eqnarray} \ddot{\theta }&=& \left[\frac{\dot{e}}{e_2} \left(\mathcal {C}^{(2)} _{s2}+ \mathcal {C} _{s2}\frac{e_2^2}{e^2}\right)- \frac{\dot{e}_2}{e}\left(\mathcal {C} _{s2}+ \mathcal {C}^{(2)} _{s2}\frac{e^2}{e_2^2}\right)\right]\cos \theta\nonumber\\ &&+\left(\mathcal {C}^{(2)} _{s2}\frac{e}{e_2}- \mathcal {C} _{s2}\frac{e_2}{e}\right)\dot{\theta }\sin \theta, \end{eqnarray} \tag{ 67 }$$

where now $\dot{e},\dot{e}_2,\dot{\theta }$ are ϕ-averaged, and thus obey the reduced equations of motion

$$\begin{equation} \dot{e} =- \mathcal {C} _{s2}e_2\sin \theta, \end{equation} \tag{ 68 }$$

$$\begin{equation} \dot{e}_2 = \mathcal {C}^{(2)} _{s2}e\sin \theta, \end{equation} \tag{ 69 }$$

$$\begin{equation} \dot{\theta } = 2 \mathcal {C}^{(2)} _{s1}-2 \mathcal {C} _{s1}+ \left(\mathcal {C}^{(2)} _{s2}\frac{e}{e_2}- \mathcal {C} _{s2}\frac{e_2}{e}\right)\cos \theta. \end{equation} \tag{ 70 }$$

When the masses of the planets are comparable, then $\mathcal {C} \sim \mathcal {C}^{(2)}$, and the first two terms on the right-hand side of Equation (70) nearly cancel. If, in addition, the planets have the appropriate values of eccentricity (e, e₂), then the final terms can also nearly cancel out. The small remaining part of the first two terms can then balance that of the second two terms, $\dot{\theta } \approx 0$, and the system can be trapped in a "θ-resonance" or apsidal resonance (for further discussion, e.g., see Chiang & Murray 2002; Murray-Clay & Chiang 2006). Numerical experiments show that the Super-Earth systems considered in this paper are not in apsidal resonance, whereas the analogs of the GJ876 system are. Further, if the planetary masses in the GJ876 systems are changed slightly, then the value of e/e₂ adjusts so that $\dot{\theta }\approx 0$.

To see how this latter adjustment works, consider the time derivative of e/e₂, i.e.,

$$\begin{equation} \frac{d}{dt}\left(\frac{e}{e_2}\right)=\frac{\dot{e}}{e_2}- \frac{e\dot{e_2}}{e_2^2}=- \mathcal {C} _{s2}\sin \theta - \mathcal {C}^{(2)} _{s2}\frac{e^2}{e_2^2}\sin \theta. \end{equation} \tag{ 71 }$$

Keep in mind that both $\mathcal {C} _{s2}$ and $\mathcal {C}^{(2)} _{s2}$ are negative. Now, suppose that θ is increasing with sin θ>0. Then e/e₂ would increase. This behavior makes the positive term $- \mathcal {C} _{s2} (e_2/e) \cos \theta$ in Equation (70) decrease in magnitude, and the negative term, $\mathcal {C}^{(2)} _{s2} (e/e_2) \cos \theta$, increase in magnitude, which would cause θ to begin to decrease. The same argument can be used to show that e/e₂ would decrease if θ became negative, which would result in θ becoming more positive.

When the system is not in apsidal resonance, the variables (θ, e, e₂) all oscillate at a frequency given by a characteristic value of $\dot{\theta }$. As shown by Equation (70), this characteristic value of $\dot{\theta }$ is of the order $\mathcal {C}$. Likewise, when the system is in apsidal resonance, so that $\dot{\theta }\approx 0$, the equation of motion for θ reduces to the approximate form

$$\begin{equation} \ddot{\theta } \approx -\left[ \mathcal {C} _{s2} \left(\mathcal {C}^{(2)} _{s2}+ \mathcal {C} _{s2}\frac{e_2^2}{e^2}\right)+ \mathcal {C}^{(2)} _{s2}\left(\mathcal {C} _{s2}+ \mathcal {C}^{(2)} _{s2} \frac{e^2}{e_2^2}\right)\right]\theta. \end{equation} \tag{ 72 }$$

In this case, the variables θ, e, and e₂ still oscillate at a characteristic frequency given by $\mathcal {C}$. This result justifies averaging over the angle ϕ, which has a characteristic oscillation frequency of $\sqrt{ \mathcal {C} }$ (keep in mind that $\sqrt{ \mathcal {C} } \gg \mathcal {C}$ because $\mathcal {C} \ll 1$). However, we cannot average the Fokker–Planck equation over θ (i.e., over ϖ and ϖ₂) in either case, because e and e₂ both vary on the same timescale and because the amplitude of this variation is substantial.

Next we consider the effect of turbulence on θ. As above, we will add a turbulent forcing term to $\dot{e}$,

$$\begin{equation} \dot{e}= - \mathcal {C} _{s2}e_2\sin \theta + \xi, \, \end{equation} \tag{ 73 }$$

where ξ has the same meaning as before. We can now write a system of first-order differential equations for θ, e, e₂,

$$\begin{equation} \dot{\theta } \equiv \omega, \end{equation} \tag{ 74 }$$

$$\begin{equation} \dot{e}=- \mathcal {C} _{s2}e_2\sin \theta +\xi, \end{equation} \tag{ 75 }$$

$$\begin{equation} \dot{e}_2= \mathcal {C}^{(2)} _{s2}e\sin \theta, \end{equation} \tag{ 76 }$$

$$\begin{eqnarray} \dot{\omega }&=& {-} \left[ \mathcal {C} _{s2} \left(\mathcal {C}^{(2)} _{s2}+ \mathcal {C} _{s2} \frac{e_2^2}{e^2}\right)+ \mathcal {C}^{(2)} _{s2} \left(\mathcal {C} _{s2}+ \mathcal {C}^{(2)} _{s2} \frac{e^2}{e_2^2}\right)\right]\cos \theta \sin \theta\nonumber\\ &&+ \left(\mathcal {C}^{(2)} _{s2}\frac{e}{e_2}- \mathcal {C} _{s2}\frac{e_2}{e}\right) \left[2 \mathcal {C}^{(2)} _{s1}-2 \mathcal {C} _{s1}\right.\nonumber\\ &&\left.+ \left(\mathcal {C}^{(2)} _{s2}\frac{e}{e_2}- \mathcal {C} _{s2} \frac{e_2}{e}\right)\cos \theta \right]\nonumber\\ &&\times \sin \theta + \frac{\xi }{e_2}\left(\mathcal {C}^{(2)} _{s2}+ \mathcal {C} _{s2}\frac{e_2^2}{e^2}\right)\cos \theta, \end{eqnarray} \tag{ 77 }$$

and find their equivalent Fokker–Planck equation:

$$\begin{eqnarray} \frac{\partial Q}{\partial t} = &-& \omega \frac{\partial Q}{\partial \theta }+ \mathcal {C} _{s2}e_2\sin \theta \frac{\partial Q}{\partial e}- \mathcal {C}^{(2)} _{s2}e \sin \theta \frac{\partial Q}{\partial e_2}+ \dot{\omega }\frac{\partial Q}{\partial \omega }\nonumber\\ &+& \frac{{ \mathcal {D} }_e}{2}\frac{\partial ^2 Q}{\partial e^2}+\frac{{ \mathcal {D} }_{e\omega }}{e_2}\left(\mathcal {C}^{(2)} _{s2}+ \mathcal {C} _{s2} \frac{e_2^2}{e^2}\right)\cos \theta \frac{\partial ^2 Q}{\partial e\partial \omega }+ \frac{{ \mathcal {D} }_\omega }{2e_2^2}\nonumber\\ &\times& \left(\mathcal {C}^{(2)} _{s2}+ \mathcal {C} _{s2}\frac{e_2^2}{e^2}\right)^2\cos ^2(\theta) \frac{\partial ^2 Q}{\partial \omega ^2} \,. \end{eqnarray}$$

Here we use Q to denote the distribution function for the apsidal resonance angle θ, which is not to be confused with the distribution function P for the original resonance angle ϕ. Equation (78) has almost the same diffusion terms (those containing the diffusion constants ${ \mathcal {D} }_e$, ${ \mathcal {D} }_{e\omega }$, and ${ \mathcal {D} }_\omega$) as the Fokker–Planck equation, Equation (65) for P (only the factors of 2 are different). This correspondence suggests that for systems starting in apsidal resonance, the bound fraction for apsidal resonance should show qualitatively the same time dependence as that of the ϕ-resonance. For example, we observe exponential decays in the bound fraction for both angles in the GJ876 system (see also Paper I). Furthermore, the effective diffusion constant for apsidal resonance is approximately four times larger than the diffusion constant for removing systems from ϕ-resonance. Systems thus tend to leave apsidal resonance four times faster, which is a trend seen in our numerical simulations.

To make further progress, we need to average over θ. This must be done carefully because all the terms in the problem vary on the same timescale as θ. We assume that the system is not in θ-resonance, and that θ = ωt, for some ω. As a result, both e and e₂ can vary with substantial amplitude about their mean values (denoted here as 〈e〉 and 〈e〉₂) in the following way:

$$\begin{equation} e = { \langle e \rangle }+ \frac{ \mathcal {C} _{s2}{ \langle e \rangle }_2}{\omega }\cos \theta \equiv { \langle e \rangle }+ \delta e \cos \theta \,, \end{equation} \tag{ 79 }$$

$$\begin{equation} e_2 = { \langle e \rangle }_2 - \frac{ \mathcal {C}^{(2)} _{s2}{ \langle e \rangle }}{\omega }\cos \theta \equiv { \langle e \rangle }_2 - \delta e_2 \cos \theta \,, \end{equation} \tag{ 80 }$$

where the second two equalities define δe and δe₂. This ansatz assumes that these perturbations of the eccentricity are smaller than the eccentricities e and e₂. However, these perturbations are nonetheless significant and, as derived below, provide the effects of interactions in this formulation.

From the original equations of motion, we see that $\omega \approx 2 \mathcal {C}^{(2)} _{s1} - 2 \mathcal {C} _{s1}$. We can now average over θ in Equation (65), and derive the simplified form for the Fokker–Planck equation:

$$\begin{eqnarray} && \frac{\partial P}{\partial t} = \frac{{ \mathcal {D} }_e}{2}\frac{\partial ^2 P}{\partial e^2}+ 2{ \mathcal {D} }_{ev}\nonumber\\ &&\times\left[\frac{ \mathcal {C}^{(2)} _{s2}}{{ \langle e \rangle }_2} \left(-\frac{\delta e_2}{2{ \langle e \rangle }_2}\right)+ \frac{ \mathcal {C} _{s2}{ \langle e \rangle }_2}{{ \langle e \rangle }^2} \left(\frac{\delta e_2}{2{ \langle e \rangle }_2}- \frac{\delta {e}}{{ \langle e \rangle }}\right)\right] \frac{\partial ^2 P}{\partial V \partial e}\nonumber \end{eqnarray}$$

$$\begin{equation} + 2{ \mathcal {D} }_v \left\lbrace \frac{1}{2}\left[ \frac{ \mathcal {C}^{(2)} _{s2}}{{ \langle e \rangle }_2} \left(1+\frac{ \mathcal {C} _{s2}{ \langle e \rangle }_2^2}{ \mathcal {C}^{(2)} _{s2}{ \langle e \rangle }^2} \right)\right]^2 + \frac{1}{2}\left(\mathcal {C}^{(2)} _r{ \langle e \rangle }\right)^2\right\rbrace \frac{\partial ^2 P}{\partial V^2} \,, \end{equation} \tag{ 81 }$$

where only the first-order terms in δe and δe₂ have been retained. In general, we can expand these formulae in a power series in δe/〈e〉 and keep higher-order terms. For simplicity, we use only the leading-order terms here. When e₂ ≪ e, we can use the result

$$\begin{equation} 1 + \frac{ \mathcal {C} _{s2}e_2^2}{ \mathcal {C}^{(2)} _{s2}e^2} \approx 1, \end{equation} \tag{ 82 }$$

and we can evaluate the integral (over θ) exactly by contour integration. In this case, we find

$$\begin{eqnarray} && \frac{\partial P}{\partial t} = \frac{{ \mathcal {D} }_e}{2}\frac{\partial ^2 P}{\partial e^2} + 2{ \mathcal {D} }_{ev} \mathcal {C}^{(2)} _{s2} \left\lbrace 1 - \left[ 1 - \left(\delta e_2 / { \langle e \rangle }_2 \right) \right]^{-1/2} \right\rbrace\nonumber\\ &&\times \frac{\partial ^2 P}{\partial V \partial e} + 2 { \mathcal {D} }_v \left(\frac{ \mathcal {C}^{(2)} _{s2}}{2{ \langle e \rangle }_2}\right)^2 {\partial ^2 P \over \partial V^2} \,. \end{eqnarray} \tag{ 83 }$$

As a consistency check, note that if we reduce Equation (81) using the approximation of Equation (82) and then evaluate Equation (83) in the limit δe₂/〈e〉₂ ≪ 1, the resulting expressions agree.

At this point, it is advantageous to average over the eccentricity e. In the Fokker–Planck equation, the second derivative terms with respect to e suggest that the eccentricity is undergoing an ordinary diffusion process, in which case its expectation value would grow like $\sqrt{t}$. In the true physical system, however, the eccentricity is bounded both from below (e ⩾ 0) and from above (e ⩽ 1). The effective upper bound is even smaller because planets that reach a sufficiently large eccentricity often experience orbit crossings and eventual ejection. For sufficiently long timescales, the eccentricity is expected to approach a nearly stationary distribution, and it becomes a reasonable approximation to average over e. We thus define the following quantities:

$$\begin{equation} A \equiv \int \frac{\partial ^2 P}{\partial e^2}\rho (e) de, \end{equation} \tag{ 84 }$$

$$\begin{equation} B \frac{\partial P}{\partial V} \equiv \int \frac{\partial ^2 P}{\partial V\partial e}\rho (e) de, \end{equation} \tag{ 85 }$$

where ρ(e) is the distribution of e. This approximation is straightforward when the original distribution P(t, e, V) is a separable function; when this is not the case, the meaning of the approximation is more complicated. In either case, after averaging, the Fokker–Planck equation becomes

$$\begin{eqnarray} &&\frac{\partial P}{\partial t} = \frac{{ \mathcal {D} }_e}{2}A + 2{ \mathcal {D} }_{ev}\nonumber\\ &&\times\,\left[\frac{ \mathcal {C}^{(2)} _{s2}}{{ \langle e \rangle }_2} \left(-\frac{\delta e_2}{2{ \langle e \rangle }_2}\right)+\frac{ \mathcal {C} _{s2}{ \langle e \rangle }_2}{{ \langle e \rangle }^2} \left(\frac{\delta e_2}{2{ \langle e \rangle }_2}-\frac{\delta {e}}{{ \langle e \rangle }}\right)\right] B\frac{\partial P}{\partial V}\nonumber\\ &&+\, 2 { \mathcal {D} }_v \left\lbrace \frac{1}{2}\left[ \frac{ \mathcal {C}^{(2)} _{s2}}{{ \langle e \rangle }_2} \left(1+\frac{ \mathcal {C} _{s2}{ \langle e \rangle }_2^2}{ \mathcal {C}^{(2)} _{s2}{ \langle e \rangle }^2}\right)\right]^2 + \frac{1}{2}\left(\mathcal {C}^{(2)} _r{ \langle e \rangle }\right)^2\right\rbrace \frac{\partial ^2 P}{\partial V^2} \,,\nonumber\\ \end{eqnarray} \tag{ 86 }$$

which is now a one-dimensional diffusion equation that can be solved analytically using standard methods.

Before proceeding further, we simplify notation by defining constants (X, Y, Z) so that the averaged version of the Fokker–Planck equation (derived above) can be written in the form

$$\begin{equation} \frac{\partial P}{\partial t} = X + Y \frac{\partial P}{\partial V} + Z\frac{\partial ^2 P}{\partial V^2} \,. \end{equation} \tag{ 87 }$$

Using standard methods, we can write the corresponding solution in the form

$$\begin{equation} P (t, V) = {1 \over 2\sqrt{\pi t Z}} \exp \left[-\frac{\left(V-Yt\right)^2}{4Zt}+Xt\right] \,. \end{equation} \tag{ 88 }$$

This solution corresponds to a distribution which is diffusing with diffusion constant Z, drifting like Yt, and decaying on the timescale 1/|X| (note that X is expected to be negative to leading order). In practice, however, we find that the quantity X is small; our numerical simulations indicate that ${ \mathcal {D} }_v / { \mathcal {D} }_e \sim 10^4$ in typical cases. This ratio of diffusion constants is the same as the ratio of the potential energy of the resonance to the potential energy of the planetary orbit (Equation (7)). The diffusion constant ${ \mathcal {D} }_e$ measures the efficacy of turbulence in changing the eccentricity, i.e., changing the orbit; the constant ${ \mathcal {D} }_v$ determines how easily turbulence can change the speed $V = {\dot{\phi }}$, which is easier to change by a factor of ∼10⁴. Because ${ \mathcal {D} }_e$ is small, the decay term in Equation (88) can be ignored on the intermediate timescales relevant to this discussion, and we set X = 0 from this point onward.

Given the solution for P(t,V), we can integrate over the bound states to find the fraction of systems in resonance as a function of time, i.e.,

$$\begin{eqnarray} { P_{\, {\rm b}} }(t) &=& \int _{-k}^k P(t,V) \, dV \, = {1 \over 2}\nonumber\\ &&\times \left[{ {\rm Erf} }\left(\frac{k+Yt}{2\sqrt{Zt}}\right)- { {\rm Erf} }\left(\frac{-k+Yt}{2\sqrt{Zt}}\right)\right], \end{eqnarray} \tag{ 89 }$$

where k is the value of V for which the oscillator begins circulating. When k ≫ Yt, the distribution function is primarily determined by diffusion, and the bound fraction takes the approximate form

$$\begin{equation} { P_{\, {\rm b}} }(t) = { {\rm Erf} }\left(\frac{k}{2\sqrt{Zt}}\right) \,, \end{equation} \tag{ 90 }$$

as in the stochastic pendulum problem. However, when Yt ≫ k, at late times we can approximate the error function with the asymptotic formula

$$\begin{equation} { {\rm Erf} }(x) \approx 1 - \frac{e^{-x^2}}{\sqrt{\pi }x} \,, \end{equation} \tag{ 91 }$$

which is valid in the limit x ≫ 1 (Abramowitz & Stegun 1970). In this limit, the fraction of bound states reduces to the form

$$\begin{eqnarray} { P_{\, {\rm b}} }(t) &\approx& \left({Z t \over \pi } \right)^{1/2} \, \exp \left(-\frac{Y^2}{4Z}t\right) \, \left[ {1 \over Yt - k} - {1 \over Yt + k} \right]\nonumber\\ &\approx& \left({Z \over \pi } \right)^{1/2} \, {2 k \over Y^2 t^{3/2}} \, \exp \left(-\frac{Y^2}{4Z}t\right) \,.\nonumber\\ \end{eqnarray} \tag{ 92 }$$

The fraction of bound states thus displays exponential decay in this limit.

In the limit where e₂ ≪ e, which holds for the GJ876 system, the constants Y and Z reduce to the forms

$$\begin{equation} Y =-{ \mathcal {D} }_{ev} \, \frac{ \mathcal {C}^{(2)} _{s2}}{{ \langle e \rangle }_2} \left(\frac{\delta e_2}{{ \langle e \rangle }_2}\right)B \qquad {\rm and} \qquad Z ={ \mathcal {D} }_v \left({ \mathcal {C}^{(2)} _{s2}\over { \langle e \rangle }_2} \right)^2 \,, \end{equation} \tag{ 93 }$$

and hence the long-term time evolution of the bound fraction is given by

$$\begin{equation} { P_{\, {\rm b}} }\sim \exp \left[ - {B^2 { \mathcal {D} }_{ev}^2 \over 4 { \mathcal {D} }_v} \left(\frac{\delta e_2}{{ \langle e \rangle }_2}\right)^2t\right] \,. \end{equation} \tag{ 94 }$$

From the definition of δe₂, we see that δe₂/〈e〉₂ will be small when $\mathcal {C}^{(2)} _{s2}/ \mathcal {C} _{s1}$ is small, i.e., when the system is not highly interactive.

We can now summarize the implications of this solution: if the system is not sufficiently interactive, then the eccentricities do not vary appreciably; in mathematical terms, the factors δe₂ and δe will be small, the parameter Y will be small, and the system will exhibit a power-law decrease in the fraction of bound states as implied by Equation (90). In the opposite limit of highly interactive systems, δe₂ and δe are nontrivial, Y ≠ 0, and the solution displays the exponential decay indicated by Equation (94) in the limit of long times. In this interactive case, since the eccentricities are varying and angular momentum is conserved, the semimajor axes of the planets will vary. As a result, the mean motions (Ω, Ω₂) will also vary and the location of the resonance will move around. This movement of the resonance thus corresponds to the "drift" in V obtained in the solution to the Fokker–Planck equation in the above analysis. When, in addition, the planets have large enough masses so that the star itself moves substantially (as in the case of GJ876), the location of the resonance moves even more.

Next we want to compare the timescales for diffusion and drift. In the absence of the drifting effect, Equation (90) shows that the bound fraction evolves on a timescale t_diffuse = k²/(4Z). Furthermore, since P_b ∼ t^−1/2, in order for turbulence to have an effect, the evolution time must take place over many diffusion times; for example, P_b ≈ 0.1 requires that t ∼ 100 t_diffuse. When the drift term becomes significant, Equation (92) shows that the bound fraction decays exponentially with a timescale t_exp = 4Z/Y². To simplify this discussion, let e and μ represent the eccentricities and mass ratios of both planets, and let δ = δe_j/〈e〉_j be the amplitude for eccentricity variations for both planets. The ratio of the two timescales is then given by

$$\begin{equation} {{ t_{\rm diffuse} }\over { t_{\rm exp} }} = \left({k Y \over 4 Z} \right)^2 \approx \, {\widetilde{B}}^2 \, \delta ^2 \, {e^6 \over \mu } \, \left({{ \mathcal {D} }_{ev} \over { \mathcal {D} }_v} \right)^2 \,, \end{equation} \tag{ 95 }$$

where we have absorbed numerical factors of order unity into the constant $\widetilde{B}$. For relatively large eccentricities, say e ∼ 0.3 near the peak of the observed distribution for extrasolar planets, the factor (e⁶/μ) is of order unity (recall also that e = 0.26 for GJ876). If the diffusion constants are comparable, then Equation (95) implies t_diffuse/t_exp ∼ δ². Thus, in order for exponential decay to be realized, the parameter δ must be relatively large, i.e., the changes in eccentricity must be comparable to the mean values 〈e〉. This condition is met when the two planets are relatively large and have comparable mass, so that $\mathcal {C} _{s2} \approx \mathcal {C}^{(2)} _{s2}$ and δ ∼ 1/2 (see Equations (79) and (80)). Because exponential decay is much faster than power-law decay, this requirement on δ is less severe than it might seem: suppose, for example, turbulence acts for 100 diffusion times, so that the bound fraction P_b ≈ 0.1 in the absence of interactions. The decay term is then roughly given by exp [ − 100t_diffuse/t_exp] ∼ exp [ − 100δ²]. Thus, the fraction of bound states can be affected for more moderately interacting systems, e.g., with values of δ ⩾ 0.1. Keep in mind, however, that the ratio of timescales depends on the other variables ($e, \mu, { \mathcal {D} }_v, { \mathcal {D} }_{ev}$) and will thus vary greatly from system to system.

To test the predictions of this model, we compare the fraction of bound states (resonant systems) from Equation (92) with the fraction derived from an ensemble of numerical integrations of the GJ876 system. The result is shown in Figure 6, which shows the bound fraction—the ratio of the number of systems in resonance to the number of systems left in orbit—as a function of time. Both the analytic model and the numerical integrations show a nearly exponential decline in the number of bound states (Figure 6 is presented as a log–linear plot so that exponential decay corresponds to a straight line). The analytic model is given in terms of error functions with dimensionless arguments of the form $x = (k \pm Y t) / 2 \sqrt{Zt}$. As a result, any fit to the numerical result corresponds to a one-parameter family of values of the parameters (k, Y, Z); the values will also depend on the choice of units. For the sake of definiteness, we have chosen k = 1. Since k is approximately the libration frequency of the resonance, the other two variables, here the diffusion constants Y and Z, are then given in terms of libration frequencies. In these units (and for time t measured in years), we find that the values Y ≈ 5.3 × 10⁻⁴ and Z ≈ 3.1 × 10⁻⁴ provide a good fit (as shown in Figure 6).

Zoom In Zoom Out Reset image size

These results are sensible, at least in order of magnitude: from the original diffusion equation, the diffusion constant Z ∼ ω²₀/t, where ω₀ is the libration frequency and t is the diffusion time. The original stochastic differential equation also defines the diffusion time, i.e., the time required for random perturbations of size η_k to change the libration speed by a total displacement of order ω₀. This timescale is given by t ∼ 9Pω²₀Ω⁻²[(Δv)/v]², where P is the period of the inner planet (in years). Putting these two results together, and inserting numerical values, we find that the diffusion constant is given by Z ∼ 4(9P)⁻¹ × 10⁻⁴ ∼ 5 × 10⁻⁴ yr⁻¹. Note that this value must be adjusted to units in which k = 1 to compare with the fitting values of our analytic theory to the numerical data; here, the libration period of GJ876 is about 9 years, so that k = 2π/P ∼ 0.7, which is close to unity, so that order-of-magnitude agreement obtains.

One key result emerging from this analysis is that turbulence leads to a power-law decrease in the fraction of resonant bound states in the regime of little interaction between planets, whereas highly interactive systems display exponential decay. As a test of these ideas, we have run the following numerical experiment: we begin with an ensemble of GJ876 system analogs, but reduce the mass of the inner planet by a factor of 10 in order to reduce the level of planet–planet interactions. All of the other system parameters are kept the same as before, and all of the systems are started in the 2:1 mean motion resonance. The level of turbulence, as set by the amplitude (Δv)/v of velocity perturbations, is also the same as before.

The resulting time evolution for the fraction of bound states is shown in Figure 7 for an ensemble of N_ens = 10³ systems. The top panel shows the early time evolution as a log–log plot, so that power-law decay corresponds to a straight line. During this early phase of evolution, planetary interactions (which are much weaker than in the true GJ876 system) do not yet have time to act and cannot lead to drift of the resonance location. As a result, the straight dotted line in the top panel provides a good fit to the numerical results, with the expected power-law slope of −1/2 (see Section 2 and Paper I). At later times, the effects of planet–planet interactions gradually accumulate, and their effect on the dynamics becomes important. As shown above, at this stage, the time evolution of the fraction of bound states turns over to an exponentially decreasing form. This behavior is illustrated in the bottom panel of Figure 7, which is presented as a log–linear plot, so that straight lines correspond to exponential decay. Here, the dotted line shows a purely exponential decay, which is expected for the late time behavior of these systems. The dashed curve in the bottom panel shows the solution from Equation (89), where the parameters are adjusted to provide a good fit. Note that the model developed herein smoothly connects the noninteractive regime, with a power-law decrease in the fraction of resonant states, with the highly interactive regime, with an exponential decrease in the fraction of bound states.

Zoom In Zoom Out Reset image size

In order to derive the simplified Fokker–Planck equation, Equation (87), and its solution for the time evolution of the fraction of bound states (Equation (89)), we have made a number of approximations, and it is useful to summarize them here: first, we averaged over the libration angle ϕ in the Fokker–Planck equation; this approximation is well known and well tested (e.g., MM04 and Paper I). At the next stage of the calculation, we adopted the ansatz of Equations (79) and (80). This approximation captures the basic behavior of the eccentricity variations for purposes of modeling the diffusion of the energy of the resonance; if we needed the full time dependence of the eccentricities, however, their full equations of motion should be retained. With this simplified description for the eccentricities, we then averaged over the apsidal angle θ = ϖ₂ − ϖ, and thereby obtained a diffusion equation in two variables (e and V). Although the resulting two-dimensional diffusion equation can be solved, it is rather cumbersome. In addition, the eccentricity values are not free to diffuse to arbitrary values, but rather are confined to the range 0 < e < e_max, where planets with e > e_max tend to experience orbit crossing and ejection. We thus average the diffusion equation over the eccentricity to obtain the simplified form of Equation (86). Although this procedure represents an uncontrolled approximation, and hence this last averaging is the least justified, the uncertainties can be encapsulated in the constants A and B. The solutions to the resulting diffusion equation can then be found exactly, and they contain additional behavior, beyond that of the solutions found in Section 2. In particular, these solutions show that the fraction of bound states can decrease exponentially when the interactions between planets are significant, and this trend agrees with the results of the corresponding numerical simulations of the full three-body problem. Finally, we note that Equation (87) is more robust than its derivation—this form of the Fokker–Planck equation represents a straightforward generalization of that obtained earlier for the simpler stochastic pendulum.

Building on earlier work (Paper I), this paper explores the effects of stochastic fluctuations on the maintenance of mean motion resonance in planetary systems. The principal result of this study is that turbulence can readily compromise mean motion resonance for the fluctuation amplitudes predicted by MHD simulations and for the solar system architectures observed in extrasolar planetary systems. Furthermore, we can understand the physical processes involved using an (primarily) analytic approach. A more specific outline of our results is given below:

The most important quantity calculated in this paper is the fraction of systems that remain in resonance as a function of time. This fraction of bound states P_b(t) is defined to be the number of systems in resonance divided by the number of systems that remain intact (without ejecting a planet). This distinction is necessary because mean motion resonance protects multiple planet systems from instability, so that planetary ejection takes place with a significant probability after a mean motion resonance is compromised. As a result, the fraction of systems remaining in resonance depends on whether the fraction is calculated relative to the original number of systems in the ensemble or the number of multiple planet systems remaining (Figure 1).

The effects of stochastic fluctuations on mean motion resonance act in qualitatively different ways in systems that are highly interactive and those which are not (compare Figure 2 with Figure 3). In systems with relatively little interaction between the planets, the fraction of bound states (systems in resonance) decreases with time as a power law, specifically P_b ∼ t^−1/2. Systems that display this type of behavior are close to the idealized, circular restricted three-body problem; the example considered in this paper contains a Jupiter in a nearly circular orbit (initially) in resonance with a Super-Earth on an interior orbit. The simple pendulum model of mean motion resonance (MD99) with turbulent forcing (Paper I) also derives from the circular restricted three-body problem and leads to this same power-law time dependence (see Equation (18)). For these noninteractive systems, ensembles of the stochastic pendulum (Sections 2.1 and 2.2), solutions to the Fokker–Planck equation (Section 2.3), and direct numerical integrations of the three-body problem (Section 2.4) all predict the same time dependence for the fraction of systems remaining in mean motion resonance. In addition, the departures of the numerical results from the expectation values are well characterized by the amplitude of root-N fluctuations (Figure 3). For this regime, we can summarize these results by writing the expected fraction of surviving resonances in the form P_b ≈ C/N_orb^1/2, where N_orb is the total number of orbits for which turbulence is active, and where the dimensionless constant C depends on the amplitude of the fluctuations and the probability ρ_ret of returning to resonance after leaving (we expect C ∼ 10–100; see Paper I). Thus, for N_orb ≈ 10⁶, we expect P_b ≈ 0.01–0.1.

In highly interactive systems, the fraction of bound states (systems in resonance) decreases more rapidly with time (see Figure 2) and eventually shows an exponential decay. In this paper, we have developed a generalized treatment for mean motion resonance that includes planetary interactions in the model equations (Section 4). This treatment initially retains six variables, and thus contains more physics than the simple one-variable pendulum model; in particular, we include terms that describe the interaction between planets due to their mutual excitation of eccentricities. In this model, the fraction of bound states P_b(t) shows two limiting regimes of behavior: for the case of minimal planet–planet interactions, or for sufficiently short timescales while the diffusion effects dominate, the bound fraction shows power-law behavior P_b ∼ t^−1/2. For highly interactive systems, or for sufficiently late times, the evolution of the distribution function is dominated by drift terms due to interactions, and the model predicts an exponential decrease in the fraction of bound states (Equation (94)). This predicted exponential behavior is in good agreement with that indicated by numerical experiments of highly interactive systems such as GJ876 (see Figures 2 and 6). For these systems, the variations in the numerical results for the bound fraction P_b(t) are somewhat larger than expected from root-N fluctuations alone; this complication is most likely due to the planetary interactions, which allow for additional degrees of freedom.

Systems that contain moderate levels of planet–planet interactions can fall in an intermediate regime, where the fraction of bound states initially decreases as a power law with P_b ∼ t^−1/2 for a substantial time interval, but eventually switches over to an exponential decay (see Figure 7). Furthermore, the duration of the initial power-law phase depends on the level of interactions in the system. During the early time evolution, interactions have little effect, and the system acts essentially as the simple stochastic pendulum; at later times the effects of interactions accumulate, and the fraction of bound states decays exponentially. These two regimes are connected at intermediate times, when P_b(t) behaves as a steeper power law in time (see Equation (92)). The model developed in Section 4 accounts for this early power-law behavior, the later exponential decay, and the smooth matching between the regimes.

Finally, for completeness, we have also studied the ramifications of including a damping term in the model equations for mean motion resonance (Section 3). For example, this damping can be driven by torques from the circumstellar disk that is responsible for planetary migration. In this case, for sufficiently short timescales, the distribution function for the ensemble of states evolves like that of the stochastic pendulum; the energy expectation value increases linearly with time (Equation (32)) and the fraction of bound states decreases as a power law P_b ∼ t^−1/2 (Equation (39)). At later times when γt ≫ 1, however, both the energy expectation value (Equation (32)) and the fraction of bound states (Equation (40)) approach asymptotic values. The analytic solutions to the Fokker–Planck equation (Section 3) are in excellent agreement with results obtained from numerical integrations of the stochastic pendulum equation with damping (see Figures 4 and 5).

This paper extends the analysis of Paper I and bolsters its main conclusion, i.e., that turbulence implies mean motion resonances in extrasolar planetary systems should be relatively rare (roughly at the level of a few percent), unless turbulence has a limited duty cycle. The fact that turbulence is capable of removing systems from resonance is not surprising. As shown in Equation (7), the potential energy associated with a bound resonant state is much smaller (typically by a factor of ∼10⁴) than the binding energy of the planet within the gravitational potential well of the star. Since planets can be ejected with relative efficiency during the early phases of solar system evolution (e.g., Rasio & Ford 1996), it makes sense that additional perturbations (here, turbulent fluctuations) can remove systems from resonance.

Given that the results of this paper (and Paper I) suggest that mean motion resonances should be rare, it is useful to compare this prediction to existing observational data. However, we first note that a detailed comparison is premature, due to the small number of systems found to date, and due to selection effects. In the current sample, 30 systems are observed to have multiple planets. As outlined in the Introduction, one system (GJ876) is found to be deep in a 2:1 resonance, whereas three other systems are either in resonance or close. If all three of these latter systems are actually in resonance, one "estimate" for the nominal rate of resonances would be 4/30 or 13%; if only GJ876 is truly in resonance, this rate becomes 1/30 or 3%. These estimates suggest that mean motion resonances are at least somewhat rare. Besides the small numbers involved, this percentage is uncertain for several reasons: one important issue is that the fraction of systems in mean motion resonance calculated here is the survival rate, i.e., the fraction of systems that start in resonance and are not removed via turbulence. A large fraction of the systems that leave resonance eject one of the planets, and hence would not be included in the number of observed multiple planet systems. Acting in the other direction, another consideration is that not all of the multiple planet systems in the current sample were ever (necessarily) in resonance. In addition, the observations are not complete, so that many of the single-planet systems in the current sample might have companions (perhaps with nearly integer period ratios), which could either add to the number of resonant systems or add to the number that are not in resonance. Our understanding of these issues will benefit from future observations and hence better statistics.

This work also shows that for systems that leave resonance, experience orbit crossing, and eject a planet, the surviving planet typically displays an eccentric orbit. Although a detailed statistical description of the resulting distribution is beyond the scope of this paper, we note that the full range of possible eccentricities is realized. This finding is consistent with previous simulations of multiple planet systems undergoing inward migration; in this setting, one of the planets is often scattered out of the system, and the remaining planets (for an ensemble of systems) are left with semimajor axes and eccentricities that fill the (a, e) plane in a manner consistent with the current observational sample (Adams & Laughlin 2003; Moorhead & Adams 2005). Preliminary numerical experiments of this process including turbulent fluctuations (Moorhead 2008) indicate that the orbital elements of the surviving planets continue to fill the (a, e) plane, but better statistics are needed (both from the simulations and for the observational sample). In any case, planetary scattering and ejection—perhaps enhanced by turbulence removing systems from resonance—provides one viable mechanism to increase the orbital eccentricity of extrasolar planets.

Extrasolar planetary systems display a great deal of diversity, and the effects of turbulence will be different for varying solar system architectures. The results of this paper show that highly interactive systems (such as GJ876) can reach a state of evolution where the fraction of resonant systems decays exponentially. Highly interactive systems are thus the least likely to be able to remain in resonance. This finding, however, makes the existence of GJ876 itself, which is observed to be deep in resonance, an even greater enigma. One way to account for the existence of the 2:1 resonance in the GJ876 system is for its disk to have had a low mass (to reduce the level of turbulent fluctuations) and/or a short survival time after planet formation (to reduce the duty cycle of the effect). However, this explanation for the survival of the resonance would make it more difficult for the disk to produce relatively large planets. Recall that both theoretical considerations (Laughlin et al. 2004) and observational findings (e.g., Johnson et al. 2007) suggest that M stars have difficulty forming planetary companions with Jovian masses.

The results of Section 4 suggest that systems with the smallest levels of interaction between planets would have the best chance of surviving in a resonant state. However, this hypothesis is not entirely true: for a given orbital spacing, for example that corresponding to a 2:1 period ratio, the level of interactions decreases as the planetary masses decrease. On the other hand, planets with the smallest masses produce the smallest gaps in their circumstellar disks, and hence experience greater levels of turbulent forcing; in other words, the disk reduction factor Γ_R due to gaps is smaller for low-mass planets (Paper I). Planetary systems with only lower-mass planets, such as the newly discovered HD 40307 system with masses of m_P = 4–9 M_E (Mayor et al. 2009), may experience greater levels of turbulent forcing than systems with larger planets, and hence would be more likely to be removed from resonance. Taken together, these theoretical results to date suggest that the systems with the best chances of maintaining mean motion resonance in the face of turbulence are those with planets of moderate mass and nearly circular orbits. The "Super-Earth" system of this paper—with a Jovian planet in an outer, nearly circular orbit and a smaller inner planet—provides one good example (see Figure 3).

The numerical simulations indicate another difference between the highly interactive and less interactive systems. After leaving a resonant state, the Super-Earth systems take a long time for their planets to experience orbit crossing and eventual planetary ejection. The GJ876-type systems eject their planets much more readily. As a result, the less interactive systems (e.g., Super-Earth systems) can stay "close" to resonance for a long time, where the period ratios are close to 2:1, but the other orbital elements do not have the proper values to be in true resonance.

This paper has made progress on our understanding of how well the simplest pendulum equation works as a model for mean motion resonances, including stochastic forcing terms, and we have generalized the model to include damping (Section 3) and planetary interactions in an approximate manner (Section 4). Nonetheless, a great deal of work remains to be done. One important issue is to understand how turbulence acts in the earliest stages of planet migration, i.e., when the planets are formed and first become locked into resonance. This issue requires much more extensive numerical work and is thus beyond the scope of this paper. Recent work (see Moorhead 2008) is starting to study how turbulence affects the early phases of migration. The issue of turbulence during planetary formation is also being considered elsewhere (e.g., Masahiro et al. 2007). In addition, we have not included the back reaction of how planets (e.g., through gap clearing) can affect the generation of turbulent fluctuations through the MRI. All of these issues provide fruitful avenues for future research.

In the coming years, as the statistics of observed multiple planet systems become sufficiently complete, the fraction of systems in resonance will provide important constraints on their formation and subsequent evolution. Since mean motion resonance can be compromised relatively easily, as shown herein, planetary systems observed in such bound states must have followed restricted historical paths (e.g., with little interaction and/or large dissipation). This type of analysis has been considered for the resonances in our Solar System (Peale 1976; Malhotra 1993) and for the GJ876 system (e.g., Lee & Peale 2002), and will soon provide interesting constraints on other extrasolar planetary systems. However, the greatest wealth of information on this topic will come when the observed sample of multiple planet systems is large enough to determine the fraction of actual mean motion resonances and the fraction of near resonances.