LONG-TERM STABILITY OF PLANETS IN THE α CENTAURI SYSTEM

We study particles orbiting α Cen A or α Cen B over a range of semimajor axes (0.2–6.0 au), as motivated from the previous study by WH97 and from the statistical fitting formula developed for planets within binaries (Holman & Wiegert 1999). In contrast to WH97, we use the hybrid symplectic method for integration with mercury6. This allows us to be more uniform with the choice of the starting timestep while maintaining a tight control on the errors in energy and angular momentum. We employ a timestep of 0.005 year step⁻¹ for all our circumstellar runs to properly handle moderately close approaches with the host star.

We analyze two separate regions of parameter space to estimate the stability of test bodies that interact through gravitational forces from the stars. In each region, we evaluate a grid of initial conditions across a range of semimajor axes considering either initially circular orbits with increments in the mutual inclination relative to the binary plane of motion or planar orbits with a range of eccentricities for the test body.

The mutual inclination regime ($a,i$) follows the methodology of WH97, where the test bodies begin on nearly circular (${e}_{o}={10}^{-6}$) circumstellar orbits starting with a semimajor axis from 0.2 to 6.0 au with an increment of 0.01 au relative to the host star. We evaluate 581 test particles per degree of mutual inclination, resulting in a grid of (181 × 581) different initial conditions that provide a high-resolution view of the system. Within this parameter space, two characterizations of motion develop that are described as either prograde or retrograde. In our usage these terms refer to the direction of orbital motion relative to the binary orbit, i.e., a prograde particle orbits its star in the same direction as the hosting star does around the center of mass. Prograde and retrograde classifications are assigned by the starting inclination of the test particle, where starting values of i less than 90° correspond to prograde and retrograde orbits begin with values of i greater than 90°. Although retrograde objects exist within the Solar System, they are expected to have undergone a more complicated early evolution with gas disk interactions (Jewitt & Haghighipour 2007), higher-order three-body effects (Lithwick & Naoz 2011; Naoz et al. 2011), and even possibly have arisen from a more drastic shift in the architecture of the Solar System giant planets (Gomes et al. 2005; Tsiganis et al. 2005). Thus, we may consider the retrograde cases to be less likely, but include them in our analyses for a more complete picture.

The second parameter space ($a,e$) prescribes the starting values of eccentricity and semimajor axis of the test bodies and assumes that they begin with a small (10⁻⁶ degree) inclination relative to the binary to allow for the freedom of motion in the third dimension. The semimajor axis values follow the same range and resolution as before, but the starting eccentricity of each test body ranges from 0.0 to 1.0 with increments of 0.01, producing a grid of (101 × 581) initial conditions. The initial configuration of the binary and phase of the test objects in these simulations also follows the prescription detailed in the previous regime.

In order to efficiently explore the transition to instability for our circumstellar runs around α Cen B on timescales up to 1 Gyr, we initially perform our runs for 10 Myr and based on those results we only continue a subset of surviving particles near the stability boundary for the full 1 Gyr. Approximate symmetry exists between the masses of the stars, so we only consider orbits around α Cen A for 100 Myr.

A different class of planet may be present that orbits both stars (circumbinary) in the system of α Cen AB around the barycenter. We also consider circumbinary test particles within the mutual inclination regime, where we evaluate test bodies from 35 to 100 au with an increment of 0.1 au. This results in a grid of (181 × 651) different initial conditions. For these integrations, we use a timestep of 2 years step⁻¹ which was adequate to control the numerical error. The computational cost of circumbinary orbits is quite low, which allows us to evaluate all of the (surviving) the test particles in that scenario for the full 1 Gyr.

Figure 1 shows the removal timescale as a function of the starting parameters of a test object. Large regions of parameter space become unstable on a relatively short ($\ll 1$ Myr) timescale. The boundary between stable and unstable starting conditions is complex. Small regions of local instability are present due to N:1 mean motion resonant interactions, where N is an integer, with the stellar perturber. As noted by WH97, the Lidov–Kozai (L–K) mechanism (Kozai 1962; Lidov 1962) is an important process in destabilizing orbits substantially inclined to the plane of the stellar orbit. Moreover, our results confirm the general findings of WH97 that the retrograde starting conditions yield a broader region of parameter space with stable orbits than do prograde initial conditions (Henon 1970).

