Long-term Stability of Tightly Packed Multi-planet Systems in Prograde, Coplanar, Circumstellar Orbits within the α Centauri AB System

B. Quarles; Jack J. Lissauer

doi:10.3847/1538-3881/aaa966

1. Introduction

The α Centauri system contains three stars. The two largest stars in the system, α Cen A and α Cen B, travel about one another on an eccentric orbit with a periapsis of ∼11 au and are broadly similar to the Sun in mass (1.133 M_⊙ and 0.972 M_⊙, Pourbaix & Boffin 2016) and luminosity (1.519 L_⊙ and 0.5 L_⊙, Thévenin et al. 2002). The loosely bound third star, α Centauri C, is much smaller and fainter and has an estimated orbital semimajor axis ∼9000 au from the center of mass (Kervella et al. 2017); α Cen C is usually referred to as Proxima Centauri because it is the nearest star to our solar system. Models of planetary accretion suggest that circumstellar planets could have formed around α Cen A and α Cen B (Quintana et al. 2002; Quintana 2003), provided the collision velocities of late stage planetesimals are not too large (Thébault et al. 2008, 2009; Thébault & Haghighipour 2015). Recent studies have investigated the necessary disk conditions that produce favorable environments for accretion to dominate over fragmentation, leading to formation of systems of multiple planets (Xie & Zhou 2008; Xie et al. 2010; Rafikov & Silsbee 2015a, 2015b; Silsbee & Rafikov 2015).

Public and scientific interest in the α Centauri AB system has increased in the past few years due to the abundance of exoplanetary detections by the Kepler mission and the recent discovery of an Earth-mass planet orbiting Proxima Centauri, Proxima b (Gillon et al. 2016). Additionally, recent observational studies have indicated the possible existence of a close-in planet α Cen B b (Dumusque et al. 2012), but scrutiny of the analysis methods have placed the claimed detection into doubt (Hatzes 2013; Rajpaul et al. 2016). Observations of the binary system continue with the possible observation (single event) of a transiting planet on a somewhat more distant orbit about α Cen B by Demory et al. (2015), and a ground-based radial velocity campaign has ruled out the presence of very massive close-in planets (Endl et al. 2015). Zhao et al. (2018) have shown that existing data places mass detection thresholds for planets in the classically defined habitable zones (HZs) at ≲50 M_⊕ for α Cen A and ≲10 M_⊕ for α Cen B. Moreover, the possibility of a space mission that will directly image any planets that may orbit within the HZs of α Cen A or α Cen B is currently being studied (Belikov & ACESat Team 2015; Belikov et al. 2015). Close conjunctions of the α Cen pair in the future may also provide a microlensing probe for the detection of planets in the system (Kervella et al. 2016). Our study investigates how the orbits of terrestrial planet systems would evolve on timescales commensurate with the estimated age of the system (∼5–7 Gyr, according to Mamajek & Hillenbrand 2008).

Gladman (1993) used Hill stability and numerical simulations to determine the minimum spacing of two Earth-mass planets orbiting a Sun-like star. He found that an interplanetary spacing of semimajor axes of $\beta \gt 2\sqrt{3}$ in units of the mutual Hill radius of the planets (Equation (1)) ensures stability of the system, whereas even slightly more closely spaced planets suffer a close encounter very rapidly. Deck et al. (2013) showed that eccentric, coplanar systems behave in a more chaotic manner (due to the overlap of first-order MMRs) that depends on the weighted eccentricity and the mass ratio of the planets.

Systems containing four or more massive bodies are fundamentally different from the three-body systems studied by Hill (1878a, 1878b, 1878c), Gladman (1993), and Deck et al. (2013) in that they have more degrees of freedom and are not well confined by integrals of motion. Studies of the planetary spacing required for stability of these many-body systems have relied primarily on numerical integrations. Chambers et al. (1996) further investigated the stability of multiple planetary systems considering up to 20 planets orbiting a Sun-like star, where the interplanetary spacing was normalized relative to the mutual Hill radius. More recent works, e.g., Smith & Lissauer 2009; Pu & Wu 2015; Obertas et al. 2017, have further investigated different aspects of planetary packing involving the variations in the total number, masses, and resonant structure of planets.

Planetary packing in binary star systems has not been extensively explored, due in part to the high dimensionality of the parameter space and the relatively few exoplanets discovered in close binary systems like α Cen.⁴ Thus, we are motivated to engage in these studies using the techniques developed for single-star systems. We build upon our previous studies of orbital stability (Quarles & Lissauer 2016) and minimization of free eccentricity (Quarles et al. 2018) of individual planets in the α Centauri AB system and focus on the number of planets that could reside within each star’s HZ.

This paper presents the results of simulations of planetary systems orbiting either star in the α Cen AB system to identify how many one Earth-mass (M_⊕) planets could survive on nested, prograde, coplanar orbits within each of the stars’ HZs. Our methods are outlined in Section 2. The results of our study are presented in Section 3. We discuss the results of our study and compare them to analogous published studies of planet packing around single stars in Section 4. In Section 5, we provide the general conclusions of our work along with the implications for future work.

2. Methodology

The numerical simulations in this work use a custom version of the mercury6 integration package that has been optimized for the study of planetary bodies within binary star systems (Chambers et al. 2002). Using this numerical tool, we investigate the stability of tightly packed planets within the α Centauri AB binary star system on a ten billion year timescale. We ignore the presence of α Cen C even though it resides in a loosely bound orbit because its pericenter distance, ∼5300 au (Kervella et al. 2017), is large compared to the sizes of both the AB binary orbit and the circumstellar HZs. Recent studies by Worth & Sigurdsson (2016) show that planetary systems near the larger stars can still form despite the possible turbulent evolution of the triple star system at earlier times. We also ignore changes in the binary stars due to stellar evolution, because we are largely sampling possible past states, where the dynamical effects due to stellar evolution are negligible.

In our study, we define long-term stability as the lack of collisions between and/or ejections of any of the planets that we numerically simulate up to the estimated age of the stellar binary system, 6 Gyr. A collision occurs when the distance between two bodies during an encounter is less than the sum of their radii. A planet is considered ejected when its distance from the host star exceeds 100 au. Almost all of our simulations are stopped due to collisions. We keep the initial orbit of the binary stars fixed while varying the number and orbital spacing of a set of coplanar Earth-mass (1 M_⊕) planets. For the binary orbit, we adopt the revised orbital parameters from Pourbaix & Boffin (2016) and utilize the initial mean anomaly prescribed in Quarles & Lissauer (2016). The dynamical outcomes do not change appreciably if one uses the parameters from Pourbaix et al. (2002), as can be seen in Figure 1. As all of the planets travel on prograde paths in the orbital plane of the binary stars, we adopt the dynamicist convention where ${i}_{\mathrm{bin}}=0^\circ$ denotes the reference plane and the longitude of pericenter is redefined as the argument of pericenter ( ${\varpi }_{\mathrm{bin}}\to {\omega }_{\mathrm{bin}}$ ).

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Comparison of the lifetimes of three-planet systems orbiting α Cen B using the older (solid curve) and the newer (dashed curve) parameters for the binary system. The interplanetary spacings β are measured in terms of the mutual Hill radius, which is slightly smaller for the newer parameters because the estimated mass of their host star has increased (see Equation (2)), thereby resulting in slightly larger values of β for a given ratio of orbital periods of neighboring 1 M_⊕ planets. The results of individual runs are connected with a solid or dashed line in order to emphasize the displacements caused by the differences in resonance locations for different stellar parameters.
Download figure:
Standard image High-resolution image

