Orbital Stability of Circumstellar Planets in Binary Systems

One of the defining features when expanding from the two-body (Kepler) problem to the three-body problem (TBP) is the emergence of chaos (see details in Mudryk & Wu 2006; Musielak & Quarles 2014), or a sensitivity to initial conditions where the future state of a system is indeterminate after a sufficient time. The possible outcomes within the TBP are broadly defined as stable (i.e., the system is bounded in space and consists of three bodies for all time) or unstable (i.e., at least one body is no longer bounded in space or a collision occurs). Notice that these definitions do not explicitly include chaos—where bundles of initial conditions can be either stable or unstable, as well as chaotic. There is a transition between stable and unstable states that has been identified within astronomical systems (Lecar et al. 1992; Smith & Szebehely 1993; Goździewski et al. 2001; Pilat-Lohinger & Dvorak 2002; Cincotta et al. 2003; Cuntz et al. 2007; Eberle et al. 2008; Quarles et al. 2011; Giuppone et al. 2012; Satyal et al. 2013) using chaotic indicators similar to the Lyapunov exponent (Benettin et al. 1980; Gonczi & Froeschle 1981; Froeschlé et al. 1997; Cincotta & Simó 1999, 2000).

Due to this transition region, we must be more precise about our definitions for stability and allow for a more probabilistic consideration. Moreover, we cannot compute the evolution of a system for all time, and thus we define a given initial condition as potentially stable if the planet survives for a predefined timescale (∼10⁵ yr) around its host star. Previous works (David et al. 2003; Fatuzzo et al. 2006; Quarles & Lissauer 2016) have shown that most of the initial conditions that will become unstable over longer timescales exist near mean motion resonances (MMRs) at the edge of stability; our approach excludes the resonant region, due to its dependence on the initial mean longitude of the planet. Unstable initial conditions are those that do not survive for the required timescale, either due to a collision of the planet with either of the host stars or because the distance to the host star exceeds 200 au (i.e., twice the largest binary semimajor axis considered).

Observations of planets in binaries are usually so limited that a full initial condition is not available. Rabl & Dvorak (1988) and Holman & Wiegert (1999) accounted for this by prescribing their stability criterion in terms of a set of observables that are the most readily obtained (e.g., binary mass ratio, binary eccentricity, semimajor axis ratio). The binary mass ratio μ (= M_B/(M_A + M_B)) and the binary eccentricity e_bin can be deduced from photometric and/or radial velocity observations. The semimajor axis ratio a_p/a_bin can be determined if the respective orbital periods are well-determined; otherwise, projected values are typically used. Due to the long orbital periods involved and small number of observations, not many other observables are typically known. As a result, Holman & Wiegert (1999) assumed that the orbital planes of the binary and planet are aligned (i.e., ${i}_{p}={i}_{\mathrm{bin}}={\omega }_{p}={{\rm{\Omega }}}_{p}={\omega }_{\mathrm{bin}}={{\rm{\Omega }}}_{\mathrm{bin}}=0^\circ$) and the planetary orbit is initially circular (i.e., e_p = 0). Additionally, the restricted TBP was implemented, where a massless test particle was used for the planet, as a matter of numerical efficiency.

A combination of parameters (μ, e_bin, and a_p/a_bin) are numerically evolved for eight equally spaced planetary phase angles, and if all eight of the trials are deemed stable, a critical semimajor axis a_c is determined. Most of our simulations follow this general approach, except we use an Earth-mass particle that gravitationally interacts with the binary and evaluate 91 initial planetary phase angles so that we can investigate the probabilistic transition region between stable and unstable orbits. The stability limit a_c is defined in our analysis as the largest semimajor axis ratio a_p/a_bin, where all 91 initial phase angles (0° ≤ λ_p ≤ 180°) survive up to our predefined timescale.

The numerical simulations in our study use a integration package that has been modified such that the orbits of planets in binaries are efficiently evaluated (mercury6; Chambers et al. 2002). This well-established code was developed for planet formation simulations within α Centauri AB (Quintana et al. 2002), where a hybrid scheme allows for switching between a symplectic integration method (i.e., Wisdom–Holman splitting (Wisdom & Holman 1991, 1992)) and an adaptive method (e.g., Bulirsch–Stoer (Press et al. 1993)). We check a subset of our runs using a newer N-body code, REBOUND, with the IAS15 integrator (Rein & Liu 2012; Rein & Spiegel 2015) to ensure the accuracy of our results.

We set up our simulations starting from a set of unitless parameters (μ, e_bin, and a_p/a_bin), which allows for scalability of our results to many different dynamical systems. Table 1 gives the range of these parameters, along with the range of inclination and mean longitudes for defining an orbit. The total mass of the binary is equal to one solar mass (M_A + M_B = 1 M_⊙), and the individual stellar masses are determined via μ (= M_B/(M_A + M_B)). The initial planetary semimajor axis a_p is 1 au and the initial value of the semimajor axis ratio is used to determine the appropriate binary semimajor axis (∼1–100 au). The initial phase of the binary begins at periastron (λ_bin = 0°), which lies on the positive x-axis (ω_bin = 0°) within the orbital plane of the binary. All of our simulations use the binary orbital plane as a reference (i_bin = Ω_bin = 0°), where the initial planetary inclination i_p is the relative angle between the planetary and binary orbital planes. Moreover, these orbital planes begin nodally aligned (Ω_p = Ω_bin).

