© ESO 2018

1 Introduction

A large portion of stellar systems are composed of multiple stars. Single stars account for approximately half, binaries a third, and triples a twelfth (Tokovinin 2014a,b). Over more than 200 yr, there have been numerous theoretical studies of the three-body problem and its restricted form. A few recent theoretical approaches to the stellar dynamics of triple systems, some incorporating numerical checks of their results, include Ford et al. (2000); Mardling & Aarseth (2001); Hamers et al. (2015) and Correia et al. (2016). Theoretical approaches can, however, have limitations in accurately representing some actual systems.

Numerical studies are a powerful tool with which to map the often discontinuous stability boundaries of systems without any restriction on their physical characteristics and configurations. Also, a number of studies have shown that, for systems that contain more than one planetary body, the orbits proposed initially were simply not dynamically feasible, which has promoted dynamical analyses as an integral part of the exoplanet discovery process (Horner et al. 2012).

A hierarchical triple stellar system consists of three bodies organized in two nested orbits – one central binary and a third star orbiting the center of mass of this binary at a larger distance. We shall call the two bodies constituting the binary stars 1 and 2, and the outer body star 3. We note that we have made no assumptions about the relative masses of the three stars. In many systems, star 3 is less massive than the central binary, but there are examples of systems with inverted mass ratios, such as HD 181068, for which star 3 is heavier than the binary (Derekas et al. 2011; Borkovits et al. 2013).

There are four possible types of planetary orbits in such a system (see Fig. 1):

• S1, S2, and S3 orbits comprise planets that orbit one of the individual stars (stars 1, 2 and 3 respectively) while being perturbed by the other stars. S1 and S2 orbits are not distinguishable and will be treated as a single S1 type;

• P1 orbits comprise circumbinary planets that orbit the central binary while being perturbed by the outer third star; and

• P2 orbits comprise circumtriple planets that orbit the whole system beyond star 3.

Planets have been found in 32 triple systems to date. Of these, S3 orbits have been found in 29 systems, S1 orbits in two and a P1 orbit in one. No planet in a circumtriple P2 orbit has yet been discovered.

Numerical studies of planets in multiple systems have focused on binaries. For example, in their key paper, Holman & Wiegert (1999) examined S-type and P-type planetary orbits in eccentric, coplanar binary systems and derived regression equations describing the stability bounds in terms of mass ratios and eccentricities. Eccentric and inclined (up to 50°) P-type orbits, for equal-mass binaries, are studied by Pilat-Lohinger et al. (2003). Work on S-type and P-type orbits is done by Musielak et al. (2005) for circular, coplanar orbits, defining the regions of stability as a function of the semi-major axis ratio between the star and planet and the mass ratio. Mudryk & Wu (2006) also analyze the numerical results of Holman & Wiegert (1999) for S-type orbits in coplanar systems and investigated the instability boundary at a higher resolution. A binary study that addressed a full range of planetary inclinations, for eccentric P-type orbits, was by Doolin & Blundell (2011), which described stability as a function of mass ratios and eccentricities.

One study that addressed retrograde planets in binary, coplanar, circular, S-type orbits was by Morais & Giuppone (2012), while Giuppone & Correia (2017) looked at S-type orbits, including retrograde ones, in a compact binary system, for various inclinations, eccentricities and orientation angles. Although none have yet been confirmed, planets in retrograde orbits within and around a triple system should be possible, similar to HAT-P-7b in its binary system.

One study of the four-body problem of planets within a triple system was by Verrier & Evans (2007). This study was confined to prograde, coplanar P-type orbits and derived regression relationships analogous to those of Holman & Wiegert (1999).

In a more general analysis, in addition to retrograde planetary orbits, it is also necessary to include highly inclined orbits of the outer star, since the resulting stellar Kozai resonance sculpts the geometry of the stable planetary region (Kozai 1962; Ford et al. 2000). It is also useful to consider retrograde motion of the outer star. For the 22 triple systems where mutual inclinations have been established, the mean is a high 63° and four of these systems are retrograde (Tokovinin 2017). One analytical study that addressed retrograde stellar orbits is by Farago & Laskar (2010). They address two limiting cases, the inner restricted problem, where one inner body (i.e., our star 2) has no mass, and the outer restricted problem, where the outer body (star 3) has negligible mass, describing the possible motions of particles in each regime and providing an expression for the boundary between the regimes.

We extended the number of parameters used in previous work to all the orbital elements of a triple’s stars and widened the ranges of both these and the relevant mass ratios, to encompass recent observations as well as potential future discoveries. For example, systems with highly inverted mass ratios, such as HD 181068, have been discovered but not yet modeled.

