EDP Sciences
Free Access
Issue
A&A
Volume 553, May 2013
Article Number A39
Number of page(s) 15
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201220482
Published online 29 April 2013

© ESO, 2013

1. Introduction

Today, nearly 100 multi-planet systems have been reported, of which roughly 1/3 possess at least one “moderate close-in planet”, that is, a planet with a semi-major axis between 0.1 and 0.5 AU1. Planets in this range are supposed to undergo significant tidal interactions, resulting in slowly modified spins and orbits. However, for the typically assumed dissipation rates for gaseous planets, the spin of moderate close-in planets reaches an equilibrium state in only a few million years, while the orbital evolution can last for the entire age of the system (Gyr timescale).

Among two-planet systems, there is a special class whose semi-major ratio a1/a2 is lower than 0.1, the so-called “hierarchical systems”. This class counts at least 20 members (1/5 of all multi-planet systems), and usually at least one of the planets’ orbits is highly eccentric1. During the formation process, the orbital eccentricities can increase through gravitational scattering (e.g. Nagasawa et al. 2008), which is simultaneously responsible for an increase of the orbits’ mutual inclination (e.g. Chatterjee et al. 2008). Evidently, smaller mass planets, that are as yet undetected, can exist in these systems, but for a semi-major axis larger than 0.1 AU the orbits usually present high values in the eccentricities, which reduces the stability areas for additional companions.

Hierarchical systems are particularly interesting from a dynamical point of view, because they can be stable with very eccentric and inclined orbits, which is an uncommon behavior. In particular, they are very interesting when the inner planet is sufficiently close to the star to undergo tidal interactions, since the final outcome of the evolution can be in a configuration that is totally different from the initial one (e.g. Wu & Murray 2003; Fabrycky & Tremaine 2007; Correia et al. 2011). Because the semi-major ratio is small, low-order mean motion resonances cannot occur, which allows us to perform analytical approximations such as averaging the orbits over the mean anomalies. In addition, tidal effects usually act over very long timescales and therefore approximate theories also allow one to accelerate the numerical simulations and to explore the parameter space much more rapidly.

Secular perturbation theories based on series expansions have been used for hierarchical systems. For low eccentricity values, the expansion of the perturbation in eccentricity series is very efficient (e.g. Wu & Goldreich 2002), but this method is not appropriate for orbits that become very eccentric. An expansion in the ratio of the semi-major axis a1/a2 is then preferred, because exact expressions can be computed for the secular system (e.g. Laskar & Boué 2010). The development to the second order in a1/a2, called the quadrupole approximation, was used by Lidov (1961, 1962) and Kozai (1962) for the restricted inner problem (the outer orbit is unperturbed). In this case, the conservation of the normal component of the angular momentum enables the inner orbit to periodically exchange its eccentricity with inclination (the so-called Lidov-Kozai mechanism). For planar problems, the series expansions in a1/a2 should be conducted to the octupole order (e.g. Marchal 1990; Ford et al. 2000; Laskar & Boué 2010), because the quadrupole approximation fails to reproduce the eccentricity oscillations (e.g. Lee & Peale 2003). Since we do not have any restrictions for the eccentricities or for the mutual inclination, we need to expand the gravitational potential in a1/a2 to the octupole order.

The ultimate stage for tidal evolution is the spin synchronization and orbital circularization (e.g. Correia 2009). Indeed, the observed mean eccentricity for planets and binary stars with a1 < 0.1 AU is close to zero within the observational limitations (e.g. Pont et al. 2011). Although tidal effects modify the spin on a much shorter timescale than they modify the orbit, synchronous rotation can only occur when the eccentricity is very close to zero: the rotation rate tends to be locked with the orbital speed at the periapsis, because tidal effects are stronger when the two bodies are closer to each other. In addition, if there is a companion body, the eccentricity oscillates (e.g. Mardling 2007; Laskar et al. 2012), and the rotation rate of the planet shows variations that follow the eccentricity. This is exactly what is observed for Mercury, whose average orbital eccentricity is around 0.2, and its rotation is captured in a 3/2 spin-orbit resonance (Correia & Laskar 2004, 2012). Therefore, spin and orbital evolution cannot be dissociated, and some unexpected behavior can be observed, such as a secular increase of the eccentricity (e.g. Correia et al. 2011, 2012).

In this paper we intend to intensify the study of hierarchical two-planet systems, in which the inner orbit undergoes tidal dissipation. We present here another counterintuitive behavior, the inclination damping, which is also a consequence of the above-mentioned eccentricity pumping when the two orbits are not coplanar. We provide a simple averaged model for the orbital and spin evolution of an exoplanet with a companion (Sect. 2), and apply it to hierarchical two-planet systems in a conservative frame (Sect. 3) and including tidal effects (Sect. 4). We finally derive some conclusions (Sect. 5), and provide an explanation for the inclination damping (Appendix A).

2. Model

We considered a system consisting of a central star of mass m0, an inner planet of mass m1, and an outer companion of mass m2. We used Jacobi canonical coordinates, with r1 being the position of m1 relative to m0, and r2 the position of m2 relative to the center of mass of m1 and m0. We assumed that the system is hierarchical, thus | r1| ≪ |r2 |. For simplicity, we express all the angles in the invariable plane of the system, i.e., the plane perpendicular to the total angular momentum (1)where L is the rotational angular momentum of the inner planet, and Gi the orbital angular momentum of each body.

The inner planet is considered an oblate ellipsoid with gravity field coefficients given by J2, rotating about the axis of maximal inertia (gyroscopic approximation), with rotation rate ω, such that (e.g. Lambeck 1988) (2) is the gravitational constant, is the radius of the planet, and k2 is the second Love number for potential (pertaining to a perfectly fluid body). We furthermore assumed that the obliquity of the planet to the orbital plane is zero (ε = 0°), that is, L and G1 are aligned. Therefore, the angle between the two orbital planes, i.e., the mutual inclination I, satisfies the relation (3)

2.1. Conservative motion