Zoom In Zoom Out Reset image size

Mean motion resonances (MMRs) extend over a range of mutual inclinations and mainly affect particles with semimajor axes larger than 2 au, whereas the L–K mechanism drives the test particles to high values of eccentricity for the mutual inclinations between the critical values, 392–1408 (Innanen et al. 1997). Our results reproduce a broad region of instability between 75° and 105° that was identified previously by WH97 where all test particles are removed on relatively short timescales. For initial parameters in the prograde ($i\lt 90^\circ$) regime of each panel, the extent of stable test particles decreases with increasing staring inclination. For the retrograde ($i\gt 90^\circ$) regime, similar features (MMRs and the L–K mechanism) are evident but present themselves differently, with larger starting semimajor axes maintaining stability despite an increased proximity to the stellar companion. The stability of test particles that orbit α Cen A follow similar trends as those initially orbiting α Cen B, albeit at slightly larger starting semimajor axis.

The results presented in Figure 1 correlate the initial starting parameters to a characteristic removal time from the system (either by ejection or collision), but the final states are dynamically evolved with larger eccentricities. We illustrate these evolved states using color-coded maps of the test particles that survive 10 Myr around α Cen A with colors representing the apastron distance, Q (Figure 2). The apastron distance evolves with time, so we show maps of both the final apastron values, Q_f, attained over the course of a simulation and the maximum values, Q_max. The values of Q_f illustrate the "phase mixing" of the oscillating eccentricities and give an estimate of likely distances that can be achieved at a given epoch.

Zoom In Zoom Out Reset image size

Most of the same spatial structures are present in Figure 2(a) as in Figure 1(a), but the color scale demonstrates the differences from the expected smooth gradient of values and the largest distances from the host star that are permitted. These differences accentuate the two previously mentioned processes that modify the eccentricity of the test particle over a simulation. For example, there are discontinuities at the border of the L–K regime for both prograde and retrograde orbits. Between the critical values of mutual inclination, 392–1408, the gradient of Q_max is nonlinear (Figure 2(b)), indicating where in the parameter space that the L–K mechanism is effective at increasing the eccentricity. In the retrograde regime, long-lived test particles can attain large Q_max values (up to ∼7 au). In contrast, Figure 2(b) illustrates further that values of Q change in time and that nearby initial conditions can exhibit significantly different variations, as shown by the sharp gradients of colors where dynamical processes are active. Figure 2(c) demonstrates test particles that attain large values of Q (≳5 au) in the retrograde regime and can remain stable over billion-year timescales.

We probe deeper into the full extent of variations due to MMRs and the L–K mechanism by plotting in Figure 3 the median value of the longitude of periastron, $\varpi =\omega \pm {\rm{\Omega }}$, of the surviving test particles after 10 Myr. The N:1 MMRs and the L–K mechanism both involve terms related to the periastron of the orbits that contribute to the calculation of an associated resonant argument, $\phi =\varpi -{\varpi }_{\star }$. We inspect this value (relative to neighboring values) to identify the appearance of possible resonances. We do not determine specifically whether these resonances are long-lived, only that libration of the resonant argument can occur over a fraction of the time during the evolution of the test particles. Figures 1(c) and 3 indicate a spatial and dynamical change in stability near the critical inclination for the L–K mechanism. In regards to the long-term stability, we look at these limiting boundaries of possible resonance and we find that the MMRs can be effective on gigayear timescales up to 30:1 (∼2 au) for the prograde cases and up to 20:1 (∼2.5 au) for the retrograde cases. The boundary for the L–K regime becomes distinct in this view and other variations appear in regions outside the regime for the L–K mechanism, possibly indicative of other secular interactions. Moreover, a large number of points outside of the L–K regime (both prograde and retrograde) tend to have a longitude of periastron nearly equal to that of the binary orbit.

Zoom In Zoom Out Reset image size

