OVERSTABLE LIBRATIONS CAN ACCOUNT FOR THE PAUCITY OF MEAN MOTION RESONANCES AMONG EXOPLANET PAIRS

As of 2013 January, more than 3000 planet candidates have been discovered by the Kepler spacecraft and over 1000 of these reside in systems that contain more than one planet (Batalha et al. 2013). Thanks to these remarkable findings, we now have a large statistical sample of multi-planet systems. Their architectures contain important clues concerning planet formation and dynamical evolution. Kepler multi-planet systems display the following characteristics.

• 1.  
  Most planets do not reside in or close to mean motion resonances as shown in Figure 1 (Fabrycky et al. 2012).
• 2.  
  There is a significant excess of planet pairs with period ratios close to but slightly larger (by 1%–2%) than that for exact resonance (see Figure 1 and Fabrycky et al. 2012).

Zoom In Zoom Out Reset image size

Several papers have been published offering explanations for the 1%–2% offset from mean motion resonance (e.g., Lithwick & Wu 2012; Batygin & Morbidelli 2013; Petrovich et al. 2013; Baruteau & Papaloizou 2013). Our investigation has a different focus. We aim to explain why mean motion resonances are rare. This is surprising because planet–disk interactions leading to orbit migration are expected to result in efficient resonance capture (see Section 2.2 and Figure 4). However, we show that for most Kepler planet pairs, resonance capture is only temporary. The reason is that librations about exact resonance are overstable and lead to passage through resonance on the eccentricity damping timescale which is about two orders of magnitude shorter than that for semimajor axis migration (e.g., Goldreich & Tremaine 1980; Ward 1988).

This paper is structured as follows. Section 2 focuses on the planar, circular, restricted three-body problem in the vicinity of a mean motion resonance including resonance capture in the presence of semimajor axis migration and eccentricity damping and the evolution within and escape from resonance. We show that together, eccentricity damping and orbit migration give rise to an equilibrium eccentricity and that librations around equilibrium are overstable for the majority of Kepler planet masses and lead to passage through resonance. We extend our model to the planar, three-body problem in Section 3. Results derived from our model are compared with properties of multi-planet systems discovered by the Kepler spacecraft in Section 4. In Section 5 we discuss reasons for the departure from exact resonance and show that our resonance model yields an excess of pairs with slightly greater than exact resonant period ratios, but that the magnitude of the offset from exact resonance is on average too small to match observations. We compare our findings with previous works that aim to explain the period offset. Our main conclusions are summarized in Section 6.

Consider a system of two planets in orbit around a host star. Assume that the outer planet moves on a fixed circular orbit. Close to a first order j + 1: j mean motion resonance, the dominant term in the inner planet's disturbing function has resonant argument ϕ such that

$$\begin{eqnarray} \phi &=&(j+1)\lambda ^{\prime }-j\lambda -\varpi, \quad \ \dot{\phi }=(j+1)n^{\prime }-jn-\dot{\varpi },\nonumber\\ \ddot{\phi }&=&-j\dot{n}-\ddot{\varpi }, \end{eqnarray} \tag{ 1 }$$

where primes distinguish the outer planet and λ, n, and ϖ denote mean longitude, mean motion, and longitude of pericenter. At conjunction, λ' = λ so ϕ is a measure of the displacement of the longitude of conjunction from the inner planet's pericenter.

Lagrange's equations of motion to first order in eccentricity for the inner planet in the vicinity of the resonance read

$$\begin{equation} \dot{n}=3j\beta \mu ^{\prime } en^2\sin \phi -\frac{n}{\tau _{n}}+p\frac{e^2n}{\tau _{e}}, \end{equation} \tag{ 2 }$$

$$\begin{equation} \dot{e}=\beta \mu ^{\prime } n\sin \phi -\frac{e}{\tau _{e}}, \end{equation} \tag{ 3 }$$

$$\begin{equation} \dot{\varpi }=-\frac{\beta \mu ^{\prime }}{e}n\cos \phi, \end{equation} \tag{ 4 }$$

$$\begin{eqnarray} \ddot{\phi }&=&{-}3 j^2 \beta \mu ^{\prime }en^2\sin \phi -\left(\frac{\beta \mu ^{\prime } n}{e}\right)^2 \cos \phi \sin \phi\nonumber\\ &&{-}\frac{\beta \mu ^{\prime } n}{e}\sin \phi \dot{\phi }+\frac{\beta \mu ^{\prime } n}{e\tau _e}\cos \phi +j\frac{n}{\tau _n}-3j\frac{e^2n}{\tau _e}, \end{eqnarray} \tag{ 5 }$$

where e and a are eccentricity and semimajor axis, μ' is the ratio of the outer planet's mass to that of the star, $\tau _n \equiv n/|\dot{n}|>0$, $\tau _e \equiv e/|\dot{e}|>0$, and β ≈ 0.8j. In Equation (2) the sign of τ_n is chosen for convergent migration. The final term in Equation (2) accounts for the contribution of eccentricity damping to changing the mean motion. For the particular case of eccentricity damping arising from energy dissipation at constant angular momentum, as applies to tides raised in a planet by its parent star, p ≈ 3. For simplicity, we assume p = 3 throughout the body of our paper but provide results applicable for general p > 0 in the Appendix. Although we will first examine the dynamics of the resonance without dissipation, we include the eccentricity damping and migration terms in Equations (2) and (3) such that we do not have to repeat the above equations later in this section.