Because we are interested in the secular behavior of the system, we averaged the equations of motion over the mean anomalies of both orbits. In the invariable plane, the averaged potential, quadrupole-level for the spin (e.g. Correia & Laskar 2010a), octupole-level for the orbits (e.g. Ford et al. 2000; Laskar & Boué 2010), and with general relativity corrections (e.g. Touma et al. 2009) is given by (4)where and ϕ is the angle between the directions of the periastrons: (9)ai is the semi-major axis (that can also be expressed using the mean motion ni), ei is the eccentricity, and ωi is the argument of the periastron. We also have β1 = m0m1/(m0 + m1), β2 = (m0 + m1)m2/(m0 + m1 + m2), , and , where ξ is the normalized moment of inertia.

The contributions to the orbits are easily obtained using the Lagrange planetary equations (e.g. Murray & Dermott 1999): (10)In addition, since the variations in e1, e2 and I are related by the conservation of the total angular momentum (Eq. (3)), we have (11)As we assumed that , if we neglect first-order terms in L/G1, we simply have (i ≠ j = 1,2): (12)Thus, and where , the constant frequencies and When I = 0°, we have ϕ = ϕ1 = ϕ2 = ω2 − ω1. Note also that the longitude of the node does not appear in the equations of motion (Eqs. (13)−(16)) because we used the invariable plane as the reference plane (Eq. (4)), and thus ΔΩ = 180°.

Table 1

Hierarchical two-planet non-resonant systems with 0.1 < a1 < 0.5 AU and a1/a2 < 0.1.

2.2. Tidal effects

In our model, we additionally considered tidal dissipation raised by the central star on the inner planet. The dissipation of the mechanical energy of tides in the planet’s interior is responsible for a time delay between the initial perturbation and the maximal deformation. Because the rheology of planets is poorly known, the exact dependence of on the tidal frequency is unknown. Several models exist (for a review see Correia et al. 2003; Efroimsky & Williams 2009), but for simplicity we adopted a model with constant , which can be made linear (Singer 1968; Mignard 1979). The contributions to the equations of motion are given by (e.g. Correia 2009; Correia et al. 2011) where and We neglected the effect of tides over the argument of the periastron, as well as the flattening of the central star. Their effect is only to add a small supplementary frequency to , similar to the contributions from the general relativity (for a complete model see Correia et al. 2011). Expression (28) for the inclination is zero, because we assumed the obliquity of the planet to be zero (ε = 0°).

Under the effect of tides alone, the equilibrium rotation rate, obtained when , is attained for (Eq. (25)) (37)Usually, , so tidal effects modify the rotation rate much faster than the orbit. It is thus tempting to replace the equilibrium rotation in expressions (26) and (27). With this simplification, one always obtains negative contributions for ȧ1 and ė1 (Correia 2009), with (40)Thus, the semi-major axis and the eccentricity can only decrease until the orbit of the planet becomes circular. However, planet-planet interactions can produce eccentricity oscillations with a period shorter than, or comparable to, the damping timescale of the spin. In that case, expression (37) is not satisfied and multi-planetary systems may show non-intuitive eccentricity evolutions, such as eccentricity pumping of the inner orbit (Correia et al. 2012).

3. Application to exoplanets

In the following sections we apply the model described in Sect. 2 to different configurations of hierarchical two-planet systems. To observe the damping effect of the mutual inclination, the spin of the planet must be fully damped, but not its orbit, that is, (see Appendix A). In addition, the damping timescale of the spin should be of the same order as the period of the eccentricity oscillations, (Eq. (A.18)). This is valid for gaseous planets roughly within 0.1 < a1 < 0.5 AU, which we call “moderate close-in planets”. In Table 1 we list all hierarchical systems known to date whose inner orbit satisfies the above condition. We focus our analysis on the HD 74156 system, but all the main results are easily extended to the remaining planetary systems.

3.1. Radius of close-in exoplanets

According to expression (A.18), the ν1 frequency (Eq. (18)) is a key parameter for the observation of the eccentricity pumping of the inner orbit and consequent damping of the inclination (by means of gx). The minimum masses and the semi-major axis are relatively well determined from the observations, so the largest incertitudes in ν1 come from the Love number k2, and particularly from the radius of the planet, which appears as a power of 5. Therefore, a correct estimate of the planetary radius is necessary to observe some effect on the inclination.

Since the radius of the planet is correlated with its mass, one solution is to adopt a constant value for the density, , and then compute the radius simply as . However, by applying this strategy to the two largest planets of the solar system, Jupiter and Saturn, we immediately see that it can give very distinct results. The density of a planet depends on many factors, such as the age of the system, the initial composition of the accretion disk, or where the planet formed in the disk. Any theoretical estimation of the radius is then subject to large incertitude, and only direct observations can give reliable values.

thumbnail Fig. 1

Radius versus the mass of the planet. We plot all known close-in planets in the range 0.1 < a1 < 0.5 AU, for which the radius was determined by the transiting method. We observe that the radius decreases with the mass in a relatively regular way.

Open with DEXTER

We used an empirical expression based on the previously observed radius of close-in planets in the range 0.1 < a1 < 0.5 AU. We found ten planets in this range1 whose radius were determined by the transiting method (Fig. 1). We observe that the radius decreases with the mass in a relatively regular way, therefore we performed a linear regression of the observational data: (41)This expression also agrees well with the solar system data, giving 1.1   RJup for Jupiter and 1.0   RSat for Saturn.

3.2. Initial conditions uncertainty

Assuming that the observational values of the minimum masses, semi-major axis, and eccentricities of the planetary systems listed in Table 1 are relatively well determined, we can use them as a starting point to study these systems.

thumbnail Fig. 2

Possible secular trajectories for the HD 74156 system seen in the (ω1,I) plane (left), and in the (ω1,e1) plane (right). We show the trajectories using the quadrupolar approximation (top), corresponding to the level curves of constant energy, and the octupolar approximation (bottom). The dashed black curves in the (ω1,I) plane correspond to the separatrix between the circulation and libration zones of ω1. The colors are preserved in all pictures, each one corresponding to a given value of the total angular momentum of the system (determined by different values of the initial mutual inclination). All trajectories are compatible to the present knowledge for this system (Table 2).

Open with DEXTER

