© ESO, 2013

1. Introduction

An Earth- to Neptune-mass planet embedded in a gaseous protoplanetary disk excites density waves in the disk (Goldreich & Tremaine 1979), and the back-reaction of these perturbations on the planet’s orbit causes it to undergo so-called Type I migration (Ward 1997). In isothermal disks Type I migration is driven only by the differential Lindblad torque and is inward-directed and fast (Tanaka et al. 2002). In radiative disks, a torque arising from material in a planet’s horseshoe region can counteract the differential Lindblad torque such that migration can be directed either inward or outward (Paardekooper & Mellema 2006; Kley & Crida 2008). There are locations within disks for which the migration is convergent, and we refer to these as convergence zones (CZs; Lyra et al. 2010; Paardekooper et al. 2011).

At a convergence zone, the positive corotation torque is in exact balance with the negative differential Lindblad torque, such that the total torque experienced by the planet is zero. It has been proposed that CZs may act to concentrate planetary embryos and build large cores (Lyra et al. 2010; Horn et al. 2012). As planets migrate toward the CZ they become trapped in mean-motion resonances (MMRs), which prevent unlimited accretion (Morbidelli et al. 2008; Sándor et al. 2011). However, collisions do occur once embryos are trapped in a long enough resonant chain to destabilize it. In addition, turbulence may act to break resonances and enhance accretion.

Bitsch & Kley (2010) show that the corotation torque is attenuated when a planet has an eccentricity such that its orbit crosses a significant fraction of the width of the horseshoe region. When two planets become trapped in resonance due to convergent migration they excite and sustain non-zero eccentricities despite eccentricity damping by the disk (e.g. Cresswell & Nelson 2008). This in turn should decrease the magnitude of the corotation torque and thus modify the balance between the differential Lindblad and corotation torques and the subsequent migration.

Here we present simulations of the convergent migration of low-mass (M = 1−10  M_⊕) planets in idealized gas disks. We use a simple model for the influence of the eccentricity on the corotation torque. We show that planets do not actually converge at the convergence zone but are instead systematically driven inward to a shifted equilibrium position with zero net torque. The location of the equilibrium position depends on the eccentricity excitation of the planets due to mutual perturbations.

The paper is laid out as follows. In Sect. 2 we present our numerical methods and assumptions. In Sect. 3 we present results for two planets. In Sect. 4 we present results for 3−10 planets. We discuss our results in Sect. 5.

2. Methods

We use an artificial convergence zone (CZ) that mimics a mass-independent CZ (whose position does not depend on the mass of the planet) at an opacity transition, see (left panel of) Fig. 1, where there is an abrupt change from outward to inward migration (see, e.g., Masset 2011). A step function was not used because the step in the “real” profile is only due to the opacity table not being smoothed. The location of the CZ was 3 AU. Inside 3 AU the total torque is positive and equal to . Outside 3 AU the total torque is − Γ₀. Here q is the ratio between the mass of the planet and the star, h is the aspect ratio that depend on the temperature profile but is typically 0.05, Σ_p, r_p, and Ω_p are respectively the surface density, the orbital distance, and the angular speed for the planet. The total torque is a sum of the differential Lindblad torque Γ_L – assumed to remain constant – and the corotation torque Γ_C. The main interest of the artificial CZ is to get rid of the very complex shape of the real profile, to only retain the CZ, and nothing more.

Bitsch & Kley (2010) show that as the eccentricity of a body increases, the structure of its horseshoe region is modified and Γ_C decreases. We implemented a simple formula that reproduces the general influence of the eccentricity on Γ_C by a simple fit to the 3D simulations of (Bitsch & Kley 2010): (1)where x_s represents the half-width of the horseshoe region divided by the semi major axis, e is the planet’s eccentricity, and our fit produced a = 0.45, b = 3.46, and c = −2.34. We define x_s as (see Paardekooper et al. 2010, Eq. (44)): (2)where γ is the adiabatic index, q the planet-to-stellar mass ratio, and h the disk aspect ratio. The right hand panel of Fig. 1 shows that our simple equation is a good match to hydrodynamical simulations of the effect of e on Γ_C, especially at small e, although the simulation data points are sparse and appear to include some random fluctuations.

Fig. 1 Left: the total torque profile of our standard disk in units of Γ0 is shown in red. The dashed blue curve shows the true torque profile for a 10  m⊕ planet near an opacity transition, calculated using the equations in Paardekooper et al. (2011). Right: decrease in the corotation torque ΓC as a function of planetary eccentricity e. We assume that the damping (0 < D < 1) of the corotation torque as a function of eccentricity is the same for isothermal and fully radiative disc. Then, we derived D from Fig. 2 of Bitsch & Kley (2010) by taking the difference between their fully radiative and isothermal case, and normalizing to 1 the “e = 0 case”.

