Highly inclined and eccentric massive planets
II. Planetplanet interactions during the disc phase
^{1} naXys, Department of Mathematics, University of Namur, 8 Rempart de la Vierge, 5000 Namur, Belgium
email: sotiris.sotiriadis@unamur.be
^{2} Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, 22100 Lund, Sweden
^{3} Laboratoire Lagrange, Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Boulevard de l’Observatoire, CS34220, 06304 Nice Cedex 4, France
^{4} Institut Universitaire de France, 103 Boulevard Saint Michel, 75005 Paris, France
Received: 9 March 2016
Accepted: 26 November 2016
Context. Observational evidence indicates that the orbits of extrasolar planets are more various than the circular and coplanar ones of the solar system. Planetplanet interactions during migration in the protoplanetary disc have been invoked to explain the formation of these eccentric and inclined orbits. However, our companion paper (Paper I) on the planetdisc interactions of highly inclined and eccentric massive planets has shown that the damping induced by the disc is significant for a massive planet, leading the planet back to the midplane with its eccentricity possibly increasing over time.
Aims. We aim to investigate the influence of the eccentricity and inclination damping due to planetdisc interactions on the final configurations of the systems, generalizing previous studies on the combined action of the gas disc and planetplanet scattering during the disc phase.
Methods. Instead of the simplistic Kprescription, our Nbody simulations adopt the damping formulae for eccentricity and inclination provided by the hydrodynamical simulations of our companion paper. We follow the orbital evolution of 11 000 numerical experiments of three giant planets in the late stage of the gas disc, exploring different initial configurations, planetary mass ratios and disc masses.
Results. The dynamical evolutions of the planetary systems are studied along the simulations, with a particular emphasis on the resonance captures and inclinationgrowth mechanisms. Most of the systems are found with small inclinations (≤ 10°) at the dispersal of the disc. Even though many systems enter an inclinationtype resonance during the migration, the disc usually damps the inclinations on a short timescale. Although the majority of the multiple systems in our simulations are quasicoplanar, ~5% of them end up with high mutual inclinations (≥ 10°). Half of these highly mutually inclined systems result from two or threebody meanmotion resonance captures, the other half being produced by orbital instability and/or planetplanet scattering. When considering the longterm evolution over 100 Myr, destabilization of the resonant systems is common, and the percentage of highly mutually inclined systems still evolving in resonance drops to 30%. Finally, the parameters of the final system configurations are in very good agreement with the semimajor axis and eccentricity distributions in the observations, showing that planetplanet interactions during the disc phase could have played an important role in sculpting planetary systems.
Key words: planets and satellites: formation / planetdisk interactions / planets and satellites: dynamical evolution and stability
© ESO, 2017
1. Introduction
The surprising diversity of the orbital parameters of exoplanets shows that the architecture of many extrasolar systems is remarkably different from that of the solar system and sets new constraints on planet formation theories (see Winn & Fabrycky 2015 for a review). The broader eccentricity distribution of the detected planets, the existence of giant planets very close to their parent star (“hot Jupiters”) and the strong spinorbit misalignment of a significant fraction (~ 40%) of them (Albrecht et al. 2012) are several characteristics that appear to be at odds with the formation of the solar system. Furthermore, some evidence of mutually inclined planetary orbits exists; for instance, the planets υ Andromedae c and d whose mutual inclination is estimated to 30° (Deitrick et al. 2015).
It seems that the interactions between the planets and their natal protoplanetary disc, at the early stages of formation, play an important role in sculpting a planetary system. The angular momentum exchange between the newborn planets (or planetary embryos) and the gaseous disc tends to shrink the semimajor axis and keeps them on coplanar and nearcircular orbits (Cresswell et al. 2007; Bitsch & Kley 2010, 2011). If the planets are massive enough to carve a gap and push the material away from their orbit, they migrate, generally inwards, with a timescale comparable to the viscous accretion rate of the gas towards the star, a phenomenon referred to as TypeII migration (Lin & Papaloizou 1986; Kley 2000; Nelson et al. 2000; and Baruteau et al. 2014 for a review). The gas also affects the eccentricities and inclinations of the embedded planets (Goldreich & Tremaine 1980; XiangGruess & Papaloizou 2013; Bitsch et al. 2013).
Several authors, using analytical prescriptions, have investigated the impact of the planetdisc interactions on the orbital arrangement of planetary systems through threedimensional (3D) Nbody simulations. Studies in the context of TypeII migration have been performed for twoplanet systems (Thommes & Lissauer 2003; Libert & Tsiganis 2009; Teyssandier & Terquem 2014) and threeplanet systems (Libert & Tsiganis 2011b). A notable outcome of these works is that, during convergent migration, the system can enter an inclinationtype resonance that pumps the mutual inclination of the planetary orbits.
After the dispersal of the disc, dynamical instabilities among the planets can result in planetplanet scattering and a subsequent excitation of eccentricities and inclinations (Weidenschilling & Marzari 1996; Rasio & Ford 1996; Lin & Ida 1997). Such a mechanism has been proposed to explain the eccentricity distribution of giant extrasolar planets (Marzari & Weidenschilling 2002; Chatterjee et al. 2008; Jurić & Tremaine 2008; Ford & Rasio 2008; Petrovich et al. 2014). The planetplanet scattering model is based on unstable initial conditions that do not take into account the imprint of the disc era (Lega et al. 2013).
More realistic approaches study the combined action of both previous mechanisms (TypeII migration and scattering). Parametric explorations via Nbody simulations (Adams & Laughlin 2003; Moorhead & Adams 2005; Matsumura et al. 2010; Libert & Tsiganis 2011a) have shown that dynamical instabilities can occur during the disc phase.
In Nbody simulations, a damping prescription depending on a scaling factor is commonly used to mimic the influence of the disc on the eccentricities of the planets, which is usually a firstorder approximation for the eccentricity damping timescale (ė/e = −K  ȧ/a ). The considered Kfactor, currently not well determined but usually estimated in the range 1−100, determines the final eccentricity and inclination distributions (see discussion in Sect. 2.3). Moreover, no damping on the planet inclination is generally included in these studies (e.g. Thommes & Lissauer 2003; Libert & Tsiganis 2009, 2011a). On the other hand, hydrodynamical simulations accurately modelize planetdisc interactions. Several works have investigated the evolution of the eccentricities of giant planets in discs using twodimensional (2D) hydrodynamical codes (e.g. Marzari et al. 2010; Moeckel & Armitage 2012; Lega et al. 2013). The influence of the disc on planetary inclinations has been studied in our accompanying paper (Bitsch et al. 2013, hereafter denoted Paper I), where we performed 3D hydrodynamical numerical simulations of protoplanetary discs with embedded high mass planets, above 1 M_{Jup}, computing the averaged torques acting on the planet over every orbit. As there is a prohibitive computational cost to following the orbital evolution for timescales comparable with the disc’s lifetime (~ 1−10 Myr) and, as a consequence, to generating large statistical studies, we have derived an explicit formula for eccentricity and inclination damping, suitable for Nbody simulations, as a function of eccentricity, inclination, planetary mass, and disc mass.
In this paper, we aim to improve the previous Nbody studies combining planetdisc interactions and planetplanet interactions by adopting the damping formulae of Paper I. In other words, we combine the speed and efficiency of an Nbody integration with a symplectic scheme, and an improved modelization of the gas effect promoted by hydrodynamical simulations. More specifically, we extend the work of Libert & Tsiganis (2011a), focusing on the orbital evolution of three giant planets in the late stage of the gas disc. In particular, we add inclination damping, which was not taken into account in their previous work; thus, comparing our results with theirs allows us to highlight the role of this phenomenon. We have performed 11 000 numerical simulations, exploring different initial configurations, planetary mass ratios and disc masses. Our goal is to investigate the influence of the eccentricity and inclination damping due to planetdisc interactions on the final configurations of planetary systems and their inclination distribution. Through our parametric analysis, we also discuss the most common threebody resonance captures occuring during the migration of the planets, as well as the most frequent mechanisms producing inclination increase.
The paper is organized as follows. In Sect. 2, we describe the setup of our numerical experiments. Typical dynamical evolutions of planetary systems during the disc phase are analyzed in detail in Sect. 3, while the orbital parameters and resonance configurations of our set of planetary systems at the dispersal of the disc are presented in Sect. 4. In Sect. 5, we address the question of the stability of the planetary systems formed in our simulations, by analyzing the effect of the longterm evolution on the architecture of the systems. Finally, our conclusions are given in Sect. 6.
2. Setup of the Nbody simulations
2.1. Physical modeling
In this work, we analyze the evolution of threeplanet systems embedded in a protoplanetary disc and evolving around a solarmass star. We consider systems with three fully formed gaseous giant planets that do not accrete any more material from the disc as they migrate inwards, towards the star. We use the symplectic integrator SyMBA (Duncan et al. 1998), which allows us to handle close encounters between the bodies by using a multiple timestep technique. Moreover, due to the interaction with the disc and the subsequent migration, some systems end up in a compact stable configuration very close to the star. In these systems, planets might have either low or high eccentricity values. For this reason, we adopt a symplectic algorithm that also has the desirable property of being able to integrate close perihelion passages with the parent star (Levison & Duncan 2000). In order to ensure high resolution for closein orbits, we set the timestep of our simulations at dt = 0.001 yr.
The migration of the planets in TypeII migration, due to angular momentum exchange with the disc, is on a similar timescale as the viscous accretion time, , where ν is the kinematic viscosity and r_{pl} is the radius of the planetary orbit. However, when the mass of the planet is comparable to the mass of the material in its vicinity, the migration rate scales with the ratio of the planetary mass over the local disc mass (Ivanov et al. 1999; Nelson et al. 2000; Crida & Morbidelli 2007). Thus, the rate for TypeII regime consists of two different cases, the discdominated case and the planetdominated case: (1)where α is the ShakuraSunyaev viscosity parameter (Shakura & Sunyaev 1973), h the disc aspect ratio, Ω the orbital frequency of the planet and is the local mass of the disc. The rate of TypeII is still under debate and may depart from the viscous time (Hasegawa & Ida 2013; Duffell et al. 2014; Dürmann & Kley 2015). However, the scope of the present paper is not to study the TypeII regime (which would require hydrodynamical simulations), but planetplanet interactions during planetary migration. Teyssandier & Terquem (2014) have pointed out that the eccentricity and inclination evolution (which is of interest for us here) of migrating giant planets in the 2:1 MMR is not affected by the TypeII timescale. Therefore, we always use Eq. (1) (like most previous similar studies), independent of e and i of the planets.
The disc parameters are set to the classical values α = 0.005 and h = 0.05. Four values of disc mass are considered in this study, being 4,8,16 and 32 M_{Jup}, in order to verify the robustness of our results. The initial surface density profile is Σ ∝ r^{0.5} and the disc’s inner and outer edges are set to R_{in} = 0.05 AU and R_{out} = 30 AU. We apply a smooth transition in the gasfree inner cavity, following Matsumoto et al. (2012), by using a hyperbolic tangent function, , where Δr = 0.001 AU. Following Libert & Tsiganis (2011a), we apply an inward migration to the outer planet only. We start the evolution of the system with fully formed planets and consider that there is not enough gas to drive migration between them^{1}. To fully mimic the interactions with the disc, we use the damping formulae of Paper I for eccentricity and inclination (see Sects. 2.2 and 2.3 for their full description).
The effects of orbital migration, eccentricity, and inclination damping are added to the acceleration of the planets, as in Papaloizou & Larwood (2000): (2)where τ_{ecc} and τ_{inc} denote the timescales for eccentricity and inclination damping, respectively (see Eqs. (3) and (4)). These rates are computed in every step of the integration and depend on the local surface density of the disc, the mass of the planet, and the eccentricity and inclination of the planet at the considered time. The damping forces were implemented in SyMBA in the same way as in Lee & Peale (2002; see their appendix for a detailed description of the method), such that it keeps the symmetry of the symplectic algorithm.
The decrease of the gas disc is implemented in two ways. In the constantmass model (hereafter denoted CM), the mass of the circumstellar disc is kept constant for 0.8 Myr, then the gas is instantly removed. We use this approach to understand, through statistical analyses, whether the inclination damping has a strong impact on the final configurations of planetary systems. A more realistic modelization of protoplanetary disc is the decreasingmass model (DM), where we decrease its mass exponentially through the evolution of the system, with a dispersal time of ~ 1 Myr (e.g., Mamajek 2009). More specifically, we neglect the interaction with the disc when dM_{Disc}/ dt < 10^{9}M_{star}/ yr. As a consequence, the surface density also evolves with time in the DM model.
During the evolution of a system, close encounters between the bodies can lead to planetplanet collisions, ejections from the system and collisions with the parent star. We treat a possible merge between two bodies as a totally plastic collision, when their distance becomes less than the sum of the radii of the planets. We assume Jupiter’s density here when computing the radius from the mass. The boundary value for planet accretion onto the star is 0.02 AU and the one for ejection from the system, 100 AU.
Fig. 1 Normalized initial and final mass distributions of our simulations compared with the massdistribution (in the range [0.65,10] M_{Jup}) of the observed giant planets with a > 0.1. 