Near resonance, e = O(μ'^1/3) and d/dt = O(μ'^2/3n). Thus the first and second terms on the right-hand side (rhs) of Equation (5) dominate over the third for small amplitude librations (i.e., ϕ = δϕ ≪ 1). In the absence of dissipation, librations of ϕ are centered on ϕ = 0, e = e₀ (e.g., Murray & Dermott 1999) where

$$\begin{equation} e_0=\frac{\beta \mu ^{\prime } n}{jn-(j+1)n^{\prime }}. \end{equation} \tag{ 6 }$$

From Equation (5), it follows that the frequency of small librations, ω, satisfies

$$\begin{equation} \frac{\omega ^2}{n^2}=3j^2\beta \mu ^{\prime } e_0+\left(\beta \mu ^{\prime }\over e_0\right)^2. \end{equation} \tag{ 7 }$$

In the absence of migration and eccentricity damping, τ_n and τ_e → ∞, Equations (2)–(4) admit two integrals⁴

$$\begin{equation} k(\phi,e^2)= \left(\frac{3}{2}j^2 e^2-\frac{\beta \mu ^{\prime }}{e}\cos \phi \right) +\frac{\dot{\phi }}{n}, \end{equation} \tag{ 8 }$$

and

$$\begin{equation} \mathcal {H}(\phi,e^2)=\left(k e^2-\frac{3}{4}j^2e^4 +2\beta \mu ^{\prime }e\cos \phi \right). \end{equation} \tag{ 9 }$$

The integral k exists because with only one resonant argument, variations of n and e are related. $\mathcal {H}$, the Jacobi constant, is also the Hamiltonian with canonically conjugate momentum, e², and coordinate, ϕ.

The topology of the phase-space as defined by $\mathcal {H}$ for fixed k changes abruptly across k = k_crit, where

$$\begin{equation} k _{{\rm crit}}= \frac{3^{4/3}}{2}(j\beta \mu ^{\prime })^{2/3}\sim j^{4/3}\mu ^{\prime 2/3}. \end{equation} \tag{ 10 }$$

For k < k_crit, there is one (stable) fixed point whereas there are three (two stable and one unstable) fixed points for k > k_crit.⁵ These distinct topologies are illustrated in Figure 2. For k > k_crit the level curve emanating from the unstable fixed point is appropriately called a separatrix because it separates the regions surrounding the two stable fixed points. The stable fixed point at

$$\begin{eqnarray} \phi &=&0 \quad {\rm and}\quad e_0=\frac{2}{j^{2/3}} \left(\frac{\beta \mu ^{\prime }}{3}\right)^{1/3}\left(k\over k_{{\rm crit}}\right)^{1/2}\nonumber\\ &=&\frac{\beta \mu ^{\prime } n}{jn-(j+1)n^{\prime }} \end{eqnarray} \tag{ 11 }$$

is present for all k and corresponds to a periodic orbit with jn > (j + 1)n'. It is located at the global maximum of $\mathcal {H}$ and owes its stability to the Coriolis acceleration in the frame rotating with angular velocity n'. The other stable fixed point at

$$\begin{eqnarray} \phi &=&\pi \quad {\rm and}\quad e_{{\rm min}}=\frac{1}{j^{2/3}}\left(\frac{\beta \mu ^{\prime }}{3}\right)^{1/3}\left(k\over k_{{\rm crit}}\right)^{1/2}\nonumber\\ &&\times(\cos \theta -3^{1/2}\sin \theta )=\frac{\beta \mu ^{\prime } n}{(j+1)n^{\prime }-jn}, \end{eqnarray} \tag{ 12 }$$

with

$$\begin{equation} \theta =\frac{1}{3}\tan ^{-1}[(k/k_{{\rm crit}})^3-1], \end{equation} \tag{ 13 }$$

only appears for k > k_crit and describes a periodic orbit with jn < (j + 1)n'. It sits at a local minimum of $\mathcal {H}$. Finally, the unstable fixed point at

$$\begin{eqnarray} \phi &=&\pi \quad {\rm and}\quad e=\frac{1}{j^{2/3}}\left(\frac{\beta \mu ^{\prime }}{3}\right)^{1/3}\left(k\over k_{{\rm crit}}\right)^{1/2}\nonumber\\ &&\times(\cos \theta +3^{1/2}\sin \theta)=\frac{\beta \mu ^{\prime }}{(j+1)n^{\prime }-jn} \end{eqnarray} \tag{ 14 }$$