We select initial conditions for the terrestrial planetary systems following the method of Smith & Lissauer (2009, 2010), where the planets begin on coplanar orbits around the host star, α Cen A (1.133 M_⊙) or α Cen B (0.972 M_⊙). The initial semimajor axes of these orbits vary based on an input spacing parameter β and follow an iterative scheme (Smith & Lissauer 2009). In this procedure, the amount of mass that lies interior to the jth planet is given by: ${\widetilde{M}}_{j}={M}_{\star }+(j-1){M}_{\oplus }$ . We then use ${\widetilde{M}}_{j}$ in a recurrence relation to assign the semimajor axis, ${a}_{j+1}$ , of the next planet. The value of ${a}_{j+1}$ depends upon the parameter β and the mutual Hill radius for two adjacent planets,

(Chambers et al. 1996), as follows:

This procedure yields systems where the initial period ratios between adjacent planets are roughly equal.

The boundaries of the habitable zone around each star can be defined relative to the empirical limits derived for the solar system or using more conservative limits that place the inner boundary so that it receives an amount of flux due to the radiation from its host star, S, equal to that intercepted by the Earth,⁵ S_⊕. We employ both definitions, as the actual conditions for planetary habitability are not generally well-defined and depend on several assumed parameters, such as atmospheric composition and evolution. The inner edge of the empirical habitable zone (1.78 S_⊕) is represented by ${r}_{\mathrm{emp}}$ and the outer edge (0.32 S_⊕) by ${r}_{\mathrm{out}}$ (Kasting et al. 1993; Kopparapu et al. 2013). Our definition of the conservative habitable zone places the inner boundary at r_⊕ (1 S_⊕) and uses the same outer edge as r_out of the empirical HZ. The periastron distance between the stars is quite large compared to r_out, so including irradiation by the stellar companion would only change the boundaries of the HZs by a small amount (Eggl et al. 2012; Kaltenegger & Haghighipour 2013), and thus we ignore heating by the companion star. The numerical values of the boundaries of the HZs relative to the host stars are given in Table 1.

Table 1. Boundary Limits around Each Star

Star	Luminosity	${r}_{\mathrm{emp}}$	T_i	${T}_{\mathrm{bin}}/{T}_{i}$	${r}_{\oplus }$	${T}_{\oplus }$	${T}_{\mathrm{bin}}/{T}_{\oplus }$	${r}_{\mathrm{out}}$	a_cr
	( ${L}_{\odot }$ )	(au)	(year)		(au)	(year)		(au)	(au)
α Cen A	1.519	0.924	0.834	95.79	1.233	1.286	62.14	2.179	2.72
α Cen B	0.500	0.530	0.391	204.23	0.707	0.603	132.55	1.250	2.60

Note. Boundary limits of the HZs around each star in terms of the distance (in au), orbital period (in year), and the period ratio relative to the binary orbital period, T_bin. The critical outer boundary of stability (in au) is also given based on Quarles et al. (2018).

Download table as: ASCII Typeset image

We begin our simulations with the first (innermost) planet having a semimajor axis, a₁, at the inner edge of the empirical or conservative habitable zone for α Cen A or α Cen B. Using this reference point, we explore systems expanding our iterative scheme relative to the first planet. The spacing in β scales with a₁, and previous studies around single stars usually set a₁ = 1 au out of convenience, but our simulations consider a significant perturbation from the secondary star, so the choice of a₁ will make a difference. The initial mean anomaly of each planet is spaced by $2\pi j\lambda$ radians = 360j λ degrees, where λ is the golden ratio ( $\lambda =(1+\sqrt{5})/2\approx 1.61803...$ ), j denotes the planet number, with j = 1 corresponding to the innermost planet, and the modulo operator is used to keep the values between 0° and 360° (i.e., ${M}_{1}=222\buildrel{\circ}\over{.} 4922$ ) .

Figure 2 illustrates how our orbital architectures are initially defined for β = 10, shown from a top-down view using the conservative (dark green) and empirical (light and dark green) habitable zone limits. Figures 2(a) and (b) show the relative sizes of the empirical HZs resulting from differences in luminosity between the stellar components and the initial locations of five-planet systems. Using the conservative definition of the habitable zone, Figure 2(c) illustrates that these regions are widely separated even at the time of periastron passage of the binary stars.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Orbital schematics of the initial conditions for planets in our simulations that initially have an orbital spacing β = 10, shown from a top-down view. The empirical HZs, 0.32 S_⊕–1.78 S_⊕, of each star are represented in green, with the dark green portions representing the conservative HZs, 0.32 S_⊕–1 S_⊕. All planets orbit in the counterclockwise direction, as do the stars. Panels (a) and (b) illustrate five-planet systems that begin the innermost planet at the inner edge of the host star’s empirical HZ (light green). Panel (c) illustrates similar conditions, but for the conservative habitable zone (dark green) and presents a broader view of the system including both stars’ HZs when the stars are at periapse.
Download figure:
Standard image High-resolution image

All of the simulations of planets orbiting α Cen A use the same timestep of 15.23 days, which samples the orbit of the innermost planet at least 20 times, and the timestep for simulations around α Cen B is chosen to be 7.15 days to account for the shorter orbital period at the inner edge of the fainter star’s empirical HZ. We iterate in the parameter β within the 2–20 range in steps of 0.05 and from 20 to 40 in steps of 0.1. This is adequate to probe the rise and turn-over of stability in most cases. Obertas et al. (2017) recently probed in β for planets around single stars at a higher resolution by drawing more than 10000 samples from a uniform distribution and found that sampling in β at the resolution that we are using is sufficient in the sense that increased resolution does not dramatically affect the results (i.e., by finding more stable cases in regions of overall unstable ones). We expect that a substantial range in β will be long-lived (t > 6 Gyr) for systems with only a few planets, and thus we stop iterating once 10 runs have reached 10 Gyr.

Our results will be dynamically shaped by outer boundary limits, most notably the outer edge of each star’s HZ, ${r}_{\mathrm{out}}$ , and the outer stability limit for a single planet on a circular orbit, a_cr (Wiegert & Holman 1997; Holman & Wiegert 1999; Quarles & Lissauer 2016; Quarles et al. 2018). These outer boundary limits are most easily defined with respect to a given distance in au, and Table 1 lists the values of a_cr in au recently found for each star by Quarles et al. (2018). However, our study warrants describing the outer boundary limits in units of β for each specified number of 1 M_⊕ planets and each choice of the inner planet’s semimajor axis. Ignoring terms of order ${m}_{j}/{M}_{\star }$ , some algebraic manipulation of Equation (9) of Obertas et al. (2017) enables us to express the value⁶ of β for a system of N planets with the inner planet at a₁ and outer planet at a_N as:

Table 2 gives the values of the parameter β for which the outermost planet would begin at either the outer edge of each star’s HZ or the outer stability limit for a single planet on a circular orbit.

Table 2. Boundary Limits in Units of β

Star	a₁	n_pl	${r}_{\mathrm{out}}$ $(\beta )$	a_cr $(\beta )$
α Cen A	${r}_{\mathrm{emp}}$	2	66.90	81.53
		3	34.95	43.59
		4	23.50	29.44
		5	17.68	22.20

	${r}_{\oplus }$	2	45.86	62.23
		3	23.39	32.29
		4	15.65	21.68
		5	11.75	16.31

α Cen B	${r}_{\mathrm{emp}}$	2	63.58	103.95
		3	33.21	59.39
		5	16.79	30.84
		6	13.45	24.79

	${r}_{\oplus }$	2	43.61	89.98
		3	22.24	49.43
		5	11.18	25.36
		6	8.95	20.36