To perform our simulations we used a modified version of the Mercury integrator (Chambers 1999). We included both the artificial convergence zone (where Γ = Γ_L + Γ_C) and Γ_C varies with the planet’s eccentricity following Eq. (1). Eccentricity and inclination damping is done using formulas from Cresswell & Nelson (2008). We assumed the presence of a gas disk with surface density profile Σ = 500⁽r/1AU^)− 1/2g   cm^-2 (used to calculate Γ₀ and the damping).

To implement the migration caused by the disk’s torque Γ, we note that , and we define an extra acceleration a_m defined as (Cresswell & Nelson 2008, Eq. (14)) where v is the planet’s velocity, and the migration time (J is the angular momentum).

In all simulations planets started on low eccentricity and low inclination orbits. We ran each simulation for three million years with a time step between 0.4 and 3 days.

3. Two-planet case

Figure 2 shows the evolution of two 1  M_⊕ planets initially on opposite sides of the CZ at 3 AU. As they approach each other, the two planets cross a series of MMRs and are eventually trapped in the 7:6 resonance. The eccentricities of the two planets reach an equilibrium between resonant excitation and damping by the disk. This equilibrium eccentricity is 0.5 times the horseshoe semi-width x_s and damps the corotation torque to roughly 80% of its zero-eccentricity value.

The planets reach a stable orbital configuration at 1.77 and 1.96 AU, with both planets interior to the pre-imposed CZ. Given their eccentricities, the innermost planet CZ is shifted to 1.95 AU while it is 1.74 for the outer (1.96 AU) planet. In this context, the shift comes from the balance between the unchanged negative Lindblad torque and the attenuated positive corotation torque. The two-planet system stabilizes at a zero net torque point – even though each planet feels a nonzero torque and neither reaches its single-planet CZ – that is shifted inward from the nominal CZ by > 1  AU.

Fig. 2 Simulation of the convergent migration of two 1  M⊕ planets toward the CZ at 3 AU including the ΓC − e feedback illustrated in the right panel of Fig. 1.

It is clear that the planets’ eccentricities – excited by planet-planet interactions – are the key factor in determining the strength of the corotation torque and the location of the effective convergence zone. For two equal-mass planets the same qualitative behavior occurs regardless of the planet mass, although the evolution for higher mass planets is faster. In addition, the eccentricity excitation varies both with the planet mass and which resonance the planets occupy: higher eccentricities imply stronger damping of Γ_C and a closer-in effective CZ.

3.1. Effect of mass ratio

We now study the case of two unequal-mass planets. Figure 3 shows the outcome of a simple experiment in which a 10  M_⊕ planet was placed at the CZ at 3AU and a second planet was placed at 4AU. The mass of the second planet was varied in a set of simulations from 0.1 to 10  M_⊕.

In Fig. 3 the outer planet remained in 3:2 resonance. In this situation, the final location of the planets is determined by their relative masses or, for this experiment, the outer planet’s mass. The more massive the outer planet, the closer in the two planets’ effective CZ. The inner planet’s eccentricity increases for a more massive outer planet, causing a corresponding decrease in Γ_C and a more drastic inward shift. Given that each planet has a different mass and eccentricity, they have different effective CZ (one per e/x_s value). The magnitude of the overall inward shift, however, is mainly determined by the most massive planet.

Figure 3 only shows a subset of the simulations run in this numerical experiment. For higher outer planet masses, the two planets were trapped in different MMRs, resulting in discontinuities in the equilibrium positions and mean eccentricity values. However, the behavior of the two planet system remained qualitatively similar.

Fig. 3 Outcome of a set of simulations that started with a 10M⊕ planet at 3 AU and a second planet of 0.1−3M⊕ at 4 AU. The panels show the equilibrium position of the planets (top) and the eccentricity normalized to the horseshoe semi-width e/xs (bottom) as a function of the outer planet’s mass.

3.2. Resonance effect

The order of MMR is also important (for an MMR of (p + q)/p, p is the order). Two planets in 3:2 resonance will have higher eccentricities than the same planets in 11:10 resonance. The simple explanation is that lower order resonances have more frequent conjunctions, hence stronger resonant perturbations (see Murray & Dermott 2000 for more details). The MMR in which a two-planet system is captured depends on the relative migration speed and the rate of eccentricity damping (e.g., Mustill & Wyatt 2011), both determined by the disk and torque profiles, and the initial positions of the planets.

