STABILITY ANALYSIS OF SINGLE-PLANET SYSTEMS AND THEIR HABITABLE ZONES

The dynamical stability of extra-solar planetary systems can constrain planet formation models, reveal commonalities among planetary systems, and may even be used to infer the existence of unseen companions. Many authors have studied the dynamical stability of our solar system and extra-solar planetary systems (see Wisdom 1982; Laskar 1989; Rasio & Ford 1996; Chambers 1996; Laughlin & Chambers 2001; Goździewski et al. 2001; Ji et al. 2002; Barnes & Quinn 2004; Ford et al. 2005; Jones et al. 2006; Raymond et al. 2009, for example). These investigations have revealed that planetary systems are close to dynamical instability, illuminated the boundaries between stable and unstable configurations, and identified the parameter space that can support additional planets.

From an astrobiological point of view, dynamically stable habitable zones (HZs) for terrestrial-mass planets (0.3 M_⊕ < M_p < 10 M_⊕) are the most interesting. Classically, the HZ is defined as the circumstellar region in which a terrestrial-mass planet with favorable atmospheric conditions can sustain liquid water on its surface (Huang 1959; Hart 1978; Kasting et al. 1993; Selsis et al. 2007, but see also Barnes et al. 2009).

Previous work (Jones et al. 2001, 2006; Menou & Tabachnik 2003; Sándor et al. 2007) investigated the orbital stability of Earth-mass planets in the HZ of systems with a Jupiter-mass companion. In their pioneering work, Jones et al. (2001) estimated the stability of four known planetary systems in the HZ of their host stars. Menou & Tabachnik (2003) considered the dynamical stability of 100 terrestrial-mass planets (modeled as test particles) in the HZs of the then known 85 extra-solar planetary systems. From their simulations, they generated a tabular list of stable HZs for all observed systems. However, that study did not systematically consider eccentricity, is not generalizable to arbitrary planet masses, and relies on numerical experiments to determine stability. A similar study by Jones et al. (2006) also examined the stability of Earth-mass planets in the HZ. Their results indicated that 41% of the systems in their sample had "sustained habitability." Their simulations were also not generalizable and based on a large set of numerical experiments which assumed that the potentially habitable planet was on a circular orbit. Most recently, Sándor et al. (2007) considered systems consisting of a giant planet with a maximum eccentricity of 0.5 and a terrestrial planet (modeled as a test particle initially in circular orbit). They used relative Lyapunov indicators and fast Lyapunov indicators to identify stable zones and generated a stability catalog, which can be applied to systems with mass ratios in the range 10⁻⁴ to 10⁻² between the giant planet and the star. Although this catalog is generalizable to massive planets between a Saturn mass and 10 M_jup, it still assumes that the terrestrial planet is on a circular orbit.

These studies made great strides toward a universal definition of HZ stability. However, several aspects of each study could be improved, such as a systematic assessment of the stability of terrestrial planets on eccentric orbits, a method that eliminates the need for computationally expensive numerical experiments, and a wide coverage of planetary masses. In this investigation, we address each of these points and develop a simple analytic approach that applies to arbitrary configurations of a giant-plus-terrestrial planetary system.

As of March 2010, 376 extra-solar planetary systems have been detected, and the majority (331, ≈88%) are single-planet systems. This opens up the possibility that there may be additional planets not yet detected, in the stable regions of these systems. According to Wright et al. (2007), more than 30% of known single-planet systems show evidence for additional companions. Furthermore, Marcy et al. (2005a) showed that the distribution of observed planets rises steeply toward smaller masses. The analyses of Wright et al. (2007) and Marcy et al. (2005a) suggest that many systems may have low-mass planets.⁴ Therefore, maps of stable regions in known planetary systems can aid observers in their quest to discover more planets in known systems.

We consider two definitions of dynamical stability. (1) Hill stability. A system is Hill stable if the ordering of planets is conserved, even if the outer-most planet escapes to infinity. (2) Lagrange stability. In this kind of stability, every planet's motion is bounded, i.e., no planet escapes from the system and exchanges are forbidden. Hill stability for a two-planet, non-resonant system can be described by an analytical expression (Marchal & Bozis 1982; Gladman 1993), whereas no analytical criteria are available for Lagrange stability so we investigate it through numerical simulations. Previous studies by Barnes & Greenberg (2006, 2007) showed that Hill stability is a reasonable approximation to Lagrange stability in the case of two approximately Jupiter-mass planets. Part of the goal of our present work is to broaden the parameter space considered by Barnes & Greenberg (2006, 2007).

