Issue 
A&A
Volume 541, May 2012



Article Number  A151  
Number of page(s)  10  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201118356  
Published online  17 May 2012 
Dynamical analysis and constraints for the HD 196885 system
^{1} Department of Physics, I3NUniversity of Aveiro, Campus Universitário de Santiago, 3810193 Aveiro, Portugal
email: cristian@ua.pt
^{2} Astronomie et Systèmes Dynamiques, IMCCECNRS UMR8028, 77 av. DenfertRochereau, 75014 Paris, France
^{3} Centro de Astrofísica, Universidade do Porto, rua das Estrelas, 4150762 Porto, Portugal
Received: 28 October 2011
Accepted: 23 March 2012
The HD 196885 system is composed of a binary star and a planet orbiting the primary. The orbit of the binary is fully constrained by astrometry, but for the planet the inclination with respect to the plane of the sky and the longitude of the node are unknown. Here we perform a full analysis of the HD 196885 system by exploring the two free parameters of the planet and choosing different sets of angular variables. We find that the most likely configurations for the planet are either nearly coplanar orbits (prograde and retrograde), or highly inclined orbits near the LidovKozai equilibrium points, i = 44° or i = 137°. Among coplanar orbits, the retrograde ones appear to be less chaotic, while for the orbits near the LidovKozai equilibria, those around ω = 270° are more reliable, where is the argument of pericenter of the planet’s orbit with respect to the binary’s orbit . From the observer’s point of view (plane of the sky) stable areas are restricted to (I_{1},Ω_{1}) ~ (65°,80°), (65°,260°), (115°,80°), and (115°,260°), where I_{1} is the inclination of the planet and Ω_{1} is the longitude of ascending node.
Key words: celestial mechanics / planets and satellites: dynamical evolution and stability / stars: individual: HD 196885 / binaries: close
© ESO, 2012
1. Introduction
About 20% of all known exoplanets have been found to inhabit multiple stellar systems (Desidera & Barbieri 2007; Mugrauer & Neuhäuser 2009). The theories of planet formation around a star that is part of a binary system pose considerable challenges. A brief review of theories can be found in Quintana et al. (2002) and Thebault (2011). However, it is very important to know the real masses and the spatial configuration of these systems to better understand the processes involved in their formation.
A few planets in binaries are found in compact systems with the semimajor axis of the binary less than 100 AU: HD 196885 (Correia et al. 2008; Chauvin et al. 2011), Gl 86 (Queloz et al. 2000; Lagrange et al. 2006), γCep (Hatzes et al. 2003; Neuhäuser et al. 2007), and HD 41004 (Zucker et al. 2004). Even for these tight systems, the poor precision of firstepoch observations together with an incomplete time span of the observations led to several bestfit solutions with almost the same residuals (e.g. Torres 2007; Correia et al. 2008). It is therefore advisable to keep this wide range of possible configurations in mind, and to combine radial velocity fits with possible formation scenarios (e.g. γCephei, Giuppone et al. 2011).
Recent observations of the HD 196885 system allowed constraining orientation of the binary orbit by combining astrometric and radial velocity observations (Chauvin et al. 2011). Thus, the real mass of the stellar companion was established at 0.45 M_{⊙}. A few numerical integrations were carried out by Chauvin et al. (2011), but it is not clear which the compatible regions of solutions for the planet around the central star are and why some individual coplanar solutions are unstable. We intend to clarify this picture by performing massive numerical integrations over the entire space of the free parameters and compare the results with analytical models. This will allow us to completely clarify all possible dynamical regimes in the HD 196885 system, and to put constraints on future observations.
2. Bestfit solution and reference angles
HD 196885 B, the stellar companion of HD 196885 A, was detected by combining imaging and spectroscopic observations as be an M1 ± IV dwarf located at 0.7′′, which corresponds to 23 AU in projected physical separation (Chauvin et al. 2007).
The planetary companion, HD 196885 Ab, was detected using radial velocity data from ELODIE, CORALIE, and CORAVEL observations spread over 14 years (Correia et al. 2008). The orbital solution for the planet gave a minimum mass of m_{Ab} sin i = 2.96 M_{Jup}, a period of P = 3.69 ± 0.03 yr, and an eccentricity of e = 0.462 ± 0.026. The authors additionally constrained the binary companion HD 196885 B with a period of 40 < P < 120 yr, a semimajor axis 14 < a < 30 AU, and a minimum mass of 0.3 < m_{B}sini < 0.6 M_{⊙}. Based on Lick observations, Fischer et al. (2009) confirmed this solution for the system. In both studies, a relatively wide range of masses and periods were found compatible with residuals for the stellar companion.
More recently, Chauvin et al. (2011) combined five astrometric measurements obtained with VLT/NACO spread over four years and all sources of radial velocity data (four sets, summarizing 187 observations) to determine the full orbit of the binary in the space and a consistent projection of planetary parameters. We show the best fit obtained by Chauvin et al. (2011) in Table 1.
The characterization of the orbits in the HD 196885 system by Chauvin et al. (2011) left only two parameters undetermined, the inclination of the planet with respect to the plane of the sky, I_{1}, and the longitude of the node, Ω_{1} (Table 1). It is then possible to cover the entire phasespace for this system by exploring these two missing angles.
Before studying the phasespace of the system, it is convenient to understand all physical configurations. A scheme of the fundamental planes for the definition of the reference angles is shown in Fig. 1. For simplicity, we reserve for the inner orbit (planet) the index 1, and for the outer orbit (stellar companion) the index 2. The central star has mass m_{0} = 1.31 M_{⊙} (Chauvin et al. 2011), the planet has mass m_{1}, and the stellar companion has mass m_{2}, with semimajor axis a_{i}, eccentricity e_{i}, mean anomaly M_{i}, argument of pericenter ω_{i}, longitude of ascending node Ω_{i}, and inclination I_{i}. We also mark some additional angles that are useful for our study: the mutual inclination i, the nodal longitude of the planet’s orbit with respect to the binary’s orbit , the difference of node longitudes ΔΩ = Ω_{1} − Ω_{2} in the plane of the sky, the argument of pericenter of the planet’s orbit with respect to the binary’s orbit , and the angle .
Fig. 1
Fundamental planes for the definition of the reference angles. The angles in black are determined directly from the observations (Table 1), while the two angles in blue (I_{1},ΔΩ) are the missing angles in the observations that prevent the full characterization of the orbits in the system. The angles in red are also undetermined, but related with the other angles. They are linked to the physical system, are therefore independent of the reference frame, and more suitable to study the dynamics. 
The dynamics of the system can be studied using different sets of the free parameters, the most important being (I_{1},ΔΩ), (), and (). The observed quantities obtained from fits to the data are deduced with respect to plane of the sky (I_{1},Ω_{1}), but the dynamics of inclined systems is more adequately described using as reference the stellar companion (the orbit 2), since it is fully determined (Table 1) and is not much disturbed by the planet.
All reference angles are interrelated, therefore whatever choice we adopt, one can easily determine the remaining angles by spherical trigonometry transformations (e.g. Smart 1965). If one studies the system for mutually inclined configurations (), we obtain for the remaining variables noting that when , we have 0° < Δω < 180° and 0° < ΔΩ < 180°, and for , we have 180° < Δω < 360° and 180° < ΔΩ < 360°.
Fig. 2
Possible secular trajectories for the HD 196885 system seen in the plane (top), and in the plane (bottom). We show the trajectories in the restricted quadrupolar approximation (left) and using nbody numerical simulations of the real system (right). The dashed black curves in a) separate zones of libration and circulation of the angle . The black crosses indicate the location of LidovKozai equilibrium points for this system. Nbody numerical integrations in the () plane b), and in the () plane d), preserve the same colors for the initial conditions in the restricted problem. 
Alternatively, if one prefers to explore the space of possible solutions starting from the observer’s point of view (I_{1},Ω_{1}), the spherical triangle defined by the sides has three known parameters ΔΩ,I_{1},I_{2} and the following relations can be used to determine the physical system:
3. Analytical model
Before studying the full massive problem, it is useful to look at the restricted inner problem (the orbit of the outer companion is fixed). This approximation is not very far from being true, because the mass of the planet is much smaller than the mass of the stellar bodies (m_{1} ≪ m_{0},m_{2}), and can therefore be seen as a “test particle” (m_{1} = 0). The restricted problem is easier to study and allows us to determine the dynamical regimes that can be expected in the HD 196885 system.
3.1. Restricted quadrupolar problem
We consider a binary star system composed of a primary (m_{0}) and a secondary (m_{2}), and a massless planet (m_{1} = 0) orbiting the primary star. The binary system’s fixed orbit with period T_{2}, semimajor axis a_{2} and eccentricity e_{2}, is the natural choice of reference frame. The planet’s osculating orbit with period T_{1}, semimajor axis a_{1} and eccentricity e_{1} has an orientation with respect to the binary system’s orbit defined by the angles i (relative inclination), (argument of the pericenter) and (longitude of the ascending node). Following Kozai (1962) and Kinoshita & Nakai (1999, 2007), we write the Hamiltonian of the planet, expand up to quadrupole order in the semimajor axis ratio a_{1}/a_{2}, and average with respect to the fast periods T_{1} and T_{2}, obtaining the secular quadrupole Hamiltonian (7)with^{1}(8)The secular motion can be described by writing Hamilton’s equations using Delaunay canonical variables where is the gravitational constant. Since the Hamiltonian (Eq. (7)) does not depend on , then the conjugate momentum H is constant. Moreover, because F (Eq. (7)) has only one degree of freedom, , we obtain secular trajectories by plotting level curves F = const., with (11)In our particular case (the HD 196885 system), we know the current value of the planet’s eccentricity (e_{10} = 0.48), but we do not know the initial i_{0} or ω_{0}. We can obtain possible secular trajectories for the planet by choosing values of i_{0} and ω_{0}, then plotting the level curves of (12)or (13)with .
Fig. 3
MEGNO chaos indicators for different sets of variables, a) (), b) (I_{1},ΔΩ), c) () with , and d) () with . Labeled lines give the mass of the planet in M_{Jup}. Orbits can be considered stable for ⟨ Y ⟩ ≤ 2. 
We show the secular trajectories described by expressions (12) and (13) in Figs. 2a and c, respectively. Because of the symmetry in the solution space we use the same color for trajectories with relative inclination i (prograde) or 180° − i (retrograde), and for librating trajectories around or . The separatrixes (dashed black curves in Fig. 2a) mark the boundary between libration and circulation of the angle and are obtained by solving the implicit equation (Kinoshita & Nakai 2007)^{2}(14)with for .
In agreement with Lidov (1961, 1962) and Kozai (1962), the regimes of secular motion consist of LidovKozai equilibrium points (with i and e_{1} fixed) and LidovKozai cycles which owing to the conservation of h (Eq. (11)) exhibit coupled oscillations in i and e_{1}. An orbit with e_{1} = 0.48 is at a LidovKozai equilibrium point if , , or . LidovKozai cycles with librating have h < 0.6 and F < F_{s} (magenta, brown, cyan, yellow, and pink orbits). LidovKozai cycles with circulating have h > 0.6 (red, blue, and green orbits) or h < 0.6 and F > F_{s} (gray, orange, and purple orbits). Moreover, highamplitude LidovKozai cycles (gray, pink and orange orbits) can reach e_{1} near unity and could thus become unstable in the full problem. In particular, an orbit that passes through i = 90° must, by conservation of h (Eq. (11)), reach e_{1} = 1 (collision orbit). An example of this is the black orbit in Fig. 2, which started at i = 90° and (although it is not obvious from Fig. 2a because Eq. (12) is undefined in this region).
Note that Figs. 2a,c are different from the standard LidovKozai diagrams that show level curves of the Hamiltonian F at constant values h = h_{0}. In our case e_{10} is fixed, but we show trajectories for different values of the initial inclination i_{0}, i.e., with different h_{0}. In particular, this explains why librating orbits in the LidovKozai regime do not encircle a LidovKozai equilibrium point (magenta, brown, cyan, yellow, and pink orbits in Fig. 2). The location of the LidovKozai equilibrium points in our diagrams occur at the current planet’s eccentricity (e_{1} = 0.48), hence circulating and librating orbits starting at nor will have eccentricity increasing (magenta, cyan, yellow and pink orbits) or decreasing (red, blue, green, and purple orbits) from this initial value (Fig. 2).
3.2. Numerical simulations
To compare the orbital behavior in the full massive problem with the restricted problem, we performed numerical integrations up to t = 30 kyr starting with the same initial conditions as in Figs. 2a,c (initial conditions from Table 1 and with i from 10 to 170 degrees, and with i from 10 to 170 degrees). We show in Figs. 2b,d these numerical integrations as well as the position of the separatrixes and the LidovKozai equilibrium points calculated with the quadrupolar approximation. We used the same colors as in Figs. 2a,c to facilitate the comparison between the two models.
We see that the numerical integrations of the full problem (Figs. 2b,d) agree well with the theoretical trajectories of the quadrupolar restricted Hamiltonian (Figs. 2a,c). The global behavior of the planet in the HD 196885 system is dominated by the dynamical regimes described above, i.e., LidovKozai equilibrium points, librating and circulating orbits. In the full problem, the secular trajectories exhibit a drift from the secular solutions of the quadrupole Hamiltonian. This behavior can be explained by including octupole and higherorder terms in the Hamiltonian (Ford et al. 2000; Laskar & Boué 2010; Lithwick & Naoz 2011). The octupole Hamiltonian is not independent of and thus h (given by Eq. (11)) is not a conserved quantity. Therefore, real secular trajectories (Figs. 2c,d) drift from the theoretical curves (Figs. 2a,c) because of changes in h.
We also expect deviations from the restricted quadrupole model solutions if the planet’s mass becomes comparable to the star’s masses, although, according to Farago & Laskar (2010), the topology of the phase space should not change. We see that this drift is more evident for orbits in the vicinity of the separatrixes and orbits that reach e_{1} ≈ 1 (nearly collision orbits). Note that the black trajectory in Figs. 2b,d is the collision orbit shown in Figs. 2a,c which is unstable after only 600 yr. Moreover, librating orbits around seem to drift more than those around , which could be due to the high planet mass solutions in this region (Sect. 4, Fig. 3).
Also, the brown orbit in Fig. 2d is not a double loop structure but two different initial conditions, one for (i < 90°) and the other for (i > 90°). In Fig. 2c both solutions coincide, while in Fig. 2d those with i < 90° have a bigger loop than those with i > 90°.
4. Dynamical analysis
Several works on dynamics in close binary systems exist. Among others, Rabl & Dvorak (1988) and Holman & Wiegert (1999), investigated the longterm stability of planets in coplanar circular orbits near one of the stars. Expressions of critical semimajor axis for the planet (prograde orbits) were derived as a function of the mass of binary components and eccentricity of the orbit. According to these approximations, prograde orbits for this system are stable with semimajor axes less than ~3.8 AU. Wiegert & Holman (1997) studied the stability of hypothetical terrestrial planets in the system αCentauri (whose semimajor axis and eccentricity are almost the same as in our system) for several mutual inclinations. However, they fixed the longitude of the node and varied the semimajor axis (in our case the semimajor axis of the planet is very well established with the radial velocity technique). The authors also noticed that highly mutually inclined orbits are unstable.
In this section we analyze the dynamics and the stability of the planetary system given in Table 1. The best choice of variables to study the global dynamics is the pair (), since it allows us to easily identify the fixed points and the different dynamical regimes (Sect. 3). However, this choice is not adequate because ω_{1} imposes constraints on , and the physical system has some restrictions in this frame (Eq. (2) or (5)). Consequently, we are left with the pair of angles (), which cover all possible configurations and are still independent of the observer.
4.1. Orbital stability
We constructed grids of initial conditions with 0.5° of resolution and each point in the grid was then numerically integrated over 30 kyr using a BurlischStoer based in an Nbody code (precision better than 10^{12}) using astrocentric and osculating variables. During the integrations we computed the averaged MEGNO chaos indicator ⟨ Y ⟩ : the regular orbits yield to ⟨ Y ⟩ ≤ 2, while higher values are indicative of chaotic motion (Cincotta & Simó 2000). The MEGNO chaos maps use a threshold that should be applied to avoid excluding stable orbits that did not converge to their theoretical value or those orbits that are weakly chaotic. Accordingly, following Maffione et al. (2011), the color scale shows “stable” orbits in blue up to ⟨ Y ⟩ ~ 2.5 (a particular choice based on integration of individual orbits for very long times and due to the characteristics of this system).
Results are shown in Fig. 3. We also show the grids in the observational frame (I_{1},ΔΩ), and in the frame () for comparison. We can see that highly mutually inclined systems (the strip within 70° < i < 110°) are highly unstable (this fact was already noticed in Wiegert & Holman (1997) for αCentauri but with a much poorer resolution).
To better understand the chaos regions, we superimposed contour level curves for the mass of the inner planet (in M_{Jup}). As we change the values of () we obtain corresponding value of I_{1} (Eq. (1)) and consequently the real mass of the planet ranges between 2.98 M_{Jup} (dashed lines) to 40 M_{Jup} (with higher values confined to little circles). Inside these circles the mass of the planet tends to infinity. Mass level curves become straight lines in the observational frame, since m_{1} ∝ 1/sinI_{1}, but the dynamical regimes are not perceptible, except instability corresponding to high masses for I_{1} = 0° and 180°.
In the frame (), all main structures become understandable: the system is unstable for high values of the mutual inclination, for high masses, but also around the separatrixes between dynamical regimes (Sect. 3). In addition, we observe that some initial conditions are not possible near highly mutual inclinations (i ~ 90°). The reason for this “forbidden” zone is that some initial conditions in the frame () are degenerate and correspond to identical solutions in the frame (), although with different values for the mass of the planet.
The size of the “forbidden” area depends on the inclination I_{2} of the stellar companion to the plane of the sky (Eqs. (1), (2)). For I_{2} = 90° there is a perfect correspondence between the frames () and (), while for I_{2} = 0° or 180°, there is a global degeneracy. For the HD 196885 system, we have I_{2} = 116.8°, that is, about 27° above 90°. As a consequence, the size of the degenerate area is about ± 27° around i = 90° in the frame (). The delimitation of this “forbidden” zone also coincides with the areas of large masses (Fig. 3c, d). In the particular case of the HD 196885 system, in the areas with overlapping solutions, one always corresponds to unstable orbits (the one that is close to the borders of the zone). Therefore, for simplicity, we can merge the two figures in a single diagram, ignoring the areas with i < 90° in Fig. 3c, and i > 90° in Fig. 3d. For the remainder of the paper we adopt this procedure in the frame ().
We used the MEGNO chaos indicator because it allows us to rapidly distinguish between regular and chaotic orbits. MEGNO is very useful in this type of studies because we do not need to check every initial condition for long periods of time. We can just follow a set of orbits in representative regions of the MEGNO map to understand how chaos leads to instability. Although the MEGNO indicator only provides indications on the regularity of the orbits, we verified that the chaotic orbits marked in red are indeed unstable in the sense that the planet reaches distances to the central star shorter than the radius of the star itself (we set the limit at 0.005 AU, which is the radius of the Sun, an overestimated quantity because for this system the central star has 1.3 R_{Sun}). The timescales of this instability depend of the chaotic region that is considered:

1.
Around the Kozai separatrices (red regions), chaotic orbitsbecome unstable on short timescales (less than ~25 000 years).

2.
When the initial mutual inclination is between 80° ≲ i ≲ 100°, the orbits are very unstable and the planet collides with the central star in less than 1000 years. As we explained in Sect. 3, an orbit that passes through i = 90° must, by conservation of h, reach e_{1} = 1 (collision orbit).

3.
Near the edge of the forbidden region: (i ~ 60°,ω < 180°) and (i ~ 118°,ω < 180°), the planets have masses > 40 M_{jup}(i_{1} ~ 0°), i.e. the strong interactions with both stars cause the planet either to collide with the central star or to escape from the system.

4.
There is a wide region of unstable orbits with high values of MEGNO inside the lefthand Kozai separatrix (55 < i < 60° and 45° ≲ ω ≲ 125°) where the collision times spans 10^{5} to 10^{7} years.
Within the Kozai libration islands around ω = 270° the stable initial conditions have eccentricities that remain always lower than 0.90. Inside the Kozai libration islands around ω = 90°, the retrograde conditions show the same behavior, while some prograde conditions only survive for times from 10^{7} to 10^{8} years. Outside the separatrix, in the prograde region, the chaotic orbits (in green at bottom left hand of the graph, where ⟨ Y ⟩ > 2.5 and 4° < i < 10°) are stable and never reach eccentricities above 0.7.
Longterm numerical simulations showed that the regions marked in blue/green in the MEGNO map do not present detectable dynamical instability, at least on timescales on the order of 10^{9} years. Owing to the high mass ratio, the effect of highorder meanmotion resonances (MMR, hereafter) may be nonnegligible (the period ratio is 20.95). Indeed, there are a multitude of highorder MMR whose overlap could be responsible for the slow chaos regions. However, these highorder MMR are difficult to identify, and their relative importance will depend on the exact orbital elements of the system. Even secular effects could be the explanation for the slow chaos regions in green. We note that chaos nearby the Kozai separatrix has already been observed in the octupole problem (Lithwick & Naoz 2011).
We prefer to conclude that the best representation of the system corresponds to a region of regular motion (the region marked in blue). Although this may seem arbitrary, it is important to recall that so far, there are no giant planets displaying significant chaotic motion in any known planetary system. Secondly, regions of regular motion are expected to be more robust with respect to additional perturbations and planetary formation.
4.2. Dynamical regimes and regular orbits
MEGNO is a very efficient tool to identify chaotic motion, but it is not suited to distinguish between different types of regular orbits, that is, orbits that are more or less unperturbed, or orbits that experience intense secular variations. To evaluate the orbital behavior, we show in Fig. 4 the amplitude of the inclination variations for the three frames, together with the separatrixes between dynamical regimes obtained in Sect. 3.
Fig. 4
Amplitude of the mutual inclination for different sets of variables, a) (), b) (I_{1},ΔΩ), c) (). Centers of the LidovKozai equilibria are marked as black circles and the separatrix is calculated with the quadrupolar approximation (Sect. 3). 
We observe that unstable regions (Fig. 3) mainly coincide with the regions of large amplitude variations of the mutual inclination. Indeed, according to expression (11) and as discussed in Sect. 3, large variations in the inclination up to i ≈ 90° can increase the eccentricity of the inner planet to very high values, which gives rise to collision with the central star.
The more regular orbits of the system are those that are close to the LidovKozai equilibrium points, and those corresponding to coplanar orbits. Also notorious is that retrograde orbits (i > 140°) are in general less perturbed, probably because close encounters last a shorter time. As we have seen before (Fig. 3), the zones of higher mass for the planet introduce some instability, but the global picture is dominated by the dynamical regimes: trajectories that are close to the separatrix and nearly collision orbits are clearly the most unstable (cf. Sect. 3).
4.3. Eccentricity of the inner planet
The eccentricity of the inner planet is a key variable for understanding the stability and evolution of the system, because high values lead to collision with the central star. Moreover, while the outer companion may destabilize the orbit of the planet, encounters with the central star may give rise to tidal effects and a subsequent secular evolution of the orbit (Sect. 4.5). In addition, among all orbital parameters listed on Table 1, the eccentricity e_{1} is the only one for which secular modifications can be observed using radial velocity data (the argument of the periastron ω_{1} also varies, but it is not directly related to the physical orbit). In Fig. 5 we show the secular period associated with the largest eccentricity oscillation of the inner planet (in years), the amplitude of these eccentricity oscillations, and the maximum eccentricity, e_{1}, attained during the integration.
Fig. 5
Dynamical indicators in the frame (). a) period associated to highest amplitude oscillation of e_{1} (in years), b) maximal amplitude of oscillation of e_{1}, c) maximum eccentricity e_{1}. 
Fig. 6
Time variation of eccentricity for the orbits in Fig. 2 showing the two secular periods. The color code is the same as in Fig. 2. Initial conditions outside the separatrix (top), and inside the separatrix (bottom). 
As discussed previously, while the behavior of highinclination orbits (i.e. LidovKozai cycles) is welldescribed with the quadrupole Hamiltonian, the secular evolution of nearly coplanar orbits must take into account the octupole Hamiltonian. In particular, the eccentricity of coplanar orbits changes secularly solely due to the octupole term (Lee & Peale 2003). The eccentricity of lowtomedium inclination orbits is modulated by the octupole and quadrupole secular periods (Krymolowski & Mazeh 1999). Moreover, the quadrupole secular period is typically shorter than the octupole secular period (Fig. 6). Octupole and higherorder terms are also responsible for making the orbits near the separatrixes chaotic Lithwick & Naoz (2011).
In Fig. 5a we see the secular period associated with the largest eccentricity, e_{1}, variations. From Fig. 5a we see that retrograde orbits have secular periods longer than direct orbits. Among stable areas, mediumamplitude LidovKozai cycles have the largest eccentricity variations (Figs. 6 and 2c,d) and periods of around 1000 yr (Figs. 6 and 5a). The eccentricity oscillations of small amplitude LidovKozai cycles occur on a 1500 yr timescale, while nearly coplanar orbits also exhibit small oscillations in e_{1} but on a timescale of around 5500 yr (Fig. 6). Although we may be able to measure eccentricity variations in the radial velocity data (for instance, by detecting a drift in the data), in practice, we will have only a short observation timespan (typically a few years). Therefore, it will be difficult to distinguish small amplitude LidovKozai cycles from nearly coplanar orbits solely based on the different secular timescales (e.g. compare the magenta and red orbits in Fig. 6). We may, however, be able to identify if the planet is on a mediumamplitude LidovKozai cycle owing to the large eccentricity variations on a short (1000 yr) timescale.
The amplitude of the eccentricity oscillations (Fig. 5b) does not differ much from the amplitude of the mutual inclination oscillations (Fig. 4c). This behavior was expectable because these two variables are correlated (Eq. (11)).
In Fig. 5c we plot the quantity log _{10}(1. − e_{1,max}), where e_{1,max} is the maximal eccentricity attained by the planet during the integration. Thus a value of log _{10}(1. − e_{1,max}) = −3 (in red) represents an eccentricity of e_{1,max} = 0.999. Lightblue regions are those with e_{1,max} ~ 0.9, darkblue for those with e_{1,max} ~ 0.7, and finally lightgray for those where e_{1,max} ~ 0.5. The initial conditions within violet region inside the KozaiLidov separatrix can survive for times from 10^{8} to 10^{9} years. Note that this confirms that the orbits marked in red in MEGNO map collide with the central star because e ~ 1 (as described in Sect. 4.1).
4.4. Different initial conditions
Because the orbital parameters in Table 1 contain some undeterminations (Chauvin et al. 2011), we also explored the stability of the system within one σ of the bestfit parameters. In particular, we aimed to test the impact of variations in the outer companion, which is the most unconstrained. For that purpose, we made grids at the upper limits a_{2} = 21.86 AU and e_{2} = 0.45, and at bottom ones a_{2} = 20.14 AU and e_{2} = 0.39 as in Fig. 3. We do not show the results because unstable (red regions) are located in the same positions (i ~ 90°, and around the Kozai separatrix). Only regions of moderate chaos (green regions) slightly varied their positions and size. However, those regions are still possible solutions. Therefore, we conclude that the global picture in Fig. 3 does not change much for different sets of initial conditions, that is, all conclusions in this paper are still valid within the errorbars of the published parameters.
We also made a grid for the restricted case (i.e. the mass of the planet m_{1} = 0) using the bestfit parameters from Table 1. The entire chaotic green regions disappeared and we were left only with regular regions (blue) and closeencounter regions (in red).
4.5. Tidal evolution
For some orbital configurations, the eccentricity of the inner planet may reach very high values and become close enough to the central star at periastron to undergo tidal effects. For an unperturbed orbit, the secular evolution of the eccentricity by tidal effect can be given by (Correia 2009) (15)with (16)and (17)where k_{2} is the second Love number, Q is the tidal dissipation factor, R is the radius of the planet, and The solution of the above equation is given by (Correia & Laskar 2010) (18)where F(e) is an implicit function of e, which converges to zero as t → + ∞. The characteristic timescale for fully dampening the eccentricity of the orbit is then τ ~ 1/K. The above expression is only valid for unperturbed orbits, but in the HD 196885 system the initial eccentricity can be seen as the maximum eccentricity attained on a Kozai cycle, e_{1,max}, since it is at this point that dissipation is maximized. Therefore, τ provides a minimal estimation of the dampening time. We can exclude the orbits for which the eccentricity can be damped during the age of the system from the possible orbital parameters of the planet, since the present eccentricity is still near 0.5 (Table 1).
In Fig. 7 we show the evolution timescale for different initial maximum eccentricities, adopting k_{2} = 0.5 and R = 1.2 R_{Jup}. We observe that for Q ~ 1000, a value similar to Jupiter, only for e_{1,max} > 0.96 the eccentricity is damped within 2 Gyr, the estimated age of the system (Correia et al. 2008). Even for Earthlike planets (Q ~ 10), only e_{1,max} > 0.92 can be dissipated in the same amount of time. Since e_{1} > 0.96 corresponds to unstable orbits, we conclude that tidal effects cannot be used to exclude any of the stable initial conditions compatible with the observational data (Table 1).
Fig. 7
Damping time for the eccentricity of the inner planet (with τ = 1/K, Eq. (17)). For Q ~ 1000, only for e_{1,max} > 0.96 the eccentricity is damped within 2 Gyr, the estimated age of the system (τ/Q = 2 × 10^{3} Gyr, red line). 
5. Formation hypothesis
From the formation point of view, the orbital excitation exerted by LidovKozai cycles implies destructive, highvelocity collisions among planetesimals, inhibiting formation of massive objects (Lissauer 1993). However, the LidovKozai equilibrium generates its own protection mechanism to achieve stable configurations when mutually inclined bodies are considered.
Thebault (2011) published a deep study of the possible scenarios of formation for the planet in the HD 196885 binary system. The author considered an axisymmetric static gas disk (no selfgravity) and estimated the impact velocities amongst a population of circumprimary planetesimals. A main conclusion was that the circumprimary disk is strongly hostile to planetesimal accretion, especially the region around 2.6 AU (the planet’s present location) where the binary perturbations induce planetesimalshattering velocities of more than 1 km s^{1}. The region around 2.6 AU is strongly hostile to planetesimal accretion, even for highly inclined orbits, and alternative solutions were proposed to justify their existence (e.g., different initial configuration of binary or disk instability). Quintana et al. (2002) studied the formation in the αCentauri system using a disk of combined large and small bodies. These authors concluded that it is possible to form planets when an inclined disk is considered, although ~ 95% of the initial mass is lost by the end of the simulation. Their results predict a wider region to form more retrograde orbits (i = 180°) than prograde ones, and that possibly large planetary embryos could be formed near the star (from 0.5 to 1.5 AU) in inclined orbits. However, because the initial semimajoraxis distribution has an upper limit set at ~ 2 AU, this may bias the final results for the allowed regions.
However, Batygin et al. (2011) addressed the formation for wide binary systems (a_{binary} ~ 1000 AU), concluding that is possible to form a single planet in an inclined orbit (LidovKozai regime), if taking into account the selfgravity of protoplanetary disk (that means planetesimals embedded in the gaseous disk). They also pointed out that the evolutionary process of formation of a planetary system at the LidovKozai equilibrium is necessarily nonunique and proposed a scenario were a multiple system might be formed protected by LidovKozai cycles and subsequent instabilities remove all the remaining planets. Unfortunately the authors did not mention to close binary systems in their approach.
The secular dynamics of small planetesimals play a fundamental role in establishing the possibility of accretional collisions in such extreme cases. The most important secular parameters are the forced eccentricity and the secular frequency, which depend on the initial conditions of the particles, as well as on the mass and orbital parameters of the secondary star (e.g. Thebault 2011; Beaugé et al. 2010). However, Giuppone et al. (2011) pointed out that for these types of compact systems, the frequencies are sometimes not satisfactorily determined from firstorder approximations (secondorder on mass are needed), and therefore, some works on formation based on these approximations should probably be reviewed.
The giant planets in the solar system, which are supposed to have migrated owing to their interaction with a disk of planetesimals (e.g. Tsiganis et al. 2005), probably experienced the same as what happened with the orbit of the planet in the HD 196885 system. Using the secular model described in Correia et al. (2011) for the full system (two stars and one planet), we performed some simple experiments on the early evolution of the system by migrating hypothetical initial orbital parameters of the planet into the present ones.
We independently tested the evolution of the semimajor axis and the eccentricity of the inner orbit, using an exponential decay (Beaugé et al. 2006; Lee et al. 2007): (19)with τ_{m} = 10 Myr. There is no variation in the orbital configuration for the semimajor axis, since the ratio a_{1}/a_{2} is a factor in the quadrupolar Hamiltonian (Eq. (8)), therefore it only modifies the timescale of the evolution. However, for the eccentricity we observe that all trajectories that begin inside the libration zone migrate into one of the LidovKozai equilibrium positions (Fig. 8). Trajectories in the circulation zone remain more or less unchanged, only the amplitude of the mutual inclination is damped. Therefore, we can assume that if the planet in the HD 196885 was able to form in libration, it will be presently observed near the LidovKozai equilibria.
Fig. 8
Early evolution of the planet when the eccentricity is damped from higher values due to the interaction with a disk of planetesimals. For each simulation, the color of the position of the planet becomes darker with time. In the circulation zone only the amplitude of the mutual inclination is damped, but for orbits starting in the libration area, the planet evolves into the LidovKozai equilibria. 
6. Discussion and conclusions
We have studied the dynamics of the HD 196885 system using the present determination of the orbital parameters (Table 1) as a starting point. We developed an insightfull analysis of the 3D space exploring the free parameters and choosing different sets of angular variables. We found that the most likely configurations for the planet in the HD 196885 system is either nearly coplanar orbits (prograde and retrograde), or highly inclined orbits near the LidovKozai equilibrium points, i − 90° = ± 47°.
Among coplanar orbits, the retrograde ones appear to be less chaotic, possibly because close encounters last shorter than in the prograde case^{3}. The orbits at the LidovKozai equilibrium around ω = 270° are more reliable than those around 90°, since the mass of the planet is smaller in the first situation. Present formation scenarios for this kind of systems are unable to distinguish between the two possibilities (coplanar or LidovKozai equilibrium). However, a simple simulation on the evolution of the initial eccentricity of the inner orbit shows that if the planet was able to form in the libration zone, it will evolve into the LidovKozai equilibrium point. We also tested the effect of tides, but they are too weak in the HD 196885 system, although they should be taken into account for tighter systems of this kind.
Although there is a wide variety of stable initial conditions for the planet in the natural frame defined by the orbit of the two stars (), from the observer’s point of view (plane of the sky) there are some restrictions. Indeed, looking at Fig. 4b we see that stable areas occur around
Ω_{1} ~ 80°  Ω_{1} ~ 260°;  
I_{1} ~ 65°  Kozai prograde  coplanar retrograde; 
Therefore, by continuing to observe the system, one will first be able to determine the dynamical regime of the planet (coplanar or LidovKozai equilibrium) and later its exact position. Moreover, it may be possible to identify medium amplitude Kozai cycles because of their large eccentricity oscillations on 1000 yr timescale. Smallamplitude Kozai cycles and nearly coplanar orbits have small eccentricity oscillations on distinct timescales (1500 yr and 4000 yr timescales, respectively). However, in practice, it may be difficult to distinguish these two configurations because we are limited to short observation timespans.
The expression for C (Eq. (7)) in Kinoshita & Nakai (2007) should have and not in the nominator.
The solution to Eq. (14) is independent of h, as we would expect, since the separatrix is an invariant curve.
An alternative explanation based on resonances’ overlap was recently given in Morais & Giuppone (2012).
Acknowledgments
We would like to express our gratitude to the referee for his detailed analysis of our paper and his comments and suggestions, which helped to improve the manuscript significantly. This work has been supported by the European Research Council/European Community under the FP7 through a Starting Grant. We also acknowledge financial support from FCTPortugal (grants PTDC/CTEAST/098528/2008 and PEstC/CTM/LA0025/2011).
References
 Batygin, K., Morbidelli, A., & Tsiganis, K. 2011, A&A, 533, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beaugé, C., Michtchenko, T. A., & FerrazMello, S. 2006, MNRAS, 365, 1160 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Beaugé, C., Leiva, A. M., Haghighipour, N., & Otto, J. C. 2010, MNRAS, 408, 503 [NASA ADS] [CrossRef] [Google Scholar]
 Chauvin, G., Lagrange, A. M., Udry, S., & Mayor, M. 2007, A&A, 475, 723 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chauvin, G., Beust, H., Lagrange, A.M., & Eggenberger, A. 2011, A&A, 528, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cincotta, P. M., & Simó, C. 2000, A&AS, 147, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M. 2009, ApJ, 704, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2010, in Exoplanets (University of Arizona Press), 534 [Google Scholar]
 Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celest. Mech. Dyn. Astron., 53 [Google Scholar]
 Correia, A. C. M., Udry, S., Mayor, M., et al. 2008, A&A, 479, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desidera, S., & Barbieri, M. 2007, A&A, 462, 345 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Farago, F., & Laskar, J. 2010, MNRAS, 401, 1189 [NASA ADS] [CrossRef] [Google Scholar]
 Fischer, D., Driscoll, P., Isaacson, H., et al. 2009, ApJ, 703, 1545 [NASA ADS] [CrossRef] [Google Scholar]
 Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Giuppone, C. A., Leiva, A. M., CorreaOtto, J., & Beaugé, C. 2011, A&A, 530, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hatzes, A. P., Cochran, W. D., Endl, M., et al. 2003, ApJ, 599, 1383 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Kinoshita, H., & Nakai, H. 1999, Celest. Mech. Dyn. Astron., 75, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Kinoshita, H., & Nakai, H. 2007, Celest. Mech. Dyn. Astron., 98, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Krymolowski, Y., & Mazeh, T. 1999, MNRAS, 304, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Lagrange, A.M., Beust, H., Udry, S., Chauvin, G., & Mayor, M. 2006, A&A, 459, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J., & Boué, G. 2010, A&A, 522, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lee, M. H., & Peale, S. J. 2003, ApJ, 592, 1201 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, M. H., Peale, S. J., Pfahl, E., & Ward, W. R. 2007, Icarus, 190, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Lidov, M. L. 1961, Iskus. sputniky Zemly (in Russian), 8, 5 [Google Scholar]
 Lidov, M. L. 1962, Plan. Space Sci., 9, 719 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J. 1993, ARA&A, 31, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Lithwick, Y., & Naoz, S. 2011, ApJ, 742, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Maffione, N. P., Darriba, L. A., Cincotta, P. M., & Giordano, C. M. 2011, Celest. Mech. Dyn. Astron., 111, 285 [NASA ADS] [CrossRef] [Google Scholar]
 Morais, M. H. M., & Giuppone, C. A. 2012, MNRAS, in press [arXiv:1204.4718] [Google Scholar]
 Mugrauer, M., & Neuhäuser, R. 2009, A&A, 494, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neuhäuser, R., Mugrauer, M., Fukagawa, M., Torres, G., & Schmidt, T. 2007, A&A, 462, 777 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Queloz, D., Mayor, M., Weber, L., et al. 2000, A&A, 354, 99 [NASA ADS] [Google Scholar]
 Quintana, E. V., Lissauer, J. J., Chambers, J. E., & Duncan, M. J. 2002, ApJ, 576, 982 [NASA ADS] [CrossRef] [Google Scholar]
 Rabl, G., & Dvorak, R. 1988, A&A, 191, 385 [NASA ADS] [Google Scholar]
 Smart, W. M. 1965, Textbook on spherical astronomy (Cambridge University Press) [Google Scholar]
 Thebault, P. 2011, Celest. Mech. Dyn. Astron., 111, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, G. 2007, ApJ, 654, 1095 [NASA ADS] [CrossRef] [Google Scholar]
 Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wiegert, P. A., & Holman, M. J. 1997, AJ, 113, 1445 [NASA ADS] [CrossRef] [Google Scholar]
 Zucker, S., Mazeh, T., Santos, N. C., Udry, S., & Mayor, M. 2004, A&A, 426, 695 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