Note. Outer edge, ${r}_{\mathrm{out}}$ , of the HZs and stability limit, a_cr, for initially circular orbits, given in units of the planetary spacing parameter β. This boundary depends on the host star, the initial semimajor axis of the innermost planet, a₁, and the number of planets in the system, n_pl.

Download table as: ASCII Typeset image

2.1. Free and Forced Eccentricity

A planet beginning on a circular orbit in the HZs of α Cen A or B rapidly reaches an eccentricity of ∼0.05 through eccentricity pumping from the binary companion, but planets assigned their forced eccentricities show smaller oscillations in the total eccentricity. The maximum eccentricity obtained in the initially circular runs is close to $2{e}_{F}$ , where e_F(a) is the eccentricity that is forced by the companion star on a planet orbiting with semimajor axis a. Apsidal precession of the pericenter also occurs in the (initially) circular runs, allowing planetary orbits to rotate from a non-crossing to a crossing orientation. These large eccentricities allow planets to strongly perturb one another, in many cases leading to collisions. Figure 3 illustrates how a three-planet system around α Cen B with β = 11.0 can evolve in eccentricity and periastron angle over a timescale of 100,000 years. The oscillation timescale appears to be ∼18,000 years, and the inner planet’s periapse angle circulates with this period. Figure 4 illustrates a similar system to that represented in Figure 3, but with assigned initial eccentricity equal to e_F(a) for each of the planets. In this case, the variations in eccentricity and periapse angles are much smaller.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Evolution of an initially *circular* three-planet system orbiting α Cen B over the first 100,000 year, where the innermost planet begins at r_emp. The top panel shows the variation of the eccentricities of each of the planets, and the middle panel illustrates the changes in the argument of periastron, ω. These orbital elements couple together in the bottom panel, which presents the eccentricity vector taken by each planet during this time interval.
Download figure:
Standard image High-resolution image

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Early evolution of a three-planet system orbiting α Cen B, where the innermost planet begins at r_emp with initial planetary *eccentricities* approximately equal to their forced eccentricity, thereby minimizing free eccentricities. The top panel shows the variation of eccentricity between the planets, and the middle panel illustrates the changes in the argument of periastron, ω. These orbital elements couple together in the bottom panel, which presents the eccentricity vectors that each planet takes during this interval of time.
Download figure:
Standard image High-resolution image

Processes during planet formation tend to drive planetary systems toward a minimum in free eccentricity. This favors circular orbits in single-star systems but can lead to eccentric orbits in both the circumbinary environment (Bromley & Kenyon 2015) and the circumstellar environment when the host star has a more distant stellar companion. These two facts motivate us to explore the difference in the smallest stable value of β both for planets with initially circular orbits and for planets with initial eccentricities chosen so that they begin with minimized free eccentricity. In order to supply a good estimate of e_F to use as the initial eccentricity for this second set of simulations, we build upon the numerical method described in Andrade-Ines et al. (2016) and the results of Quarles et al. (2018) for prograde bodies. We use the results from Quarles et al. (2018) that includes a piecewise-quadratic formula, which adequately provides a forced eccentricity e_F as a function of the starting semimajor axis for the jth planet a_j in our standard “eccentric” simulations. This formula is applied to orbits around each star in the binary system using the form of ${e}_{F}={C}_{1}{a}_{j}^{2}+{C}_{2}{a}_{j}+{C}_{3}$ , whose coefficients are given in Table 3 for each star. Note that the difference between coefficients for the two stellar components are small because the stars have similar masses.

Table 3. Coefficients for Forced Eccentricity

	Star	C₁	C₂	C₃
$a\lt 2\,\mathrm{au}$	α Cen A	−0.007	0.044	−0.002
	α Cen B	−0.009	0.047	−0.003
$a\geqslant 2\,\mathrm{au}$	α Cen A	−0.027	0.123	−0.080
	α Cen B	−0.030	0.130	−0.085

Note. Coefficients for determining the forced eccentricity, e_F, for planets on prograde orbits as a quadratic function in starting semimajor axis around each stellar component (Quarles et al. 2018).

Download table as: ASCII Typeset image

2.2. Two-, Three-, and Five-planet Systems

We perform a series of simulations considering two-, three-, and five-planet systems around either star in the α Cen system. The inner planet in each of these runs begins at a semimajor axis that corresponds to the inner edge of either the empirical HZ ( ${r}_{\mathrm{emp}}$ ) or the conservative HZ ( ${r}_{\oplus }$ ) of its host star, and we only iterate in β until the first ten runs that survive 10 Gyr are reached. The length of our simulations extend beyond our criterion for stability (6 Gyr) to account for extremely late instabilities (i.e., those that can occur in the future), as the solar system can be considered unstable as per our definition (Batygin & Laughlin 2008; Laskar & Gastineau 2009; Batygin et al. 2015) and there remains some uncertainty in the age of the α Cen AB system (Mamajek & Hillenbrand 2008). In addition, our choice in the initial eccentricity of the planets may affect our results, so we explore a similar set of simulations where each of the planets begins with an eccentricity close to the forced eccentricity at its location in the binary star system.

In most cases, our three- and five-planet runs broadly explore up to β = 40 to show the full range of outcomes. Lifetimes generally increase with β and then decrease once the starting semimajor axis of the outermost planet approaches the stability limit for single-planet systems in the α Cen binary star system (Wiegert & Holman 1997; Holman & Wiegert 1999; Quarles & Lissauer 2016), where the specific values of β are given in Table 2.

We plot the initial semimajor axis of the outermost planet, a₃ or a₅, as a function of β in Figure 5(a), to show how the single-planet stability limit affects the exploration in β. The points in Figure 5(a) are delineated by color (red or blue) representing the host star, by symbol (star or triangle) denoting the number of planets in the system, and by the innermost semimajor axis beginning at the inner edge of the empirical (open) or conservative (filled) HZ of each star. The stability limit (red solid line) for five planets (red open stars) in the conservative HZ orbiting α Cen A is crossed at β ≈ 21, whereas the three-planet (red open triangles) cases crosses the stability boundary at β ≈ 41 (See Table 2 for more precise values.). The outer edge of the HZ (dashed) is also indicated, showing that a range of values in β may be allowed through stability without all the planets residing within the respective HZ of the host star. Some of our simulations include planets exterior to the HZ, making our numerical results applicable to any choice for the outer boundary of the HZ.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** The starting semimajor axis of the outermost planet for three (triangles), four (squares), and five (stars) planets orbiting α Cen A or three (triangles), five (stars), and six planets (hexagons) orbiting α Cen B relative to the spacing parameter β, where the starting semimajor axis begins at the inner edge of the empirical (open symbols) or conservative (filled symbols) HZ. The outer edge of the HZ for each stellar host is represented by the dashed (red and blue) horizontal lines. The solid (red and blue) horizontal lines denote the stability limits for single planets on prograde orbits about α Cen A and α Cen B, respectively, using those established in Quarles et al. (2018).
Download figure:
Standard image High-resolution image

2.3. Four- and Six-planet Systems

One of the primary goals of our simulations is to determine how many planets can reside within the HZs of α Cen A and B. Therefore, in addition to our general study of two-, three-, and five-planet systems, we also perform a limited study of four-planet systems with a forced eccentricity around α Cen A and six-planet systems with a forced eccentricity in α Cen B. The motivation for including these specific systems is largely informed by our results in Section 3.1.2, as discussed in the first paragraph of Section 3.2. We have limited our investigation to β = 25, as the peak in lifetime typically occurs at a smaller value of β than this. We show in Figure 5(b) the location of the outer planet as a function of β and also delineate the stability limit and the outer edge of each star’s respective HZ, ${r}_{\mathrm{out}}$ . (See Table 2 for specific values.)

