LidovKozai stability regions in the α Centauri system
^{1} Universidad Nacional de Córdoba, Observatorio Astronómico, IATE, Laprida 854, 5000 Córdoba, Argentina
email: cristian@oac.unc.edu.ar
^{2} CIDMA, Departamento de Física, Universidade de Aveiro, Campus de Santiago, 3810193 Aveiro, Portugal
^{3} ASD, IMCCECNRS UMR 8028, Observatoire de Paris, 77 Av. DenfertRochereau, 75014 Paris, France
Received: 24 January 2017
Accepted: 26 June 2017
The stability of planets in the α Centauri AB stellar system has been studied extensively. However, most studies either focus on the orbital plane of the binary or consider inclined circular orbits. Here, we numerically investigate the stability of a possible planet in the α Centauri AB binary system for Stype orbits in an arbitrary spatial configuration. In particular, we focus on inclined orbits and explore the stability for different eccentricities and orientation angles. We show that large stable and regular regions are present for very eccentric and inclined orbits, corresponding to libration in the LidovKozai resonance. We additionally show that these extreme orbits can survive over the age of the system, despite the effect of tides. Our results remain qualitatively the same for any compact binary system.
Key words: binaries: close / methods: numerical / celestial mechanics / planetary systems
© ESO, 2017
1. Introduction
The nearest neighbour to our solar system, α Centauri stellar system, always captured the attention of diverse studies in astronomy. The compact binary system α Cen A and α Cen B has an orbital period of 79 yr and has a very eccentric orbit, thus challenging the formation theories of potential existing planets. Planetary accretion models suggest that Stype planets could have formed within the α Cen system (Quintana et al. 2002, 2007) provided that the collision velocities of latestage planetesimals are not too large (Thébault et al. 2008, 2009; Thebault & Haghighipour 2014). However, the accretional collisions that form planets in compact binary systems is a complicate mechanism, which depends on the initial conditions of the particles and on the mass and orbital parameters of the secondary star (e.g. Beaugé et al. 2010).
The announcement of an Earthmass planet candidate in a 3.24 day orbit α Cen B b (Dumusque et al. 2012) and the tentative detection of a transiting planet on a more distant orbit by Demory et al. (2015) put the spotlight again on this system. Nonetheless, Hatzes (2013) concluded that the presence of the activity signal from the star may boost the velocity amplitude to values comparable to the planet signature. Recently, a groundbased radial velocity campaign has ruled out the presence of massive closein planets (Endl et al. 2015), although a Jupiter (or less massive) distant planet may exist but has not been detected because the time span of observations is not long enough yet.
Early works from Benest (1988) and Wiegert & Holman (1997) examined the stability of planetary orbits in the α Centauri system with the logical CPU limitations at the time. More recently, AndradeInes & Michtchenko (2014) and Quarles & Lissauer (2016) have studied the stability regions of this system again. However, all these studies focused on nearly coplanar orbits or inclined circular orbits. Moreover, they have chosen very particular orientation angles for the orbit of the planet, which result in a reduced exploration of the phase space.
In this paper, we study the stability of an additional planet covering all the orbital parameters. This allows us to uncover previously unnoticed stability regions. We describe the methods in Sect. 2. Results are presented in Sects. 3–5, and comparisons with other binary systems are discussed in Sect. 6. Finally, our conclusions are presented in Sect. 7
2. Methods
The orbit of the α Cen AB binary is fully constrained by astrometry. Owing to the previous observational constraints, we assumed that if the system is hosting a planet, its mass should be lower than the mass of Jupiter (Endl et al. 2015). Since the two binary stars have similar masses, we studied the motion around the more massive star (α Cen A) using the orbital plane of the binary as a reference frame.
Semimajor axis, eccentricity, and masses for the binary systems studied in this paper.
Each planetary orbit can be described by six orbital elements: three “actions”, i.e. the semimajor axis a, eccentricity e, and inclination J, with respect to the binary orbital plane; and three “conjugated” angles, the mean longitude λ, longitude of pericentre ϖ, and longitude of the ascending node Ω, respectively, measured from the direction of the pericentre of the binary.
Previous works have always focused on the action variables (a,e,J) and these works arbitrarily fix the angles, usually at zero (λ = 0,ϖ = 0,Ω = 0). This strategy is understandable. Since the amplitude of the interactions during close encounters depends on action variables, the main features are captured this way. However, for some particular choices of the angles it is also possible to avoid close encounters, namely when they are involved in resonances (mean motion or secular). A full exploration of the phase space thus requires the inspection of its six free parameters. In particular, we need to explore the pairs of conjugated elements, i.e. (a,λ) to identify regions with meanmotion resonances; and (e,ϖ) or (J,Ω) to study secular resonances.
We solved the threebody equations of motion numerically with a BurlischStoer integrator with double precision and tolerance 10^{12}. We stopped the integrations when the distance of the planet to one of the stars is lower than one stellar radius or when the planet is ejected from the system, identifying the time for any of these situations as “disruption times”. The integration time is 2 × 10^{5} yr, which corresponds to several periods of the secular variations.
We analysed the stability of a test planet in a Stype orbit (a<a_{B}) for a wide variety of configurations. We constructed stability maps integrating the system on a regular 2D mesh of initial conditions for a pair of orbital parameters, while the remaining four parameters are initially set at nominal values. Unstable orbits during the integration time, which either collide with one of the stars or escape from the system, are identified in the figures in white.
For each initial condition we computed the mean exponential growth of nearby orbits (MEGNO) value ⟨ Y ⟩ because this value can identify chaotic orbits in less CPU time than other indicators (Cincotta & Simó 2000; Maffione et al. 2011). We compared this value with other chaos indicators (Lyapunov exponent and frequency analysis) and the chaotic regions coincide. However, MEGNO cannot give a precise representation of the structure of a resonance, as it only differentiates regular (⟨ Y ⟩ ~ 2, blue regions) from chaotic orbits (⟨ Y ⟩ ≫ 2, red regions).
To study the structure of the secular resonances, we used the amplitude of maximum variation of the eccentricity of the planet attained during the integrations, (1)The Δe indicator is an extremely useful tool to map the resonant structure in Nbody problems (Ramos et al. 2015); however Δe is not a measure of chaotic motion. Abrupt changes in Δe are often traces for the presence of resonances, while regions with large variations in Δe are more sensitive to perturbations, thus are very likely chaotic (e.g. Giuppone et al. 2012; Martí et al. 2013).
In Fig. 1 we compare the two indicators used in this work for different values of the initial mutual inclinations. We plot results for ⟨ Y ⟩ and Δe values integrated over 10^{5} yr. For visual representation, the MEGNO unstable orbits are identified at the top of the scale of ⟨ Y ⟩. We can see that the minimum variation of Δe corresponds to regions where the orbits are regular ⟨ Y ⟩ ~ 2. There are some regions with chaotic orbits (2 < ⟨ Y ⟩ < 16), that we checked to remain stable over 10^{8} yr with almost the same value of Δe. In Fig. 2 we show the evolution of two bounded chaotic orbits with the initial conditions at J = 60° and J = 63°. They undergo slow diffusion in the orbital elements and still survived over 10^{9} yr. Thus, the MEGNO values larger than 2 are somehow related to this diffusion and not necessarily to unstability, i.e. “stable chaos” as defined by Milani & Nobili (1992). Large MEGNO values give us an estimation of the longterm stability, while Δe measures the amplitude of the orbital secular variations.
Fig. 1 Stability indicators for different initial mutual inclination with a = 2 au, e = 0.9, ω = 90°, and λ = Ω = 0. The left scale corresponds to the MEGNO chaos indicator, ⟨ Y ⟩ (red crosses), while the right scale corresponds to the Δe indicator (green circles). 