A striking observation is that the eccentricity of at least one of the planets can be very high. Because hierarchical systems exhibit high eccentricity values, the mutual inclinations can also be very high (e.g. Chatterjee et al. 2008). However, currently their mutual inclinations are unknown, not only because we are unable to determine the inclination with respect to the plane of the sky, Ii (and hence the true masses), but also because we are unable to determine the longitude of the nodes in the plane of the sky, Ωi: (42)The only partial exception is the system HD 38529, for which I2 ≈ 48°, estimated using astrometric measurements from the Hubble Space Telescope (Benedict et al. 2010). The planetary masses m1 and m2 given in Table 1 correspond to the minimum masses (assuming I1 = I2 = 90°), except for HD 38529, where the masses were estimated using I1 = I2 = 48°.

According to expression (3), without knowing the mutual inclination of these systems it is impossible to determine the total angular momentum H, which is a constant of the motion. Even if we assume that the system is coplanar and prograde (cosI = 1), the total angular momentum is undetermined because the true planetary masses appearing in the expressions of G1 and G2 are also unknown (except for HD 38529). As a consequence, the present dynamics of these systems can be considerably different, depending on the true H value.

Moreover, although the argument of the periastron of these planets is known in the plane of the sky (angle in Table 1), the arguments of the periastron in the invariable plane of the system (ω1 and ω2) are also unknown because they depend on Ii and Ωi (e.g. Giuppone et al. 2012): (43)where is the argument of the periastron in the plane of the sky. For ω2 we have an identical relation, the only difference is that ω1 is measured from the ascending node, while ω2 is measured from the descending node.

Therefore, apart from the semi-major axes, eccentricities and minimum masses, currently there are few constraints on the remaining orbital parameters of hierarchical systems. Because we are only concerned with the tidal evolution of the mutual inclination, we adopted the masses listed in Table 1 for the numerical simulations. The only free parameters are then the mutual inclination and the arguments of the periastron, whose uncertainty are related to the lack of knowledge of the longitude of the node in the plane of the sky (Ω2 − Ω1).

Using the quadrupolar approximation for gravitational interactions, the potential energy (Eq. (4)) is independent of the ω2, and therefore e2 is constant (Eq. (10)). The dynamics of the system is then fully described by the couples (I,ω1) or (e1,ω1). Adopting the minimum masses (Table 1), we show in Fig. 2a, b some possibilities for the HD 74156 system for different mutual inclinations (corresponding to different H values). For nearly coplanar systems (I < 20°), only small variations are observed for e1 and I. However, as the inclination increases, the dynamics of the system is considerably perturbed by the presence of Lidov-Kozai cycles (Lidov 1961, 1962; Kozai 1962). In this regime, we can observe significant exchanges between the inclination and the eccentricity of the inner orbit. In some cases the eccentricity can reach values much higher than today, and thus enhance the tidal dissipation.

Figure 2a, b is different from the standard Lidov-Kozai diagrams that show level curves of the quadrupole Hamiltonian at fixed values G1   cosI, because for HD 74156 the initial eccentricity of the inner planet is already fixed at e1 = 0.64 (Table 1). Instead, we followed the procedure in Giuppone et al. (2012) and show the trajectories for different values of the initial mutual inclination I, i.e., we varied the total angular momentum H. This explains why librating orbits in the Lidov-Kozai regime do not encircle a Lidov-Kozai equilibrium point (which occurs at the current eccentricity value).

The impact of the initial uncertainty on the ω1 value is shown in Fig. 3. Depending on this value, the observed eccentricity can correspond to a maximum or to a minimum for an identical total angular momentum H (Fig. 3a). Moreover, two trajectories may be in circulation or in libration (Fig. 3b). It is therefore very important to completely explore the phase space of the initial conditions (I,ω1) to capture the global dynamics of hierarchical two-planet systems.

3.3. Octupole contribution

So far, we restricted our analysis to the quadrupolar gravitational interactions, because they are mainly responsible for the inclination variations. However, for the hierarchical systems listed in Table 1, the range of semi-major ratios is 0.03 < a1/a2 < 0.1, meaning that octupolar interactions cannot be neglected. Indeed, although the impact of octupolar terms on the eccentricity variations is weaker than that of the quadrupole terms, octupolar interactions are strong enough to produce secular drifts when combined with tidal effects (Correia et al. 2012).

In the planar prograde case (cosI = 1), the potential energy (Eq. (4)) only depends on ϕ = ω2 − ω1, and therefore I is constant (Eq. (12)). The dynamics of the system is then fully described by the couples (ϕ,e1) or (ϕ,e2). Adopting the minimum masses, the angle ϕ in the plane of the sky listed in Table 1 corresponds to the angle (ω2 − ω1 + 180°) in the invariant plane of the system, because for I1 = I2 = 90° we have and (Eq. (43)). In Fig. 4a we show the expected eccentricity variations for the HD 74156 system in this unique situation for which the system is fully characterized.

Table 2

Stability analysis of the HD 74156 system for different sets (ω1,I) of initial conditions (Fig. 2).

thumbnail Fig. 3

Possible secular trajectories for the HD 74156 system seen in the (ω1,e1) plane for two different values of the initial mutual inclination, I0 = 40° (top), and I0 = 50° (bottom). We show the trajectories using the quadrupolar approximation, corresponding to the level curves of constant energy. Each one corresponds to an initial value of the argument of the periastron ω1, ranging from 0° (black) to 90° (red) with a step of 10°. The dashed black curves correspond to the observed eccentricity of the planet, that is, it gives the initial condition for ω1. All trajectories are compatible with the current knowledge for this system (Table 1).

Open with DEXTER

thumbnail Fig. 4

Possible secular trajectories for the HD 74156 system seen in the (ϕ,e1) plane for I = 0° (top) and in the (ω1,e1) plane for I = 35° (bottom). All trajectories are compatible with the current knowledge for this system (Table 1). We show the trajectories using the octupolar approximation (blue) and direct numerical simulations (red). In a) the blue path corresponds to the level curves of constant energy for coplanar orbits. The dot marks the present position of the planet. In b) we additionally show the trajectories using the quadrupolar approximation (green), which corresponds to the level curves of constant energy.

Open with DEXTER