We simulated test bodies across a different portion of parameter space, where the initial states began nearly planar with the binary (${i}_{o}={10}^{-6}$ degree) but a full range of initial eccentricities were allowed. Figure 4 shows stability diagrams of this region of parameter space, arranged by timescale and the host star in a similar manner as those in Figure 1. This region of parameter space shows a sharp boundary between the stable and unstable runs. Not surprisingly, the maximum stable semimajor axis decreases when the initial eccentricity, e_o, of the test body is increased; we find that the rate of decrease is slow for ${e}_{o}\lt 0.8$ and rapid for higher eccentricity. Above a starting eccentricity of 0.8, an actual planet would be affected by other interactions (i.e., tides, oblateness, General Relativity) at periastron q ($=a(1-e)$) and/or significant perturbations from the stellar companion at apastron Q. Although we have not considered the effects of tides or General Relativity and our choice of a constant timestep compromises the accuracy of our integrations when q is small, we see the effects of large perturbations at apastron destabilizing the orbits of most test particles with eccentricities greater than 0.8.

Zoom In Zoom Out Reset image size

Within this region of parameter space, we probed how the stability boundary varied with respect to the initial eccentricity of the test body by considering the largest stable starting semimajor axis for a given initial eccentricity. We define these points as the "stability boundary" up to the full simulation time. This procedure was used in previous studies of planetary stability in binary systems (Dvorak 1986; Rabl & Dvorak 1988; Holman & Wiegert 1999), but they have focused on the eccentricity of the stars. Using this stability boundary, we found that the variation of the largest stable periastron distance q was approximately linear with respect to the planetary eccentricity (see also Popova & Shevchenko 2012), which allowed us to apply Monte Carlo methods to determine the best fitting parameters for the slope m and y-intercept b of a linear function (i.e., $y={mx}+b$) within a parameter space of eccentricity and the value of q. We fitted a line to the stable starting values of q, using emcee (Foreman-Mackey et al. 2013), allowing us to transform into our input variables of semimajor axis and eccentricity via the following:

$$\begin{eqnarray}&&e={mq}+b.\end{eqnarray} \tag{ 1 }$$

For low-mass bodies orbiting α Cen A, we find $m=-0.356\pm 0.012$ and $b=0.913\pm 0.011$ and are able to produce a fit to the instability boundary through some algebra to obtain

$$\begin{eqnarray}&&e=\displaystyle \frac{(-0.356\pm 0.012)a+0.913\pm 0.011}{1+(0.356\mp 0.012)a}.\end{eqnarray} \tag{ 2 }$$

The same type of fitting can be applied for test bodies orbiting α Cen B, where we find $m=-0.367\pm 0.012$ and $b=0.890\pm 0.012$. In this case the largest stable semimajor axis for an initially circular orbit (${a}_{o}=-b/m$) is slightly smaller than for bodies that orbit α Cen A. Also, we note that we find the largest stable semimajor axis for an initially planar, circular prograde orbit around α Cen A to be ∼2.56 au after 100 Myr of simulation time, which is larger than the value of 2.34 au from WH97 and suggests that the erosion of initial states for this case is not as severe over long integration timescales as previously indicated.

The removal times do not provide a full dynamical picture of the parameter space, thus we produce views in Figure 5 illustrating the variations in both the maximum and final values of Q. Figure 5(a) demonstrates that Q_max extends smoothly from low to high starting semimajor axes apart from the effects of MMRs near the stability limit. The 20:1 and 15:1 MMRs produce clusters of "lucky" particles that likely survive for long timescales due to their initial conditions and orbital phase relative to the epoch of binary periastron. Figure 5(b) shows timescale variations indicative of oscillation in Q that have become out of phase. The white curves show that Q changes little for ${Q}_{o}\lt 2.5\;{\rm{au}}$. Note that some initial conditions can achieve a value of Q up to $\approx 5\;{\rm{au}}$.

Zoom In Zoom Out Reset image size

