Highly inclined and eccentric massive planets
I. Planetdisc interactions
^{1} University of NiceSophia Antipolis / CNRS / Observatoire de la Côte d’Azur, Laboratoire Lagrange UMR 7293, Boulevard de l’Observatoire, BP4229, 06304 Nice Cedex 4, France
email: bertram.bitsch@oca.eu
^{2} NaXys, Department of Mathematics, University of Namur, 8 rempart de la Vierge, 5000 Namur, Belgium
^{3} Observatoire de Lille (LALIMCCE), CNRSUMR 8028, 1 impasse de l’Observatoire, 59000 Lille, France
Received: 30 August 2012
Accepted: 11 May 2013
Context. In the solar system, planets have a small inclination with respect to the equatorial plane of the Sun, but there is evidence that in extrasolar systems the inclination can be very high. This spinorbit misalignment is unexpected, as planets form in a protoplanetary disc supposedly aligned with the stellar spin. It has been proposed that planetplanet interactions can lead to mutual inclinations during migration in the protoplanetary disc. However, the effect of the gas disc on inclined giant planets is still unknown.
Aims. In this paper we investigate planetdisc interactions for planets above 1 M_{Jup}. We check the influence of three parameters: the inclination i, eccentricity e, and mass M_{p} of the planet. This analysis also aims at providing a general expression of the eccentricity and inclination damping exerted on the planet by the disc.
Methods. We perform threedimensional numerical simulations of protoplanetary discs with embedded highmass planets on fixed orbits. We use the explicit/implicit hydrodynamical code NIRVANA in 3D with an isothermal equation of state.
Results. We provide damping formulae for i and e as a function of i, e, and M_{p} that fit the numerical data. For highly inclined massive planets, the gap opening is reduced, and the damping of i occurs on timescales of the order of 10^{4} deg/year·M_{disc}/(0.01 M_{⋆}) with the damping of e on a smaller timescale. While the inclination of low planetary masses (<5 M_{Jup}) is always damped, large planetary masses with large i can undergo a Kozaicycle with the disc. These Kozaicycles are damped through the disc in time. Eccentricity is generally damped, except for very massive planets (M_{p} ~ 5 M_{Jup}) where eccentricity can increase for low inclinations. So the dynamics tends to a final state: planets end up in midplane and can then, over time, increase their eccentricity as a result of interactions with the disc.
Conclusions. The interactions with the disc lead to damping of i and e after a scattering event of highmass planets. If i is sufficiently reduced, the eccentricity can be pumped up because of interactions with the disc. If the planet is scattered to high inclination, it can undergo a Kozaicycle with the disc that makes it hard to predict the exact movement of the planet and its orbital parameters at the dispersal of the disc.
Key words: accretion, accretion disks / planets and satellites: formation / hydrodynamics / planetdisk interactions
© ESO, 2013
1. Introduction
In the solar system, the orbits of all the planets are nearly coplanar (within 4 degrees, except for Mercury). The ecliptic (the plane of the Earth’s orbit) is also close to the equatorial plane of the Sun: the spinorbit misalignment is only β_{⊕} = 7.5°. The low inclination of the massive planets with respect to the ecliptic is normally taken as an indication that planets form within a flattened protoplanetary disc, itself closely aligned with the stellar equator. The newly discovered Kepler30 system (SanchisOjeda et al. 2012) is even flatter, and confirms this view. However, exoplanets with strong spinorbit misalignment have been detected (e.g. β > 50°; Moutou et al. 2011a,b; Hébrard et al. 2011; Simpson et al. 2011). Considering that the plane of the past protoplanetary disc should be identical to the present stellar equator^{1}, the orbital plane of these planets must have been changed by some mechanism.
One process generally invoked to explain inclined orbits is scattering by multiple planets in the system after the protoplanetary disc has dissipated (e.g. Marzari & Weidenschilling 2002; Chatterjee et al. 2008; Jurić & Tremaine 2008). These works assume that unstable crowded systems are formed, and undergo planetplanet scattering after a relatively short time when the gas nebula dissipates. However, recent work suggests that unstable systems reach instability while still embedded in the gas disc (Lega et al. 2013). A second process is planetplanet interactions during migration in the protoplanetary disc (Thommes & Lissauer 2003; Libert & Tsiganis 2009, 2011a,b). During the gasdriven migration, the system can enter an inclinationtype resonance or the resonant configuration becomes unstable as the resonance excites the eccentricities of the planets and planetplanet scattering sets in. All this affirms the need for a better understanding of the interactions between giant planets and a gaseous protoplanetary disc when the orbit of the former is highly inclined with respect to the midplane of the later. Here, we study this phenomenon, in detail.
Tanaka & Ward (2004) have shown in linear studies that the inclination of a lowmass planet embedded in a disc is exponentially damped by planetdisc interactions for any nonvanishing inclination. Such results are formally valid only for i ≪ H/r. However, numerical simulations of more highly inclined planets have shown that the exponential damping might be valid up to i ≈ 2H/r. If the planet has an even greater inclination, the damping rates deviate from being exponential and it can be fitted by a di/dt ∝ i^{2} function (Cresswell et al. 2007; Bitsch & Kley 2011). However, for highmass planets, the linear regime is no longer valid. Marzari & Nelson (2009) considered Joviantype planets on inclined and eccentric orbits. They find highly inclined and eccentric planets with Jovian masses lose their inclination and eccentricity very quickly (on a timescale of the order of 10^{3} years) when entering the disc again (when i < H/r). Since a highly inclined planet is only slightly disturbed by the accretion disc (and vice versa), this kind of planet is only able to open a gap in the disc when the inclination drops to i < 10.0°.
Planetdisc interactions also influence the eccentricity of embedded planets, as has been shown by Goldreich & Tremaine (1980). It has been suggested, by performing linear analysis, that the planetary eccentricity can be increased through planetdisc interaction under some conditions (Goldreich & Sari 2003; Sari & Goldreich 2004; Moorhead & Adams 2008). They estimate that eccentric Lindblad resonances can cause eccentricity growth for gapforming planets. However, numerical simulations show that eccentricity in the disc is damped for a variety of masses (Cresswell et al. 2007; Moorhead & Ford 2009; Bitsch & Kley 2010).
For veryhighmass planets, on the other hand an eccentric instability in the disc can arise (Kley & Dirksen 2006). In turn, this eccentric disc can possibly increase the planetary eccentricity (Papaloizou et al. 2001; D’Angelo et al. 2006). However, this process can only explain the eccentricity of very massive (≈5–10 M_{Jup}) planets. XiangGruess & Papaloizou (2013) have recently studied the interactions between Jupitermass planets and circumstellar discs as well. However, they did not consider planets on eccentric orbits and they were using SPH simulations, while we use a gridbased code.
In this paper, we investigate planetdisc interactions for planets above 1 M_{Jup}, considering different inclination and eccentricity values. Our analysis also aims at deriving a formula for the change of eccentricity and inclination due to planetdisc interactions, in order to study the longterm evolution of systems with massive planets. Indeed, longterm evolution studies of planetary systems cannot be done with hydrodynamical simulations, as the computation time is too long, and NBody codes that consider the gravitational effects only are used. A correct damping rate of eccentricity and inclination is needed in order to simulate the evolution correctly. This study will be the topic of our Paper II.
We use isothermal threedimensional (3D) simulations to determine the change of inclination and eccentricity due to planet disc interactions. In Sect. 2 we describe the numerical methods used, as well as the procedure to calculate the forces acting on the embedded planets to determine di/dt and de/dt. In Sect. 3 we show di/dt and de/dt as a function of inclination i and eccentricity e, and provide fitting formulae. Additionally an observed oscillatory behaviour is discussed in this section. The implications for singleplanet systems are shown in Sect. 4.
2. Physical modelling
The protoplanetary disc is modelled as a 3D, nonselfgravitating gas whose motion is described by the NavierStokes equations. We use the code Nirvana (Ziegler & Yorke 1997; Kley et al. 2001), which uses the FARGOalgorithm (Masset 2000) and was described in our earlier work on planets on inclined orbits (Bitsch & Kley 2011). We note that the use of the FARGOalgorithm may not be straight forward in the case of highly inclined planets (i = 75.0°). Our test simulations, however, show that this algorithm can also be used in highly inclined planets, see Appendix. A. Here we treat the disc as a viscous medium in the locally isothermal regime. We do not use radiation transport, as we focus here on highmass planets that open a gap inside a disc, where the effects of heating and cooling of the disc are much less important than for lowmass planets (Kley et al. 2009). A more detailed description of the used code can be found in Kley et al. (2009).
2.1. Smoothing of the planetary potential
An important issue in modelling planetary dynamics in discs is the gravitational potential of the planet since this has to be artificially smoothed to avoid singularities. While in two dimensions a potential smoothing takes care of the otherwise neglected vertical extension of the disc, in 3D simulations the most accurate potential should be used. As the planetary radius is much smaller than a typical grid cell, and the planet is treated as a point mass, a smoothing of the potential is required to ensure numerical stability.
In Kley et al. (2009) two different kinds of planetary potentials for 3D discs have been discussed. The first is the classic ϵ_{sm}potential (1)Here M_{P} is the planetary mass, and d =  r − r_{P}  denotes the distance of the disc element to the planet. This potential has the advantage that it leads to very stable evolutions when the parameter ϵ_{sm} is a significant fraction of the Roche radius. The disadvantage is that for smaller ϵ_{sm}, which would yield a higher accuracy at larger d, the potential becomes very deep at the planetary position. Additionally, the potential differs from the exact 1/r potential even for medium to larger distances d from the planet.
To resolve these problems at small and large d simultaneously, the following cubicpotential has been suggested (Klahr & Kley 2006; Kley et al. 2009) (2)The construction of the planetary potential is such that for distances larger than r_{sm} the potential matches the correct 1/r potential. Inside this radius (d < r_{sm}) it is smoothed by a cubic polynomial. This potential has the advantage of exactness outside the specified distance r_{sm}, while being finite inside.
For 1 M_{Jup} and 5 M_{Jup} we use the cubic potential with r_{sm} = 0.8R_{H}. For the 10 M_{Jup} planet, we use the ϵ_{sm}potential with r_{sm} = 0.8R_{H}, with the Hill radius R_{H} given by (3)where a_{p} is the semi major axis of the planet, and M_{⋆} is the mass of the central star.
As the planetary mass increases, so does the amount of material accumulated near the planet. In order to resolve the gradients of density in that region correctly, a much higher resolution is required. Therefore, we change the cubic potential to the ϵ_{sm}potential for the 10 M_{Jup} planet. For the torque acting on the planets, the consequences are minimal, as we use a torque cutoff function in the Hill sphere of the planet, as described below. Additional information regarding the smoothing length can be found in Appendix A.
2.2. Initial setup
The 3D (r,θ,φ) computational domain consists of a complete annulus of the protoplanetary disc centred on the star, extending from r_{min} = 0.2 to r_{max} = 4.2 in units of r_{0} = a_{Jup} = 5.2 AU, where we put the planet. The planet is held on a fixed orbit during the evolution. The eccentricity of the planet can be e_{0} = 0.0, e_{0} = 0.2, or e_{0} = 0.4. We use 390 × 48 × 576 active cells for the simulations with 1 M_{Jup} and 260 × 32 × 384 active cells for 5 M_{Jup} and 10 M_{Jup}. This resolution is sufficient, as we still resolve the horseshoe width with a few grid cells for all planetary masses. The horseshoe width is defined for large planets as (Masset et al. 2006), where q is the planetstar mass ratio. Tests regarding the numerical resolution can be found in Appendix A.
In the vertical direction, the annulus extends 7° below and above the disc’s midplane, meaning 83° < θ < 97°. Here θ denotes the polar angle of our spherical polar coordinate system measured from the polar axis, therefore the midplane of the disc is at θ = 90.0°. We use closed boundary conditions in the radial and vertical directions. In the azimuthal direction, periodic boundary conditions are used. The central star has one solar mass M_{∗} = M_{⊙}, and the total disc mass inside [r_{min},r_{max}] is M_{disc} = 0.01M_{⊙}. The aspect ratio of the disc is H/r = 0.05. We use an α prescription of the viscosity, where (Shakura & Sunyaev 1973) with α = 0.005 ; Ω_{K} is the Kepler frequency ; denotes the isothermal sound speed, P the pressure, ρ the volume density of the gas, and H = c_{s}/Ω.
The models are initialised with constant temperatures on cylinders with a profile T(s) ∝ s^{1} with s = rsinθ. This yields a constant ratio of the disc’s vertical height H to the radius s. The initial vertical density stratification is given approximately by a Gaussian (4)Here, the density in the midplane is ρ_{0}(r) ∝ r^{1.5} which leads to a Σ(r) ∝ r^{− 1/2} profile of the vertically integrated surface density. In the radial and θdirection we set the initial velocities to zero, while for the azimuthal component the initial velocity u_{φ} is given by the equilibrium of gravity, centrifugal acceleration and the radial pressure gradient. This corresponds to the equilibrium configuration for a purely isothermal disc with constant viscosity. However, as the massive planets in the disc start to open gaps, the density and surface density profile get distorted.
2.3. Calculation of forces
To determine the change of orbital elements for planets on fixed inclined orbits, we follow Burns (1976) and compute the forces as described in Bitsch & Kley (2011). The gravitational torques and forces acting on the planet are calculated by integrating over the whole disc, where we apply a tapering function to exclude the inner parts of the Hill sphere of the planet. Specifically, we use the smooth (Fermitype) function (5)which increases from 0 at the planet location (d = 0) to 1 outside d ≥ R_{H} with a midpoint f_{b} = 1/2 at d = bR_{H}, i.e. the quantity b denotes the torquecutoff radius in units of the Hill radius. This torquecutoff is necessary to avoid large, probably noisy contributions from the inner parts of the Roche lobe and to disregard material that is possibly gravitationally bound to the planet (Crida et al. 2009). Here we assume b = 0.8, as a change in b did not influence the results significantly (Kley et al. 2009).
If a small disturbing force dF (given per unit mass) due to the disc is acting on the planet, the planet changes its orbit. This small disturbing force dF may change the planetary orbit in size (semimajor axis a), eccentricity e, and inclination i. The inclination i gives the angle between the orbital plane and an arbitrary fixed plane, which is in our case the equatorial plane (θ = 90°), which corresponds to the midplane of the disc. Only forces lying in the orbit plane can change the orbit’s size and shape, while these forces cannot change the orientation of the orbital plane. In Burns (1976) the specific disturbing force is written as (6)where each e represents the relevant orthogonal component of the unit vector. The perturbing force can be split into these components: R is radially outwards along r; T is transverse to the radial vector in the orbit plane (positive in the direction of motion of the planet); and N is normal to the orbit planet in the direction R × T.
Burns (1976) finds for the change of inclination (7)where the numerator is the component of the torque which rotates the specific angular momentum H = r × ṙ about the line of nodes (and which thereby changes the inclination of the orbital plane). The specific angular momentum H is defined as (8)The angle ξ is related to the true anomaly f by f = ξ − ω, with ω being the argument of periapsis and ξ describes the angle between the line of nodes and the planet on its orbit around the star. For the case of circular orbits, the argument of periapsis ω is zero.
The change of eccentricity is given by Burns (1976) as (9)where ϵ is the eccentric anomaly, which is given by (10)With this set of equations, we are able to calculate the forces acting on planets on fixed orbits and determine di/dt and de/dt.
3. Planets on inclined and eccentric orbits
In this section we investigate the changes of the planetary orbit due to planetdisc interactions. The planets are put in fixed orbits with inclinations ranging from i_{0} = 1.0° to i_{0} = 75°, with a total of ten different inclinations. For each inclination we also adopt three different eccentricities, which are e_{0} = 0.0, e_{0} = 0.2 and e_{0} = 0.4.
We note that the orbit of highly inclined planets is not embedded completely in the hydrodynamical grid, since the grid is only extended up to 7° above and below midplane. However, the density distribution in the vertical direction follows a Gaussian profile and for an aspect ratio of 0.05 we are at about 2.5σ at 7° so that the contribution of the gas can be neglected at larger θ.
3.1. Gaps in discs
The criterion for gap opening depends on the viscosity, the pressure, and the planetary mass (Crida et al. 2006). Giant planets (M ≳ 0.5 M_{Jup}) are generally massive enough to split the disc. However, the inclination of a giant planet plays a very important role in opening a gap as well, as can be seen in Fig. 1, where we display the surface density profile of discs with embedded 10 M_{Jup} planets on different inclinations.
Fig. 1 Surface density for disc simulations with 10 M_{Jup} planets in circular and eccentric orbits with different inclinations. The surface density is plotted after 400 planetary orbits. The evolution has reached an equilibrium state, meaning that the surface density does not change in time any more. 