Open with DEXTER 
Initial disc mass, integration time, and number of simulations for each disc modelization.
In our simulations, gas giant planets are initially on coplanar and quasicircular orbits (e ∈ [0.001,0.01] and i ∈ [0.01°,0.1°]). The initial semimajor axis of the inner body is fixed to 3 AU and the middle and outer ones follow uniform distributions in the intervals a_{2} ∈ [4.3,6.5] AU and a_{3} ∈ [8,16] AU, respectively. These initial semimajor axis distributions provide a wide range of possible resonance captures. We choose randomly initial planetary masses from a loguniform distribution that approximately fits with the observational data^{2} in the interval [0.65,10]M_{Jup} (black and blue curves in Fig. 1). In Table 1, we describe the different sets of simulations we performed for both disc modelizations and different disc masses. In total, 4000 simulations were realized for the CM model and 7000 for the DM model. The computational effort required for our investigation was ~1.8 × 10^{5} computational hours.
2.2. Eccentricity and inclination damping formulae
The damping rates for eccentricity and inclination used in the present work originate from the accompanying Paper I, where the evolutions of the eccentricity and inclination of massive planets (≥ 1 M_{Jup}) in isothermal protoplanetary discs have been studied with the explicit/implicit hydrodynamical code NIRVANA in 3D. The forces from the disc, acting onto the planet kept on a fixed orbit, were calculated, and a change of de/ dt and di/ dt was thereby determined.
Eccentricity and inclination are mostly damped by the interactions with the disc. For highly inclined massive planets, the damping of i occurs on smaller timescales than the damping of e. The only exception is lowinclination planets with a sufficient mass (≥ 4−5 M_{Jup}), for which the interactions of the planet with the disc result in an increase of the planet eccentricity. As a result, for singleplanet systems, the dynamics tend to damp the planet towards the midplane of the disc, in a circular orbit in the case of a lowmass planet and in an orbit whose eccentricity increases over time due to the interactions with the disc in the case of a highmass planet and a sufficiently massive gas disc. The increase of the planetary eccentricity for massive planets may appear surprising, but it is a wellestablished phenomenon, further discussed in Sect. 2.3.
The damping formulae depend on the planet mass M_{P} (in Jupiter masses), the eccentricity e_{P} and the inclination i_{P} (in degrees). The eccentricity damping function is given by (3)where i_{D} = M_{p}/ 3 degrees, M_{⋆} is the mass of the star in solar mass, and with coefficients
The damping function for inclination is given, in degrees per orbit, by (4)with coefficients
Let us highlight the limitations of the implemented formulae. These expressions have been derived from fitting the results of hydrodynamical simulations for planets between 1 and 10 M_{Jup} (to cover the range of giant planets) with an eccentricity smaller than 0.65 (the expression for coefficient c_{i} is clearly not valid for e_{P} > 2/3). Practically, when the eccentricity of a planet exceeds the limit value during the integration, the square root of a negative number must be computed, so we instead give a small number (10^{5}) to the problematic factor (2−3e_{P}) of the inclination damping formula. Moreover, systems with an overly massive planet (>10M_{Jup}), and following a merging event, will be disregarded from our parametric analysis.
An additional limitation concerning the planetary inclinations is that it is rather unclear if an inclined planet will continue to migrate on a viscous accretion timescale. Thommes & Lissauer (2003) stated that the interaction of a planet away from the midplane with the gas disc is highly uncertain. Furthermore, it has been shown that Kozai oscillations with the disc govern the longterm evolution of highly inclined planets (e.g., Terquem & Ajmia 2010; Teyssandier et al. 2013; Paper I), which makes it hard to predict the exact movement of the planet and its orbital parameters at the dispersal of the disc. We note that the study of the exact TypeII timescales due to planetdisc interactions is beyond the scope of the present work; therefore, for simplicity and despite these limitations, the viscous accretion timescale is adopted independently of the planetary inclination in our simulations. However, as we will see in Sects. 4.3 and 4.7, reaching such high inclination during the disc phase is rather occasional, starting from a coplanar system in interaction with the disc. We note also that no inclination damping is applied for planets with inclinations below 0.5°.
In Eqs. (3) and (4), the damping (and excitation) formulae scale with the mass of the disc. It should be noted that the total mass of the disc is actually an irrelevant parameter: whatever the power law chosen, the integration of Σ(r) from r = 0 to + ∞ diverges. Hence, the total mass in the grid actually depends at least as much on the boundaries as on the local surface density. But the damping is done locally by the disc in the neighborhood of the planet. Therefore, the M_{disc} parameter in Eqs. (3) and (4) is actually the disc mass between 0.2 and 2.5 a_{Jup}. It should be noted that this was not the total mass of gas included in the grid, because the grid was extended to 4.2 a_{Jup}. In the present paper, for convenience, the notation M_{Disc} refers to the total mass of the disc, instead of the local one. Therefore, we have scaled the damping formulae accordingly, taking M_{disc} = M_{Disc} × (2.5/4.2)^{3/2} = 0.4 M_{Disc}.
The previous formulae have been derived from hydrodynamical simulations with fixed disc parameters; for instance, the shape of the density profile, the viscosity, and the aspect ratio. A change of values of these parameters in the present study would not be consistent. We do not aim to realize an exhaustive study of the free parameters of the hydrodynamical simulations, but, with realistic parameters, analyze in detail the dynamical interactions between the planets in the disc.
The final state of singleplanet systems in a protoplanetary disc has been studied in detail in Paper I. In the present work, we consider systems of multiple planets that excite each other’s inclination during the migration in the gas disc while suffering from the strong damping influence of the disc. When two planets share a common gap (i.e., mutual distance <4 R_{Hill}), we reduce the effect of the damping by a factor of 2, because each planet only interacts with either the inner or the outer disc, that is half the mass of gas a single planet would interact with. For reference, this applies for a long time scale in only 2% of our simulations, since the planets are initially on wellseparated orbits and never evolve in a resonance with close enough orbital separation to share a common gap according to our criterion. When three planets share a gap (which applies for a long time scale in only <1% of our simulations), the middle one is hardly in contact with the gas disc; therefore, we reduce the damping rates by a factor of 100 for the middle planet. Choosing 100 is relatively arbitrary, but the basic idea is that the middle planet’s eccentricity and inclination are almost not damped by the inner and outer disc because they are too far away, and the gas density inside the gap is at most a hundredth of the unperturbed surface density. Taking 50 or 500 would not significantly change the evolution of the system: the middle planet is excited freely by the two others.
2.3. Discussion on damping timescales
We stress that the final results of our simulations depend on the recipe chosen for the damping of the eccentricity and inclination of the planets. Here we review the benefits and limitations of the Kprescription and the prescription of Paper I.
In Nbody simulations, a Kprescription for the eccentricity damping timescale, (5)is commonly used to mimic the influence of the disc on the eccentricities of the planets. For instance, Lee & Peale (2002) reproduced the GJ 876 system using a fixed ratio K = 100 when migration and eccentricity damping are applied to the outer planet only, and K = 10 when applied to the two planets. Concerning the formation of inclined systems, Thommes & Lissauer (2003), Libert & Tsiganis (2009), and Teyssandier & Terquem (2014) performed extensive parametric studies for different K values and orbital parameters of the planets. They showed that planets captured in a meanmotion resonance during TypeII migration can undergo an inclinationtype resonance, if eccentricity damping is not too efficient. Libert & Tsiganis (2011a,b) generalized these results to threeplanet systems, studying the establishment of threeplanet resonances (similar to the Laplace resonance in the Galilean satellites) and their effects on the mutual inclinations of the orbital planes of the planets. These works concluded that the higher the value of K, the smaller the final eccentricities and the less the number of mutually inclined systems.
However, since the actual value of K is estimated in the wide range 1−100, no eccentricity or inclination distributions can be inferred from these studies, the final system configurations being strongly dependent on the value of K considered in the simulations. Moreover, the Kprescription approach is too simplistic, as highlighted in several works. For instance, Crida et al. (2008) showed that for gap opening planets, the ratio K is not constant, but depends on the eccentricity of the planet (see their Fig. 3, fourth panel). Note also that the Kprescription has been developed by Tanaka & Ward (2004), using firstorder, linear calculations which apply to lowmass planets only. For giant planets, it appears more appropriate to use another recipe for eccentricity damping, such as the one provided in Paper I, based on hydrodynamical simulations.
In hydrodynamical simulations (included in Paper I), planetary eccentricity growth has sometimes been observed. This is a result of the disc becoming eccentric, and is a longterm effect. The simulations of Kley & Dirksen (2006) stated that a 3Jupiter mass planet can excite the eccentricity of the disc. Reciprocally, planetary eccentricities could be excited by the disc in some circumstances (Papaloizou et al. 2001; Goldreich & Sari 2003; Ogilvie & Lubow 2003). Simulations by D’Angelo et al. (2006) investigated the longterm evolution of Jupiterlike planets in protoplanetary discs and found eccentricity growth for planets even smaller than 3Jupiter masses on timescales of several thousand orbits. The growth of the eccentricity happens when the 2:4, 3:5, or 4:6 outer Lindblad resonance becomes dominant (Teyssandier & Ogilvie 2016; Duffell & Chiang 2015). This process therefore depends on the depth and width of the gap in the disc caused by the planet, hence on the planet mass and disc parameters (e.g., Crida et al. 2006). Duffell & Chiang (2015) computed a gap opening parameter from works of different authors, including Paper I, and showed that in all these studies, eccentricity damping is similarly observed for and eccentricity growth for (their Fig. 9).
The formulae for eccentricity and inclination damping derived in Paper I and used in the present work depend on the eccentricity, inclination, and mass of the planet, as well as on the migration rate via the surface density profile Σ (local mass disc M_{disc} in Eqs. (3) and (4)). Previous works using the Kprescription argued that the final system configurations are strongly related to the damping chosen (e.g., Thommes & Lissauer 2003; Libert & Tsiganis 2009, 2011b; Teyssandier & Terquem 2014). Instead of the Kfactor dependency, the key parameter of the formulae used in this work is M_{Disc}. To check the robustness of our results, we consider four different disc masses in our simulations, being 4,8,16, and 32 M_{Jup}, or equivalently four different migration timescales (see Eq. (1)).
For consistency, let us give an evaluation of the (initial) ratio of the migration timescale to the eccentricity damping timescale (τ_{II}/τ_{ecc}) observed in our simulations. The eccentricity damping is in the order of 10^{4}−10^{5} yr (depending on the mass of the planet, see Fig. 4 of Paper I), while TypeII migration timescale is of the order of 10^{5} yr (see Eq. (1)), which would result in a K value of 1−10. This is smaller than the K value assumed in Lee & Peale (2002), but in good agreement with measures by Crida et al. (2008). Note that a preliminary study, with a similar threeplanet initial setup as in the present work but making use of the Kprescription for the eccentricity damping (no inclination damping considered), was achieved in Libert & Tsiganis (2011b) for this range of K values. Inclination damping timescale τ_{inc} is in the same range as τ_{ecc}, depending again on the mass of the planet. So the ratio τ_{II}/τ_{inc} is in the range 1−10 throughout the evolution of the system.
Fig. 2 Inclination excitation due to planetplanet scattering (DM model). After the hyperbolic ejection of the inner planet, the system consists of two planets with a large orbital separation, large eccentricities and slightly inclined orbits. The planetary masses are m_{1} = 1.98, m_{2} = 8.35, m_{3} = 3.8 M_{Jup}. 