The objective of this work was to provide a broader mapping of the regions of planetary stability in generalized triples. This resulted in 24 semi-empirical models describing the stability bounds for planets in each possible combination of orbital motions in triples and for a range of stellar orbital elements and mass ratios.

2 Methods

2.1 Selection of parameter space

Our aim was to determine the orbital stability of both the host stars and their planets, where collisions between the stars or the ejection of one of the three stars (typically the least massive body) does not occur over secular timescales that are very long compared to the orbital periods. Our treatment of secular perturbations was based on classical Newtonian dynamics and assumed that all three bodies are point masses that do not interact other than through gravitation and do not evolve in any way. Relativistic and tidal effects, as well as spin, were not addressed. The triple system’s nomenclature is shown in the schematic illustration in Fig. 1.

There are two S-type orbits, as the S1 and S2 orbits are equivalent, and two P-type orbits, one circumbinary and one circumtriple. The stable region for circumbinary (P1) orbits will have an inner and outer bound, while that for circumtriple (P2) orbits will have an inner bound only. Circumstellar (S1 and S3) regions of stability will have their outer edges bounded. The objective is to determine the outer bound a_io of the stable regions for S1, S3, and P1 orbits and the inner bound a_oi of the P2 orbits. The following subsections list the selected parameter ranges for the stellar configurations.

Fig. 1 Triple system with S-type and P-type orbits.

2.1.1 Semi-major axis ratio a₂ /a₁

The lower limit was selected to be just inside the stability criterion for triples as given by Mardling & Aarseth (2001). The smallest ratio used was approximately three, lower than the smallest ratio of 4.1 found in the Eggleton & Tokovinin (2008) catalog of 285 triples and well above the a₂∕a₁ ≫ 2 requirement of the HJS symplectic integrator used. An upper limit of 100 was used; most of the semi-major axis ratios for triple systems fall within this range.

2.1.2 Inner mass ratio μ₁ = m₂ /(m₁ + m₂)

The range for this ratio is between zero and one. However, only the range 0 < μ₁ < 0.5 needs to be studied, as mass ratios of μ₁ and (1 − μ₁) are equivalent, other than for a 180° change in the longitude of the ascending node. An upper limit of one and lower limits of 0.001 and 0.1 were used for P and S orbits, respectively.

2.1.3 Outer mass ratio μ₂ = m₃/(m₁ + m₂)

Outer mass ratio definitions vary widely. In the Tokovinin (2014a) survey, the largest mass ratio (defined there as m₃ ∕m₁) found was 8.9 (for HIP 29860). However, 94% of mass ratios are below 2 and 99% are below 5. Taking an upper limit of five and assuming that the masses of the binary pair are broadly comparable, the equivalent value for our ratio definition is 2.5. A mass ratio greater than 1 implies an inverted system, where the outermost star is more massive than the aggregate inner binary. An example is the triple system HD 181068, which has a mass ratio of ~1.7. The lowest mass ratio found in the survey was 0.07. We therefore used ratios ranging from 0.2 or 0.001, for P and S orbits, respectively, up to 2.5.

2.1.4 Inner and outer eccentricities e₁ and e₂

The only data on inner and outer stellar eccentricities in triples is from Sterzik & Tokovinin (2002). The mean eccentricity of all orbits is high at 0.39 and in most (70%) of the systems the inner orbit is more eccentric than the outer orbit. The difference in eccentricities within these systems ranges from effectively 0 to 0.57. A study of 222 Kepler triples finds outer eccentricities spanning the full range, with a broad peak in the middle of the range and an unexplained narrow peak near e₂ ≃ 0.28 (Borkovits et al. 2015). The integrations, therefore, covered eccentricities in both inner and outer orbits ranging from 0 to 0.7–0.9 to cover fully this parameter range, with the higher limit required when Kozai resonance occurs.

2.1.5 Inner and outer inclinations i₁ and i₂

The inner inclination i₁ was always set to zero, so the outer inclination i₂ was the mutual inclination i₂ − i₁ of the outer orbit relative to the inner orbit. For the low-inclination simulations, small ranges of mutual inclination were selected, of 0°–60° for prograde orbits and 120°–180° for retrograde orbits. In the high-inclination integrations, the full range of 0°–180° was used. Borkovits et al. (2015) found that the distribution of mutual inclinations for 62 Kepler triples had a large peak at 0°–5°, indicating close-to-coplanar configurations, with a significant 38% portion of the systems in a secondary peak centred at 40°, suggesting Kozai effects.

2.1.6 Other orbital elements

The outer star’s longitude of ascending node Ω₂ and argument of periapsis ω₂ were both varied from 0° to 360° while those for the inner star were set to zero.

2.1.7 Test particles