We performed a set of 100 simulations (1 Myr each) with two 3M_⊕ planets that were each initially placed randomly between 1 and 10 AU, with the same CZ at 3 AU as before. Figure 4 shows that in all cases the planets are trapped in MMRs between 11:10 and 3:2 (resonant order from 2 to 10). As expected, the excited eccentricity decreases for higher order resonances and leads to a smaller inward shift of the planetary system. The magnitude of the inward shift ranges from 0.2 to 1.5 AU. In two simulations the planets started so close to each other that they were trapped in co-orbital (1:1) resonance. In these cases the eccentricities remained very low, and both planets migrated to the nominal (classical unshifted) CZ.

Fig. 4 Final semimajor axes (top) and eccentricities (bottom) of two 3M⊕ planets trapped in different MMRs. For a 3M⊕ planet, the semiwidth of the horseshoe region is about 0.014.

Fig. 5 Two example simulations with 5 planets, one that is quite stable even if short resonance breaking occurs (top panel) and one that exhibits sustained chaotic behavior (bottom).

4. Evolution with >2 planets

We now turn our attention to the case of many planets in a disk. We ran ten simulations each with two, three, five, and ten planets of 3M_⊕ each. The planets were placed randomly between 1 and 10 AU and were given low initial eccentricities and inclinations. As before, each simulation was evolved for 3 Myr in a static disk.

For the three-planet case we find three different outcomes. First, the most likely outcome is that the three planets can be trapped in a resonant chain and migrate inward together to a zero net torque zone. This zone is typically between 2 and 2.5 AU. The eccentricities of the three planets are not the same: the middle planet is generally more excited. In the second most likely outcome, two planets are trapped in MMR and shifted inward, but a third, outer planet is too far away to be trapped in a resonant chain. The outer planet remains at the nominal CZ at 3 AU, while the inner planets end up at a new shifted zero torque zone. Finally, in some cases collisions occur between two planets and the system reverts to the two planet, unequal-mass case seen in Fig. 3.

For the case of five or ten planets the situation is more complex. The five-planet systems become trapped in chains of MMRs and migrate inward to an equilibrium position with zero net torque. However, perturbations between planets add a stochastic element to the planets’ behavior. Even the most stable systems include periods during which the planets all drift radially in the same direction. These periods are triggered by the breaking of a resonance between two planets that propagates as a perturbation through the whole system. The amplitude and frequency of these chaotic periods vary from case to case. For example, in the top simulation from Fig. 5 the resonant chain undergoes several small hiccups but the drift is small compared with the typical inter-planetary spacing. In contrast, the perturbations in the simulation in the bottom panel of Fig. 5 are much stronger. For instance, we consider the interval between roughly 1.1 and 1.3 Myr in the bottom case. At 1.12 Myr the two outer planets become trapped in 4:3 MMR. They migrate inward due to the resonant eccentricity increase. This perturbation propagates to the inner system, and the eccentricity increase causes the planets to migrate inward all together. Five thousand years later, the outer planets leave and then re-enter the 4:3 MMR, again perturbing the system. Finally, at 1.13 Myr the outer two planets leave the 4:3 MMR for good. As they leave the MMR, the two planets act like isolated bodies and migrate toward the nominal convergence zone with near zero eccentricities. Without the net negative torque supplied by the outer two planets, the inner planets reactively migrate outward toward an effective three-planet-system zero torque zone. This is the cause of the abrupt outward migration of all five planets. However, the two outer planets quickly find themselves in the 5:4 MMR. Over the subsequent 0.15 Myr interval, they are periodically trapped in 5:4 MMR but outward migration continues because most of the time they are not in MMR and their eccentricity is low. The system-wide outward migration stops at 1.335 Myr when the outer planets cross over the 5:4 MMR and are trapped in 6:5 MMR. This configuration stabilize the system, excites the outer planets’ eccentricities, and causes the whole system to migrate inward again, which signals the end of this particular metaperturbation.

The rest of the evolution is composed of the same kinds of perturbations. Perturbations come from planets entering or leaving resonances, then propagate between planets. When leaving a resonance, planets’ eccentricities usually drop, causing outward migration. Conversely, being trapped in resonance excites a sustained eccentricity that leads to inward migration. Via these perturbations, entire systems of planets undergo modest-scale chaotic migration due to the difficulty of maintaining resonant chains over long times. Each of the five planet systems we simulated remained stable, in that no collisions occurred, but the amplitude of the chaotic migration varied from case to case; the two examples from Fig. 5 show the extremes. The ten-planet simulations were even more chaotic, and collisions did occur.

