Stability Limits of Circumbinary Planets: Is There a Pile-up in the Kepler CBPs?

Our simulations use a modified scheme within the popular mercury integration package (Chambers et al. 2002) that has been designed for the efficient simulation of circumbinary systems. This modification allows for the integration of the inner binary orbit and an outer planet at different timescales while preserving the symplectic nature of the integration method. As a result, we find that the largest integration step for the planet to be ∼2.5% of the planetary Keplerian period. Our numerical scheme stops the simulation when an instability event occurs, which we define as an intersection with the binary orbit or when the radial distance of the planet to the more massive star exceeds 10 au. Using 10 au from the more massive star as a distance cutoff is justified because the planets begin with small semimajor axes, which rules out such a large apastron distance, and planets that reach this distance are likely to exceed the respective escape velocity.

A majority of our simulations use ideal initial conditions for the binary orbit, where the binary semimajor axis is 1 au and the total mass (M_A + M_B) of the stellar components is 1 ${{\rm{M}}}_{\odot }$. Our runs consider a range of binary mass ratios (μ = M_B/(M_A + M_B)) from 0.01 to 0.5 in steps of 0.01 and include one additional case, μ = 0.001, for a total of 51 steps. The eccentricity of the binary orbit varies from 0.0 to 0.80 in steps of 0.01. Most of our integrations begin the binary orbit at periastron (λ_bin = 0°) because previous investigations (HW99) have shown this assumption to produce a conservative estimate for the stability limit. However, we do investigate a small subset of runs to quantify how beginning the binary at apastron (λ = 180°) would change our results.

In order to compare with previous results (i.e., HW99), we perform integrations to determine a critical semimajor axis a_c for a Jupiter-mass planet that begins on an initially coplanar, circular orbit. We define the stability limit as the the critical semimajor axis ratio a_c in units of the binary semimajor axis a_bin and measure it by the smallest planetary semimajor axis where a planet is stable for all choices of initial mean anomaly or phases (i.e., the lower critical orbit, Dvorak et al. 1989) relative to the host binary orbit. This definition is motivated by our models of gas giant CBPs that employ migration from a larger distance through interactions with the disk (Pierens & Nelson 2008; Paardekooper et al. 2012; Kley & Haghighipour 2014, 2015; Meschiari 2014; Bromley & Kenyon 2015). Such studies have shown that gas disk migration around tight binaries can occur in a similar manner as in single star systems, where gas drag acts to circularize planetary orbits. After the gas disk dissipates, the binary excites the eccentricities of close-in exoplanets leading to scattering events or expulsion from the system (e.g., Kley & Haghighipour 2015; Silsbee & Rafikov 2015; Thebault & Haghighipour 2015; Vartanyan et al. 2016).

In order to determine a_c consistently, given the above definition, we perform simulations over a grid of orbital parameters, where the total simulation time per integration is 10⁵ binary orbits. Our grid of planetary orbital parameters, for each combination of the stellar μ and e_bin, vary the semimajor axis ratio (a_p/a_bin) from 1.01 to 5.0 in steps of 0.01 and the initial planetary mean anomaly from 0° to 180° in steps of 2°.

HW99 performed their simulations using a similar definition for the stability limit but only include eight initial phases (0°–315°) for each test particle. Our simulations take advantage of a symmetry with respect to the initial phase and increase the number of trial phases because some initial values can become unstable in between the 45° increments employed by HW99. A higher resolution is necessary to ensure that our definition for the stability limit is reliable.

For very small values of both μ and e_bin, our simulations approach conditions consistent with Hill stability (Szebehely & McKenzie 1981; Gladman 1993). This type of analysis identifies the gravitational radius of influence that the secondary mass has on the tertiary mass through the Hill radius R_H = a(μ/3)^1/3.

Parameterization using the Hill radius has also been used in the stability of planetary systems around single stars, but has been modified slightly to include the average semimajor axis between adjacent pairs of planets and the mass interior to the outermost body (Chambers et al. 1996). Using this formalism, we measure the dynamical fullness of the system. The mathematical definitions of the mutual Hill Radius, R_H,m, and dynamical spacing, β, between the k,k + 1 planets with identical mass, m, are:

$$\begin{eqnarray}&&{R}_{H,m}=\displaystyle \frac{1}{2}({a}_{k}+{a}_{k+1}){\left(\displaystyle \frac{2m}{3{M}_{T}}\right)}^{1/3},\,\mathrm{and}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\beta =\displaystyle \frac{{a}_{k+1}-{a}_{k}}{{R}_{H,m}},\end{eqnarray} \tag{ 2 }$$