All Tables
All Figures
Fig. 1
Fundamental planes for the definition of the reference angles. The angles in black are determined directly from the observations (Table 1), while the two angles in blue (I_{1},ΔΩ) are the missing angles in the observations that prevent the full characterization of the orbits in the system. The angles in red are also undetermined, but related with the other angles. They are linked to the physical system, are therefore independent of the reference frame, and more suitable to study the dynamics. 

In the text 
Fig. 2
Possible secular trajectories for the HD 196885 system seen in the plane (top), and in the plane (bottom). We show the trajectories in the restricted quadrupolar approximation (left) and using nbody numerical simulations of the real system (right). The dashed black curves in a) separate zones of libration and circulation of the angle . The black crosses indicate the location of LidovKozai equilibrium points for this system. Nbody numerical integrations in the () plane b), and in the () plane d), preserve the same colors for the initial conditions in the restricted problem. 

In the text 
Fig. 3
MEGNO chaos indicators for different sets of variables, a) (), b) (I_{1},ΔΩ), c) () with , and d) () with . Labeled lines give the mass of the planet in M_{Jup}. Orbits can be considered stable for ⟨ Y ⟩ ≤ 2. 

In the text 
Fig. 4
Amplitude of the mutual inclination for different sets of variables, a) (), b) (I_{1},ΔΩ), c) (). Centers of the LidovKozai equilibria are marked as black circles and the separatrix is calculated with the quadrupolar approximation (Sect. 3). 

In the text 
Fig. 5
Dynamical indicators in the frame (). a) period associated to highest amplitude oscillation of e_{1} (in years), b) maximal amplitude of oscillation of e_{1}, c) maximum eccentricity e_{1}. 

In the text 
Fig. 6
Time variation of eccentricity for the orbits in Fig. 2 showing the two secular periods. The color code is the same as in Fig. 2. Initial conditions outside the separatrix (top), and inside the separatrix (bottom). 

In the text 
Fig. 7
Damping time for the eccentricity of the inner planet (with τ = 1/K, Eq. (17)). For Q ~ 1000, only for e_{1,max} > 0.96 the eccentricity is damped within 2 Gyr, the estimated age of the system (τ/Q = 2 × 10^{3} Gyr, red line). 

In the text 
Fig. 8
Early evolution of the planet when the eccentricity is damped from higher values due to the interaction with a disk of planetesimals. For each simulation, the color of the position of the planet becomes darker with time. In the circulation zone only the amplitude of the mutual inclination is damped, but for orbits starting in the libration area, the planet evolves into the LidovKozai equilibria. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.