Planets were represented by massless test particles. Their initial orbital elements ranged between selected lower and upper limits – for eccentricity these were 0–0.9, for inclination 0°–90° or 90°–180° for prograde or retrograde orbits, respectively, and for ω, Ω and mean anomaly M they were 0°–360°. Most initial orbital elements were randomly generated from uniform distributions. However, for inclination i₂ the direction of angular momentum should be uniform over the celestial sphere, which requires a uniform distribution in cos i₂ between [−1, 1].

A log scale was used for semi-major axes. The particles’ semi-major axes ranged between 0.02 AU and 0.9a₁ for S1 orbits and 0.02 AU and fa₂ for S3 orbits, where f varied between 0.07 and 0.60 depending on the stellar configuration. For P1 and P2 orbits, the lower limit was 0.01 AU while the outer limit was set conservatively at twice the distance at which the 5:1 mean motion resonance occurred. The planet orbiting closest to its central body is PSR 1719-14 b, at 0.0044 AU, while Kepler-42c has the closest orbit to a “normal” star, at 0.006 AU. The distance at which a test particle was stopped as being too close to the central body was set at 0.002 AU. The minimum semi-major axis used must be larger than this, and so 0.01–0.02 AU was selected for all orbit types. The furthest planet from its host star yet discovered is HIP 77900 b, at 3200 AU. The distance at which a test particle was assumed to have escaped from the central body was set at 10 000 AU.

2.2 Numerical methods

2.2.1 Integrator

We used thesymplectic HJS (Hierarchical Jacobi Symplectic) integrator developed by Beust (2003) from the Swift code (Levison & Duncan 1994, 2013). It is designed specifically to allow the integration of multiple stellar systems and to handle hierarchical systems of any number and structure, provided the hierarchy is preserved along the integration. The HJS algorithm has proved itself a valuable tool in many studies of hierarchical systems. Most triples are hierarchical, simply because if they were not, the system is likely to be unstable and fragment into a binary and an ejected third star.

2.2.2 Computational parameters

An integration time step of 1∕20 of the orbital period of the inner binary was normally used in the integrations, and varied to verify that results were not affected by numerical errors. Symplectic integration schemes usually ensure energy conservation with 10⁻⁶ –10⁻⁸ relative accuracy using this time step size (Levison & Duncan 1994; Beust 2003; Verrier & Evans 2007). We found an overall fractional change in the system energy ΔE∕E₀ of ~10⁻⁷ over a 10⁵ yr integration, although this varied widely depending on the triple configuration being integrated. Orbital instabilities tended to occur quickly – if an orbit did not become unstable in as little as 10–100 yr, it was usually stable up to 10⁵ yr. The fact that the stellar systems we investigated were largely compact was helpful, as 10⁵ yr was often equivalent to many hundreds of thousands of orbits of the outer star.

The secular evolution of orbital parameters was sampled for each of the configurations used, and an integration time of 10⁵ yr was found to be sufficiently long in most cases. Each batch of integrations included at least one run on the least compact configuration with an integration time of 10⁶ or 10⁷ yr, to confirm this.

The number of planetary test particles used varied from 1000 to 10 000, with higher numbers being used when the rate of instability and hence removal of test particles was high. For coplanar configurations, the initial test particle cloud took the form of a disk aligned with the plane of the inner binary, while for configurations where the outer star had a high inclination the test particle disk was aligned with the invariable plane. A procedure for automatically identifying the edges of the test particle cloud at the end of an integration and measuring its semi-major axis was incorporated into the HJS algorithm.

3 Results and discussion

3.1 Orbit types P1 and P2

3.1.1 Results

For inner (P1) and outer (P2) orbits, the respective outer and inner edges of the cleared area of the test particle cloud were standardized by taking them as ratios of the semi-major axis of the outer star, a₂. These standardized dependent variables a_io /a₂ and a_oi /a₂ delineate the bounds of the stable regions and are termed the inner and outer critical semi-major axis ratios respectively.

The mean critical semi-major axis ratio for these inner and outer orbits, for the eight possible combinations of orbital motion, are shown in Table 1. These are the average ratios found over all the combinations used in the parameter space. For P1 planetary orbits, the mean critical ratio is materially different (36%) for prograde and retrograde planetary orbits. Both these ratios are slightly (~2%) larger for retrograde stellar orbits. For P2 orbits, the mean critical ratio is similarly (−33%) different for prograde and retrograde planetary motions. For retrograde stellar orbits, these ratios are slightly (~3%) smaller. The difference made by the orbital direction of the outer star is small for both prograde and retrograde planetary orbits, but is statistically meaningful for inner retrograde planetary orbits and outer prograde planetary orbits.