lies on a saddle point. The stable and unstable fixed points at ϕ = π bifurcate at k = k_crit where e_crit = (βμ'/3j²)^1/3. These values correspond to the boundary between stable and unstable linear oscillations around ϕ = π, e = βμ'n/((j + 1)n' − jn) as can be verified from Equation (5). In anticipation of later discussion, we make note of the asymmetry of the stable fixed points on opposite sides of exact resonance.

Zoom In Zoom Out Reset image size

2.1. Planet–Disk Interactions

Up to this point, migration and eccentricity damping were merely parameterized by their respective timescales τ_n and τ_e. For dissipation due a planet's interaction with a protoplanetary disk, τ_n and τ_e are given by

$$\begin{equation} \frac{1}{\tau _{n}} \equiv \frac{1}{n}\frac{dn}{dt} \sim \mu \mu _d \left(\frac{a}{h}\right)^2n, \end{equation} \tag{ 15 }$$

$$\begin{equation} \frac{1}{\tau _{e}} \equiv \frac{1}{e} \frac{de}{dt} \sim \mu \mu _d \left(\frac{a}{h}\right)^4n, \end{equation} \tag{ 16 }$$

where μ_d = Σ_da²/M_* is the disk to star mass ratio and h is the disk's scale height. Here we have assumed that the planet's mass is too small for it to clear a gap in the disk.

The migration rate is based on the balance of torques at principal Lindblad resonances located interior and exterior to the planet (Goldreich & Tremaine 1980; Ward 1986). Eccentricity damping in a gapless disk is dominated by first order, coorbital Lindblad resonances (Ward 1988; Artymowicz 1993). These come in pairs with the individual members approximately equally displaced interior and exterior from the planet's orbit. Consequently, each pair damps the planet's orbital energy at nearly constant angular momentum. However, remote first order Lindblad resonances along with corotation resonances also contribute to eccentricity excitation and damping. Thus the coefficient p in Equation (2) deviates from 3. For example, Tanaka & Ward (2004) obtain p ∼ 1.3 from their 3D planet–disk model. Overstable librations require that the final term in Equation (2) be positive as shown by Equation (A2).

2.2. Capture Into Resonance

Next we consider conditions under which the inner planet would be captured into the 2:1 resonance if it were migrating outward at a rate ${|\dot{n}|}/n=1/\tau _n$. Capture would be guaranteed provided the orbital eccentricity of the inner body far from the resonance were sufficiently small and the rate of outward migration sufficiently slow. Let us assume that the former condition is satisfied. Then upon approach to resonance, the inner planet would stay close to the stable fixed point at ϕ = 0 and e₀ = βμ'n/(jn − (j + 1)n'). As the resonance is approached, ω first decreases and then increases (see Figure 3). Capture is most problematic during passage through

$$\begin{equation} \frac{\omega _{{\rm min}}}{n}=\frac{3^{4/3}}{2^{1/3}}(j\beta \mu ^{\prime })^{2/3}\sim (j^2\mu ^{\prime })^{2/3}, \end{equation} \tag{ 17 }$$

at which point

$$\begin{equation} e_0(\omega _{{\rm min}})=\left(2\beta \mu ^{\prime }\over 3j^2\right)^{1/3}\sim \left(\mu ^{\prime }\over j\right)^{1/3}, \end{equation} \tag{ 18 }$$

and

$$\begin{equation} \frac{\Delta n}{n}\equiv 1-\frac{(j+1)n^{\prime }}{jn}=\left(3\over 2j\right)^{1/3}(\beta \mu ^{\prime })^{2/3}\sim j^{1/3}\mu ^{\prime 2/3}. \end{equation} \tag{ 19 }$$

To cross the resonance width, Δn/n ∼ j^1/3μ'^2/3 given by Equation (19) while migrating at rate $|\dot{n}|/n=1/\tau _n$ takes Δt ∼ j^1/3μ'^2/3τ_n. Capture would be assured provided ω_minΔt ∼ j^5/3μ'^4/3nτ_n ≫ 1. Careful analytical and numerical calculations sharpen this result to

$$\begin{equation} j^{5/3}\mu ^{\prime 4/3}n\tau _n\ge 2.5. \end{equation} \tag{ 20 }$$

For smaller τ_n, temporary capture could occur at larger separations from resonance where ω > ω_min but the resonance lock would be broken before ω_min were reached. With τ_n given by Equation (15), the condition of certain capture becomes

$$\begin{equation} \mu ^{\prime } \gtrsim \frac{\mu _d^3}{j^5}\left(\frac{a}{h}\right)^6. \end{equation} \tag{ 21 }$$

Zoom In Zoom Out Reset image size

Having defined sufficiently slow migration, we do the same for sufficiently small initial eccentricity. The key is to make use of the adiabatic invariant

$$\begin{equation} {\rm AI}=\oint _C \, d\phi \,e^2. \end{equation} \tag{ 22 }$$

For motion on the separatrix at k = k_crit, it can be shown that AI_sep = 12π(βμ'/3j²)^2/3. Thus if well before a planet encounters the resonance, its orbital eccentricity were smaller than

$$\begin{equation} e_{{\rm cap}}=\left({\rm AI}_{{\rm sep}}\over 2\pi \right)^{1/2}=1.7\left(\beta \mu ^{\prime } \over j^2\right)^{1/3}\sim \left(\mu ^{\prime } \over j\right)^{1/3}, \end{equation} \tag{ 23 }$$

then at k = k_crit, $\mathcal {H}$ would be larger than ${\mathcal {H}}_{{\rm sep}}$ and capture would be guaranteed.

Conditions for resonance capture appropriate to planets migrating in a protoplanetary disk are displayed in Figure 4 as a function of the local disk and planet masses. The local disk mass, μ_d, is that which resides within a factor of two of the planet's orbital radius. For a total disk mass equal to 10⁻² M_* distributed with a radial surface density profile ∝a⁻¹, the local disk mass at 0.5 AU is two orders of magnitude smaller than the local disk mass at 50 AU. Figure 4 indicates that resonance capture will only occur for local disk masses less than 10⁻³ M_*. Thus in a typical disk with total mass 10⁻² M_* and an a⁻¹ surface density profile, only planets within about 10 AU would be caught into resonance. Figure 4 further suggests planets with masses between those of Earth and Neptune which orbit within 1 AU have the best chances to capture smaller planets.

Zoom In Zoom Out Reset image size

To a significant extent, the current situation regarding mean motion resonances among pairs of exoplanets is analogous to that which pertained to mean motion resonances in the satellite systems of Jupiter and Saturn more than half a century ago. A early paper by Roy & Ovenden (1954) called attention to the overabundance of orbital resonances among the satellite systems of Jupiter and Saturn. Goldreich (1965) showed that slow convergent migration driven by tides raised in a planet would lead to the formation of mean motion resonances stabilized by the transfer of angular momentum from the inner to the outer satellite. Early work emphasized slow migration because it is appropriate to cases in which tidal torques drive the expansion of satellite orbits. The probability of capture into resonance for slow convergent migration was first solved by Yoder (1979). Subsequently, simpler derivations emphasizing the role of adiabatic invariance were provided in Henrard (1982) and Borderies & Goldreich (1984). Recent investigations focused on the maximum migration rate permitting capture. The analytic derivation by Friedland (2001) is in close agreement with results from numerical experiments (Quillen 2006; Ogihara & Kobayashi 2013).

2.3. Evolution in Resonance in the Presence of Migration

2.4. Effects of Eccentricity Damping

Evolution in resonance is very different in the presence of eccentricity damping. This is because as the planet evolves deeper into resonance the associated increase in its eccentricity can be balanced by a decrease due to eccentricity damping. Equations (2) and (3) imply the existence of an equilibrium eccentricity and resonance phase angle given by

$$\begin{equation} e_{{\rm eq}}= \left(\frac{\tau _{e}}{3(j+1)\tau _{n}}\right)^{1/2} \end{equation} \tag{ 24 }$$

and

$$\begin{equation} \sin \phi _{{\rm eq}} = \frac{e_{{\rm eq}}}{\beta \mu ^{\prime }\tau _{e} n}. \end{equation} \tag{ 25 }$$

No such equilibrium exists in the presence of migration alone.

In what follows we show that evolution in resonance with eccentricity damping and migration has three possible outcomes which depend upon the values of βμ' and τ_e/τ_n.

• 1.  
  For

  $$\begin{equation} \mu ^{\prime } > \frac{j}{\sqrt{3}(j+1)^{3/2}\beta }\left(\tau _{e}\over \tau _{n}\right)^{3/2}, \end{equation} \tag{ 26 }$$

  the planet is permanently trapped in resonance at the fixed point ϕ = ϕ_eq and e = e_eq. An example is shown in the upper panel of Figure 5.
• 2.  
  For

  $$\begin{eqnarray} \frac{j^2}{8\sqrt{3}(j+1)^{3/2}\beta }\left(\tau _{e}\over \tau _{n}\right)^{3/2} &<& \mu ^{\prime } < \frac{j}{\sqrt{3}(j+1)^{3/2}\beta }\nonumber\\ &\times &\left(\tau _{e}\over \tau _{n}\right)^{3/2}, \end{eqnarray} \tag{ 27 }$$

  the planet is permanently caught in resonance and its libration amplitude saturates at a finite value. The middle panel of Figure 5 displays an example.
• 3.  
  For

  $$\begin{equation} \mu ^{\prime } < \frac{j^2}{8\sqrt{3}(j+1)^{3/2}\beta }\left(\tau _{e}\over \tau _{n}\right)^{3/2}, \end{equation} \tag{ 28 }$$

  the planet is caught in resonance, but escapes on timescale τ_e. An example is provided in the lower panel of Figure 5.

Zoom In Zoom Out Reset image size

Note that the boundary between permanent capture and temporary capture leading to escape shifts monotonically to smaller μ'∝j^−1/2 with increasing j. Moreover, the intermediate range of μ' corresponding to finite amplitude librations shrinks with increasing j and does not exist for j ⩾ 8.

Figure 6 summarizes the three different outcomes for the evolution in resonance for various planet to star mass ratios and as a function of τ_e/τ_n for a 2:1 mean motion resonance.

Zoom In Zoom Out Reset image size

Below we derive analytic results for different evolutions in resonance.

2.4.1. Instability: Transition from Damping to Growth

Linearizing Equations (1) for $\dot{\phi }$, (2) for $\dot{n}$, and (3) for $\dot{e}$ around the equilibrium given by ϕ_eq and e_eq and setting d/dt = s yields a system of three homogeneous linear equations in the variables δϕ, δn and δe. The determinant of the coefficient matrix is a cubic polynomial in s with real coefficients whose roots are the eigenvalues. A standard exercise in algebra reveals that one root is negative and that the other two are a complex conjugate pair. For slow migration and eccentricity damping, ωτ_e ≫ 1 and ωτ_n ≫ 1, the imaginary parts of the complex roots are essentially ±iω evaluated at e₀ = e_eq and the real part is well-approximated by

$$\begin{equation} s=\frac{1}{\tau _{e}}\left(\frac{n}{\omega _{{\rm eq}}}\right)^2\left(3j\beta \mu ^{\prime }e_{{\rm eq}}-\left(\frac{\beta \mu ^{\prime }}{e_{{\rm eq}}}\right)^2\right). \end{equation} \tag{ 29 }$$

Thus the evolution in resonance switches from damping to growth at e_eq = (βμ'/3j)^1/3. Since e_eq = (τ_e/(3(j + 1)τ_n))^1/2, the libration amplitude grows if

$$\begin{equation} \mu ^{\prime }<\frac{j}{\sqrt{3}(j+1)^{3/2}\beta }\left(\tau _{e}\over \tau _{n}\right)^{3/2}. \end{equation} \tag{ 30 }$$

Furthermore, the growth timescale from Equation (29) is of order τ_e. Thus, escape from resonance, if it occurs, will happen on a timescale τ_e.

At a basic level, overstability can arise because e₀ corresponds to a maximum of $\mathcal {H}$ and eccentricity damping causes the system to slide downhill. It is necessary to dig a little deeper to understand the origin of the transition from damping to overstability. Once again, an adiabatic invariant comes into play, this time involving small librations around ϕ = 0, e = e₀.⁶ A simple calculation using Equation (22) yields

$$\begin{equation} {\rm AI}_{e_0}=2\pi \beta \mu ^{\prime }\frac{n}{\omega _0}e_0\Phi ^2, \end{equation} \tag{ 31 }$$

where Φ is the amplitude of the libration and ω₀ its frequency. For e₀ ≪ μ'^1/3, ω₀/n ∼ μ'/e₀ whereas for e₀ ≫ μ'^1/3, ω₀/n ∼ (μ'e₀)^1/2. By itself, convergent migration results in a monotonic increase of e₀. The adiabatic invariant then implies that the libration amplitude decreases as e₀ increases. But for e₀ ≪ μ'^1/3, $\Phi \propto e_0^{-1}$ whereas for e₀ ≫ μ'^1/3, $\Phi \propto e_0^{-1/4}$. Eccentricity damping leads to an increase in Φ which migration can overcome for e₀ ≪ μ'^1/3 but not for e₀ ≫ μ'^1/3.

2.4.2. Escape from Resonance

Following capture in resonance, escape involves damping to the stable fixed point at ϕ = π, e = e_min. This can only occur if k > k_crit and involves crossing the inner branch of the separatrix which leads to damped circulation about the enclosed fixed point. As shown at the beginning of Section 2, the separatrix exists for e₀ ⩾ 2(βμ'/3j²)^1/3. Equating e_eq = (τ_e/(3(j + 1)τ_n))^1/2 to e₀ demonstrates that resonance escape can occur provided

$$\begin{equation} \mu ^{\prime } < \frac{j^2}{8\sqrt{3}(j+1)^{3/2}\beta }\left(\tau _{e}\over \tau _{n}\right)^{3/2}. \end{equation} \tag{ 32 }$$

Different outcomes for evolution in resonance as a function of μ' and τ_e/τ_n are illustrated in Figure 7 for the particular case of a 2:1 resonance. Plotted points show results from direct numerical integrations of the resonant equations of motion and the solid and dashed lines present our analytic results. The blue circles, which mark the transition from damping to growth, are indistinguishable for τ_nn = 10⁵ and τ_nn = 10⁶ and are in accord with our analytic limit (solid blue line). The black squares and red diamonds show the transition from finite amplitude growth to escape from resonance for τ_nn = 10⁵ and τ_nn = 10⁶, respectively. Although the numerical results for τ_nn = 10⁶ are in close agreement with our analytic expression (dashed red line), those for τ_nn = 10⁵ fall well below the analytic limit for μ' ≲ 10⁻⁴. The latter is not surprising. These combinations of τ_n and μ' do not conform to the basic assumption of our analytic treatment, namely that damping only makes a small perturbation to the dissipation free motion. In particular, we suspect that the problem in this case is that the eccentricity damping timescale is comparable to the libration period.

Zoom In Zoom Out Reset image size

So far we have assumed that the more massive outer planet moves on a fixed circular orbit and that m₁/m₂ ≪ 1.⁷ In this case there is only one resonant argument. Here we investigate a scenario in which both the inner and outer planet are in resonance. We still assume that the outer planet with mass m₂ is the more massive one thus ensuring that type-I migration is convergent (Equation (15)).

3.1. Two-planet System with Convergent Migration and Eccentricity Damping

Consider two planets that orbit their host star in the vicinity of a 2:1 mean motion resonance. The first order resonant terms in the disturbing function for a 2:1 mean motion resonance for the inner and outer planets are, respectively,

$$\begin{equation} \mathcal {R}_1=\alpha a_1^2\mu _2 n_1^2\left(-be_1\cos \phi _1+ce_2 \cos \phi _2\right) \end{equation} \tag{ 33 }$$

and

$$\begin{equation} \mathcal {R}_2=a_2^2\mu _1 n_2^2\left(-be_1\cos \phi _1+ce_2 \cos \phi _2\right), \end{equation} \tag{ 34 }$$

where α = a₁/a₂ = 0.630, b = 1.190, and c = 0.428 and ϕ₁ = 2λ₂ − λ₁ − ϖ₁ and ϕ₂ = 2λ₂ − λ₁ − ϖ₂. In the presence of migration and eccentricity damping

$$\begin{eqnarray} \dot{n}_1&=&{-}\frac{3}{a_1^2}\frac{\partial \mathcal {R}_1}{\partial \lambda _1} +\frac{n_1}{\tau _{n1}}+\frac{3n_1e_1^2}{\tau _{e1}}\nonumber\\ &=&3\alpha \mu _2n_1^2(be_1\sin \phi _1-ce_2\sin \phi _2)+\frac{n_1}{\tau _{n1}}+\frac{3n_1e_1^2}{\tau _{e1}}, \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} \dot{n}_2&=&{-}\frac{3}{a_2^2}\frac{\partial \mathcal {R}_2}{\partial \lambda _2} +\frac{n_2}{\tau _{n2}}+\frac{3n_2e_2^2}{\tau _{e2}}\nonumber\\ &=&{-}6 \mu _1n_2^2(be_1\sin \phi _1-ce_2\sin \phi _2)+\frac{n_2}{\tau _{n2}}+\frac{3n_2e_2^2}{\tau _{e2}}, \nonumber\\ \end{eqnarray} \tag{ 36 }$$

$$\begin{equation} \dot{e}_1=-\frac{1}{n_1a_1^2e_1}\frac{\partial \mathcal {R}_1}{\partial \varpi _1}-\frac{e_1}{\tau _{e1}}=\alpha b\mu _2 n_1\sin \phi _1-\frac{e_1}{\tau _{e1}}, \end{equation} \tag{ 37 }$$

$$\begin{equation} \dot{e}_2=-\frac{1}{n_2a_2^2e_2}\frac{\partial \mathcal {R}_2}{\partial \varpi _2}-\frac{e_2}{\tau _{e2}}=-c\mu _1 n_2\sin \phi _2-\frac{e_2}{\tau _{e2}}, \end{equation} \tag{ 38 }$$

$$\begin{equation} \dot{\varpi }_1=-\frac{1}{n_1a_1^2e_1}\frac{\partial e_1}{\partial \varpi _1}=- \alpha b\mu _2 n_1\frac{\cos \phi _1}{e_1}, \end{equation} \tag{ 39 }$$

$$\begin{equation} \dot{\varpi }_2=-\frac{1}{n_2a_2^2e_2}\frac{\partial e_2}{\partial \varpi _1}=c\mu _1 n_2\frac{\cos \phi _2}{e_2}. \end{equation} \tag{ 40 }$$

Since we are dealing with a 2:1 mean motion resonance, ϕ₁ = 2λ₂ − λ₁ − ϖ₁ and ϕ₂ = 2λ₂ − λ₁ − ϖ₂, which simplify to ϕ₁ = σ − ϖ₁ and ϕ₂ = σ − ϖ₂ where

$$\begin{equation} \dot{\sigma }= 2n_2 -n_1. \end{equation} \tag{ 41 }$$

After substituting for ϕ₁ and ϕ₂, Equations (35)–(41) yield a set of seven first order differential equations in the seven independent variables σ, n₁, n₂, e₁, e₂, ϖ₁, and ϖ₂.

If the planets are caught in an exact 2:1 mean motion resonance, 2n₂ − n₁ = 0 and $\dot{\varpi }_1=\dot{\varpi }_2$ which implies

$$\begin{equation} \frac{e_1}{e_2}= 2\frac{\alpha b}{c}\frac{\mu _2}{\mu _1}. \end{equation} \tag{ 42 }$$

Furthermore, in the presence of eccentricity damping, the librations of the resonant argument are no longer centered exactly at ϕ₁ = 0 and ϕ₂ = π but are offset from 0 and π such that

$$\begin{equation} \sin \phi _1 = \frac{e_1}{\alpha b\tau _{e1} \mu _2 n_1},\ \quad \sin \phi _2 = -\frac{e_2}{c\tau _{e2} \mu _1 n_2}. \end{equation} \tag{ 43 }$$

Assuming that τ_e1/τ_e2 = τ_n1/τ_n2 = μ₂/μ₁, as expected for dissipation due to interactions with a protoplanetary disk (see Equations (15) and (16)), the condition 2n₂ − n₁ = 0 implies

$$\begin{equation} e_1= \left(\frac{\tau _{e1}}{6\tau _{n2}}\right)^{1/2} \left(\frac{1-\mu _1/\mu _2}{(1+\mu _1/(2\alpha \mu _2))(1+(c/b)^2/(4\alpha))}\right)^{1/2}. \end{equation} \tag{ 44 }$$

Figures 8 and 9 display the results of the numerical integration of Equations (35)–(41). Figure 8 shows the evolution of the semimajor axis of the two planets. Convergent migration leads to temporary capture into 2:1 mean motion resonance, but escape from this resonance occurs after a few eccentricity damping timescales. Figure 9 displays the corresponding eccentricity evolution of the two planets (solid lines) and the expected equilibrium eccentricity calculated analytically in Equations (42) and (44) (dashed lines). The numerical and analytical results are in very good agreement.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the limit that μ₁/μ₂ ≪ 1, Equation (44) simplifies to

$$\begin{equation} e_1= \left(\frac{\tau _{e1}}{6\tau _{n2}}\right)^{1/2} \left(\frac{1}{(1+(c/b)^2/(4\alpha))}\right)^{1/2} \simeq 0.98 \left(\frac{\tau _{e1}}{6\tau _{n2}}\right)^{1/2} \end{equation} \tag{ 45 }$$

and the two-planet resonance problem reduces to a one-planet resonance problem in which the inner planet migrates toward the outer one at the rate τ_n2 such that

$$\begin{equation} \dot{n}_1 -2 \dot{n}_2 = 3\alpha b \mu _2n_1^2e_1\sin \phi _1-\frac{n_1}{\tau _{n2}}+\frac{3n_1e_1^2}{\tau _{e1}}. \end{equation} \tag{ 46 }$$

The rhs of the above equation is identical to Equation (2) except the migration rate is that of the outer planet rather than that of the inner one. The equations of motion for $\dot{e}_1$ and $\dot{\varpi }_1$ remain unchanged as given by Equations (37) and (39), respectively. Taking the ratio of Equations (15) and (16) yields

$$\begin{equation} \frac{\tau _{e1}}{\tau _{n2}} \sim \left(\frac{h}{a}\right)^2 \left(\frac{\mu _2}{\mu _1} \right), \end{equation} \tag{ 47 }$$

which evaluates to about 0.01(μ₂/μ₁) for typical disk scale heights h/a ≈ 0.1. The results of this section demonstrate that a simplified treatment of mean motion resonance based on a single active planet is capable of revealing the essential features of the more complex two-planet dynamics, provided τ_n is set equal to $\tau _{n_2}$ and τ_e is set to $\tau _{e_1}$.

3.2. Comparison with Results of other Numerical Integrations

Our investigation implies that pairs of Kepler planets found close to resonance are comprised of those that were permanently captured and others that were temporarily captured but were caught in resonance when the protoplanetary disk dispersed. Given plausible ranges of μ and τ_e/τ_n, it is not surprising that these total only several percent of the pairs in multi-planet systems. This argument rests on the assumption that semimajor axes are largely frozen during and subsequent to disk dispersal.

Figures 10 and 11 plot the more massive planet to star mass ratio, μ', as a function of the planet pair mass ratio, m_</m_>, for all Kepler two-planet systems known as of 2013 June. The subscripts "<" and ">" refer to the smaller and larger of the two planets, respectively. These ratios were calculated based on the assumption that each planet has a mean density of 2 g cm⁻³. Planet pairs above the dashed red lines are candidates for permanent capture in the 2:1 (Figure 10) and 3:2 (Figure 11) resonance. Resonant arguments for those above the solid red lines would be damped whereas those located below it would be overstable. Most of the planet pairs lie below the dashed red lines. These might have undergone temporary capture but would have escaped from resonance on timescale $\tau _{e_<}$. The fate of the pairs between the dashed and solid lines is less certain. For many, the eccentricity damping timescale, $\tau _{e_<}$, is so short that escape would be probable as found in our numerical integrations summarized in Figure 7. Given that planet pairs which ultimately escape resonance only spend a time $\tau _{e_<}$ in each resonance but about a time $\tau _{n_>}$ between successive resonances, the fraction of temporarily captured planet pairs in or near mean motion resonances is $\tau _{e_<}/\tau _{n_>} \sim 0.01 \mu _>/\mu _<$ which evaluates to 3% for the Kepler two-planet sample. Perhaps a comparable or even slightly greater percent might have been permanently captured.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 shows that peaks associated with period ratios close to 2:1 and 3:2 are systematically displaced to larger values by an order of 1%–2%. The direction of displacement is almost certainly due to the asymmetry that requires convergent migration for capture in resonance. However, its magnitude is less easily accounted for. This is illustrated in Figure 12 which shows that to match a 1% average offset from exact 2:1 resonance requires a permanently captured test particle to be paired with a planet for which μ₂ ∼ 4 × 10⁻⁴. More generally, test particles paired in mean motion resonances with massive planets would have period ratios offset from exact resonance by ${\sim} j^{1/3}\mu _2^{2/3}$ (cf. Equation (19)). Petrovich et al. (2013) arrive at similar mass requirements and scalings from a model in which planets grow in mass at a prescribed rate without orbital migration or eccentricity damping. Period offsets for planets temporarily caught in resonance that might remain close to resonance after the disk dissipates would be smaller still. These considerations strongly suggest that the departure from exact resonance observed in the Kepler data occurred either during or after the protoplanetary disk disappeared.

Zoom In Zoom Out Reset image size

Retreat from resonance could occur as a result of eccentricity damping provided it were slow on the timescale of small amplitude librations about the fixed point at ϕ = 0, e = e₀. Several mechanisms come to mind. Prominent among them are tides raised by the star in the planets as suggested by Lithwick & Wu (2012) and Batygin & Morbidelli (2013), dynamical friction due to a residual particle disk, and the excitation of apsidal waves. Proposals involving tidal damping of eccentricity (Lithwick & Wu 2012; Batygin & Morbidelli 2013) have been critically examined by Lee et al. (2013) who conclude that with reasonable estimates of the tidal Love number, k, and dissipation factor, Q, tidal damping cannot account for the majority of offsets from resonance. Recent work by Baruteau & Papaloizou (2013) finds that interactions between a planet and the wake of a companion can reverse convergent migration and significantly increase the period ratio away from the resonant value. They suggest that this may help to account for the diversity of period ratios in Kepler multiple planet systems.

If a planetesimal disk contains a significant mass, interactions between the planets and the disk may lead to significant damping of orbital eccentricities. The resonant structure of the Kuiper Belt and the "late veneer" found on Earth, Moon, and Mars provide evidence for a leftover planetesimal disk in our solar system that was still present at the end of planet formation (Schlichting et al. 2012).

In some instances, interactions with a particle disk can be considered in analogy to those with a gaseous disk. If the disk were composed of small bodies, collisions might maintain the ratio of its thickness to radius at a much lower value than typical of a gaseous disk. This could lower τ_e/τ_n and perhaps provide sufficient eccentricity damping to match the observed offsets from exact resonance. However, this proposal is plagued by uncertainty. Gaps are easily opened in disks with small velocity dispersion and the absence of material corotating with the planet would eliminate eccentricity damping associated with coorbital Lindblad resonances. Coorbital Lindblad resonances dominate eccentricity damping during type I migration.

By virtue of their slow detuning away from exact resonance, apsidal waves launched at a secular resonance can propagate in low optical depth disks where waves associated with other Lindblad resonances cannot (Goldreich & Tremaine 1978; Ward & Hahn 1998). That an apsidal wave transports negative angular momentum follows from the fact that it excites the eccentricities of disk particles while hardly affecting their energies. As a consequence, the excitation of the apsidal wave damps the planet's orbital eccentricity. There is no difficulty conjuring up reasonable scenarios in which the apsidal wave is responsible for damping a planet's orbital eccentricity by a few orders of magnitude within the timescale that the particle disk might survive destruction or accretion onto planets and the central star. However, just as for damping by dynamical friction, uncertainties abound. Foremost among them are the mass and lifetime of the particle disk, the particle size distribution, and the radial dependence of the apsidal precession rate. We leave these issues for a future investigation.

Our investigation is predicated on the assumption that the orbits of the Kepler planets were largely determined by processes that operated within their protoplanetary disks and that they have undergone little modification since it disappeared. That τ_e/τ_n, which depends on temperature but is independent of disk surface density, plays a central role in the outcome of evolution in resonance is consistent with, but far from, proof that this assumption is valid. Given this starting point, we argue that the small fraction of Kepler pairs found close to mean motion resonance is compatible with standard estimates for the rates of orbital migration and eccentricity damping due to planet–disk interactions. Figure 4 shows that capture into first order mean motion resonance during convergent migration must have been a common occurrence given the observed orbital parameters and estimated planet masses. However, as the result of eccentricity damping, permanent capture in resonance was rare. Most captures were only temporary with durations of the order of τ_e which is short in comparison to the timescale τ_n for migration between neighboring resonances. The temporary nature of resonance capture is due to the overstability of librations of the resonant argument about the equilibrium given by Equation (24). The paucity of resonances among Kepler pairs should therefore not be taken as evidence for in situ planet formation or the disruptive effects of disk turbulence (Rein 2012).

The resonance model and our analytic solutions presented here only include the lowest order term in eccentricity. The assumption of small eccentricities seems to be justified for the majority of Kepler systems given that estimates for typical disk parameters yield equilibrium eccentricities smaller than 0.1.

The overstability criterion described in Section 2.4.1 is the major technical accomplishment of our paper. Permanent resonance capture is possible if equilibrium occurs at large enough separation from resonance so that a separatrix is not present. As indicated by Figures 10 and 11, only the most massive Kepler planets are likely to be permanently captured in resonance. For the others, overstable librations lead to passage through resonance and to the resumption of unimpeded convergent migration. In this context, the drop off in the number of pairs with period ratios below 1.3 seen in Figure 1 presents a puzzle, but it may, at least partly, be due to chaos caused by resonance overlap (Wisdom 1980; Deck et al. 2013).

H.S. gratefully acknowledges support from the Wade Fund. We thank Scott Tremaine, Jing Luan, and Glen Stewart for helpful comments that led to an improved manuscript.

Equilibrium eccentricities and conditions for overstable librations were derived in the body of our paper under the assumption that eccentricity is damped at constant angular momentum such that p = 3 in Equation (2). Here, we relax this assumption and provide expressions for the equilibrium eccentricity and conditions for overstable librations appropriate for a general value of p > 0. Analogous expressions to the ones derived in Section 2.4 are as follows. Equation (24) for the equilibrium eccentricity generalizes to

$$\begin{equation} e_{\rm eq} = \left(\frac {\tau_e} {(3j+p)\tau_n}\right)^{1/2}. \end{equation} \tag{ A1 }$$

Equation (29) for the growth rate of overstability is replaced by

$$\begin{equation} s = \frac{1}{\tau_e}\left(\frac{n}{\omega_{\rm eq}}\right)^2\left(jp\beta\mu^{\prime}e_{\rm eq}-\left(\frac{\beta\mu^\prime}{e_{\rm eq}}\right)^2\right). \end{equation} \tag{ A2 }$$

Thus overstable librations require p > 0. Equations (A1) and (A2) imply that a planet is permanently trapped in resonance and its librations damped provided

$$\begin{equation} \mu^\prime>\frac{jp}{(3j+p)^{3/2}\beta}\left(\frac{\tau_e}{\tau_n}\right)^{3/2}. \end{equation} \tag{ A3 }$$

For

$$\begin{eqnarray} &&\frac{j^2p}{8(3j+p)^{3/2}\beta}\left(\frac{\tau_e}{\tau_n}\right)^{3/2}<\mu^\prime<\frac{jp}{(3j+p)^{3/2}\beta}\left(\frac{\tau_e}{\tau_n}\right)^{3/2},\nonumber\\ \end{eqnarray} \tag{ A4 }$$

the planet is permanently caught in resonance and its libration amplitude saturates at a finite value. Lastly, for

$$\begin{equation} \mu^\prime<\frac{j^2p}{8(3j+p)^{3/2}\beta}\left(\frac{\tau_e}{\tau_n}\right)^{3/2}, \end{equation} \tag{ A5 }$$

the planet is caught in resonance but then escapes on timescale τ_e.