EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 541 (May 2012) A&A, 541 (2012) A151 Full HTML
Free Access
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/0004-6361/201118356
Published online 17 May 2012

© 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 first-epoch observations together with an incomplete time span of the observations led to several best-fit 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. Best-fit 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 mAb   sin   i = 2.96   MJup, 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 < mBsini < 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.

Table 1

Orbital parameters for the HD 196885 system obtained by Chauvin et al. (2011) for JD = 2 455 198.

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, I1, and the longitude of the node, Ω1 (Table 1). It is then possible to cover the entire phase-space for this system by exploring these two missing angles.

Before studying the phase-space 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 m0 = 1.31   M⊙ (Chauvin et al. 2011), the planet has mass m1, and the stellar companion has mass m2, with semimajor axis ai, eccentricity ei, mean anomaly Mi, argument of pericenter ωi, longitude of ascending node Ωi, and inclination Ii. 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 .

thumbnail 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 (I1,ΔΩ) 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 (I1,ΔΩ), (), and (). The observed quantities obtained from fits to the data are deduced with respect to plane of the sky (I1,Ω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 inter-related, 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°.

thumbnail 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 n-body 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 Lidov-Kozai equilibrium points for this system. N-body 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 (I1,Ω1), the spherical triangle defined by the sides has three known parameters ΔΩ,I1,I2 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 (m1 ≪ m0,m2), and can therefore be seen as a “test particle” (m1 = 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 (m0) and a secondary (m2), and a massless planet (m1 = 0) orbiting the primary star. The binary system’s fixed orbit with period T2, semimajor axis a2 and eccentricity e2, is the natural choice of reference frame. The planet’s osculating orbit with period T1, semimajor axis a1 and eccentricity e1 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 a1/a2, and average with respect to the fast periods T1 and T2, obtaining the secular quadrupole Hamiltonian (7)with1(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 (e10 = 0.48), but we do not know the initial i0 or ω0. We can obtain possible secular trajectories for the planet by choosing values of i0 and ω0, then plotting the level curves of (12)or (13)with .

thumbnail Fig. 3

MEGNO chaos indicators for different sets of variables, a) (), b) (I1,ΔΩ), c) () with , and d) () with . Labeled lines give the mass of the planet in MJup. 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 Lidov-Kozai equilibrium points (with i and e1 fixed) and Lidov-Kozai cycles which owing to the conservation of h (Eq. (11)) exhibit coupled oscillations in i and e1. An orbit with e1 = 0.48 is at a Lidov-Kozai equilibrium point if , , or . Lidov-Kozai cycles with librating have h < 0.6 and F < Fs (magenta, brown, cyan, yellow, and pink orbits). Lidov-Kozai cycles with circulating have h > 0.6 (red, blue, and green orbits) or h < 0.6 and F > Fs (gray, orange, and purple orbits). Moreover, high-amplitude Lidov-Kozai cycles (gray, pink and orange orbits) can reach e1 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 e1 = 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 Lidov-Kozai diagrams that show level curves of the Hamiltonian F at constant values h = h0. In our case e10 is fixed, but we show trajectories for different values of the initial inclination i0, i.e., with different h0. In particular, this explains why librating orbits in the Lidov-Kozai regime do not encircle a Lidov-Kozai equilibrium point (magenta, brown, cyan, yellow, and pink orbits in Fig. 2). The location of the Lidov-Kozai equilibrium points in our diagrams occur at the current planet’s eccentricity (e1 = 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 Lidov-Kozai 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., Lidov-Kozai 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 higher-order 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 e1 ≈ 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 long-term 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 Burlisch-Stoer based in an N-body 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 (I1,ΔΩ), 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 MJup). As we change the values of () we obtain corresponding value of I1 (Eq. (1)) and consequently the real mass of the planet ranges between 2.98 MJup (dashed lines) to 40 MJup (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 m1 ∝ 1/sinI1, but the dynamical regimes are not perceptible, except instability corresponding to high masses for I1 = 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 I2 of the stellar companion to the plane of the sky (Eqs. (1), (2)). For I2 = 90° there is a perfect correspondence between the frames () and (), while for I2 = 0° or 180°, there is a global degeneracy. For the HD 196885 system, we have I2 = 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 RSun). 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 e1 = 1 (collision orbit).

  • 3.

    Near the edge of the forbidden region: (i ~ 60°,ω < 180°) and (i ~ 118°,ω < 180°), the planets have masses  > 40   Mjup(i1 ~ 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 left-hand Kozai separatrix (55 < i < 60° and 45° ≲ ω ≲ 125°) where the collision times spans 105 to 107 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 107 to 108 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.

Long-term numerical simulations showed that the regions marked in blue/green in the MEGNO map do not present detectable dynamical instability, at least on time-scales on the order of 109 years. Owing to the high mass ratio, the effect of high-order mean-motion resonances (MMR, hereafter) may be non-negligible (the period ratio is 20.95). Indeed, there are a multitude of high-order MMR whose overlap could be responsible for the slow chaos regions. However, these high-order 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.

thumbnail Fig. 4

Amplitude of the mutual inclination for different sets of variables, a) (), b) (I1,ΔΩ), c) (). Centers of the Lidov-Kozai 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 Lidov-Kozai 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 e1 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, e1, attained during the integration.