Open with DEXTER 
Fig. 2 Evolution of initial conditions from Fig. 1 at J = 60° (left) and J = 63° (right). Slow diffusion is present in the orbital elements. MEGNO values are proportional to such diffusion (bottom), while Δe measures the maximal amplitude eccentricity variations (top). Both orbits survived for 10^{9} yr. 

Open with DEXTER 
3. Actionangle maps
3.1. Meanmotion resonances
We begin our quest for stability regions by varying the pair (a,λ), which allow us to identify the presence of possible meanmotion resonances. We consider coplanar prograde (J = 0°) and retrograde orbits (J = 180°). We additionally set ϖ = 0 and Ω = 0, which correspond to orbits with aligned pericentres, as these orbits are among the most favourable for meanmotion resonances to occur.
Fig. 3 Stability maps of ⟨ Y ⟩ in the plane (a,λ) for some eccentricity values with Ω = ω = 0°, J = 0° (left) and J = 180° (right). From top to bottom the initial e = 0,0.3,0.6, and 0.9. Vertical lines indicate the main meanmotion resonances. 

Open with DEXTER 
In Fig. 3 we show the stability maps for different values of the initial eccentricity, e = 0,0.3,0.6 and 0.9. We extend the semimajor axis up to a = 6 au, since the Hill radius of α Cen A is ~ 6.4 au (e.g. Marchal & Bozis 1982). The initial conditions with MEGNO values ⟨ Y ⟩ ~ 2 are regular, orbits with 2 < ⟨ Y ⟩ < 16 show very small diffusion in orbital elements, and orbits with ⟨ Y ⟩ ≳ 16, identified as red regions in the dynamical maps, show high diffusion and eventually collide with one of the stars.
As in all previous studies, we observe that in the prograde case, stable orbits are only possible for a ≲ 3 au (Holman & Wiegert 1999). Stability for 3 ≲ a ≲ 6 au is not possible because of resonance overlap, which leads to chaotic motions in these regions (Wisdom 1980). Although stability slightly depends on the initial λ value, we see that it is not possible to trap a planet in a low order meanmotion resonance with the companion star. However, for small eccentricities we can observe some small resonant islands for the 15:1 and the 16:1 meanmotion resonances.
Marzari & Gallina (2016) have used frequency map analysis to study the stability of a planet in a binary system with e_{B} = 0 and e_{B} = 0.4. These authors have shown that in the case of an eccentric binary, for a given semimajor axis the orbit of the planet can be regular or chaotic, depending on the initial mean longitude, λ. However, their initial λ values were chosen randomly, so it was not possible to understand the origin of this behavior exactly. In Fig. 3 we clearly see that for λ close to 0° and 180° stable resonant islands exist at the middle of chaotic regions, which allow different stability regimes for the same semimajor axis value.
For retrograde orbits we observe that stability is possible for larger values of the semimajor axis. In particular, for e = 0.3 stability is possible up to 6 au, very close to the Hill sphere. We additionally observe that capture in lower order meanmotion resonances, such as 7:1 or 6:1, is also possible. Indeed, retrograde orbits in binary systems are more stable than the prograde orbits because of a different structure of meanmotion resonance overlaps (see Morais & Giuppone 2012). We hence conclude that the inclination value is a very important parameter that shapes the stability in binary systems.
3.2. Secular resonances
Fig. 4 Stability maps of ⟨ Y ⟩ in the plane (ω,J) for some semimajor axis and eccentricity values with λ = Ω = 0. From left to right the initial a = 2 au, a = 2.42 au, and a = 3 au. From top to bottom the initial e = 0, 0.3, 0.6, and 0.9. 