The accuracy of our simulations is controlled by choosing a time step that is 5% of the planetary orbital period T_p (∝(1 − μ)^−1/2). Stability studies of α Cen AB (Wiegert & Holman 1997; Quarles & Lissauer 2016) have shown that the region of stability does not change appreciably until i_p is greater than 40°, and the extent of stability increases for retrograde orbits. There is an intermediate region (40° ≳ i_p ≳ 140°) where the Lidov–Kozai (LK) mechanism (Kozai 1962; Lidov 1962) drives large eccentricity oscillations that can reduce stability zones unless the planet and binary are apsidally misaligned by 90° (Giuppone & Correia 2017). We probe a wide range of planetary inclinations (0°, 30°, 60°, 85°, 95°, 120°, 150°, and 180°) for circular binaries, and then we limit our investigation of eccentric binaries to four planetary inclinations: 0°, 30°, 45°, and 180°. The mass ratio μ is varied uniformly in steps of 0.01 from 0.01 to 0.99, where we also evaluate the special cases of 0.001 and 0.999. The binary eccentricity e_bin is sampled in 0.01 steps from 0 to 0.8, while the semimajor axis ratio a_p/a_bin ranges from 0.01 to 0.8 in steps of 0.001. While the initial phase of the binary remains fixed (λ_bin = 0°), the mean longitude of the planet is varied in 2° steps from 0° to 180°.

The range of parameters that we explore is quite broad—and necessarily so, in order to provide the most general results. Our previous study of circumbinary planets (Quarles et al. 2018b) exploited a symmetry in the initial mean longitude, and allowed for the computations to be performed efficiently. To evaluate whether similar symmetries exist, we find the maximum eccentricity a planet attains when beginning from a circular orbit. Figure 1 illustrates the maximum eccentricity that a coplanar planet attains for stable configurations as a function of the initial semimajor axis ratio a_p/a_bin and mean longitude λ_p, where the white regions represent initial conditions that are unstable on a 500,000 yr timescale. The panels are labeled to show the mass ratio μ and stability limit a_c determined. As the mass ratio increases (top to bottom), the stability limit a_c decreases and the maximum eccentricity increases by an order of magnitude. Additionally, the first column shows a symmetry in the mean longitude that appears more strongly with increasing μ. The growth of the maximum eccentricity is more pronounced when the binary eccentricity e_bin increases (left to right). Indeed, these trends are known from secular studies (Andrade-Ines et al. 2016; Andrade-Ines & Eggl 2017, and references therein), and further indicates the accuracy of our numerical simulations. These trends justify our numerical setup (see Section 2.2) and particular restrictions in mean longitude when surveying the full parameter space.

Zoom In Zoom Out Reset image size

Other dynamical effects, such as the Lidov–Kozai (LK) mechanism (Kozai 1962; Lidov 1962), can arise when considering inclined orbits and affect how much we can exploit particular symmetries. In Figure 2, we perform similar simulations as in Figure 1, but for inclined orbits with a constant mass ratio (μ = 0.1). The planetary inclination i_p and determined stability limit a_c are given as tuples in the lower left of each panel. The difference between a planet inclined by 30° relative to the coplanar case is relatively minor, where variations occur mainly at large semimajor axis ratio that depend on the initial mean longitude. When the planetary inclination is increased to 45°, the stability limit a_c decreases more substantially relative to the coplanar case and the typical maximum eccentricity is ∼0.5 across most semimajor axis ratios. Finally, our retrograde (i_p = 180°) runs demonstrate a ∼25% increase in the stability limit a_c relative to the respective coplanar runs. This is roughly what one would expect from previous studies of the restricted TBP, where stable particles fill a larger portion of the Hill radius (Henon 1970). Overall, Figures 1 and 2 demonstrate that the range of mean longitude required to determine the stability limit can be reduced to 0–180° independent of the initial planetary inclination i_p, and it is possible that a separate condition can be placed on the maximum eccentricity, where ${e}_{\max }\lesssim 0.5$.

Zoom In Zoom Out Reset image size

One of the oldest problems studied within orbital dynamics is the circular restricted three-body problem (CRTBP), where the two more massive primaries begin on circular orbits and the third body is not massive enough to substantially alter the orbit of the primaries; see Musielak & Quarles (2014) for a review. Due to the restriction on the binary orbit's eccentricity, an integral arises that reduces some of the complexity. The solution to the integral, or Jacobi constant C_J in the CRTBP, provides a natural means to orient our analysis. The so-called zero velocity contour (ZVC) for a given mass ratio μ and semimajor axis ratio a_p/a_bin describes a topological boundary that limits the trajectory of a planet and can be related to the Jacobi constant (Eberle et al. 2008).

We use this formalism to explore the changes in the trajectory for a single initial condition (μ = 0.3, a_p/a_bin = 0.369, and λ_p = 180°) with respect to a changing planetary inclination. Figure 3 shows each trajectory over one orbit of the binary using inertial or static coordinates (left column) and rotated coordinates (right column), where the rotation rate is equal to the binary mean motion. The nominal location of the N:1 MMRs, in terms of semimajor axis ratio, is ${a}_{p}/{a}_{\mathrm{bin}}={((1-\mu )/{N}^{2})}^{1/3}$, which places the initial condition between the 4:1 and 3:1 MMRs. For the coplanar case in Figure 3(a), we trace the trajectory of the planet to nearly complete three orbits within one binary orbit. This results in a quasi-periodicity for the trajectory in the rotated coordinates (Figure 3(b)) and is bounded by the ZVC. When the planet is inclined by 30° in Figure 3(c), it follows a similar trajectory within a tilted plane and small differences appear due to the projection onto the X–Y plane. The planetary orbit is more eccentric (not just in appearance due to projection) once the inclination is increased to 45° in Figure 3(e). A retrograde (i_p = 180°) planet executes the most regular trajectory, where almost five orbits are now completed within a single binary orbit (Figure 3(g)) and the trajectory is periodic in rotated coordinates (Figure 3(h)). From these example trajectories, we expect that the stability results for planetary inclinations up to 30° to not differ too much from those determine for a coplanar orbit. The increased eccentricity from the Lidov–Kozai mechanism will introduce a dependence on the initial mean longitude, due to apsidal precession, and retrograde orbits will be stable to significantly larger semimajor axis ratios.

