Free Access
Issue
A&A
Volume 555, July 2013
Article Number A124
Number of page(s) 13
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201220310
Published online 11 July 2013

© 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 spin-orbit 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 Kepler-30 system (Sanchis-Ojeda et al. 2012) is even flatter, and confirms this view. However, exo-planets with strong spin-orbit 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 equator1, 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 planet-planet 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 planet-planet interactions during migration in the protoplanetary disc (Thommes & Lissauer 2003; Libert & Tsiganis 2009, 2011a,b). During the gas-driven migration, the system can enter an inclination-type resonance or the resonant configuration becomes unstable as the resonance excites the eccentricities of the planets and planet-planet 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 low-mass planet embedded in a disc is exponentially damped by planet-disc interactions for any non-vanishing 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 high-mass planets, the linear regime is no longer valid. Marzari & Nelson (2009) considered Jovian-type 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 time-scale of the order of 103 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°.

Planet-disc 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 planet-disc interaction under some conditions (Goldreich & Sari 2003; Sari & Goldreich 2004; Moorhead & Adams 2008). They estimate that eccentric Lindblad resonances can cause eccentricity growth for gap-forming 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 very-high-mass 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  MJup) planets. Xiang-Gruess & Papaloizou (2013) have recently studied the interactions between Jupiter-mass 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 grid-based code.

In this paper, we investigate planet-disc interactions for planets above 1 MJup, considering different inclination and eccentricity values. Our analysis also aims at deriving a formula for the change of eccentricity and inclination due to planet-disc interactions, in order to study the long-term evolution of systems with massive planets. Indeed, long-term evolution studies of planetary systems cannot be done with hydrodynamical simulations, as the computation time is too long, and N-Body 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 three-dimensional (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 single-planet systems are shown in Sect. 4.

2. Physical modelling

The protoplanetary disc is modelled as a 3D, non-self-gravitating gas whose motion is described by the Navier-Stokes equations. We use the code Nirvana (Ziegler & Yorke 1997; Kley et al. 2001), which uses the FARGO-algorithm (Masset 2000) and was described in our earlier work on planets on inclined orbits (Bitsch & Kley 2011). We note that the use of the FARGO-algorithm 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 high-mass planets that open a gap inside a disc, where the effects of heating and cooling of the disc are much less important than for low-mass 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 MP is the planetary mass, and d =  | r − rP | 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 cubic-potential has been suggested (Klahr & Kley 2006; Kley et al. 2009) (2)The construction of the planetary potential is such that for distances larger than rsm the potential matches the correct 1/r potential. Inside this radius (d < rsm) it is smoothed by a cubic polynomial. This potential has the advantage of exactness outside the specified distance rsm, while being finite inside.

For 1 MJup and 5 MJup we use the cubic potential with rsm = 0.8RH. For the 10 MJup planet, we use the ϵsm-potential with rsm = 0.8RH, with the Hill radius RH given by (3)where ap 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 MJup planet. For the torque acting on the planets, the consequences are minimal, as we use a torque cut-off 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 rmin = 0.2 to rmax = 4.2 in units of r0 = aJup = 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 e0 = 0.0, e0 = 0.2, or e0 = 0.4. We use 390 × 48 × 576 active cells for the simulations with 1 MJup and 260 × 32 × 384 active cells for 5 MJup and 10 MJup. 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 planet-star 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 [rmin,rmax] is Mdisc = 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 = cs/Ω.

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 (Fermi-type) function (5)which increases from 0 at the planet location (d = 0) to 1 outside d ≥ RH with a midpoint fb = 1/2 at d = bRH, i.e. the quantity b denotes the torque-cutoff radius in units of the Hill radius. This torque-cutoff 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 (semi-major 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 planet-disc interactions. The planets are put in fixed orbits with inclinations ranging from i0 = 1.0° to i0 = 75°, with a total of ten different inclinations. For each inclination we also adopt three different eccentricities, which are e0 = 0.0, e0 = 0.2 and e0 = 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 MJup) 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 MJup planets on different inclinations.

thumbnail Fig. 1

Surface density for disc simulations with 10 MJup 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 time-scale 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 − z-plane for the disc’s density for 10 MJup 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.

thumbnail Fig. 2

Density (in g/cm3) of a r − θ-slice through the disc at the azimuth of an embedded 10 MJup planet on a fixed circular inclined orbit with i0 = 1.0° (top), i0 = 20° (middle), and i0 = 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 MJup 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.

thumbnail Fig. 3

Eccentricity (top) and inclination (bottom) of the disc with a 10 MJup planet influencing the disc structure after 400 planetary orbits.

Open with DEXTER

The inclination of the disc for the i0 = 1° planets is greater mostly around the planet’s location (at r = 1.0aJup) 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 non-inclined.

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.0aJup) compared to the case of low inclinations. However, the outer parts of the disc show a non-zero inclination (which has a tendency to be larger for larger planetary eccentricities), which was not visible for the low-inclination planets. Additionally, they show a small peak of inclination at r ≈ 1.25aJup.

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 planet-disc 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 MJup, 5 MJup, and 10 MJup on orbits with different inclinations and eccentricities. The change of de/dt has been studied in the past for coplanar low-mass planets (Cresswell et al. 2007; Bitsch & Kley 2010) and for high-mass planets (Papaloizou et al. 2001; Kley & Dirksen 2006).