Open with DEXTER 
The main frequency involved in the angle λ is the orbital mean motion, n, thus this angle varies rapidly. The angles Ω and ϖ vary in a much longer timescale owing to the presence of the binary companion, hence the name secular. Since we are studying a threebody problem, owing to the conservation of the total angular momentum there is a single frequency associated with the precession of the line of the nodes, . Moreover, since most of the angular momentum is on the binary orbit, the precession frequency associated with the pericentre of the binary is almost zero.
Thus, the longitude of the pericentre is mainly driven by a single frequency , which corresponds to the precession rate of the pericentre of the planet. Indeed, in the restricted problem, g and s are the only secular frequencies in the system. As a consequence, secular resonances can only occur when g = s, which is usually known in the literature by the LidovKozai resonance (Lidov 1962; Kozai 1962). The particular geometry of this almost restricted threebody problem allows us to explore the phasespace more rapidly. Instead of using the angles Ω and ϖ separately, we can adopt the argument of the pericentre ω = ϖ − Ω, for which resonances occur when . Therefore, all the significant information on secular resonances can be captured by a (ω,e) or (ω,J) diagram.
In Fig. 4 we show the stability maps in the plane (ω,J) for different values of the initial eccentricity (e = 0, 0.3, 0.6 and 0.9), and three different values of the semimajor axis, corresponding to three different stability regions: a = 2 au, which places the planet inside a stable region for prograde orbits; a = 2.42 au, corresponding to the 22:1 meanmotion resonance and near the unstable region; and a = 3 au, already in a chaotic region for prograde orbits (see Fig. 3).
In Fig. 4 we also observe that polar orbits (J ~ 90°) are always unstable. Owing to the conservation of the orbital angular momentum, the following quantity is conserved:(2)As a consequence, the eccentricity of polar orbits can reach values very close to unity, which may place the planet outside the Hill sphere or in collision with the star.
For initial circular orbits (e = 0) stability is only possible for nearly coplanar orbits (J ≲ 30° or J ≳ 150°). Prograde orbits (J ≲ 30°) become unstable for a ~ 3 au, in conformity with the results shown in Fig. 3, except for a small zone of aligned and antialigned orbits (ω ≈ 0° and ω ≈ 180°, respectively). However, as also shown in Fig. 3, the stability region is extended beyond this value of the semimajor axis for retrograde orbits (J ≳ 150°). Some additional chaotic structures can also be seen for some ω values, probably due to secondary nonlinear resonances.
It is often assumed that circular orbits provide an upper limit for stability with a given semimajor axis, since the minimal distance between the planet and the perturber decreases with the eccentricity. However, we observe that as we increase the initial eccentricity, the coplanar regions become indeed less stable, but new stable islands emerge in the region within 30° ≲ J ≲ 150°. An interesting result is that for a> 2 au stable prograde coplanar orbits are no more possible for moderate eccentricities, but stability can still be achieved in these islands for very high values of eccentricity and mutual inclination (Fig. 4).
The stability islands at high inclinations are centred at ω = 90° and ω = 270° and correspond to the secular LidovKozai resonances. For the restricted problem, the resonant motion is possible whenever (Lidov 1962; Kozai 1962) (3)We then conclude that resonant motion is only possible for J_{c} ≤ J ≤ π − J_{c}, with the critical inclination (corresponding to zero eccentricity). The equality above also corresponds to exact resonance. For a given initial eccentricity e, the equilibrium inclination J_{r} at exact resonance is then given by(4)For instance, for e = 0.6 we have J_{r} = 51.7° or J_{r} = 128.3° and for e = 0.9 we have J_{r} = 70.3° or J_{r} = 109.7°. For trajectories in libration around this equilibrium, the minimum value of the inclination for prograde orbits lies in the interval J_{c} ≤ J<J_{r} and the maximum inclination in retrograde orbits lies in π − J_{c} ≥ J>π − J_{r}. Since for initial circular orbits (e = 0) we have J_{c} = J_{r}, the entire region occupied by the LidovKozai resonance is unstable. However, as we increase the initial eccentricity, the equilibrium mutual inclination for prograde (retrograde) orbits moves to higher (lower) values and stable resonant regions appear in the vicinity of ω = 90° and ω = 270°.
In Fig. 1 we show the value of the chaotic indicators for a vertical line with ω = 90° in the bottom left panel of Fig. 4 (a = 2 au and e = 0.9). We observe that the most regular region is obtained for the libration regions of the LidovKozai resonance. Regular coplanar retrograde orbits are still possible in this case, but they already present eccentricities very close to instability.
4. New actionaction maps
Maps involving the three actions (a,e,J) were extensively explored in previous studies (e.g. AndradeInes & Michtchenko 2014; Quarles & Lissauer 2016). However, they usually fix ω = 0°, which corresponds to a region of the phasespace that is always outside the libration zone of the LidovKozai resonance (Fig. 4). As we just saw in previous section, stable regions for high eccentricity and mutual inclination are near the centre of libration, which is placed at ω = 90° or ω = 270°. Therefore, it is better to fix ω at one of these two values to capture the resonant regions in actionaction maps.
Fig. 5 Stability maps of ⟨ Y ⟩ in the plane (a,J) for some eccentricity values with λ = Ω = 0° and ω = 0° (left) or ω = 90° (right). From top to bottom the initial eccentricity is e = 0, 0.3, 0.6, and 0.9. The horizontal grey lines give the centre of libration of the LidovKozai resonance (Eq. (4)). 