Open with DEXTER 
Furthermore, let us note that while the damping on the inclinations is either not considered or chosen in an arbitrary way in previous Nbody studies, the formulae for eccentricity and inclination damping used in this work have been derived consistently. In particular, it means that no additional free parameter has to be introduced for the inclination damping (as it would be the case for the Kprescription where two Kfactors should be considered). For this reason, the formulae of Paper I are appropriate for a study of the impact of the inclination damping on the formation of noncoplanar systems. In the following sections, we identify the dynamical mechanisms producing inclination increase and describe the eccentricity and inclination distributions of planetary systems obtained when adopting these damping formulae.
Fig. 3 Inclination excitation due to inclinationtype resonance in a Laplace 1:2:4 meanmotion resonance (DM model). The planetary masses are m_{1} = 2.66, m_{2} = 1.44 and m_{3} = 0.73 M_{Jup}. 

Open with DEXTER 
3. Typical dynamical evolutions
This section shows two typical evolutions of our simulations in the DM model, characterized by inclination excitation. They illustrate well the constant competition between the planetplanet interactions that excite the eccentricity and inclination of giant planets on initially circular and coplanar orbits, and the planetdisc interactions that damp the same quantities. Two dynamical mechanisms producing inclination increase during the disc phase are discussed in the following: planetplanet scattering and inclinationtype resonance. The frequency of these dynamical mechanisms will be discussed in Sect. 4.7.
The planetplanet scattering mechanism for planetary masses of m_{1} = 1.98, m_{2} = 8.35 and m_{3} = 3.8 M_{Jup} is illustrated in Fig. 2. While the outer planet migrates, m_{2} and m_{3} are captured in a 1:3 meanmotion resonance, at approximately 2 × 10^{5} yr, as indicated by the libration of both resonant angles θ_{1} = λ_{2}−3λ_{3} + 2ω_{2} and θ_{2} = λ_{2}−3λ_{3} + 2ω_{3}, with λ being the mean longitude and ω the argument of the pericenter. Following the resonance capture, the eccentricities are excited and the two planets migrate as a pair. When m_{2} approaches a 3:7 commensurability with the inner planet, the whole system is destabilized, leading to planetplanet scattering. Indeed, heavy planets usually lead to dynamical instability before the establishment of a threebody meanmotion resonance (Libert & Tsiganis 2011b). The inner less massive body is then ejected from the system at 4 × 10^{5} yr, leaving the remaining planets in (slightly) inclined orbits with large eccentricity variations and large orbital separation. As shown here, an excitation in inclinations can result from the planetplanet scattering phase. Let us note here that only during the brief chaotic phase before the ejection of the inner planet, did all the bodies share a common gap, and a reduced rate for eccentricity and inclination damping was applied (see Sect. 2.2).
Figure 3 shows an example of the establishment of a threebody resonance and the subsequent inclination excitation. In this example, the three planetary masses are m_{1} = 2.66 M_{Jup}, m_{2} = 1.44 M_{Jup} and m_{3} = 0.73 M_{Jup}. The two inner planets are initially in a 1:2 meanmotion resonance. The convergent migration of the outer planet leads to the establishment of a threebody meanmotion resonance, the Laplace 1:2:4 resonance (since n_{2}:n_{1} = 1:2 and n_{3}:n_{2} = 1:2, the multipleplanet resonance is labeled as n_{3}:n_{2}:n_{1} = 1:2:4), at approximately 2 × 10^{4} yr. The libration of the resonant angle φ_{L} = λ_{1}−3λ_{2} + 2λ_{3} is shown in the bottom panel of Fig. 3. From then on, the three planets migrate together while in resonance (they never share a common gap) and their eccentricities rapidly increase (up to 0.4 for e_{2}), despite the damping exerted by the disc. Let us recall that there is no eccentricity excitation from the disc for lowmass planets. When the eccentricities are high enough, the system enters an inclinationtype resonance, at approximately 1 × 10^{5} yr: the angles and (where Ω is the longitude of the ascending node) start to librate and a rapid growth of the inclinations is observed. However, the strong damping exerted by the disc on the planets leads the planets back to the midplane and the system is temporarily outside the inclinationtype resonance. On the other hand, the eccentricities of the planets keep increasing as migration continues, in such a way that the system reenters the inclinationtype resonance at ~3 × 10^{5} yr, with both critical angles and in libration and a new sudden increase in inclination. This time the inclination values is maintained for a long time, since the damping exerted by the gas disc is weaker due to the exponential decay of the disc mass. On the contrary, the weak inclination damping contributes in an appropriate way to the longterm stability of the planetary system by keeping the eccentricity and inclination values constant. This example shows that resonant dynamical interactions between the planets during the disc phase can also form noncoplanar planetary systems. However, when the capture in the inclinationtype resonance happens very rapidly, the gas disc is so massive (discdominated case) that it forces the planet to get back to the midplane in a very short timescale (~5 × 10^{4} yr). This is the case for the large majority of our simulations (see Sect. 5).
4. Parameter distributions and resonance configurations at the dispersal of the disc
This section describes the main characteristics of the systems formed by our scenario combining migration in the disc and planetplanet interactions. On the one hand, we aim to see whether or not the parameters of the systems formed in our simulations are consistent with the observations, showing how important the role of planetplanet interactions during the disc phase is in the formation of planetary systems. On the other hand, our objective is to study the impact of the inclination damping on the final configurations (Sects. 4.1 and 4.2), especially on the formation of noncoplanar systems and the resonance captures (Sects. 4.3–4.6).
The systems are followed on a timescale slightly longer than the disc’s lifetime, that is 1.2 × 10^{6} yr for the 4000 simulations of the CM model and 1.4 × 10^{6} yr for 5800 simulations of the DM model (see Table 1). This timescale is appropriate for the study of the influence of the discplanet and planetplanet interactions during the disc phase. The systems emerging from the disc phase would most probably experience orbital rearrangement in the future due to planetplanet interactions. The issue of the longterm evolution and stability of the systems will be addressed in Sect. 5.
Let us begin with an overview of our simulations. As the eccentricities of the planets increase due to resonance capture, most of the systems evolve towards close encounters between the bodies, leading to ejections, mergings and/or accretions on the star. The percentages of these three possible outcomes are given in Table 2, for both disc models. In the DM model, ejections are preferred over collisions among the planets and with the central star.
Table 3 shows the final number of planets at the end of the integration time. The constantmass model preferably forms systems with a small number of planets, since migration is efficient during all the simulations (discdominated case favored) and pushes the planets closer to one another and closer to the parent star. However, the more realistic decreasingmass scenario shows that systems with two or three planets are the most usual outcomes: 50% of our systems consist of two planets at the end of the integration time and 40% of the systems are still composed of three planets.
Percentages of discard events for both disc models at the dispersal of the disc.
Number of planets in the systems of our simulations at the end of the integration time, for both disc models.
Fig. 4 Normalized semimajor axis distributions in both disc modelizations for all the initial disc masses (blue and green lines), and in the observed giant planet population (with a> 0.1 AU and M_{p} ∈ [0.65,10] M_{Jup}, black line). The bin size is Δlog (a) = 0.2. 