3. Results

The results of our simulations are plotted in Figures 6–9 and 11. We use the common convention of computing the lifetime of the initial configuration (which we plot on a logarithmic scale) as a function of the planetary spacing parameter β. We delineate our results using a color-code, where black dots signify those runs with a₁ set to the inner edge of the empirical habitable zone, ${r}_{\mathrm{emp}}$ , and red dots represent runs with a₁ set to an orbit where the planet intercepts the same flux of radiation as does the Earth, r_⊕, i.e., the inner edge of our conservative HZ. We refer to a configuration as “stable” when its lifetime extends beyond 6 Gyr. To guide the eye, we have indicated the stable runs by highlighting them with upwards vertical ticks, retaining the color-code to signify the location of the innermost planet, and used open symbols to designate those runs that survived the entire 10 Gyr simulated.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Lifetimes of initially circular two (a), three (b), and five (c) planet systems that orbit α Cen A are shown as a function of the interplanetary spacing parameter β. Black (red) points represent those systems that begin at the inner edge of the empirical (conservative) habitable zone. Open points represent systems that survived for the entire 10 Gyr time interval simulated. Simulations where the survival time t exceeds 6 Gyr are highlighted with vertical ticks in the space beyond 10 Gyr. A vertical (gray) line marks the two-planet single-star stability limit, $\beta =2\sqrt{3}$ .
Download figure:
Standard image High-resolution image

**Figure 6.** Lifetimes of initially circular two (a), three (b), and five (c) planet systems that orbit α Cen A are shown as a function of the interplanetary spacing parameter β. Black (red) points represent those systems that begin at the inner edge of the empirical (conservative) habitable zone. Open points represent systems that survived for the entire 10 Gyr time interval simulated. Simulations where the survival time t exceeds 6 Gyr are highlighted with vertical ticks in the space beyond 10 Gyr. A vertical (gray) line marks the two-planet single-star stability limit, $\beta =2\sqrt{3}$ .
Download figure:
Standard image High-resolution image

3.1. Two-, Three-, and Five-planet Systems

These runs explore the lifetime of nested coplanar orbits of two-, three-, or five-Earth-mass planets orbiting either α Cen A or α Cen B. The stability of two-planet systems has been well-studied around single Sun-like stars (Gladman 1993; Chambers et al. 1996), but has not been extensively explored in the case of binary star systems. The presence of a second star changes the character of the problem because the gravitational potential of the many-body system can be substantially different from one in which only three bodies are present, allowing for more degrees of freedom in the system and thereby preventing systems from becoming Hill stable.

The lifetime verses spacing plots share many characteristics with analogous plots for systems of three or more planets around single stars (Smith & Lissauer 2009; Obertas et al. 2017). For closely spaced systems, the logarithm of system lifetime increases roughly in proportion to planetary spacing, with shorter lifetimes near strong mean motion resonances and scatter caused by chaos. However, systems tend to be shorter-lived than comparable systems orbiting single stars. And for configurations in which the outermost planet is strongly perturbed by the binary companion, system lifetime can flatten out with increased spacing and then turn over, with lifetime decreasing as planetary spacing is increased further.

3.1.1. Circular Orbits

We plot the lifetimes of initially circular, coplanar systems around α Cen A and α Cen B in Figures 6 and 7, respectively. Overall, Figures 6 and 7 show several common features. The transition to stability occurs at substantially higher values of β than in the case of single stars. The additional small perturbations from the secondary star cause the location of the transition to stability to depend on the starting semimajor axis of the innermost planet ( ${r}_{\oplus }$ or ${r}_{\mathrm{emp}}$ ). As is the case for multi-planet systems orbiting single stars (Smith & Lissauer 2009; Obertas et al. 2017), interplanetary perturbations due to mean motion resonances between the planets can reduce the lifetimes of systems by up to several orders of magnitude, resulting in stable islands in β between these resonances. For β > 6, simulations with the innermost planet beginning near the inner edge of the empirical HZ ( ${r}_{\mathrm{emp}}$ ) typically survive for longer times than those runs with the innermost planet initially near the inner edge of the conservative HZ ( ${r}_{\oplus }$ ). Also, systems of planets orbiting α Cen A are typically shorter-lived at a given value of β > 5 than systems with the same number of planets that intercept the same radiation flux around α Cen B, because the difference between the luminosities of the stars is much larger than the difference in their masses.

Figure 6(a) illustrates that two-planet systems around α Cen A with initially circular orbits are much more stable if ${a}_{1}={r}_{\mathrm{emp}}$ than if ${a}_{1}={r}_{\oplus }$ . Stable orbits are possible for relatively tightly spaced (β = 9.55) pairs of planets near the inner boundary of the empirical HZ. In contrast, much more widely spaced orbits (β > 21) are necessary for stability by our definition (survival for >6 Gyr) if the innermost planet begins at 1 Earth flux, although several planet pairs that are much more closely spaced are long-lived, with the system at β = 10 surviving for >4.5 Gyr.

The three-planet runs in Figure 6(b) show stable configurations are achieved only for ${r}_{\mathrm{emp}}$ beginning at β ≈ 23, and all of the three-planet r_⊕ runs undergo collisions between the planets within 1 Gyr. Finally, the five-planet systems in Figure 6(c) become unstable on significantly shorter timescales, where the longest lifetimes are less than 100 Myr. This is largely due to outermost planet beginning close to or beyond the outer stability limit even for the smallest values of β for which three-planet systems are stable with the innermost planet located at ${r}_{\mathrm{emp}};$ see Figure 5(a) and Table 2.

The initially circular, coplanar systems around α Cen B whose lifetimes are displayed in Figure 7 show a behavior intermediate between the systems around α Cen A displayed in Figure 6 and planetary systems orbiting the Sun. The minimum separation for two planets to be stable is β ∼ 6.5 (Figure 7). Three planets become stable at larger values (β ∼ 10) for the r_emp cases as shown in Figure 7(b), somewhat larger than the β ∼ 6.4 minimum spacing of comparable systems around single stars found by Smith & Lissauer (2009). For the r_⊕ runs, stable systems occur when β > 18. Only two of the five-planet systems are long-lived in Figure 7(c) and both have the starting semimajor axis of the innermost planet located at r_emp. We note that even though stable configurations exist, only four of the five planets reside interior to the outer boundary of our defined habitable zone.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Lifetime of two (a), three (b), and five (c) initially circular planet systems that orbit α Cen B as a function of the interplanetary spacing parameter β. See caption to Figure 6 for an explanation of the symbols used on this plot.
Download figure:
Standard image High-resolution image

3.1.2. Orbits with Minimal Free Eccentricity

A single planet initially on a circular orbit about an isolated star will remain on the same orbit indefinitely (neglecting stellar evolution and gravitational radiation), but an eccentric binary companion in a system like α Cen will excite a forced eccentricity. Therefore, we simulate systems of two, three, or five planets around each of the stellar components with a forced eccentricity, e_F, that depends on the starting semimajor axis of each planet in a simulation. All of the planets begin with a common argument of periastron that is aligned with the binary orbit, so that ${\rm{\Delta }}\varpi =0^\circ$ .