With increasing mutual inclination, we are left with four free parameters (e1,ω1,e2,ω2), and it becomes impossible to capture the dynamics of the systems in a two-dimensional plot. Nevertheless, we can perform numerical simulations of the equations of motion (Eqs. (13)−(16)) to understand how the octupolar terms modify the quadrupolar approximation (Fig. 2). We observe that the main effect is to add some diffusion around the quadrupolar trajectories. The diffusion is more pronounced for orbits in circulation around the separatrix (35° < I < 145°). The stability of the orbits can be measured with a frequency analysis (Laskar 1990, 1993). We determined the precession frequency g and g′ of the argument of the pericenter ω1 over two consecutive time intervals of length T = 5 Myr. In Table 2 we compute the difference D = |g − g′|/g, which is a measure of the chaotic diffusion of the trajectories (Correia et al. 2005; Couetdic et al. 2010). It should be close to zero for a regular solution, and values with D > 10-6 correspond to chaotic motion. This is observed for all trajectories in circulation close to the separatrix.

In Fig. 4b we plot simultaneously the eccentricity evolution obtained with the two approximations for I0 = 35°, and compare it with direct numerical simulations. We conclude that 1) the quadrupolar approximation correctly describes the average dynamics in inclined hierarchical systems; 2) the octupolar approximation is essential to derive a more realistic behavior and obtain results similar to direct numerical simulations.

4. Tidal evolution

We now include the effect of tides described in Sect. 2.2 to the conservative equations of the motion (Eqs. (13)−(16)), and perform some numerical simulations. In all simulations we adopt for the innermost planet , k2 = 1/2, and a dissipation time lag  s. For HD 74156 this dissipation is equivalent to (Eq. (31)), which is comparable to the value estimated for Saturn (Lainey et al. 2012). In addition, we always set ω2 = 180° and  day. The impact of the initial ω2 value can be obtained by adjusting a different value for ω1. Similarly, the initial rotation rate is not a critical initial parameter, since tidal effects quickly bring the rotation near to the equilibrium value (Eq. (37)).

4.1. Effect of the spin

In Fig. 5 we show some examples for the evolution of the HD 74156 planetary system, in three different situations. The radius of the inner planet is estimated to be (Eq. (41)), and we initially assume I0 = 40° and ω1 = 0°.

In a first experiment, we only consider tidal effects on the orbit (Eqs. (38), (39)), as it is often done in previous studies. That is, since the rotation of close-in planets evolves very fast, we assume that the spin is locked in its equilibrium position (Eq. (37)). We observe that the eccentricity and the semi-major axis of the inner orbit slowly decrease, while the mutual inclination and the eccentricity of the outer orbit only oscillate around a constant mean value (Fig. 5a). Therefore, we conclude that in this case the only effect of tides is to circularize the inner orbit in a timescale longer than the age of the system, a well-know result in the literature (e.g. Correia & Laskar 2010b).

In a second experiment, we neglect the effects on the orbit and we only consider the effect on the rotation (Eq. (25)). This situation corresponds to the opposite of the previous one, and it is not realistic, but it allows to highlight the importance of not neglecting the rotation rate evolution. Indeed, although there is no direct dissipative contribution to the eccentricity or to the inclination, we observe that these two parameters undergo significant variations (Fig. 5b), the eccentricity of the inner planet rising almost up to 1. The pumping effect on the eccentricity due to the spin excitation was reported in the planar case by Correia et al. (2012), for which only octupolar terms are important. In the non-planar case, the pumping effect it is even more pronounced, since quadrupole order terms additionally contribute. In addition, because the angular momentum is mainly exchanged between the inner planet eccentricity and the inclination, while the first increases, the second decreases. In Appendix A we provide the full explanation for this effect in the frame of the quadrupolar approximation. We also observe that the eccentricity of the outer planet is slightly damped during this process, because of the octupole order interactions (Correia et al. 2012).

thumbnail Fig. 5

Long-term evolution of the HD 74156 system for different tidal models with I0 = 40°: only dissipation in the orbit is considered a); only dissipation in the spin is considered b); full model c). We plot the mutual inclination I (top), the eccentricities e1 (blue) and e2 (green) (middle), and the semi-major axis a1 (bottom).

Open with DEXTER

thumbnail Fig. 6

Long-term evolution of the HD 74156 system with I0 = 40° for different values of the argument of the periastron, ω1 = 30°a); 50°b); and 70°c). We plot the mutual inclination I (top), the eccentricities e1 (blue) and e2 (green) (middle), and the semi-major axis a1 (bottom).

Open with DEXTER

thumbnail Fig. 7

Long-term evolution of the HD 74156 system with I0 = 50° for different values of the argument of the periastron, ω1 = 30°a); 60°b); and 90°c). As a function of time, we plot the mutual inclination I (top) and the eccentricities e1 (blue) and e2 (green) (middle). As a function of ω1, we plot again the mutual inclination, where the color of each dot becomes darker with time (bottom).

Open with DEXTER

Finally, since orbital and spin evolution cannot be dissociated, we integrate the full set of equations for the tidal evolution (Eqs. (25)−(28)) (Fig. 5c). We observe that the initial behavior of the system is identical to the situation without dissipation on the orbit (Fig. 5b), with a significant damping of the mutual inclination. However, as the eccentricity increases, the inner planet comes closer to the star at periastron, and tidal effects on the orbit become stronger. As a consequence, the semi-major axis decreases and the damping effect on the eccentricity (Eq. (27)) overrides the pumping drift (Eq. (A.18)). At this point, the inclination damping is less efficient, and it ceases when the pumping drift is completely gone. The system ultimately evolves into a circular orbit as usual, but in a considerable much shorter timescale and it is left with a final lower mutual inclination (~15°).

4.2. Effect of the argument of periastron

In Sect. 3.2 we have seen that the initial choice of the argument of the periastron of the inner planet, ω1, plays an important role on how the present eccentricity is changing (Fig. 3). Similarly, it also changes the initial trend of the inclination: for increasing eccentricity e1 the inclination decreases, and vice-versa. In the example from previous section (Fig. 5c), we used ω1 = 0°, that is, we assumed that the observed value of the eccentricity (e1 = 0.64) is a minimum, and the inclination I = 40° a maximum (Fig. 3a).

In order to test the impact of the initial argument of the periastron in the HD 74156 system, in Fig. 6 we plot its evolution for different initial ω1 values. We observe similar behavior as before for all situations, the only significant difference being the evolution timescale. Until ω1 < 50° this timescale increases, because the present eccentricity is no longer a minimum value. As a consequence, the average value of the eccentricity oscillations is shifted down, and tidal dissipation is less effective, since at the periastron the inner planet is farther from the star.

