<figure><img src="image_1.png" /></figure>

# Migration of planets into and out of mean motion resonances in protoplanetary discs: overstability of capture and non-linear eccentricity damping

Wenrui Xu,<sup>1★</sup> Dong Lai<sup>2</sup> and Alessandro Morbidelli<sup>3</sup>

Accepted 2018 August 31. Received 2018 August 30; in original form 2018 May 18

## ABSTRACT

A number of multiplanet systems are observed to contain planets very close to mean motion resonances (MMRs), although there is no significant pile-up of precise resonance pairs. We present theoretical and numerical studies on the outcome of capture into first-order MMRs using a parametrized planet migration model that takes into account non-linear eccentricity damping due to planet–disc interaction. This parametrization is based on numerical hydrodynamical simulations and is more realistic than the simple linear parametrization widely used in previous analytic studies. We find that non-linear eccentricity damping can significantly influence the stability and outcome of resonance capture. In particular, the equilibrium eccentricity of the planet captured into MMRs becomes larger, and the captured MMR state tends to be more stable compared to the prediction based on the simple migration model. In addition, when the migration is sufficiently fast or/and the planet mass ratio is sufficiently small, we observe a novel phenomenon of eccentricity overshoot, where the planet’s eccentricity becomes very large before settling down to the lower equilibrium value. This can lead to the ejection of the smaller planet if its eccentricity becomes too large during the overshoot. This may help explain the intra-system mass uniformity observed in compact multiplanet systems and the lack of a low-mass planet companion of hot Jupiters when compared to warm Jupiters.

**Key words:** methods: analytical – celestial mechanics – planets and satellites: dynamical evolution and stability – planets and satellites: formation.

## 1 INTRODUCTION

The *Kepler* mission has discovered thousands of exoplanets, many of which are in multiplanet systems (Batalha et al. 2013; Coughlin et al. 2016). The period ratio distribution of the *Kepler* planets shows a significant excess of planet pairs with period ratio near (predominantly first-order) mean motion resonances (MMRs) (Fabrycky et al. 2014). This excess of planets near (or in) MMRs, together with the discovery of several resonant chain systems, such as Kepler-223 (Mills et al. 2016) and TRAPPIST-1 (Luger et al. 2017), suggests that resonance capture during disc-driven migration can be common. However, the MMR capture rate predicted using a relatively ‘clean’ migration model is much higher than the observed occurrence rate of MMRs. This discrepancy is often explained by the disruption of MMRs by physical processes after the resonance capture, including instability of the captured state during disc-driven

migration (Goldreich & Schlichting 2014; Deck & Batygin 2015; Delisle, Correia & Laskar 2015; Xu & Lai 2017), tidal dissipation in planets (Lithwick & Wu 2012; Batygin & Morbidelli 2013; Delisle, Laskar & Correia 2014), late-time dynamical instability (Pu & Wu 2015; Izidoro et al. 2017), outward (divergent) migration due to planetesimal scattering (Chatterjee & Ford 2015), and disc turbulence (Batygin & Adams 2017). Regardless of whether MMRs are maintained or destroyed by any of these processes, it is important to recognize that MMRs, even if temporarily maintained, play a significant role in the early evolution of planetary systems and can profoundly shape their final architectures.

A majority of the studies on the outcome of MMR capture (such as the impact of MMR on the orbital parameters of planets and the stability of the resonance) include the effect of disc-driven migration using a simple parametrized migration model, the most commonly used being that given by Goldreich & Tremaine (1980). The choice of this parametrized migration model makes the equation of motion of the system relatively simple, which is ideal for long-term numerical integrations or analytical studies. However, this model

<footnote n="1"><sup>1</sup>Department of Astrophysical Sciences, Princeton University, Princeton NJ 08544, USA</footnote>

<footnote n="2"><sup>2</sup>Department of Astronomy, Center for Astrophysics and Planetary Science, Cornell University, Ithaca NY 14853, USA</footnote>

<footnote n="3"><sup>3</sup>Laboratoire Lagrange, Université Cote d'Azur, Observatoire de la Cote d'Azur, CNRS, CS 34229, F-06304 Nice, France</footnote>

<footnote>* E-mail: wenruix@princeton.edu</footnote>
only works well for small eccentricities ($e \lesssim H/r$, the aspect ratio of the disc). As we show in this paper, the eccentricities of the planets near the MMR can often lie in the regime where the Goldreich & Tremaine (1980) result is no longer valid. This can impact the outcome of the resonance capture. There are also a number of studies that include more realistic migration models, such as those using parametrized forcing in N-body integration (e.g. Terquem & Papaloizou 2007; Migaszewski 2015) or using self-consistent hydrodynamics (e.g. Kley et al. 2005; Papaloizou & Szuszkiewicz 2005; Crida, Sándor & Kley 2008; Zhang et al. 2014; André & Papaloizou 2016). However, these studies tend to focus on explaining the behaviours of particular systems and do not survey a sufficiently large parameter space to obtain various possible outcomes. The goal of our paper is to remedy this situation. In particular, we generalize previous analyses (Goldreich & Schlichting 2014; Deck & Batygin 2015; Delisle et al. 2015; Xu & Lai 2017) by adopting a more realistic parametrization for the migration and eccentricity damping, and examine how different model parameters affect the outcome of the MMR capture.

This paper is organized as follows. Section 2 summarizes the parametrizations for the rates of orbit decay and eccentricity damping due to planet-disc interactions. In Section 3, we consider the simple case when one of the planets is massless and study how different parametrizations can affect the outcome of MMR capture. We find that using the more realistic migration model can sometimes cause the ejection of the small planet, but otherwise tends to increase the stability of the resonance. In Section 4, we study the more realistic case when both planets have finite masses. While most of the results from Section 3 can be generalized, we also observe several new phenomena that arise only when both planets have finite masses. In addition to analytical calculations, we use 3-body integrations to validate our results. We conclude in Section 5 and discuss how our results affect the architecture of multiplanet systems.

## 2 PARAMETRIZATIONS OF THE RATES OF ORBIT DECAY AND ECCENTRICITY DAMPING

Consider a small planet undergoing type I migration in a gaseous disc. At low eccentricity, the rates of orbit decay and eccentricity damping due to planet-disc interaction are approximately given by (Goldreich & Tremaine 1980)

$$ \frac{\dot{a}}{a} = -\frac{1}{T_m} - \frac{2pe^2}{T_e}, \qquad (1) $$

$$ \frac{\dot{e}}{e} = -\frac{1}{T_e}, \qquad (2) $$