In this investigation, we explore the stability of hypothetical 1 M_⊕ and 10 M_⊕ planets in the HZ and in the presence of giant and super-Earth planets. We consider non-zero initial eccentricities of terrestrial planets and find that a modified version of the Hill stability criterion adequately describes the Lagrange stability boundary. Furthermore, we provide an analytical expression that identifies the Lagrange stability boundary of two-planet, non-resonant systems.

Utilizing these boundaries, we provide a catalog of fractions of HZs that are Lagrange stable for terrestrial-mass planets in all the currently known single-planet systems. This catalog can help guide observers toward systems that can host terrestrial-size planets in their HZ.

The plan of our paper is as follows. In Section 2, we discuss the Hill and Lagrange stability criteria, describe our numerical methods, and present our model of the HZ. In Section 3, we present our results and explain how to identify the Lagrange stability boundary for any system with one ⩾10 M_⊕ planet and one ⩽10 M_⊕ planet. In Section 4, we apply our results to some of the known single-planet systems. Finally, in Section 5, we summarize the investigation, discuss its importance for observational programs, and suggest directions for future research.

According to Marchal & Bozis (1982), a system is Hill stable if the following inequality is satisfied:

$$\begin{eqnarray} &&{-} \frac{2 M}{G^{2} M^{3}_{\star }} c^{2} h > 1 + 3^{4/3} \frac{m_{1} m_{2}}{m^{2/3}_{3} (m_{1} + m_{2})^{4/3}} +\cdots, \end{eqnarray} \tag{ 1 }$$

where M is the total mass of the system, G is the gravitational constant, M_⋆ = m₁m₂ + m₂m₃ + m₃m₁, c is the total angular momentum of the system, h is the total energy, and m₁, m₂, and m₃ are the masses of the planets and the star, respectively. We call the left-hand side of Equation (1) β and the right-hand side β_crit. If β/β_crit > 1, then a system is definitely Hill stable, if not the Hill stability is unknown.

Studies by Barnes & Greenberg (2006, 2007) found that for two Jupiter-mass planets, if β/β_crit ≳ 1 (and no resonances are present), then the system is Lagrange stable. Moreover, Barnes et al. (2008a) found that systems tend to be packed if β/β_crit ≲ 1.5 and not packed when β/β_crit ≳ 2. Barnes & Greenberg (2007) pointed out that the vast majority of two-planet systems are observed with β/β_crit < 1.5 and hence are packed. Recently, Kopparapu et al. (2009) proposed that the HD 47186 planetary system, with β/β_crit = 6.13, the largest value among known, non-controversial systems that have not been affected by tides,⁵ may have at least one additional (terrestrial mass) companion in the HZ between the two known planets.

To determine the dynamically stable regions around single-planet systems, we numerically explore the mass, semimajor axis, and eccentricity space of model systems, which cover the range of observed extra-solar planets. In all the models (listed in Table 1), we assume that the hypothetical additional planet is either 1 M_⊕ or 10 M_⊕ and consider the following massive companions (which we presume are already known to exist): (1) 10 M_jup, (2) 5.6 M_jup, (3) 3 M_jup, (4) 1.77 M_jup, (5) 1 M_jup, (6) 1.86 M_saturn, (7) 1 M_saturn, (8) 56 M_⊕, (9) 30 M_⊕, (10) 17.7 M_⊕, and (11) 10 M_⊕. Most simulations assume that the host star has the same mass, effective temperature (T_eff), and luminosity as the Sun. Orbital elements such as the longitude of periastron ϖ are chosen randomly before the beginning of the simulation (Equation (1) only depends weakly on them). For "known" Saturns and super-Earths, we fix the semimajor axis a at 0.5 AU (and the HZ is exterior) or at 2 AU (the HZ is interior). For super-Jupiter and Jupiter mass, a is fixed either at 0.25 AU or at 4 AU. These choices allow at least part of the HZ to be Lagrange stable. Although we choose configurations that focus on the HZ, the results should apply to all regions in the system.