In the quadrupolar approximation, the behavior described above should be maintained up to ω1 < 90°, for which the observed eccentricity is a maximum and the inclination a minimum (Fig. 3a). However, the fact that the initial inclination increases when the eccentricity decreases has a strong implication when including octupolar terms (Fig. 2): for high inclinations the trajectories are closer to the separatrix, which results in a higher oscillation of the eccentricity. Thus, for ω1 > 60° we observe that the evolution timescale is reduced again, since the inner orbit eccentricity is allowed to reach much higher values than those predicted by the quadrupolar approximation.

Up to now, we have been considering an initial mutual inclination I0 = 40°. For I0 ≲ 40°, there is only one dynamical regime for the HD 74256 planetary system, consisting of trajectories in circulation around the Lidov-Kozai equilibira (Fig. 3a). However, for higher values of the initial inclination we can also observe the libration regime (Fig. 3b). In Fig. 7 we show some numerical simulations with initial I0 = 50° using different values for the argument of the periastron.

For ω1 = 30° the inner planet is still in circulation (Fig. 7a), so we observe identical behavior for the eccentricity and inclination as in the case with I0 = 40° (Fig. 6). At the bottom of Fig. 7 we plot the inclination as a function of the argument of the periastron (I,ω1). We plot a dot each 105 yr and its color becomes darker with time. There we can clearly see that the planet is always in circulation, and that the amplitude of the oscillations is damped with time.

In the remaining two situations shown in Fig. 7, the inner planet is in libration around the Lidov-Kozai equilibrium located at ω1 = 90°. The effect of tides is to slowly increase the amplitude of both the eccentricity and inclination. As a consequence, the orbit of the planet will cross the separatrix of the libration zone and start to circulate as in the previous case.

For ω1 = 60° (Fig. 7b) the planet already starts close to the separatrix, so the initial oscillations are higher and it takes only about 75 Myr to cross it. For ω1 = 90° (Fig. 7c) the planet is placed close to the Lidov-Kozai equilibrium, so it takes much longer to reach the separatrix. In both situations, just after the transition of dynamical regime, the eccentricity reaches a very high value close to unity. Therefore, tidal effects with the central star become very strong and the final evolution is rapid: the semi-major axis decreases, and the inner orbit becomes circular. However, in a more realistic simulation where we integrate the full equations of motion and take into account the bodies dimensions, the planet most likely collides with the star. In both situations the system is either destroyed, or its configuration completely modified from the initial situation.

In order to better understand the variation of the evolution timescale with the initial choice of the argument of the pericenter, in Fig. 8 we plot the circularization time (e1 < 0.01) as a function of ω1. The circularization time is more or less equivalent to the inclination damping time, since the final stages of the evolution are very fast. For I0 = 40° the timescale decreases around the Lidov-Kozai equilibria, because the system is the circulation regime. For I0 = 50°, the system is in libration, but it is most likely destroyed after it crosses the separatrix.

Giuppone et al. (2012) also studied the evolution of planets inside the circulation zone of Lidov-Kozai equilibria. They performed some numerical simulations using the quadrupolar approximation and damping of the inner orbit eccentricity due to the presence of a primordial disk. They concluded that the planet stays in libration and migrates into the Lidov-Kozai equilibrium position, which is exactly the contrary that we observed here. Since the eccentricity of the inner orbit is also damped in our model, these results appear somehow contradictory. Therefore, we performed one simulation where the eccentricity is damped, but the semi-major axis is held constant. In this unrealistic situation, one observe that the planet migrates into the equilibrium like in Giuppone et al. (2012). As a consequence, it seems that there is no inconsistency between the two models, but it becomes clear that the semi-major axis evolution plays an important role in destabilizing the Lidov-Kozai equilibria. It appears that it cannot be neglected in future studies on the migration of the initial orbits as in Giuppone et al. (2012).

thumbnail Fig. 8

Circularization time (e1 < 0.01) for the HD 74156 system with I0 = 40° (red) and I0 = 50° (blue) using different values of the initial argument of the pericenter. For I0 = 40° the timescale decreases around the Lidov-Kozai equilibria, because the system is the circulation regime. For I0 = 50°, the system is in libration, but it is most likely destroyed after it crosses the separatrix.

Open with DEXTER

4.3. Constraints for the mutual inclination

In previous sections, we saw that in mutually inclined hierarchical two-planet systems there is a significant increase in the eccentricity of the inner planet’s orbit. As a result, the tidal dissipation is enhanced when the planet is at the periastron, and the system evolves faster into an equilibrium configuration.

In Fig. 9 we show the evolution of the HD 74156 system for three different values of the initial inclination I0 = 15°, 30° and 45°. As expected, when we increase the mutual inclination, the evolution timescale decreases. For I0 = 15° the eccentricity and the semi-major take more than 10 Gyr to be completely damped, while for I0 = 45° the system is fully evolved only after 100 Myr. Moreover, for I0 = 15° there is almost no effect on the mutual inclination, we only observe some amplitude damping when the eccentricity is decreased to low values, because gravitational perturbations no longer force the inclination. On the contrary, for I0 = 30°, the pumping effect on the eccentricity is already present, and hence we observe a significant reduction of the final mutual inclination.

thumbnail Fig. 9

Long-term evolution of the HD 74156 system for different values of the initial inclination, I0 = 15°a); 30°b); and 45°c). We plot the mutual inclination I (top), the eccentricities e1 (blue) and e2 (green) (middle), and the semi-major axis a1 (bottom).

Open with DEXTER

In Fig. 10 we show the same kind of evolution as before, but for an initial retrograde orbit with initial inclinations I0 = 165°, 150° and 135°. In this case, the evolution of the system does not differ much from the prograde situation, the only significant difference is that the inclination is damped to high values close to 180° (coplanar system with a retrograde orbit). For I0 = 165° the inclination is more or less conserved and the eccentricity is damped over 10 Gyr, while for lower values of the initial inclination, the inclination is damped and the system evolves in much shorter timescales.

thumbnail Fig. 10

Long-term evolution of the HD 74156 system for different values of the initial inclination, I0 = 165°a); 150°b); and 135°c). We plot the mutual inclination I (top), the eccentricities e1 (blue) and e2 (green) (middle), and the semi-major axis a1 (bottom).