Open with DEXTER 
Fig. 6 Stability maps of ⟨ Y ⟩ in the plane (a,J) for e = 0.3 and ω = 0° with λ = Ω = 0. Left panel: prograde orbits and right panel: retrograde orbits. The vertical labels indicate the nominal position of the N:1 meanmotion resonances. 

Open with DEXTER 
4.1. Semimajor axis versus mutual inclination
In Fig. 5 we show the stability maps in the plane (a,J) for different values of the initial eccentricity e = 0, 0.3, 0.6, and 0.9. We fix λ = Ω = 0, and ω = 0° (left column) or ω = 90° (right column) to compare better with previous studies. We can see that for both ω values the region around polar orbits is very chaotic and splits the regions corresponding to prograde (J< 90°) and retrograde (J> 90°) orbits. Generally, collisions and/or ejections occur in less than 5 × 10^{4} yr for nearby orbits. We scale our integrations time such that more distant orbits are integrated at least 5 × 10^{4} planetary periods.
For initial circular orbits (e = 0) there is no big difference between the two ω values. Indeed, in the restricted quadrupolar problem, the Hamiltonian only depends on the product e^{2}cos2ω (e.g. Kozai 1962; Giuppone et al. 2012), so the initial value of ω does not change the energy at the quadrupole order, which is the dominating term. We also observe that retrograde coplanar orbits are stable for larger values of semimajor axis than prograde coplanar orbits (as also noted by Holman & Wiegert 1999; Quarles & Lissauer 2016). Indeed, for compact binary systems the retrograde planets are stable up to distances closer to the perturber than prograde planets owing to the higher order of overlap of nearby resonances (Morais & Giuppone 2012).
As we increase the initial eccentricity, the coplanar orbits remain more stable in the case ω = 0°, as they correspond to aligned orbits (see Giuppone et al. 2013). However, for inclined orbits we observe that a new stability region appears for maps with ω = 90°, in a strip for mutual inclinations given by expression (4), corresponding to libration in the LidovKozai resonance. We also observe a chaotic strip clearly delimiting the coplanar and resonant regions corresponding to the separatrix of this resonance.
In Fig. 5 we see that the eccentricity does not limit the stability in binary systems, provided that we change the mutual inclination following the resonant equilibrium points (Eq. (4)). Moderate eccentricities can also facilitate stability for coplanar orbits. Indeed, for e = 0.3 and ω = 0° (aligned orbits) we observe that stable regions exist beyond 3 au and 6 au for prograde and retrograde orbits, respectively. In Fig. 6 we zoom in on the coplanar regions for these initial conditions. We superimposed the nominal location of the N:1 meanmotion resonances to better understand the structures present in these regions.
We observe that stability islands are associated with meanmotion resonances between the planet and the stellar companion α Cen B. The last stable resonances correspond to the 15:1 for prograde orbits and 6:1 for retrograde orbits. Lower order resonances are not possible because they lie outside the Hill sphere of α Cen A. The only real limitation for stability is thus the semimajor axis; for the remaining orbital parameters, the stability can always be achieved at some particular combinations.
4.2. Eccentricity versus mutual inclination
In Fig. 7 we show the stability maps in the plane (e,J) for different values of the semimajor axis a = 0.55, 1.5, 2.0, 3.0, and 4.0 au. We fix λ = Ω = 0 and ω = 90° such that the LidovKozai resonance is visible. We also plot the curve corresponding to the centre of this resonance, given by expression (4) obtained in the frame of the restricted quadrupolar approximation.
Fig. 7 Stability maps in the plane (e,J) with λ = Ω = 0 and ω = 90° for initial semimajor axis from left to right: a = 0.55, 1.5, 2.0, 3.0, and 4.0 au. We show the MEGNO ⟨ Y ⟩ (top) and the Δe (bottom) stability indicators. The grey curve gives the centre of libration of the LidovKozai resonance (Eq. (4)), while the red vertical line gives the limit for tidal stability (Eq. (8)). 

