DISK–PLANETS INTERACTIONS AND THE DIVERSITY OF PERIOD RATIOS IN KEPLER'S MULTI-PLANETARY SYSTEMS

Figure 2 displays typical outcomes of our hydrodynamical simulations of the Kepler-46 system. Here, λ_{i, o} and ω_{i, o} denote the mean longitude and longitude of periastron, with the subscripts i and o referring to the inner and outer planets, respectively.

Zoom In Zoom Out Reset image size

The top row shows results with M_inner = 0.4 M_J and Σ₀ = 3 × 10⁻⁴. As the planets progressively clear a partial gap around their orbit, convergent migration does not occur at a constant rate. In particular, the inward migration of the outer planet slightly accelerates early due to a stage of type-III migration (Masset & Papaloizou 2003) before slowing down upon approaching the 2:1 MMR with the inner planet. The planets are locked in the 2:1 MMR from about 300 orbits onward, as evidenced by the libration of the resonant angles associated with this resonance. The resonance capture increases the inner planet's eccentricity to about 0.035. The increase in the outer planet's eccentricity is much more modest, reaching about 0.005, even though the mass ratio of the planets is near unity. This modest increase in the outer planet's eccentricity is a result of the contribution of the indirect term nearly canceling the contribution from the direct term to the term that is first order in the eccentricities in the development of the disturbing function, which governs the response of the outer planet at a 2:1 MMR (see, e.g., Murray-Clay & Chiang 2005; Papaloizou & Terquem 2010). Damping of the eccentricities follows from the gravitational interaction with the background gas disk. The ratio of orbital periods increases from 2.02 to 2.15 in only 1,500 orbits, while the eccentricities are damped. The physical origin of this divergent evolution is discussed in Section 2.2.2. Resonant coupling is maintained throughout this divergent evolution, the width of the 2:1 MMR increasing as the eccentricities decrease.

The second row of panels in Figure 2 displays results for M_inner = 0.6 M_J and Σ₀ = 6 × 10⁻⁴. This case shows that convergent migration may stall with an orbital period ratio close to the observed value (≈1.7). The accelerating inward migration of the outer planet is particularly visible in the first 100 orbits. This runaway migration is fueled by the dense gas region located between the two planet gaps. This dense material is continuously funneled beyond the outer planet's orbit through horseshoe streamlines. The density between the two planet gaps therefore decreases, and so does the migration rate of the outer planet. Convergent migration then proceeds at a slower pace and the planets merge their gap, thereby forming a common gap in the disk. This is illustrated in Figure 3, which compares the disk's surface density at 400 orbits for the two aforementioned simulations. While in the former case a narrow ring of gas is left between the two planet gaps, the latter case shows that both planets end up in a common gap.

Zoom In Zoom Out Reset image size

In the second row of panels in Figure 2, we see that the orbital period ratio and the eccentricities remain approximately stationary from about 400 orbits. The time-averaged eccentricity of the inner planet is close to 0.01; that of the outer planet is ≈0.025. These values are consistent, albeit slightly larger than those inferred from observations (see Table 1). The critical angles show that the planets are not locked in the 5:3 MMR.

The third row of panels in Figure 2 shows M_inner = 0.2 M_J and Σ₀ = 8 × 10⁻⁴. Here, convergent migration is rapid enough that the planets approximately reach the nominal location of their 3:2 MMR (the minimum value of the period ratio is 1.508 at 170 orbits). The planet eccentricities increase to about 0.03 (inner planet) and 0.015 (outer planet). Subsequent damping of the eccentricities occurs along with a slight divergent evolution of the planets. The rate of divergent evolution is much reduced compared with the 2:1 MMR case shown in the top row of panels in Figure 2. We analyze this result further in Section 2.2.2. No steady state has been reached at the end of the simulation. For this simulation, as for all simulations of the Kepler-46 system, results are shown until the inner planet reaches r ∼ 0.4, below which the proximity of the grid's inner edge (at r = 0.2) starts affecting the planet's orbital evolution.