Zoom In Zoom Out Reset image size

In the CRTBP, some regions of the parameter space are more chaotic; these appear once the ZVC opens at the Lagrange point L₁, which lies between the stars (Quarles et al. 2011). Additional chaos can arise when the other collinear Lagrange points no longer lie within the forbidden region, due to changes in the mass ratio μ or the semimajor axis ratio a_p/a_bin. The origin of chaos in the CRTBP lies in the overlap of the N:1 MMRs with the binary, once the planetary eccentricity is substantially excited by a resonance (Murray & Dermott 1999; Mardling 2008). We take a wide view of the possible initial conditions (μ, a_p/a_bin, and λ_p) to produce probabilistic maps for a range of planetary inclinations in Figure 4. The probability for stability is calculated simply as the percentage of simulations (out of 91 values for λ_p) that survive for 500,000 yr for a given mass ratio and semimajor axis ratio. As a result, we identify regions that are stable (black), unstable (white), and quasi-stable (light or dark gray). The quasi-stable region is split into two regimes where either 1–18 trials survive (1%–20%; light gray) or 19–90 trials survive (20%–99%; dark gray) to delineate between chaotic zones and possible resonant islands that strongly depend on λ_p. In Figure 4, each panel also marks the critical boundaries (red curves) for which the ZVC opens for a collinear Lagrange point (L₁, L₂, or L₃) as detailed in Eberle et al. (2008) using the respective Jacobi constant (C₁, C₂, or C₃).

Zoom In Zoom Out Reset image size

Eberle et al. (2008) and Quarles et al. (2011) explored the coplanar case (Figure 4(a)) in detail using the Jacobi constant and maximum Lyapunov component, respectively, for planets that orbit the more massive primary (μ ≤ 0.5). However, Figure 4(a) shows that using the Jacobi constant C₁ as a stability criterion agrees well for all μ and only a few stable (black) regions exist to the right of the solid red curve, C₁. Between the curves for C₁ and C₃, there are initial conditions that can be stable depending on the initial mean longitude due to MMRs (Morais & Giuppone 2012; Morais & Namouni 2013), as indicated by the quasi-stable (light and dark gray) region, and Quarles et al. (2011) identified this regime as chaotic using the maximum Lyapunov exponent. All of the initial values to the right of C₃ (dashed red curve) are unstable due to large excitations of the planetary eccentricity by the perturbing star and the opening of all the collinear Lagrange points in the ZVC.