Open with DEXTER 
Fig. 5 Same as Fig. 4 for the different initial disc masses. 

Open with DEXTER 
Hereafter, we describe the orbital characteristics of the planetary systems formed by our mechanism. Semimajor axes, eccentricities, and inclinations are discussed in Sects. 4.1–4.3, respectively. Section 4.5 analyzes the mutual Hill separation of the orbits, while threebody resonances are studied in Sect. 4.6. Unless otherwise stated, our analysis relates to the planetary systems of the DM model with a > 0.1 AU (because tidal/relativistic effects are not included in our work) and mass in the interval [0.65,10] M_{Jup}, as they emerge from the disc phase (at 1.4 × 10^{6} yr), namely 13 036 planets in 5644 systems. Comparisons are made with the exoplanets of the observational data that suffer from the same limitations in semimajor axis and mass.
4.1. Semimajor axes
We first investigate the effect of the gas disc on the final distribution of semimajor axis. Figure 4 shows the normalized semimajor axis distribution (in logarithmic scale), for all the systems of the DM model as they emerge from the gas phase (green dashed line). For completeness, the CM model is also added to the plot (blue dashed line) in order to evaluate the efficiency of the orbital migration. Indeed, the efficiency of the migration mechanism depends on the amount of surrounding gas in the vicinity of the planet, and the planets end up with smaller semimajor axis in the CM modelization, as expected. The black solid line represents the observed giant planet population (we have excluded the planets with a < 0.1 AU and mass out of the interval [0.65,10] M_{Jup}). Unlike the CM model, the DM semimajor axis distribution and the one of the observations follow a similar trend.
Fig. 6 Cumulative eccentricity distributions at the dispersal of the disc. The black lines correspond to the observations. Top panel: the dependence of the distribution of all systems on the number of planets. Bottom panel: the dependence of the distribution of multiplanetary systems on the initial disc mass. 