The final period ratio obtained in all our simulations is displayed in Figure 4 as a function of Σ₀ for the three masses assumed for the inner planet. By final period ratio, we mean either the period ratio obtained at the end of the simulation (after between 1500 and 2000 planet orbits), or that obtained when the inner planet's orbital radius falls below r = 0.4. Figure 4 highlights the three major outcomes of our simulations: (1) capture in 2:1 MMR generally followed by rapid divergent evolution, (2) capture in 3:2 MMR followed by slow divergent evolution, and (3) convergent migration stalling with a period ratio between 1.6 and 1.7, depending on the disk's density distribution inside the planets' common gap. Several simulations can reproduce the observed period ratio of Kepler-46b and Kepler-46c. The best agreement with observations is obtained for M_inner = 0.6 M_J, our upper mass value. This mass would imply a mean density for Kepler-46b about 20% larger than Jupiter's (about 1.6 g cm⁻³).

Zoom In Zoom Out Reset image size

The results of hydrodynamical simulations in Section 2.2.1 show that the orbital period ratio of two partial gap-opening planets may increase from a near resonant value. This divergent evolution is reminiscent of the late evolution of short-period planet pairs through star–planet tidal interactions. In late evolution cases, tidal circularization causes a slow divergent evolution of the orbits (Papaloizou & Terquem 2010; Papaloizou 2011; Lithwick & Wu 2012; Batygin & Morbidelli 2013), a mechanism known as tidally driven resonant repulsion (Lithwick & Wu 2012). Tidally driven resonant repulsion is based on energy dissipation of the planet pair while its total angular momentum remains constant (Papaloizou 2011).

It is possible that there is an analogous resonant repulsion mechanism due to disk–planets interactions. To address this possibility, it is important to first look at the evolution of two embedded planets that do not interact gravitationally. In Figure 5, we compare the results of two hydrodynamical simulations: one in which the gravitational interaction between the planets is switched off (both the direct and the indirect interaction terms are discarded), and another for which the planet–planet interaction is fully accounted. The inner planet's mass is M_inner = 0.6M_J and the initial disk surface density at r = 1 is Σ₀ = 2 × 10⁻⁴. The other disk and planet parameters are those described in Section 2.1. When planet–planet interactions are on, capture into the 2:1 MMR is followed by rapid divergent evolution and damping of the eccentricities. This is similar to the case presented in the top row of panels in Figure 2. We see that a similar divergent evolution is also obtained when planet–planet interactions are switched off, while the planets are not resonantly coupled (the resonant angles do not librate) and their eccentricities take vanishingly small stationary values. This comparison highlights that when disk–planet interactions are considered, eccentricity damping is not necessarily responsible for the planets' divergent evolution. The divergent evolution of a planet pair may instead be driven by changes to the disk–planet interaction brought about through the proximity of the planets, such as the development of an interaction between each planet and the wake of its companion. We argue below that the damping of the eccentricities can actually be a consequence of divergent evolution driven by wake-planet interactions.

Zoom In Zoom Out Reset image size

Before describing the process of wake-planet interactions, we briefly report that in some other simulations not illustrated here, with planet–planet interactions switched off, a tiny increase in the orbital eccentricities is observed when the orbital period ratio takes near resonant values. Eccentricities may reach a few × 10⁻³ before being damped by the disk. Although this is a small effect, it suggests that planets could be weakly resonantly coupled through interaction with the wakes and/or the circumplanetary disk of their companions. In our simulations, the mass of the circumplanetary disks does not exceed one percent of the planet's mass.