After the planet is inclined 30° relative to the binary orbit (Figure 4(b)), the stable (black) region does not change much. Panels (c) and (d) of Figure 4 illustrate significant differences from the coplanar cases due to the LK mechanism, where a large periodic eccentricity oscillation occurs that erodes much of the parameter space between the C₁ and C₃ curves. The stable regime for Figure 4(d) lies significantly below the C₁ curve, demonstrating the effectiveness of the LK mechanism. Some resonant regions (i.e., N:1 MMRs) can remain stable for longer periods, but must maintain a special configuration (e.g., similar to Pluto's 3:2 MMR with Neptune).

Panels (e)–(h) of Figure 4 show how the stability regions vary for retrograde orbits (${i}_{p}^{\mathrm{retro}}=180^\circ -{i}_{p}$), where the stable retrograde initial conditions (black) in Figure 4(e) represent a much larger fraction of the parameter space and do not agree with the Jacobi constant stability criterion from Eberle et al. (2008). The extent of the quasi-stable regime (light or dark gray) also increases, but it reaches a limit at the 3:2 MMR with the binary orbit. Figure 4(f) can be largely constrained with the Jacobi constant, and appears very similar in extent to Figure 4(b). Panels (g) and (h) in Figure 4 are similar to their prograde counterparts (Figure 4(c) and (d)), but the stabilizing effects from starting in retrograde increase the extent of the quasi-stable region at high μ and a_p/a_bin.

Another special circumstance can arise when the planet begins with a significant orbital velocity out of the binary plane so that the Coriolis force is minimized. In this scenario, the initial orbital velocity is calculated assuming the lower-mass star is the host (i.e., high μ), but the high semimajor axis ratio a_p/a_bin starts the planet within the Hill sphere of the more massive star. As a result, the planet can be captured into N:1 MMRs by the primary. Although this setup is allowed (i.e., no physical laws broken), it could be very unlikely to occur within known pathways for planet formation, but we highlight it for completeness. Figure 5 illustrates the evolution of a particular initial condition in rotated coordinates that can lead to this type of capture, where the planet begins either coplanar (Figure 5(a)) or highly inclined (Figure 5(c)). To confirm a capture into the 3:1 MMR, we plot the evolution of the resonance angle ϕ_3:1, where there are points in the coplanar case (Figure 5(b)) that are ill-defined (i.e., vertical jumps) and the highly inclined simulation (Figure 5(d)) shows libration.

Zoom In Zoom Out Reset image size

Even in the simplified case of the CRTBP, there is not a simple expression to define the likelihood that a planetary orbit will be stable. The most straightforward prescription to estimate stability is N-body simulations on a case-by-case basis, but this process can be greatly expedited or avoided using our results. For a nearly coplanar planet within a circular binary (e.g., transiting an eclipsing binary), the stability limit can be estimated using the Jacobi constant at the first Lagrange point C₁. We provide an algorithm to determine this curve as follows:

• 1.  
  Choose a value of μ from a list.
• 2.  
  Approximate through a power series or numerically determine the Lagrange point L₁ using a given μ and a root finding function.
• 3.  
  Use the location of L₁ to determine the Jacobi constant C₁ using the equation for the zero velocity surface ($2U={C}_{J};$ see Chap. 3 of Murray & Dermott 1999)
• 4.  
  Use the value of C₁ from (3) and numerically solve Equation (14) from Eberle et al. (2008) using a root finding function.
• 5.  
  Go back to (1), until the end of the list is reached.

Table 2 lists values for the semimajor axis ratio and mass ratio, where openings of the ZVC occur at L₁ or L₃ using the above algorithm, and these values correspond to the red curves (C₁ or C₃, respectively) in Figure 4. For a moderately inclined planet within a binary, similar curves provide a good approximation for the stability limit, but a coplanar retrograde orbit can significantly expand the stability limit as shown by C₁^retro. Finally, we provide Python tools on GitHub to query the data set represented in Figure 4, and there is a routine that returns the probability for stability through grid interpolation for a given combination of μ and a_p/a_bin.

Planets orbiting a single star of a circular binary represent only a subset of all the possibilities—and thus, we extend our method to investigate binary systems with a significant eccentricity (e_bin ≤ 0.8). This extra dimension expands the volume of our parameter space by almost two orders of magnitude, and thus we perform our numerical simulations with computational efficiency in mind. In Section 3.1, we explore eight different mutual inclinations because the eccentricity forcing on the planet is mainly due to MMRs. For eccentric binaries, we limit our investigation to only four mutual inclinations (i_p = 30°, 45°, and 180°), as the LK mechanism greatly reduces the stability limit within the 40–140° range of inclination (Quarles & Lissauer 2016; Quarles et al. 2018a).

Our coplanar simulations (i_p = 0°) use the uniform stepping described in Section 2.2 to determine the stability boundary a_c. To make the exploration of the inclined cased more efficient, we use the method of Lam & Kipping (2018), where they focused on a range in semimajor axis ratio (i.e., window) that would most likely contain the stability limit. To do this, we use the results from the coplanar runs to determine an appropriate trial window for the other inclined cases (i_p = 30°, 45°, and 180°). We then numerically investigate for longer timescales (500,000 yr) using semimajor axis ratios within $0.5{a}_{c}^{0}\mbox{--}1.5{a}_{c}^{0}$, where a_c⁰ is the stability limit in the respective prograde, coplanar simulation. If all the simulations in this range $0.5{a}_{c}^{0}\mbox{--}1.5{a}_{c}^{0}$ are unstable, then we explore a range of smaller semimajor axis ratios ($0.001\leqslant {a}_{p}/{a}_{\mathrm{bin}}\leqslant 0.5{a}_{c}^{0}$). This procedure was sufficient to cover regions of parameter space that we expect to either be truncated (e.g., Lidov–Kozai mechanism for ${i}_{p}=45^\circ$) or expanded (e.g., retrograde orbits for i_p = 180°).

Following Holman & Wiegert (1999), we take these results and compare to previous studies (Rabl & Dvorak 1988) that prescribe a stability limit for equal mass stars (μ = 0.5). Figure 6 illustrates how our results compare for the special case of equal mass stars with a coplanar planet, the range of stability for two regimes of μ (red and blue), and the changes to the median stability limit within each regime for each planetary inclination. Rabl & Dvorak (1988) performed simulations for only 300 binary periods and explored up to e_bin = 0.6. Thus, we expect for deviations to appear for high binary eccentricities (0.6 ≤ e_bin ≤ 0.8 where the expression from Rabl & Dvorak (1988) is extrapolated as seen in Figure 6(a). There is also a small difference between the curve from Holman & Wiegert (1999) (black dashed) and our results (solid gray) in the high binary eccentricity region, as well as in the e_bin ≤ 0.1 region. These differences are due to the factor of 10 increase in sampling that we perform, and amount to slight changes in the stability limit for a coplanar planet. The shaded regions (red and blue) in Figure 6(a) demonstrate the extent to which the stability limit can vary as a function of μ, where μ = 0.01 follows the upper red boundary and μ = 0.99 is much flatter along the lower blue boundary. As in the case for the CRTBP, we can see that the stability limit is not symmetric about the equal mass case.

Zoom In Zoom Out Reset image size

The planetary inclination also plays a role in determining the stability limit a_c, where Figure 6(b) shows these differences through quadratic fits to the median stability limit for each e_bin. The i_p = 0° (solid) and i_p = 30° (dashed) curves are nearly identical in Figure 6(b), and there is not a strong dependence on μ when split into the two regimes (red and blue). As a result, the stability limit for coplanar planets will be similar enough, in most cases, to be valid up to ∼40°, until the Lidov–Kozai mechanism becomes active. The same will be true for retrograde planetary inclinations ($140^\circ \leqslant {i}_{p}\lt 180^\circ$) relative the exactly retrograde case (i_p = 180°). Sufficiently inclined planets (i_p = 45°; dashed–dotted in Figure 6(b)) where the Lidov–Kozai mechanism is active have a lower median stability limit, due to eccentricity excitation, and it is largely independent of μ. On the other hand, planets in retrograde orbits (dotted in Figure 6(b)) have a larger median stability limit.

Table 3 provides coefficients and uncertainties for the quadratic fitting formula (${a}_{c}/{a}_{\mathrm{bin}}={c}_{1}+{c}_{2}{e}_{\mathrm{bin}}+{c}_{3}{e}_{\mathrm{bin}}^{2}$) from Rabl & Dvorak (1988), Holman & Wiegert (1999), and using the stability limits determined through our simulations. The difference between our coplanar results (i_p = 0°) and the previous studies is that we break the fitting formula into two regimes for μ, as motivated by the CRTBP and Figure 6. If we find the average value for each coefficient (c₁ − c₃), then our results are consistent with the previous works. Our c₁ and c₂ coefficients for i_p = 30° are similar to our coplanar coefficients. The coefficients for the i_p = 45° fitting formula depend more strongly on the binary eccentricity due to the Lidov–Kozai mechanism. The coefficients for the retrograde fitting formulas are larger than the coplanar, prograde formulas by more than 30%. Moreover, each set of coefficients reveal that the stability limit could be up to two times bigger for planets that orbit the more massive primary over the less massive secondary star.

Table 3.  Critical Semimajor Axis as a Function of the Binary Eccentricity e_bin

	ip	c1 ± σ1	c2 ± σ2	c3 ± σ3
$0.1\leqslant \mu \leqslant 0.9\dagger$	0°	0.262 ± 0.006	−0.254 ± 0.017	0.060 ± 0.027
$0.1\leqslant \mu \leqslant {0.9}^{\star }$	0°	0.274 ± 0.008	−0.338 ± 0.045	0.051 ± 0.055
$0.01\leqslant \mu \leqslant 0.5$	0°	0.363 ± 0.001	−0.492 ± 0.006	0.129 ± 0.005
$0.51\leqslant \mu \leqslant 0.99$	0°	0.186 ± 0.001	−0.193 ± 0.003	0.001 ± 0.003
$0.01\leqslant \mu \leqslant 0.5$	30°	0.346 ± 0.002	−0.464 ± 0.006	0.117 ± 0.005
$0.51\leqslant \mu \leqslant 0.99$	30°	0.198 ± 0.001	−0.243 ± 0.002	0.043 ± 0.002
$0.01\leqslant \mu \leqslant 0.5$	45°	0.247 ± 0.005	−0.487 ± 0.018	0.268 ± 0.016
$0.51\leqslant \mu \leqslant 0.99$	45°	0.213 ± 0.002	−0.441 ± 0.009	0.252 ± 0.008
$0.01\leqslant \mu \leqslant 0.5$	180°	0.479 ± 0.002	−0.647 ± 0.008	0.168 ± 0.008
$0.51\leqslant \mu \leqslant 0.99$	180°	0.298 ± 0.001	−0.378 ± 0.006	0.072 ± 0.006

Note. The listed coefficients (${c}_{1}-{c}_{3}$ and uncertainties (${\sigma }_{1}-{\sigma }_{3}$) from $\dagger$Rabl & Dvorak (1988) and ${}^{\star }$Holman & Wiegert (1999) use a quadratic fitting function, ${a}_{c}/{a}_{\mathrm{bin}}={c}_{1}+{c}_{2}{e}_{\mathrm{bin}}+{c}_{3}{e}_{\mathrm{bin}}^{2}$. We use the same function in this work, but we split it between two domains in μ for each planetary inclination i_p.

Download table as:  ASCIITypeset image

An accurate and reliable fitting formula should include a dependence on the mass ratio μ, given the asymmetries we have described in the CRTBP and in Figure 6. Holman & Wiegert (1999) assumed that the stability limit depends linearly on μ and quadratically for e_bin. We follow the same prescription and provide updated values for the coefficients ${c}_{1}-{c}_{6}$, using our much larger data set in Table 4 along with the previous values from Holman & Wiegert (1999) for comparison. Moreover, coefficients are listed for inclined planets (i_p = 30°, 45°, and 180°), where we see that the coefficients for the coplanar and 30° formulas are similar. When the planet is initially inclined to 45°, the c₃ − c₆ terms dominate due to the strong dependence on e_bin from the Lidov–Kozai mechanism. The c₁, c₃, and c₅ coefficients for retrograde are larger than the respective coefficients for i_p = 0° formula, which indicates that retrograde stability depends more strongly on e_bin, due to a lower magnitude forcing on the planet from the stellar companion (Henon 1970).

Table 4.  Coefficients for the Critical Semimajor Axis

	ip	${c}_{1}\pm {\sigma }_{1}$	${c}_{2}\pm {\sigma }_{2}$	${c}_{3}\pm {\sigma }_{3}$	${c}_{4}\pm {\sigma }_{4}$	${c}_{5}\pm {\sigma }_{5}$	${c}_{6}\pm {\sigma }_{6}$
Previous Work${}^{\star }$	0°	0.464 ± 0.006	−0.380 ± 0.010	−0.631 ± 0.034	0.586 ± 0.061	0.650 ± 0.041	−0.198 ± 0.074
This work	0°	0.501 ± 0.002	−0.435 ± 0.003	−0.668 ± 0.009	0.644 ± 0.015	0.152 ± 0.011	−0.196 ± 0.019
This work	30°	0.485 ± 0.002	−0.405 ± 0.003	−0.684 ± 0.009	0.603 ± 0.015	0.190 ± 0.011	−0.182 ± 0.019
This work	45°	0.428 ± 0.002	−0.318 ± 0.003	−1.128 ± 0.011	0.987 ± 0.020	0.839 ± 0.014	−0.825 ± 0.024
This work	180°	0.617 ± 0.001	−0.457 ± 0.002	−0.787 ± 0.008	0.586 ± 0.014	0.163 ± 0.010	−0.128 ± 0.017

Note. The coefficients (${c}_{1}-{c}_{6}$) and uncertainties (${\sigma }_{1}-{\sigma }_{6}$) from ${}^{\star }$Holman & Wiegert (1999) are listed using the fitting formula, ${a}_{c}/{a}_{\mathrm{bin}}={c}_{1}+{c}_{2}\mu \,+{c}_{3}{e}_{\mathrm{bin}}+{c}_{4}\mu {e}_{\mathrm{bin}}+{c}_{5}{e}_{\mathrm{bin}}^{2}+{c}_{6}\mu {e}_{\mathrm{bin}}^{2}$. We use the same function in this work, but also fit for a range of planetary inclination i_p.

Download table as:  ASCIITypeset image

The fitting formulas in Section 3.2 are dependent on the quality and breadth of the numerical simulations that underlie their derivation. They are reliable in characterizing a population of exoplanets, but they can also be inaccurate when describing particular systems. For instance, we can compare the case when e_bin = 0 in Table 4 to the results in Table 2. The estimate for the critical semimajor axis a_c from Holman & Wiegert (1999) and our own newly updated estimation is different by more than 10% for μ < 0.05. As a result, we suggest a different approach using a lookup table or interpolation map that relies less on statistical averaging, which occurs in deriving a single fitting formula.

Figure 7 shows⁷ the critical semimajor axis ratio a_c/a_bin (color-coded) over the full range of mass ratio μ and binary eccentricity e_bin. Table 4 shows only small differences in the fitting formulas for i_p = 0° and i_p = 30°, where this trend is visually apparent in Figure 7(a) and (b). In addition, panels (a) and (b) illustrate the asymmetry of the stability limit with respect to the mass ratio μ, where for μ > 0.5 the typical value is ${a}_{c}/{a}_{\mathrm{bin}}\lesssim 0.2$. As a result, we find the volume for a stable planet for the primary star to host a nearly planar (i_p ≤ 30°) planet in a low mass ratio (μ < 0.3) binary to be typically two times larger than for the lower-mass secondary (μ > 0.7). The propensity for stable orbits is substantially reduced for a more inclined (i_p = 45°) planet when e_bin ≳ 0.2 and the critical semimajor axis becomes more symmetric about μ = 0.5 (Figure 7(c)). A retrograde planet (Figure 7(d)) can be stable to much larger semimajor axis ratios and for a larger fraction of binary parameters (μ and e_bin), which has been shown to some degree in several other investigations (Henon 1970; Benest 1988b; Wiegert & Holman 1997; Quarles & Lissauer 2016).

Zoom In Zoom Out Reset image size

The results in Figure 7 illustrates the variations of the critical semimajor axis with respect to the mass ratio μ and eccentricity e_bin, but observations of binary star populations, especially with Sunlike primaries, are not uniform with respect to these binary parameters (Raghavan et al. 2010; Moe & Di Stefano 2017). Therefore, a statistical study is needed to gain further insights from our results using Monte Carlo approaches. David et al. (2003) and Fatuzzo et al. (2006) performed a statistical study with the binary periastron distance as a stability criterion, and using the contemporary stability formula (Holman & Wiegert 1999) and a survey of binary stars (Duquennoy & Mayor 1991). We follow their approach, but use our our determination of the critical semimajor axis and more up-to-date results from binary star surveys (Raghavan et al. 2010; Moe & Di Stefano 2017). This is important for observations because, when planets are discovered orbiting a star with a distant stellar companion, e.g., K2-288b (Feinstein et al. 2019), the binary period and eccentricity are often unknown, but the planetary semimajor axis and/or the projected binary semimajor axis can be determined (see Table 5 for a list of known planet hosting binary stars). Moreover, observations of young binary stars can reveal the stability limit through dust and gas emission around either star (e.g., Alves et al. 2019), but this is also a projected distance.

Table 5.  Mass Ratio, Semimajor Axis, and Stability Limits for Known Planet-hosting Binaries with $5\leqslant {a}_{\mathrm{bin}}\lt 500\,\mathrm{au}$

System	MA	MB	μ	abin	${a}_{c}^{A}$,a	${a}_{c}^{A}$,b	${a}_{c}^{B}$,a	${a}_{c}^{B}$,b
	(${M}_{\odot }$)	(${M}_{\odot }$)		(au)	(abin)	(abin)	(abin)	(abin)
HD 109749	1.100	0.780	0.415	490	0.046	0.299	0.040	0.235
HD 133131	0.950	0.930	0.495	360	0.043	0.267	0.042	0.263
HD 106515	0.910	0.880	0.492	345	0.043	0.267	0.042	0.261
Kepler-108	1.300	0.960	0.425	327	0.045	0.295	0.040	0.238
WASP-77	1.002	0.710	0.415	306	0.046	0.299	0.040	0.235
KELT-2	1.317	0.780	0.372	295	0.047	0.321	0.037	0.222
HD 114729	0.930	0.253	0.214	282	0.053	0.379	0.031	0.169
Kepler-14	1.510	1.390	0.479	280	0.043	0.272	0.042	0.257
HD 27442	1.200	0.750	0.385	240	0.047	0.314	0.038	0.227
TrES-2	0.980	0.509	0.342	232	0.048	0.338	0.037	0.213
HD 212301	1.270	0.350	0.216	230	0.053	0.380	0.032	0.170
HD 16141	1.010	0.286	0.221	220	0.052	0.380	0.032	0.171
HD 189733	0.800	0.200	0.200	216	0.054	0.377	0.031	0.164
HD 217786	1.020	0.160	0.136	155	0.057	0.435	0.027	0.141
HD 142	1.200	0.590	0.330	138	0.048	0.346	0.036	0.209
HD 114762	0.840	0.138	0.141	132	0.057	0.432	0.028	0.143
HD 195019	1.060	0.700	0.398	131	0.046	0.307	0.040	0.230
WASP-2	0.890	0.480	0.350	106	0.047	0.334	0.037	0.215
HD 19994	1.340	0.350	0.207	100	0.054	0.378	0.031	0.167
HD 177830	1.450	0.230	0.137	97	0.057	0.434	0.027	0.142
Gliese 15	0.380	0.150	0.283	93	0.051	0.378	0.035	0.194
Kepler-296	0.500	0.330	0.398	80	0.046	0.307	0.040	0.230
GJ 3021	0.900	0.130	0.126	68	0.057	0.440	0.026	0.137
K2-288	0.520	0.330	0.388	54.8	0.047	0.312	0.038	0.227
HD 120136	1.400	0.400	0.222	45	0.052	0.380	0.032	0.172
WASP-11	0.820	0.340	0.293	42	0.050	0.369	0.035	0.197
K2-136	0.740	0.100	0.119	40	0.057	0.442	0.026	0.135
HD 164509	1.130	0.420	0.271	37	0.051	0.383	0.034	0.190
HD 41004	0.700	0.400	0.363	23	0.048	0.326	0.037	0.219
HD 196885	1.330	0.550	0.293	23	0.050	0.369	0.035	0.197
HD 4113	0.990	0.550	0.357	23	0.048	0.330	0.037	0.217
GJ 86	0.800	0.490	0.380	21	0.047	0.317	0.038	0.225
γ Cep	1.180	0.320	0.213	20	0.053	0.379	0.031	0.169
HD 8673	1.400	0.400	0.222	10	0.052	0.380	0.032	0.172
Kepler-420	0.990	0.700	0.414	5.3	0.046	0.300	0.039	0.235

Notes. Critical semimajor axis a_c calculated using interpolation of Figure 7(a), where planets orbit star A for μ and star B for $1-\mu$. The stellar parameters are from an online catalog of planets in binaries maintained by Richard Schwarz (https://www.univie.ac.at/adg/schwarz/multiple.html).

^a ${e}_{\mathrm{bin}}=0.8.$ ^b ${e}_{\mathrm{bin}}=0.0.$

Download table as:  ASCIITypeset image

Our approach is to determine a probability distribution function (PDF) for the stability of planets in binaries as a function of the critical semimajor axis ratio a_c/a_bin, where the sum of the binary masses is equal to a solar mass (i.e., M_A + M_B = 1 M_⊙). The binary period (and semimajor axis) follows a log-normal probability distribution (${p}_{\mathrm{log}P}\propto {e}^{-{(\mathrm{log}P-\zeta )}^{2}/2{\sigma }^{2}}$), where ζ = 5.03 and σ = 2.28 (Raghavan et al. 2010). This probability distribution in Raghavan et al. (2010) has a broad range (−2 ≤ log P ≤ 10, in days), where we limit our analysis to a smaller range in binary period (4 ≤ log P ≤ 7) that corresponds to a reasonable range in binary semimajor axis (10 au ≤ a_bin ≤ 1000 au).

Observers typically use the binary mass quotient q (=M_B/M_A, where M_B ≤ M_A), which algebraically relates to the dynamical mass ratio μ ($=q/(1+q)$). From the mass quotient q, the range for μ is limited to 0.5 and our exploration of planets orbiting the secondary star (μ > 0.5) is accomplished by subtracting the μ samples from unity (i.e., 1 − μ). We use a probability distribution (${p}_{q}\propto {q}^{\gamma }$) for q derived by Moe & Di Stefano (2017) that is a broken power law, which has γ₁ = 0.3 ± 0.4 for 0.1 ≤ q ≤ 0.3 and γ₂ = −0.5 ± 0.3 for 0.3 ≤ q ≤ 1. An excess twin fraction F_twin = 0.1 ± 0.3 is added to the PDF for q ≥ 0.95; see Figure 2 of Moe & Di Stefano (2017). A probability distribution (${p}_{e}\propto {e}_{\mathrm{bin}}^{\eta }$) for the binary eccentricity follows a single power law, which has η = 0.4 ± 0.3 (Moe & Di Stefano 2017). We note that the exponents used in these power laws are dependent on the range of sampled binary period log P (Moe & Di Stefano 2017) and adjustments to the exponents are necessary if one expands our analysis to a broader range in binary period.

Using these probability distributions, we employ Monte Carlo integration to determine the fraction of observed binaries with a critical semimajor axis larger than a given a_c/a_bin using the data in panels (a), (c), and (d) of Figure 7. We do not use Figure 7(b) because of its strong similarity to Figure 7(a). Figure 7(a) illustrates that all binaries within the domain will have a critical semimajor axis larger than 0.001 and very few binaries will have a a_c/a_bin > 0.5 where this is independent of the weighting given by the probability distributions we choose. As a result, we integrate a_c/a_bin = 0.001–0.01 using 0.001 sized bins, and a_c/a_bin = 0.01–1.0 using 0.01 sized bins. Figure 10 shows the PDF (Figure 10(a)–(c)) and cumulative distribution function (CDF; Figure 10(d)–(f)) resulting from the Monte Carlo integration considering planetary inclinations of 0°, 45°, and 180°. For planets orbiting either the primary (star A; black) or the secondary (star B; gray) star, there are many binaries with a_c/a_bin > 0.05, and hence both curves in Figure 10(a)–(c) have the highest probability within this range. Planets that orbit the secondary constitute a larger fraction of binaries with a small critical semimajor axis, and thus the probability is slightly larger. This trend reverses for a_c/a_bin ≳ 0.12, as the stability limit for the primary extends to larger distances compared to the secondary. Sufficiently inclined planets (e.g., i_p = 45°) have their stability limits truncated due to the large growth of eccentricity from the Lidov–Kozai mechanism, and stable orbits are possible for predominantly low critical semimajor axis ratios. Retrograde orbits extend to larger critical semimajor axis ratios (Henon 1970; Quarles & Lissauer 2016), and the PDF is flatter.

Zoom In Zoom Out Reset image size

The uncertainty in the binary semimajor axis a_bin from observations can be quite large, where discoveries from direct imaging may only have an estimate for the projected semimajor axis. From such situations, a CDF is more useful, where the probability can be calculated from a difference of values on the CDF and a single value describes the fraction of binaries with a stability limit less than a given a_c/a_bin. For example, Figure 10(d) shows that half of all binaries within our domain will have a critical semimajor axis for coplanar planets that is less than 0.06 (star B) or 0.078 (star A). The slope in the CDF for the Lidov–Kozai regime (Figure 10(e)) is steeper, which shifts the median to smaller a_c/a_bin. Moreover, the slope is nearly identical between whichever star the planet is orbiting (star A or B). The shallower slope for retrograde planets (Figure 10(f)) has a median critical semimajor axis closer to 0.1 and extends to larger values of a_c/a_bin.

The black and gray curves in Figure 10(d)–(f) are fit well by a modified exponential distribution:

$$\begin{eqnarray}&&F(\xi )=1-{e}^{-{c}_{1}{\xi }^{2}-{c}_{2}\xi },\end{eqnarray} \tag{ 1 }$$

where ξ = a_c/a_bin. The red curves in Figure 10(d)–(f) mark sample curves within the 1σ of the best-fit values for F(ξ), where the largest differences occur in the top quintile as marked by the respective inset panels. David et al. (2003) used a similar function, except their formulation depends on the binary periastron distance. Table 6 provides values for the coefficients (c₁ and c₂), along with their uncertainties using the Levenberg–Marquardt algorithm in the curve fit function within the scipy module. In addition to the coefficients, we provide the the critical semimajor axis values at 50% and 99.7% quantiles.

Table 6.  Cumulative Distribution Function Coefficients for the Critical Semimajor Axis

	ip	${c}_{1}\pm {\sigma }_{1}$	${c}_{2}\pm {\sigma }_{2}$	50%	99.7%
Star A	0°	28.96 ± 0.75	6.65 ± 0.09	0.078	0.348
Star B	0°	54.50 ± 1.62	8.31 ± 0.14	0.060	0.259
Star A	45°	15.11 ± 2.30	13.28 ± 0.23	0.049	0.321
Star B	45°	20.69 ± 3.28	13.82 ± 0.30	0.047	0.292
Star A	180°	15.41 ± 0.37	4.82 ± 0.06	0.107	0.477
Star B	180°	22.97 ± 0.62	5.88 ± 0.08	0.088	0.391

Note. The coefficients (c₁ and c₂) and uncertainties (${\sigma }_{1}$ and ${\sigma }_{2}$) using a fitting formula for the CDF, $F(\xi )=1-{e}^{-{c}_{1}{\xi }^{2}-{c}_{2}\xi }$ where $\xi ={a}_{c}/{a}_{\mathrm{bin}}$.

Download table as:  ASCIITypeset image

NASA's TESS will observe nearly 85% of the sky over the course of its primary two-year mission. The spacecraft has four cameras that provide a $24^\circ \times 96^\circ$ field of view, and a photometric precision between  ∼100 ppm for apparent magnitude ${I}_{C}\sim 4$, and $\sim {10}^{5}\ \mathrm{ppm}$ for ${I}_{C}\sim 18$ (Ricker et al. 2015). It is expected that TESS will discover thousands of transiting exoplanets around bright, nearby stars (Sullivan et al. 2015; Bouma et al. 2017; Barclay et al. 2018), including hundreds of planets smaller than $2{{\rm{R}}}_{\oplus }$. A significant number of stars (some of which will have transiting planets) observed by TESS will be in multiple systems, i.e., from $\sim 20 \% \ \mathrm{for}\ {M}_{\mathrm{prim}}\,\lt 0.1\ {M}_{\odot }$ to $\sim 75 \% \ \mathrm{for}\ {M}_{\mathrm{prim}}\gt 1.4\ {M}_{\odot }$, and spanning semimajor axes from a few to hundreds of au (Sullivan et al. 2015). In combination with data from Gaia, by virtue of being bright targets (and thus favorable to follow-up observations) those of the multiple star systems hosting S-type planet candidates will be amenable to measurements of the mass ratio. The results we present here will be directly applicable to the stability of these systems.

Recently, Winters et al. (2019) uncovered a transiting planet in a triple star system, LTT 1445ABC, where star A (M_A = 0.257 ${\text{}}{M}_{\odot }$) orbits a nearly equal-mass stellar binary (M_B = 0.215 ${\text{}}{M}_{\odot }$ and M_C = 0.161 ${\text{}}{M}_{\odot }$) and the separation from star A to the BC barycenter is ∼34 au. The nearly Earth-sized exoplanet LTT 1445Ab (${R}_{p}\approx 1.38{R}_{\oplus }$) orbits very quickly (5.358 days, or a_p = 0.038 au) around LTT 1445A, which is well within the stability limit, assuming that the orbit of LTT 1445A is nearly circular around the barycenter of LTT 1445BC. However, the A–BC orbit is likely eccentric, considering observations over the last few decades that indicate significant astrometric variation (Mason et al. 2001; Dieterich et al. 2012; Winters et al. 2019, and references therein), and a noncircular orbit may affect the stability of the observed exoplanet. Using the combined mass for LTT 1445BC (M_BC = 0.376 M_⊙), the dynamical mass ratio of the A–BC pair is μ = 0.594, given that the planet is orbiting the less massive component of the pair. The semimajor axis a_p is quite small (0.038 au), where the planetary orbital stability may be compromised for a large eccentricity e_bin of the A–BC orbit.

If we consider planets within the A–BC plane, then the stability limit varies between 1.38 and 8.15 au depending on the value of e_bin, the eccentricity of star A around the BC barycenter. Even for a high value for e_bin, the stability limit is still ∼36 times farther from the host star than the observed planetary semimajor axis. Therefore, we are confident that LTT 1445Ab is on a stable orbit, as the planet is very far from the stability limit under reasonably extreme assumptions on the shape of the LTT 1445 A–BC stellar orbit. If TESS observes the system during the extended mission and detects potential transits of additional planets on longer orbits, reducing the BC binary to its barycenter may need to be revised and the eccentricity of star A around the BC barycenter properly evaluated. We do not investigate nonplanar conditions (i_p = 45°) where the stability limit can be further reduced, because the astrometric measurements show that the LTT 1445ABC system is nearly planar and the planet is transiting—which implies that the planetary orbit is likely nearly aligned with the LTT 1445ABC orbital plane.