Open with DEXTER 
We observe that for small values of the semimajor axis (a< 2 au), almost all configurations are stable, except for those in near polar orbits. However, for the semimajor axis of roughly 3 au, most configurations become unstable. As in previous figures, we see that only a small retrograde zone subsists for small eccentricities, together with the LidovKozai regions. In particular, there is a perfect agreement between the theoretical prediction given by expression (4) and the stable regions with high inclination.
In Fig. 7 (bottom) we also show the Δe stability indicator (Eq. (1)) to get a clearer view of the secular dynamics in the α Cen system. We observe that the libration resonant areas are the most stable structures in the system. Libration regions are present for closein orbits (a< 1 au) together with the stable coplanar regions, but the former subsist for more distant semimajor axes (a ~ 3 au), while the coplanar regions are no longer stable. The stability regions are also larger for high values of eccentricity and inclination, since the libration zone is more extended.
5. Tidal evolution
For some orbital configurations, the eccentricity of the planet may reach very high values and become close enough to the central star at periastron to undergo tidal effects. In that case, the semimajor axis and the eccentricity will decrease and the final configuration can be completely different from the initial configuration.
Fig. 8 Tidal evolution in the plane (e,J) obtained with a secular model for the same initial conditions from Fig. 7. The colour index gives the variation between the initial and final semimajor axes, log _{10}(Δa/a), and white stands for unstable orbits. The initial semimajor axes from left to right are a = 0.55, 1.5, and 2.0 au. The black vertical lines give the tidal stability limit (Eq. (8)). 