where M_T = M_A + M_B + m and represents the total mass interior to the outermost body. Our analysis uses this representation to evaluate whether a planet of identical mass could be placed at a_c in addition to the observed CBP. Others have investigated planet packing for CBPs and determined that values of β = 5–7 would represent the minimum dynamical spacing necessary near the stability limit (Kratter & Shannon 2014; Andrade-Ines & Robutel 2018). We use this formalism to measure the dynamical fullness of each system, where dynamically full systems relative to the stability limit are potential evidence for a pile-up in the CBPs.

We use numerical simulations to identify the stability limit for a Jupiter-mass planet in an initially circular, coplanar orbit around a range of binary parameters. The results of such simulations can vary with the assumed binary orbit. Therefore, we first identify the range of variation we can expect when the binary begins at either periastron (λ_bin = 0°) or apastron (λ_bin = 180°).

Figure 1 illustrates how the maximum eccentricity (color-coded) of the planet varies with respect to the initial semimajor axis ratio (a_p/a_bin) and planetary mean anomaly when the binary begins at periastron. Each subplot varies the binary stellar parameters (μ, e_bin), where the respective values are given in the upper right corner. Additionally, the stability limit a_c is identified by a horizontal cyan line, and the value is given in the lower left corner. The white space denotes regions of parameter space that are unstable on the timescale of 10⁵ binary orbits.

Zoom In Zoom Out Reset image size

The stability limit in the top row of Figure 1 increases as the binary eccentricity, e_bin, increases. There are also increases in a_c as μ increases (i.e., starting from the top row and going down a given column). However, the largest value in a_c does not occur at both the largest μ and e_bin combination, rather when e_bin is large (0.5) and μ is modest (0.1–0.3). Also, the stability islands—that depend on mean motion resonances—are symmetric about 180^° in the planetary mean anomaly. Deck et al. (2013) found similar results when investigating first-order resonance overlap in close two planet systems around a single star, where a similar dynamical environment exists. The symmetry justifies our choice to investigate only from 0° to 180° in the initial planetary mean anomaly in the more computationally expensive portion of our study (See Section 3).

In contrast, Figure 2 demonstrates how the stability limit a_c changes when the binary begins at apastron (λ_bin = 180°). Similar trends are present and most changes in a_c between respective subplots are small (<0.1). The most drastic change occurs when μ = 0.5 and e_bin = 0.5, where the difference in a_c is 0.37. In order to produce conservative (and possibly more reliable) results, we begin the binary at periastron for simulations used in the rest of our analysis.

Zoom In Zoom Out Reset image size

We perform a multitude of simulations⁵ (∼150 million) to improve the accuracy of the stability limit, a_c, for CBPs. We determine the a_c for a given combination of binary parameters, μ and e_bin, using a grid of simulations (e.g., Figures 1 and 2), where we limit initial planetary mean anomaly to 180^° (see Section 2.4).

From these results, we analyze how a_c varies at a given binary eccentricity, e_bin, as a function of the binary mass ratio. Figure 3 shows these results where the color-code represents the binary mass ratio, μ, and the smallest value μ = 0.001 is excluded because it is significantly flatter relative to the rest of the points. When including this broad range of mass ratio, the value of a_c, can typically vary by ∼0.5–1.0 a_bin. This variation changes with the binary eccentricity and the median value is not proportional to the mean value. We have over-plotted the median value (black points) with error bars indicating the upper and lower extremes in Figure 3 to illustrate how the stability limit is affected by the binary mass ratio, μ.

Zoom In Zoom Out Reset image size

Most of the variation in the lower bound occurs for μ < 0.1. If points with μ < 0.1 were excluded, then a polynomial function could approximately reflect a_c statistically. HW99 included results where μ = 0.1–0.5 and determined a quadratic polynomial to be appropriate, although their selection on μ was to exclude a regime that can be modeled using Hill stability.

In order to make a fair comparison with HW99, we plot in Figure 4 the median values of the stability limit, a_c, using error bars to indicate the total range (maximum/minimum values) at a given binary eccentricity e_bin. In addition to the data points (blue), we also include the curves using the respective coefficients from Dvorak et al. (1989) (solid black), HW99 (dashed black), and those determined from our simulations (solid red). The coefficients and uncertainties are provided in Table 1, where the values for Dvorak et al. (1989) and HW99 are both quoted from HW99 due to the possible errors in labeling noted by HW99.