Planetary systems around α Cen A with minimized initial free eccentricity are stable for much tighter configurations than those begun on circular orbits, as can be seen by comparing Figures 8 and 9 to Figures 6 and 7. The most tightly packed long-lived (>6 Gyr) three-planet runs occur at β ∼ 11 for e_o = e_F, in contrast to β ∼ 23 for e_o = 0. In Figure 8(c), a region exists between 16 and 21 in β that allows for five-planet configurations with ${a}_{1}={r}_{\mathrm{emp}}$ to survive past 1 Gyr, although they do not last up to 6 Gyr to meet our stability criterion. As for systems of planets with e_o = 0, some regions in β are dominated by chaos or mean motion resonances. However, lifetimes of systems with a forced eccentricity are less sensitive to the starting semimajor axis of the innermost planet a₁, as shown by a smaller deviation between runs at a constant β.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Lifetime of two (a), three (b), and five (c) initially eccentric planet systems that orbit α Cen A, with eccentricity chosen to minimize free eccentricity, as a function of the interplanetary spacing parameter β. See caption to Figure 6 for an explanation of the symbols used on this plot.
Download figure:
Standard image High-resolution image

For two planets initially orbiting α Cen B, Figure 9(a) shows the minimum long-lived β is nearly the same value (∼6.7) as compared to Figure 7(a). But the long-lived three-planet runs (vertical ticks), appear at a smaller value in β for minimized free eccentricity (Figure 9(b)) than those for initially circular orbits (Figure 7(b)), with the difference being much larger for systems beginning at r_⊕ than for those with inner planet at r_emp. This is also the case in the five-planet runs, but to a much stronger degree, with a difference of ∼10 in the smallest stable β for ${a}_{1}={r}_{\mathrm{emp}}$ , and five planets can even survive within the conservative HZ.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** Lifetime of two (a), three (b), and five (c) initially eccentric planet systems that orbit α Cen B as a function of the interplanetary spacing parameter β. See caption to Figure 6 for an explanation of the symbols used on this plot.
Download figure:
Standard image High-resolution image

Our algorithm for choosing eccentricities of each planet independently leads to different eccentricities for each planet and thus to non-zero initial relative eccentricities of the planets. One might think that starting all of the planets with the same eccentricity, set by our estimate of the forced eccentricity of the outermost (and most perturbed) planet, might lead to longer system lifetimes. We thus re-ran the simulations of five planets systems with the inner planet orbiting α Cen A at ${r}_{\mathrm{emp}}$ with all planetary eccentricities initially equal to the forced eccentricity of the outermost planet. As shown in Figure 10, on average, this prescription typically resulted in slightly shorter system lifetimes.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Comparison of lifetimes of systems wherein the initial eccentricities of all the planets are set to the forced eccentricity at the semimajor axis of the outermost planet, ${e}_{F}({a}_{5})$ , (blue) to those with initial eccentricity of each planets is set to the forced eccentricity at its own semimajor axis, ${e}_{F}({a}_{j})$ , (black; plotted in panel (a) using the same data used for the black points in Figure 8(c)) for five planets initially orbiting α Cen A. Panel (a) illustrates the similarities between the runs, and panel (b) shows at an expanded vertical scale the pointwise differences between the runs, with black points representing values of β where the system with ${e}_{F}({a}_{j})$ is longer lived and blue points representing values of β where the system with ${e}_{F}({a}_{5})$ survives longer. Horizontal red lines are plotted to guide the eye and differentiate between regions of random scatter and those that could be more significant. A vertical (gray) line marks the two-planet single-star stability limit, $\beta =2\sqrt{3}$ .
Download figure:
Standard image High-resolution image

**Figure 10.** Comparison of lifetimes of systems wherein the initial eccentricities of all the planets are set to the forced eccentricity at the semimajor axis of the outermost planet, ${e}_{F}({a}_{5})$ , (blue) to those with initial eccentricity of each planets is set to the forced eccentricity at its own semimajor axis, ${e}_{F}({a}_{j})$ , (black; plotted in panel (a) using the same data used for the black points in Figure 8(c)) for five planets initially orbiting α Cen A. Panel (a) illustrates the similarities between the runs, and panel (b) shows at an expanded vertical scale the pointwise differences between the runs, with black points representing values of β where the system with ${e}_{F}({a}_{j})$ is longer lived and blue points representing values of β where the system with ${e}_{F}({a}_{5})$ survives longer. Horizontal red lines are plotted to guide the eye and differentiate between regions of random scatter and those that could be more significant. A vertical (gray) line marks the two-planet single-star stability limit, $\beta =2\sqrt{3}$ .
Download figure:
Standard image High-resolution image

3.2. Four- and Six-planet Systems

Our three-planet simulations around α Cen A that minimize free eccentricity show stable values of β with the semimajor axis of the outermost planet occurring well inside the outer limit of the HZ and the stability limit for single planets (Figure 8(b)), suggesting that four-planet systems might be stable within the HZ of α Cen A. Similarly, our five-planet simulations (Figure 9(c)) suggest that six planets might be stable within the HZ of α Cen B. Thus, we also explore the stability of four-planet systems around α Cen A and six-planet systems around α Cen B. We note that the outer planet is located at the single-planet stability limit for a lower value of β when the number of planets is larger (see Figure 5(b) and Table 2).

Figure 11(a) shows that four-planet configurations beginning with ${r}_{\mathrm{emp}}$ can be stable around α Cen A starting near β = 13.6, well inside the outer edge of the HZ. However, the lifetimes of systems with ${a}_{1}={r}_{\oplus }$ are typically lower than those with ${a}_{1}={r}_{\mathrm{emp}}$ for $\beta \gtrsim 9$ . Few simulations systems with ${a}_{1}={r}_{\oplus }$ survive beyond 1 Gyr and none last long enough to be deemed stable by our definition. None of the five-planet runs in Figure 8(c) are stable for either choice of the innermost planet’s semimajor axis.

Figure 11. Refer to the following caption and surrounding text. — **Figure 11.** Lifetime of four (a) and six (b) initially eccentric planetary systems that orbit α Cen A or α Cen B, respectively, as a function of the interplanetary spacing parameter β. See caption to Figure 6 for an explanation of the symbols used on this plot.
Download figure:
Standard image High-resolution image

The inner edge of the HZ of α Cen B is much closer to the host star, possibly allowing for more planets to be stable within it. However, the width of the HZ is smaller than around α Cen A, albeit only slightly smaller when measured in units of β for Earth-mass planets, and the outer edge of α Cen B’s HZ is well interior to the largest stable prograde orbit of a single planet around α Cen B found by Quarles & Lissauer (2016). More planets could stably orbit the star, but not all the planets would necessarily fit within the boundaries of our prescribed HZs. In a similar fashion to Figure 11(a), Figure 11(b) indicates that six-planet systems in the ${r}_{\mathrm{emp}}$ runs are stable starting at β = 13.50, which places the outermost planet just beyond ${r}_{\mathrm{out}}(\beta )\,=\,13.45$ (see Table 2), i.e., the outer planet orbits just exterior to the outer edge of the HZ. None of the six-planet systems that began with ${a}_{1}={r}_{\oplus }$ survive up to our 6 Gyr threshold for stability, and those that lasted the longest had multiple planets orbiting exterior to the HZ.

3.3. Mean Motion Resonances

The build-up of interplanetary perturbations near mean motion resonances can destabilize planetary systems, and putative planets in orbit about α Cen A or α Cen B would also experience kicks from the stellar companion. We focus on the lifetimes of systems around α Cen A as a function of the spacing in period ratios of adjacent planets. We find similar results for our runs around α Cen B, but because our prescription starts the innermost planet at the inner edge of the host star’s HZ, which lies much closer to α Cen B than for α Cen A, strong effects due to the stellar companion only occur at much larger values of β for a given planet multiplicity. Figures 12–14 display local instabilities near mean motion resonances. The contrast is much stronger in systems with minimized free eccentricity; for instance, the instability near the 6:5 resonance appears much sharper in Figures 13(b) and (c) compared to Figures 12(b) and (c). Systems with minimized free eccentricity are typically longer lived for a given value of β and planetary multiplicity, allowing more time for resonant perturbations to build up.