Open with DEXTER 
We further analyze the DM model semimajor axis distribution, investigating the influence of the initial disc mass in Fig. 5. The different panels present the normalized semimajor axes distribution (in logarithmic scale) for initial disc mass of 4, 8, 16, and 32 M_{Jup}, respectively. For more massive discs, the migration mechanism is more efficient and the planets evolve closer to the parent star. Let us note that the best agreement between the simulated and observed distributions is obtained for the 16 M_{Jup} disc (bottom left panel of Fig. 5).
4.2. Eccentricities
We find good agreement between the eccentricities of the planetary systems obtained by our simulations and the detected ones, as shown in Fig. 6, where both cumulative eccentricity distributions are displayed (top panel). We see that the eccentricities observed in our simulations are well diversified and present the same general pattern as the observed eccentricities. In particular, the similarity of the curves for eccentricities smaller than 0.35 is impressive. However, our model underproduces highly eccentric orbits: only ~9.7% of the planets have an eccentricity higher than 0.5 in our simulations. We highlight that the modelization of the damping formulae is not accurate for e> 2/3 (Sect. 2). Also, the orbital elements of our systems are considered immediately after the dispersal of the disc (1.4 × 10^{6} yr) and we show in Sect. 5 that orbital adjustments due to planetplanet interactions can still come into play afterwards.
Furthermore, it is well known that part of the high eccentricities reported in the observations could originate from additional mechanisms not considered here, such as the gravitational interactions with a distant companion (binary star or massive giant planet). These highly eccentric planets are mostly in single planet systems, at least in the observational data for a distance of 30 AU from the star. For this reason, in the bottom panel of Fig. 6, we focus on the eccentricities of (final) multiplanetary systems only. We first observe that removing the singleplanet systems reduces the disagreement at high eccentricities, in particular from 0.35 up to 0.55, higher eccentricities not being correctly modelized by our study.
The cumulative eccentricity distributions for the four initial disc masses considered in this study are also shown in the bottom panel of Fig. 6 (dashed lines). The lower the mass of the disc, the lower the eccentricities of the multiplanetary configurations observed in the simulations. The initial disc mass of 16M_{Jup} gives the best approximation of the observational distribution. To put this in perspective, let us add that the mean value of the total mass of the planetary systems at the beginning of the simulations is 10.5M_{Jup}. Since our approximation well matches the observations for small and medium eccentricities, this would imply that the observed extrasolar systems are consistent with an initial disc mass comparable to the total mass of the planets initially formed in the disc.
The correlation between semimajor axis and eccentricity is examined in the left panel of Fig. 7. In this analysis, we consider the 16 M_{Jup} disc simulations, since the previous figures suggest that this initial disc mass presents the best curve fittings for both the semimajor axis and the eccentricity independently. The left panel of Fig. 7 shows the matching between our simulations and the observations in the semimajor axis vs. eccentricity graph. The results are in agreement with the study of Matsumura et al. (2010), combining Nbody dynamics with hydrodynamical disc evolution.
Finally, the mean eccentricity for each system that ends up in a multiplanetary configuration in our simulations is shown in the right panel of Fig. 7 as a function of the total mass of the system. For the systems consisting of two planets at the dispersal of the disc, the higher the total mass of the system, the higher the mean eccentricity. Indeed, the eccentricities are excited by meanmotion resonance capture, planetplanet scattering and interactions with the disc. These mechanisms are more efficient for massive planets. In particular, concerning the last one, the disc, instead of damping the planetary eccentricities, can induce eccentricity excitation (Papaloizou et al. 2001; D’Angelo et al. 2006) for highmass planets (M_{p} ≳ 5.5 M_{Jup}). This feature is included in the damping formulae derived in Paper I. For the threeplanet systems, the mean eccentricities are smaller, since planetplanet scattering did not take place and only the other two mechanisms operated.
Fig. 7 Left panel: semimajor axis vs. eccentricity for our simulations with M_{Disc} = 16 M_{Jup} (green dots) and the observations (black dots). Right panel: total mass vs. mean eccentricity for each multiplanetary configuration of our 16 M_{Jup} disc simulations. Black dots correspond to threeplanet systems and green dots to twoplanet systems. 

Open with DEXTER 
We do not discuss here the eccentricities of the CM model as no significant difference was observed. The same holds true for the inclinations.
Fig. 8 Normalized inclination distribution for all the planets with i> 1° at the dispersal of the disc. Top panel: the dependence of the distribution on the initial disc mass. Bottom panel: the dependence of the distribution on the number of planets in the (final) systems. The bin size is Δlog (i) = 0.2. 

Open with DEXTER 
Fig. 9 Mutual inclination of the twoplanet systems (in logarithmic scale) as a function of the total mass of the planets in the system for the four different disc masses. 