Open with DEXTER 
Clearly, a lower inclination produces a much wider and deeper gap inside the disc. For larger inclinations, the gap opening is reduced, as the planet spends less and less time inside the disc to push material away from its orbit. Additionally, eccentric planets open up gaps that are less deep than their circular counter parts. This effect is very important for the damping of inclination and eccentricity, as an open gap inside the disc prolongs the damping timescale of inclination (Bitsch & Kley 2011) and of eccentricity (Bitsch & Kley 2010). Gap opening also indicates that linear analysis of the situation is no longer applicable.
In Fig. 2 we present slices in the x − zplane for the disc’s density for 10 M_{Jup} planets on inclinations of 1°, 20°, and 75° degrees. The inclinations correspond to those shown in the surface density plot (Fig. 1). Clearly the depth of the gap shown in the surface density is reflected in the 2D plots. Additionally, the density structures show no effects at the upper and lower boundaries because of boundary conditions, indicating that an opening angle of 7° is sufficient for our simulations.
Fig. 2 Density (in g/cm^{3}) of a r − θslice through the disc at the azimuth of an embedded 10 M_{Jup} planet on a fixed circular inclined orbit with i_{0} = 1.0° (top), i_{0} = 20° (middle), and i_{0} = 75.0° (bottom). The planet is at its lowest point in orbit (lower culmination) at the time of the snapshot, which was taken after 400 planetary orbits. We note the slightly different colour scale for each plot. The black line indicates the midplane of the grid to which the inclination of the disc is measured (see Sect. 3.2). 