Wake-planet interactions. The results of simulations shown in Figure 5 make clear that disk-driven divergent evolution of a planet pair is not necessarily driven by eccentricity damping and it may continue to operate even when the eccentricities are very small. We indicate how the interaction between a planet and the wake of its companion could act as a driving mechanism for divergent evolution. The inner planet's outer wake carries positive fluxes of angular momentum and energy. A fraction of this angular momentum flux can be deposited in the co-orbital region of the outer planet through the dissipation of shocks, and then transferred to the outer planet via an effective positive corotation torque (see Podlewska-Gaca et al. 2012). Associated with the angular momentum flux, F_J, deposited in the co-orbital region of the outer planet is an energy flux, F_E = n_iF_J, where n_i denotes the inner planet's mean motion (or angular velocity). From this energy flux, n_oF_J is used to keep the outer planet and its co-orbital region nearly circular (see Figure 5; here, n_o is the outer planet's mean motion). The remaining energy flux, (n_i − n_o)F_J, is dissipated in the outer planet's co-orbital region. Reciprocally, a fraction of the negative inward flux of angular momentum carried by the inner wake of the outer planet can be deposited in the inner planet's co-orbital region and subsequently transmitted to the inner planet. Again, this exchange is associated with energy dissipation in the disk at a rate equal to the absolute value of the product of the flux of angular momentum transmitted to the inner planet by the outer one and the difference in the planets angular velocities. This dissipation occurs, as is usual in disk–planet interactions, through the tidally excited density waves steepening into shocks (e.g., Goodman & Rafikov 2001). Numerically, shocks are handled through a combination of numerical diffusion and the application of an artificial viscosity. As a locally isothermal equation of state is adopted in our simulations, this dissipated energy is effectively lost to the system through cooling. Though the total orbital energy of the two planets is dissipated due to the wake-planet interactions described above, their total angular momentum remains unaltered. As studied in the context of star–planet tidal interactions, energy dissipation at constant angular momentum leads to the divergent evolution of a planet pair (Papaloizou 2011; Lithwick & Wu 2012; Batygin & Morbidelli 2013).

We derive an equation satisfied by the period ratio of the two planets as a direct consequence of considering the evolution of their total orbital energy and angular momentum. We denote by a_i, e_i, and M_i, the semi-major axis, eccentricity, and mass of the inner planet (with the subscript changed to o denoting the same quantities for the outer planet), respectively. The total angular momentum of the two planets is

$$\begin{equation} J = J_{\rm i} + J_{\rm o} = M_{\rm i} \sqrt{GM_{\star }a_{\rm i} \left(1-e_{\rm i}^2\right)} + M_{\rm o} \sqrt{GM_{\star }a_{\rm o} \left(1-e_{\rm o}^2\right)}, \end{equation} \tag{ 1 }$$

which gives, correct to second order in the eccentricities,

$$\begin{equation} \dot{J} \equiv \frac{dJ}{dt} = J_{\rm i} \left(\frac{\dot{a}_{\rm i}}{2a_{\rm i}} - e_{\rm i} \dot{e}_{\rm i}\right) + J_{\rm o} \left(\frac{\dot{a}_{\rm o}}{2a_{\rm o}} - e_{\rm o} \dot{e}_{\rm o}\right). \end{equation} \tag{ 2 }$$

Since we are looking for changes in the planets' semi-major axes and eccentricities that occur on similar timescales, we may therefore discard the eccentricity terms in the expressions for $\dot{J},$ J_i, and J_o. Denoting the torques exerted by the disk on the inner and outer planets by Γ_i and Γ_o, respectively, conservation of angular momentum reads

$$\begin{equation} \dot{J} = J_{\rm i} \frac{\dot{a}_{\rm i}}{2a_{\rm i}} + J_{\rm o} \frac{\dot{a}_{\rm o}}{2a_{\rm o}} = \Gamma _{\rm i} + \Gamma _{\rm o}. \end{equation} \tag{ 3 }$$

The total orbital energy is given by

$$\begin{equation*} E = -\frac{GM_{\star }M_{\rm i}}{2a_{\rm i}} -\frac{GM_{\star }M_{\rm o}}{2a_{\rm o}}, \end{equation*}$$

and conservation of energy reads

$$\begin{eqnarray} \dot{E} & = & \frac{GM_{\star }M_{\rm i}}{2a_{\rm i}^2} \dot{a}_{\rm i} + \frac{GM_{\star }M_{\rm o}}{2a_{\rm o}^2} \dot{a}_{\rm o}\nonumber \\ & = & n_{\rm i} \Gamma _{\rm i} + n_{\rm o} \Gamma _{\rm o} - \frac{GM_{\star }M_{\rm i} e_{\rm i}^2}{a_{\rm i} \tau _{\rm c,i}} - \frac{GM_{\star }M_{\rm o} e_{\rm o}^2}{a_{\rm o} \tau _{\rm c,o}} - |\dot{E}|_{\rm wake},\nonumber\\ && \end{eqnarray} \tag{ 4 }$$

where $n_{\rm i,o} = \sqrt{\vphantom{A^A}\smash{{{GM_{\star } / a^3_{\rm i,o}}}}}$ is the planet's mean motion. The right-hand side of Equation (4) contains three contributions. The first arises from the energy dissipation in the disk resulting from the action of the torques Γ_{i, o}. Since these are negative, this energy is removed from the orbital motion. The second originates from orbital circularization due to disk–planets interactions. The quantities τ_{c, i} and τ_{c, o} are the circularization times for the planets that arise from the disk torque and energy dissipation acting on each planet through the effects of its own wake (Nelson & Papaloizou 2002; Papaloizou & Szuszkiewicz 2005). The third and final contribution, $-|\dot{E}|_{\rm wake}$, stems from the energy dissipation due to wake-planet interactions. We remark that wake dissipation in the co-orbital region of the outer planet contributes (n_i − n_o)|F_J| to $|{\dot{E}}|_{\rm wake}$, in addition to a corresponding contribution arising from the co-orbital region of the inner planet, as indicated above. Terms of the order of $e_{\rm i,o}^2$ are retained above even though the eccentricities are small because it is assumed that τ_{i, o} are small enough so that dissipation through circularization could be as important as, or even dominate, the other effects under some circumstances.

Combining Equations (3) and (4), we find

$$\begin{eqnarray} \frac{1}{3}\frac{n_{\rm o}}{n_{\rm i}}\frac{d}{dt}\left(\frac{n_{\rm i}}{n_{\rm o}} \right) &=& \left(\frac{\Gamma _{\rm o}}{J_{\rm o}} - \frac{\Gamma _{\rm i}}{J_{\rm i}} \right) + \frac{J_{\rm i} + J_{\rm o}}{J_{\rm i} J_{\rm o} (n_{\rm i} - n_{\rm o})} \nonumber \\ && \times\left(n_{\rm i} J_{\rm i} \frac{e_{\rm i}^2}{\tau _{\rm c,i}} + n_{\rm o} J_{\rm o} \frac{e_{\rm o}^2}{\tau _{\rm c,o}} + |\dot{E}|_{\rm wake} \right). \end{eqnarray} \tag{ 5 }$$

When wake-planet interactions are negligible ($|\dot{E}|_{\rm wake}=0$), Equation (5) shows that far from resonance, where the eccentricities are very small, the planets' period ratio n_i/n_o decreases as a consequence of convergent migration, provided that Γ_o/J_o < Γ_i/J_i. As expected, the latter condition is that the inward migration rate of the outer planet should exceed the inward migration rate of the inner planet. As the eccentricities increase as resonance is approached, a steady state may be reached in which the planets' period ratio and eccentricities take stationary values. If the eccentricity damping timescales are sufficiently short, disk-driven circularization of the orbits can reverse convergent migration, thus the planets' period ratio can increase.

Equation (5) highlights that wake-planet interactions may also lead to divergent evolution of a planet pair. Divergent evolution driven by wake-planet interactions does not require the planets to be resonantly coupled, as already inferred from the simulations shown in Figure 5. Similar divergent evolution of a non-resonant planet pair embedded in a disk, comprised of a Jupiter and a super-Earth, was reported by Podlewska-Gaca et al. (2012). When divergent evolution is primarily due to wake-planets interactions, as the period ratio of the planet pair increases away from a resonant value, the distance to resonance increases, and the forced eccentricities therefore decrease (e.g., Papaloizou 2011). This means that orbital circularization is a natural consequence of divergent evolution driven by wake-planet interactions. This is a major difference between the divergent evolutions driven by wake-planet interactions and by orbital circularization (through the disk or stellar tides): eccentricity damping is the consequence of the former mechanism and the cause of the latter. With both mechanisms, the divergent evolution of a resonant planet pair maintains the libration of the resonant angles by damping the eccentricities.

We call disk-driven resonant repulsion the divergent evolution of a resonant planet pair embedded in a disk, following the terminology introduced by Lithwick & Wu (2012) for the divergent evolution of a resonant planet pair due to star–planet tidal interactions. We stress that disk-driven resonant repulsion can have two origins: (1) disk circularization of the orbits and (2) wake-planet interactions. Disk-driven repulsion comprises (1) circularization-driven repulsion and (2) wake-driven repulsion. The results of hydrodynamical simulations shown above indicate that for the disk and planet parameters that we have considered, wake-planet repulsion dominates circularization-driven repulsion. This point will be further illustrated in Section 3.2.1.

When is wake-driven repulsion important? We have shown that the period ratio of a planet pair embedded in a disk may increase due to angular momentum transfer between the planets through their wakes. We separate the planet pair into a donor and a recipient. The bigger the donor's mass, the stronger its wake, and the easier for the wake to deposit angular momentum in the co-orbital region of the recipient planet through the dissipation of shocks. This typically requires the dimensionless parameter q/h³ of the donor planet to exceed unity (q is the planet-to-star mass ratio and h the disk's aspect ratio). The donor planet can thus be a partial or a deep gap-opening planet. The angular momentum that is deposited by the donor's wake is at least partly transferred to the recipient planet through an effective corotation torque (Podlewska-Gaca et al. 2012). If the recipient's mass is too big, the density in its co-orbital region will be too low for this effective corotation torque to matter. The recipient could be a partial gap-opening planet or even a type-I migrating planet. These considerations show why divergent evolution is not expected for two type-I migrating planets (the donor's wake being too weak), or for two type-II migrating planets (the recipient's gap being too deep). We have checked this with hydrodynamical simulations, not illustrated here.

