Tidal damping of the mutual inclination in hierarchical systems
^{1}
Departamento de Física, I3NUniversidade de Aveiro,
Campus de Santiago,
3810193
Aveiro,
Portugal
email:
correia@ua.pt
^{2}
ASD, IMCCE – CNRS, UMR8028, Observatoire de Paris,
UPMC, 77 Av.
DenfertRochereau, 75014
Paris,
France
^{3}
Department of Astronomy and Astrophysics, University of
Chicago, 5640 S. Ellis
Ave, Chicago,
IL
95064,
USA
Received:
2
October
2012
Accepted:
18
February
2013
Hierarchical twoplanet systems, in which the inner body’s semimajor axis is between 0.1 and 0.5 AU, usually present high eccentricity values, at least for one of the orbits. As a result of the formation process, one may expect that planetary systems with high eccentricities also have high mutual inclinations. However, here we show that tidal effects combined with gravitational interactions damp the initial mutual inclination to modest values in timescales that are shorter than the age of the system. This effect is not a direct consequence of tides on the orbits, but it results from a secular forcing of the inner planet’s flattening. We then conclude that these hierarchical planetary systems are unlikely to present very high mutual inclinations, at least as long as the orbits remain outside the LidovKozai libration areas. The present study can also be extended to systems of binary stars and to planetsatellite systems.
Key words: celestial mechanics / planets and satellites: dynamical evolution and stability / planetstar interactions
© ESO, 2013
1. Introduction
Today, nearly 100 multiplanet systems have been reported, of which roughly 1/3 possess at least one “moderate closein planet”, that is, a planet with a semimajor axis between 0.1 and 0.5 AU^{1}. 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 closein 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 twoplanet systems, there is a special class whose semimajor ratio a_{1}/a_{2} is lower than 0.1, the socalled “hierarchical systems”. This class counts at least 20 members (1/5 of all multiplanet systems), and usually at least one of the planets’ orbits is highly eccentric^{1}. 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 semimajor 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 semimajor ratio is small, loworder 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 semimajor axis a_{1}/a_{2} 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 a_{1}/a_{2}, 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 socalled LidovKozai mechanism). For planar problems, the series expansions in a_{1}/a_{2} 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 a_{1}/a_{2} 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 a_{1} < 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 spinorbit 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 twoplanet 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 abovementioned 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 twoplanet 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 m_{0}, an inner planet of mass m_{1}, and an outer companion of mass m_{2}. We used Jacobi canonical coordinates, with r_{1} being the position of m_{1} relative to m_{0}, and r_{2} the position of m_{2} relative to the center of mass of m_{1} and m_{0}. We assumed that the system is hierarchical, thus  r_{1} ≪ r_{2} . 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 G_{i} the orbital angular momentum of each body.
The inner planet is considered an oblate ellipsoid with gravity field coefficients given by J_{2}, 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 k_{2} 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 G_{1} 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, quadrupolelevel for the spin (e.g. Correia & Laskar 2010a), octupolelevel 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)a_{i} is the semimajor axis (that can also be expressed using the mean motion n_{i}), e_{i} is the eccentricity, and ω_{i} is the argument of the periastron. We also have β_{1} = m_{0}m_{1}/(m_{0} + m_{1}), β_{2} = (m_{0} + m_{1})m_{2}/(m_{0} + m_{1} + m_{2}), , 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 e_{1}, e_{2} and I are related by the conservation of the total angular momentum (Eq. (3)), we have (11)As we assumed that , if we neglect firstorder terms in L/G_{1}, 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°.
Hierarchical twoplanet nonresonant systems with 0.1 < a_{1} < 0.5 AU and a_{1}/a_{2} < 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 semimajor axis and the eccentricity can only decrease until the orbit of the planet becomes circular. However, planetplanet 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 multiplanetary systems may show nonintuitive 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 twoplanet 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 < a_{1} < 0.5 AU, which we call “moderate closein 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 closein 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 g_{x}). The minimum masses and the semimajor axis are relatively well determined from the observations, so the largest incertitudes in ν_{1} come from the Love number k_{2}, 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.
Fig. 1 Radius versus the mass of the planet. We plot all known closein planets in the range 0.1 < a_{1} < 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 closein planets in the range 0.1 < a_{1} < 0.5 AU. We found ten planets in this range^{1} 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 R_{Jup} for Jupiter and 1.0 R_{Sat} for Saturn.
3.2. Initial conditions uncertainty
Assuming that the observational values of the minimum masses, semimajor 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.
Fig. 2 Possible secular trajectories for the HD 74156 system seen in the (ω_{1},I) plane (left), and in the (ω_{1},e_{1}) 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, I_{i} (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 I_{2} ≈ 48°, estimated using astrometric measurements from the Hubble Space Telescope (Benedict et al. 2010). The planetary masses m_{1} and m_{2} given in Table 1 correspond to the minimum masses (assuming I_{1} = I_{2} = 90°), except for HD 38529, where the masses were estimated using I_{1} = I_{2} = 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 G_{1} and G_{2} 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 I_{i} 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 semimajor 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 e_{2} is constant (Eq. (10)). The dynamics of the system is then fully described by the couples (I,ω_{1}) or (e_{1},ω_{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 e_{1} and I. However, as the inclination increases, the dynamics of the system is considerably perturbed by the presence of LidovKozai 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 LidovKozai diagrams that show level curves of the quadrupole Hamiltonian at fixed values G_{1} cosI, because for HD 74156 the initial eccentricity of the inner planet is already fixed at e_{1} = 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 LidovKozai regime do not encircle a LidovKozai 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 twoplanet 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 semimajor ratios is 0.03 < a_{1}/a_{2} < 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 (ϕ,e_{1}) or (ϕ,e_{2}). 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 I_{1} = I_{2} = 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.
Fig. 3 Possible secular trajectories for the HD 74156 system seen in the (ω_{1},e_{1}) plane for two different values of the initial mutual inclination, I_{0} = 40° (top), and I_{0} = 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 
Fig. 4 Possible secular trajectories for the HD 74156 system seen in the (ϕ,e_{1}) plane for I = 0° (top) and in the (ω_{1},e_{1}) 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 (e_{1},ω_{1},e_{2},ω_{2}), and it becomes impossible to capture the dynamics of the systems in a twodimensional 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 I_{0} = 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 , k_{2} = 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 I_{0} = 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 closein planets evolves very fast, we assume that the spin is locked in its equilibrium position (Eq. (37)). We observe that the eccentricity and the semimajor 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 wellknow 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 nonplanar 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).
Fig. 5 Longterm evolution of the HD 74156 system for different tidal models with I_{0} = 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 e_{1} (blue) and e_{2} (green) (middle), and the semimajor axis a_{1} (bottom). 

Open with DEXTER 
Fig. 6 Longterm evolution of the HD 74156 system with I_{0} = 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 e_{1} (blue) and e_{2} (green) (middle), and the semimajor axis a_{1} (bottom). 

Open with DEXTER 
Fig. 7 Longterm evolution of the HD 74156 system with I_{0} = 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 e_{1} (blue) and e_{2} (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 semimajor 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 e_{1} the inclination decreases, and viceversa. In the example from previous section (Fig. 5c), we used ω_{1} = 0°, that is, we assumed that the observed value of the eccentricity (e_{1} = 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 I_{0} = 40°. For I_{0} ≲ 40°, there is only one dynamical regime for the HD 74256 planetary system, consisting of trajectories in circulation around the LidovKozai 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 I_{0} = 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 I_{0} = 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 10^{5} 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 LidovKozai 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 LidovKozai 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 semimajor 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 (e_{1} < 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 I_{0} = 40° the timescale decreases around the LidovKozai equilibria, because the system is the circulation regime. For I_{0} = 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 LidovKozai 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 LidovKozai 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 semimajor 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 semimajor axis evolution plays an important role in destabilizing the LidovKozai equilibria. It appears that it cannot be neglected in future studies on the migration of the initial orbits as in Giuppone et al. (2012).
Fig. 8 Circularization time (e_{1} < 0.01) for the HD 74156 system with I_{0} = 40° (red) and I_{0} = 50° (blue) using different values of the initial argument of the pericenter. For I_{0} = 40° the timescale decreases around the LidovKozai equilibria, because the system is the circulation regime. For I_{0} = 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 twoplanet 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 I_{0} = 15°, 30° and 45°. As expected, when we increase the mutual inclination, the evolution timescale decreases. For I_{0} = 15° the eccentricity and the semimajor take more than 10 Gyr to be completely damped, while for I_{0} = 45° the system is fully evolved only after 100 Myr. Moreover, for I_{0} = 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 I_{0} = 30°, the pumping effect on the eccentricity is already present, and hence we observe a significant reduction of the final mutual inclination.
Fig. 9 Longterm evolution of the HD 74156 system for different values of the initial inclination, I_{0} = 15°a); 30°b); and 45°c). We plot the mutual inclination I (top), the eccentricities e_{1} (blue) and e_{2} (green) (middle), and the semimajor axis a_{1} (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 I_{0} = 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 I_{0} = 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.
Fig. 10 Longterm evolution of the HD 74156 system for different values of the initial inclination, I_{0} = 165°a); 150°b); and 135°c). We plot the mutual inclination I (top), the eccentricities e_{1} (blue) and e_{2} (green) (middle), and the semimajor axis a_{1} (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 innerorbit’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 k_{2}Δt = 100 s, and different initial values for I and ω_{1}. All trajectories that circularize the innerorbit 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 k_{2}Δ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 semimajor 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 k_{2}Δt = 100 s, but adopting a_{1} = 0.2 AU and e_{1} = 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 semimajor axis could have been higher, providing that the inner orbit eccentricity was also higher (for instance a_{1} = 0.25 AU and e_{1} = 0.4).
4.4. Effect of the dissipation rate
In the former sections we have been using k_{2}Δt = 100 s, or, in terms of Qfactor, (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 hotJupiters also adopted similar values for Δt (e.g. Fabrycky & Tremaine 2007; Correia et al. 2011).
However, the Qfactor 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 hotJupiters 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 Qfactor 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 Qvalues 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 k_{2}/Q ≈ 10^{6}, 10^{5}, 10^{4}, respectively), starting with I_{0} = 50°. We observe that, although the evolution timescale is longer for higher values (as expected), the inclination damping is still present.
Fig. 11 Circularization time (e_{1} < 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 
Fig. 12 Possible past evolution of the HD 190360 system with I_{0} = 43°, a_{1} = 0.2 AU, and e_{1} = 0.25. We plot the mutual inclination I (top), the eccentricities e_{1} (blue) and e_{2} (green) (middle), and the semimajor axis a_{1} (bottom). 

Open with DEXTER 
In Sect. 4.3 we saw that when we increase the initial mutual inclination I_{0}, the evolution timescale decreases very fast. For instance, in the case of the HD 74156 system, for I_{0} = 40° the inner orbit becomes fully damped after 1 Gyr, for I_{0} = 50° it takes about 100 Myr, and for I_{0} = 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.
Fig. 13 Possible past evolution of the HD 74156 system with I_{0} = 50°, using different dissipation rates. We plot the mutual inclination I for k_{2}Δt = 100 s (top), k_{2}Δt = 10 s (middle), and k_{2}Δt = 1 s (bottom). 

Open with DEXTER 
5. Conclusion
Many twoplanet 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 semimajor axis 0.1 < a_{1} < 0.5 AU (we called these “moderate closein 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 phaseshifted (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 innerorbit’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 LidozKozai 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 PNPCNRS, CS of Paris Observatory, PICS05998 FrancePortugal 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/CTEAST/098528/2008, SFRH/BSAB/1148/2011, PEstC/CTM/LA0025/2011).
References
 Bean, J. L., McArthur, B. E., Benedict, G. F., & Armstrong, A. 2008, ApJ, 672, 1202 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Benedict, G. F., McArthur, B. E., Bean, J. L., et al. 2010, AJ, 139, 1844 [NASA ADS] [CrossRef] [Google Scholar]
 Butler, R. P., Wright, J. T., Marcy, G. W., et al. 2006, ApJ, 646, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M. 2009, ApJ, 704, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M. 2011, in IAU Symp. 276, eds. A. Sozzetti, M. G. Lattanzi, & A. P. Boss, 287 [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2004, Nature, 429, 848 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2010a, Icarus, 205, 338 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2010b, in Exoplanets (University of Arizona Press), 534 [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2012, ApJ, 751, L43 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Laskar, J., & Néron de Surgy, O. 2003, Icarus, 163, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Udry, S., Mayor, M., et al. 2005, A&A, 440, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celest. Mech. Dynamic. Astron., 111, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Boué, G., & Laskar, J. 2012, ApJ, 744, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Couetdic, J., Laskar, J., Correia, A. C. M., Mayor, M., & Udry, S. 2010, A&A, 519, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Efroimsky, M., & Williams, J. G. 2009, Celest. Mech. Dynamic. Astron., 104, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298 [NASA ADS] [CrossRef] [Google Scholar]
 Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Giguere, M. J., Fischer, D. A., Howard, A. W., et al. 2012, ApJ, 744, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Giuppone, C. A., Morais, M. H. M., Boué, G., & Correia, A. C. M. 2012, A&A, 541, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goldreich, P., & Soter, S. 1966, Icarus, 5, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, B. M. S. 2010, ApJ, 723, 285 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, B., Greenberg, R., & Barnes, R. 2008, ApJ, 678, 1396 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Lainey, V., Arlot, J.E., Karatekin, Ö., & van Hoolst, T. 2009, Nature, 459, 957 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lainey, V., Karatekin, Ö., Desmars, J., et al. 2012, ApJ, 752, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Lambeck, K. 1988, Geophysical geodesy: the slow deformations of the earth Lambeck (Oxford, UK: Clarendon Press; New York: Oxford University Press) [Google Scholar]
 Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1993, Phys. D Nonlinear Phenomena, 67, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Boué, G. 2010, A&A, 522, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J., Boué, G., & Correia, A. C. M. 2012, A&A, 538, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lee, M. H., & Peale, S. J. 2003, ApJ, 592, 1201 [NASA ADS] [CrossRef] [Google Scholar]
 Levrard, B., Correia, A. C. M., Chabrier, G., et al. 2007, A&A, 462, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lidov, M. L. 1961, Iskus. sputniky Zemly [in Russian], 8, 5 [Google Scholar]
 Lidov, M. L. 1962, Planet. Space Sci., 9, 719 [NASA ADS] [CrossRef] [Google Scholar]
 Marchal, C. 1990, The ThreeBody Problem (Amsterdam: Elsevier) [Google Scholar]
 Mardling, R. A. 2007, MNRAS, 382, 1768 [NASA ADS] [Google Scholar]
 Meschiari, S., Laughlin, G., Vogt, S. S., et al. 2011, ApJ, 727, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Mignard, F. 1979, Moon and Planets, 20, 301 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge University Press) [Google Scholar]
 Naef, D., Mayor, M., Beuzit, J. L., et al. 2004, A&A, 414, 351 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nagasawa, M., Ida, S., & Bessho, T. 2008, ApJ, 678, 498 [NASA ADS] [CrossRef] [Google Scholar]
 Pilyavsky, G., Mahadevan, S., Kane, S. R., et al. 2011, ApJ, 743, 162 [NASA ADS] [CrossRef] [Google Scholar]
 Pont, F., Husnoo, N., Mazeh, T., & Fabrycky, D. 2011, MNRAS, 414, 1278 [NASA ADS] [CrossRef] [Google Scholar]
 Ségransan, D., Udry, S., Mayor, M., et al. 2010, A&A, 511, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Singer, S. F. 1968, Geophys. J. R. Astron. Soc. , 15, 205 [CrossRef] [Google Scholar]
 Touma, J. R., Tremaine, S., & Kazandjian, M. V. 2009, MNRAS, 394, 1085 [NASA ADS] [CrossRef] [Google Scholar]
 Vogt, S. S., Butler, R. P., Marcy, G. W., et al. 2005, ApJ, 632, 638 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, J. T., Upadhyay, S., Marcy, G. W., et al. 2009, ApJ, 693, 1084 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y., & Goldreich, P. 2002, ApJ, 564, 1024 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y., & Murray, N. 2003, ApJ, 589, 605 [NASA ADS] [CrossRef] [Google Scholar]
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 highobliquity states are possible, they are unlikely for closein planets (Levrard et al. 2007), so we only expect these planets to be trapped in lowobliquity Cassini states, and hence dI/dt ≈ 0.
When tidal effects are combined with planetary perturbations, some counterintuitive behavior may appear (e.g. Wu & Murray 2003; Mardling 2007; Correia et al. 2012). In particular, it has been shown that for moderate closein 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 e_{1} 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 semimajor axis and the mean motion are thus constant, as is the eccentricity of the outer orbit and G_{2}.
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, e_{1}, and I. Let x = x_{0} + δx, where x_{0} is the solution of (37), e_{1} = e_{10} + δe_{1}, and I = I_{0} + δI. In the following, δI is expressed as a function of δe_{1} 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 e_{1} presents periodic variations around an equilibrium value e_{10}, with amplitude Δe and frequency 2g. Since g_{x}δx,g_{e}δe_{1} ≪ g, the above solution for the eccentricity can be adopted as the zerothorder 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 g_{x}Δx,g_{e}Δ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 (g_{x}Δx,g_{e}Δ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 G_{2} ≫ G_{1} we obtain from expression (A.3) that for I_{0} < 90°, and for I_{0} > 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 nonlinearized 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
Hierarchical twoplanet nonresonant systems with 0.1 < a_{1} < 0.5 AU and a_{1}/a_{2} < 0.1.
All Figures
Fig. 1 Radius versus the mass of the planet. We plot all known closein planets in the range 0.1 < a_{1} < 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 
Fig. 2 Possible secular trajectories for the HD 74156 system seen in the (ω_{1},I) plane (left), and in the (ω_{1},e_{1}) 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 
Fig. 3 Possible secular trajectories for the HD 74156 system seen in the (ω_{1},e_{1}) plane for two different values of the initial mutual inclination, I_{0} = 40° (top), and I_{0} = 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 
Fig. 4 Possible secular trajectories for the HD 74156 system seen in the (ϕ,e_{1}) plane for I = 0° (top) and in the (ω_{1},e_{1}) 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 
Fig. 5 Longterm evolution of the HD 74156 system for different tidal models with I_{0} = 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 e_{1} (blue) and e_{2} (green) (middle), and the semimajor axis a_{1} (bottom). 

Open with DEXTER  
In the text 
Fig. 6 Longterm evolution of the HD 74156 system with I_{0} = 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 e_{1} (blue) and e_{2} (green) (middle), and the semimajor axis a_{1} (bottom). 

Open with DEXTER  
In the text 
Fig. 7 Longterm evolution of the HD 74156 system with I_{0} = 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 e_{1} (blue) and e_{2} (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 
Fig. 8 Circularization time (e_{1} < 0.01) for the HD 74156 system with I_{0} = 40° (red) and I_{0} = 50° (blue) using different values of the initial argument of the pericenter. For I_{0} = 40° the timescale decreases around the LidovKozai equilibria, because the system is the circulation regime. For I_{0} = 50°, the system is in libration, but it is most likely destroyed after it crosses the separatrix. 

Open with DEXTER  
In the text 
Fig. 9 Longterm evolution of the HD 74156 system for different values of the initial inclination, I_{0} = 15°a); 30°b); and 45°c). We plot the mutual inclination I (top), the eccentricities e_{1} (blue) and e_{2} (green) (middle), and the semimajor axis a_{1} (bottom). 

Open with DEXTER  
In the text 
Fig. 10 Longterm evolution of the HD 74156 system for different values of the initial inclination, I_{0} = 165°a); 150°b); and 135°c). We plot the mutual inclination I (top), the eccentricities e_{1} (blue) and e_{2} (green) (middle), and the semimajor axis a_{1} (bottom). 

Open with DEXTER  
In the text 
Fig. 11 Circularization time (e_{1} < 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 
Fig. 12 Possible past evolution of the HD 190360 system with I_{0} = 43°, a_{1} = 0.2 AU, and e_{1} = 0.25. We plot the mutual inclination I (top), the eccentricities e_{1} (blue) and e_{2} (green) (middle), and the semimajor axis a_{1} (bottom). 

Open with DEXTER  
In the text 
Fig. 13 Possible past evolution of the HD 74156 system with I_{0} = 50°, using different dissipation rates. We plot the mutual inclination I for k_{2}Δt = 100 s (top), k_{2}Δt = 10 s (middle), and k_{2}Δt = 1 s (bottom). 

Open with DEXTER  
In the text 