Open with DEXTER 
3.2. Change of the disc structure
It has been known since several years that massive planets are able not only to open up a gap in the disc, but are also able to change the shape of the whole disc by turning it eccentric (Papaloizou et al. 2001; Kley & Dirksen 2006). Additionally, the inclination of the disc will change due to the interactions with the inclined planet. In this section, we discuss the impact of a massive planet on the eccentricity and inclination of the disc.
In Fig. 3 we display the eccentricity (top) and inclination (bottom) of the disc interacting with a 10 M_{Jup} planet with different inclinations (1° and 75°) and eccentricities. The calculations for deriving the eccentricity and inclination of the disc can be found in Appendix B.
For low planetary inclinations, the influence of the planet on the eccentricity of the protoplanetary disc is greater than on high planetary inclinations, simply, because the planet is closer to midplane and can therefore influence the eccentricity of the disc more strongly by pushing the material away. The eccentricity increase of the disc is stronger for planets in circular orbits than for planets that are already in an eccentric orbit. For highly inclined planets, the situation is reversed. The disc is most eccentric for planets that are already in an eccentric orbit and the disc is less eccentric for planets in circular orbits. Additionally, the eccentricity of the disc is highest close to the planet and drops with distance from the planet, independent of the inclination of the planet.
Fig. 3 Eccentricity (top) and inclination (bottom) of the disc with a 10 M_{Jup} planet influencing the disc structure after 400 planetary orbits. 

Open with DEXTER 
The inclination of the disc for the i_{0} = 1° planets is greater mostly around the planet’s location (at r = 1.0a_{Jup}) because the influence of the planet is strongest there. The inclination of the disc can be larger than the inclination of the planet. This is possible because the planet opens a gap inside the disc and pushes the material away from the planet (Fig. 1, top), which can also be seen in the 2D density configuration (top of Fig. 2). In the outer parts of the disc the disc remains noninclined.
For planets with high inclinations the situation is slightly different than for planets with low inclinations. The maximum of inclination is lower and there is no distinct maximum of inclination visible inside the planetary orbit (r < 1.0a_{Jup}) compared to the case of low inclinations. However, the outer parts of the disc show a nonzero inclination (which has a tendency to be larger for larger planetary eccentricities), which was not visible for the lowinclination planets. Additionally, they show a small peak of inclination at r ≈ 1.25a_{Jup}.
3.3. Change of orbital parameters
3.3.1. Eccentricity
As stated in Sect. 2.3, the forces acting on a planet on a fixed orbit can be calculated and then used to determine a rate of change for the inclination and eccentricity. The damping rates are taken when the planetdisc interactions are in an equilibrium state and do not change on average any more. The damping given by Eq. (7) varies strongly within the time of an orbit and slightly from one orbit to an other. Thus, we averaged it over 40 planetary orbits.
In Fig. 4 we present the change of eccentricity de/dt for planets of 1 M_{Jup}, 5 M_{Jup}, and 10 M_{Jup} on orbits with different inclinations and eccentricities. The change of de/dt has been studied in the past for coplanar lowmass planets (Cresswell et al. 2007; Bitsch & Kley 2010) and for highmass planets (Papaloizou et al. 2001; Kley & Dirksen 2006).
Fig. 4 Change of eccentricity de/dt for planets with 1 M_{Jup}, 5 M_{Jup}, and 10 M_{Jup} with different eccentricities. Points are results from numerical simulations, while lines indicate the fitting of the data. The 1 M_{Jup} planets have been evolved for 200 planetary orbits, the 5 M_{Jup} and 10 M_{Jup} planets have been evolved for 400 orbits. The forces used to calculate the data points have been averaged over 40 planetary orbits for all simulations. 