The above considerations also help in understanding the different rates of divergent evolution obtained in Figure 2. The opening of a partial gap for both planets is a favorable situation for wake-planet interactions to play a significant role. In the top row of panels in Figure 2, where the planets become locked in the 2:1 MMR, the high-density region between the two gaps (see upper panel in Figure 3) causes significant deposition of angular momentum in the planets co-orbital region and, therefore, efficient divergent evolution. In contrast, as seen in the bottom row of panels in Figure 2, where the planets merge their gap and become locked in the 3:2 MMR, the moderate disk density left between the planets can only trigger slow divergent evolution. In the middle row of panels, where the period ratio stalls near 1.7, a steady state is reached with background convergent migration balanced by wake-planet divergent evolution.

Evolution toward a steady state. When the planets are widely separated, they are expected to interact with the disk independently of each other and they can then enter a convergent migration phase. As they approach each other, both planet–planet interactions and wake-planet interactions, where the planets interact with each others wakes, can occur. Both these processes involve exchanges of orbital energy and angular momentum between the planets. Planet–planet interactions can result in resonant trapping associated with the growth of eccentricities, while wake-planet interactions act to produce divergent migration without the growth of eccentricities. From Equation (5) we see that both circularization- and wake-driven resonant repulsions act together to increase the period ratio. If wake-driven repulsion dominates after the eccentricities have damped to low values, with $|{\dot{E}}|_{{\rm wake}}$ remaining constant, the period ratio should increase linearly with time; this result is roughly what is obtained in the simulation presented in the top row of panels in Figure 2. However, in the longer term, $|{\dot{E}}|_{{\rm wake}}$ should decrease with time as the planets separate, which indicates that a steady state should ultimately be reached with convergent migration driven by background disk torques balancing wake-driven repulsion. Wake-planet interactions cause a non-linear evolution of the disk's density profile near the planets, which may affect the time evolution of the background torques responsible for convergent migration. At least in the simulation without the planet–planet interaction of Figure 5, a slow decrease in the background torques could explain why the period ratio decreases and then increases, instead of reaching a stationary non-resonant value. Another possibility is that convergent migration and disk-driven repulsion behave differently as the planet pair moves closer to the star, if the torques responsible for both mechanisms have different radial dependences. These aspects require further study.