thumbnail Fig. 4

Change of eccentricity de/dt for planets with 1 MJup, 5 MJup, and 10 MJup with different eccentricities. Points are results from numerical simulations, while lines indicate the fitting of the data. The 1 MJup planets have been evolved for 200 planetary orbits, the 5 MJup and 10 MJup 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 (i0 < 10°) the damping of eccentricity is stronger than for larger inclinations in the case of 1 MJup. The maximum damping rate is also dependent on the initial eccentricity e0, where a larger e0 provides a faster damping. The damping of eccentricity is reduced significantly for larger inclinations i0 > 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 (i0 < 10°) and low eccentricities (e0 < 0.2), the 5 MJup 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 (e0 = 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 i0 > 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 MJup), we observe an eccentricity increase for low planetary inclinations for all non-zero eccentricities. But for larger inclinations (i0 > 15°), the eccentricity is damped again. The largest value of damping is at i0 ≈ 30° − 50°, depending on the planet’s eccentricity and is then reduced for higher inclinations again, following the trend described for the 5 MJup 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 MP > 20 MJup, while our simulations indicate it clearly already for MP > 5 MJup (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 high-mass 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.5aP and a large planetary eccentricity (the i0 = 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 MJup the inclination is damped for all initial inclinations. For increasing inclinations with i0 < 15° (smaller for increasing eccentricity), the damping of inclination increases. This increase is nearly linear, as has been shown for low-mass planets in theory (Tanaka & Ward 2004) and in numerical simulations (Cresswell et al. 2007; Bitsch & Kley 2011). The rates of inclination damping for zero-eccentricity planets are comparable to those stated in Xiang-Gruess & Papaloizou (2013).

thumbnail Fig. 5

Change of inclination di/dt for planets with 1 MJup, 5 MJup, and 10 MJup with different eccentricities. di/dt is in degrees per orbit at the planet’s location rP = 1.0aJup. Points are results from numerical simulations, while lines indicate the fitting of the data. The 1 MJup planets have been evolved for 200 planetary orbits, the 5 MJup and 10 MJup 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 i0 > 15°, the damping rate of the inclination is a decreasing function of inclination ; this is consistent with the planet-disc interaction being weaker when the planet spends more time farther from the midplane.

For 5 MJup the damping of inclination is almost the same as for the 1 MJup planet, but with a maximum at i0 ≈ 20°. However, there is a significant difference for high inclined and low eccentric planets: the inclination is not damped if i0 > 50°, but it increases for e0 < 0.1. This behaviour will be discussed in Sect. 3.4.

The 10 MJup planet shows the same general behaviour as the 5 MJup planet, but the inclination increase already sets in at , depending on e0. Still, no inclination increase is observed in the high eccentricity simulations (e0 = 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 high-mass planets due to interactions with the disc has been studied in Lubow & Ogilvie (2001). They state that the 1:3 mean-motion resonance also acts to increase inclination. This resonance is at rr = 2.08rP, which clearly is not inside an open gap in the case of i0 = 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 i0 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. Short-term 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 MJup and 10 MJup with an inclination of i0 = 40° and i0 = 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).

thumbnail Fig. 6

Evolution of inclination of 5 MJup and 10 MJup planets with an initial inclination of i0 = 40.0° (top) and i0 = 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 low-mass 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 vice-versa is much larger.

3.4.2. Long-term evolution

The long-term 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.

thumbnail Fig. 7

Long-term 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 MJup, e0 = 0.0, i0 = 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 MJup), 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 (e0 = 0.2 or e0 = 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 (i0 <   40° in the restricted problem of Kozai 1962). We therefore also display a planet with i0 = 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 e0 = 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 low-mass 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 MJup(e = 0.0): the eccentricity value reached during the second cycle of inclination increase (at 450 orbits) is higher for the 5 MJup planet, as expected.

The effect of Kozai oscillations between a disc and planet was also stated in Xiang-Gruess & Papaloizou (2013) however, Xiang-Gruess & Papaloizou (2013) were not able to resolve a full Kozai cycle, probably because their mass-ratio between planet and disc is smaller than in our case. This shows that for i0 > 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 long-term evolution of the orbital parameters. This phenomenon can be of crucial importance for the study of the fate of planets scattered on high-inclination 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 MP, the eccentricity eP, and inclination iP. The inclination iP 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 i0 < 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 iP for small iP, and a decreasing power law of iP for large iP. 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 iP, , and for large iP, . 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 Mdisc/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 iP tends to zero. A pure increasing power law of iP is inappropriate here. The damping function will be given by (12)where iD is a small inclination so that for iP ≈ 0, . We are using degrees in Eq. (12), where is the planet mass in Jupiter masses. For small eccentricities, it is expected that eP/(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 high-mass planets. The damping and excitation of eP 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 i0 = 0° high mass-planets. 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 MP < 5 MJup 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 iP centred on 0°, so that this term is not affected for small iP 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 iP, di/dt should be close to linear in iP, 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 ci 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 eP and iP. In Fig. 8 the di/dt for different inclinations and eccentricities for 5 MJup and 10 MJup according to the formulae is presented. In Fig. 9 the de/dt for the same two planetary masses is plotted.

thumbnail Fig. 8

Top: di/dt for a 5 MJup 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 MJup. 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 (iP ≈ 15°). The inclination damping rate indicates that planets that are scattered during the gas disc phase in orbits with moderate inclination (iP < 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 high-mass cases, but the increase of eccentricity declines with increasing eccentricity and inclination. Additionally, the threshold of eP and iP 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.

thumbnail Fig. 9

Top: de/dt for a 5 MJup 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 MJup. 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 single-planet 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 Kozai-oscillation 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 MJup planet with i0 = 75° and e0 = 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 Kozai-oscillations 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 non-zero eccentricity. This non-zero 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 104 years. A typical damping rate of di/dt = 0.01 deg/orbit indicates that the planet will lose ≈8.5° of inclination in 104 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 non-zero eccentricity. This is illustrated by the blue lines in Fig. 10 that correspond to calculated trajectories of 10 MJup 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 high-mass planets (MP ≥ 1 MJup) 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 MJup the damping rate of e and i is highest for only very small inclinations (i0 ≈ 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 low-inclination planets with a sufficient mass (MP > 4 − 5 MJup). 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 (MP ≈ Mdisc, in our case for MP > 5 MJup) on high initial inclinations (i0 > 40°). In the long-term evolution of the planet, eccentricity can increase, while inclination is damped and vice-versa. The planet undergoes a Kozai-cycle 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.

thumbnail Fig. 10

Evolution of e and i of the 10 MJup planet (shown in Fig. 7) with i0 = 75° and e0 = 0.4 in the e-i 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 MJup planets from the fitting formulae.

Open with DEXTER

In Sect. 3.5 we provided formulae for di/dt and de/dt for high-mass planets, which we fitted to the numerical hydrodynamical simulations. The formulae can now be used to calculate the long-term evolution of planetary systems during the gas phase of the disc with N-Body 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 long-term evolution of planets is hard to predict).

In the end, the planet’s inclination will be damped to zero. Low-mass planets (MP < 4 − 5 MJup) will end up in circular orbits in the midplane of the disc, while higher mass planets (MP > 5 MJup) 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 long-term evolution of multi-body systems will be studied in a future paper.


1

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

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 FN (which has been averaged over one orbit) of planets with i = 75° on circular orbits. The two simulations shown feature different time-step 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 FN 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 ≈ aP.

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 rsm. A smaller smoothing length rsm 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 MJup planet.

In Fig. A.2 we present the inclination damping for 5 MJup planets in circular orbits for two different smoothing lengths, rsm = 0.8 and rsm = 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 i0 = 55° case with a smoothing length of rsm = 0.5. This has also been observed for the 10 MJup 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 MJup planet. For i0 = 75° the inclination seems to increase for models of planets in fixed orbits for both 5 MJup and 10 MJup. We therefore conclude that the general trend is conserved regardless of the chosen planetary potential and smoothing length.

thumbnail Fig. A.1

FN for 5 MJup planets with i = 75° on circular orbits. The planet’s feature different time steps, as indicated in the plot. FN has been averaged over one running orbit.

Open with DEXTER

thumbnail Fig. A.2

Change of inclination di/dt for 5 MJup and 10 MJup 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

thumbnail Fig. A.3

FN for a 10 MJup planet with i = 3° in a e0 = 0.4 orbit. FN 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 FN has been averaged over one running orbit. Keep in mind that FN has been averaged over 40 orbits to determine the change of inclination in the end. The simulations feature a 10 MJup 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 FN 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 mass-weighted 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 semi-major 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 Lspec. = r × v is the specific angular momentum of each grid cell. To get an estimate of the eccentricity of the disc, we make a mass-weighted average of the eccentricity of each grid cell (averaged in azimuthal and polar coordinates) in order to get edisc(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 Lspec. 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 mass-averaged specific angular momentum (averaged in polar and azimuthal direction) (B.9)where the subscript c denotes the grid cell number, and mc the corresponding mass of the grid cell. Now we can compute the angle between Lav. and the z-axis, 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 high-mass planets, as it is given by Ag 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 ed is the disc’s eccentricity and rP the planetary distance to star. The integral basically gives the disc mass, which is comparable to the planet’s mass, but as ed is much smaller than eP (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 ed = 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 = 2ePγ(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

thumbnail Fig. 1

Surface density for disc simulations with 10 MJup 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
thumbnail Fig. 2

Density (in g/cm3) of a r − θ-slice through the disc at the azimuth of an embedded 10 MJup planet on a fixed circular inclined orbit with i0 = 1.0° (top), i0 = 20° (middle), and i0 = 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
thumbnail Fig. 3

Eccentricity (top) and inclination (bottom) of the disc with a 10 MJup planet influencing the disc structure after 400 planetary orbits.

Open with DEXTER
In the text
thumbnail Fig. 4

Change of eccentricity de/dt for planets with 1 MJup, 5 MJup, and 10 MJup with different eccentricities. Points are results from numerical simulations, while lines indicate the fitting of the data. The 1 MJup planets have been evolved for 200 planetary orbits, the 5 MJup and 10 MJup 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
thumbnail Fig. 5

Change of inclination di/dt for planets with 1 MJup, 5 MJup, and 10 MJup with different eccentricities. di/dt is in degrees per orbit at the planet’s location rP = 1.0aJup. Points are results from numerical simulations, while lines indicate the fitting of the data. The 1 MJup planets have been evolved for 200 planetary orbits, the 5 MJup and 10 MJup 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
thumbnail Fig. 6

Evolution of inclination of 5 MJup and 10 MJup planets with an initial inclination of i0 = 40.0° (top) and i0 = 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
thumbnail Fig. 7

Long-term 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
thumbnail Fig. 8

Top: di/dt for a 5 MJup 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 MJup. We note the different colour coding as the change is dependent on the planetary mass.

Open with DEXTER
In the text
thumbnail Fig. 9

Top: de/dt for a 5 MJup 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 MJup. 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
thumbnail Fig. 10

Evolution of e and i of the 10 MJup planet (shown in Fig. 7) with i0 = 75° and e0 = 0.4 in the e-i 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 MJup planets from the fitting formulae.

Open with DEXTER
In the text
thumbnail Fig. A.1

FN for 5 MJup planets with i = 75° on circular orbits. The planet’s feature different time steps, as indicated in the plot. FN has been averaged over one running orbit.

Open with DEXTER
In the text
thumbnail Fig. A.2

Change of inclination di/dt for 5 MJup and 10 MJup 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
thumbnail Fig. A.3

FN for a 10 MJup planet with i = 3° in a e0 = 0.4 orbit. FN has been averaged over one running orbit.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.