Open with DEXTER 
For an unperturbed orbit, the secular evolution of the eccentricity by tidal effect using a linear dissipation model can be given by (Correia 2009) (5)with (6)and (7)where is the initial mean motion, k_{2} is the second Love number, Q is the tidal dissipation factor, m is the mass of the planet, and R its radius.
The solution of the above equation is given by (Correia & Laskar 2010) (8)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_{0}. Orbits with τ smaller than the age of the system can be excluded from the stability diagrams since the planet will not stay at the original semimajor axis because the orbit of the planet evolves. This does not mean that the planet is necessarily destroyed, only that it evolves into a different value of the semimajor axis and/or eccentricity.
The time τ depends on many uncertain parameters, so it is not easy to place a clear limit for tidal stability. In particular, τ should be different for rocky and gaseous planets, since rocky bodies usually dissipate energy more efficiently. Indeed, rocky planets in the solar system present k_{2}/Q ~ 10^{2} − 10^{3}, while for gaseous planets k_{2}/Q ~ 10^{4} − 10^{5} (Yoder 1995). However, the mass and radius of a gaseous planet is in general 10^{2} and 10 times larger than mass and the radius of a rocky planet, respectively. When replacing all these values in expression (7) we get similar values for τ for both types of planets.
In Fig. 7 we trace a vertical red line corresponding to the solution of Eq. (8) for a timescale τ smaller than 1 Gyr assuming a Jupiterlike planet with k_{2}/Q = 1.1 × 10^{5} (Lainey et al. 2009); the solutions with higher initial eccentricities should not be considered, since they evolve in a period of time shorter than the age of the system. As expected, we observe that for smaller values of semimajor axis we have to exclude more configurations, because tides are stronger and the orbits evolve faster. However, for a = 1.5 au, we only need to exclude orbits with e> 0.96. For the LidovKozai resonance, these eccentricities correspond to an equilibrium mutual inclination 78°<J< 102°, which is also unstable in the absence of tides. We thus conclude that tidal effects only need to be taken into account for closein planets (a< 1.5 au). In particular, they do not disturb the orbits at the edge of stability (a> 2 au). Therefore, the stability islands observed at high eccentricities and inclinations remain a possibility to find planets in closein binaries.
The above equations are only valid for unperturbed orbits, but the initial eccentricity can be seen as the maximal eccentricity over a cycle, so τ provides a minimal estimation of the dampening time. Indeed, for orbits inside the LidovKozai resonant region, the orbital damping may drive the planet into the exact resonance (Giuppone et al. 2012), so the eccentricity can stay very high. A complete analysis requires a study that combines tidal effects with orbital forcing. Adopting the secular tidal model^{1} from Correia et al. (2016), we have run simulations for three different semimajor axes with tides using the same initial conditions from Fig. 7. In Fig. 8 we show the final evolution of the semimajor axis after 5 Gyr. The results show that the theoretical estimation given by expression (8) is relatively accurate and can be used to put constraints on the tidal evolution. In addition, Fig. 8 also shows that the main chaotic structures are captured by the secular octupolar model.
6. Other compact binary systems
Fig. 9 Stability maps of ⟨ Y ⟩ in the plane (a/a_{B},J) for some eccentricity values, with λ = Ω = 0 and ω = 90° for three different compact binary systems. From left to right: α Centauri, HD 196885, and HD 41004 are shown. From top to bottom, the initial eccentricity is e = 0, 0.3, 0.6, and 0.9. 