The above simulations of the Kepler-46 system show that the planets' period ratio can take a range of values, depending on the convergent migration rate during the early evolution of the system. The limited duration of the simulations (up to a few thousand orbits) and the fact that planets still migrate at the end of the simulations, however, raise the question of the long-term evolution of the system. Even when the inward migration of the planets is stalled,² further evolution of the period ratio can be maintained by disk-driven resonant repulsion. A steady state may a priori be achieved after depletion of the protoplanetary disk. Disk dispersal typically occurs between 10⁶ and 10⁷ yr after formation of the central star. Photoevaporation, driven by extreme ultraviolet radiation from the star, opens a gap in the disk at separations of typically a few astronomical units (AU) for Sun-like stars (e.g., Owen et al. 2010; Alexander & Pascucci 2012). For planets below that separation, the background disk will be depleted over a viscous timescale at the orbital separation where photoevaporation occurs. Assuming a local viscous alpha parameter of a few × 10⁻³, viscous draining of the inner disk could operate in a few 10⁴ yr, that is typically a few 10⁵ orbital periods at the present location of Kepler-46b and Kepler-46c.

To assess the impact of disk evaporation on our results, we restart two simulations adopting a simple exponential decay of the surface density profile. This is done by solving $\partial _t \Sigma = -(\overline{\Sigma } - \Sigma _{\rm target})/\tau _{\rm evap}$, in addition to the hydrodynamical equations, where $\overline{\Sigma }$ is the azimuthally averaged density profile at restart time (from which evaporation switches on), $\Sigma _{\rm target} = 10^{-3}\overline{\Sigma }$ is an arbitrarily small density profile in a steady state (taken to be not zero for numerical convenience), and τ_evap is the evaporation timescale. We consider three short evaporation timescales: 100, 500, and 2000 orbits, simulations being restarted at 1000 orbits. Results are shown in Figure 6.