We explore dynamical stability by performing a large number of N-body simulations, each with a different initial condition. For the known planet, we keep a constant and vary its initial eccentricity, e, from 0 to 0.6 in steps of 0.05. We calculate β/β_crit from Equation (1), by varying the hypothetical planet's semimajor axis and initial eccentricity. In order to find the Lagrange stability boundary, we perform numerical simulations along a particular β/β_crit curve, with Mercury (Chambers 1999), using the hybrid integrator. We integrate each configuration for 10⁷ yr, long enough to identify unstable regions (Barnes & Quinn 2004). The time step was small enough that energy is conserved to better than 1 part in 10⁶. A system is considered Lagrange unstable if the semimajor axis of the terrestrial-mass planet changes by 15% of the initial value or if the two planets come within 3.5 Hill radii of each other.⁶ In total, we ran ∼70, 000 simulations which required ∼35, 000 hr of CPU time.

We use the definition of the "eccentric habitable zone" (EHZ; Barnes et al. 2008b), which is the HZ from Selsis et al. (2007), with 50% cloud cover, but assumes that the orbit-averaged flux determines surface temperature (Williams & Pollard 2002). In other words, the EHZ is the range of orbits for which a planet receives as much flux over an orbit as a planet on a circular orbit in the HZ of Selsis et al. (2007).

3.1. Jupiter-mass Planet with Hypothetical Earth-mass Planet

In Figures 1 and 2, we show representative results of our numerical simulations from the Jupiter-mass planet with hypothetical Earth-mass planet case discussed in Section 2. In all panels of Figures 1 and 2, the blue squares and red triangles represent Lagrange stable and unstable simulations, respectively, the black circle represents the "known" planet and the shaded green region represents the EHZ. For each case, we also plot β/β_crit contours calculated from Equation (1). In any given panel, as a increases, the curves change from all unstable (all red triangles) to all stable (all blue squares), with a transition region in between.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We designate a particular β/β_crit contour as τ_s, beyond which (larger values) a hypothetical terrestrial-mass planet is stable for all values of a and e, for at least 10⁷ yr. We tested that τ_s is the first β/β_crit (close to the known massive planet) that is completely stable (only blue squares). For β/β_crit curves below τ_s, all or some locations along those curves may be unstable; hence, τ_s is a conservative representation of the Lagrange stability boundary. Similarly, we designate τ_u as the largest value of β/β_crit for which all configurations are unstable. Therefore, the range τ_u < β/β_crit < τ_s is a transition region, where the hypothetical planet's orbit changes from unstable (τ_u) to stable (τ_s). Typically, this transition occurs over 10⁻³β/β_crit. Although Figures 1 and 2 only show the curves in this transition region, we performed many more integrations at larger and smaller values of β/β_crit, but exclude them from the plot to improve readability. For all cases, all our simulations with β/β_crit > τ_s are stable, and all with β/β_crit < τ_u are unstable.

These figures show that the Lagrange stability boundary significantly diverges from Hill stability boundary, as the eccentricity of the known Jupiter-mass planet increases. Moreover, τ_s is more or less independent (within 0.1%) of whether the Jupiter-mass planet lies at 0.25 AU or at 4 AU. If an extra-solar planetary system is known to have a Jupiter-mass planet, then one can calculate β/β_crit over a range of a and e, and those regions with β/β_crit > τ_s are stable. We show explicit examples of this methodology in Section 4.

We also consider host star masses of 0.3 M_☉ and perform additional simulations. We do not show our results here, but they indicate that the mass of the star does not effect stability boundaries.

3.2. Lagrange Stability Boundary as a Function of Planetary Mass and Eccentricity

In this section, we consider the broader range of "known" planetary masses discussed in Section 2 and listed in Table 1. Figures 1 and 2 show that as the eccentricity of the "known" planet increases, τ_s and τ_u appear to change monotonically. This trend is apparent in all our simulations and suggests that τ_s and τ_u may be described by an analytic function of the mass and eccentricity of the known planet. Therefore, instead of plotting the results from these models in a–e space, as shown in Figures 1 and 2, we identified these analytical expressions that relate τ_s and τ_u to mass m₁ and eccentricity e₁ of the known massive planet. Although these fits were made for planets near the host star, these fits should apply in all cases, irrespective of its distance from the star. In the following equations, the parameter x = log [m₁], where m₁ is expressed in Earth masses and y = e₁. The stability boundaries for systems with hypothetical 1 M_⊕ and 10 M_⊕ mass companion are