Open with DEXTER

From the observation of Figs. 9 and 10 we then conclude that mutual inclinations closer to 90° speed up the final evolution of the system. Since most hierarchical systems listed in Table 1 (except HD 190360 and HD 38529) still present substantial values for the inner-orbit’s eccentricity, we then expect that their mutual inclinations are not extremely high. For all those systems we run several numerical simulations starting with the present initial conditions from Table 1, adopting k2Δt = 100 s, and different initial values for I and ω1. All trajectories that circularize the inner-orbit in less than ~10 Gyr can then be ruled out, while those not showing significant modifications can be retained as possible representations of the real system. Therefore, we are able to set some constraints for the maximal mutual inclination of each system, whose limitations are listed in Table 1. In Fig. 11 we show the outcome of these simulations for the HD 74156 system, which corresponds to a summary of the more detailed evolutions shown in previous Figures. Orbits with 20° < I < 150° circularize the system in less than 10 Gyr, so they can be discarded.

When we run the same kind of simulations for the HD 38529, we observe that the eccentricity is damped very quickly, even for coplanar orbits. One possibility is that we overestimated the dissipation. However, even if we adopt k2Δt = 10 s, that is, one order of magnitude lower than for the remaining planets, the system still circularizes in a timescale shorter than the age of the system. Another possibility is to suppose that the inner planet semi-major axis was higher in the past. This hypothesis can also be extended to the HD 190360 system, for which the inner orbit is already circularized, but it may have had a higher eccentricity value in the past. In order to test this scenario, for all planetary systems in Table 1, we run several numerical simulations for different initial values for I and ω1, keeping k2Δt = 100 s, but adopting a1 = 0.2 AU and e1 = 0.25 as initial values, instead of the current values. In Fig. 12 we plot an example for the HD 190360 system. By modifying the initial conditions, we are able to reproduce the present observations. Note that the HD 190360 system is older than the HD 38529 one (Table 1), so both systems may have undergone an identical evolution, but they are observed at different stages. The initial semi-major axis could have been higher, providing that the inner orbit eccentricity was also higher (for instance a1 = 0.25 AU and e1 = 0.4).

4.4. Effect of the dissipation rate

In the former sections we have been using k2Δt = 100 s, or, in terms of Q-factor, (Eq. (31)), which is similar to the present value measured for Saturn (Lainey et al. 2012). We adopted this value mainly for a better comparison with the previous paper on the planar case (Correia et al. 2012). Other works on the tidal evolution of hot-Jupiters also adopted similar values for Δt (e.g. Fabrycky & Tremaine 2007; Correia et al. 2011).

However, the Q-factor of planets is unknown, and may vary by some orders of magnitude (e.g. Goldreich & Soter 1966). Indeed, the value measured for Jupiter appears to differ by a factor of ten, (Lainey et al. 2009). On the other hand, statistical studies on the observed eccentricity distribution of hot-Jupiters give (e.g. Jackson et al. 2008; Hansen 2010). Part of the problem is that the nature of tidal dissipation in these planets is still poorly understood. In addition, the Q-factor is frequency dependent (i.e., model dependent), and therefore sometimes it is difficult to translate from one system to another. In order to test the robustness of the inclination damping, here we test our model with lower dissipation rates (higher Q values).

According to the tidal equations (Eqs. (25)−(28)), the evolution timescales are linearly proportional to , so higher Q-values delay the final evolution of a system. In Fig. 13 we plot the evolution of the HD 74156 system for three different dissipation values,  s (which is equivalent to k2/Q ≈ 10-6,   10-5,   10-4, respectively), starting with I0 = 50°. We observe that, although the evolution timescale is longer for higher values (as expected), the inclination damping is still present.

thumbnail Fig. 11

Circularization time (e1 < 0.01) for the HD 74156 system using different values of the initial mutual inclination and argument of the pericenter. Since the estimated age of the system is several Gyr (Table 1), we can rule out mutual inclinations within 20° < I < 150°.

Open with DEXTER

thumbnail Fig. 12

Possible past evolution of the HD 190360 system with I0 = 43°, a1 = 0.2 AU, and e1 = 0.25. We plot the mutual inclination I (top), the eccentricities e1 (blue) and e2 (green) (middle), and the semi-major axis a1 (bottom).

Open with DEXTER

In Sect. 4.3 we saw that when we increase the initial mutual inclination I0, the evolution timescale decreases very fast. For instance, in the case of the HD 74156 system, for I0 = 40° the inner orbit becomes fully damped after 1 Gyr, for I0 = 50° it takes about 100 Myr, and for I0 = 60° only 10 Myr. Therefore, for different dissipation rates we can still put constraints on the mutual inclination of hierarchical systems, the only consequence is that as we increase , the maximal mutual inclination that one can expect to observe also increases.

thumbnail Fig. 13

Possible past evolution of the HD 74156 system with I0 = 50°, using different dissipation rates. We plot the mutual inclination I for k2Δt = 100 s (top), k2Δt = 10 s (middle), and k2Δt = 1 s (bottom).

Open with DEXTER

5. Conclusion

Many two-planet systems have been reported in hierarchical configurations. For most of these systems the mutual inclinations are unknown. However, since the orbital eccentricities are typically high, we may expect that the formation mechanism that increased the eccentricities also increased the mutual inclinations. Very often the innermost planet in these systems is close enough to the star to undergo tidal dissipation, which can pose constraints on the final evolution.

Here we have studied a particular subgroup of hierarchical planetary systems in which the inner planet’s semi-major axis 0.1 < a1 < 0.5 AU (we called these “moderate close-in planets”). This range is important to ensure that the spin of the inner planet is fully evolved, but not its orbit. Using an averaged secular model that takes into account gravitational interactions up to octupole order, we showed that for many initial conditions the mutual inclination is damped to relatively low values (I ~ 15° for HD 74156) on timescales shorter than the age of the system (less than one Gyr).