Figure 12. Refer to the following caption and surrounding text. — **Figure 12.** Data shown in Figure 6, but the horizontal axis has been re-scaled to measure the planetary spacing in terms of the ratio of the initial orbital periods of adjacent pairs of planets. The locations of first- and second-order mean motion resonances of adjacent planets are marked along the top axis. See caption to Figure 6 for an explanation of the symbols used.
Download figure:
Standard image High-resolution image

Figure 13. Refer to the following caption and surrounding text. — **Figure 13.** Data shown in Figure 8, but the horizontal axis has been re-scaled to measure the planetary spacing in terms of the ratio of the initial orbital periods of adjacent pairs of planets. The locations of first- and second-order mean motion resonances of adjacent planets are marked along the top axis. See caption to Figure 6 for an explanation of the symbols used.
Download figure:
Standard image High-resolution image

Figures 12(b) and (c) show small dips in stability occurring exterior to the 7:5 resonance between adjacent planets. These dips are due to the high-order N:1 mean motion resonances (N ≥ 15) of the outer planet with the stellar companion. See Quarles & Lissauer (2016) for a discussion of the effects of these N:1 resonances on the stability of single planets in the α Cen system. This effect is strongest when the planet begins beyond ∼2 au and depends sensitively on the starting semimajor axis of the planet. Figure 14 dramatically shows how the other binary component destabilizes planetary systems in which the outermost planet orbits near one of the high-order N:1 mean motion resonances with the stellar orbit.

Figure 14. Refer to the following caption and surrounding text. — **Figure 14.** Data shown in Figure 11, but the horizontal axis has been re-scaled to measure the planetary spacing in terms of the ratio of the initial orbital periods of adjacent pairs of planets. The locations of first-order mean motion resonances of adjacent planets are marked in black along the top axis. The vertical red lines show the locations of the N:1 mean motion resonances of the outermost planet with the binary orbit for the r_⊕ runs, for which they are stronger at a given value of β, with odd values of N corresponding to these resonances shown above the plots. See caption to Figure 6 for an explanation of the symbols used.
Download figure:
Standard image High-resolution image

4. Discussion

4.1. Effect of Free Eccentricity

Starting planets with their forced values of eccentricity significantly improves how tightly planets can be packed on nested prograde orbits in the HZs of α Cen A and α Cen B. The maximum value of forced eccentricity that we impart for a single planet is ∼0.06, which occurs when the outer planet begins near 2 au from the host star. However, the most tightly packed stable systems typically occur where the outermost planet begins with a semimajor axis less than 1.5 au, and for these configurations our method for forced eccentricity assigns nearly the same value ( ${e}_{j}=0.02\mbox{--}0.03$ ) to each of the planets.

Our method chooses the initial eccentricity of each planet, e_j, individually as a function of its semimajor axis, a_j. We find that this methodology usually produces systems with slightly longer lifetimes than setting the initial eccentricity of all of the planets equal to the forced eccentricity of the outermost planet (Figure 10). This is likely stems from starting each planet near its secular equilibrium solution, which keeps the planetary orbits aligned with the binary orbit (Andrade-Ines et al. 2016; Andrade-Ines & Eggl 2017).

4.2. Power Law Fits to Lifetimes of Closely Packed Systems

Previous studies of planet packing in simulations with four or more massive bodies (planets + star) have sought to fit the logarithm of the lifetimes of the systems to a linear function of the form:

in order to evaluate the slope toward stability relative to the spacing parameter β (Chambers et al. 1996; Smith & Lissauer 2009; Obertas et al. 2017). Their results are summarized in Table 4. We apply this type of analysis independently to each of the 28 sets of our runs whose instability times are plotted in Figures 6–9 and 11. We compare our results for differing sets of system parameters to one another and to the single-star studies by Smith & Lissauer (2009) and Obertas et al. (2017), because the types of planetary systems that they simulated are similar to those being evaluated in our study. Chambers et al. (1996) used the same methodology, but because of their more limited computer resources they did not fit over the same range in β and they also used different planetary masses.

Table 4. Summary Table of Fitted Coefficients from Previous Works

References	n_pl	b	${\sigma }_{b}$	c	${\sigma }_{c}$	$f(2\sqrt{3})$
Chambers et al. (1996)	3	1.176	0.051	−1.663	0.274	2.411
Smith & Lissauer (2009)	3	1.496	…	−3.142	…	2.040

Chambers et al. (1996)	5	0.765	0.03	−0.030	0.192	2.620
Smith & Lissauer (2009)	5	1.012	…	−1.686	…	1.820
Obertas et al. (2017)	5	0.951	…	−1.202	…	2.092

Note. Summary of coefficients (b and c) for linear fits to Equation (4) from previous studies around single stars with respect to the initial number of planets (n_pl) in each system. The range in planetary spacing considered for these fits is $2\sqrt{3}\lt \beta \leqslant 8.3$ . We also provide the value of $\mathrm{log}\ t$ when $\beta =2\sqrt{3}$ (i.e., $f(2\sqrt{3})$ ). The ellipses denote when uncertainties are not available in the respective work. Uncertainties are not available for any of the previous estimates of $f(2\sqrt{3})$ because the values of the parameters b and c are correlated.

Download table as: ASCII Typeset image

Table 5 provides a summary of our results in terms of the key parameters that describe the linear trends at small β for each of our sets of simulations. The mass of α Cen B is less than that of α Cen A, but the luminosity of α Cen B is substantially lower, which shrinks the width of the associated HZ and shifts the inner edge of each HZ closer to α Cen B. This results in a much steeper slope and indicates that planets orbiting in the HZ of α Cen B are less perturbed by the long-range gravitational effects of α Cen A. Note that our results should not be used to extrapolate beyond the fitted region because the binary companion becomes increasingly important as β increases, so the systems deviate from the power law dependence found for planetary systems orbiting a single star. A more sophisticated fitting function is necessary to capture system lifetimes at high β values (β ≳ 9).

Table 5. Summary Table of Fitted Coefficients

Star	e_o	a₁	n_pl	b	${\sigma }_{b}$	c	${\sigma }_{c}$	${b}^{{\prime} }$	${\sigma }_{{b}^{{\prime} }}$	${c}^{{\prime} }$	${\sigma }_{{c}^{{\prime} }}$	$f(2\sqrt{3})$
α Cen A	circular	${r}_{\mathrm{emp}}$	2	0.845	0.066	1.549	0.379	0.854	0.065	4.442	0.152	4.476
			3	0.646	0.040	0.754	0.229	0.626	0.038	3.062	0.091	2.992
			5	0.379	0.028	1.560	0.160	0.385	0.029	2.840	0.073	2.873

		${r}_{\oplus }$	2	0.484	0.066	2.822	0.388	0.487	0.064	4.483	0.162	4.499
			3	0.351	0.038	2.121	0.217	0.359	0.031	3.331	0.079	3.337
			5	0.274	0.023	1.938	0.134	0.258	0.020	2.950	0.046	2.887

	eccentric	${r}_{\mathrm{emp}}$	2	0.953	0.072	1.212	0.392	0.952	0.073	4.527	0.162	4.513
			3	0.823	0.040	−0.306	0.221	0.810	0.036	2.607	0.081	2.545
			4	0.723	0.034	−0.079	0.196	0.705	0.034	2.456	0.079	2.426
			5	0.629	0.030	0.244	0.165	0.625	0.027	2.446	0.067	2.423

		${r}_{\oplus }$	2	0.799	0.092	1.599	0.489	0.793	0.091	4.373	0.188	4.367
			3	0.623	0.034	0.872	0.192	0.627	0.036	2.997	0.088	3.030
			4	0.588	0.038	0.712	0.210	0.558	0.032	2.879	0.077	2.749
			5	0.590	0.028	0.496	0.161	0.571	0.028	2.588	0.067	2.540