$$\begin{eqnarray} \tau _\mathrm{j} &=& c_{1}+ \frac{c_{2}}{x} + c_{3} y + \frac{c_{4}}{x^{2}} + c_{5} y^{2} + c_{6} \frac{y}{x} + \frac{c_{7}}{x^{3}} \nonumber\\ &&+ c_{8} y^{3} + c_{9} \frac{y^{2}}{x} + c_{10} \frac{y}{x^{2}}, \end{eqnarray} \tag{ 2 }$$

where j = s, u indicates stable or unstable and the coefficients for each case are given in Table 2.

Table 2. Best-fit Properties for Equation (2)

Coefficients	1 M⊕		10 M⊕
	τs	τu	τs	τu
c1	1.0018	1.0098	0.9868	1.0609
c2	−0.0375	−0.0589	0.0024	−0.3547
c3	0.0633	0.04196	0.1438	0.0105
c4	0.1283	0.1078	0.2155	0.6483
c5	−1.0492	−1.0139	-1.7093	−1.2313
c6	−0.2539	−0.1913	-0.2485	−0.0827
c7	−0.0899	−0.0690	-0.1827	−0.4456
c8	−0.0316	−0.0558	0.1196	−0.0279
c9	0.2349	0.1932	1.8752	0.9615
c10	0.2067	0.1577	-0.0289	0.1042
R2	0.996	0.997	0.931	0.977
σ	0.0065	0.0061	0.0257	0.0141
Max. dev.	0.08	0.08	0.15	0.05

Download table as:  ASCIITypeset image

The coefficients in the above expression were obtained by finding a best-fit curve to our model data that maximizes the R² statistic,

$$\begin{eqnarray} R^{2} &=& 1 - \frac{ \sum _{i}^n \big(\tau ^{\rm model}_{i} - \tau ^{\rm fit}_{i}\big)^2}{\sum _{i}^n \big(\tau ^{\rm model}_{i} - \overline{\tau ^{\rm model}}\big)^2}, \end{eqnarray} \tag{ 3 }$$

where τ^model_i is the ithmodel value of τ from numerical simulations, τ^fit_i is the corresponding model value from the curve fit, $\overline{\tau ^{\rm model}}$ is the average of all the τ_model values, and n = 572 is the number of models (including mass, eccentricity, and locations of the massive planet). Values close to 1 indicate a better quality of the fit. In Figure 3, the top panels (a) and (b) show contour maps of τ_s as a function of log [m₁] and e₁ between model data (solid line) and best fit (dashed line). The R² values for 1 M_⊕ companion (Figure 3(a)) and 10 M_⊕ companion (Figure 3(b)) are 0.99 and 0.93, respectively, for τ_s. In both the cases, the model and the fit deviate when the masses of both the planets are near terrestrial mass. Therefore, our analysis is most robust for more unequal mass planets. The residuals between the model and the predicted τ_s values are also shown in Figure 3(c) (1 M_⊕ companion) and Figure 3(d) (10 M_⊕ companion). The standard deviation of these residuals is 0.0065 and 0.0257 for 1 M_⊕ and 10 M_⊕, respectively, though the 1 M_⊕ case has an outlier which does not significantly effect the fit. The maximum deviation is 0.08 for 1 M_⊕ and 0.15 for 10 M_⊕ cases.

Zoom In Zoom Out Reset image size

The expression given in Equation (2) can be used to identify Lagrange stable regions (β/β_crit > τ_s) for terrestrial-mass planets around stars with one known planet with e ⩽ 0.6 and may provide an important tool for the observers to locate these planets.⁷ Once a Lagrange stability boundary is identified, it is straightforward to calculate the range of a and e that is stable for a hypothetical terrestrial-mass planet, using Equation (1). In the next section, we illustrate the applicability of our method for selected observed systems.