Open with DEXTER 
We have seen that the LidovKozai resonance is an important mechanism that allows stable regions with high eccentricity and high inclination in the α Cen system. Although no planets are known for this system (Hatzes 2013; Endl et al. 2015; Rajpaul et al. 2016), other similar compact binary systems exist for which planets have been reported in very eccentric orbits such as HD 196885 b and HD 41004 Ab (see Table 1). We may then wonder how the stability regions are modified for the different mass ratios of these systems. As for α Cen, the stability in the two other binary systems has already been studied before (e.g. Funk et al. 2015), but they focus on prograde and nearly coplanar systems (J< 60°).
In Fig. 9 we show the stability maps in the plane (a,J) with ω = 90° for the three compact binary systems listed in Table 1, which include α Cen. These maps are the same as shown in Fig. 5, but the semimajor axis of the planet is normalised by the semimajor axis of the binary to enable a better comparison between the different systems. Actually, the semimajor axis of the binary does not change much between the systems, but the mass ratios are M_{B}/M_{A} = 0.85, 0.33, and 0.57 for α Cen, HD 196885, and HD 41004, respectively. We observe that the results are qualitatively the same with the LidovKozai regions located at same places. The only difference is that stability can be obtained for more distant semimajor axis ratios, since the Hill sphere of the main star is larger. We hence conclude that the stability maps drawn for α Cen are very general and can be used as reference for other compact binary systems.
7. Conclusions
In this paper we have numerically investigated the stability of Stype planetary orbits in the α Centauri system. In particular, we studied the stability on inclined orbits for high eccentricities and various orientation angles.
The Hill radius of α Cen A is ~ 6.4 au, but stability for coplanar prograde orbits can only be achieved for a< 3 au owing to meanmotion resonances overlap (Fig. 3). We have shown that nearly coplanar retrograde orbits (J> 150°) with moderate eccentricity (e ~ 0.3) can extend the stability regions beyond 6 au, very close to the limits of the Hill sphere (Fig. 5).
We have also shown that an exhaustive study of the stability regions cannot be restricted to the action variables (a,e,J). The conjugated angles are also important, in particular the argument of the pericentre ω = ϖ − Ω. For simplicity, previous studies usually set ω = 0°, but this choice limits the stability regions at very high inclinations. For 39.2°<J< 140.8° large stable regions appear located around ω = 90° and ω = 270°, corresponding to libration in the LidovKozai resonance.
As the eccentricity increases, the mutual inclination at the centre of the LidovKozai resonance also increases. The LidovKozai resonant region is thus the most stable region for planets in eccentric orbits. It persists for high inclinations, but also for semimajor axes close to 4 au. As for coplanar orbits, the retrograde regions of this resonance 100°<J< 140.8° are also more stable than the prograde regions 39.2°<J< 80°. For very eccentric orbits (e> 0.9), tidal effects can also modify the LidovKozai equilibrium, but only for closein planets (a< 1.5 au).
In this paper we focused on the stability of Stype orbits in the α Centauri system. Nevertheless, our results remain qualitatively the same for any compact binary system with significant eccentricity (e_{B}> 0.4). For binaries in nearly circular orbits, low order resonance capture is possible and the global picture may be different (Marzari & Gallina 2016).
Finally, we may wonder about the reliability of forming planets in very eccentric and inclined orbits in binary systems. Indeed, at present no planets have been found in such configurations. However, it appears to be possible to trap single planets at ~ 2 au in LidovKozai configurations in compact binary systems (a_{B} ~ 20 au) when tides are considered, through a close flyby of a background star (Martí & Beaugé 2012, 2015).
This model uses the octupolar nonrestricted approximation for the orbital interactions, general relativity corrections, the quadrupolar approximation for the spins and the viscous linear model for tides. Although in Correia et al. (2016) the authors apply their model to study Ptype circumbinary orbits, it is also valid to study Stype orbits as in (Correia et al. 2011).
Acknowledgments
The authors acknowledge financial support from CIDMA strategic project UID/MAT/04106/2013. The computations were performed at the BlaFis cluster at the University of Aveiro.
References
 AndradeInes, E., & Michtchenko, T. A. 2014, MNRAS, 444, 2167 [NASA ADS] [CrossRef] [Google Scholar]
 Beaugé, C., Leiva, A. M., Haghighipour, N., & Otto, J. C. 2010, MNRAS, 408, 503 [NASA ADS] [CrossRef] [Google Scholar]
 Benest, D. 1988, A&A, 206, 143 [NASA ADS] [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., Boué, G., & Laskar, J. 2016, Celest. Mech. Dyn. Astron., 126, 189 [NASA ADS] [CrossRef] [Google Scholar]
 Demory, B.O., Ehrenreich, D., Queloz, D., et al. 2015, MNRAS, 450, 2043 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Pepe, F., Lovis, C., et al. 2012, Nature, 491, 207 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Endl, M., Bergmann, C., Hearnshaw, J., et al. 2015, Int. J. Astrobiol., 14, 305 [CrossRef] [Google Scholar]
 Funk, B., PilatLohinger, E., & Eggl, S. 2015, MNRAS, 448, 3797 [NASA ADS] [CrossRef] [Google Scholar]
 Giuppone, C. A., Morais, M. H. M., Boué, G., & Correia, A. C. M. 2012, A&A, 541, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Giuppone, C. A., Morais, M. H. M., & Correia, A. C. M. 2013, MNRAS, 436, 3547 [NASA ADS] [CrossRef] [Google Scholar]
 Hatzes, A. P. 2013, ApJ, 770, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Lainey, V., Arlot, J.E., Karatekin, Ö., & van Hoolst, T. 2009, Nature, 459, 957 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lidov, M. L. 1962, Planet. Space Sci., 9, 719 [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]
 Marchal, C., & Bozis, G. 1982, Celest. Mech., 26, 311 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Martí, J. G., & Beaugé, C. 2012, A&A, 544, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martí, J. G., & Beaugé, C. 2015, Int. J. Astrobiol., 14, 313 [CrossRef] [Google Scholar]
 Martí, J. G., Giuppone, C. A., & Beaugé, C. 2013, MNRAS, 433, 928 [NASA ADS] [CrossRef] [Google Scholar]
 Marzari, F., & Gallina, G. 2016, A&A, 594, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Milani, A., & Nobili, A. M. 1992, Nature, 357, 569 [NASA ADS] [CrossRef] [Google Scholar]
 Morais, M. H. M., & Giuppone, C. A. 2012, MNRAS, 424, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Quarles, B., & Lissauer, J. J. 2016, AJ, 151, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Quintana, E. V., Lissauer, J. J., Chambers, J. E., & Duncan, M. J. 2002, ApJ, 576, 982 [NASA ADS] [CrossRef] [Google Scholar]
 Quintana, E. V., Adams, F. C., Lissauer, J. J., & Chambers, J. E. 2007, ApJ, 660, 807 [NASA ADS] [CrossRef] [Google Scholar]
 Rajpaul, V., Aigrain, S., & Roberts, S. 2016, MNRAS, 456, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Ramos, X. S., CorreaOtto, J. A., & Beaugé, C. 2015, Celest. Mech. Dyn. Astron., 123, 453 [NASA ADS] [CrossRef] [Google Scholar]
 Thebault, P., & Haghighipour, N. 2014, ArXiv eprints [arXiv:1406.1357] [Google Scholar]
 Thébault, P., Marzari, F., & Scholl, H. 2008, MNRAS, 388, 1528 [NASA ADS] [CrossRef] [Google Scholar]
 Thébault, P., Marzari, F., & Scholl, H. 2009, MNRAS, 393, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Wiegert, P. A., & Holman, M. J. 1997, AJ, 113, 1445 [NASA ADS] [CrossRef] [Google Scholar]
 Wisdom, J. 1980, AJ, 85, 1122 [NASA ADS] [CrossRef] [Google Scholar]
 Yoder, C. F. 1995, in Global Earth Physics: A Handbook of Physical Constants (Washington D.C: American Geophysical Union), 1–31 [Google Scholar]
 Zucker, S., Mazeh, T., Santos, N. C., Udry, S., & Mayor, M. 2004, A&A, 426, 695 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Semimajor axis, eccentricity, and masses for the binary systems studied in this paper.
All Figures
Fig. 1 Stability indicators for different initial mutual inclination with a = 2 au, e = 0.9, ω = 90°, and λ = Ω = 0. The left scale corresponds to the MEGNO chaos indicator, ⟨ Y ⟩ (red crosses), while the right scale corresponds to the Δe indicator (green circles). 

Open with DEXTER  
In the text 
Fig. 2 Evolution of initial conditions from Fig. 1 at J = 60° (left) and J = 63° (right). Slow diffusion is present in the orbital elements. MEGNO values are proportional to such diffusion (bottom), while Δe measures the maximal amplitude eccentricity variations (top). Both orbits survived for 10^{9} yr. 

Open with DEXTER  
In the text 
Fig. 3 Stability maps of ⟨ Y ⟩ in the plane (a,λ) for some eccentricity values with Ω = ω = 0°, J = 0° (left) and J = 180° (right). From top to bottom the initial e = 0,0.3,0.6, and 0.9. Vertical lines indicate the main meanmotion resonances. 

Open with DEXTER  
In the text 
Fig. 4 Stability maps of ⟨ Y ⟩ in the plane (ω,J) for some semimajor axis and eccentricity values with λ = Ω = 0. From left to right the initial a = 2 au, a = 2.42 au, and a = 3 au. From top to bottom the initial e = 0, 0.3, 0.6, and 0.9. 

Open with DEXTER  
In the text 
Fig. 5 Stability maps of ⟨ Y ⟩ in the plane (a,J) for some eccentricity values with λ = Ω = 0° and ω = 0° (left) or ω = 90° (right). From top to bottom the initial eccentricity is e = 0, 0.3, 0.6, and 0.9. The horizontal grey lines give the centre of libration of the LidovKozai resonance (Eq. (4)). 

Open with DEXTER  
In the text 
Fig. 6 Stability maps of ⟨ Y ⟩ in the plane (a,J) for e = 0.3 and ω = 0° with λ = Ω = 0. Left panel: prograde orbits and right panel: retrograde orbits. The vertical labels indicate the nominal position of the N:1 meanmotion resonances. 

Open with DEXTER  
In the text 
Fig. 7 Stability maps in the plane (e,J) with λ = Ω = 0 and ω = 90° for initial semimajor axis from left to right: a = 0.55, 1.5, 2.0, 3.0, and 4.0 au. We show the MEGNO ⟨ Y ⟩ (top) and the Δe (bottom) stability indicators. The grey curve gives the centre of libration of the LidovKozai resonance (Eq. (4)), while the red vertical line gives the limit for tidal stability (Eq. (8)). 

Open with DEXTER  
In the text 
Fig. 8 Tidal evolution in the plane (e,J) obtained with a secular model for the same initial conditions from Fig. 7. The colour index gives the variation between the initial and final semimajor axes, log _{10}(Δa/a), and white stands for unstable orbits. The initial semimajor axes from left to right are a = 0.55, 1.5, and 2.0 au. The black vertical lines give the tidal stability limit (Eq. (8)). 

Open with DEXTER  
In the text 
Fig. 9 Stability maps of ⟨ Y ⟩ in the plane (a/a_{B},J) for some eccentricity values, with λ = Ω = 0 and ω = 90° for three different compact binary systems. From left to right: α Centauri, HD 196885, and HD 41004 are shown. From top to bottom, the initial eccentricity is e = 0, 0.3, 0.6, and 0.9. 

Open with DEXTER  
In the text 