where $T_m$ and $T_e$ are independent of $e$ and $T_m \sim T_e h^{-2}$, with $h \equiv H(r)/r$ ($H$ is the disc's scale height). Analytic calculations of Chiang & Goldreich (1997) give $h \sim 0.05$ at $\sim 1$ au, with weak dependence on $a$. Here we assume $h$ is constant, for simplicity. The parameter $p$ characterizes the coupling between orbit decay and eccentricity damping; here we take $p = 1$, which corresponds to eccentricity damping that conserves angular momentum. This is the parametrized migration model used in most studies of MMR capture (Goldreich & Schlichting 2014; Deck & Batygin 2015; Delisle et al. 2015; Xu & Lai 2017).

However, this migration model is accurate only for small eccentricities, $e \lesssim h$. For larger eccentricities, hydrodynamic simulations (Cresswell et al. 2007; Cresswell & Nelson 2008) show that the

orbit decay rate and eccentricity damping rate both decrease. As an empirical fit to the numerical results, $T_m$ and $T_e$ are functions of $e/h$ given by (based on equations 11 and 13 of Cresswell & Nelson 2008)

$$ T_m = T_{m,0} \frac{1 + (e/2.25h)^{1.2} + (e/2.84h)^6}{1 - (e/2.02h)^4}, \qquad (3) $$

$$ T_e = T_{e,0} \left( 1 - 0.14 \frac{e^2}{h^2} + 0.06 \frac{e^3}{h^3} \right), \qquad (4) $$

$$ \text{with } T_{m,0} \equiv \frac{t_{\text{wave}}}{2.7 + 1.1\beta} h^{-2}, \quad T_{e,0} \equiv \frac{t_{\text{wave}}}{0.78}. \qquad (5) $$

Here we assume that the disk has a density profile $\Sigma(r) \propto r^{-\beta}$; we adopt $\beta = 0$ (i.e. a disk with uniform surface density) unless otherwise specified. The time-scale $t_{\text{wave}}$ is given by (Tanaka & Ward 2004)

$$ t_{\text{wave}} = \frac{M_*^2}{\Sigma a^2 m} h^4 \Omega^{-1}, \qquad (6) $$

with $\Omega$ being the angular velocity of the unperturbed disk.

In this paper we compare two different migration models/ parametrizations: the ‘simple’ model, with $T_m = T_{m,0}$ and $T_e = T_{e,0}$ independent of $e$, and the ‘realistic’ model, with $T_m$, $T_e$ given by equations (3) and (4). The two models produce identical $T_e$, $T_m$ for $e/h \ll 1$,<fnref n="1" /> but can give very different orbit decay and eccentricity damping rates when $e/h$ is large. This is illustrated in Fig. 1. In particular, for the realistic migration model, the eccentricity damping rate scales as $e^{-3}$ when $e/h \gg 1$.

Planet-disc interaction is a very complex problem, and is strongly influenced by the structure of the protoplanetary disc (see reviews in Baruteau et al. 2014; Nelson 2018). The migration model we adopt (based on the simulations of Cresswell et al. 2007; Cresswell & Nelson 2008) assumes a simple disc profile and includes only hydrodynamic interactions. Nevertheless, we expect the non-linear eccentricity damping (i.e. the decrease of damping strength for $e \gtrsim h$ a few $h$) to be one of the most important differences between reality and the simple migration model, and this feature is fully captured in the ‘realistic’ migration model we use.

## 3 OUTCOME OF MMR CAPTURE: MASSLESS INNER PLANET

To gain some analytical understanding of the general problem of MMR capture with planets of comparable mass, in this section we consider a simpler case: a planet with negligible mass ($m$) perturbed by an outer massive planet ($m'$) on a circular orbit near a first-order $j$: $j+1$ MMR. To this end, we take $T'_{m'}$, the orbit decay time-scale of the outer planet, to be a free parameter. This allows us to explore how the equilibrium eccentricity (of the inner planet), which is determined by the net convergent migration rate, affects the outcome of the MMR capture. In reality, both planets undergo migration. For $m \ll m'$ and type I migration, we expect $T_m$, $T_e h^{-2} \gg T'_{m'}$. The results in this section should qualitatively illustrate how the outcomes of MMR capture are affected when the realistic migration model is applied (see Section 4).

<footnote n="1">¹However, the derivatives of $T_e$, $T_m$ with respect to $e/h$ produced by the two models are still significantly different for very small $e/h$, which can cause some issue when computing the stability of the equilibrium state at small $e/h$. See more discussion in Section 3.2.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Orbit decay and eccentricity damping rates given by the simple migration model (eccentricity-independent $T_e$, $T_m$, black dashed curves) and the realistic migration model (eccentricity-dependent $T_e$, $T_m$, red solid curves). The difference becomes prominent when $e/h \gtrsim 3$.</figcaption></figure>

In this section, we also assume that the planet-disc interaction is weak, so $T_m$, $T_e$, $T'_m$ (note that $T'_e$ is irrelevant since the outer planet's orbit is always circular) are much greater than the time-scale of libration, $T_{res}$, given by

$$T_{\text{res}} \approx 0.8 j^{-4/3} (\mu')^{-2/3} \frac{2\pi}{n}, \quad (7)$$

where $n$ is the mean motion of the inner planet and $\mu' = m'/M_\star$ the mass ratio between the outer planet and the star. (For the exact definition of $T_{res}$, see equation B6 in appendix B of Xu & Lai 2017.) For simplicity, we assume that $T_{e,0}$, $T_{m,0}$, $T'_m$ remain constant (i.e. their variations due to the evolution of the planets' semi-major axes are ignored).

## 3.1 Existence of equilibrium

We first study the eccentricity at the equilibrium state (and whether such an equilibrium state exists). Near a first-order $j: j + 1$ MMR, the resonant motion conserves

$$\alpha_0 \equiv \alpha(1 + je^2), \quad (8)$$

where $\alpha = a/a' < 1$ is the semi-major axis ratio. When the system undergoes convergent migration, the inner planet can be captured into the resonance. It reaches an equilibrium state when $d\alpha_0/dt = 0$, which corresponds to

$$T_{m,\text{eff}}^{-1} = 2(j + 1)e^2 T_e^{-1}. \quad (9)$$

Here $T_{m,\text{eff}}$ (which may depend on $e$) is the effective convergent migration rate given by $T_{m,\text{eff}}^{-1} = T_m^{-1} - T_m^{-1}$. Note that when the

outer planet is much more massive it should migrate much faster than the inner planet, so $T_{m,\text{eff}} \approx T'_m$.

For the simple migration model with constant $T_e = T_{e,0}$ and $T_m = T_{m,0}$, the equilibrium always exists, with the corresponding eccentricity given by (Goldreich & Schlichting 2014)

$$e_{\text{eq},0} = \sqrt{\frac{T_{e,0}}{2(j+1)T_{m,\text{eff},0}}} \approx \sqrt{\frac{T_{e,0}}{2(j+1)T'_m}}, \quad (10)$$

where $T_{m,\text{eff},0}^{-1} \equiv T_m^{-1} - T_{m,0}^{-1}$.

However, for the realistic migration model with eccentricity-dependent $T_e$ and $T_m$, the right-hand side of equation (9) has a finite maximum value because for $e \gtrsim$ a few $h$, $e^2 T_e^{-1} \propto e^{-1}$ decreases as $e$ increases. Therefore, the equilibrium may not exist when the outer planet's migration is too fast. The maximum value of $e^2 T_e^{-1}$ occurs at $e \approx 3h$; thus the equilibrium ceases to exist when

$$T_m^{-1} \gtrsim 2(j+1)T_{e,0}^{-1}(3h)^2. \quad (11)$$

Fig. 2 gives an example of the evolution of the system when the equilibrium of resonance capture does not exist. This result is obtained by doing 3-body integrations using the MERCURY code (Chambers 1999), with $a' = 1$ au and $T_{e,0} = 10T_{\text{res}}$ [$T_{\text{res}}$ is given by equation (7)]. Since the equilibrium does not exist, the planet's eccentricity grows unboundedly as the system goes deeper into the resonance.

## 3.2 Stability of capture

The migration model can also affect the stability of the captured (equilibrium) state.

For the simple migration model, the stability of the equilibrium state has been studied by Goldreich & Schlichting (2014). Under the assumption that planet-disc interaction is weak, the behaviour of the system depends only on the ratio $\mu'/e_{\text{eq},0}^3$, where $\mu' = m'/M_\star$ and $e_{\text{eq},0}$ is the previously defined equilibrium eccentricity. The equilibrium is stable when the outer planet is sufficiently massive (with $\mu' \gtrsim e_{\text{eq},0}^3$); in this case the resonant angle librates with small amplitude. For $\mu' \sim e_{\text{eq},0}^3$, the libration amplitude saturates at a finite value, and the system stays in resonance. For a less massive outer planet (with $\mu' \lesssim e_{\text{eq},0}^3$), the equilibrium state is overstable (i.e. the amplitude of libration increases with time) and the system eventually escapes from resonance.

For the realistic migration model, however, the stability of the equilibrium state depends on not only $\mu'/e_{\text{eq},0}^3$ but also $e_{\text{eq},0}/h$; the latter parameter characterizes how significantly the system is affected by including the eccentricity dependence in the migration model.

Fig. 3 plots the regimes of different behaviours in the $e_{\text{eq},0}/h - \mu'/e_{\text{eq},0}^3$ parameter space for a 2:3 MMR when the realistic migration model is applied. We integrate the equation of motion derived from the resonance Hamiltonian (see e.g. appendix B of Xu & Lai 2017), and include the dissipative terms associated with migration and eccentricity damping.<fnref n="2" /> For each data point, we choose $T_{e,0}/T_{m,0}$ to produce the desired $e_{\text{eq},0}/h$, then choose $\mu'$ to produce the desired $\mu'/e_{\text{eq},0}^3$. We integrate the system until it passes the resonance (the period ratio becomes significantly lower than the resonant value), its eccentricity blows up ($e/h$ keeps increasing and

<footnote n="2">²Direct integration of the equation of motion is necessary because the outcome when the equilibrium state is overstable (whether the libration saturates at a finite amplitude, or the system eventually escapes the resonance) cannot be obtained from linear stability analysis of the equilibrium state.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Evolution of a system near 2:3 MMR with $h = 0.05$, $T'_m = 10T_{e,0}$, and the outer planet mass $\mu' = m'/M_* = 10^{-3}$. Left-hand panel: period ratio. Centre panel: eccentricity of the inner planet. The black dashed line shows $e = 3h$ for reference. Right-hand panel: resonant angle. The equilibrium does not exist and the planet's eccentricity grows unboundedly as the system goes deeper into the resonance. This increasing eccentricity should ultimately cause ejection or collision of the planet.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Regimes of different behaviours in the $e_{eq,0}/h - \mu'/e_{eq,0}^3$ parameter space for a 2:3 MMR. The two dashed lines mark the analytical estimates for the boundary between libration with small amplitude, libration with finite amplitude and escape given by Goldreich & Schlichting (2014) for the simple migration model.</figcaption></figure>

goes beyond 20), the libration amplitude stops changing (change by <5 per cent for the last half of the simulation), or the libration amplitude decreases to small values (<0.1 rad). This divides the parameter space into four possible outcomes/behaviours:

(i) When $e_{eq,0}$ is larger than $2.4h$, the equilibrium state of resonance capture does not exist because the eccentricity damping is too weak to balance the eccentricity excitation due to resonant interaction, and the planet's eccentricity $e$ grows unboundedly until the system becomes unstable (red diamonds in Fig. 3).

(ii)–(iv) When $e_{eq,0}$ is small enough to allow the existence of an equilibrium state, this equilibrium can be stable or overstable. When it is stable, the system exhibits small libration around the equilibrium state with the libration amplitude converging to zero (purple circles in Fig. 3). When it is overstable, the system can either end up in a stable state with a finite libration amplitude (green

triangles) or exit the resonance with damped eccentricity (orange squares). Only these three behaviours are possible in the simple migration model.

Although Fig. 3 refers to the 2:3 MMR, we find that the results for other first-order MMRs are qualitatively similar.

Three-body simulations (see below) show that the results obtained from the resonant Hamiltonian in Fig. 3 are qualitatively correct, with tolerable error for the boundaries between different behaviours. Note that the boundaries between the last three behaviours (stable libration with finite and small amplitude, and escape) depend sensitively on the migration model, since the stability of the equilibrium is affected by the derivatives of $T_e$ and $T_m$.<fnref n="3" />

Fig. 3 differs from the result based on the simple migration model (e.g. Goldreich & Schlichting 2014) in several aspects. First, as noted above, there exists a new regime where the planet's eccentricity can grow unboundedly because of the decrease of the eccentricity damping rate for $e \gtrsim h$. Secondly, near the boundary of this 'eccentricity blowing up' regime ($2.3 \lesssim e_{eq,0}/h \lesssim 2.45$), the stable finite-amplitude libration regime occupies a large parameter space; in particular, the system can stay in resonance with a finite-amplitude libration even when $\mu'/e_{eq,0}^3$ is as small as 0.6 (by contrast, the simple migration model would predict the system escape from the resonance due to overstability). Thirdly, the boundaries between the different regimes, even at $e_{eq,0}/h \lesssim 2$ (for which $T_m$ and $T_e$ deviate little from the simple model), are significantly distorted due to the use of the more realistic migration model, showing that these boundaries are indeed sensitive to the migration model (and disc parameters). Note that for low eccentricity ($e_{eq,0}/h \lesssim 1$), the derivative of $T_m$ with respect to $e/h$ for the realistic migration model converges to the simple migration model very slowly [because of the $(e/h)^{1.2}$ term in the fitting formula of $T_m$ (equation 3)]. This is visible in the upper panel of Fig. 1: Even at $e/h$ very close to 0, the slopes of $\dot{a}/a$ for the two models are significantly different. We believe that this difference arises because Cresswell & Nelson (2008) slightly overfitted their result for small $e/h$, and in reality the two models should converge faster. Since the stability of the

<footnote n="3">³One can see this by considering how the stability of the equilibrium point is calculated. The stability is determined by the eigenvalues of a matrix with entries of the form $\partial_\lambda(dy/dt)$, where x, y can be either a or e. These entries depend not only on the values of $T_e$ and $T_m$ but also on their derivatives with respect to a or e.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Evolution of a system captured into the 2:3 MMR with $e_{eq,0}/h = 2.3$ and $\mu'/e_{eq,0}^3 = 1$. The four panels plot the inner planet's resonant angle, eccentricity, the two planets' period ratio, and the inner planet's trajectory in the $ecos\theta - esin\theta$ phase space after reaching equilibrium. The system ends up in a stable state with a finite libration amplitude.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Same as Fig. 4, but with $e_{eq,0}/h = 2.5$ and $\mu'/e_{eq,0}^3 = 1$. Equilibrium no longer exists and the planet's eccentricity grows unboundedly. The system eventually becomes unstable and the inner planet gets ejected shortly after $e$ reaches ~0.8.</figcaption></figure>

equilibrium is sensitive to the derivatives of $T_e$ and $T_m$, we do not expect our model using equation (3) to recover the analytical result of Goldreich & Schlichting (2014) for small $e_{eq}$, and the result for $e_{eq}/h < 1$ is not shown in Fig. 3. In future studies, using a fitting formula that produces the correct behaviour of both $T_e$, $T_m$ and their derivatives for small $e/h$ can help improve the accuracy of the stability prediction for $e/h \lesssim 1$.

Figs 4–7 show the behaviour of the system in each regime depicted in Fig. 3. Similar to Fig. 2, these results are obtained by doing 3-body integrations using the MERCURY code, with $h = 0.025$, $a' = 1$ au, and $T_{e,0} = 10T_{\text{res}}$. The other parameters of the system can be solved to match the given $e_{eq,0}/h$ and $\mu'/e_{eq,0}^3$ values. In practice,

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Same as Fig. 4, but with $e_{eq,0}/h = 2$ and $\mu'/e_{eq,0}^3 = 1$. The planet escapes from the resonance due to the overstability of the equilibrium state. Note that soon after the planet exits the 2:3 MMR, it gets captured into a 3:4 MMR.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Same as Fig. 4, but with $e_{eq,0}/h = 2$ and $\mu'/e_{eq,0}^3 = 3.5$. The equilibrium state is stable and the libration of the resonant angle stays small.</figcaption></figure>

to avoid having the planets migrate too far inward during the integration (which will make it necessary to choose a much smaller time-step to account for the planet's short orbital period), we fix the outer planet and let the inner planet's semi-major axis increase at the rate $\dot{a}/a = -1/T_m - 2e^2/T_e + 1/T'_m$ – note that this parametrized treatment is necessary because the overstability time-scale of the equilibrium state can be $\gtrsim 10T'_m$ in many cases.

## 4 OUTCOME OF MMR CAPTURE: TWO MASSIVE PLANETS

To apply our results to realistic systems, it is important to study the case where both planets have finite masses. As we will show in this section, the perturbation on the more massive planet from the smaller planet can qualitatively affect the outcome of resonance
capture even when the mass ratio is very small. We will also discuss the effect of a strong eccentricity damping rate and non-adiabatic evolution due to fast migration.

## 4.1 Existence and location of equilibrium

Consider two planets near a $j: j + 1$ MMR, with both planets having finite masses. Let the inner (outer) planet have mass $m_1$ ($m_2$) and semi-major axis $a_1$ ($a_2$).<fnref n="4" /> The Hamiltonian of the system, to first order in eccentricity and with all non-resonant terms averaged out, is given by

$$H = -\frac{GM_*m_1}{2a_1} - \frac{GM_*m_2}{2a_2} - \frac{Gm_1m_2}{a_2} (f_{j+1,27} e_1 \cos \theta_1 + \tilde{f}_{j+1,31} e_2 \cos \theta_2). \quad (12)$$

Here $\alpha = a_1/a_2$, and $f_{m,n}$ are functions of $\alpha$ (evaluated at $\alpha_0 \equiv [j/(j+1)]^{2/3}$) given in appendix B of Murray & Dermott (1999), with $\tilde{f}_{j+1,31} \equiv f_{j+1,31} - \delta_{j,1}2\alpha_0$. The interaction between the two planets conserves the total angular momentum

$$\mathcal{L} \equiv \Lambda_1 \sqrt{1 - e_1^2} + \Lambda_2 \sqrt{1 - e_2^2}, \quad (13)$$

where $\Lambda_i = m_i \sqrt{GM_*a_i}$. The Hamiltonian (12) also admits a second constant of motion (Michtchenko & Ferraz-Mello 2001),

$$\mathcal{K} \equiv \frac{j+1}{j} \Lambda_1 + \Lambda_2. \quad (14)$$

Combining the two ($\mathcal{L}$ and $\mathcal{K}$) produces a conserved quantity $\eta$, given by

$$\begin{aligned} \eta &\equiv -2(q\alpha_0^{-1} + 1) \left[ \frac{\mathcal{L}}{\mathcal{K}} - \left( \frac{\mathcal{L}}{\mathcal{K}} \right)_{\alpha=\alpha_0, e_i=0} \right] \\ &\approx \frac{q}{j\alpha_0^{1/2}(q\alpha_0^{-1} + 1)} (\alpha - \alpha_0) + \alpha_0^{1/2} q e_1^2 + e_2^2, \end{aligned} \quad (15)$$

where $q = m_1/m_2$ is the mass ratio, and $\alpha_0 = [j/(j+1)]^{2/3}$ is the semi-major axis ratio at resonance. In the second line of equation (15), we have expanded the result to the lowest order in $(\alpha - \alpha_0)$ and $e_2^2$. The parameter $\eta$ characterizes how deep the system is inside the resonance when captured: For larger $\eta$, the system is deeper inside the resonance, and the fixed point (libration centre) of the system corresponds to larger eccentricities.

Consider the evolution of $e_i$, $\varpi_1 - \varpi_2$, and $\eta$. At the equilibrium state, $\varpi_1 - \varpi_2$ is constant because the resonant angles $\theta_i \equiv (j + 1)\lambda_2 - j\lambda_1 - \varpi_i$ are constant; $\eta$, which is a function of $\alpha$ and $e_i$, should also be constant because $\alpha$ and $e_i$ are constant. Therefore, $e_i$ and $\theta_i$ at the equilibrium state can be solved from the following equations<fnref n="5" />:

$$\frac{de_1}{dt} = -\mu_2 n_1 \alpha_0 f_{j+1,27} \sin \theta_1 - \frac{e_1}{T_{e,1}} = 0, \quad (16)$$

$$\frac{de_2}{dt} = -\mu_1 n_2 \tilde{f}_{j+1,31} \sin \theta_2 - \frac{e_2}{T_{e,2}} = 0, \quad (17)$$

$$\frac{d(\varpi_1 - \varpi_2)}{dt} \propto f_{j+1,27} e_2 \cos \theta_1 - \tilde{f}_{j+1,37} e_1 q \alpha_0^{1/2} \cos \theta_2 = 0, \quad (18)$$

$$\begin{aligned} \frac{d\eta}{dt} = & \frac{q\alpha_0^{1/2}}{j(q\alpha_0^{-1} + 1)} \left( \frac{1}{T_{m,2}} - \frac{1}{T_{m,1}} + \frac{2e_1^2}{T_{e,1}} - \frac{2e_2^2}{T_{e,2}} \right) \\ & - q\alpha_0^{1/2} \frac{2e_1^2}{T_{e,1}} - \frac{2e_2^2}{T_{e,2}} = 0. \end{aligned} \quad (19)$$

Note that since $\eta$ is conserved in the absence of dissipation (planet-disc interaction), equation (19) only includes contributions from planet-disc interactions. Equation (19) can be interpreted physically as convergent migration tending to push the system deeper into resonance (i.e. increasing $\eta$ and eccentricities) while eccentricity damping (from planet-disc interaction) counters the effect of migration. Equilibrium is reached ($\eta$ ceases to evolve) when migration and eccentricity damping balance each other.

A caveat of the above method is that the Hamiltonian (equation 12) and the equations for the equilibrium (equation 16–19) only include the lowest-order terms in eccentricities; i.e., we have effectively assumed $e_1, e_2 \ll 1$. For realistic systems, when $e_1$ becomes large, higher-order secular couplings may affect the result. In fact, as we will see in Section 4.2, the results can be quite different from the predictions of this Hamiltonian model when $e_1$ approaches order unity.

### 4.1.1 Weak eccentricity damping

First consider the case when the eccentricity damping is weak, i.e. $\mu_2 n_1 \gg e_1/T_{e,1}$ and $\mu_1 n_2 \gg e_2/T_{e,2}$. In this case, $|\cos \theta_i| \approx 1$, and equation (18) gives $e_1/e_2$ at the equilibrium. Note that $e_1/e_2 \sim q^{-1} = m_2/m_1$ and is independent of $T_{e,i}$ and $T_{m,i}$. With $e_1/e_2$ known, and with $T_{e,i}$, $T_{m,i}$ as a function of $e_i$ (see Section 2), we can solve equation (19) to obtain $e_i$ at the equilibrium.

Fig. 8 shows the equilibrium eccentricities of the two planets calculated using the above method. For the simple migration model, equations (18) and (19) give $e_{1,eq} \sim q^{-1/2}h$ and $e_{2,eq} \sim q^{1/2}h$. At equilibrium, the $e_1^2$ terms and the $e_2^2$ terms in equation (19) are comparable when $q \lesssim 1$: The eccentricity terms in the first line of (19) are comparable to or smaller than the corresponding eccentricity terms in the second line when $q \lesssim 1$, and $qe_1^2/T_{e,1} \sim e_2^2/T_{e,1}$ given that $T_{e,1}/T_{e,2} \sim q^{-1}$.

For the realistic migration model, the result is similar to that of the simple migration model when $q$ is relatively large ($q \gtrsim 0.15$, 0.04, and 0.03 for the 1:2, 2:3, and 3:4 MMR, respectively). When $q$ is smaller, however, $e_{1,eq}$ exceeds $\sim 3h$ and the damping rate $T_{e,1}^{-1}$ is reduced. Therefore, the equilibrium eccentricities of the planets are larger compared to those of the simple migration model in order to satisfy equation (19). For 1:2 MMR, the equilibrium eccentricity exhibits a step-like feature in the left-hand panel of Fig. 8. This is due to a bifurcation where as $q$ decreases, a pair of solutions to $d\eta/dt = 0$ ceases to exist (see Fig. 9). Such bifurcation does not happen for 2:3 and 3:4 MMR, since $d\eta/dt$ is monotonic for them.

The equilibrium always exists when both planets have finite masses, although it may correspond to $e_{1,eq} \gtrsim 1$, which implies that the smaller planet can be ejected due to instability before reaching the equilibrium. This is very different from the ‘massless inner planet’ case considered in Section 3, where the equilibrium state may not exist. The reason of such a difference is that while the eccentricity of the smaller planet can exceed $3h$, the eccentricity of the more massive planet always remains well below $3h$, so the eccentricity damping from the more massive planet is able to balance the migration, ensuring the existence of an equilibrium state.

<footnote n="4">⁴The notation is different from Section 3 in order to emphasize the fact that both planets have finite masses.</footnote>

<footnote n="5">⁵Another method is to directly solve for the equilibrium state by linking the evolution of all quantities to that of $R \equiv a_2/a_1$, and imposing that $\dot{R} = 0$, all the while considering the torques exerted by the disc on the planets that result from the migration model (Pichierri et al. 2018, in preparation). Our approach makes it easier to analyse how the equilibrium eccentricities are affected by using different migration models.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Equilibrium eccentricities versus mass ratio $q = m_1/m_2$ of the two planets captured into a first-order MMR in a disc with uniform surface density ($\beta = 0$ in $\Sigma \propto r^{-\beta}$). Note that the $e_i/h$ value depends only on $q$ and is independent of the total mass of the planets and $h$. Only the $q < 1$ region is plotted since convergent migration requires $q \lesssim 1$. The blue (red) curve shows the eccentricity of the inner (outer) planet. The solid (dashed) curve shows the eccentricity for the realistic (simple) migration model. These results are calculated under the assumption that eccentricity damping is weak.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** $d\eta/dt$ for 1:2 MMR as a function of $e_1/h$ at three different $q$ values before, at, and after the bifurcation that creates the step-like feature in Fig. 8. For each $q$, the equilibrium is marked by a dot. The parameters and assumptions are the same as those in Fig. 8, and $e_2$ is chosen such that equation (18) is satisfied.</figcaption></figure>

The equilibrium eccentricities also depend on the density profile of the disc, which is characterized by the parameter $\beta$ [assuming that the disc has $\Sigma(r) \propto r^{-\beta}$; note that we adopt $\beta = 0$ everywhere else in this paper]. Fig. 10 shows that the equilibrium eccentricities of the planets depend weakly on $\beta$.

### 4.1.2 Effect of strong eccentricity damping

When $q$ is small, the resonant perturbation from the inner planet is no longer much stronger than the eccentricity damping of the outer planet, and the second term in equation (17) can no longer be ignored. In this regime, the equilibrium eccentricities can be significantly affected when the realistic migration model is applied.

For a sufficiently small $q$ (which gives a large $e_1/h$), the terms proportional to $e_1^2$ in equation (19) are negligible, so $e_{2,eq} \sim q^{1/2}h$ can be determined directly from equation (19) and are independent of the strength of the eccentricity damping. Meanwhile, equations (17) and (18) suggest that for smaller $\mu_1$ (or larger $1/T_{e,2}$), $|\sin\theta_2|$ increases, $|\cos\theta_2|$ decreases, and $e_{1,eq}/e_{2,eq}$ increases. In particular, when $\mu_1$ is sufficiently small (i.e. $\mu_1 n_2 \sim e_{2,eq}/T_{e,2} \sim q^{1/2}h/T_{e,2}$), $|\cos\theta_2| \to 0$ and $e_{1,eq}/e_{2,eq}$ diverge. Since $e_{2,eq}$ is finite, this means

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Equilibrium eccentricities of the two planets captured into a 1:2 MMR in discs with different $\beta$. The blue (red) curves show the eccentricities of the inner (outer) planet. The solid, dashed, and dotted curves are the eccentricities for $\beta = 0, 1, \text{and } -1$, respectively. The results are similar for different values of $\beta$.</figcaption></figure>

that $e_{1,eq}$ diverges (i.e. ejection or collision of the smaller planet should happen before the equilibrium is reached).

Fig. 11 (based on numerical calculations of the equilibrium eccentricities) demonstrates this effect. For given $\mu_1 + \mu_2$ and $n_1, n_2$, the critical $q$ at which $e_{1,eq}$ diverges is related to the characteristic eccentricity damping rate by $q_{crit} \propto t_{w,0}^{-2}$, where $t_{w,0}$ is a time-scale characterizing the migration and eccentricity damping defined as $t_{wave}$ (see equation 6) evaluated at $m = m_1 + m_2$ and $a = a_2$. Note that $t_{w,0}$ is determined by the disc parameters, and is comparable to $T_{e,0}$ of the larger planet. This scaling for $q_{crit}$ can be explained as follows: The eccentricity of the smaller planet $e_{1,eq}$ diverges when $\cos\theta_2 \to 0$ according to equation (18). When $\cos\theta_2 \to 0$ (and $\sin\theta_2 \to 1$), equation (17), together with the fact that $e_{2,eq} \sim q^{1/2}h$, gives (assuming $q$ is small)

$$ \mu_1 n_2 \sim e_{2,eq}/T_{e,2} \Rightarrow q(\mu_1 + \mu_2)n_2 \sim q^{1/2}h/t_{w,0}, \quad (20) $$ which then gives $q_{crit} \propto t_{w,0}^{-2}$.
<figure><img src="image_12.png" /></figure>

**Figure 11.** Equilibrium eccentricities of the two planets with $\mu_1 + \mu_2 = 10^{-3}$, $M_* = 1M_\odot$, and $a_2 = 1$ au captured into the 1:2 MMR in a disc with uniform surface density, for different strengths of eccentricity damping. The eccentricity damping rate is characterized by $t_{w,0}$, defined as $t_{wave}$ (see equations 3–6) evaluated at $m = m_1 + m_2$ and $a = a_2$. The blue (red) curves show the eccentricity of the inner (outer) planet. Different line styles correspond to different $t_{w,0}$, with the solid curves corresponding to very slow migration ($t_{w,0} \to \infty$).

## 4.2 Three-body simulations: effects of non-linear eccentricities and non-adiabatic evolution

We now use 3-body simulations to check our semi-analytical results obtained in the previous subsection. This is necessary since the Hamiltonian (equation 12) assumes that the eccentricities are small, which may lead to non-trivial errors when $e_1$ attains large values. In addition, it is useful to use 3-body integrations to investigate at which point and for what reason(s) the inner planet becomes dynamically unstable at high eccentricities.

Figs 12 and 13 compare the 3-body integration results for the 1:2 and 2:3 MMRs using MERCURY with our analytical results. Forcing due to planet-disc interaction is implemented as described in Cresswell & Nelson (2008) to agree with equations (3) and (4). Overall, the 3-body integration results agree with our analytical results, showing the general trend that the equilibrium eccentricities increase (compared to the simple migration model) for small $q$. However, there are several important effects that the semi-analytical linear theory fails to capture, and we discuss these effects below.

### 4.2.1 Effect of high-order coupling at large $e_1$

Fig. 12 and the $t_{w,0} = 100$ yr curves<fnref n="6" /> in Fig. 13 show that $e_{1,eq}$ is smaller than the analytical prediction when $e_{1,eq} \sim 1$. This is likely due to the higher-order secular coupling between the planets; such coupling prevents $e_1$ from reaching unity while $e_2$ remains finite. As a result, the divergence of $e_1$ due to finite eccentricity damping (discussed in Section 4.1.2) does not occur in real systems. (The ejection of the inner planet for small $q$ depicted in Figs 12 and 13 is due to eccentricity overshoot, a phenomenon we will discuss next.)

<figure><img src="image_13.png" /></figure>

**Figure 12.** Equilibrium eccentricities of two planets with $\mu_1 + \mu_2 = 10^{-3}$, $M_* = 1M_\odot$, and initial $a_2 = 1$ au captured into the 1:2 MMR in a disc with uniform surface density and $h = 0.025$. The black curves show our analytical results, as given in Fig. 8. The 3-body integration results for $t_{w,0} = 10$ yr (100 yr) are shown in crosses (saltires). For $t_{w,0} = 10$ yr, $e_1$ ($e_2$) is marked by the blue (red) curve; for $t_{w,0} = 100$ yr, $e_1$ ($e_2$) is marked by the green (yellow) curve. When the inner planet is ejected (or collides with the other planet or the star), we set $e_1 = \infty$ and $e_2 = 0$.

<figure><img src="image_14.png" /></figure>

**Figure 13.** Same as Fig. 12, but for the 2:3 MMR. The behaviour of the system is slightly different; see the text for more discussion.

### 4.2.2 Effect of non-adiabatic evolution: eccentricity overshoot

For sufficiently slow migration, the evolution of the system is adiabatic (i.e. the evolution of $\eta$, the 'resonance depth' parameter, is sufficiently slow so that the system stays close to the libration centre as the libration centre moves in the phase space) and the eccentricities of both planets should slowly increase until they reach the equilibrium values. In this case, the equilibrium eccentricities are the maximum eccentricities that the planets can reach. However, when $q$ is small or when migration is fast (i.e. $t_{w,0}$ is small), the growth of $e_2$ is too slow, and the initial evolution of $e_1$ is similar to the restricted problem studied in Section 3: Due to the inefficient eccentricity damping, $\eta$ and $e_1$ both keep increasing, and $e_1$ can easily grow beyond the equilibrium value. The growth of $e_1$ stops only when it becomes so large that the secular interaction between

<footnote n="6">⁶The other curves in Fig. 13 will be discussed in Section 4.2.3.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** An example where the eccentricity $e_1$ overshoots to a large value before reaching equilibrium, for the system in Fig. 12 with $t_{w,0} = 10$ yr and $q = 0.0125$. The different panels show the two planets' eccentricities ($e_1$, $e_2$), the period ratio, and the resonant angles ($\theta_1 = 2\lambda_2 - \lambda_1 - \omega_1$, $\theta_2 = 2\lambda_2 - \lambda_1 - \omega_2$). In this example, $e_1$ reaches a maximum value of 0.935 before decreasing to the equilibrium value. Note that for a slightly smaller $q$, the inner planet will have $e_1 \to 1$ and become unstable during the overshoot.</figcaption></figure>

the planets forces $e_2$ to increase. Since eccentricity damping of $e_2$ is still efficient, this stops the system from going deeper into the resonance (i.e. stops $\eta$ from further increasing). Eventually, the system will reach equilibrium, provided that the smaller planet has not become dynamically unstable during the high-$e_1$ phase.

Fig. 14 shows an example. Before the system reaches equilibrium, the eccentricity $e_1$ first overshoots to a very large value, then decreases back to the equilibrium value. When $q$ is smaller (or when the migration is faster), the inner planet will be ejected because it reaches $e_1 \to 1$ during this overshooting phase. This is the reason for the ejection of the smaller planet at low $q$ in Figs 12 and 13.

The condition of significant eccentricity overshoot (which may lead to planet ejection) can be estimated as follows. In order for significant eccentricity overshoot to happen, the time it takes for $e_1$ to reach $O(1)$ when $e_2$ is small compared to the equilibrium value should be no larger than the time it takes for $e_2$ to reach the equilibrium value. The first time-scale is of the order of $T_{m,2}$, since for $e_1$ to reach $O(1)$ requires the angular momentum of the system to change by $O(1)$. Using equation (17) (assuming $|\sin\theta_2| \sim 1$ and ignoring the second term) and the fact that $e_{2,eq} \sim q^{1/2}h$ (see Section 4.1.1), the second time-scale is of the order of $P_2\mu_1^{-1}q^{1/2}h$, where $P_2$ is the period of the outer planet. Note that these two estimations are very crude since we dropped too many order unity factors. We can improve the estimate by including a constant factor (which may depend on $\beta$ but no other parameter) such that the example in Fig. 14 lies close to the threshold of significant eccentricity overshoot. After this calibration, the condition of significant eccentricity overshoot becomes

$$T_{m,2} \lesssim 30P_2\mu_2^{-1}q^{-1/2}h. \qquad (21)$$

This predicts that the critical $q$ below which the planet is ejected is approximately $\propto T_{m,2}^{-2}$ (or $t_{w,0}^{-2}$) when $q$ is small, in agreement with Fig. 12. Equation (21) does not assume any specific migration model for the outer planet. If the outer planet undergoes type I migration, equation (21) can be rewritten [using equations (5) and

(6)] as
$$q \lesssim 3.6 \left(\frac{h}{0.05}\right)^{-2} \left(\frac{\Sigma}{1700 \text{ g} \cdot \text{cm}^{-2}}\right)^{2} \left(\frac{a_2}{1 \text{ au}}\right)^{4} \left(\frac{M_{\star}}{M_{\odot}}\right)^{-2}. \qquad (22)$$

This condition can be easily satisfied in typical protoplanetary discs, since the RHS is often $>1$ (while convergent migration already requires $q \lesssim 1$). For type-I migration, the condition of eccentricity overshoot depends on the mass ratio and disc density, but does not directly depend on the mass of the planets. Note that in addition to equation (22), eccentricity overshoot also requires the equilibrium eccentricity of the less massive planet to be $\gtrsim 3h$, so that eccentricity damping on it becomes negligible; this corresponds to $q \lesssim 0.15$, 0.04, and 0.03 for 1:2, 2:3, and 3:4 MMR, respectively, when $\beta = 0$.

It is worth noting that significant eccentricity overshoot is a phenomenon unique to the realistic migration model. For the simple migration model, since the eccentricity damping of the inner planet is efficient (i.e. $e_1^2/T_{e,1}$ always increases as $e_1$ increases), the system will cease to go deeper into the resonance once the $e_1^2$ terms in equation (19) can balance the migration; this corresponds to an insignificant eccentricity overshoot.

### 4.2.3 Effect of non-adiabatic evolution: bifurcation of the equilibrium state

In Fig. 13, we observe that the equilibrium eccentricity of the small planet increases abruptly when $q$ goes below $q \simeq 0.02$ (0.005) for $t_{w,0} = 10$ yr (20 yr); at a somewhat smaller $q$ the system becomes unstable. It is likely that this abrupt change corresponds to a bifurcation, with the equilibrium states before and after the bifurcation corresponding to two different fixed points of the system. One possible reason for this bifurcation is that the finite migration rate, together with the more realistic migration model, affects the stability of the fixed points. This different equilibrium state with a higher equilibrium eccentricity is not captured by our analytical result. For this new equilibrium state, we observe less eccentricity overshoot, and the scaling relation (21) for overshoot condition no longer holds.

As $t_{w,0}$ increases, the intermediate region where the system reaches this different equilibrium state with high eccentricity shrinks; when $t_{w,0}$ is sufficiently large, the system always becomes unstable (due to eccentricity overshoot) before the bifurcation happens and this intermediate region disappears.

## 4.3 Stability of capture

Similar to the case when the smaller planet is massless (Section 3), using the realistic migration model affects the stability of MMR capture. We observe that when the equilibrium eccentricity is $\gtrsim$ a few $h$, the system tends to be more stable compared to the prediction of the simple migration model. Since it is difficult to do a thorough survey of the parameter space, we illustrate this by an example. Figs 15 and 16 show the different outcomes of a 1:2 MMR capture when the simple migration model (eccentricity-independent $T_e$ and $T_m$) and the realistic migration model (eccentricity-dependent $T_e$ and $T_m$) are used. For the simple migration model, the equilibrium state is overstable, and the system eventually escapes the resonance. For the realistic migration model, the equilibrium state becomes stable (the eccentricity at the equilibrium also increases compared to the simple migration model).
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Outcome of the 1:2 MMR capture for the system depicted in Fig. 12 with $q = 0.01$ and $t_{w,0} = 100$ yr, using the simple migration model. We see that the system escapes the resonance at $t \approx 60$ kyr.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Same as Fig. 15, except that the realistic migration model is used. The equilibrium state is stable. Note that the equilibrium eccentricity is also increased compared to Fig. 15.</figcaption></figure>

Moreover, all numerical examples summarized in Figs 12 and 13 (except those cases where the inner planet is ejected) have stable equilibrium states. This suggests that for planets undergoing type I migration, the capture into a first-order MMR is stable for typical disk configurations if we use the realistic migration model. By contrast, if we use the simple migration model for the systems in Figs 12 and 13, the equilibrium state becomes unstable for $q \lesssim 0.1$.

Deck & Batygin (2015) have previously carried out an extensive study on the stability of the equilibrium state of first-order MMRs for general planet mass ratios. Their analysis was entirely based on the simple migration model. They found a region of the parameter space leading to overstability and proposed a criterion for overstability of the equilibrium state. Since an overstable system tends to evolve to an adjacent MMR equilibrium state, they concluded that

the overstability of the equilibrium state cannot fully explain the observed paucity of resonant pairs in the Kepler sample.

However, the overstability criterion of Deck & Batygin (2015) cannot be directly generalized to the realistic migration model (with eccentricity-dependent $T_m$, $T_e$). This is because the stability of the equilibrium state depends on both $\dot{a}/a$, $\dot{e}/e$ at the equilibrium and their partial derivatives with respect to the eccentricity. Although it is possible to tune the parameters of the simple migration model ($T_{e,0}$, $T_{m,0}$ for each planet and $p$) to obtain $\dot{a}/a$, $\dot{e}/e$ and $\partial(\dot{a}/a)/\partial e$ that locally match the values for the realistic migration model near the equilibrium, the local value of $\partial(\dot{e}/e)/\partial e$ in general cannot be matched by tuning the parameters of the simple migration model. Still, if we oversimplify the problem by plugging the local values of $T_e$, $T_m$ at the equilibrium into the overstability criterion of Deck & Batygin (2015), the stability does tend to increase compared to the simple migration model (with $T_e = T_{e,0}$, $T_m = T_{m,0}$) when $e_1 \gtrsim$ a few $h$. This is mainly because $T_e$ of the inner planet for the realistic migration model is larger than that for the simple migration model, which pushes the system away from the instability zone (see figs 2 and 3 of Deck & Batygin 2015). Ultimately, when $e_1$ is large enough so that $q^2 T_{e,1} > T_{e,2}$, the equilibrium state should become stable regardless of the total mass of the planets. For our migration model, this happens when $(e_1/h) \gtrsim q^{-2/3}$. This increased stability agrees with the trend shown in Fig. 3 around $e_{eq,0} \sim 2.4h$. Note that this is only an intuitive explanation of our finding of the increased stability and cannot serve as a rigorous analysis.

# 5 SUMMARY AND DISCUSSION

## 5.1 Summary of key results

In this paper, we have carried out theoretical and numerical studies on the outcomes of first-order MMR capture for planets undergoing convergent type I migration. Unlike previous works (Goldreich & Schlichting 2014; Deck & Batygin 2015; Delisle et al. 2015; Xu & Lai 2017), which adopted a simple migration model where the eccentricity damping rate and orbit decay rate [$T_e^{-1}$ and $T_m^{-1}$, respectively; see equations (1) and (2)] are independent of the planet's eccentricity, we consider a more realistic model for $T_e$ and $T_m$ which captures their non-linear eccentricity dependence when the eccentricity exceeds $\sim h$ (where $h \equiv H/r$ is the aspect ratio of the disc). We find that this more realistic migration model can significantly affect the outcomes of MMR capture and lead to several new dynamical behaviours.

First, the equilibrium eccentricities of planets captured into the MMR can be larger by a factor of a few than those predicted by the simple migration model (which assumes eccentricity-independent $T_e$, $T_m$). This arises because when $e \gtrsim 3h$, eccentricity damping becomes weaker and the system migrates deeper into the resonance before reaching equilibrium. When the inner planet is massless, the equilibrium state no longer exists if the equilibrium eccentricity predicted using the simple migration model is $\gtrsim 3h$, and the planet's eccentricity grows and eventually becomes unstable (Section 3.1). For general planet mass ratios (Section 4), the more massive planet's eccentricity stays below $3h$, and the eccentricity damping of this more massive planet ensures the existence of the equilibrium state. However, the equilibrium eccentricity is larger than the prediction using the simple migration model when the mass ratio $q = m_1/m_2$ is sufficiently small (Section 4.1). This increase in eccentricity is very significant for the 1:2 MMR, and less significant for other first-order MMRs (see Fig. 8). For typical disc parameters, the critical mass
ratio below which such increase occurs is around 0.03–0.15 (see Figs 8–9).

Secondly, the stability of the equilibrium state can be strongly affected by the migration model. Our analytical calculation and parameter survey for the case when the inner planet is massless (Section 3) show that the equilibrium state becomes more stable when the equilibrium eccentricity is $e_{eq} \gtrsim 2h$ (Section 3.2; see Fig. 3). This increased level of stability of the MMR is also seen when both planets have finite masses (Section 4.3). In particular, for realistic disc configurations, the simple migration model predicts that the equilibrium state is unstable for small *q*, while the realistic migration model predicts that the equilibrium state is virtually always stable (provided that the small planet does not suffer dynamical ejection at high eccentricities; see below).

Another new phenomenon we have found is that when the migration is fast and/or the inner planet's mass is sufficiently small, the eccentricity growth of the more massive planet (due to the resonant perturbation from the inner, smaller planet) becomes too slow; this causes the eccentricity of the smaller planet to overshoot the equilibrium value before the system reaches the equilibrium state (Section 4.2.2; see Fig. 13). Such an overshoot can be very significant and may cause the smaller planet to be ejected at high eccentricities even when the equilibrium eccentricity is modest. The condition of eccentricity overshoot is given in equation (21), and can be simplified to equation (22) if the outer planet undergoes type I migration – this condition can be satisfied for typical protoplanetary discs.

Overall, using the more realistic migration model tends to increase the equilibrium eccentricities of planets captured in MMRs and makes the equilibrium state less prone to overstability. However, when migration is sufficiently fast (or the small planet has too small a mass), it also causes the ejection of the smaller planet during eccentricity overshoot – this behaviour is much less significant when the simple migration model is used. All of these can affect the ways in which MMRs shape planetary system architecture.

## 5.2 Implications for multiplanet system architecture

### 5.2.1 Occurrence of MMRs and mass ratio in multiplanet systems

For planets with similar masses ($q \sim 1$), since the equilibrium eccentricities of the planets captured into MMRs are usually small, previous results concerning the stability of MMRs remain valid (Deck & Batygin 2015; Delisle et al. 2015; Xu & Lai 2017). observationally, short-period, compact systems (e.g. Kepler multiplanet systems) are found to have small eccentricities (Van Eylen et al. 2018) and to exhibit intra-system uniformity in mass and radius (Millholland, Wang & Laughlin 2017; Weiss et al. 2018). Therefore, the occurrence of MMR in most observed systems can be explained without invoking a more realistic migration models.

For a smaller mass ratio ($q \lesssim 0.1$), however, the maximum eccentricity that the smaller planet can reach is much larger when the realistic migration model (with eccentricity-dependent $T_e, T_m$) is applied (compared to the results obtained with eccentricity-independent $T_e, T_m$) due to the increased equilibrium eccentricity and eccentricity overshoot. The large eccentricity can lead to the ejection of the smaller planet (when its eccentricity approaches unity) or make it scatter with a third planet in the system (if their orbits cross). This tends to reduce the multiplicity of the system when it initially hosts a pair of convergently migrating planets with a small mass ratio. Therefore, it is possible that the unexpected intra-system mass and radius uniformity in Kepler multiplanet sys-

tems (Millholland et al. 2017) is partially caused by the removal of small planets during MMR encounters.

### 5.2.2 Loneliness of hot Jupiters

The eccentricity overshoot phenomenon (which occurs when the mass ratio is small and the migration is sufficiently fast) provides an efficient way of removing super-Earth companions of fast-migrating giant planets. This may help explain the loneliness (the lack of low-mass planet neighbours) of hot Jupiters (Huang, Wu & Triaud 2016) if they are formed through disc-driven migration.<fnref n="7" /> In this picture, hot Jupiters arrived at their current locations through fast type II migration, with the migration time-scale much less than the disc lifetime.<fnref n="8" /> If they had any inner low-mass companion (a super-Earth), it could be removed when captured into an MMR with the Jupiter during its migration due to the instability caused by eccentricity overshoot. On the other hand, warm Jupiters do commonly have low-mass-planet companions (Huang et al. 2016). This may be explained by their slow migration rates: During such slow migration, their low-mass companions do not suffer eccentricity overshoot and therefore are kept in safety upon capture into MMRs. Note that the rate of type II migration is sensitive to the property of the disc, especially its viscosity (Ward 1997). Thus, in this scenario, whether a system forms hot Jupiters (without low-mass companions) or warm Jupiters (with low-mass companions) simply reflects the different disc properties and the resulting different migration history of giant planets.

Of course, hot Jupiters may also form by high-eccentricity migration, in which the eccentricity of a giant planet is excited by distant stellar or planetary companions, followed by tidal circularization and orbital decay (e.g. Dawson & Johnson 2018). In this scenario, the loneliness of hot Jupiters can be naturally explained because a giant planet undergoing high-amplitude eccentricity oscillations can easily eject smaller planets interior of its initial orbit.<fnref n="9" /> Our discussion here does not aim to prove or disprove any particular formation scenario; we simply argue that one should not rule out the disc-driven (low-eccentricity) migration scenario using solely the loneliness of hot Jupiters.

## ACKNOWLEDGEMENTS

We thank M. Lambrechts, E. Lega, and G. Pichierri for useful discussions and comments. This work has been supported in part by NASA grant NNX14AG94G and NSF grant AST-1715246. WX thanks the undergraduate research fellowship from the Hopkins Foundation (Summer 2017). DL thanks the Laboratoire Lagrange OCA for hospitality where this work started.

<footnote n="7">⁷Giant planets should undergo type II instead of type I migration. Although in this paper we have focused on type I migration models, our results should also be reasonably accurate if the more massive planet undergoes type II migration, since the eccentricity dependence of the more massive planet's migration and eccentricity damping rates does not play an important role in our analysis.</footnote>

<footnote n="8">⁸Type II migration can be fast enough to push a Jupiter to the disc inner edge before the gas disperses (Hasegawa & Ida 2013). Some of these hot Jupiters may not fall into the star probably because their migration stops when reaching the inner edge of the disc or when the inner part of the disc induces outward migration (e.g. Lega et al. 2015).</footnote>

<footnote n="9">⁹However, in this scenario, warm Jupiters with an inner companion should be formed through a different channel.</footnote>
REFERENCES

André Q., Papaloizou J. C. B., 2016, MNRAS, 461, 4406

Baruteau C. et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. University of Arizona Press, Tucson, p. 667

Batalha N. M. et al., 2013, ApJS, 204, 24

Batygin K., Adams F. C., 2017, AJ, 153, 120

Batygin K., Morbidelli A., 2013, AJ, 145, 1

Chambers J. E., 1999, MNRAS, 304, 793

Chatterjee S., Ford E. B., 2015, ApJ, 803, 33

Chiang E. I., Goldreich P., 1997, ApJ, 490, 368

Coughlin J. L. et al., 2016, ApJS, 224, 12

Cresswell P., Nelson R. P., 2008, A&A, 482, 677

Cresswell P., Dirksen G., Kley W., Nelson R. P., 2007, A&A, 473, 329

Crida A., Sándor Z., Kley W., 2008, A&A, 483, 325

Dawson R. I., Johnson J. A., 2018, preprint (arXiv:1801.06117)

Deck K. M., Batygin K., 2015, ApJ, 810, 119

Delisle J.-B., Laskar J., Correia A. C. M., 2014, A&A, 566, A137

Delisle J.-B., Correia A. C. M., Laskar J., 2015, A&A, 579, A128

Fabrycky D. C. et al., 2014, ApJ, 790, 146

Goldreich P., Schlichting H. E., 2014, AJ, 147, 32

Goldreich P., Tremaine S., 1980, ApJ, 241, 425

Hasegawa Y., Ida S., 2013, ApJ, 774, 146

Huang C., Wu Y., Triaud A. H., 2016, ApJ, 825, 98

Izidoro A., Ogihara M., Raymond S. N., Morbidelli A., Pierens A., Bitsch B., Cossou C., Hersant F., 2017, MNRAS, 470, 1750

Kley W., Lee M. H., Murray N., Peale S. J., 2005, A&A, 437, 727

Lega E., Morbidelli A., Bitsch B., Crida A., Szulágyi J., 2015, MNRAS, 452, 1717

Lithwick Y., Wu Y., 2012, ApJ, 756, L11

Luger R. et al., 2017, Nature Astronomy, 1, 0129

Michtchenko T. A., Ferraz-Mello S., 2001, Icarus, 149, 357

Migaszewski C., 2015, MNRAS, 453, 1632

Millholland S., Wang S., Laughlin G., 2017, ApJ, 849, L33

Mills S. M., Fabrycky D. C., Migaszewski C., Ford E. B., Petigura E., Isaacson H., 2016, Nature, 533, 509

Murray C. D., Dermott S. F., 1999, Solar System Dynamics. Cambridge Univ. Press

Nelson R. P., 2018, Handbook of Exoplanets, in Deeg H. J. Belmonte J. A., eds, Belmonte, Springer

Papaloizou J. C. B., Szuszkiewicz E., 2005, MNRAS, 363, 153

Pu B., Wu Y., 2015, ApJ, 807, 44

Tanaka H., Ward W. R., 2004, AJ, 602, 388

Terquem C., Papaloizou J. C. B., 2007, ApJ, 654, 1110

Van Eylen V. et al., 2018, preprint (arXiv:1807.00549)

Ward W. R., 1997, Icarus, 126, 261

Weiss L. M. et al., 2018, AJ, 155, 48

Xu W., Lai D., 2017, MNRAS, 468, 3223

Zhang X., Li H., Li S., Lin D. N. C., 2014, ApJ, 789, L23

This paper has been typeset from a *TEX/LATEX* file prepared by the author.