The expressions for τ_s given in Section 3.2 can be very useful in calculating the parts of HZs that are stable for all currently known single-planet systems. In order to calculate this fraction, we used orbital parameters from the Exoplanet Data Explorer maintained by the California Planet Survey consortium⁸ and selected all 236 single-planet systems in this database with masses in the range 10 M_jup–10 M_⊕ and e ⩽ 0.6.

Table 3 lists the properties of the example systems that we consider in Sections 4.1–4.4 along with the orbital parameters of the known companions and stellar mass. The procedure to determine the extent of the stable region for a hypothetical 1 M_⊕ and 10 M_⊕ is as follows. (1) Identify the mass (m₁) and eccentricity (e) of the known planet. (2) Determine τ_s from Equation (2) with coefficients from Table 2. (3) Calculate β/β_crit over the range of orbits (a and e) around the known planet using Equation (1). (4) The Lagrange stability boundary is the β/β_crit = τ_s curve.

4.1. Rho CrB

As an illustration of the internal Jupiter + Earth case, we consider the Rho CrB system. Rho CrB is a G0V star with a mass similar to the Sun, but with greater luminosity. Noyes et al. (1997) discovered a Jupiter-mass planet orbiting at a distance of 0.23 AU with low eccentricity (e = 0.04). Since the current inner edge of the circular HZ of this star lies at 0.90 AU, there is a good possibility for terrestrial planets to remain stable within the HZ. Indeed, Jones et al. (2001) found that stable orbits may be prevalent in the present day circular HZ of Rho CrB for Earth-mass planets.

Figure 4(a) shows the EHZ (green shaded) assuming a 50% cloud cover in the a–e space of Rho CrB. The Jupiter-mass planet is the blue filled circle. Corresponding τ_s, values calculated from Equation (2) for 1 M_⊕ companion (0.998, dashed magenta line) and 10 M_⊕ companion (1.009, black solid line) are also shown. These two contours represent the stable boundary beyond which an Earth mass or super-Earth will remain stable for all values of a and e (cf. Figures 1(a) and (b)). The fraction of HZ (FHZ) that is stable for 1 M_⊕ is 72.2% and for 10 M_⊕ is 77.0%. Therefore, the Lagrange stable region is larger for a larger terrestrial planet. We conclude that the HZ of rho CrB can support terrestrial-mass planets, except for very high eccentricity (e > 0.6).

Zoom In Zoom Out Reset image size

These results are in agreement with the conclusion of Jones et al. (2001) and Menou & Tabachnik (2003), who found that a planet with a mass equivalent to the Earth–moon system, when launched with e = 0 within the HZ of Rho CrB, can remain stable for ∼10⁸ yr. They also varied the mass of Rho CrB b up to 8.8 M_jup and still found that the HZ is stable. Our models also considered systems with 3  M_jup, 5  M_jup, and 10  M_jup and our results show that even for these high masses, if the initial eccentricity of the Earth-mass planet is less than 0.3, then it is stable.

To show the detectability of a 10 M_⊕ planet, we have also drawn an RV contour of 1 m s⁻¹ (red curve), which indicates that a 10 M_⊕ planet in the HZ is detectable. A similar contour for an Earth-mass planet is not shown because the precision required is extremely high.

4.2. HD 164922

4.3. GJ 674

4.4. HD 7924

4.5. Fraction of Stable HZ

For astrobiological purposes, the utility of τ_s is multi-fold. Not only it is useful in identifying stable regions within the HZ of a given system, but it can also provide (based on the range of a and e) what FHZ is stable. We have calculated this fraction for all single-planet systems in the Exoplanet Data Explorer as of 2010 March 25. The distribution of fractions of currently known single-planet systems is shown in Figure 5 and tabulated in Table 4. A bimodal distribution can be clearly seen. Nearly 40% of the systems have more than 90% of their HZ stable and 38% of the systems have less than 10% of their HZ stable. The total FHZ that is stable in all known single-planet systems is ∼50%. Note that if we include systems with masses >10 M_jup and also e > 0.6 (which tend to have a ∼ 1 AU; Wright et al. 2009), the distribution will change and there will be relatively fewer HZs that are fully stable.

Zoom In Zoom Out Reset image size

Table 4. Lagrange Stable (τ_s) and Unstable (τ_u) Boundaries, and the Corresponding FHZ Stable for Terrestrial-mass Planets in Known Single-planet Systems