We quantify the full extent of our simulations by two measures: the rate of survival/loss as a function of simulation time and the largest distance a bound test particle can achieve over the course of a simulation. Figure 6 demonstrates how the number of test particles remaining N_r (Figure 6(a)) and the number lost N_l (Figure 6(b)) change with time (logarithmic scale). These results in Figure 6 are derived from simulations of test particles initially orbiting α Cen B, and similar results are found when considering α Cen A as the host star.

Zoom In Zoom Out Reset image size

In Figure 6(a), the results have been grouped into prograde (black), retrograde (red), and eccentric (green) as previously defined in Figures 1 and 4. The prograde and retrograde runs each start with 52,290 test particles, whereas the eccentric parameter space begins with 58,681. Of these initial states, the eccentric simulations retain more particles than the prograde, which indicates that inclination has a more adverse effect on stability than does eccentricity.

Figure 6(b) illustrates the time at which test bodies are removed for the initial nearly circular orbit runs, which are divided into prograde and retrograde groupings. The corresponding histogram for the eccentric runs looks similar, but the main peak is shifted slightly to shorter times. For each of the ranges of inclination that we examine, the peak in the bodies lost per logarithmic bin in time occurs near 1000 years. The percentage (of the initial population) lost between 10 and 100 Myr is ∼3%, but only about half as many particles are lost in the next decade (100 Myr–1 Gyr). This amounts to about ∼1500 test bodies being lost between 100 Myr and 1 Gyr of simulation time along the transition region of stability in the inclination parameter space, or around 7 test bodies per degree of mutual inclination. But, differences occur when we consider the magnitude of the peak percentage lost with respect to the L–K regime (dashed) and those outside (dotted) the region. The percentage lost beyond 100 Myr differs in the retrograde regime for initial inclinations above 140° when compared to those below 140°, implying that retrograde test particles near the MMRs become unstable on a longer timescale.

Observers would like to know how far from the host star a typical test body can remain stable. Figure 7 illustrates how the apastron distance Q differs between the prograde (black) and retrograde (red) cases in Figure 7(a). In this view, we see that most of the surviving test particles exist within 3 au (prograde) and 5 au (retrograde). Almost all of the stable prograde test particles remain within 4 au. In contrast, the falloff in particle numbers with distance in the retrograde case is more gradual with some stable particles spending time beyond 6 au and a few reaching a distance of ∼7 au. For the eccentric runs (Figure 7(b)), the bulk of surviving particles are within 3 au as is the case for the inclined, prograde runs; however, the eccentric runs extend ∼0.25 au farther in Q_f and ∼0.5 au in Q_max. This becomes important because the extent of the values of Q gives information as to how far away from a host star one should potentially look for circumstellar planets in α Cen. Also, we note that some of the largest values of Q_max represent particles between 2 and 3 au, where MMRs aid in the excitation of eccentricity.

Zoom In Zoom Out Reset image size

Figure 8 shows our results for planets on circumbinary orbits. As with circumstellar orbits, the extent of the stable region for retrograde orbits is larger than that for prograde orbits. However, inclinations around 90° are typically more stable than lower inclinations, in sharp contrast to the trend for circumstellar case (Section 3.1). This is likely due to a slightly different interaction with the quadrupole moment (Ćuk & Gladman 2005), which couples the variation of inclination with the eccentricity of the inner binary. As a result, initial conditions between ∼40°–140° of mutual inclination can periodically switch from prograde to retrograde. Also, the effects of external MMRs, which can cause instabilities through chaos and resonance overlap (Chirikov 1979; Wisdom 1980; Mudryk & Wu 2006), are clearly evident in Figure 8.

Zoom In Zoom Out Reset image size

We have illustrated the regions of phase space where planets may reside within the α Centauri system, but to find such worlds we also consider how their orbits appear to an observer in our Solar System through a projection onto the sky plane. Planets on circumbinary orbits could appear anywhere in the sky in the general vicinity of α Cen, but those on circumstellar orbits would only appear close to their stellar host. Figure 9 shows the sky projected view of the binary orbit in both the astrocentric (dashed) and barycentric (red and blue) coordinate frames. Because this is a projected view, we note that the points of binary apastron and periastron do not fall on the extreme points of the dashed line. Thus, we have included a projection of the major axis for the astrocentric frame centered on α Cen A.