The critical semi-major axis ratios were then regressed against the parameters discussed in the previous section using linear regression to extract semi-empirical relationships of the form $$\frac{a_{\text{io}}}{a_2},\frac{a_{\text{oi}}}{a_2}=C+b_1{\text{μ}}_1+b_2{\mu }_2+b_3{e}_1+b_4{e}_2+b_5{i}_2+b_6{\Omega }_2+b_7{\text{ω}}_2$$(1)

for the outer bounds of S1, S3 and P1 orbits, and the inner bounds of P2 orbits, respectively, where C and b_i are the regression constant and coefficients. In all the regressions of P-type and S-type orbits, the univariate relationships between thecritical ratios and the independent variables were linear, with one exception – the relationship between a_io ∕a₂ and μ₂ in S3 orbits.

The results of the regressions are tabulated in Table 2. In this and subsequent tables of regression results, the various parameters are defined as follows: σ–standard deviation of critical ratios; R²–coefficient of determination, F–F–statistic for overall significance, SE–standard error and MAPE–mean average percentage error, all being for the overall regression fit; t–t–statistic for individual coefficients and N-number of data points.

The various regressions resulted in 24 regression constants and 192 coefficients. The semi-major axis ratio a is included as an error check as the critical ratio scales with it and its coefficient should be zero.

Examining the constants, for circular, coplanar outer stellar orbits, the inner and outer stability boundaries for prograde planetary orbits are found at around 0.4 and 2.4 times the distance of the outer star, respectively, and for retrograde planetary orbits, at around 0.6 and 1.5 times, respectively. The parameters i₂, Ω₂ and ω₂ have a negligible influence on the critical ratio in these low-inclination cases.

For P1 orbits, the dominant influence on the critical semi-major axis ratio is only the outer star’s eccentricity, particularly for retrograde planetary orbits. For P2 orbits this effect is even stronger, for both planetary motions. In the case of a retrograde outer star these influences are largely unchanged for P1 orbits, but for P2 orbits the inner mass ratio becomes important for a prograde planet.

The outer star dominates the regions of stability in a triple system, with its eccentricity having by far the largest influence. The configuration of the inner binary has little effect on either inner or outer orbits. This results from the fact that, as a consequenceof the Mardling stability limit, the outer star is sufficiently far away that the inner binary effectively resembles a single point mass. These conclusions apply to both prograde and retrograde planetary orbits. However, the greater stability of retrograde planetary orbits results in inner and outer bounds that are closer to the outer star compared with the prograde case.

The regression coefficient of the outer eccentricity e₂ is large – for highly eccentric orbits of the outer star, the critical semi-major axis ratio can expand by over 80% for prograde planetary orbits and more than double for retrograde orbits. The inner bound can shrink by a quarter for prograde planetary orbits and by over 80% for retrograde orbits.

The effect of the direction of motion of the outer star on the planetary stability bounds is shown in Table 3. The significance level is shown for the mean critical ratio, and differences in the regression constant are shown for comparison. Although the absolute differences in mean critical ratio are small, two are statistically significant.

The effect of the direction of motion of planets on their stability bounds is shown in Table 4. The absolute differences are all large, averaging 37% for P1 orbits and 31% for P2 orbits, with the signs being consistent with the greater stability of retrograde orbits.

3.1.2 Stability bounds and the outer star’s eccentricity

Since the eccentricity of the outer star is by far the most influential variable on both the inner and outer planetary stability bounds, a series of integrations was run to examine the relationship between these two variables. Typical relationships for the outer bounds are illustrated in Fig. 2, for one stellar configuration.

For the outer bound the critical ratio should be an increasing function of outer eccentricity e₂. This is true for all four cases shown, although for the two retrograde planet cases the critical ratio flattens out for e₂ ≳ 0.6. There are discontinuities in the critical ratio, resulting from resonances, in each case. This is most visible in the prograde star and prograde planet case, with gaps occurring at eccentricities of 0.10, where the critical ratio jumps from 2.37 to 2.70; and at 0.39, with the ratio undergoing a step change from 2.82 to 3.05. These instabilities are far less pronounced in the two retrograde planet cases.

A retrograde outer body allows the outer bound for a prograde planet to move substantially closer to it, with a critical ratio at e₂ = 0 of ~1.3, compared with ~2.1 for a prograde outer star. However, this difference diminishes with increasing outer eccentricity, with both critical ratios converging toward 3.2 as this eccentricity approaches unity.

Typical relationships for the inner bounds are illustrated in Fig. 3, for the same stellar configuration. For the inner bound, the critical ratio is expected to be a decreasing function of outer eccentricity. This is the general trend in each case, albeit with high scatter, which reflects that the inner orbit bounds are always more diffuse than the outer orbit bounds. This is a result of the greater sparsity of surviving test particles compared with outer orbits, sinceinner orbits tend to be more chaotic. No stable inner orbits are found for outer eccentricities above 0.7.