Open with DEXTER 
4.3. Inclinations
Fig. 10 Initial and final Hill neighbor separation normalized distributions (D_{Hill} = (Δa_{i,j}) /R_{Hill,ai,aj}) of the systems of our simulations. The left column displays the twoplanet systems and the right column the inner and outer pairs of the threeplanet systems. Each panel represents a different initial disc mass. The bin size is ΔD_{Hill} = 1.2. 

Open with DEXTER 
In Fig. 8 we present the normalized inclination distribution (black solid line), as a function of the initial disc mass (top panel), and the number of planets in the (final) systems (bottom panel). Only the planets with i > 1° are represented. A clear outcome from the top panel is that the smaller the initial disc mass, the larger the fraction of planets with higher inclinations in the timescale of disc’s lifetime. For a 4M_{Jup} disc, ~ 24% of the planets have i > 1°, while this percentage drops dramatically to approximately 7% for a 32 M_{Jup} disc. It is simply explained by the fact that a more massive disc exerts a stronger damping on the inclination of the planet.
The normalized inclination distribution (when i > 1°) also strongly depends on the final number of bodies in the system after the gas phase, as shown in the bottom panel of Fig. 8. As expected, higher inclinations are typically found for twoplanet systems, since these systems have undergone a scattering event during the gas phase. Moreover, regarding the mutual inclination of twoplanet systems, there is no correlation between the total mass of the system and the mutual inclination of the planetary orbits. In Fig. 9 we show that, for the four disc masses considered here, mutual inclinations can be pumped regardless of the total mass of the systems.
4.4. Hot Jupiters
Concerning hot Jupiters, ~ 2% of the planets are found with semimajor axes in the range [0.02,0.2] AU after the dispersal of the disc. We remind the reader that our model is not accurate for planets with a < 0.1 AU, since we do not consider tidal/relativistic effects. However, our simulations seem to indicate that planetplanet interactions during migration in the protoplanetary disc can produce hot Jupiters on eccentric orbits while the migration of single planets would leave them circular. In our results, ~ 24% of these closein planets have eccentricities higher than 0.2. This suggests that planetary migration, including planetplanet interactions, is a valid mechanism to produce both circular and eccentric hot Jupiters.
Moreover, it appears difficult to excite their inclinations larger than 10° with respect to the midplane of the disc (only ~ 12% of them reach inclination higher than 1° in our simulations). In this frame, the observed misalignment of approximately half of the hot Jupiters with respect to the stellar equator could be the result of a misalignment between the old disc midplane and the present stellar equator (see for instance Crida & Batygin 2014).
4.5. Mutual Hill separation
The mutual Hill radius of two planets is defined, by Gladman (1993), as: (6)and we express the Hill neighbor separation of two planetary orbits as D_{Hill} = (a_{j}−a_{i}) /R_{Hill,ai,aj}. In Fig. 10, we show the initial and final Hill neighbor separations of the systems of our simulations. Systems ending up with twoplanets are shown in the left column, while the right column shows both inner and outer pairs of the threeplanet systems with different colored dashed lines. Each panel corresponds to a different initial value of the mass of the disc.
It is clear that the threeplanet configurations are more compact compared with systems that suffered scattering events. Also, the distribution of the inner pair is slightly wider than the one of the outer pair since only the outer planet is migrating inwards due to the interaction with the disc in our simulations. Concerning the Hill neighbor separation distribution of the twoplanet systems, we see that the higher the initial disc mass, the larger the mean separation.
Threebody meanmotion resonances at the dispersal of the disc for both disc modelizations.
4.6. Threebody resonances
As previously mentioned, 40% of the systems (in our DM model simulations) end up in a threeplanet configuration (Table 3). Depending on the initial separation of the planets, the systems are generally captured in two or threebody meanmotion resonances during the migration phase. These resonances can either survive until the end of the disc phase or be disrupted in an instability phase. For the systems consisting of three planets at the end of the disc phase, ~10% of them are not in meanmotion resonance, ~ 25% are in a twobody meanmotion resonance and ~ 65% are in a threebody meanmotion resonance. This high percentage (~ 90%) of resonant systems is also observed in Matsumura et al. (2010). Table 4 shows that half of the threebody resonant systems are in a 1:2:4 Laplace configuration. The n_{3}:n_{2}:n_{1} = 1:3:6 and n_{3}:n_{2}:n_{1} = 2:5:10 resonances are the second most common configurations. These percentages being the same in both disc modelizations shows that unlike the semimajor axis distribution, the establishment of the resonant threebody configurations is not affected by the decrease of the gas disc.
Dynamical history of highly mutually inclined systems. The first column shows all the possible outcomes.
Fig. 11 Examples of several scenarios producing mutual inclination excitation: a) planetplanet scattering followed by meanmotion resonance; b) orbital instability; c) orbital instability followed by threebody resonance; and d) threebody resonance followed by orbital instability. 

Open with DEXTER 
4.7. Dynamical history of highly mutually inclined systems
Here we focus on systems formed on highly mutually inclined orbits (≥ 10°), which represent ~3% of our simulations at the dispersal of the disc. Mechanisms producing inclination increase have been identified in Sect. 3, namely inclinationtype resonance and planetplanet scattering. We now aim at understanding the frequency of each mechanism, keeping in mind that the occurrence of both mechanims is possible during the longterm evolution in the disc phase.
By carefully studying the dynamical evolution of the highly mutually inclined systems formed in our simulations, we have found seven possible dynamical histories and the frequency of each one is given in Table 5 (third column). Concerning the systems finally composed of two planets at the dispersal of the disc (1.4×10^{6} yr), three scenarios could have happened: either planetplanet scattering induces the ejection/collision of a planet and excites the mutual inclinations of the remaining planets (similarly to Fig. 2), the inclinations produced by an ejection/collision of a body due to planetplanet scattering are rapidly damped by the disc and a meanmotion resonance capture of the two bodies induces the increase of the inclinations (Fig. 11, panel a), or a threebody inclinationtype resonance excites the inclinations and consequently leads to planetplanet scattering and the ejection of a planet (similarly to Fig. 12). Mutual inclination is always observed in threebody systems following four different histories: either inclination increase is produced by orbital instability of the orbits (without planet ejection) (Fig. 11, panel b), the orbital elements excited by orbital instability are damped until a threebody resonance capture and a resonant excitation of the inclinations (Fig. 11, panel c), the system evolves smoothly in a threebody inclinationtype resonance (similarly to Fig. 3), or the phase of threebody inclinationtype resonance is followed by a destabilization of the orbits (Fig. 11, panel d).
Considering the last mechanism as the source of the inclination excitation, we see that half of the highly mutually inclined systems of our simulations result from two or threebody meanmotion resonance captures, the other half being produced by orbital instability and/or planetplanet scattering. This emphasizes the importance of meanmotion resonance captures during the disc phase, on the final 3D configurations of planetary systems.
5. Longterm evolution
Section 4 focuses on the distributions of the orbital elements of the planetary systems considered immediately after the dispersal of the disc. As previously discussed, orbital adjustments due to planetplanet interactions can occur on a longer timescale. An example is given in Fig. 12, showing the destabilization of a system in a 1:2:6 resonance at the end of the disc phase. At ~ 38 Myr, the middle planet is ejected from the system and the two surviving bodies are left in wellseparated and stable orbits.
In this section, we aim to investigate whether or not the longterm evolution of planetary systems produces significant changes on the final distribution of the orbital elements discussed hereabove. To study the longterm evolution of planetary systems, we ran two additional sets of simulations for 100 Myr: 400 systems for an initial disc mass of 8 M_{Jup} and 800 systems for 16 M_{Jup}, both with an exponential decay of the mass disc (DM model, with the same dispersal time of ~1 Myr). As expected, many systems of our longterm simulations are destabilized after millions of years, and the percentages given by Table 3 on the number of planets in the final configurations have changed significantly. There is a clear tendency towards systems with fewer planets on a long timescale.
Fig. 12 Destabilization of a threebody resonance on a long timescale. While the system is locked in a 1:2:6 resonance at the dispersal of the gas disc, a planetplanet scattering event finally takes place and the middle planet is ejected from the system at ~ 38 Myr. The planetary masses are m_{1} = 0.87, m_{2} = 1.39 and m_{3} = 8.43 M_{Jup}. The initial mass of the disc is 16 M_{Jup}. 