Without planetary perturbations and for zero obliquity there is no effect from tides on the evolution of the inclination (Eq. (28)). The inclination damping is thus not a direct consequence of tidal effects on the orbits. The key element is a inner planet’s eccentricity oscillation at a secular timescale similar to the synchronization time of its spin. Indeed, the rotation (and thus the flattening (Eq. (2))) of the planet is driven by its eccentricity variations (Eq. (25)). In response to these excitations, the rotation is phase-shifted (Eq. (A.14)) and the lag tends to pump the eccentricity (Eq. (A.18)). Because the total angular momentum must be conserved, the increase in the inner orbit’s eccentricity is accompanied by a subsequent reduction of the mutual inclination. When the eccentricity pumping ceases, the inclination damping also stops.

For high mutual inclination values, quadrupolar gravitational exchanges with the eccentricity are more efficient, and so is the eccentricity pumping. As a consequence, the inner orbit’s eccentricity reaches higher values, tidal effects are enhanced at the periastron, and the system evolves on shorter timescales. A strong inclination damping is then often followed by a fast circularization of the inner orbit. Since most of the observed hierarchical systems still present substantial values of the inner-orbit’s eccentricity after several Gyr, we expect that they cannot have very high mutual inclinations. In particular, we are able to set constraints on the highest mutual inclination in these systems.

The evolution timescale also depends on the argument of the periastron in the invariant plane of the system, which is unknown at present for most systems. Indeed, the uncertainty on this parameter is much higher than the uncertainty on the dissipation time lag Δt. If ω1 is in circulation, the evolution can be extremely fast (a few million years) for values close to ω1 = 90° or ω1 = 270°. However, for high mutual inclinations ω1 can be in Lidoz-Kozai libration for values close to ω1 = 90° or ω1 = 270°. In that case, the evolution timescale can be delayed to several Gyr, but the equilibrium is unstable and broken when the system crosses the separatrix with the circulation regime. After that, the inner planet is most likely lost.


Acknowledgments

We wish to thank the referee Rosemary Mardling, who made very valuable suggestions. We acknowledge support by PNP-CNRS, CS of Paris Observatory, PICS05998 France-Portugal program, the European Research Council/European Community under the FP7 through a Starting Grant, and Fundacão para a Ciência e a Tecnologia, Portugal (grants PTDC/CTE-AST/098528/2008, SFRH/BSAB/1148/2011, PEst-C/CTM/LA0025/2011).

References

Appendix A: Inclination damping

According to expression (28), the direct contribution of tidal effects to the inclination can be neglected for low obliquities. We adopted ε = 0° for simplicity, but in a more realistic situation the obliquity is expected to be trapped in a Cassini state (e.g. Correia et al. 2011). Although high-obliquity states are possible, they are unlikely for close-in planets (Levrard et al. 2007), so we only expect these planets to be trapped in low-obliquity Cassini states, and hence dI/dt ≈ 0.

When tidal effects are combined with planetary perturbations, some counter-intuitive behavior may appear (e.g. Wu & Murray 2003; Mardling 2007; Correia et al. 2012). In particular, it has been shown that for moderate close-in planets, the gravitational perturbations of a distant coplanar companion combined with tides can increase the eccentricity of the inner orbit to very high values (Correia et al. 2012). Because of the total angular momentum conservation, the eccentricity of the outer orbit is simultaneously decreased.

For coplanar systems (sinI = 0), the octupolar terms in the expansion of the potential (Eq. (4)) are the main responsible for the eccentricity variations of the inner orbit (Eq. (13)). However, for mutually inclined systems, the quadrupolar terms become dominating, and we can neglect the octupole contribution: (A.1)In addition, the contributions to the outer orbit eccentricity vanish (ė2 ≈ 0), meaning that the angular momentum exchanges only occur between e1 and I (Eq. (3)). Thus, if somehow the quadrupole interactions are able to pump the inner orbit eccentricity as in the planar case (Correia et al. 2012), one can expect to observe a decrease in the mutual inclination of the system.

To understand this mechanism, we can simplify equations of motion without loss of generality. Since we assume that the orbits are not coplanar, we retain only the quadrupolar terms for the gravitational perturbations, that is, the conservative motion can be described solely by expression (A.1) and the first four terms in expression (15): (A.2)We can neglect tidal effects on orbital quantities (Eqs. (26), (27)), which is justified since (Eq. (30)). The only contribution of tides is then to the rotation rate (Eq. (25)). The semi-major axis and the mean motion are thus constant, as is the eccentricity of the outer orbit and G2.

For simplicity, we average expression (A.2) over ω1, and linearize the set of equations of motion in the vicinity of the averaged values of x, e1, and I. Let x = x0 + δx, where x0 is the solution of (37), e1 = e10 + δe1, and I = I0 + δI. In the following, δI is expressed as a function of δe1 using the angular momentum conservation (Eq. (3)): (A.3)Then, the set of equations of motion ((A.1), (A.2), and (25)) reduces to with where , and .

At first order, the precession of the periastron is constant , and the eccentricity is simply given from expression (A.4) as (A.13)where Δe = A/2g, and ω1 = gt + ω10. That is, the eccentricity e1 presents periodic variations around an equilibrium value e10, with amplitude Δe and frequency 2g. Since gxδx,geδe1 ≪ g, the above solution for the eccentricity can be adopted as the zeroth-order solution of the system of Eqs. (A.4)−(A.6). With this approximation, the equation of motion of δx (A.6) becomes that of a driven harmonic oscillator whose steady state solution is (A.14)with , and . The rotation rate thus presents an oscillation identical to the eccentricity (Eq. (A.13)), but with smaller amplitude and delayed by an angle φ (see Correia 2011). Using the above expression in Eq. (A.5) and integrating, gives for the periastron (A.15)Finally, substituting in expression (A.4) and using the approximation gxΔx,geΔe ≪ g gives (A.16)or, combining the two products of periodic functions, (A.17)The two middle terms in the above equation can be neglected since they are periodic and have a very small amplitude (gxΔx,geΔe ≪ g). However, the last term in sinφ is constant and it adds a small drift to the eccentricity, (A.18)Note that the phase lag φ between the eccentricity (Eq. (A.13)) and the rotation variations (Eq. (A.14)) is essential to obtain a drift on the eccentricity. From expressions (A.7) and (A.8) we have A ~ ν21 ~ g, while from expressions (A.11) and (A.12) we get . Consequently, (A.19)and (A.20)that is, (A.21)Thus, the drift vanishes when φ = 0 or 90 degrees, i.e., for strong dissipation (), where , or for weak dissipation (), where , respectively. The pumping effect on the eccentricity is then maximized when φ = 45°, for , with .