Fig. 2 Outer stability bound for P2 orbits as a function of outer star eccentricity, for a = 100 AU, μ1 = μ2 = 0.5, e1 = i1 = Ω2 = ω2 = 0 and i2 = 0° and 180°.

Fig. 3 Inner stability bound for P1 orbits as a function of outer star eccentricity, for a = 100 AU, μ1 = μ2 = 0.5, e1 = i1 = Ω2 = ω2 = 0 and i2 = 0° and 180°.

3.1.3 Comparison with observed P1 and P2 orbits

The only P1orbit found in a triple system is in HW Virginis, with a semi-major axis ratio of 0.37 (Beuermann et al. 2012). No circumtriple P2 orbits have been discovered to date. The planet with the smallest P1 orbit in a binary, in terms of absolute semi-major axis, is Kepler-47b (Orosz et al. 2012), with a semi-major axis ratio (here relative to the binary) of a_io ∕a₁ = 3.54. The smallest semi-major axis ratio of a planet’s orbit is for Kepler-16b, with a_io ∕a₁ = 3.14 (Doyle et al. 2011). Both lie well outside the smallest mean critical semi-major axis ratios of ~0.1 found in thesimulations, so the possibility exists of finding planets in even smaller orbits.

3.1.4 P1 and P2 orbits in triples compared with binaries

To highlight how the P1 and P2 planetary stability bounds in a triple differ from those in a binary, the integrations were repeated after reducing the triple to a binary by condensing the inner binary into a single body. The results of this approximation for the mean critical semi-major axis ratios for the prograde stellar case are compared in Table 5.

Both the inner and outer mean critical semi-major axis ratios move toward the central star, for both prograde and retrograde planetary orbits. While these movements are statistically significant, they are nevertheless small, confirming that the influence of the outer star is dominant. The regression coefficients are compared in Table 6.

In all cases, aside from the constant the only variable of major influence is the eccentricity of the outer star, e₂. Generally, the effect of the outer star’s eccentricity is far larger on the outer stability bounds than on the inner stability bounds. For inner orbits, its influence is larger in binaries than in triples, for both prograde and retrograde planetary orbits. For outer orbits, its effect in triples and binaries is essentially the same for retrograde orbits, but for prograde orbits it has a larger influence in triples than in binaries. Compared with binaries, the influence of the outer star’s eccentricity is significantly smaller for inner prograde orbits and materially larger for outer prograde orbits; there are no similar differences for retrograde orbits. Overall, the relatively small differences between triples and binaries result from the Mardling stability limit for triples, which precludes them from becoming too compact.

3.2 Orbit types P1 and P2 in highly inclined triple systems

3.2.1 Vertical characteristics of the stability region

Here we examine the stable planetary region for large inclinations of the outer star relative to the invariable plane of the triple system. Once the orbit of the outer star is no longer coplanar with the orbit of the inner binary, in other words, the mutual inclination of the two orbits is no longer zero, the stable planetary region will also no longer be coplanar but will extend vertically, with its shape being sculpted primarily by the outer star. An example is illustrated in Fig. 4, where the x–y plane is aligned with the invariable plane. The outer limit of the test particle cloud is of no relevance here, since this is simply where it has been truncated.

Examining the paths of individual test particles over time yields the following observations. Firstly, their orbits librate, with the argument of periapsis ω oscillating. Secondly, the semi-major axis of their orbits, a, remains effectively constant. Thirdly, their orbits remain circular, with their eccentricity staying very close to zero. Fourthly, their orbits vary in inclination. In an integration they are initially coplanar with the invariable plane and then become more inclined over time, eventually oscillating around a final non-zero inclination with a period much longer than those of the other orbital elements. The typical annular vertical “chimney” shown in Fig. 4 therefore consists of test particles in circular orbits of constant semi-major axis but varying inclination. Fifthly, some test particles develop Kozai resonance and some of these eventually undergo orbit flipping into retrograde motion.

There isa limit to the maximum inclination of planetary orbits. The number of stable orbital bounds falls exponentially with increasing inclination, as shown in Fig. 5, where the lines end when there are no more stable orbital bounds.

Fig. 4 Position in space of surviving test particles after 100 kyr. P1 and P2 prograde orbits, with a = 40 AU, e1 = e2 = 0, i1 = 0°, i2 = 65°, μ1 = 0.5, μ2 = 2.0.

Fig. 5 Number of stable P1 and P2 prograde orbits remaining after 1 Myr as a function of initial planetary inclination, for various triple configurations.

3.2.2 Planetary stability bounds and the outer star’s inclination

Both prograde and retrograde stellar orbits were used, and planetary motions were prograde. In practical terms, we know that for real bodies, as opposed to test particles, the Kozai critical inclination is larger than the theoretical 39°, for example Grishin et al. (2017). Also, if the ratio of the initial angular momentum of the outer orbit to that of the inner orbit is ≳4, then significant Kozai resonance usually occurs, for example Beust et al. (2012).