Open with DEXTER 
For low inclinations (i_{0} < 10°) the damping of eccentricity is stronger than for larger inclinations in the case of 1 M_{Jup}. The maximum damping rate is also dependent on the initial eccentricity e_{0}, where a larger e_{0} provides a faster damping. The damping of eccentricity is reduced significantly for larger inclinations i_{0} > 20°. As soon as the planet is no longer embedded in the disc, the damping reduces, as it is most efficient when the planet is inside the disc and not high above or below the disc for most of its orbit.
For low inclinations (i_{0} < 10°) and low eccentricities (e_{0} < 0.2), the 5 M_{Jup} planet opens up a large gap inside the disc. As the planet opens a gap inside the disc the damping is reduced because there is less material close to the planet to damp its orbit. For large initial eccentricities (e_{0} = 0.4), an increase of eccentricity is observable for low planetary inclinations. But for higher inclinations, the damping of eccentricity increases as well, because the planet does not open up such a deep gap (Fig. 1). However, for i_{0} > 40° the damping of eccentricity becomes smaller again, because the planet spends less and less time in the midplane of the disc where most of the disc material is, which is responsible for damping.
For even higher masses (10 M_{Jup}), we observe an eccentricity increase for low planetary inclinations for all nonzero eccentricities. But for larger inclinations (i_{0} > 15°), the eccentricity is damped again. The largest value of damping is at i_{0} ≈ 30° − 50°, depending on the planet’s eccentricity and is then reduced for higher inclinations again, following the trend described for the 5 M_{Jup} planet.
For large planets with low inclinations, the eccentricity of the planet can rise, which has been observed in Papaloizou et al. (2001) and Kley & Dirksen (2006). Papaloizou et al. (2001) stated that if the planet opens up a large gap, the m = 2 spiral wave at the 1:3 outer eccentric Lindblad resonance becomes dominant (because the order 1 resonances lie inside the gap) and induces eccentricity growth. However, they found an eccentricity increase only for M_{P} > 20 M_{Jup}, while our simulations indicate it clearly already for M_{P} > 5 M_{Jup} (Fig. 4, bottom). The differences between their 2D simulations and our 3D simulations can be the cause of the change in the required planetary mass for eccentricity growth.
Additionally, by embedding a highmass planet inside a disc, the disc can become eccentric, as shown in Sect. 3.2. The disc’s eccentricity is dependent on the planet’s inclination and slightly dependent on its eccentricity as well (see Fig. 3).
It seems that the coupling between a large disc eccentricity at r ≈ 1 − 1.5a_{P} and a large planetary eccentricity (the i_{0} = 1° with e = 0.4 case) results in a large force on the planet. This effect is increased as the planet in an eccentric orbit opens a small gap leaving more material at that location. This leads then to a greater increase of eccentricity for highly eccentric planets, compared to those with small eccentricity.
3.3.2. Inclination
In Fig. 5 we present the rate of change of inclination di/dt, presented in degrees per orbit, for planets with different masses and different eccentricities. For 1 M_{Jup} the inclination is damped for all initial inclinations. For increasing inclinations with i_{0} < 15° (smaller for increasing eccentricity), the damping of inclination increases. This increase is nearly linear, as has been shown for lowmass planets in theory (Tanaka & Ward 2004) and in numerical simulations (Cresswell et al. 2007; Bitsch & Kley 2011). The rates of inclination damping for zeroeccentricity planets are comparable to those stated in XiangGruess & Papaloizou (2013).
Fig. 5 Change of inclination di/dt for planets with 1 M_{Jup}, 5 M_{Jup}, and 10 M_{Jup} with different eccentricities. di/dt is in degrees per orbit at the planet’s location r_{P} = 1.0a_{Jup}. Points are results from numerical simulations, while lines indicate the fitting of the data. The 1 M_{Jup} planets have been evolved for 200 planetary orbits, the 5 M_{Jup} and 10 M_{Jup} planets have been evolved for 400 orbits. The forces used to calculate the data points have been averaged over 40 planetary orbits for all simulations. 

Open with DEXTER 
For i_{0} > 15°, the damping rate of the inclination is a decreasing function of inclination ; this is consistent with the planetdisc interaction being weaker when the planet spends more time farther from the midplane.
For 5 M_{Jup} the damping of inclination is almost the same as for the 1 M_{Jup} planet, but with a maximum at i_{0} ≈ 20°. However, there is a significant difference for high inclined and low eccentric planets: the inclination is not damped if i_{0} > 50°, but it increases for e_{0} < 0.1. This behaviour will be discussed in Sect. 3.4.
The 10 M_{Jup} planet shows the same general behaviour as the 5 M_{Jup} planet, but the inclination increase already sets in at , depending on e_{0}. Still, no inclination increase is observed in the high eccentricity simulations (e_{0} = 0.4). We also want to stress here that the damping rate significantly increases with increasing planetary eccentricity for all planetary masses.
The increase of inclination for highmass planets due to interactions with the disc has been studied in Lubow & Ogilvie (2001). They state that the 1:3 meanmotion resonance also acts to increase inclination. This resonance is at r_{r} = 2.08r_{P}, which clearly is not inside an open gap in the case of i_{0} = 75° (see Fig. 1). However, the resonances closer to the planet (1:2 and 2:3) are also not completely inside the gap, so that there should be some damping effects, but the damping of inclination through these resonances is weaker than the increase from the 1:3 resonance because in total the inclination increases for high inclined planets (Fig. 5). Lubow & Ogilvie (2001) also used small i_{0} for their calculations in order to apply linear theory, which does not apply for large inclinations. The situation for our high inclination planets might therefore be completely different from their calculations.
3.4. Moving planets in discs
3.4.1. Shortterm evolution
In order to verify the results of inclination and eccentricity change, we present in this section simulations of planets evolving freely in the disc. The planets are moving because of the influences of the discs forces. We present here several interesting cases for planets with high inclinations. The first case is for 5 M_{Jup} and 10 M_{Jup} with an inclination of i_{0} = 40° and i_{0} = 75° in circular orbits with an evolution time of 80 planetary orbits. In Fig. 6 the evolution of inclination with time is presented for the two different planets and inclinations. The evolution is nearly identical, as was predicted by the measured forces for the planets on fixed orbits, which is shown by the solid lines (rates from Fig. 5).
Fig. 6 Evolution of inclination of 5 M_{Jup} and 10 M_{Jup} planets with an initial inclination of i_{0} = 40.0° (top) and i_{0} = 75.0° (bottom) in circular orbits. The simulations have been restarted with a moving planet after the disc was evolved for a fixed planet for 400 planetary orbits. The time index has been reset to zero and the lines correspond to the expected damping rates from Fig. 5. 

Open with DEXTER 
One should be aware, however, that by keeping the planet in a fixed orbit, angular momentum in the system is not conserved because, for example, the inclination of the disc is rising (see Fig. 3) while the planet remains in a fixed orbit. The effect of conserving angular momentum is not a problem for lowmass planets, where the measured forces match perfectly with the inclination damping rates for moving planets (Bitsch & Kley 2011), but for big planets of several Jupiter masses this can lead to small differences because the angular momentum transfer from disc to planet and viceversa is much larger.
3.4.2. Longterm evolution
The longterm evolution of planets with different inclinations and eccentricities is displayed in Fig. 7. At the beginning of the evolution, the change of inclination and eccentricity matches those presented in Figs. 4 and 5 for planets in fixed orbits. However, the evolution after the initial orbits is quite different from what was expected by the previous simulations.
Fig. 7 Longterm evolution of planets with different inclinations, eccentricities, and masses in discs. The simulations are started from the equilibrium structures with fixed planets, where the planets are then released and allowed to move freely inside the disc. The top plot features the inclination of the planets, while the bottom plot shows the eccentricity of the planets. In the beginning the changes of eccentricity and inclination match the ones displayed in Figs. 4 and 5. 

