© 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 S-type planets could have formed within the α Cen system (Quintana et al. 2002, 2007) provided that the collision velocities of late-stage 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 Earth-mass 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 ground-based radial velocity campaign has ruled out the presence of massive close-in 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, Andrade-Ines & 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.

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 mean-motion resonances; and (e,ϖ) or (J,Ω) to study secular resonances.

We solved the three-body equations of motion numerically with a Burlisch-Stoer 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⁵ yr, which corresponds to several periods of the secular variations.

We analysed the stability of a test planet in a S-type 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 N-body 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⁵ 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⁸ 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⁹ 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 long-term 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).

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 109 yr.

3. Action-angle maps

3.1. Mean-motion resonances

We begin our quest for stability regions by varying the pair (a,λ), which allow us to identify the presence of possible mean-motion 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 mean-motion 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 mean-motion resonances.

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 mean-motion resonance with the companion star. However, for small eccentricities we can observe some small resonant islands for the 15:1 and the 16:1 mean-motion 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 mean-motion 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 mean-motion 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.

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 three-body 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 Lidov-Kozai resonance (Lidov 1962; Kozai 1962). The particular geometry of this almost restricted three-body problem allows us to explore the phase-space 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 mean-motion 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 anti-aligned 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 non-linear 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 Lidov-Kozai 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 Lidov-Kozai 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 Lidov-Kozai resonance. Regular coplanar retrograde orbits are still possible in this case, but they already present eccentricities very close to instability.

4. New action-action maps

Maps involving the three actions (a,e,J) were extensively explored in previous studies (e.g. Andrade-Ines & Michtchenko 2014; Quarles & Lissauer 2016). However, they usually fix ω = 0°, which corresponds to a region of the phase-space that is always outside the libration zone of the Lidov-Kozai 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 action-action 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 Lidov-Kozai resonance (Eq. (4)).

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 mean-motion resonances.

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⁴ yr for nearby orbits. We scale our integrations time such that more distant orbits are integrated at least 5 × 10⁴ 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²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 Lidov-Kozai 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 mean-motion resonances to better understand the structures present in these regions.

We observe that stability islands are associated with mean-motion 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 Lidov-Kozai 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 Lidov-Kozai resonance (Eq. (4)), while the red vertical line gives the limit for tidal stability (Eq. (8)).

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 Lidov-Kozai 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 close-in 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)).

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₂ 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₀. 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₂/Q ~ 10^-2 − 10^-3, while for gaseous planets k₂/Q ~ 10^-4 − 10^-5 (Yoder 1995). However, the mass and radius of a gaseous planet is in general 10² 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 Jupiter-like planet with k₂/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 Lidov-Kozai 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 close-in 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 close-in 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 Lidov-Kozai 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¹ 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/aB,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.

We have seen that the Lidov-Kozai 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 Lidov-Kozai 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 S-type 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 mean-motion 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 Lidov-Kozai resonance.

As the eccentricity increases, the mutual inclination at the centre of the Lidov-Kozai resonance also increases. The Lidov-Kozai 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 Lidov-Kozai equilibrium, but only for close-in planets (a< 1.5 au).

In this paper we focused on the stability of S-type 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 Lidov-Kozai configurations in compact binary systems (a_B ~ 20 au) when tides are considered, through a close fly-by of a background star (Martí & Beaugé 2012, 2015).

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

All Tables

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).
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 109 yr.
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 mean-motion resonances.
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.
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 Lidov-Kozai resonance (Eq. (4)).
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 mean-motion resonances.
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 Lidov-Kozai resonance (Eq. (4)), while the red vertical line gives the limit for tidal stability (Eq. (8)).
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)).
In the text

	Fig. 9 Stability maps of ⟨ Y ⟩ in the plane (a/aB,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.
In the text