Since we used a relatively massive outer star, the integrations were split, not by the theoretical Kozai critical inclination, but into what were experimentally determined to actually be non-Kozai and Kozai regimes, that is, those of relatively low inclination whose nodes circulate and those of high inclination where the nodes librate about 90° or 270°. For the integrations, the sample was first split into prograde and retrograde stellar cases. Each case was then split into subsets where Kozai resonance did or did not occur. Figure 6 shows the critical semi-major axis ratios for the inner and outer planetary orbits against the “absolute” inclination i₂ (i.e., i₂ for prograde orbits and 180° − i₂ for retrograde orbits), for both Kozai and non-Kozai regimes. As usual, there is more scatter for inner orbits than outer orbits.

For inner orbits, the critical semi-major axis ratio has almost no dependence on the outer star’s inclination for progradestellar orbits, and only a weak one for retrograde orbits. Retrograde stellar orbits allow the stable planetary region to extend further out (with a regression constant of 0.69 versus 0.46), but as inclination approaches 90°, this effect decreases. (The regression lines for inner and outer orbits should ideally coincide at an inclination of 90°.)

For outer orbits, the critical ratio also has minimal dependence on inclination for prograde stellar orbits. However, it can be seen in the bottom panel of Fig. 6 that this low average dependence is made up of two parts. For low inclinations the critical ratio is relatively constant at 1.9–2.0, but as inclinations increase the critical ratio begins to decline, in other words, the orbital stability bound begins to move inwards. This makes intuitive sense because once Kozai resonance begins, as the outer star’s inclination increases, the eccentricity of the inner binary decreases, which allows this to occur. This divergence begins at around 45° and is probably a manifestation of the start of Kozai resonance. This suggests that for real bodies which may also be relatively large, Kozai resonance is more likely to begin when i₂ ≳ 45° than for i₂ > 39°.

For retrograde stellar orbits, dependence on i₂ is much stronger. In this case planetary orbits can again approach much closer to the outer star (with a regression constant of 1.2 versus 1.9), but this increased stability declines as inclination rises to 90°. Also of interest are the regions of stability interspersed with gaps of instability.

The mean semi-major axis ratios for the non-Kozai and Kozai cases are shown in Table 7. Generally, with a retrograde outer star the inner orbits move outwards toward this star and the outer orbits move inwards, reflecting the greater stability of retrograde orbits. For a prograde outer star, the existence of Kozai resonance makes no difference to the critical ratios of either theinner or outer bounds, which remain virtually identical. For a retrograde outer star, however, the bounds are very different. When Kozai resonance occurs, the inner bound contracts inwards, while the outer bound moves outwards, and by a larger proportional amount. For both inner and outer orbits, these movements are in a direction opposite to that usually induced by a retrograde stellar orbit. This shows, unsurprisingly, that Kozai resonance increases planetary instability.

At the 5% significance level, the critical semi-major axis ratios for planetary orbits under non-Kozai and Kozai regimes do not differ for a prograde outer body, but for the retrograde case they are substantially different, by an absolute 20–30%. The number of stable inner orbits found under Kozai resonance with a prograde outer star is small.

The regression data for these four cases is shown in Table 2. None of the orbital parameters have any significant influence on the critical semi-major axis ratio, with the possible exception of μ₁ in the prograde star/Kozai case. This stability bound is therefore effectively represented by the constant. The data for the constant in the four cases is shown in Table 8.

For a prograde outer star, Kozai resonance results in the constant being 44% smaller for inner orbits and 22% larger for outer orbits. A retrograde outer body does not result in a significant change in either constant.

Fig. 6 Critical semi-major axis ratios versus outer star inclination – panel a: inner stability bound; panel b: outer stability bound.

3.3 Orbit type S1

3.3.1 Results

The S1 and S2 orbits around stars 1 and 2, respectively, are interchangeable, so only one case is examined, denoted S1. The mean critical semi-major axis ratios for S1 orbits, for the four possible combinations of orbital motion, are shown in Table 9.

For S1 planetary orbits in triples, the mean critical semi-major axis ratio is in the range 0.180–0.197. S1 planetary orbits on average extend ~8% further out when the outer star is in a retrograde orbit. For prograde outer stellar orbits, retrograde planetary orbits extend slightly (3%) further out while for retrograde stellar orbits they are virtually the same. The results of the corresponding regressions are shown in Table 2. The key determinants of the critical semi-major axis ratio for planetary S1 orbits are only the inner binary’s mass ratio and eccentricity, and for a retrograde outer star the influence of eccentricity is much weaker.