Open with DEXTER 
In the 10 M_{Jup}, e_{0} = 0.0, i_{0} = 75° case, the inclination initially increases slightly with a rate that corresponds to the predicted rate (see also Fig. 6). At the same time, the eccentricity of the planet increases and after about 250 orbits it reaches e ≈ 0.25. This eccentricity corresponds to inclination damping (Fig. 5), which is what happens in the evolution of the planet: the inclination drops. However, the eccentricity still increases at the same time, which was not predicted by the analysis of planets in fixed orbits (Fig. 4). The eccentricity then rises until a peak of e ≈ 0.9, where it starts to drop again. At the same time, the inclination decrease stops and the inclination starts to rise again. As soon as the eccentricity has dropped to e ≈ 0.4, the inclination starts to decrease again.
This exchange of inclination and eccentricity is representative of the Kozai mechanism, introduced initially to describe the evolution of a highly inclined asteroid perturbed by Jupiter (Lidov 1962; Kozai 1962). A similar Kozai mechanism affects the orbits of highly inclined planets with respect to a disc (Terquem & Ajmia 2010; Teyssandier et al. 2013). For inclinations above a critical value, the gravitational force exerted by the disc on the planet produces Kozai cycles where the eccentricity of the planet can be pumped to high values, in antiphase with its inclination. We note that the Kozai mechanism is visible in the given computation time because of the high mass values considered in our study (5 and 10 M_{Jup}), comparable to the total mass of the disc (0.01M_{⊙}). Indeed high masses induce faster dynamical evolution.
When the planet starts at a larger initial eccentricity (e_{0} = 0.2 or e_{0} = 0.4), the general behaviour is similar as can be seen in Fig. 7, but the first Kozai cycle occurs earlier. Circular orbits at high inclination constitute an unstable equilibrium of the secular dynamics, so the evolution at zero initial eccentricity remains for a while close to the separatrix associated with the equilibrium (Libert & Henrard 2007). The Kozai effect does not act for initial inclinations smaller than a critical value (i_{0} < 40° in the restricted problem of Kozai 1962). We therefore also display a planet with i_{0} = 40°, and show that the eccentricity and inclination oscillations are significantly reduced.
Even if eccentricity can be pumped to high values, the Kozai mechanism only postpones the alignment with the disc and the circularization of the orbit induced by damping forces of the disc on the planet (given in Figs. 4 and 5). As clearly shown by the evolution of the planet with e_{0} = 0.2, Kozai cycles repeat with reduced intensity. The drop of inclination is much larger than the raise of inclination, after the eccentricity increase/decrease cycle. These results are in agreement with Teyssandier et al. (2013), showing that lowmass planets would remain on inclined and eccentric orbits over the disc lifetime, while higher mass planets would align and circularize. We also illustrate in Fig. 7 the influence of the planet mass by considering a planet of 5 M_{Jup}(e = 0.0): the eccentricity value reached during the second cycle of inclination increase (at 450 orbits) is higher for the 5 M_{Jup} planet, as expected.
The effect of Kozai oscillations between a disc and planet was also stated in XiangGruess & Papaloizou (2013) however, XiangGruess & Papaloizou (2013) were not able to resolve a full Kozai cycle, probably because their massratio between planet and disc is smaller than in our case. This shows that for i_{0} > 40°, the measure of the forces on a planet on a fixed orbit is not relevant. In this case, damped Kozai oscillations will govern the longterm evolution of the orbital parameters. This phenomenon can be of crucial importance for the study of the fate of planets scattered on highinclination orbits.
3.5. Fitting for e and i
In Figs. 5 and 4 we provided the change of di/dt and de/dt for different planetary masses. In these plots, lines indicate a fit for these data points. We now present the fitting formulae, which depend on the planet mass M_{P}, the eccentricity e_{P}, and inclination i_{P}. The inclination i_{P} used in the presented formulae is given in degrees, as is the resulting di/dt. As discussed in the previous section, these formulae are only relevant for i_{0} < 40° where no complex cycles are observed. Therefore, in fitting our parameters we have ignored the data points corresponding to high inclinations, in particular the ones showing inclination increase. This applies to the fitting of inclination and eccentricity.
As can be seen in the figures, the results of the numerical simulations are all but smooth. Therefore, one should not expect the fit to be very accurate with simple functions. However, our goal is to catch the big picture, and to provide an acceptable order of magnitude of the effect of the disc on the inclination and eccentricity. In log scale, the data appear relatively close to an increasing power law of i_{P} for small i_{P}, and a decreasing power law of i_{P} for large i_{P}. Therefore, we base our fits on the general form for the damping rates (11)where b is positive and d is negative. This way, for small i_{P}, , and for large i_{P}, . The coefficients a, b, c, and d depend on the planet mass and eccentricity, and are fitted to the data as follows. The damping rate also has to be linear dependent on the disc mass M_{disc}/M_{⋆}, as our simulations linearly scale with the gas density.
3.5.1. Eccentricity
We do not want (de/dt) to tend to zero when i_{P} tends to zero. A pure increasing power law of i_{P} is inappropriate here. The damping function will be given by (12)where i_{D} is a small inclination so that for i_{P} ≈ 0, . We are using degrees in Eq. (12), where is the planet mass in Jupiter masses. For small eccentricities, it is expected that e_{P}/(de/dt) = τ_{e} is constant. This makes the coefficient a proportional to . We find that de/dt is well fitted by the above general form using the coefficients (13)The second degree polynomial function of in the expression of a is just a refinement, its value being between 30 and 50 for . We note, however, that it is negative for so this expression only applies for , but this covers the range of giant planets. To describe the change of eccentricity we add a second function , which describes the eccentricity increase for highmass planets. The damping and excitation of e_{P} are two different mechanisms that add on the planet, and one of them finally dominates, setting the sign of de/dt. The fits in Fig. 4 are the added functions.
For we use the result of Papaloizou et al. (2001) who calculated the eccentricity excitation for i_{0} = 0° high massplanets. Our calculation is presented in Appendix C and gives (14)Then, we find that this excitation decreases with i as a Gaussian function, finally making (15)In principle planets with M_{P} < 5 M_{Jup} and e < 0.3 do not show any signs of eccentricity increase and the Gaussian function should not be added in that case. However, the function is designed to scale with the planetary mass, so that lower mass planets are not affected by it. The change of eccentricity is then given by the sum of ℱ_{e} and .
3.5.2. Inclination
In the case of inclination damping data, we notice that the decreasing power law dominates actually before the intersection with the increasing power law; thus, we multiply the term in our general formula by a Gaussian function of i_{P} centred on 0°, so that this term is not affected for small i_{P} but vanishes more quickly than is natural. It allows our fitting formula to catch the peak of damping in inclination observed around 5 to 20 degrees in Fig. 5. For small i_{P}, di/dt should be close to linear in i_{P}, so the coefficient b should be close to 1. The damping function for inclination ℱ_{i} is then given, in degrees per orbit, by (16)We note that the expression for coefficient c_{i} is clearly not valid for e > 2/3.
From our formulae for de/dt and di/dt we can now estimate how the eccentricity and inclination of a planet will evolve for all e_{P} and i_{P}. In Fig. 8 the di/dt for different inclinations and eccentricities for 5 M_{Jup} and 10 M_{Jup} according to the formulae is presented. In Fig. 9 the de/dt for the same two planetary masses is plotted.
Fig. 8 Top: di/dt for a 5 M_{Jup} planet with different eccentricities and inclinations. The values of di/dt have been determined with the formula given in Sect. 3.5. Bottom: same, but for 10 M_{Jup}. We note the different colour coding as the change is dependent on the planetary mass. 