The most important factor in determining the amplitude of chaotic oscillations is the order of the resonances. Low-order resonances maintain higher eccentricities and are less stable because they are sensitive to eccentricity variations. On the other hand, high-order resonances maintain lower eccentricities, and eccentricity perturbations are less important. For instance, in the system in the top panel of Fig. 5, perturbations are quickly smoothed out, whereas in the system in the bottom panel the frequency of perturbations is high enough that the system does not tend to a stably evolving configuration. Therefore, in this context a more compact system of planets is more stable than an extended system, which is the opposite of the pure gravitational situation (Marchal & Bozis 1982).

5. Discussion

We have shown that planets do not actually converge in convergence zones (CZs). Instead, embryos rapidly migrate toward the CZ and are trapped in chains of MMRs. This causes their eccentricities to increase and remain high enough to attenuate the corotation torque. The zero-torque equilibrium point of the resonant chain of planets is determined by the sum of each individual planet torque (sum of attenuated corotation torque and unattenuated differential Linbdlad torque). In practice, the effective zero net torque zone is shifted inward and is most strongly determined by the CZ of the most massive planet in the resonant chain. It is not a true CZ because each planet feels a different CZ (depending on its eccentricity).

The inward shift exists because planets’ eccentricities are sustained by resonant perturbations. The amplitude of a planet’s eccentricity is a result of the competition between resonant excitation and eccentricity damping by the disk. For sufficiently high sustained eccentricities, an entire resonant chain of planets can migrate all the way in to the inner edge of the disk. Changing the disk properties may thus change the typical sustained eccentricities by affecting the eccentricity damping timescale. However, changing the disk properties affects other properties of the system such as the torque profile. In changing the properties of the disk, it is not straightforward to assess the effect on the planets’ evolution, as planets may be trapped in different resonances with different eccentricity excitations and even different stability criteria.

The CZ depends on the disk parameters such as the viscosity, temperature and surface density profiles (e.g. Paardekooper et al. 2011). Here we have used a disk profile that is motivated by complex models but is nonetheless artificial. Even if the results shown use a particular model, they are robust when we vary the torque shape in function of the orbital distance, as long as a CZ zone exists. In a more realistic disk we would expect a few differences. First, there may be multiple CZs in the same disk owing to different effects (Lyra et al. 2010; Hasegawa & Pudritz 2011). Second, mass-dependent convergence zones may exist in the outer part of the disk, where this mechanism should be less efficient because planets do not migrate to the same location in the disk and may not produce the resonant chains essential for our mechanism to occur (Cossou et al., in prep.). Third, as the disk dissipates, so does the torque profile and the location of CZs (Lyra et al. 2010; Horn et al. 2012). Finally, turbulence is thought to be common in real protoplanetary disks (Armitage 2011). Although turbulence does not affect the long-term evolution of a single planet in a radiative disk (Pierens et al. 2012), we expect it to modify resonant capture and eccentricity evolution (see Pierens et al. 2011).

Acknowledgments

We thank B. Bitsch, F. Selsis, F. Hersant, and an anonymous referee for stimulating discussions. We acknowledge funding from the CNRS’s PNP program and the Conseil Régional d’Aquitaine. Computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l’Adour.

References

All Figures

Fig. 1 Left: the total torque profile of our standard disk in units of Γ0 is shown in red. The dashed blue curve shows the true torque profile for a 10  m⊕ planet near an opacity transition, calculated using the equations in Paardekooper et al. (2011). Right: decrease in the corotation torque ΓC as a function of planetary eccentricity e. We assume that the damping (0 < D < 1) of the corotation torque as a function of eccentricity is the same for isothermal and fully radiative disc. Then, we derived D from Fig. 2 of Bitsch & Kley (2010) by taking the difference between their fully radiative and isothermal case, and normalizing to 1 the “e = 0 case”.

In the text

	Fig. 2 Simulation of the convergent migration of two 1  M⊕ planets toward the CZ at 3 AU including the ΓC − e feedback illustrated in the right panel of Fig. 1.
In the text

	Fig. 3 Outcome of a set of simulations that started with a 10M⊕ planet at 3 AU and a second planet of 0.1−3M⊕ at 4 AU. The panels show the equilibrium position of the planets (top) and the eccentricity normalized to the horseshoe semi-width e/xs (bottom) as a function of the outer planet’s mass.
In the text

	Fig. 4 Final semimajor axes (top) and eccentricities (bottom) of two 3M⊕ planets trapped in different MMRs. For a 3M⊕ planet, the semiwidth of the horseshoe region is about 0.014.
In the text

	Fig. 5 Two example simulations with 5 planets, one that is quite stable even if short resonance breaking occurs (top panel) and one that exhibits sustained chaotic behavior (bottom).
In the text