α Cen B	circular	${r}_{\mathrm{emp}}$	2	1.494	0.144	−0.920	0.694	1.480	0.145	4.294	0.228	4.255
			3	1.017	0.050	−1.345	0.267	0.996	0.049	2.234	0.102	2.178
			5	0.758	0.039	−0.587	0.211	0.742	0.038	2.084	0.087	2.039

		${r}_{\oplus }$	2	1.227	0.148	0.259	0.695	1.249	0.143	4.471	0.209	4.509
			3	0.825	0.045	−0.249	0.244	0.833	0.048	2.581	0.109	2.609
			5	0.616	0.039	0.193	0.215	0.611	0.041	2.325	0.094	2.327

	eccentric	${r}_{\mathrm{emp}}$	2	1.371	0.136	−0.177	0.655	1.417	0.133	4.499	0.210	4.572
			3	0.931	0.044	−0.789	0.236	0.919	0.042	2.474	0.092	2.436
			5	0.744	0.037	−0.575	0.204	0.792	0.038	1.891	0.083	2.002
			6	0.764	0.031	−0.795	0.173	0.737	0.034	1.888	0.079	1.852

		${r}_{\oplus }$	2	1.653	0.151	−1.822	0.695	1.653	0.150	3.915	0.209	3.904
			3	0.823	0.043	−0.343	0.235	0.796	0.045	2.593	0.101	2.508
			5	0.679	0.034	−0.116	0.191	0.707	0.037	2.146	0.089	2.236
			6	0.706	0.032	−0.447	0.179	0.697	0.033	2.063	0.074	1.999

Note. Summary of coefficients (b and c) and uncertainties ( ${\sigma }_{b}$ and ${\sigma }_{c}$ ) for linear fits to $\mathrm{log}\ t=f(\beta )=b\beta +c$ for each of our runs, where these values depend on the host star, initial eccentricity (e_o), innermost semimajor axis (a₁), and initial number of planets (n_pl) in the system. A second set of coefficients ( ${b}^{{\prime} }$ and ${c}^{{\prime} }$ ) and uncertainties ( ${\sigma }_{{b}^{{\prime} }}$ and ${\sigma }_{{c}^{{\prime} }}$ ) are given where a translation ( ${\beta }_{2\sqrt{3}}=\beta -2\sqrt{3}$ ) is applied. We also provide the value of $\mathrm{log}\ t$ when $\beta =2\sqrt{3}$ (i.e., $f(2\sqrt{3})$ ). The range in β considered for these fits are values greater than $2\sqrt{3}$ and include up to the lesser of 8.3 and ${\beta }^{\dagger }$ .

Download table as: ASCII Typeset image

In addition to performing fits analogous to those of previous studies, we also fit our data to a version of Equation (4) shifted to place the origin where the separation between planetary orbits is equal to the critical value for two planets orbiting an isolated star:

where ${\beta }_{2\sqrt{3}}\equiv \beta -2\sqrt{3}$ and the fitted coefficients in the shifted coordinate system are ${b}^{{\prime} }$ and ${c}^{{\prime} }$ . In order to produce Table 5, we used the python module emcee (Foreman-Mackey et al. 2013) to determine the best-fitting parameters and their associated uncertainties σ. For most of our fits, we include the approximately linear regime, $2\sqrt{3}\leqslant \beta \leqslant 8.3$ ( $0\leqslant {\beta }_{2\sqrt{3}}\,\leqslant 8.3-2\sqrt{3}\approx 4.836$ ), in order to compare as fairly as possible with other works. However for the two-planet runs, the smallest value of β that produces a stable configuration, ${\beta }^{\dagger }$ , is typically less than 8.3, and when this is the case we use ${\beta }^{\dagger }$ as the outer limit for the fit instead. Note that in all cases ${b}^{{\prime} }\approx b$ and ${c}^{{\prime} }\approx f(2\sqrt{3})\equiv 2\sqrt{3}b+c$ .

Figure 15(a) displays the results of the fits for the coefficients b and c for the various sets of planetary systems that we modeled herein as well as those studied by Smith & Lissauer (2009) and Obertas et al. (2017). A strong anti-correlation between b and c is apparent, especially if we restrict our attention to systems with a given number of planets. This trend, combined with the similarity between the lifetimes of systems with planetary separations close to $\beta =2\sqrt{3}$ , motivated us to introduce ${b}^{{\prime} }$ and ${c}^{{\prime} }$ . The differences in the slope of the anti-correlation is not much different between three- and five-planet systems, analogous to the limiting effects on stability due to planet multiplicity around single stars (Chambers et al. 1996). On the other hand, planet multiplicity in α Cen is limited by the outer boundary of stability for single-planet systems in binaries (Holman & Wiegert 1999).

Figure 15. Refer to the following caption and surrounding text. — **Figure 15.** Comparison of the slopes and intercepts of the logarithmic fits given in Equations (4) and (5). Symbol shapes show the number of planets in the system, using the code shown above the plots. Planetary systems orbiting α Cen A are shown in red, systems orbiting α Cen B use blue, and the black points are for single-star planetary systems (Smith & Lissauer 2009; Obertas et al. 2017). Open and filled symbols correspond to systems that begin at the inner edge of the empirical (r_emp) or conservative HZ (r_⊕), respectively. Small pink dots are used to denote systems with initially eccentric planetary orbits. The numerical values of the points in panels a and b are given in Tables 4 and 5, respectively. Panel a shows the traditional fitting for the coefficients b and c, where the cyan lines in this panel show the best linear fit to our results for two-, three-, and five-planet systems. Panel b shows fits for the slope (b′) and y-intercept (c′) with β shifted by $2\sqrt{3}$ .
Download figure:
Standard image High-resolution image

Figure 15(b) shows an analogous plot that displays b′ and c′ for the various sets of planetary systems that we modeled herein, where b and $f(2\sqrt{3})$ are used for the fits in previous studies (Smith & Lissauer 2009; Obertas et al. 2017) that did not fit their data for b′ and c′. The range in the value of c′ in Figure 15(b) is much less than the range in the value of c in Figure 15(a), with the two-planet systems all having very similar values of c′. Transforming to c′ demonstrates the similarity in lifetimes of systems with a given number of planets at $\beta =2\sqrt{3}$ (see Figures 6–9), whereas c estimates system lifetimes substantially outside the fitted region (i.e., β = 0). Systems with 3–6 planets typically have smaller values of b′ and c′ than for the two-planet systems, and still have an anti-correlation between the parameters plotted, albeit a much weaker one than is present in Figure 15(a).

Considering two-, three-, and five-planet systems separately, we fit the (anti-)correlation between b and c and show our result as the three cyan lines in Figure 15(a). We use the slopes of these lines to determine the shift in β required to completely remove the correlation between b and c. We refer to this shift, whose value is ${3.56}_{-0.34}^{+0.42}$ for two-planet systems, ${5.10}_{-0.42}^{+0.43}$ for three-planet systems and ${5.24}_{-0.34}^{+0.36}$ for five-planet systems, as ξ. The traditional logarithmic fits of systems with the same number planets orbiting different stars or with different initial eccentricities have similar lifetimes near β = ξ, which implies that the local interplanetary interactions dominate over perturbations from the stellar companion for β ≲ ξ.