Zoom In Zoom Out Reset image size

We compare the size, location and shape of the distribution of stable particles on inclined prograde orbits about α Cen A at different phases of the binary orbit. The parameters given in Table 2 represent different ellipses determined through the statistical covariance of the test particles for the prograde, inclined simulations. We produce statistical distributions of the test particles by co-adding them across 10 consecutive binary periods at each of 16 different phases. This process of co-adding helps fill out the area on the sky for which planets might occupy without infringing upon structures introduced by the binary interaction. We then compute the covariances that encloses ∼99% of the particles on the sky plane for a given phase (or mean anomaly) of the binary. The results in Table 2 show that the parameters of the distribution do not vary that much across phases of the binary orbit. Thus we give the co-added distributions across 10 binary periods but sampled at a single binary phase, apastron, in Figure 10. Figure 10 shows views centered on α Cen A using our resulting distribution of stable particles after 10 Myr of their evolution in the eccentric (Figure 10(a)), prograde (Figure 10(b)), and retrograde (Figure 10(c)) runs.

Zoom In Zoom Out Reset image size

Table 2.  Distribution of the Test Particle Survivors

Mean Anom.	$\bar{x}$	${\lambda }_{x}$	$\bar{y}$	${\lambda }_{y}$	θ
(degree)	(au)	(au)	(au)	(au)	(degree)
5.43	0.041	1.053	−0.045	0.525	24.58
72.54	0.027	1.048	−0.021	0.525	24.77
107.56	0.034	1.042	−0.041	0.529	24.54
128.07	0.034	1.045	−0.037	0.522	24.72
142.40	0.033	1.047	−0.032	0.525	25.08
153.70	0.035	1.052	−0.036	0.527	24.51
163.36	0.032	1.044	−0.038	0.526	24.83
172.14	0.032	1.043	−0.048	0.526	24.78
180.55	0.032	1.042	−0.037	0.526	24.67
188.98	0.032	1.044	−0.047	0.526	24.54
197.84	0.032	1.046	−0.045	0.528	24.67
207.66	0.032	1.043	−0.036	0.527	24.72
219.26	0.031	1.045	−0.047	0.525	24.53
234.17	0.031	1.043	−0.037	0.526	24.56
255.90	0.037	1.048	−0.048	0.527	24.84
294.04	0.031	1.049	−0.044	0.526	24.73

Note. Statistical properties in the distribution of test particle survivors (prograde, inclined around α Cen A) on the sky co-added over 10 binary orbits at 16 different mean anomalies that are equally spaced in time. The values in $\bar{x},\bar{y}$ represent the centers of each distribution and the values ${\lambda }_{x},{\lambda }_{y}$ are the eigenvalues of the covariance matrix from the respective distributions. The angle θ corresponds to the angle of the largest eigenvector relative to the positive y-axis of the sky coordinates.

Download table as:  ASCIITypeset image

Figure 10(a) shows the distribution of particles along with a green ellipse with a semimajor axis ${\lambda }_{x}$ and semiminor axis ${\lambda }_{y}$ from Table 2 at ∼180°. We find some structure imposed by the binary orbit onto each distribution through an enhancement (region where stability is more favored) of points on the side between the stars when they are at apastron. Finally, the eccentric runs begin planar with the binary (and do not evolve significantly in inclination), so that the resulting distribution of particles looks like a tilted disk when viewed on the sky.

Figures 10(b) and (c) follow the same procedure but use the results from the prograde and retrograde runs, respectively. These results are naturally more inclined than those in Figure 10(a) and hence they extend further from the binary plane. On the sky plane this extension increases in the northwest and southeast directions. In Figure 10(b), we further decompose these results and show those (in black) with initial mutual inclination below the critical value for the L–K mechanism as a distinct population relative to the distribution (gray) within the L–K regime. These two regimes cover very different areas on the sky, but some ambiguity can exist for orbits ∼1 au from the host star because they can occupy the same area on the sky. Figure 10(c) shows the results using the retrograde runs. In this case a much larger portion of the sky can be covered, with the area enclosed by the green ellipse increasing by ∼50%.