3.3.2 Comparison with observed S1 orbits

The two S1 orbits in triple systems found to date, HD 126614 A b and HD 2638 b, have semi-major axis ratios of 0.0649 and 0.00172, respectively. There is some uncertainty whether a third, Fomalhaut b, is a planet or a dust cloud or disk.

Among the closest S1 orbits found in a binary are 0.7 AU for OGLE-2013-BLG-0341L b, with a semi-major axis ratio a_io ∕a₁ of 0.058–0.041 (Gould et al. 2014); 1.09 AU for HD 59686 b, with a_io ∕a₁ ~ 0.080 (Trifonov et al. 2018) and 0.382 AU for KOI-1257 b, giving a_io∕a₁ ~ 0.072 (Santerne et al. 2014). Although some of our integrations showed stable bounds with semi-major axis ratios as low as 0.015, they also extended to almost 0.8, with the greatest concentration in the region 0.03–0.25, so there appears to be significant scope for planets to be found further out from their host stars.

3.3.3 S1 orbits in triples compared with binaries

The triple was reduced to a binary by setting the outer star’s mass to a negligible value and the integrations were repeated. The results of this approximation for the mean critical semi-major axis ratios for the prograde stellar case are compared in Table 10.

The outer star of a triple has a significant constraining influence on S1-type orbits compared with the binary case. In triples the added influence of the outer star makes S1 orbits move inwards, with the mean critical semi-major axis ratio reducing by over 30%, from ~0.270 to ~0.183, while the difference in critical ratio between prograde and retrograde planetary orbits shrinks from 12% to 2%. The regressions for the binary case are compared with those for the triple case in Table 11.

The results for prograde and retrograde planetary orbits in triples and binaries are qualitatively similar, in that the only variables of influence are the inner mass ratio μ₁ and eccentricity e₁. This effect is not attributable to Kozai, as the mutual inclination has no influence. However, the effect of the mass ratio is significantly larger in triples than in binaries, while that of eccentricity (and the constant) are effectively the same.

3.4 Orbit type S3

3.4.1 Results

The mean critical semi-major axis ratios for S3 orbits, for the four possible combinations of orbital motion, are shown in Table 12.

For S3 planetary orbits in triples, the critical semi-major axis ratio is in the range 0.287–0.361 (compared with 0.180–0.196 for S1 orbits) with wide ranges around these values depending on the parameters of the triple, and is essentially independentof the direction of motion of the outer star. The critical semi-major axis ratio is around 25% larger for retrograde planetary orbits than for prograde orbits, significantly larger than the 8% difference found for S1 orbits.

3.4.2 Comparison with observed S3 orbits

The smallestand largest critical semi-major axis ratios a_io∕a₂ calculated for planets in S3 orbits, based on data in Wagner et al. (2016), are 0.00061 and 0.265, respectively. Recently discovered KELT-4 Ab appears to have a very small a_io ∕a₂ of 0.000132 (Eastman et al. 2016), as does Proxima Centauri b’s 3.3 × 10⁻⁶ (Anglada-Escudé et al. 2016). Although the integrations found some stable bounds with semi-major axis ratios as low as 0.009 for a few stellar configurations, there were very few critical ratios lying below 0.1 and most ranged from 0.1 to 0.6, with a few as high as 0.9. This suggests that there may be planets at much greater distances from their outer host stars than discovered to date.

3.4.3 S3 orbits in triples compared with binaries

For S3 orbits the triple was reduced to a binary by merging the inner binary into a single star, as for the P1/P2 case, after which the integrations were repeated. The results of this approximation for the prograde stellar case are compared in Table 13.

The S3 orbits in a triple system have a smaller mean critical semi-major axis ratio than in a binary, since the inner binary in a triple has a larger effect on the planet than a single body of equal mass at the same distance. The average difference in S3 orbits for binaries and triples is around 5%, which is much smaller than the average 50% difference found for S1 orbits and more comparable with the average (absolute) 9% difference for P1/P2 orbits. The regressions for the binary case are compared with those for the triple case in Table 14.

Aside from the constant, the only orbital element of importance is e₂, the eccentricity of the outer star containing the S3 orbit. The next largest influence, that of the outer mass ratio μ₂, is far smaller. The effect of the mass ratio μ₁ is smaller still, unlike in the S1 case, where it was of equal quantitative importance to e₂. These coefficients are virtually identical for prograde and retrograde planetary motions as well as the direction of motion of the outer star.

3.5 Comparison with previous work