Open with DEXTER 
Fig. 13 Left panel: normalized eccentricity distributions of the 16 M_{Jup} disc simulations for two integration timescales; 1.4 × 10^{6} and 1 × 10^{8} yr. The bin size is Δe = 0.05. Right panel: normalized semimajor axis distributions of the 16 M_{Jup} disc simulations, for the same two integration timescales. The bin size is Δlog (a) = 0.2. 

Open with DEXTER 
Fig. 14 Normalized inclination distributions for two integration timescales (1.4 × 10^{6} and 1 × 10^{8} yr) and two initial disc masses (8 and 16 M_{Jup}). Only the inclinations higher than 1° are represented. The bin size is Δlog (i) = 0.2. 

Open with DEXTER 
Figure 13 shows no significant change on the semimajor axis and eccentricity distributions on the 100 Myr timescale. However, the inclinations have considerably increased on a longer timescale, as appears clearly in Fig. 14. Just after the disc phase, ~ 10% of the planets have inclinations higher than 1° and this percentage has almost doubled at 100 Myr (~ 25% for the 8 M_{Jup} and ~ 17% for the 16 M_{Jup}).
Concerning the mutual inclinations, ~ 3% of the threeplanet systems are highly mutually inclined (I_{mut}> 10°) at 1.4 × 10^{6} yr and this percentage remains approximately the same for the longterm simulations. The situation is quite different for the twoplanet systems. There are approximately 2% of highly mutually inclined systems at the dispersal of the gas disc and ~ 7% at 100 Myr. As a result, ~ 5% of the multiple systems of our simulations have high mutual inclinations (I_{mut}> 10°).
The destabilization of highly mutually inclined systems is also common when considering the longterm evolution on 100 Myr. The same observation has recently been pointed out by Barnes et al. (2015), concerning planetary systems in meanmotion resonance with mutual inclinations. Referring to Table 5 (fourth column), we see that the percentage of highly mutually inclined systems still evolving in resonance drops to 30% (instead of 50% at the dispersal of the disc).
6. Conclusions
In this study we followed the orbital evolution of three giant planets in the late stage of the gas disc. Our scenario for the formation of planetary systems combines TypeII migration, with the consistent eccentricity and inclination damping of Paper I (Bitsch et al. 2013), and planetplanet scattering. The results shown in this work are naturally attached to the disc parameters of the hydrodynamical simulations of Paper I from which the damping formulae are issued. Our parametric experiments consisted of 11 000 numerical simulations, considering a variety of initial configurations, planet mass ratios, and disc masses. Moreover, two modelizations of the gas disc were taken into account: the constantmass model (no gas dissipation during 0.8 Myr) and the decreasingmass model (gas exponential decay with an efolding timescale of 1 Myr). The first case leads to more merging, more migration, and less ejections of planets.
We focused on the impact of eccentricity and inclination damping on the final configuration of planetary systems. We have shown that the eccentricities are already welldiversified at the dispersal of the disc, despite the strong eccentricity damping exerted by the gas disc, and subsequent inclination increase is possible. Concerning the inclinations, in contrast with previous works that did not include inclination damping, we found that most of the planets end up in the midplane of the disc (i.e., in quasicoplanar orbits with i< 1°), showing the efficiency of the inclination damping. One should keep in mind that planets formed in the disc midplane could appear inclined with respect to the stellar equatorial plane if these two planes differ. Needless to say, the higher the initial disc mass, the smaller the inclinations of the planets at the dispersal of the disc. Nevertheless, in multiple systems, inclinationtype resonance and planetplanet scattering events during/after the gas phase can produce inclination excitation of the inclinations. Approximately 5% of highly mutually inclined systems (I_{mut}> 10°) have been formed in our scenario. In future observations, this percentage could help to discriminate between the formation scenarios.
The dynamical mechanisms producing inclination increase were identified, namely inclinationtype resonance, orbital instability, and/or planetplanet scattering. We have shown that resonance captures play an important role in the formation of mutually inclined systems. Indeed, half of these systems originate from two or threebody meanmotion resonance captures. The longterm evolution of the systems was been investigated, showing that destabilization of the resonant systems is common. However, 30% of the systems still evolve in resonance after 100 Myr.
As a byproduct of our study, we found a very good agreement between our simulations and the observed population of extrasolar systems, in particular for the semimajor axis and eccentricity distributions. Although a full exploration of the parameter space and a real population synthesis study are far beyond the scope of this paper, this agreement suggests strongly that planetplanet interactions during the migration inside the protoplanetary disc could account for most of the eccentricity excitation observed among exoplanets.
When all the planets migrate and do not share a common gap, the rate for TypeII regime given by Eq. (1) leads to divergent migration (instead of convergent migration like in the present work) and the same phenomenons of resonance capture and eccentricity/inclination excitation are only temporarily observed (see e.g. Libert & Tsiganis 2011b).
Acknowledgments
This work was supported by the Fonds de la Recherche ScientifiqueFNRS under Grant No. T.0029.13 (“ExtraOrDynHa” research project). Computational resources were provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.FNRS) under Grant No. 2.5020.11. We thank the referee for his constructive suggestions that have improved the manuscript.
References
 Adams, F. C., & Laughlin, G. 2003, Icarus, 163, 290 [NASA ADS] [CrossRef] [Google Scholar]
 Albrecht, S., Winn, J. N., Johnson, J. A., et al. 2012, ApJ, 757, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, R., Deitrick, R., Greenberg, R., Quinn, T. R., & Raymond, S. N. 2015, ApJ, 801, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., Crida, A., Paardekooper, S.J., et al. 2014, Protostars and Planets VI, 667 [Google Scholar]
 Bitsch, B., & Kley, W. 2010, A&A, 523, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., & Kley, W. 2011, A&A, 530, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Crida, A., Libert, A.S., & Lega, E. 2013, A&A, 555, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chatterjee, S., Ford, E., Matsumura, S., & Rasio, F. 2008, ApJ, 686 [Google Scholar]
 Cresswell, P., Dirksen, G., Kley, W., & Nelson, R. P. 2007, A&A, 473, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crida, A., & Batygin, K. 2014, A&A, 567, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crida, A., & Morbidelli, A. 2007, MNRAS, 377, 1324 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., Sándor, Z., & Kley, W. 2008, A&A, 483, 325 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 D’Angelo, G., Lubow, S. H., & Bate, M. R. 2006, ApJ, 652, 1698 [NASA ADS] [CrossRef] [Google Scholar]
 Deitrick, R., Barnes, R., McArthur, B., et al. 2015, ApJ, 798, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Duffell, P. C., & Chiang, E. 2015, ApJ, 812, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Duffell, P. C., Haiman, Z., MacFadyen, A. I., D’Orazio, D. J., & Farris, B. D. 2014, ApJ, 792, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067 [NASA ADS] [CrossRef] [Google Scholar]
 Dürmann, C., & Kley, W. 2015, A&A, 574, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ford, E., & Rasio, F. 2008, ApJ, 686, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Gladman, B. 1993, Icarus, 106, 247 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Sari, R. 2003, ApJ, 585, 1024 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hasegawa, Y., & Ida, S. 2013, ApJ, 774, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Ivanov, P., Papaloizou, J. C. B., & Polnarev, A. 1999, MNRAS, 307, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Jurić, M., & Tremaine, S. 2008, ApJ, 686 [Google Scholar]
 Kley, W. 2000, MNRAS, 314, L47 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W., & Dirksen, G. 2006, A&A, 447, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596 [NASA ADS] [CrossRef] [Google Scholar]
 Lega, E., Morbidelli, A., & Nesvorný, D. 2013, MNRAS, 431, 3494 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., & Duncan, M. J. 2000, AJ, 120, 2117 [NASA ADS] [CrossRef] [Google Scholar]
 Libert, A.S., & Tsiganis, K. 2009, MNRAS, 400, 1373 [NASA ADS] [CrossRef] [Google Scholar]
 Libert, A.S., & Tsiganis, K. 2011a, MNRAS, 412, 2353 [NASA ADS] [CrossRef] [Google Scholar]
 Libert, A.S., & Tsiganis, K. 2011b, Celest. Mech. Dyn. Astron., 111, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D. N. C., & Ida, S. 1997, ApJ, 477, 781 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. C. B. 1986, ApJ, 309, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Mamajek, E. E. 2009, in AIP Conf. Ser. 1158, eds. T. Usuda, M. Tamura, & M. Ishii, 3 [Google Scholar]
 Marzari, F., & Weidenschilling, S. J. 2002, Icarus, 156, 570 [NASA ADS] [CrossRef] [Google Scholar]
 Marzari, F., Baruteau, C., & Scholl, H. 2010, A&A, 514, 4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matsumoto, Y., Nagasawa, M., & Ida, S. 2012, Icarus, 221, 624 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumura, S., Thommes, E., Chatterjee, S., & Rasio, F. 2010, ApJ, 714, 194 [NASA ADS] [CrossRef] [Google Scholar]
 Moeckel, N., & Armitage, P. 2012, MNRAS, 419, 366 [NASA ADS] [CrossRef] [Google Scholar]
 Moorhead, A. V., & Adams, F. C. 2005, Icarus, 178, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Ogilvie, G. I., & Lubow, S. H. 2003, ApJ, 587, 398 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., & Larwood, J. 2000, MNRAS, 315, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Petrovich, C., Tremaine, S., & Rafikov, R. 2014, ApJ, 786, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Rasio, F. A., & Ford, E. B. 1996, Science, 274, 954 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388 [NASA ADS] [CrossRef] [Google Scholar]
 Terquem, C., & Ajmia, A. 2010, MNRAS, 404, 409 [NASA ADS] [Google Scholar]
 Teyssandier, J., & Ogilvie, G. I. 2016, MNRAS, 458, 3221 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssandier, J., & Terquem, C. 2014, MNRAS, 443, 568 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssandier, J., Terquem, C., & Papaloizou, J. C. B. 2013, MNRAS, 428, 658 [NASA ADS] [CrossRef] [Google Scholar]
 Thommes, E., & Lissauer, J. 2003, ApJ, 597, 566 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J., & Marzari, F. 1996, Nature, 384, 619 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409 [NASA ADS] [CrossRef] [Google Scholar]
 XiangGruess, M., & Papaloizou, J. C. B. 2013, MNRAS, 431, 1320 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Initial disc mass, integration time, and number of simulations for each disc modelization.