Open with DEXTER 
Figure 8 clearly indicates that the damping rate of inclination is highest for planets with a large eccentricity that are moderately inclined above the midplane (i_{P} ≈ 15°). The inclination damping rate indicates that planets that are scattered during the gas disc phase in orbits with moderate inclination (i_{P} < 40°), would lose their inclination well within the gas dispersal of the disc.
As already indicated in Fig. 4, the eccentricity is always damped for high inclinations. For high planetary masses, the eccentricity of the planet can increase for low planetary inclinations because of interactions with disc. We find an eccentricity increase for both highmass cases, but the increase of eccentricity declines with increasing eccentricity and inclination. Additionally, the threshold of e_{P} and i_{P} for which eccentricity can increase is larger for higher mass planets, which is indicated by the white line in Fig. 9 that represents the transition from eccentricity increase to decrease. Below the line eccentricity increases, above the line eccentricity decreases.
Fig. 9 Top: de/dt for a 5 M_{Jup} planet with different eccentricities and inclinations. The values of de/dt have been determined with the formula given in Sect. 3.5. Bottom: same, but for 10 M_{Jup}. The white line in the figure indicates the transition between eccentricity increase and eccentricity damping. Below the white line, the eccentricity increases, above the line eccentricity decreases. We note the different colour coding as the change is dependent on the planetary mass. 

Open with DEXTER 
4. Application to singleplanet systems
The movement of a single planet in the disc can only be predicted if i < 40° and e < 0.65 as the planet would undergo a Kozaioscillation for larger i. Additionally, the fitting formula might not be totally accurate for e > 0.5, since our simulations only cover an eccentricity space of up to e = 0.4. In Fig. 10 the trajectory of the 10 M_{Jup} planet with i_{0} = 75° and e_{0} = 0.4, which was shown in Fig. 7 is displayed. This illustrates that the movement of the planet is a complex process as long as the Kozaioscillations are still operational, but as soon as i < 40°, the planet loses inclination, which is then not converted back into eccentricity. The planet is damped towards midplane on a nonzero eccentricity. This nonzero eccentricity will actually hold in this case (see Sect. 3.3.1).
A typical damping rate of de/dt = 0.001 /orbit would suggest that the planet will lose ≈0.085 in eccentricity in the period of 10^{4} years. A typical damping rate of di/dt = 0.01 deg/orbit indicates that the planet will lose ≈8.5° of inclination in 10^{4} years.
The important parameters for the evolution of the orbit of a planet are the damping timescales τ_{e} = e/(de/dt) and τ_{i} = i/(di/dt). We find that e/ℱ_{e} is much smaller than i/ℱ_{i} for i > 10 − 20°, depending on the planet mass and eccentricity. Thus, planets scattered on highly inclined orbits will follow a certain pattern. While the inclination is damped slowly and still high, the eccentricity is damped to zero. After the inclination is damped further, the eccentricity of the planet can rise because of interactions with the disc (if e is below the white line in Fig. 10). Finally the inclination is damped to zero and the planet remains with a nonzero eccentricity. This is illustrated by the blue lines in Fig. 10 that correspond to calculated trajectories of 10 M_{Jup} planets.
Nevertheless, this suggests that at the time of the disc dispersal, the favoured endstate for the planet’s evolution is an eccentric orbit in midplane of the disc. This implies that the scattering process of inclined planets must have taken place after the gas is depleted or gone.
5. Summary
We have presented the evolution of eccentricity e and inclination i of highmass planets (M_{P} ≥ 1 M_{Jup}) in isothermal protoplanetary discs. The planets have been kept on fixed orbits around the host star, and the forces from the disc acting onto the planet have been calculated. By using these forces, a change of de/dt and di/dt has been determined.
Inclination and eccentricity are in general damped by the interactions with the disc. For 1 M_{Jup} the damping rate of e and i is highest for only very small inclinations (i_{0} ≈ 3°), while the maximal damping rate is shifted to larger inclinations for more massive planets. As the more massive planets carve deeper gaps inside the disc, the damping interactions with the disc are reduced. But for larger inclinations, the planet can feel the full damping potential of the disc and is therefore damped in e and i at a faster rate.
There are two exceptions. The first is for lowinclination planets with a sufficient mass (M_{P} > 4 − 5 M_{Jup}). In this case, the interactions of the planet with the disc result in an increase of eccentricity of the planet, which has already been observed and studied (Papaloizou et al. 2001; Kley & Dirksen 2006). However, our 3D results predict an increase of eccentricity for lower planetary masses than the previous studies.
The second exception arises for massive planets (M_{P} ≈ M_{disc}, in our case for M_{P} > 5 M_{Jup}) on high initial inclinations (i_{0} > 40°). In the longterm evolution of the planet, eccentricity can increase, while inclination is damped and viceversa. The planet undergoes a Kozaicycle with the disc, but in time the oscillations of the planet in e and i diminish, as e and i get damped by the disc at the same time. The planet will end up in midplane through the interactions with the disc.
Fig. 10 Evolution of e and i of the 10 M_{Jup} planet (shown in Fig. 7) with i_{0} = 75° and e_{0} = 0.4 in the ei plane (black line). The background is the extended formula of the fit for de/dt and the white line marks the transition between eccentricity increase and damping as in Fig. 8. The blue lines indicate calculated trajectories of 10 M_{Jup} planets from the fitting formulae. 