The only empirical work on triples to compare with is that of Verrier & Evans (2007). Our prograde P1 mean critical ratio of 0.383 is 18%smaller than the 0.466 from that study, while our P2 mean critical ratio of 2.94 is close to their 2.92. The reason for the P1 discrepancy is not clear. Although a larger parameter space and larger clouds of test particles were used in our study, the boundary ofP1 orbits is always much more tenuous than that for P2 orbits, and the selection of a density cutoff to define its edge is necessarily arbitrary. So while results should be consistent within a study, a difference of this magnitude across studies is quite possible for P1 orbits.

A largerbut still sparse area of overlap exists between our “triple-reduced-to-binary” cases and previous work on binaries. Regression constants from this study are compared with those of previous empirical studies, where provided, in Fig. 7.

For prograde S1 and P1 orbits, our results are very close to those in previous work. There have been no empirical studies on S1 retrograde orbits. For P1 retrograde orbits, there are two other results; ours coincides with the lower one, from Doolin & Blundell (2011). The critical semi-major axis ratios found in previous research are shown in Fig. 8, together with our ratios and regression constants for comparison.

For binary S1 prograde planetary orbits, previous results for the critical semi-major axis ratio varied widely, ranging from 0.22 to 0.46 and averaging 0.38. Our corresponding mean critical ratio of 0.26 and regression constant of 0.36 are within this range. For P1 prograde orbits, the range is 1.60–3.85, the study by Holman & Wiegert (1999) being the lowest, and with an average of 2.44. Our mean critical ratio and constant of 2.70 and 2.71, respectively, are 11% higher than this.

Two empirical studies have addressed retrograde planetary orbits in binaries. The one for S1 orbits finds a critical ratio of 0.60, with our results for the mean critical ratio and constant, of 0.29 and 0.38, respectively, being around 45% lower. The study for P1 orbits, by Doolin & Blundell (2011), finds critical ratios in the range 1.3–2.7, comparable with our results that have a mean critical ratio of 1.69 and a regression constant of 1.25. Despite the small sample sizes, the above results generally compare well. The one relatively large difference, for retrograde S1 orbits, suggests further investigation.

Fig. 7 Comparison of results for binaries - regression constants.

Fig. 8 Comparison of results for binaries - critical semi-major axis ratios and regression constants.

4 Conclusions

We have investigated numerically the long-term stability of planets in a large number of possible orbits around the stars in a triple system and provided a generalized mapping of the regions of stability. This was done for prograde and retrograde motion of the planets and for the outer body of the triple, as well as for highly inclined orbits of the outer body. The constructionof multiple regression equations resulted in 24 semi-empirical models, one for each type of orbital configuration.

The greatest influences on a planet’s critical semi-major axis ratio are as follows: for P1 and P2 orbits, the eccentricity of the outer body and, to a far lesser extent, the inner mass ratio; for S1 orbits, the inner binary’s eccentricity and inner mass ratio; for S3 orbits, the outer body’s eccentricity and outer mass ratio. Further details can be found in Busetti (2018).

We extended the number of parameters used to all relevant orbital elements of the triple’s stars and their mass ratios. We also expanded these parameters to wider ranges that will accommodate recent and possibly future observational discoveries. The investigation of how the regions of planetary stability in triples differed from those in binaries enabled comparison with, and confirmation of, some previous studies. It also contributed some new data points.

These generalized results can be useful in the investigation of observed systems, providing a fast method of determining their stability bounds within the large parameter space that results from observational uncertainties. The relationships expressed in the regression models can also be used to guide searches for planets in triple systems and to select candidates for surveys of triple systems. The geometry of the stable zone indicates not only where to look for planets but also the most appropriate technique to search for them.

References

All Tables

All Figures

	Fig. 1 Triple system with S-type and P-type orbits.
In the text

	Fig. 2 Outer stability bound for P2 orbits as a function of outer star eccentricity, for a = 100 AU, μ1 = μ2 = 0.5, e1 = i1 = Ω2 = ω2 = 0 and i2 = 0° and 180°.
In the text

	Fig. 3 Inner stability bound for P1 orbits as a function of outer star eccentricity, for a = 100 AU, μ1 = μ2 = 0.5, e1 = i1 = Ω2 = ω2 = 0 and i2 = 0° and 180°.
In the text

	Fig. 4 Position in space of surviving test particles after 100 kyr. P1 and P2 prograde orbits, with a = 40 AU, e1 = e2 = 0, i1 = 0°, i2 = 65°, μ1 = 0.5, μ2 = 2.0.
In the text

	Fig. 5 Number of stable P1 and P2 prograde orbits remaining after 1 Myr as a function of initial planetary inclination, for various triple configurations.
In the text

	Fig. 6 Critical semi-major axis ratios versus outer star inclination – panel a: inner stability bound; panel b: outer stability bound.
In the text

	Fig. 7 Comparison of results for binaries - regression constants.
In the text

	Fig. 8 Comparison of results for binaries - critical semi-major axis ratios and regression constants.
In the text