thumbnail Fig. 5

Dynamical indicators in the frame (). a) period associated to highest amplitude oscillation of e1 (in years), b) maximal amplitude of oscillation of e1, c) maximum eccentricity e1.

thumbnail 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 high-inclination orbits (i.e. Lidov-Kozai cycles) is well-described 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 low-to-medium 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 higher-order 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, e1, variations. From Fig. 5a we see that retrograde orbits have secular periods longer than direct orbits. Among stable areas, medium-amplitude Lidov-Kozai 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 Lidov-Kozai cycles occur on a 1500 yr timescale, while nearly coplanar orbits also exhibit small oscillations in e1 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 Lidov-Kozai 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 medium-amplitude Lidov-Kozai 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. − e1,max), where e1,max is the maximal eccentricity attained by the planet during the integration. Thus a value of log 10(1. − e1,max) =  −3 (in red) represents an eccentricity of e1,max = 0.999. Light-blue regions are those with e1,max ~ 0.9, dark-blue for those with e1,max ~ 0.7, and finally light-gray for those where e1,max ~ 0.5. The initial conditions within violet region inside the Kozai-Lidov separatrix can survive for times from 108 to 109 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 best-fit 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 a2 = 21.86 AU and e2 = 0.45, and at bottom ones a2 = 20.14 AU and e2 = 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 m1 = 0) using the best-fit parameters from Table 1. The entire chaotic green regions disappeared and we were left only with regular regions (blue) and close-encounter 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 k2 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 time-scale 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, e1,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 time-scale for different initial maximum eccentricities, adopting k2 = 0.5 and R = 1.2   RJup. We observe that for Q ~ 1000, a value similar to Jupiter, only for e1,max > 0.96 the eccentricity is damped within 2 Gyr, the estimated age of the system (Correia et al. 2008). Even for Earth-like planets (Q ~ 10), only e1,max > 0.92 can be dissipated in the same amount of time. Since e1 > 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).

thumbnail Fig. 7