Zoom In Zoom Out Reset image size

The upper panels are obtained with M_inner = 0.4 M_J and Σ₀ = 3 × 10⁻⁴, for which a rapid divergent evolution occurs after capture in the 2:1 MMR. We see that the evolution is frozen out at a few evaporation timescales after restart, the planets reaching stationary eccentricities and semi-major axes. Similarly, the lower panels in Figure 6 show the results for M_inner = 0.6 M_J and Σ₀ = 6 × 10⁻⁴, for which convergent migration stalls at about the observed period ratio in the Kepler-46 system. In that case, including disk evaporation stalls the migration of both planets, as expected, but it does not significantly affect the planets' eccentricities and period ratio. The final period ratio remains close to the observed value.

As highlighted in Figure 1, Kepler's multi-planetary systems show a clear tendency for planet pairs near resonances to feature period ratios slightly greater than strict commensurability. Divergent evolution of a short-period resonant planet pair due to star–planet tidal interactions could partly explain this feature (Lithwick & Wu 2012; Batygin & Morbidelli 2013). Although the efficiency of tidal dissipation remains largely uncertain, it is thought that efficient tidal circularization requires the inner planet to orbit its star in less than a few days. As can be seen in Figure 1, many of Kepler's candidate systems have an inner planet beyond 10 days, for which tidal circularization is unlikely to have caused significant resonant repulsion (see Appendix A). For these systems we propose that resonant repulsion driven by disk–planets interactions could have played a prominent role in determining the observed period ratios.

The vast majority of the planets in Kepler's multi-planetary systems have physical radii between 0.1 and 0.3 Jupiter radii, which places them in the super-Earth to Neptune-mass range. From the mass–radius diagram of Kepler candidates followed up by radial velocity, the median mass of a 0.2 Jupiter-radius planet is about 10 to 15 Earth-masses.³ Based on generally considered temperatures and viscosities in protoplanetary disks (corresponding to h = 0.05 and α ∼ a few × 10⁻³), such planets are expected to experience type-I migration. Recent N-body experiments by Rein (2012) have shown that, in contrast with Kepler's data, convergent migration of type-I migrating planets yields very little departure from strict commensurability. He showed that the inclusion of stochastic forces, in addition to type-I migration, could nicely reproduce the observed distribution of period ratios. Such agreement has been obtained assuming that the initial period ratio distribution is the observed one, and taking the same convergent migration timescale and the same amplitude of stochastic forces for all of Kepler's multiple systems, which remains uncertain given the expected diversity of physical properties in protoplanetary disks. Still, Rein's (2012) results indicate that type-I migration alone may not be able to account for Kepler's diversity of period ratios.

Based on our results of hydrodynamical simulations for the Kepler-46 system, we propose that some of Kepler's planet candidates in multi-planetary systems could have opened partial gaps in their parent disk, thereby escaping the type-I migration regime. This could be the case if, for instance, planets formed and/or migrated in regions of low turbulent activity (dead zones; see, e.g., Fleming & Stone 2003) or if the disk's aspect ratio takes smaller values than commonly adopted. Aspect ratios ≲ 3% are likely typical of the short orbital separations at which Kepler's candidate systems are detected. To address this possibility, we carry out hydrodynamical simulations with two planets in the super-Earth mass range. The inner planet's mass, M_inner = 15 M_⊕, the outer planet's mass, M_outer = 13 M_⊕, and the central star is one Solar mass. The planet-to-primary mass ratios are therefore q_inner = 4.4 × 10⁻⁵ and q_outer = 4 × 10⁻⁵, that is, an order-of-magnitude smaller than those of Kepler-46b (median value) and Kepler-46c. To get similar gap depths in the super-Earths and Kepler-46 simulations, we adopt values of h and α such that the dimensionless parameters in the gap-opening criterion of Crida et al. (2006) take the same values. For the super-Earth's simulations, this yields h = 0.023 and α = 2.3 × 10⁻³. All other simulation parameters are otherwise identical, except the grid resolution, which we increase to 600 × 1200.