Figures 15(a) and (b) both include eight points representing our results for two-, three-, and five-planet systems. These represent a parameter space in which the star being orbited, assumed initial planetary eccentricities, and orbital distance of the inner planet each can take two values. For a given number of planets, single-star systems have the steepest slopes. The initially circular systems orbiting around α Cen A with innermost planet at r_⊕ (filled red points lacking pink dots) have the flattest slopes because they are subjected to the largest eccentricity forcing by the companion star. In general, for the same number of planets, the less perturbed systems (open symbols with pink dots) tend to have a steeper slope than more perturbed systems. With other parameters held fixed, systems with more planets have shallower slopes. For very closely spaced planets (small values of β), perturbations from the companion star are not significant because of shorter system lifetimes and the smaller semimajor axis of the outermost planet. The lifetimes of such closely spaced systems vary more closely with the orbital periods of the inner planets (Table 1) than with the strength of forcing by the stellar companion.

5. Conclusions

We have numerically integrated tightly packed systems containing from 2–6 planets, each being one Earth-mass (1 M ${}_{\oplus }$ ), on prograde orbits in (and in some cases slightly exterior to) the HZs of our stellar neighbors, α Centauri A and B, for up to 10 Gyr. The orbits of the planets are initially uniformly spaced in units of the mutual Hill radius of the planets, ${R}_{{H}_{j,j+1}}$ , the value of which is given by Equation (1). Table 6 summarizes our results and lists the smallest difference in initial orbital semimajor axes of neighboring planets in units of ${R}_{{H}_{j,j+1}}$ that survives longer than 6 Gyr (the nominal age of the α Centauri system), ${\beta }^{\dagger };$ the number of those planets initially within the designated HZ, n_HZ; and the 10th value of orbital separation parameter β for which the system survived for the entire 10 Gyr simulated, ${\beta }^{\ddagger}$ .

Table 6. Summary of Parameters for the Most Tightly Packed Systems

Star	e_o	a₁	n_pl	${\beta }^{\dagger }$	n_HZ	${\beta }^{\ddagger}$
α Cen A	circular	${r}_{\mathrm{emp}}$	2	9.55	2	12.10
			3	22.80	3	29.30
			5	…	…	…

		${r}_{\oplus }$	2	21.30	2	24.50
			3	…	…	…
			5	…	…	…

	eccentric	${r}_{\mathrm{emp}}$	2	7.95	2	9.75
			3	10.75	3	13.15
			4	13.60	4	16.60
			5	…	…	…
		${r}_{\oplus }$	2	7.45	2	8.60
			3	12.90	3	17.15
			4	…	…	…
			5	…	…	…

α Cen B	circular	${r}_{\mathrm{emp}}$	2	6.75	2	8.60
			3	9.95	3	11.90
			5	22.40	4	…

		${r}_{\oplus }$	2	6.30	2	7.35
			3	18.30	3	21.50
			5	…	…	…

	eccentric	${r}_{\mathrm{emp}}$	2	6.65	2	8.25
			3	8.95	3	11.25
			5	12.65	5	14.20
			6	13.50	5	15.80

		${r}_{\oplus }$	2	6.35	2	7.25
			3	9.90	3	12.10
			5	13.90	4	16.20
			6	…	…	…

Note. Summary of most tightly packed stable (lifetime >6 Gyr) value of β ( ${\beta }^{\dagger }$ ), number of planets in the HZ at ${\beta }^{\dagger }$ (n_HZ), and the 10th value of β from our simulations that remained stable until our integrations were terminated at 10 Gyr ( ${\beta }^{\ddagger}$ ). These values depend on the host star, initial eccentricity (e_o), innermost semimajor axis (a₁), and initial number of planets (n_pl) in the system. Boldface values indicate when the initial number of planets in the HZ is less than the initial number of planets in the system. The ellipses represent sets of systems that survived for $\lt 6\,\mathrm{Gyr}$ for all values of β that we simulated ( ${\beta }^{\dagger }$ ) or for which we did not find 10 configurations that survived up to 10 Gyr ( ${\beta }^{\ddagger}$ ).

Download table as: ASCII Typeset image

Perturbations from the companion star substantially reduce the stability of closely spaced multi-planet systems orbiting within the HZs of α Centauri A and B relative to that of comparably spaced planets orbiting single stars. The most closely spaced two-planet system that we found to survive for at least 6 Gyr was in orbit around α Cen B and initially had β = 6.3, which is much greater than the minimum value of $\beta =2\sqrt{3}\approx 3.464$ required for Hill stability around single stars. The companion star induces a “forced” eccentricity on circumstellar planetary orbits within a binary system, and starting planets with orbital eccentricities close to their forced eccentricities generally allowed for more closely packed systems than for planets with initially circular orbits.

We found that three 1 M_⊕ planets can stably orbit around α Cen A within the conservative HZ and four planets can stably orbit within its optimistic (empirical) HZ. For α Cen B, four planets can survive in the conservative HZ, while five can remain in its optimistic HZ. We also found cases where one additional planet could orbit slightly exterior to the HZ of α Cen B, so it is quite possible that with a higher resolution grid in β and/or a slightly different choice for initial planetary eccentricities, five planets might survive in the conservative HZ and six could remain in its optimistic HZ. (See Table 2 for HZ limits in terms of β.) Spacing planets non-uniformly may also allow more planets to orbit stably within these stars’ HZs, as could more exotic configurations such as retrograde trajectories relative to the binary orbit and co-orbital (e.g., Trojan) planets, but we leave studies of these configurations to future work.

Overall, at least nine planets could orbit within HZs of the stars α Centauri A and B. However, more studies are needed to fully understand the possible formation environment for these types of systems and which types of initial architectures would be more typical. Our results demonstrate some of the various dynamical effects possible within the lifetime of the binary system. If all of the planets initially orbit well interior to the single-planet stability limit found by Quarles & Lissauer (2016), then perturbations from the companion star have only a slight destabilizing effect on three- and five-planet systems; however, they are still fundamentally important to two-planet systems. The regime near the outer stability limit exhibits strong destabilizing perturbations of the outermost planet in the system by long-term external forcing due to N:1 resonances with the stars’ mutual orbit. Fortuitously, multi-planet systems can survive in large regions of the HZs of both stars, and these regions should be the focus of efforts to detect many potential worlds orbiting around our closest stellar neighbors.

The authors thank the anonymous referee for a careful and thoughtful review that resulted in improving the clarity and quality of the manuscript. B.Q. gratefully acknowledges support by an appointment to the NASA Postdoctoral Program at the Ames Research Center, administered by Oak Ridge Associated Universities through a contract with NASA. Some of this project was performed at the OU Supercomputing Center for Education & Research (OSCER) at the University of Oklahoma (OU). We thank R. Belikov and N. Kaib for helpful comments on the manuscript.

Long-term Stability of Tightly Packed Multi-planet Systems in Prograde, Coplanar, Circumstellar Orbits within the α Centauri AB System

Author notes

Notes

Article metrics

Permissions

Share this article

Dates

Abstract

1. Introduction

2. Methodology

2.1. Free and Forced Eccentricity

2.2. Two-, Three-, and Five-planet Systems

2.3. Four- and Six-planet Systems

3. Results

3.1. Two-, Three-, and Five-planet Systems

3.1.1. Circular Orbits

3.1.2. Orbits with Minimal Free Eccentricity

3.2. Four- and Six-planet Systems

3.3. Mean Motion Resonances

4. Discussion

4.1. Effect of Free Eccentricity

4.2. Power Law Fits to Lifetimes of Closely Packed Systems

5. Conclusions

Footnotes