Number of planets in the systems of our simulations at the end of the integration time, for both disc models.
Threebody meanmotion resonances at the dispersal of the disc for both disc modelizations.
Dynamical history of highly mutually inclined systems. The first column shows all the possible outcomes.
All Figures
Fig. 1 Normalized initial and final mass distributions of our simulations compared with the massdistribution (in the range [0.65,10] M_{Jup}) of the observed giant planets with a > 0.1. 

Open with DEXTER  
In the text 
Fig. 2 Inclination excitation due to planetplanet scattering (DM model). After the hyperbolic ejection of the inner planet, the system consists of two planets with a large orbital separation, large eccentricities and slightly inclined orbits. The planetary masses are m_{1} = 1.98, m_{2} = 8.35, m_{3} = 3.8 M_{Jup}. 

Open with DEXTER  
In the text 
Fig. 3 Inclination excitation due to inclinationtype resonance in a Laplace 1:2:4 meanmotion resonance (DM model). The planetary masses are m_{1} = 2.66, m_{2} = 1.44 and m_{3} = 0.73 M_{Jup}. 

Open with DEXTER  
In the text 
Fig. 4 Normalized semimajor axis distributions in both disc modelizations for all the initial disc masses (blue and green lines), and in the observed giant planet population (with a> 0.1 AU and M_{p} ∈ [0.65,10] M_{Jup}, black line). The bin size is Δlog (a) = 0.2. 

Open with DEXTER  
In the text 
Fig. 5 Same as Fig. 4 for the different initial disc masses. 

Open with DEXTER  
In the text 
Fig. 6 Cumulative eccentricity distributions at the dispersal of the disc. The black lines correspond to the observations. Top panel: the dependence of the distribution of all systems on the number of planets. Bottom panel: the dependence of the distribution of multiplanetary systems on the initial disc mass. 

Open with DEXTER  
In the text 
Fig. 7 Left panel: semimajor axis vs. eccentricity for our simulations with M_{Disc} = 16 M_{Jup} (green dots) and the observations (black dots). Right panel: total mass vs. mean eccentricity for each multiplanetary configuration of our 16 M_{Jup} disc simulations. Black dots correspond to threeplanet systems and green dots to twoplanet systems. 

Open with DEXTER  
In the text 
Fig. 8 Normalized inclination distribution for all the planets with i> 1° at the dispersal of the disc. Top panel: the dependence of the distribution on the initial disc mass. Bottom panel: the dependence of the distribution on the number of planets in the (final) systems. The bin size is Δlog (i) = 0.2. 

Open with DEXTER  
In the text 
Fig. 9 Mutual inclination of the twoplanet systems (in logarithmic scale) as a function of the total mass of the planets in the system for the four different disc masses. 

Open with DEXTER  
In the text 
Fig. 10 Initial and final Hill neighbor separation normalized distributions (D_{Hill} = (Δa_{i,j}) /R_{Hill,ai,aj}) of the systems of our simulations. The left column displays the twoplanet systems and the right column the inner and outer pairs of the threeplanet systems. Each panel represents a different initial disc mass. The bin size is ΔD_{Hill} = 1.2. 

Open with DEXTER  
In the text 
Fig. 11 Examples of several scenarios producing mutual inclination excitation: a) planetplanet scattering followed by meanmotion resonance; b) orbital instability; c) orbital instability followed by threebody resonance; and d) threebody resonance followed by orbital instability. 

Open with DEXTER  
In the text 
Fig. 12 Destabilization of a threebody resonance on a long timescale. While the system is locked in a 1:2:6 resonance at the dispersal of the gas disc, a planetplanet scattering event finally takes place and the middle planet is ejected from the system at ~ 38 Myr. The planetary masses are m_{1} = 0.87, m_{2} = 1.39 and m_{3} = 8.43 M_{Jup}. The initial mass of the disc is 16 M_{Jup}. 

Open with DEXTER  
In the text 
Fig. 13 Left panel: normalized eccentricity distributions of the 16 M_{Jup} disc simulations for two integration timescales; 1.4 × 10^{6} and 1 × 10^{8} yr. The bin size is Δe = 0.05. Right panel: normalized semimajor axis distributions of the 16 M_{Jup} disc simulations, for the same two integration timescales. The bin size is Δlog (a) = 0.2. 

Open with DEXTER  
In the text 
Fig. 14 Normalized inclination distributions for two integration timescales (1.4 × 10^{6} and 1 × 10^{8} yr) and two initial disc masses (8 and 16 M_{Jup}). Only the inclinations higher than 1° are represented. The bin size is Δlog (i) = 0.2. 

Open with DEXTER  
In the text 