Figure 7 displays the results of two simulations. In the top panels, the unperturbed disk's surface density is Σ₀ = 8 × 10⁻⁵. We see that the planets lock themselves in the 3:2 MMR from about 1300 orbits. Divergent evolution proceeds with a decrease in the planets' eccentricities. After reaching exactly 1.5 at about 1500 orbits, the period ratio increases to 1.54 at 5000 orbits. As illustrated by the screenshot of the disk surface density, the proximity of the planets allows their wakes to penetrate into each others coorbital region, where they can deposit part of their angular momentum flux. As we show in Section 2.2.2, transfer of angular momentum between planets through wake-planet interactions will result in a loss rate of the total orbital energy of the two planets. This has the consequence of increasing the period ratio and decreasing the eccentricities. In the bottom panels of Figure 7, the unperturbed disk surface density parameter is decreased to Σ₀ = 3 × 10⁻⁵. Here, convergent migration due to the background disk torques is slow enough to lock the planets in the 2:1 MMR, but, this time, the resonant capture does not lead to divergent evolution. Instead, a steady state is reached in which both the eccentricities and the orbital period ratio take stationary values. This steady state is a consequence of negligible wake-planet interactions in this case. As illustrated in the bottom-right panel of Figure 7, the planets are sufficiently far from each other that their wakes primarily transfer energy and angular momentum to regions of the disk that are not co-orbital with the planets.

Zoom In Zoom Out Reset image size

Our three-body simulations solve the equations of motion for the two planets and the central star using a standard technique. Planets are assumed to be coplanar. Disk–planets interactions are incorporated by applying appropriate dissipative forces (for details, see Papaloizou 2011). Equation (5) shows that in order to model the planet's evolution, we need to specify (1) the rate of convergent migration before the onset of disk-driven repulsion, (2) the damping rate of the eccentricities, (3) the efficiency of wake-driven divergent evolution, and (4) the method of modeling of disk dispersal. The efficiency of wake-driven divergent evolution remains uncertain at this stage, since the term $|\dot{E}|_{\rm wake}$ in the right-hand side of Equation (5), which models energy dissipation due to wake-planet interactions, should depend sensitively on the planet masses, their mutual separation and distance from the central star, and several disk quantities (including its temperature and viscosity). A systematic study of this dependence goes beyond the scope of this paper, and is therefore left for future work. To do simple modeling and to minimize the number of free parameters, we do not attempt to parameterize the $|\dot{E}|_{\rm wake}$ term. Instead, we model disk-driven repulsion as circularization-driven repulsion only. In Equation (5), we compensate the absence of the $|\dot{E}|_{\rm wake}$ term in our three-body integrations by considering circularization timescales τ_{c, i} and τ_{c, o}, which are shorter than in the hydrodynamical simulations. Damping of the planets' eccentricity is modeled as an exponential decay with a constant damping timescale to orbital period ratio, which we take as a free parameter (specified below).

To model the convergent migration prior to resonant repulsion, we adopt a simple prescription that reproduces the results of hydrodynamical simulations. Denoting by R, the ratio of semi-major axes (R = a_outer/a_inner), we found good agreement using

$$\begin{equation} R(t) = R(t=0) \times 10^{f} \end{equation} \tag{ 6 }$$

with

$$\begin{equation} f = \log _{10}\left[ \frac{{\rm (pr)}^{2/3}}{R(t=0)} \right] \times \frac{(t / \tau)^3 }{\sqrt{1+(t/\tau)^6}}, \end{equation} \tag{ 7 }$$

where pr and τ are fitting parameters. The expression in Equation (6) tends to a constant value for t ≫ τ, which mimics the fact that convergent migration is found to nearly stall before the onset of resonant repulsion. In Equation (7), pr denotes the period ratio at which convergent migration approximately stalls, which occurs at t = τ. The good agreement between this migration prescription and our hydrodynamical simulations is illustrated in Figure 8. In the three-body simulations, disk-driven migration is applied to the outer planet only.

Zoom In Zoom Out Reset image size