Open with DEXTER 
In Sect. 3.5 we provided formulae for di/dt and de/dt for highmass planets, which we fitted to the numerical hydrodynamical simulations. The formulae can now be used to calculate the longterm evolution of planetary systems during the gas phase of the disc with NBody codes. However, we recommend not using the fitting formula, if the planetary eccentricity is e > 0.65 and if i > 40° (because of the Kozai interactions, a fit that can be used for the longterm evolution of planets is hard to predict).
In the end, the planet’s inclination will be damped to zero. Lowmass planets (M_{P} < 4 − 5 M_{Jup}) will end up in circular orbits in the midplane of the disc, while higher mass planets (M_{P} > 5 M_{Jup}) will pump their eccentricity to larger values because of interactions with the disc. This implies that the scattering process of inclined planets must have taken place after the gas is well depleted.
The influence of the gaseous protoplanetary disc on the inclination is also of crucial importance, if multiple planets are present in the disc that excite each other’s inclination during their migration (Libert & Tsiganis 2009, 2011a). The influence of the disc on the longterm evolution of multibody systems will be studied in a future paper.
This is generally accepted, but is actually the subject of debate (see e.g. Cébron et al. 2011; Batygin 2012).
Acknowledgments
B. Bitsch has been sponsored through the Helmholtz Alliance Planetary Evolution and Life. The work of A.S. Libert is supported by an FNRS Postdoctoral Research Fellowship. The calculations were performed on systems of the Computer centre (ZDV) of the University of Tübingen and systems operated by the ZDV on behalf of bwGRiD, the grid of the Baden Württemberg state. We thank the referee Willy Kley for his useful and helpful remarks that improve the paper.
References
 Batygin, K. 2012, Nature, 418 [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]
 Burns, J. A. 1976, Am. J. Phys., 44, 944 [NASA ADS] [CrossRef] [Google Scholar]
 Cébron, D., Moutou, C., Le Bars, M., Le Gal, P., & Farès, R. 2011, in Detection and Dynamics of Transiting Exoplanets, 11 [Google Scholar]
 Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580 [NASA ADS] [CrossRef] [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., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., Baruteau, C., Kley, W., & Masset, F. 2009, A&A, 502, 679 [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]
 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]
 Hébrard, G., Ehrenreich, D., Bouchy, F., et al. 2011, A&A, 527, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jurić, M., & Tremaine, S. 2008, ApJ, 686, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Klahr, H., & Kley, W. 2006, A&A, 445, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., & Dirksen, G. 2006, A&A, 447, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., D’Angelo, G., & Henning, T. 2001, ApJ, 547, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Lega, E., Morbidelli, A., & Nesvorny, D. 2013, MNRAS, in press [Google Scholar]
 Libert, A.S., & Henrard, J. 2007, Icarus, 191, 469 [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]
 Lidov, M. L. 1962, Planetary and Space Science, 719 [Google Scholar]
 Lubow, S. H., & Ogilvie, G. I. 2001, ApJ, 560, 997 [NASA ADS] [CrossRef] [Google Scholar]
 Marzari, F., & Nelson, A. F. 2009, ApJ, 705, 1575 [NASA ADS] [CrossRef] [Google Scholar]
 Marzari, F., & Weidenschilling, S. J. 2002, Icarus, 156, 570 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Masset, F. S., D’Angelo, G., & Kley, W. 2006, ApJ, 652, 730 [NASA ADS] [CrossRef] [Google Scholar]
 Moorhead, A. V., & Adams, F. C. 2008, Icarus, 193, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Moorhead, A. V., & Ford, E. B. 2009 [arXiv:0904.3336] [Google Scholar]
 Moutou, C., Díaz, R. F., Udry, S., et al. 2011a, A&A, 533, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moutou, C., Hébrard, G., & Bouchy, F. 2011b, in EPSCDPS Joint Meeting 2011, 596 [Google Scholar]
 Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 SanchisOjeda, R., Fabrycky, D., Winn, J. N., Barclay, T., et al. 2012, Nature, 487, 449 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Sari, R., & Goldreich, P. 2004, ApJ, 606, L77 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Simpson, E. K., Pollacco, D., Cameron, A. C., Hébrard, G., et al. 2011, MNRAS, 414, 3023 [NASA ADS] [CrossRef] [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., Terquem, C., & Papaloizou, J. 2013, MNRAS, 428, 658 [NASA ADS] [CrossRef] [Google Scholar]
 Thommes, E. W., & Lissauer, J. J. 2003, ApJ, 597, 566 [NASA ADS] [CrossRef] [Google Scholar]
 XiangGruess, M., & Papaloizou, J. C. B. 2013, MNRAS, 431, 1320 [NASA ADS] [CrossRef] [Google Scholar]
 Ziegler, U., & Yorke, H. 1997, Comput. Phys. Comm., 101, 54 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Additional information on numerics
In principle, a fast vertical movement (more than one grid cell per timestep) through the grid could cause problems with the Fargo algorithm, as Fargo shifts the grid cells for several cells azimuthally and the effects of the planet on the gas might get corrupted. In Fig. A.1 we present the evolution of the normal component of the disturbing force F_{N} (which has been averaged over one orbit) of planets with i = 75° on circular orbits. The two simulations shown feature different timestep lengths. For the simulation with larger timestep, the planet moves through about one vertical grid cell in each time step. For the shorter timestep, three timesteps are needed to cover the vertical extent of one grid cell. The evolution of F_{N} seems to be identical, indicating that the length of the timestep is not of crucial importance here also because we use a rotating frame so that the planet is on a fixed position inside the numerical grid where the Fargo algorithm does not shift grid cells for r ≈ a_{P}.
The smoothing of the planetary potential is crucial for avoiding singularities at the planet’s location. In Sect. 2.1 the numerical potential for the planets was introduced. Of crucial importance here is the smoothing length r_{sm}. A smaller smoothing length r_{sm} leads to a deeper planetary potential. This in turn leads to a larger accumulation of mass at the planet’s location, but this increase in density near the planet can be very high for large planets, especially in the isothermal case. This increase of density near the planet is unphysical, as normally the temperature and pressure gradients should stop the accumulation of gas at some point, which is not possible in the isothermal case. In this situation, the density can become so large that the gradients of density near the planet become too steep and the timestep inside the code collapses down to very small values, which makes an integration over several orbits impossible. We therefore return to the ϵpotential for the 10 M_{Jup} planet.
In Fig. A.2 we present the inclination damping for 5 M_{Jup} planets in circular orbits for two different smoothing lengths, r_{sm} = 0.8 and r_{sm} = 0.5. Changing the planetary potential seems to influence the damping of inclination by up to ± 15%, but the general trend is the same. Even with a deeper planetary potential, the inclination of a planet seems to increase for large initial inclinations. The main difference seems to be that no inclination increase can be observed for the i_{0} = 55° case with a smoothing length of r_{sm} = 0.5. This has also been observed for the 10 M_{Jup} planet where the difference between the depth of the two potentials is supposed to be stronger (as we change from the ϵ to the cubic potential), but the trend is the same as for the 5 M_{Jup} planet. For i_{0} = 75° the inclination seems to increase for models of planets in fixed orbits for both 5 M_{Jup} and 10 M_{Jup}. We therefore conclude that the general trend is conserved regardless of the chosen planetary potential and smoothing length.
Fig. A.1 F_{N} for 5 M_{Jup} planets with i = 75° on circular orbits. The planet’s feature different time steps, as indicated in the plot. F_{N} has been averaged over one running orbit. 

Open with DEXTER 
Fig. A.2 Change of inclination di/dt for 5 M_{Jup} and 10 M_{Jup} planets in circular orbits for two different smoothing length for the planetary potential. The cubic potential is used for all displayed simulations, unless stated otherwise. 

Open with DEXTER 
Fig. A.3 F_{N} for a 10 M_{Jup} planet with i = 3° in a e_{0} = 0.4 orbit. F_{N} has been averaged over one running orbit. 

Open with DEXTER 
In order to find the sufficient numerical resolution for our simulations of inclination damping, we have performed several resolution tests. In Fig. A.3 we present the results of these tests. The plotted quantity F_{N} has been averaged over one running orbit. Keep in mind that F_{N} has been averaged over 40 orbits to determine the change of inclination in the end. The simulations feature a 10 M_{Jup} planet with i = 3°, so it is well embedded inside the disc. The numerical resolution of the grid has been changed from 260 × 32 × 384 to 390 × 48 × 576. As the crucial force F_{N} for inclination damping gives the same results for both resolutions, we use the lower resolution for our simulation with confidence.
Appendix B: Eccentricity and inclination of the disc
To determine the eccentricity and inclination of the disc, we take a massweighted average of the eccentricity of all grid cells. To compute the eccentricity we take the specific total energy (in mass units) (B.1)where v is the velocity vector of a given grid cell and the radial component towards the grid cell. The total energy is given by (B.2)where a is the semimajor axis towards the grid cell. With that, we can compute a(B.3)With a we can now compute the eccentricity e of each grid cell: (B.4)where L_{spec.} = r × v is the specific angular momentum of each grid cell. To get an estimate of the eccentricity of the disc, we make a massweighted average of the eccentricity of each grid cell (averaged in azimuthal and polar coordinates) in order to get e_{disc}(r): (B.5)where m_{θφ} is the mass of the grid cell.
Because we use spherical coordinates r, θ, φ for the inclination, we have to transform L_{spec.} first into Cartesian coordinates in order to calculate the mass average. This has to be done because each product L = r × v is given in a different local coordinate system of each grid cell, but for the average the angular momentum vectors should always be in the same coordinate frame. The angular momentum vector is given in the two coordinate systems by (B.6)where (B.7)with the angles θ and φ of the grid cell, which differ for each grid cell. This gives us for L in Cartesian coordinates (B.8)For the inclination of the disc we now take a massaveraged specific angular momentum (averaged in polar and azimuthal direction) (B.9)where the subscript c denotes the grid cell number, and m_{c} the corresponding mass of the grid cell. Now we can compute the angle between L_{av.} and the zaxis, which gives us the averaged inclination at each ring of the disc.
Appendix C: Increase of eccentricity
We follow Papaloizou et al. (2001) to calculate the maximum value of eccentricity increase for highmass planets, as it is given by A_{g} in Eq. (15). In Papaloizou et al. (2001) the increase of eccentricity is calculated through the growth rates of the modes of the Lindblad resonance, which is given by (C.1)where (C.2)where e_{d} is the disc’s eccentricity and r_{P} the planetary distance to star. The integral basically gives the disc mass, which is comparable to the planet’s mass, but as e_{d} is much smaller than e_{P} (see Fig. 3), the term concerning the disc eccentricity is much smaller than the term concerning the planetary eccentricity. We therefore choose to neglect it in our estimate of the eccentricity increase. We then get (C.3)which leads to (C.4)As also (C.5)where we set e_{d} = 0.0 and because we are only interested in a first order estimate of the eccentricity increase from the disc. With we find for ė_{P} = 2e_{P}γ(C.6)which gives the increase of eccentricity for a planet orbiting in the midplane of the disc, which fits quite well with the results of our simulations (Fig. 4).
All Figures
Fig. 1 Surface density for disc simulations with 10 M_{Jup} planets in circular and eccentric orbits with different inclinations. The surface density is plotted after 400 planetary orbits. The evolution has reached an equilibrium state, meaning that the surface density does not change in time any more. 

Open with DEXTER  
In the text 
Fig. 2 Density (in g/cm^{3}) of a r − θslice through the disc at the azimuth of an embedded 10 M_{Jup} planet on a fixed circular inclined orbit with i_{0} = 1.0° (top), i_{0} = 20° (middle), and i_{0} = 75.0° (bottom). The planet is at its lowest point in orbit (lower culmination) at the time of the snapshot, which was taken after 400 planetary orbits. We note the slightly different colour scale for each plot. The black line indicates the midplane of the grid to which the inclination of the disc is measured (see Sect. 3.2). 

Open with DEXTER  
In the text 
Fig. 3 Eccentricity (top) and inclination (bottom) of the disc with a 10 M_{Jup} planet influencing the disc structure after 400 planetary orbits. 

Open with DEXTER  
In the text 
Fig. 4 Change of eccentricity de/dt for planets with 1 M_{Jup}, 5 M_{Jup}, and 10 M_{Jup} with different eccentricities. Points are results from numerical simulations, while lines indicate the fitting of the data. The 1 M_{Jup} planets have been evolved for 200 planetary orbits, the 5 M_{Jup} and 10 M_{Jup} planets have been evolved for 400 orbits. The forces used to calculate the data points have been averaged over 40 planetary orbits for all simulations. 

Open with DEXTER  
In the text 
Fig. 5 Change of inclination di/dt for planets with 1 M_{Jup}, 5 M_{Jup}, and 10 M_{Jup} with different eccentricities. di/dt is in degrees per orbit at the planet’s location r_{P} = 1.0a_{Jup}. Points are results from numerical simulations, while lines indicate the fitting of the data. The 1 M_{Jup} planets have been evolved for 200 planetary orbits, the 5 M_{Jup} and 10 M_{Jup} planets have been evolved for 400 orbits. The forces used to calculate the data points have been averaged over 40 planetary orbits for all simulations. 

Open with DEXTER  
In the text 
Fig. 6 Evolution of inclination of 5 M_{Jup} and 10 M_{Jup} planets with an initial inclination of i_{0} = 40.0° (top) and i_{0} = 75.0° (bottom) in circular orbits. The simulations have been restarted with a moving planet after the disc was evolved for a fixed planet for 400 planetary orbits. The time index has been reset to zero and the lines correspond to the expected damping rates from Fig. 5. 

Open with DEXTER  
In the text 
Fig. 7 Longterm evolution of planets with different inclinations, eccentricities, and masses in discs. The simulations are started from the equilibrium structures with fixed planets, where the planets are then released and allowed to move freely inside the disc. The top plot features the inclination of the planets, while the bottom plot shows the eccentricity of the planets. In the beginning the changes of eccentricity and inclination match the ones displayed in Figs. 4 and 5. 

Open with DEXTER  
In the text 
Fig. 8 Top: di/dt for a 5 M_{Jup} planet with different eccentricities and inclinations. The values of di/dt have been determined with the formula given in Sect. 3.5. Bottom: same, but for 10 M_{Jup}. We note the different colour coding as the change is dependent on the planetary mass. 

Open with DEXTER  
In the text 
Fig. 9 Top: de/dt for a 5 M_{Jup} planet with different eccentricities and inclinations. The values of de/dt have been determined with the formula given in Sect. 3.5. Bottom: same, but for 10 M_{Jup}. The white line in the figure indicates the transition between eccentricity increase and eccentricity damping. Below the white line, the eccentricity increases, above the line eccentricity decreases. We note the different colour coding as the change is dependent on the planetary mass. 

Open with DEXTER  
In the text 
Fig. 10 Evolution of e and i of the 10 M_{Jup} planet (shown in Fig. 7) with i_{0} = 75° and e_{0} = 0.4 in the ei plane (black line). The background is the extended formula of the fit for de/dt and the white line marks the transition between eccentricity increase and damping as in Fig. 8. The blue lines indicate calculated trajectories of 10 M_{Jup} planets from the fitting formulae. 

Open with DEXTER  
In the text 
Fig. A.1 F_{N} for 5 M_{Jup} planets with i = 75° on circular orbits. The planet’s feature different time steps, as indicated in the plot. F_{N} has been averaged over one running orbit. 

Open with DEXTER  
In the text 
Fig. A.2 Change of inclination di/dt for 5 M_{Jup} and 10 M_{Jup} planets in circular orbits for two different smoothing length for the planetary potential. The cubic potential is used for all displayed simulations, unless stated otherwise. 

Open with DEXTER  
In the text 
Fig. A.3 F_{N} for a 10 M_{Jup} planet with i = 3° in a e_{0} = 0.4 orbit. F_{N} has been averaged over one running orbit. 

Open with DEXTER  
In the text 