© 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 spin-orbit 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 semi-major axis and keeps them on coplanar and near-circular 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 Type-II 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; Xiang-Gruess & Papaloizou 2013; Bitsch et al. 2013).

Several authors, using analytical prescriptions, have investigated the impact of the planet-disc interactions on the orbital arrangement of planetary systems through three-dimensional (3D) N-body simulations. Studies in the context of Type-II migration have been performed for two-planet systems (Thommes & Lissauer 2003; Libert & Tsiganis 2009; Teyssandier & Terquem 2014) and three-planet systems (Libert & Tsiganis 2011b). A notable outcome of these works is that, during convergent migration, the system can enter an inclination-type resonance that pumps the mutual inclination of the planetary orbits.

After the dispersal of the disc, dynamical instabilities among the planets can result in planet-planet 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 planet-planet 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 (Type-II migration and scattering). Parametric explorations via N-body 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 N-body 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 first-order approximation for the eccentricity damping timescale (ė/e = −K | ȧ/a |). The considered K-factor, 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 planet-disc interactions. Several works have investigated the evolution of the eccentricities of giant planets in discs using two-dimensional (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 N-body simulations, as a function of eccentricity, inclination, planetary mass, and disc mass.

In this paper, we aim to improve the previous N-body studies combining planet-disc interactions and planet-planet interactions by adopting the damping formulae of Paper I. In other words, we combine the speed and efficiency of an N-body 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 planet-disc interactions on the final configurations of planetary systems and their inclination distribution. Through our parametric analysis, we also discuss the most common three-body 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 set-up 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 long-term evolution on the architecture of the systems. Finally, our conclusions are given in Sect. 6.

2. Set-up of the N-body simulations

2.1. Physical modeling

In this work, we analyze the evolution of three-planet systems embedded in a protoplanetary disc and evolving around a solar-mass 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 time-step 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 close-in orbits, we set the time-step of our simulations at dt = 0.001 yr.

The migration of the planets in Type-II migration, due to angular momentum exchange with the disc, is on a similar timescale as the viscous accretion time, ${\mathit{t}}_{\mathit{\nu }}\mathrm{~}{\mathit{r}}_{\mathrm{pl}}^{\mathrm{2}}\mathit{/}\mathit{\nu }$, 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 Type-II regime consists of two different cases, the disc-dominated case and the planet-dominated case: ${\mathit{\tau }}_{\mathrm{II}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\alpha }}^{-1}{\mathit{h}}^{-2}{\mathrm{\Omega }}_{\mathrm{pl}}^{-1}\mathrm{\times }\begin{array}{c}\max \end{array}\left(\mathrm{1}\mathit{,}\frac{{\mathit{M}}_{\mathrm{p}}}{\mathrm{(}\mathrm{4}\mathit{\pi }\mathit{/}\mathrm{3}\mathrm{)}\mathrm{\Sigma }\mathrm{\left(}{\mathit{r}}_{\mathrm{pl}}\mathrm{\right)}{\mathit{r}}_{\mathrm{pl}}^{\mathrm{2}}}\right)\mathit{,}$(1)where α is the Shakura-Sunyaev viscosity parameter (Shakura & Sunyaev 1973), h the disc aspect ratio, Ω the orbital frequency of the planet and $\mathrm{(}\mathrm{4}\mathit{\pi }\mathit{/}\mathrm{3}\mathrm{)}\mathrm{\Sigma }\mathrm{\left(}{\mathit{r}}_{\mathrm{pl}}\mathrm{\right)}{\mathit{r}}_{\mathrm{pl}}^{\mathrm{2}}$ is the local mass of the disc. The rate of Type-II 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 Type-II regime (which would require hydrodynamical simulations), but planet-planet 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 Type-II 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 gas-free inner cavity, following Matsumoto et al. (2012), by using a hyperbolic tangent function, $\tanh \frac{\mathrm{(}\mathit{r}\mathrm{-}{\mathit{R}}_{\mathrm{in}}\mathrm{)}}{\mathrm{\Delta }\mathit{r}}$, 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¹. 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): ${{a}}_{\mathrm{disc}}\mathrm{=}\mathrm{-}\frac{{v}}{{\mathit{\tau }}_{\mathrm{II}}}\mathrm{-}\mathrm{2}\frac{\mathrm{(}{v}\mathrm{·}{r}\mathrm{)}{r}}{{\mathit{r}}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{ecc}}}\mathrm{-}\mathrm{2}\frac{\mathrm{(}{v}\mathrm{·}{k}\mathrm{)}{k}}{{\mathit{\tau }}_{\mathrm{inc}}}\mathit{,}$(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 constant-mass 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 decreasing-mass 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^-9M_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 planet-planet 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 mass-distribution (in the range [0.65,10] MJup) of the observed giant planets with a > 0.1.

In our simulations, gas giant planets are initially on coplanar and quasi-circular orbits (e ∈ [0.001,0.01] and i ∈ [0.01°,0.1°]). The initial semi-major axis of the inner body is fixed to 3 AU and the middle and outer ones follow uniform distributions in the intervals a₂ ∈ [4.3,6.5] AU and a₃ ∈ [8,16] AU, respectively. These initial semi-major axis distributions provide a wide range of possible resonance captures. We choose randomly initial planetary masses from a log-uniform distribution that approximately fits with the observational data² 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⁵ 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 low-inclination 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 single-planet systems, the dynamics tend to damp the planet towards the midplane of the disc, in a circular orbit in the case of a low-mass planet and in an orbit whose eccentricity increases over time due to the interactions with the disc in the case of a high-mass planet and a sufficiently massive gas disc. The increase of the planetary eccentricity for massive planets may appear surprising, but it is a well-established 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 $\begin{array}{ccc}\frac{\mathrm{d}\mathit{e}}{\mathrm{d}\mathit{t}}\mathrm{\left(}{\mathit{M}}_{\mathrm{P}}\mathit{,}{\mathit{e}}_{\mathrm{P}}\mathit{,}{\mathit{i}}_{\mathrm{P}}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{{\mathit{M}}_{\mathrm{disc}}}{\mathrm{0.01}{\mathit{M}}_{\mathit{\star }}}{\mathrm{(}\mathit{a}\mathrm{(}{\mathit{i}}_{\mathrm{P}}\mathrm{+}{\mathit{i}}_{\mathrm{D}}{\mathrm{)}}^{\mathrm{-}\mathrm{2}\mathit{b}}\mathrm{+}\mathit{c}{{\mathit{i}}_{\mathrm{P}}^{\mathrm{-}\mathrm{2}\mathit{d}}}^{\mathrm{)}}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\\ & & \end{array}$(3)where i_D = M_p/ 3 degrees, M_⋆ is the mass of the star in solar mass, and with coefficients

$\begin{array}{ccc}& & {\mathit{a}}_{\mathit{e}}\mathrm{\left(}{\mathit{M}}_{\mathrm{P}}\mathit{,}{\mathit{e}}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}\mathrm{80}\hspace{0.17em}{\mathit{e}}_{\mathrm{P}}^{-2}\hspace{0.17em}\exp \mathrm{(}\mathrm{-}{\mathit{e}}_{\mathrm{P}}^{\mathrm{2}}{\mathit{M}}_{\mathrm{p}}\mathit{/}{\mathrm{0.26}}^{\mathrm{)}}{\mathrm{15}}^{{\mathit{M}}_{\mathrm{p}}}\hspace{0.17em}\mathrm{(}\mathrm{20}\mathrm{+}\mathrm{11}{\mathit{M}}_{\mathrm{p}}\mathrm{-}{{\mathit{M}}_{\mathrm{p}}^{\mathrm{2}}}^{\mathrm{)}}\\ & & {\mathit{b}}_{\mathit{e}}\mathrm{\left(}{\mathit{M}}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}\mathrm{0.3}{\mathit{M}}_{\mathrm{p}}\\ & & {\mathit{c}}_{\mathit{e}}\mathrm{\left(}{\mathit{M}}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}\mathrm{450}\mathrm{+}{\mathrm{2}}^{{\mathit{M}}_{\mathrm{p}}}\\ & & {\mathit{d}}_{\mathit{e}}\mathrm{\left(}{\mathit{M}}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{1.4}\mathrm{+}\sqrt{{\mathit{M}}_{\mathrm{p}}}\mathit{/}\mathrm{6.}\end{array}$The damping function for inclination is given, in degrees per orbit, by $\begin{array}{ccc}\frac{\mathrm{d}\mathit{i}}{\mathrm{d}\mathit{t}}\mathrm{\left(}{\mathit{M}}_{\mathrm{P}}\mathit{,}{\mathit{e}}_{\mathrm{P}}\mathit{,}{\mathit{i}}_{\mathrm{P}}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{{\mathit{M}}_{\mathrm{disc}}}{\mathrm{0.01}{\mathit{M}}_{\mathit{\star }}}[{\mathit{a}}_{\mathit{i}}\hspace{0.17em}{\left(\frac{{\mathit{i}}_{\mathrm{P}}}{{\mathrm{1}}^{\mathrm{◦}}}\right)}^{\mathrm{-}\mathrm{2}{\mathit{b}}_{\mathit{i}}}\exp \mathrm{(}\mathrm{-}\mathrm{(}{\mathit{i}}_{\mathrm{P}}\mathit{/}{\mathit{g}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{)}\\ & & \end{array}$(4)with coefficients

$\begin{array}{ccc}& & {\mathit{a}}_{\mathit{i}}\mathrm{\left(}{\mathit{M}}_{\mathrm{P}}\mathit{,}{\mathit{e}}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}\mathrm{1.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{3}{\mathit{e}}_{\mathrm{P}}\mathrm{)}{{\mathit{M}}_{\mathrm{p}}}^{\mathrm{3}}\\ & & {\mathit{b}}_{\mathit{i}}\mathrm{\left(}{\mathit{M}}_{\mathrm{P}}\mathit{,}{\mathit{e}}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}{\mathit{M}}_{\mathrm{p}}{\mathit{e}}_{\mathrm{P}}^{\mathrm{2}}\mathit{/}\mathrm{10}\\ & & {\mathit{c}}_{\mathit{i}}\mathrm{\left(}{\mathit{M}}_{\mathrm{P}}\mathit{,}{\mathit{e}}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}\mathrm{1.2}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}\mathit{/}[\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{3}{\mathit{e}}_{\mathrm{P}}\mathrm{)}\mathrm{(}\mathrm{5}\mathrm{+}{\mathit{e}}_{\mathrm{P}}^{\mathrm{2}}\mathrm{(}{\mathit{M}}_{\mathrm{p}}\mathrm{+}\mathrm{2}{\mathrm{)}}^{\mathrm{3}}\mathrm{)}]\\ & & {\mathit{d}}_{\mathit{i}}\mathrm{\left(}{\mathit{e}}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{3}\mathrm{+}\mathrm{2}{\mathit{e}}_{\mathrm{P}}\\ & & {\mathit{g}}_{\mathit{i}}\mathrm{\left(}{\mathit{M}}_{\mathrm{P}}\mathit{,}{\mathit{e}}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{3}{\mathit{M}}_{\mathrm{p}}\mathit{/}\mathrm{(}{\mathit{e}}_{\mathrm{P}}\mathrm{+}\mathrm{0.001}\mathrm{)}}\mathrm{\times }{\mathrm{1}}^{\mathrm{◦}}\mathit{.}\end{array}$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 long-term 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 Type-II timescales due to planet-disc 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 single-planet 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 well-separated 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 K-prescription and the prescription of Paper I.

In N-body simulations, a K-prescription for the eccentricity damping timescale, $\frac{\mathit{ė}}{\mathit{e}}\mathrm{=}\mathrm{-}\mathit{K}\left|\hspace{0.17em}\frac{\mathit{ȧ}}{\mathit{a}}\right|\mathit{,}$(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 mean-motion resonance during Type-II migration can undergo an inclination-type resonance, if eccentricity damping is not too efficient. Libert & Tsiganis (2011a,b) generalized these results to three-planet systems, studying the establishment of three-planet 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 K-prescription 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 K-prescription has been developed by Tanaka & Ward (2004), using first-order, linear calculations which apply to low-mass 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 long-term effect. The simulations of Kley & Dirksen (2006) stated that a 3-Jupiter 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 long-term evolution of Jupiter-like planets in protoplanetary discs and found eccentricity growth for planets even smaller than 3-Jupiter 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 K-prescription 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 K-factor 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⁴−10⁵ yr (depending on the mass of the planet, see Fig. 4 of Paper I), while Type-II migration timescale is of the order of 10⁵ 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 three-planet initial set-up as in the present work but making use of the K-prescription 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 planet-planet 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 m1 = 1.98, m2 = 8.35, m3 = 3.8 MJup.

Furthermore, let us note that while the damping on the inclinations is either not considered or chosen in an arbitrary way in previous N-body 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 K-prescription where two K-factors 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 non-coplanar 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 inclination-type resonance in a Laplace 1:2:4 mean-motion resonance (DM model). The planetary masses are m1 = 2.66, m2 = 1.44 and m3 = 0.73 MJup.

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 planet-planet interactions that excite the eccentricity and inclination of giant planets on initially circular and coplanar orbits, and the planet-disc interactions that damp the same quantities. Two dynamical mechanisms producing inclination increase during the disc phase are discussed in the following: planet-planet scattering and inclination-type resonance. The frequency of these dynamical mechanisms will be discussed in Sect. 4.7.

The planet-planet scattering mechanism for planetary masses of m₁ = 1.98, m₂ = 8.35 and m₃ = 3.8 M_Jup is illustrated in Fig. 2. While the outer planet migrates, m₂ and m₃ are captured in a 1:3 mean-motion resonance, at approximately 2 × 10⁵ yr, as indicated by the libration of both resonant angles θ₁ = λ₂−3λ₃ + 2ω₂ and θ₂ = λ₂−3λ₃ + 2ω₃, 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₂ approaches a 3:7 commensurability with the inner planet, the whole system is destabilized, leading to planet-planet scattering. Indeed, heavy planets usually lead to dynamical instability before the establishment of a three-body mean-motion resonance (Libert & Tsiganis 2011b). The inner less massive body is then ejected from the system at 4 × 10⁵ 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 planet-planet 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 three-body resonance and the subsequent inclination excitation. In this example, the three planetary masses are m₁ = 2.66 M_Jup, m₂ = 1.44 M_Jup and m₃ = 0.73 M_Jup. The two inner planets are initially in a 1:2 mean-motion resonance. The convergent migration of the outer planet leads to the establishment of a three-body mean-motion resonance, the Laplace 1:2:4 resonance (since n₂:n₁ = 1:2 and n₃:n₂ = 1:2, the multiple-planet resonance is labeled as n₃:n₂:n₁ = 1:2:4), at approximately 2 × 10⁴ yr. The libration of the resonant angle φ_L = λ₁−3λ₂ + 2λ₃ 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₂), despite the damping exerted by the disc. Let us recall that there is no eccentricity excitation from the disc for low-mass planets. When the eccentricities are high enough, the system enters an inclination-type resonance, at approximately 1 × 10⁵ yr: the angles ${\mathit{\theta }}_{{\mathit{i}}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{=}\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\mathrm{-}\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathrm{\Omega }}_{\mathrm{1}}$ and ${\mathit{\theta }}_{{\mathit{i}}_{\mathrm{3}}^{\mathrm{2}}}\mathrm{=}\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\mathrm{-}\mathrm{4}{\mathit{\lambda }}_{\mathrm{3}}\mathrm{+}\mathrm{2}{\mathrm{\Omega }}_{\mathrm{2}}$ (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 inclination-type resonance. On the other hand, the eccentricities of the planets keep increasing as migration continues, in such a way that the system re-enters the inclination-type resonance at ~3 × 10⁵ yr, with both critical angles ${\mathit{\theta }}_{{\mathit{i}}_{\mathrm{1}}^{\mathrm{2}}}$ and ${\mathit{\theta }}_{{\mathit{i}}_{\mathrm{3}}^{\mathrm{2}}}$ 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 long-term 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 non-coplanar planetary systems. However, when the capture in the inclination-type resonance happens very rapidly, the gas disc is so massive (disc-dominated case) that it forces the planet to get back to the midplane in a very short timescale (~5 × 10⁴ 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 planet-planet 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 planet-planet 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 non-coplanar 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⁶ yr for the 4000 simulations of the CM model and 1.4 × 10⁶ yr for 5800 simulations of the DM model (see Table 1). This timescale is appropriate for the study of the influence of the disc-planet and planet-planet interactions during the disc phase. The systems emerging from the disc phase would most probably experience orbital rearrangement in the future due to planet-planet interactions. The issue of the long-term 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 constant-mass model preferably forms systems with a small number of planets, since migration is efficient during all the simulations (disc-dominated case favored) and pushes the planets closer to one another and closer to the parent star. However, the more realistic decreasing-mass 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.

Fig. 4 Normalized semi-major 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 Mp ∈ [0.65,10] MJup, black line). The bin size is Δlog (a) = 0.2.

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

Hereafter, we describe the orbital characteristics of the planetary systems formed by our mechanism. Semi-major 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 three-body 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⁶ 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 semi-major axis and mass.

4.1. Semi-major axes

We first investigate the effect of the gas disc on the final distribution of semi-major axis. Figure 4 shows the normalized semi-major 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 semi-major 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 semi-major 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 multi-planetary systems on the initial disc mass.

We further analyze the DM model semi-major axis distribution, investigating the influence of the initial disc mass in Fig. 5. The different panels present the normalized semi-major 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⁶ yr) and we show in Sect. 5 that orbital adjustments due to planet-planet 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) multi-planetary systems only. We first observe that removing the single-planet 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 multi-planetary 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 semi-major 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 semi-major axis and the eccentricity independently. The left panel of Fig. 7 shows the matching between our simulations and the observations in the semi-major axis vs. eccentricity graph. The results are in agreement with the study of Matsumura et al. (2010), combining N-body dynamics with hydrodynamical disc evolution.

Finally, the mean eccentricity for each system that ends up in a multi-planetary 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 mean-motion resonance capture, planet-planet 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 high-mass planets (M_p ≳ 5.5 M_Jup). This feature is included in the damping formulae derived in Paper I. For the three-planet systems, the mean eccentricities are smaller, since planet-planet scattering did not take place and only the other two mechanisms operated.

Fig. 7 Left panel: semi-major axis vs. eccentricity for our simulations with MDisc = 16 MJup (green dots) and the observations (black dots). Right panel: total mass vs. mean eccentricity for each multi-planetary configuration of our 16 MJup disc simulations. Black dots correspond to three-planet systems and green dots to two-planet systems.

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.

Fig. 9 Mutual inclination of the two-planet systems (in logarithmic scale) as a function of the total mass of the planets in the system for the four different disc masses.

4.3. Inclinations

Fig. 10 Initial and final Hill neighbor separation normalized distributions (DHill = (Δai,j) /RHill,ai,aj) of the systems of our simulations. The left column displays the two-planet systems and the right column the inner and outer pairs of the three-planet systems. Each panel represents a different initial disc mass. The bin size is ΔDHill = 1.2.

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 two-planet systems, since these systems have undergone a scattering event during the gas phase. Moreover, regarding the mutual inclination of two-planet 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 semi-major 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 planet-planet 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 close-in planets have eccentricities higher than 0.2. This suggests that planetary migration, including planet-planet 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: ${\mathit{R}}_{\mathrm{Hill}\mathit{,}{\mathit{a}}_{\mathit{i}}\mathit{,}{\mathit{a}}_{\mathit{j}}}\mathrm{=}{\left(\frac{{\mathit{m}}_{\mathit{i}}\mathrm{+}{\mathit{m}}_{\mathit{j}}}{\mathrm{3}{\mathit{M}}_{\mathrm{star}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\left(\frac{{\mathit{a}}_{\mathit{i}}\mathrm{+}{\mathit{a}}_{\mathit{j}}}{\mathrm{2}}\right)\mathit{,}$(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 two-planets are shown in the left column, while the right column shows both inner and outer pairs of the three-planet 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 three-planet 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 two-planet systems, we see that the higher the initial disc mass, the larger the mean separation.

4.6. Three-body resonances

As previously mentioned, 40% of the systems (in our DM model simulations) end up in a three-planet configuration (Table 3). Depending on the initial separation of the planets, the systems are generally captured in two- or three-body mean-motion 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 mean-motion resonance, ~ 25% are in a two-body mean-motion resonance and ~ 65% are in a three-body mean-motion resonance. This high percentage (~ 90%) of resonant systems is also observed in Matsumura et al. (2010). Table 4 shows that half of the three-body resonant systems are in a 1:2:4 Laplace configuration. The n₃:n₂:n₁ = 1:3:6 and n₃:n₂:n₁ = 2:5:10 resonances are the second most common configurations. These percentages being the same in both disc modelizations shows that unlike the semi-major axis distribution, the establishment of the resonant three-body configurations is not affected by the decrease of the gas disc.

Fig. 11 Examples of several scenarios producing mutual inclination excitation: a) planet-planet scattering followed by mean-motion resonance; b) orbital instability; c) orbital instability followed by three-body resonance; and d) three-body resonance followed by orbital instability.

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 inclination-type resonance and planet-planet scattering. We now aim at understanding the frequency of each mechanism, keeping in mind that the occurrence of both mechanims is possible during the long-term 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⁶ yr), three scenarios could have happened: either planet-planet 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 planet-planet scattering are rapidly damped by the disc and a mean-motion resonance capture of the two bodies induces the increase of the inclinations (Fig. 11, panel a), or a three-body inclination-type resonance excites the inclinations and consequently leads to planet-planet scattering and the ejection of a planet (similarly to Fig. 12). Mutual inclination is always observed in three-body 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 three-body resonance capture and a resonant excitation of the inclinations (Fig. 11, panel c), the system evolves smoothly in a three-body inclination-type resonance (similarly to Fig. 3), or the phase of three-body inclination-type 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 three-body mean-motion resonance captures, the other half being produced by orbital instability and/or planet-planet scattering. This emphasizes the importance of mean-motion resonance captures during the disc phase, on the final 3D configurations of planetary systems.

5. Long-term 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 planet-planet 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 well-separated and stable orbits.

In this section, we aim to investigate whether or not the long-term evolution of planetary systems produces significant changes on the final distribution of the orbital elements discussed hereabove. To study the long-term 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 long-term 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 three-body resonance on a long timescale. While the system is locked in a 1:2:6 resonance at the dispersal of the gas disc, a planet-planet scattering event finally takes place and the middle planet is ejected from the system at ~ 38 Myr. The planetary masses are m1 = 0.87, m2 = 1.39 and m3 = 8.43 MJup. The initial mass of the disc is 16 MJup.

Fig. 13 Left panel: normalized eccentricity distributions of the 16 MJup disc simulations for two integration timescales; 1.4 × 106 and 1 × 108 yr. The bin size is Δe = 0.05. Right panel: normalized semi-major axis distributions of the 16 MJup disc simulations, for the same two integration timescales. The bin size is Δlog (a) = 0.2.

Fig. 14 Normalized inclination distributions for two integration timescales (1.4 × 106 and 1 × 108 yr) and two initial disc masses (8 and 16 MJup). Only the inclinations higher than 1° are represented. The bin size is Δlog (i) = 0.2.

Figure 13 shows no significant change on the semi-major 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 three-planet systems are highly mutually inclined (I_mut> 10°) at 1.4 × 10⁶ yr and this percentage remains approximately the same for the long-term simulations. The situation is quite different for the two-planet 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 long-term evolution on 100 Myr. The same observation has recently been pointed out by Barnes et al. (2015), concerning planetary systems in mean-motion 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 Type-II migration, with the consistent eccentricity and inclination damping of Paper I (Bitsch et al. 2013), and planet-planet 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 constant-mass model (no gas dissipation during 0.8 Myr) and the decreasing-mass model (gas exponential decay with an e-folding 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 well-diversified 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 quasi-coplanar 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, inclination-type resonance and planet-planet 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 inclination-type resonance, orbital instability, and/or planet-planet 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 three-body mean-motion resonance captures. The long-term 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 by-product of our study, we found a very good agreement between our simulations and the observed population of extrasolar systems, in particular for the semi-major 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 planet-planet interactions during the migration inside the protoplanetary disc could account for most of the eccentricity excitation observed among exoplanets.

Acknowledgments

This work was supported by the Fonds de la Recherche Scientifique-FNRS 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

All Tables

All Figures

	Fig. 1 Normalized initial and final mass distributions of our simulations compared with the mass-distribution (in the range [0.65,10] MJup) of the observed giant planets with a > 0.1.
In the text

	Fig. 2 Inclination excitation due to planet-planet 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 m1 = 1.98, m2 = 8.35, m3 = 3.8 MJup.
In the text

	Fig. 3 Inclination excitation due to inclination-type resonance in a Laplace 1:2:4 mean-motion resonance (DM model). The planetary masses are m1 = 2.66, m2 = 1.44 and m3 = 0.73 MJup.
In the text

	Fig. 4 Normalized semi-major 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 Mp ∈ [0.65,10] MJup, black line). The bin size is Δlog (a) = 0.2.
In the text

	Fig. 5 Same as Fig. 4 for the different initial disc masses.
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 multi-planetary systems on the initial disc mass.
In the text

	Fig. 7 Left panel: semi-major axis vs. eccentricity for our simulations with MDisc = 16 MJup (green dots) and the observations (black dots). Right panel: total mass vs. mean eccentricity for each multi-planetary configuration of our 16 MJup disc simulations. Black dots correspond to three-planet systems and green dots to two-planet systems.
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.
In the text

	Fig. 9 Mutual inclination of the two-planet systems (in logarithmic scale) as a function of the total mass of the planets in the system for the four different disc masses.
In the text

	Fig. 10 Initial and final Hill neighbor separation normalized distributions (DHill = (Δai,j) /RHill,ai,aj) of the systems of our simulations. The left column displays the two-planet systems and the right column the inner and outer pairs of the three-planet systems. Each panel represents a different initial disc mass. The bin size is ΔDHill = 1.2.
In the text

	Fig. 11 Examples of several scenarios producing mutual inclination excitation: a) planet-planet scattering followed by mean-motion resonance; b) orbital instability; c) orbital instability followed by three-body resonance; and d) three-body resonance followed by orbital instability.
In the text

	Fig. 12 Destabilization of a three-body resonance on a long timescale. While the system is locked in a 1:2:6 resonance at the dispersal of the gas disc, a planet-planet scattering event finally takes place and the middle planet is ejected from the system at ~ 38 Myr. The planetary masses are m1 = 0.87, m2 = 1.39 and m3 = 8.43 MJup. The initial mass of the disc is 16 MJup.
In the text

	Fig. 13 Left panel: normalized eccentricity distributions of the 16 MJup disc simulations for two integration timescales; 1.4 × 106 and 1 × 108 yr. The bin size is Δe = 0.05. Right panel: normalized semi-major axis distributions of the 16 MJup disc simulations, for the same two integration timescales. The bin size is Δlog (a) = 0.2.
In the text

	Fig. 14 Normalized inclination distributions for two integration timescales (1.4 × 106 and 1 × 108 yr) and two initial disc masses (8 and 16 MJup). Only the inclinations higher than 1° are represented. The bin size is Δlog (i) = 0.2.
In the text