We display, in Figure 9, the results of a three-body simulation with pr = 1.50 and τ = 1000 orbits, which mimics the convergent migration obtained in the hydrodynamical simulation of Section 3.1.3 with Σ₀ = 8 × 10⁻⁵ (see the right panel of Figure 8). In the top-left panel of Figure 9, only the semi-major axis of the outer planet varies in the first 1500 orbits as our prescription for disk-driven migration is applied to the outer planet. Once resonant repulsion sets in near 1500 orbits, the semi-major axis of the inner planet decreases slightly and the outer planet's increases slightly. The eccentricity damping timescale here is 20 planet orbits; this rather short timescale is found to give very similar divergent evolutions in the three-body and hydrodynamical simulations (see the top-right panel in Figure 9). At 5000 orbits, the period ratio is 1.53 in the three-body run, and 1.54 in the hydrodynamical run. This small difference can be qualitatively explained by the slightly larger pace of the repulsion in the hydrodynamical simulation than modeled in the three-body run. Furthermore, the eccentricity peaks near the 2:1 and 3:2 MMRs are about a factor of two to three smaller in the three-body simulation (compare this with the eccentricities of the hydrodynamical run shown in the second upper panel of Figure 7). The reason for this difference is most likely a reflection of the longer actual circularization time in the hydrodynamical run than that adopted in the three-body integration, consistent with the operation of wake-planet repulsion, in addition to circularization-driven repulsion as discussed above.

Zoom In Zoom Out Reset image size

We also incorporate a simple model for disk dispersal, which occurs between 10⁴ and 10⁶ orbits. Since planets are assumed to be already formed at the beginning of the simulations, this range of timescales is meant to account for different planet formation timescales. The longest time (10⁶ orbits) corresponds to a case where both planet formation and capture into resonance following convergent migration occur rapidly. The shortest one (10⁴ orbits) stands for the opposite. Once disk dispersal switches on, the eccentricity damping timescale is progressively increased as an exponential growth with characteristic timescale τ_evap = 2000 orbits. This timescale is shorter than the typical viscous draining timescale of a disk inside the photoevaporation radius. As it essentially ceases to operate beyond this point, the effectiveness of disk-driven resonant repulsion is nearly entirely determined by the time from the beginning of the simulations at which dispersal switches on.

We display in Figure 10 the results of a series of three-body simulations using the model described in Section 3.2.1. Recall that M_inner = 15 M_⊕, M_outer = 13 M_⊕, and M_⋆ = M_☉, which can be thought of as being representative of Kepler's sample of multi-planetary systems (the planets' mass ratio is arbitrary, it takes the same value as in the hydrodynamical simulations in Section 3.1.3). Planets are initiated with a period ratio of 2.4. The prescription for convergent migration given at Equation (6) is used with pr = 1.50 and τ = 1000 orbits. This value of pr is meant to explore the diversity of period ratios that can be achieved by disk-driven resonant repulsion away from the 3:2 MMR. The eccentricity damping timescale varies from 20 to 200 orbits, and the time at which disk dispersal switches on from 10⁴ to 10⁶ orbits (see Section 3.2.1). Given the duration of disk-driven resonant repulsion, our results have very little dependence on the "convergent migration timescale," τ. Simulations were carried out over 1.2 × 10⁶ orbits, and the period ratio time-averaged over the last 1000 orbits of the simulations is indicated by the color bar. For a given eccentricity damping timescale, the later disk dispersal occurs, the longer planets' eccentricities are damped, and, therefore, the larger the final period ratio. At a given time prior to disk dispersal, the shorter the eccentricity damping timescale, the larger the final period ratio, again as expected. For the shortest damping timescales that we consider, period ratios can reach ∼1.75. With such short timescales, planets can evolve continuously from the 3:2 MMR to the 2:1 MMR (that is, the 2:1 resonant angles start to librate when the 3:2 resonant angles no longer librate). Had we taken even shorter damping timescales or delayed disk dispersal even more, resonant repulsion would have led to period ratios exceeding 2.

Zoom In Zoom Out Reset image size

The main conclusion that can be drawn from Figure 10 is that partial gap-opening super-Earths may experience significant disk-driven resonant repulsion away from the nominal 3:2 MMR. Under favorable circumstances (particularly efficient repulsion maintained over a particularly long timescale), period ratios can easily exceed nominal resonant values. Although we consider fixed planet masses for illustration purposes, we believe that disk-driven resonant repulsion is a generic mechanism that may have occurred among Kepler's multi-planetary systems.