We therefore conclude that the quadrupolar perturbation of an inclined companion enhances the effect of the octupolar terms in the planar case (Correia et al. 2012). However, since the angular momentum is exchanged with the inclination, by assuming G2 ≫ G1 we obtain from expression (A.3) that for I0 < 90°, and for I0 > 90°. That is, the mutual inclination is damped as long as the inner orbit eccentricity is pumped.

The main difference when we consider the full non-linearized problem is that the drift in the eccentricity (Eq. (A.18)) cannot grow indefinitely. Indeed, when the eccentricity reaches high values, the drift vanishes (Fig. 5b). Moreover, tidal effects are also enhanced for high eccentricities and counterbalance the drift (Eq. (27)). Hence, the drift in the eccentricity is never permanent, neither is the inclination damping, although they can last for the entire age of the system.

All Tables

Table 1

Hierarchical two-planet non-resonant systems with 0.1 < a1 < 0.5 AU and a1/a2 < 0.1.

Table 2

Stability analysis of the HD 74156 system for different sets (ω1,I) of initial conditions (Fig. 2).

All Figures

thumbnail Fig. 1

Radius versus the mass of the planet. We plot all known close-in planets in the range 0.1 < a1 < 0.5 AU, for which the radius was determined by the transiting method. We observe that the radius decreases with the mass in a relatively regular way.

Open with DEXTER
In the text
thumbnail Fig. 2

Possible secular trajectories for the HD 74156 system seen in the (ω1,I) plane (left), and in the (ω1,e1) plane (right). We show the trajectories using the quadrupolar approximation (top), corresponding to the level curves of constant energy, and the octupolar approximation (bottom). The dashed black curves in the (ω1,I) plane correspond to the separatrix between the circulation and libration zones of ω1. The colors are preserved in all pictures, each one corresponding to a given value of the total angular momentum of the system (determined by different values of the initial mutual inclination). All trajectories are compatible to the present knowledge for this system (Table 2).

Open with DEXTER
In the text
thumbnail Fig. 3

Possible secular trajectories for the HD 74156 system seen in the (ω1,e1) plane for two different values of the initial mutual inclination, I0 = 40° (top), and I0 = 50° (bottom). We show the trajectories using the quadrupolar approximation, corresponding to the level curves of constant energy. Each one corresponds to an initial value of the argument of the periastron ω1, ranging from 0° (black) to 90° (red) with a step of 10°. The dashed black curves correspond to the observed eccentricity of the planet, that is, it gives the initial condition for ω1. All trajectories are compatible with the current knowledge for this system (Table 1).

Open with DEXTER
In the text
thumbnail Fig. 4

Possible secular trajectories for the HD 74156 system seen in the (ϕ,e1) plane for I = 0° (top) and in the (ω1,e1) plane for I = 35° (bottom). All trajectories are compatible with the current knowledge for this system (Table 1). We show the trajectories using the octupolar approximation (blue) and direct numerical simulations (red). In a) the blue path corresponds to the level curves of constant energy for coplanar orbits. The dot marks the present position of the planet. In b) we additionally show the trajectories using the quadrupolar approximation (green), which corresponds to the level curves of constant energy.

Open with DEXTER
In the text
thumbnail Fig. 5

Long-term evolution of the HD 74156 system for different tidal models with I0 = 40°: only dissipation in the orbit is considered a); only dissipation in the spin is considered b); full model c). We plot the mutual inclination I (top), the eccentricities e1 (blue) and e2 (green) (middle), and the semi-major axis a1 (bottom).

Open with DEXTER
In the text
thumbnail Fig. 6

Long-term evolution of the HD 74156 system with I0 = 40° for different values of the argument of the periastron, ω1 = 30°a); 50°b); and 70°c). We plot the mutual inclination I (top), the eccentricities e1 (blue) and e2 (green) (middle), and the semi-major axis a1 (bottom).

Open with DEXTER
In the text
thumbnail Fig. 7

Long-term evolution of the HD 74156 system with I0 = 50° for different values of the argument of the periastron, ω1 = 30°a); 60°b); and 90°c). As a function of time, we plot the mutual inclination I (top) and the eccentricities e1 (blue) and e2 (green) (middle). As a function of ω1, we plot again the mutual inclination, where the color of each dot becomes darker with time (bottom).

Open with DEXTER
In the text
thumbnail Fig. 8

Circularization time (e1 < 0.01) for the HD 74156 system with I0 = 40° (red) and I0 = 50° (blue) using different values of the initial argument of the pericenter. For I0 = 40° the timescale decreases around the Lidov-Kozai equilibria, because the system is the circulation regime. For I0 = 50°, the system is in libration, but it is most likely destroyed after it crosses the separatrix.

Open with DEXTER
In the text
thumbnail Fig. 9

Long-term evolution of the HD 74156 system for different values of the initial inclination, I0 = 15°a); 30°b); and 45°c). We plot the mutual inclination I (top), the eccentricities e1 (blue) and e2 (green) (middle), and the semi-major axis a1 (bottom).

Open with DEXTER
In the text
thumbnail Fig. 10

Long-term evolution of the HD 74156 system for different values of the initial inclination, I0 = 165°a); 150°b); and 135°c). We plot the mutual inclination I (top), the eccentricities e1 (blue) and e2 (green) (middle), and the semi-major axis a1 (bottom).

Open with DEXTER
In the text
thumbnail Fig. 11

Circularization time (e1 < 0.01) for the HD 74156 system using different values of the initial mutual inclination and argument of the pericenter. Since the estimated age of the system is several Gyr (Table 1), we can rule out mutual inclinations within 20° < I < 150°.

Open with DEXTER
In the text
thumbnail Fig. 12

Possible past evolution of the HD 190360 system with I0 = 43°, a1 = 0.2 AU, and e1 = 0.25. We plot the mutual inclination I (top), the eccentricities e1 (blue) and e2 (green) (middle), and the semi-major axis a1 (bottom).

Open with DEXTER
In the text
thumbnail Fig. 13

Possible past evolution of the HD 74156 system with I0 = 50°, using different dissipation rates. We plot the mutual inclination I for k2Δt = 100 s (top), k2Δt = 10 s (middle), and k2Δt = 1 s (bottom).

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.