Zoom In Zoom Out Reset image size

Table 1.  Coefficients for the Critical Semimajor Axis Using e_bin

	C1	C2	C3
Dvorak et al. (1989)	2.37	2.76	−1.04
HW99	2.278+0.008−0.008	${3.824}_{-0.33}^{+0.33}$	−1.71${}_{-0.10}^{+0.10}$
this work	${2.170}_{-0.017}^{+0.017}$	${4.017}_{-0.10}^{+0.10}$	−1.75${}_{-0.14}^{+0.14}$

Note. The coefficients (and uncertainties) for C₁, C₂, and C₃ from previous studies are listed that use a quadratic fitting function ignoring μ, ${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 use the maximum, median, and minimum values of a_c (e.g., Figure 4).

Download table as:  ASCIITypeset image

Upon inspection of Figure 4 and Table 1, we reaffirm the previous results, where most of our coefficients overlap (within errors) with those of HW99. However, both fits are applicable at a statistical distribution level and not very accurate individually due to the effective marginalization over the binary mass ratio, μ. In Figure 4, we also mark the expected locations of the mean motion resonances the planet would encounter with the binary orbit, which act to destabilize CBPs (Mudryk & Wu 2006). Doolin & Blundell (2011), Quarles & Lissauer (2016), and others have shown through large parameter space studies that these resonances produce unstable gaps and stability islands can exist at locations approximately half-way between the resonances.

Another method utilized by HW99 is to allow both the binary semimajor axis, μ, and the binary eccentricity, e_bin, to vary as quadratic functions. We perform a similar approach (using all our data) and provide our results alongside those determined by HW99 in Table 2. The reduced chi-square statistic is provided using the HW99 coefficients as well as our own. We provide an additional fitting where we make the replacement of $\mu \,\to \,{\mu }^{1/3}$ as motivated by stability studies using the Hill radius (i.e., planet packing around a single star). Both of our fittings produce a lower chi-square statistic than HW99, although our Fit 2 is likely to be biased in that a large portion of our simulations (∼40%) have a binary mass ratio within the Hill regime.

Table 2.  Coefficients for the Critical Semimajor Axis

	C1	C2	C3	C4	C5	C6	C7	${\chi }_{\nu }^{2}$
HW99	${1.60}_{-0.04}^{+0.04}$	${5.10}_{-0.05}^{+0.05}$	−2.22${}_{-0.11}^{+0.11}$	${4.12}_{-0.09}^{+0.09}$	−4.27${}_{-0.17}^{+0.17}$	−5.09${}_{-0.11}^{+0.11}$	${4.61}_{-0.36}^{+0.36}$	2015.97a
Fit 1	${1.48}_{-0.01}^{+0.01}$	${3.92}_{-0.06}^{+0.06}$	−1.41${}_{-0.06}^{+0.06}$	${5.14}_{-0.10}^{+0.10}$	${0.33}_{-0.19}^{+0.19}$	−7.95${}_{-0.15}^{+0.15}$	−4.89${}_{-0.44}^{+0.44}$	876.25
Fit 2	${0.93}_{-0.02}^{+0.02}$	${2.67}_{-0.08}^{+0.08}$	−0.25${}_{-0.06}^{+0.06}$	${3.72}_{-0.06}^{+0.06}$	${2.25}_{-0.12}^{+0.12}$	−2.72${}_{-0.05}^{+0.05}$	−4.17${}_{-0.15}^{+0.15}$	450.55b

Notes. The coefficients and uncertainties for ${C}_{1}-{C}_{7}$ from HW99 are listed using the fitting formula, ${a}_{c}/{a}_{\mathrm{bin}}={C}_{1}+{C}_{2}{e}_{\mathrm{bin}}+{C}_{3}{e}_{\mathrm{bin}}^{2}+{C}_{4}\mu +{C}_{5}{e}_{\mathrm{bin}}\mu \,+{C}_{6}{\mu }^{2}+{C}_{7}{e}_{\mathrm{bin}}^{2}{\mu }^{2}$. We perform two separate fits (Fit 1 and Fit 2) using all of our data and list the resulting reduced chi-square value, ${\chi }_{\nu }^{2}$.

^aThis value was calculated using the coefficients listed in HW99 and our larger data set. ^bThis fit modifies the equation where $\mu \to {\mu }^{1/3}$ in order to better match the form with the Hill radius when e_bin and μ are small.

Download table as:  ASCIITypeset image

Although we find good agreement statistically with HW99, the final result will represent CBPs at the population level and there are not enough detections made thus far to justify a completely statistical treatment. Therefore, we suggest a different approach, which is to think of a stability surface (i.e., two-dimensional) rather than a stability limit. In this interpretation, we can obtain much higher accuracy at an individual system level through grid interpolation of our results.

Figure 5 illustrates how our data set can be used to make such a map.⁶ Each point is color-coded to the stability limit, a_c, determined through a smaller grid of simulations (e.g., Figure 7). Additionally, in Figure 5, we have over-plotted (white dots) the locations corresponding to the stellar parameters of the Kepler CBPs. The smallest value of a_c is 1.31 and is located where one would expect. Interestingly, the largest value of a_c is 4.49 and is not produced considering the largest value of μ that we consider. HW99 also observed a similar feature (a_c = 4.2–4.3), but their range in e_bin and resolution did not allow them to identify this location accurately.

Zoom In Zoom Out Reset image size

The Kepler mission has uncovered nine CBP systems, whose stellar and planetary properties vary widely, and comparing them statistically in terms of a stability limit may not be reliable. Li et al. (2016) examined how the mutual inclination of CBPs relative to the binary orbital plane would affect the probability of observing a pile-up of the Kepler CBPs at the stability limit. Their study demonstrated that different conclusions can be drawn when Kepler-1647 is and is not included in the sample of CBPs, and more systems need to be observed in order to distinguish between their two scenarios.

As a result, we compare each of the Kepler CBPs at a system-by-system level using the values of the critical semimajor axis, a_c determined in Section 3.1. We also note that our analysis represents a conservative estimate of the stability limit as our determined limits for a_c could decrease with an increased mutual inclination of the CBP (Doolin & Blundell 2011; Li et al. 2016) or the stellar binary begins closer to apastron rather than periastron (See Section 2.4). Table 3 summarizes the observationally determined stellar masses and orbital parameters of each of the known Kepler CBPs.

Table 3.  Stellar Parameters for the Kepler CBPs

	MA (${M}_{\odot }$)	MB (${M}_{\odot }$)	μ	abin (au)	ebin	ω (deg.)	MA (deg.)	Ref.
Kepler-16	0.6897	0.20255	0.2270	0.22431	0.15944	263.464	188.888	Doyle et al. (2011)
Kepler-34	1.0479	1.0208	0.4934	0.22882	0.52087	71.437	228.760	Welsh et al. (2012)
Kepler-35	0.8877	0.8094	0.4769	0.17617	0.1421	89.1784	2.9021	Welsh et al. (2012)
Kepler-38	0.949	0.249	0.208	0.1469	0.1032	268.68	181.32	Orosz et al. (2012b)
Kepler-47	0.957	0.342	0.263	0.08145	0.0288	226.253	310.818	Orosz et al. (2012a)
Kepler-64	1.528	0.408	0.211	0.1744	0.2117	219.7504	251.558	Schwamb et al. (2013)
Kepler-413	0.820	0.5423	0.398	0.10148	0.0365	279.54	169.5328	Kostov et al. (2014)
Kepler-453	0.944	0.1951	0.171	0.185319	0.0524	263.05	187.7059	Welsh et al. (2015)
Kepler-1647	1.2207	0.9678	0.4422	0.1276	0.1602	300.5442	139.0749	Kostov et al. (2016)

Note. The stellar parameters (M_A, M_B, μ, a_bin, e_bin, ω, and MA) of the Kepler CBPs are listed. The definitions of these orbital parameters carry their usual meaning from the exoplanet literature.

Download table as:  ASCIITypeset image

To measure the proximity of the Kepler CBPs to our determined stability limit, we first determine a_c through a grid interpolation of Figure 5 using the μ and e_bin values given in Table 3. The result of this interpolation for each CBP is given in Table 4. The observed planetary semimajor axis a_p is also provided along with a measure of the percentage difference between a_p and a_c. The comparison using percent difference shows that some CBPs are much closer to a_c than others, where the average difference between a_p and a_c is ∼42%. For systems that are much lower than 42%, we initially classify them to be at the stability limit, and those that are much higher are not at the stability limit.

Another method for determining the proximity to the stability limit uses formalisms from planet packing studies (e.g., Kratter & Shannon 2014) and requires the calculation of the mutual Hill radius, R_H,m (See Section 2.3). For this calculation, we propose that planets classified to reside at the stability limit should not allow for an additional equal-mass planet to exist on an interior orbit at our determined a_c. Along with R_H,m, we also determine the dynamical spacing, β_c, between an equal-mass planet at a_c relative to the observed planet at a_p in Table 4.

Kratter & Shannon (2014) determined that stability is possible with β_c = 5–7, where we define in this analysis that β_c ≤ 7 does not allow for an interior equal-mass planet to exist at a_c. Using this criterion, we find that five out of nine CBP systems (55%) do allow for an interior equal-mass planet. However, the previous study (Kratter & Shannon 2014) did not take the binary eccentricity into account, and we perform a limited suite of N-body simulations to confirm the above estimate for the Kepler CBPs.

In these simulations, we introduce an equal-mass planet with a semimajor axis between a_bin and q_p ( = a_p(1–e_p)) with steps of 0.001 au. The binary can induce a forced eccentricity on the inner planet (e.g., Mudryk & Wu 2006). As a result, we choose the initial eccentricity vectors of both planets to be aligned (ω = 0°) with the binary orbit and vary the magnitude of the eccentricity vector from 0.0 to 0.50 in steps of 0.01. We follow a similar relative phase setup from Gladman (1993), where each planet pair starts 180^° out of phase from one another and the inner planet begins at periastron (MA = 0°). In order to identify robust regions of stability, these simulations are integrated up to 500 million orbits for a planet at a_c (see Table 4 for values of T_c). The integration step is adjusted for each simulation at 2.5% of the initial Keplerian period for the inner planet.

Our full N-body simulations justify our criterion, β_c > 7, to allow for additional equal-mass planets to stably orbit within five of the Kepler CBP systems. Figure 6 illustrates the relative distribution of the Kepler CBPs through their respective values of β_c as concentric circles, where the origin denotes the location that is exactly at the stability limit. In this schematic, the dynamical separation, β_c, from the stability limit does not appear to cluster at any particular value. If we choose the inner planet mass to be Earth-like, then our values of β_c would increase by ∼2^1/3 and potentially allow for an additional planet in Kepler-35. Note: we emphasize that we are not confirming the existence of any planets interior to the observed Kepler CBPs. If such planets do exist, then they must be on sufficiently inclined orbits at the present epoch to have avoided detection.

Zoom In Zoom Out Reset image size

The Kepler CBPs are a snapshot of the larger population of CBPs, where we want to investigate them at a system-by-system level. We identify the ways that the initial conditions could alter our determination of the critical semimajor axis ratio, a_c, by choices in the: initial phase of the binary orbit, initial phase of the planetary orbit, or the initial eccentricity of the planetary orbit.

We examine the stability of systems similar to the Kepler CBPs in the binary parameters (μ, e_bin) using the results from our simulations in Figure 5. Figure 7 illustrates the variation of stability over 10⁵ binary orbits with respect to variations in the initial semimajor axis ratio and mean anomaly of the planetary orbit, when the binary begins at periastron (λ_bin = 0°). We emphasize these assumptions because the actual Kepler CBPs will likely not adhere to them and shifts in the initial phase may be necessary for 1:1 comparisons.

Zoom In Zoom Out Reset image size

In Figure 7, most systems are not strongly dependent on the choice of the planetary mean anomaly (except Kepler-34, Kepler-38, and Kepler-64) and our definition of stability (see Section 2.2) appears to be robust when stability islands exist at specific ranges in the mean anomaly of the planetary orbit. Comparing the values of a_c at these points are consistent (within ∼1%) to those given in Table 4, after multiplying by a_bin, that were determined through a grid interpolation.

We go beyond our ideal setup (excluding Kepler-1647, see Kostov et al. (2016) for a stability map) that makes assumptions on the binary and planetary orbit. For this, we evaluate the variation of stability considering the actual host binary orbit (see Table 3) with a range of initial eccentricity (0–0.5 in steps of 0.01) and semimajor axis (a_bin–1.5 au in steps of 0.001 au) for a Jupiter-mass planet. The planet begins along the reference node so that ω = Ω = MA = 0°. We plot the initial conditions that are stable for at least 100,000 years in Figure 8 using a color-code, the location (green dot) of the observed Kepler CBP parameters, and over-plot the approximate stability boundary for three methods: our Fit 1 (cyan, see Table 2), our interpolation (yellow, see Table 4), and HW99 (violet, see Table 2) using mean values where applicable. The stable initial conditions are color-coded based upon the range of eccentricity (Δe = e_max−e_min) a planet attains over the simulation time on a base-10 logarithmic scale. Ramos et al. (2015) and Giuppone & Correia (2017) have used a similar metric because it highlights dynamical regions affected by resonant interactions, where our definition differs from theirs by a factor of two for clarity. Our results from Section 3.2 investigating whether interior planets could be possible are also shown as gray squares, where those simulations used planet pairs more similar to the actual mass of the Kepler CBPs.

Zoom In Zoom Out Reset image size

For the stability boundary, we assume that a critical pericenter distance, q_c, exists for an eccentric orbit that corresponds approximately to the critical semimajor axis for a circular orbit (i.e., Popova & Shevchenko 2016). The mathematical expression that we use to approximate the critical eccentricity, e_c, is:

$$\begin{eqnarray}&&{e}_{c}=0.8\left(1-\displaystyle \frac{{a}_{c}}{{a}_{p}}\right),\end{eqnarray} \tag{ 3 }$$

where a_c is the critical semimajor axis (in au) derived via each method, and a_p is the planetary semimajor axis in au. The 0.8 factor in our equation is arbitrary, but we found that using this value consistently improves the fit of the upper boundary (high values of e_p) of stability for most cases. Our interpolation method (yellow curve) in Figure 8 typically agrees well with the innermost stable circular orbit (within ∼1%). The other two methods (Fit 1 and HW99) also agree within their error limits, although a substantial fraction within the error range exists in a region of unstable parameter space.

With a sample of only 10 planets, one of which is already an outlier in terms of orbital separation (Kepler-1647b), interpretations of the available data—such as the proposed pile-up of CBPs at the dynamical stability limit—may be affected by observational bias. Increasing the sample by a factor of two would be very useful. Increasing it by a factor of 10 would be fantastic, and with the help of the tools we develop here, would enable comprehensive statistical studies for or against the potential pile-up.

An order of magnitude increase in the number of known CBPs will be indeed possible with the TESS mission by using a novel method for CBP detection based on the occurrence of multiple transits during the same conjunction. This method has already been demonstrated for the case of two such transits of Kepler-1647b, where Kostov et al. (2016) estimated a planet period within 5% of the true period by combining radial velocity measurements with transit timing-independently from the full photodynamical solution of the systems. Similar transits can easily occur within the 30-day all-sky observing window of TESS.

TESS will observe the entire sky for at least 30 days, and continuously measure the brightness of ∼20 million stars (including ∼500,000 eclipsing binaries) brighter than R ∼ 15 with mmag precision (Sullivan et al. 2015). Based on the CBP results from Kepler, we expect the CBP yield of TESS to be a few hundred planets similar to Kepler's (V. B. Kostov et al. 2018, in preparation). The tools and methods we describe here will be directly applicable for both estimating the stability of each new CBP candidate during the initial detection phase, as well as for detailed dynamical investigations of the entire sample after the comprehensive photodynamical characterization of all planets. For example, our method would allow rapid identification (<1 s) of the likelihood that a candidate is a false positive—and thus immediately guide follow-up efforts—based on stability criteria and dynamical packing. In addition, if any TESS CBP system exhibits extra transits, not associated with either the binary or the detected planets, by applying the methodology presented here, we will be able to rule out the orbital parameter space available to additional planets in the respective systems.

We perform a multitude of simulations into order to determine the most general and reliable stability limit given a set of binary parameters (μ, e_bin). The results of these simulations are available through GitHub⁷ and Zenodo.⁸ The GitHub repository contains scripts to identify the stability limit, a_c, and reproduce the figures contained in this paper using Matplotlib (Hunter 2007; Droettboom et al. 2016). The determination of a_c is not limited to the binary parameters used to make Figure 5, but can be interpolated using routines from Scipy (Jones et al. 2001) in Python or other programming languages (Press et al. 1992).

The full data set is available as a compressed tar archive on Zenodo. The archive contains text files that are delineated by the assumed binary parameters (μ, e_bin) in the file names. Python scripts to manipulate the data set without extracting all the files are available in the GitHub repository. Each comma delimited file in the archive lists the results of a given simulation, where the columns are the initial semimajor axis ratio, the initial planetary phase in degrees, the maximum planetary eccentricity attained, the minimum eccentricity attained, and the collision/escape time in years. For initial conditions that survived the full simulation time, a value of 10⁵ years is reported the final column.