System	m1(Mjup)	a (AU)	e	τs	τu	FHZ	τs	τu	FHZ
				(1 M⊕)	(1 M⊕)	(1 M⊕)	(10 M⊕)	(10 M⊕)	(10 M⊕)
HD 142b	1.3057	1.04292	0.26	0.9347	0.9323	0.000	0.9552	0.9320	0.000
HD 1237	3.3748	0.49467	0.51	0.7407	0.7401	0.213	0.7549	0.7450	0.000
HD 1461	0.0240	0.06352	0.14	0.9920	0.9780	0.976	1.0200	0.9200	0.959
WASP-1	0.9101	0.03957	0.00	1.0022	0.9990	0.990	1.0200	0.9980	0.991
HIP 2247	5.1232	1.33884	0.54	0.7138	0.7111	0.000	0.7490	0.7387	0.000

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

We have empirically determined the Lagrange stability boundary for a planetary system consisting of one terrestrial-mass planet and one massive planet, with initial eccentricities less than 0.6. Our analysis shows that for two-planet systems with one terrestrial-like planet and one more massive planet, Equation (2) defines Lagrange stable configurations and can be used to identify systems with HZs stable for terrestrial-mass planets. Furthermore, in Table 4 we provide a catalog of exoplanets, identifying the FHZ that is Lagrange stable for terrestrial-mass planets. A full version of the table is available in the electronic edition of the journal.⁹

In order to identify stable configurations for a terrestrial planet, one can calculate a stability boundary (denoted as τ_s in Equation (2)) for a given system (depending on the eccentricity and mass of the known planet) and calculate the range of a and e that can support a terrestrial planet, as shown in Section 4. For the transitional region between unstable and stable (τ_u < β/β_crit < τ_s), a numerical integration should be made. Our results are in general agreement with previous studies (Menou & Tabachnik 2003; Jones et al. 2006; Sándor et al. 2007), but crucially our approach does not (usually) require a large suite of N-body integrations to determine stability.

We have only considered two-planet systems, but the possibility that the star hosts more currently undetected planets is real and may change the stability boundaries outlined here. However, the presence of additional companions will likely decrease the size of the stable regions shown in this study. Therefore, those systems that have fully unstable HZs from our analysis will likely continue to have unstable HZs as more companions are detected (assuming that the mass and orbital parameters of the known planet do not change with these additional discoveries). The discovery of an additional planet outside the HZ that destabilizes the HZ is also an important information.

As more extra-solar planets are discovered, the resources required to follow-up grows. Furthermore, as surveys push to lower planet masses, time on large telescopes is required, which is in limited supply. The study of exoplanets seems poised to transition to an era in which systems with the potential to host terrestrial-mass planets in HZs will be the focus of surveys. With limited resources, it will be important to identify systems that can actually support a planet in the HZ. The parameter τ_s can therefore guide observers as they hunt for the grand prize in exoplanet research, an inhabited planet.

Although the current work focuses on terrestrial-mass planets, the same analysis can be applied to arbitrary configurations that cover all possible orbital parameters. Such a study could represent a significant improvement on the work of Barnes & Greenberg (2007). The results presented here show that β/β_crit = 1 is not always the Lagrange stability boundary, as they suggested. An expansion of this research to a wider range of planetary and stellar masses and larger eccentricities could provide an important tool for determining the stability and packing of exoplanetary systems. Moreover, it could reveal an empirical relationship that describes the Lagrange stability boundary for two-planet systems. As new planets are discovered in the future, the stability maps presented here will guide future research on the stability of extra-solar planetary systems.

R.K. gratefully acknowledges the support of National Science Foundation Grants PHY 06-53462 and PHY 05-55615, and NASA Grant NNG05GF71G, awarded to the Pennsylvania State University. R.B. acknowledges funding from NASA Astrobiology Institute's Virtual Planetary Laboratory lead team, supported by NASA under cooperative agreement NNH05ZDA001C. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. The authors acknowledge the Research Computing and Cyberinfrastructure unit (http://rcc.its.psu.edu) of Information Technology Services at the Pennsylvania State University for providing HPC resources and services that contributed to the research results reported in this paper.