Damping time for the eccentricity of the inner planet (with τ = 1/K, Eq. (17)). For Q ~ 1000, only for e1,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 Lidov-Kozai cycles implies destructive, high-velocity collisions among planetesimals, inhibiting formation of massive objects (Lissauer 1993). However, the Lidov-Kozai 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 self-gravity) 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 planetesimal-shattering 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 semimajor-axis 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 (abinary ~ 1000 AU), concluding that is possible to form a single planet in an inclined orbit (Lidov-Kozai regime), if taking into account the self-gravity of proto-planetary 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 Lidov-Kozai equilibrium is necessarily non-unique and proposed a scenario were a multiple system might be formed protected by Lidov-Kozai 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 first-order approximations (second-order 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 a1/a2 is a factor in the quadrupolar Hamiltonian (Eq. (8)), therefore it only modifies the time-scale of the evolution. However, for the eccentricity we observe that all trajectories that begin inside the libration zone migrate into one of the Lidov-Kozai 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 Lidov-Kozai equilibria.

thumbnail 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 Lidov-Kozai 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 insight-full analysis of the 3-D 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 Lidov-Kozai 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 case3. The orbits at the Lidov-Kozai 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 Lidov-Kozai 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 Lidov-Kozai 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°;
I1 ~ 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 Lidov-Kozai 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. Small-amplitude 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.

1

The expression for C (Eq. (7)) in Kinoshita & Nakai (2007) should have and not in the nominator.

2

The solution to Eq. (14) is independent of h, as we would expect, since the separatrix is an invariant curve.

3

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 FCT-Portugal (grants PTDC/CTE-AST/098528/2008 and PEst-C/CTM/LA0025/2011).

References

  1. Batygin, K., Morbidelli, A., & Tsiganis, K. 2011, A&A, 533, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Beaugé, C., Michtchenko, T. A., & Ferraz-Mello, S. 2006, MNRAS, 365, 1160 [Google Scholar]
  3. Beaugé, C., Leiva, A. M., Haghighipour, N., & Otto, J. C. 2010, MNRAS, 408, 503 [NASA ADS] [CrossRef] [Google Scholar]
  4. Chauvin, G., Lagrange, A. M., Udry, S., & Mayor, M. 2007, A&A, 475, 723 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Chauvin, G., Beust, H., Lagrange, A.-M., & Eggenberger, A. 2011, A&A, 528, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Cincotta, P. M., & Simó, C. 2000, A&AS, 147, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Correia, A. C. M. 2009, ApJ, 704, L1 [NASA ADS] [CrossRef] [Google Scholar]
  8. Correia, A. C. M., & Laskar, J. 2010, in Exoplanets (University of Arizona Press), 534 [Google Scholar]
  9. Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celest. Mech. Dyn. Astron., 53 [Google Scholar]
  10. Correia, A. C. M., Udry, S., Mayor, M., et al. 2008, A&A, 479, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Desidera, S., & Barbieri, M. 2007, A&A, 462, 345 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Farago, F., & Laskar, J. 2010, MNRAS, 401, 1189 [NASA ADS] [CrossRef] [Google Scholar]
  13. Fischer, D., Driscoll, P., Isaacson, H., et al. 2009, ApJ, 703, 1545 [NASA ADS] [CrossRef] [Google Scholar]
  14. Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385 [Google Scholar]
  15. Giuppone, C. A., Leiva, A. M., Correa-Otto, J., & Beaugé, C. 2011, A&A, 530, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Hatzes, A. P., Cochran, W. D., Endl, M., et al. 2003, ApJ, 599, 1383 [NASA ADS] [CrossRef] [Google Scholar]
  17. Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [NASA ADS] [CrossRef] [Google Scholar]
  18. Kinoshita, H., & Nakai, H. 1999, Celest. Mech. Dyn. Astron., 75, 125 [NASA ADS] [CrossRef] [Google Scholar]
  19. Kinoshita, H., & Nakai, H. 2007, Celest. Mech. Dyn. Astron., 98, 67 [Google Scholar]
  20. Kozai, Y. 1962, AJ, 67, 591 [Google Scholar]
  21. Krymolowski, Y., & Mazeh, T. 1999, MNRAS, 304, 720 [NASA ADS] [CrossRef] [Google Scholar]
  22. Lagrange, A.-M., Beust, H., Udry, S., Chauvin, G., & Mayor, M. 2006, A&A, 459, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Laskar, J., & Boué, G. 2010, A&A, 522, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Lee, M. H., & Peale, S. J. 2003, ApJ, 592, 1201 [NASA ADS] [CrossRef] [Google Scholar]
  25. Lee, M. H., Peale, S. J., Pfahl, E., & Ward, W. R. 2007, Icarus, 190, 103 [NASA ADS] [CrossRef] [Google Scholar]
  26. Lidov, M. L. 1961, Iskus. sputniky Zemly (in Russian), 8, 5 [Google Scholar]
  27. Lidov, M. L. 1962, Plan. Space Sci., 9, 719 [NASA ADS] [CrossRef] [Google Scholar]
  28. Lissauer, J. J. 1993, ARA&A, 31, 129 [Google Scholar]
  29. Lithwick, Y., & Naoz, S. 2011, ApJ, 742, 94 [NASA ADS] [CrossRef] [Google Scholar]
  30. Maffione, N. P., Darriba, L. A., Cincotta, P. M., & Giordano, C. M. 2011, Celest. Mech. Dyn. Astron., 111, 285 [NASA ADS] [CrossRef] [Google Scholar]
  31. Morais, M. H. M., & Giuppone, C. A. 2012, MNRAS, in press [arXiv:1204.4718] [Google Scholar]
  32. Mugrauer, M., & Neuhäuser, R. 2009, A&A, 494, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Neuhäuser, R., Mugrauer, M., Fukagawa, M., Torres, G., & Schmidt, T. 2007, A&A, 462, 777 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Queloz, D., Mayor, M., Weber, L., et al. 2000, A&A, 354, 99 [NASA ADS] [Google Scholar]
  35. Quintana, E. V., Lissauer, J. J., Chambers, J. E., & Duncan, M. J. 2002, ApJ, 576, 982 [NASA ADS] [CrossRef] [Google Scholar]
  36. Rabl, G., & Dvorak, R. 1988, A&A, 191, 385 [NASA ADS] [Google Scholar]
  37. Smart, W. M. 1965, Text-book on spherical astronomy (Cambridge University Press) [Google Scholar]
  38. Thebault, P. 2011, Celest. Mech. Dyn. Astron., 111, 29 [NASA ADS] [CrossRef] [Google Scholar]
  39. Torres, G. 2007, ApJ, 654, 1095 [NASA ADS] [CrossRef] [Google Scholar]
  40. Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  41. Wiegert, P. A., & Holman, M. J. 1997, AJ, 113, 1445 [NASA ADS] [CrossRef] [Google Scholar]
  42. 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

Table 1

Orbital parameters for the HD 196885 system obtained by Chauvin et al. (2011) for JD = 2 455 198.

In the text

All Figures

thumbnail 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 (I1,ΔΩ) 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
thumbnail 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 n-body 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 Lidov-Kozai equilibrium points for this system. N-body 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
thumbnail Fig. 3

MEGNO chaos indicators for different sets of variables, a) (), b) (I1,ΔΩ), c) () with , and d) () with . Labeled lines give the mass of the planet in MJup. Orbits can be considered stable for  ⟨ Y ⟩  ≤ 2.
In the text
thumbnail Fig. 4

Amplitude of the mutual inclination for different sets of variables, a) (), b) (I1,ΔΩ), c) (). Centers of the Lidov-Kozai equilibria are marked as black circles and the separatrix is calculated with the quadrupolar approximation (Sect. 3).
In the text
thumbnail Fig. 5

Dynamical indicators in the frame (). a) period associated to highest amplitude oscillation of e1 (in years), b) maximal amplitude of oscillation of e1, c) maximum eccentricity e1.
In the text
thumbnail 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
thumbnail Fig. 7

Damping time for the eccentricity of the inner planet (with τ = 1/K, Eq. (17)). For Q ~ 1000, only for e1,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
thumbnail 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 Lidov-Kozai equilibria.
In the text

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

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

Initial download of the metrics may take a while.