Dynamical stability of the Gliese 581 exoplanetary system

Abstract

Using numerical methods, we investigate the dynamical stability of the Gliese 581 exoplanetary system. The system is known to harbour four planets (b–e). The existence of another planet (g) in the liquid water habitable zone of the star is debated after the latest analyses of the radial velocity (RV) measurements. We integrated the four- and five-planet model of Vogt et al. with initial circular orbits. To characterize stability, the maximum eccentricity was used that the planets reached over the time of the integrations and the Lyapunov characteristic indicator and relative Lyapunov indicator to identify chaotic motion. Since circular orbits in the RV fits seem to be a too strong restriction and the true orbits might be elliptic, we investigated the stability of the planets as a function of their eccentricity. The integration of the circular four-planet model shows that it is stable on a longer time-scale for even an inclination i = 5°, i.e. high planetary masses. A fifth planetary body in the four-planet model could have a stable orbit between the two super-Earth-sized planets c and d, and beyond the orbit of planet d, although another planet would likely only be stable on a circular or near-circular orbit in the habitable zone of the star. Gliese 581 g in the five-planet model would have a dynamically stable orbit, even for a wider range of orbital parameters, but its stability is strongly dependent on the eccentricity of planet d. The low-mass planet e, which quickly became unstable in eccentric models, remains stable in the circular four-planet model, but the stable region around its initial semimajor axis and eccentricity is rather small. The stability of the inner planets e and c is dependent on the eccentricity of the Neptune-sized planet b. The outermost planet d is far away from the adjacent planet c to considerably influence its stability; however, the existence of a planet between the two super-Earth planets c and d constrains its eccentricity.

1 INTRODUCTION

Based on high-precision radial velocity (RV) measurements acquired using the HARPS (Bonfils et al. 2005; Udry et al. 2007; Mayor et al. 2009) and HIRES (Vogt et al. 2010) spectrographs, an exoplanetary system was discovered around the red dwarf star Gliese 581 recently. With new RV measurements over time and different analyses of the individual and combined RV data sets, the number of planets revealed in the planetary system has varied between 3 and 6: Gliese 581 e, b, c are certain, the planetary signal of Gliese 581 d became questionable in light of a recent analysis (Baluev 2013), while the existence of f and g is heavily disputed (Anglada-Escude & Dawson 2010; Forveille et al. 2011; Gregory 2011; Tuomi 2011; Makarov, Berghea & Efroimsky 2012; Tadeu dos Santos et al. 2012; Tuomi & Jenkins 2012; Vogt, Butler & Haghighipour 2012; Hatzes 2013).

The first four planets (b, c, d, and e) were all discovered in the HARPS data sets by Bonfils et al. (2005), Udry et al. (2007) and Mayor et al. (2009), respectively. Gliese 581 b is a Neptune-mass planet on a 5.36-d orbit, Gliese 581 c and d are two super-Earth-sized planets with orbital periods of 12.9 and 67 d, and Gliese 581 e is a low-mass planet closest to the star with an orbital period of 3.15 d. The super-Earth planets in the system are located on the two edges of the liquid water habitable zone, and while the greenhouse effect of the atmosphere would make planet c too warm and therefore unable to host liquid water, high concentrations of carbon dioxide or other greenhouse gases would be sufficient to keep planet d from freezing out at the cold edge of the habitable zone (Selsis et al. 2007).

Dynamical stability and evolution of the three-planet system (Gliese 581 b, c and d) of Udry et al. (2007) was investigated by Beust et al. (2008) for different inclinations and was found almost always dynamically stable, even in close to pole-on configurations. The semimajor axes remained extremely stable, the eccentricities underwent only small amplitude variations over the 10⁸ yr integration time. In this three-planet model, both planet c and d had slightly elliptic orbits (e_c = 0.16, e_d = 0.12), and they showed significant chaotic behaviour based on their Ljapunov exponents. Zollinger & Armstrong (2009) investigated the dynamical possibility of a fourth planet. An ∼ 2.5 M_⊕ planet with a semimajor axis between 0.11 and 0.21 au (between planet c and d) and eccentricity below 0.25 would not have compromised the overall stability of the system, although the fourth planet was later found inside the orbit of planet b at 0.028 au. The four-planet model including the then discovered planet e was integrated by Mayor et al. (2009) for 10⁸ yr, with the eccentricity of planet e and b fixed at 0, of planet c and d greater than 0 (e_c = 0.17, e_d = 0.39). This four-planet solution remained only stable for i ≥ 30°, because planet e escaped the system after a few Myr for lower inclinations. Allowing planet e's eccentricity to be 0.1 tightened even more the minimum inclination constraint. The addition of a fifth ∼ 1–10 M_⊕ planet to this four-planet model within 0.3 au quickly introduced dynamical instability to the system, mainly as a consequence of the highly eccentric orbits of the two super-Earth planets c and d (Tóth, Nagy & Dobos 2011).

Combining the 4.3-yr HARPS set with the 11-yr set of HIRES RVs, Vogt et al. (2010) announced the discovery of planet f and g orbiting the star with periods of 433 and 36.6 d. Both planet f and g were indicated in the combined RV data sets using a Keplerian fit with (forced) circular orbits. The quality of the fit could be improved only when the eccentricities of the 67-and 36.6-d planets’ orbits were allowed to float, with these two planets being in secular resonance. Gliese 581 g with a mass of ∼ 3 M_⊕ would be a rocky planet in the middle of the habitable zone of the system and might be habitable for a wider range of atmospheric conditions.

Bayesian reanalysis of both the combined and individual data sets by Gregory (2011) and Tuomi (2011) confirmed only four clear planetary signals (planet e, b, c and d) and higher probability for the existence of planet f (with an orbital period of ∼400 d). Both studies found that the eccentricities for three of the four orbits are consistent with zero, the orbit of planet d is elliptical (e_d ≃ 0.4), similarly to the earlier four-planet solution of Mayor et al. (2009).

Anglada-Escude, López-Morales & Chambers (2010) noted how solutions fitting RV data sets with a single-planet eccentric orbit can hide two planets in circular resonant orbits. This is because there is a degeneracy between the resonant and eccentric solutions, as their Keplerian motion equations are identical up to the first order in the eccentricity. In a subsequent study, Anglada-Escude & Dawson (2010) showed that the first eccentricity harmonic of Gliese 581 d (∼33.5 d) coincides with a yearly alias of the newly reported planet g (∼33.2 d), thus, the high eccentricity of planet d can partially absorb the signal from planet g. Nevertheless, based on statistical tests they concluded that the presence of planet g is well supported by the available RV data. Tadeu dos Santos et al. (2012) similarly concluded that the existence of the 36-d planet g depends on the eccentricity of the 67-d planet d and its detection requires the assumption that all planets are on circular orbits. The signal of planet f was found, but only in the threshold of their confidence level with a period of ∼455 d.

Forveille et al. (2011) released an additional set of HARPS RV measurements and analysed the then total 7-yr data set. Their four-planet Keplerian-fits, with either fixed or freely floating eccentri-cities, revealed no significant residual signals after identifying four, therefore no support for the two additional planets g and f. The additional measurements revised the mass of planet d down to 6 M_⊕, making a rocky composition more likely. Vogt et al. (2012) then reanalysed this data set and warned that allowing the eccentricities to float, and in particular the eccentricity of planet e to rise, leads to instability and therefore to highly unphysical Keplerian models. None of the N-body simulations of the eccentric Keplerian fit of Forveille et al. (2011) remained dynamically stable on longer time-scales due to high eccentricity of planet e. The four-planet all-circular interacting model of Vogt et al. (2012) on the other hand remained dynamically stable for 20 Myr. Furthermore, it offered confirmative support for a fifth planetary signal near 32–33 d, which could be planet g at its 36-d yearly alias period.

The recent analyses of Baluev (2013) and Tuomi & Jenkins (2012) revealed that the RV data contain a significant correlated noise component (red noise), which was not treated by any previous astrostatistical methods, as they assumed that the measurement errors were statistically independent, i.e. the noise was uncorrelated (white). Both studies thoroughly analysed the individual and combined RV sets using a red-noise model. Baluev (2013) found three robust planetary signals (b, c, e), but the significance of planet d's signal dropped with the red-noise model and its existence will therefore require independent confirmation. Tuomi & Jenkins (2012), on the other hand, identified four clear planets (b–d) based on the two data sets using a red-noise model. Interestingly, the HARPS data alone showed a fifth signal with a period of 190 d, but they concluded against its planetary origin, because it was not detected in the combined RV data sets.

Using submillimetre wavebands of the Herschel Space Observatory, Lestrade et al. (2012) have spatially resolved a cold debris disc around Gliese 581, extending radially from 25 ± 12 to 60 au. The debris disc was determined to have an inclination within the range 30° < i < 70°. Assuming that the disc mid-plane and the planetary orbits are coplanar, this range of inclination makes the masses of the planets no more than ∼1.6 times their measured minimum masses by RV measurements. The known planets around the star within 0.3 au cannot dynamically perturb the disc sufficiently over the life-span of the star (2–8 Gyr), which suggests that another undiscovered planet may exist further out, keeping the comet belt replenished.

Since different statistical methods and noise modelling might lead to different orbital solutions in the future, the debate over the number of planets in the Gliese 581 exoplanetary system is far from over yet. Moreover, Gliese 581 is still the target of exoplanet surveys and the addition of new high-precision RV data, additional or already discovered periodic signals may become statistically significant and prove to be of planetary origin. Nevertheless, dynamical stability of current models is important to investigate, especially in the light of the debate over the orbital eccentricities of the planets.

Here, we will present a numerical investigation into the dynamical stability of the Gliese 581 system, with special focus on stability as a function of eccentricity. Apart from how the eccentricities were treated, the orbital parameters and planetary masses did not change considerably after the four-planet solution of Vogt et al. (2012), so we look into the long-term dynamical stability of this RV fit. From the initially circular orbits, we expect that this model is the most stable dynamically; therefore, we search for stable orbits for a fifth planetary body in the inner part of the system, inside the orbit of planet d. (Here, we exclude the region of long-period orbits, such as planet f in the six-planet model of Vogt et al. (2010) would have. Due to the long distance, the inner planets would not influence their stability; conversely, their low masses would not affect considerably the stability of the inner orbits.) We will see that there are stable regions between planet c and d for a fifth planet in the system, and this region corresponds to the liquid water habitable zone of the star. Given this potential, we investigate the long-term stability of the five-planet model of Vogt et al. (2012), which includes a low-mass planet (the unconfirmed planet g) in this zone. The main objective of our study is to set dynamical constraints on the orbits of these RV fits. Using stability maps established from well-known stability indicators, we will determine the limits of the orbital eccentricities of the planets.

2 DYNAMICAL MODELS AND METHODS

We investigated the dynamical stability of the Gliese 581 planetary system by using numerical integrations of the planetary orbits. The model parameters were taken from the four- and five-planet fits of Vogt et al. (2012), these are presented in Table 1 and 3 for reference. For the mass of Gliese 581, we took 0.31 M_⊙. We made the assumption of coplanarity of all orbits, and in case of the circular four-planet model and with the addition of a fifth test planet, we checked the dynamical stability for different orbital inclinations (to the line of sight), otherwise sin i = 1 (an edge-on system) was assumed. We confined our study area between 0.01 and 0.41 au. The stability of the following models was investigated:

• five-body problem consisting of the star and the four planets (b–e) to check the long-term stability of the four-planet system,

• restricted six-body problem consisting of the star, four planets (b–e) and a massless fifth planet to search for stable regions for an additional planetary body,

• six-body problem consisting of the star and five planets to verify the dynamical stability of the five-planet system including planet g as presumed by Vogt et al. (2012).

The numerical integrations were performed using the method of Lie-integration with an adaptive step-size, which is a very fast and precise integration method due to the recurrence of the Lie-terms (Hanslmeier & Dvorak 1984; Pál & Süli 2007). In order to characterize the stability of the models, we used the Lyapunov characteristic indicator (LCI), the relative Lyapunov indicator (RLI) and the maximum eccentricity method (MEM). The LCI, the finite time approximation of the maximal Lyapunov Exponent, is a well-known chaos indicator of a dynamical system, it estimates the exponential divergence rate of infinitesimally close trajectories in the phase space. The RLI is the difference of the Lyapunov indicators of two very close trajectories and provides an effective tool for detection of chaotic behaviour during short integration times (Sándor et al. 2004). The MEM provides information about the evolution of the orbit during the integration time through the largest value of the eccentricity of the planets and indicates close encounters and escapes from the region of motion as well (e.g. Dvorak et al. 2003; Süli, Dvorak & Freistetter 2005; Nagy, Süli & Érdi 2006). While the chaos indicators LCI and RLI characterize the structure of the phase space, dynamical stability is explicitly described by the maximum eccentricity. Since chaos does not necessarily mean dynamical instability, rather a sensitivity to initial conditions, we perform long-term integrations of the four- and five-planet models as well.

3 STABILITY OF THE four-PLANET SYSTEM

3.1 Basic stability of the four-planet model

First, we investigated the stability of the circular four-planet model of Vogt et al. (2012) in detail, which was found to be stable for at least 100 kyr by the authors. The analysis of RV variations is able to constrain the mass (M) of exoplanets by a lower limit, since only the quantity M sin i is determined, where i is the unknown orbital inclination. Below a given value of i – or, equivalently above given masses for the planets – we expect the system to become unstable as the dynamical interactions can increase with lower inclinations, i.e. higher planetary masses. Although the debris disc around the star narrows down the inclination to 30° < i < 70° possibly for the planets as well (Lestrade et al. 2012), this range is not necessarily where the planetary orbits lie. Therefore, to check for which range of inclination the system has this basic stability, we integrated the circular four-planet model (Table 1) over 1 Myr, varying the inclination by ▵i = 5° from i = 90° (edge on) to i = 5° (almost pole on).

In all cases, the variations of the semimajor axes and eccentricities of the four planets are not significant. The initially circular four-planet system remains stable for 1 Myr down to i = 5° (almost pole-on) and coplanar orbits. Dynamical stability with i ≥ 5° suggests that the mass of each planet can be ∼12.5 times of its minimum mass. For Gliese 581 e, b, c, and d those limits are 21, 183, 62, and 60.2 M_⊕.

3.2 Stability of planet e

The low-mass planet e has little effect on the stability of the more massive planets in the system, but an increase in its eccentricity quickly leads to orbital crossings with planet b or ejection from the system (as it was already pointed out by Mayor et al. 2009; Vogt et al. 2012). To characterize planet e's orbit in the circular four-planet model, we carried out integrations with the following grid of initial semimajor axis and eccentricity: 0.01 ≤ a ≤ 0.041 au with steps ▵a = 0.003 au and 0 ≤ e ≤ 1 with steps ▵e = 0.01, for a time period of 5000P_e, where P_e is the orbital period of planet e.

The maximum eccentricity that planet e reached over the time of the integration is plotted on the a-e stability map in Fig. 1. Planet e's orbit is stable (small e_max), although the vicinity around its initial a_e = 0.028, e_e = 0 values represents a stable island on the map, so relatively small perturbations could lead to instability and thus collision or ejection. Having varied only the eccentricity e_e, the LCI became suddenly high above e_e = 0.22 (data not shown). This indicates chaos, and also points to the highest eccentricity for planet e's orbit to remain stable.

We also performed integrations with different starting positions (l) for planet e varying l by 2° between l = 0°and 360°. The LCI values (data not shown) remain low during the time of the integration for each case, therefore planet e's stability does not depend on its orbital position in the model.

4 STABILITY OF THE five-PLANET SYSTEM

4.1 Stable regions for a fifth planet

We integrated the four-planet system with an additional massless body in order to find regions of possible orbital stability for a fifth planet. The test planet's initial parameters were varied, in one set the semimajor axis within the range 0.01 ≤ a ≤ 0.41 au with steps of ▵a = 0.001 au and the eccentricity 0 ≤ e ≤ 1 with steps of ▵e = 0.01, and in another set for the same semimajor axis range and e = 0, the starting position 0 ≤ l ≤ 360° with steps of ▵l = 2°. The time of the integrations was 10 000 P_test for calculating the LCI, where P_test is the orbital period of the test planet. The LCI values of each integration for the grid of different a, e are plotted in Fig. 2(a), for a, l in Fig. 3. The e_max that the test planet reached over 5000 P_test was calculated for the same grid of a, e, and for eight different starting orbital positions with steps ▵l = 45° between l = 0° and 315°. The average of these eight e_max values (Nagy et al. 2006) are represented on one map in Fig. 2(b).

The orbits of a fifth fictitious planet are chaotic in the vicinity of the innermost planets e and b, but extensive stable regions exist between the larger planets c and d, and outside the orbit of planet d, with relatively low eccentricities (Fig. 2a). The starting position of the fictitious planet only plays a role in its stability around the innermost planets and its orbit is less chaotic in the outer part of the system (Fig. 3). The a-e stability map with the maximum eccentricity essentially shows the same boundaries for stable regions, in fact the average e_max values confirm that irrespective of the starting position of the test planet, its orbit can be stable in the region between the super-Earth planets and beyond them (Fig. 2b). However, unstable configurations can be observed in both of these regions on the map (see further details about the habitable zone in Section 5.2). This is in agreement with the results of Tadeu dos Santos et al. (2012), who gave dynamical constraints for the existence of an additional planet between planet c and d.

The integration of the four-planet model including a fifth test planet was also carried out for different inclinations to the line of sight and for different degrees of the argument of pericentre (ω). The integrations ran for 5000 P_test for inclinations i = 30°, 50°, 70° – in the range where the inclination of the debris disc around the star lies – while the argument of pericentre was varied within the range 0 ≤ w ≤ 360° with steps of ▵ω = 90°. Neither of these parameters changed significantly the possible stable regions in the system for a fifth planet based on the maximum eccentricity of the test planet, in fact the resulting a-e stability maps (see Fig. S1 in the Supplementary Material online) resemble very closely that in Fig. 2(b) in each case.

On the stability maps in Figs 2 and 3, there are small stable islands visible sporadically, where both the e_max and the LCI have low values at the end of the integration. These point to mean-motion resonances (MMRs) between the test planet and the other planets. While no MMR could be identified in the inner part of the system (due to the even parametrization of the stability maps), several MMRs were found involving the test planet and the two super-Earth planets around the habitable zone. The main resonances are indicated in Fig. 2(b) and explained in Table 2. MMRs enhance stability and create stability islands in these region, because planet c and d are assumed to have zero or very low eccentricity orbits (Yoshikawa 1989; Beust & Morbidelli 1996). Nevertheless, none of these low-order MMRs overlap, which could lead to instability in the planetary system.

4.2 Stability of the five-planet model including planet g

Even though the existence of planet g, the fifth planet candidate in the star's classical liquid water habitable zone, is debated and might as well be a false periodic signal caused by noise, we investigated its possibility from a dynamical point of view. Therefore, the five-planet model of Vogt et al. (2012, as a reference see Table 3) containing planet g was integrated for 1 Myr to verify its long-term stability. All five planets remained stable during the long integration time and their orbital parameters did not change significantly, so this five-planet solution of the RV data with the 32-d planet g can be considered dynamically stable.

5 STABILITY AS A FUNCTION OF ECCENTRICITY

5.1 Planets e, b, c

The orbits of the planets in the Gliese 581 planetary system might be elliptic as a consequence of dynamical interactions. Highly eccentric orbits can cause instability though, we derived therefore limits for the eccentricities using a-e stability maps. First, we investigated the situation of the inner planets e, b, and c. We have assumed that given its higher mass, planet b is likely to be more stable in comparison to the adjacent planets, so we looked at the stability of planet e for a range of initial e_e and e_b, and the stability of planet c for a range of initial e_b and e_c (in this case planet e was left out of the system). The maximum eccentricity was calculated over a time period of 5000P_e and 5000P_c, respectively. The eccentricities ranged until the elliptical orbit of one planet crossed the other one's circular orbit.

The maximum eccentricity, plotted in Fig. 4 as a function of initial e_e and e_b, indicates that planet e's eccentricity is coupled to the more massive planet b's eccentricity. If e_b ≃ 0, the limit of planet e's eccentricity is e_e ≤ 0.18. If planet b's eccentricity is above zero but relatively low, then planet e remains stable below e_e ≤ 0.2, which is approximately the same upper limit as the one based on the LCI (in Section 3.2). Small stable regions exist with eccentricities above e_e = 0.3 too (Fig. 4), though these isolated configurations are less likely to remain stable for a longer time.

Looking at the stability of planet c as a function of e_b and e_c in Fig. 5, it is apparent that planet c puts less constraint on the eccentricity of planet b (which is then mainly given by the smaller planet e). While in case of circular or near-circular orbits for planet b, the upper limit of planet c's eccentricity is e_c ≤ 0.32.

5.2 Planets c, d and a planet in the habitable zone

The stability of the outer planet d is less influenced by planet c because of the large distance between them, but a planet between planet c and d in the habitable zone would put constraints on both planet c and d's eccentricity. For this reason, the stability of planet c and d is investigated with the 32-d planet g in the system or in general with another planet in the habitable zone.

We ran two set of integrations to characterize the stability of planet g's orbit as a function of the adjacent planets’ eccentri-cities. In one set, for a time of 5000P_g, the eccentricity of the adjacent planets c and d were varied until the point when their orbits cross each other. The LCI of planet g in Fig. 6 shows that planet g's orbit can be considered stable in the rectangular region on the map where the eccentricity of planet c stays below e_c ≤ 0.32 and of planet d below e_d ≤ 0.28. These eccentricity limits mean r_per = 0.049 au, r_ap = 0.096 au pericentre and apocentre distances for planet c, and r_per = 0.158 au, r_ap = 0.282 au for planet d. The apocentre of the orbit of planet c leaves a little bit more space for planet g (at 0.13 au) than the pericentre of the orbit of planet d. Interestingly, however, planet c then causes instability in the orbit of planet g even in case of moderately close encounters, even though planet c's mass is lower than that of planet d. The reason behind this may be that planet c's orbital period is about the fifth of the period of planet d, thus planet c simply approaches planet g more frequently, causing more perturbation. Referring back to the previous section, planet c's eccentricity as high as e_c ≃ 0.32 would not actually be the only cause of the instability in planet g's orbit. When e_c ≃ 0.32, the adjacent, more massive planet b already considerably affects the stability of planet c (Fig. 5) and this is what may ultimately lead to the instability of planet g.

Since planet d's real eccentricity could be greater than zero based on the RV fits, as well as the true orbital period and thus semimajor axis of planet g can be different, in another set of integrations, also for a time of 5000P_g, we varied a_g in the range between the neighbouring planets, and e_d until planet d's orbit would cross the orbit of planet c. Not surprisingly, as Fig. 7 shows, the more elliptic the orbit of planet d, the further away planet g has to be from it in order to remain stable. At the largest distance, planet g's orbit would be stable in case of e_d ≈ 0.55, but such an elliptic orbit would leave only a very narrow region in the system for the fifth planet. Eccentricity as high as e_d ≃ 0.4 that was suggested in previous models for planet d (Mayor et al. 2009; Forveille et al. 2011; Gregory 2011; Tuomi 2011), would limit planet g's orbit to be in a narrow region between ∼ 0.08 and 0.12 au (Fig. 7). As Tadeu dos Santos et al. (2012) has already pointed it out, stability of planet g is therefore strongly dependent on the eccentricity of the adjacent planet d. More generally, based on the map in Fig. 7, it can be stated that this dependence is true for the stability of a low-mass planet orbiting in the habitable zone between planet c and d. Also, several two- or three-body MMRs between planet g and planet c or d were determined in the region (Fig. 2b), which are indicated here by the small vertical chaotic regions observed on the map in Fig. 7.

Finally, we investigated the general case of a planet in the habitable zone by checking for which range of eccentricity a planet would remain stable here. We integrated the all-circular four-planet model with the addition of a massless test planet and varied its eccentricity. The integration ran for a time of 5000 P_test and for eight different starting orbital positions with steps ▵l = 45° between l = 0° and 315°. This time we calculated the RLI of the test planet as a function of the semimajor axis and eccentricity. The average of these eight RLI values is plotted on the a-e stability map in Fig. 8. In general, the dynamical stability of another planet between planets c and d greatly depends on the eccentricity, and a planet would be only stable in this region on circular or nearly circular orbits. Even this very small stable region close to zero eccentricity is not continuous, probably as a consequence of the above-mentioned MMRs between the fifth planet and planets c and d.

6 CONCLUSIONS

In the process of fitting multiplanet models to precise RV data, it is a question of choice whether to allow the orbital eccentricities to float freely or to hold them fixed at zero. Comparing the four-planet models of Mayor et al. (2009), Forveille et al. (2011) and Vogt et al. (2012) in regards to their long-term dynamical stability, eccentricity plays a great role (see also the discussion of Vogt et al. 2012). When eccentricities are fixed at zero, the degree of freedom in fitting the RV curve is reduced, and such a model is more likely to describe a stable system. As the case of the Gliese 581 planetary system shows, freely floating eccentricities can lead to solutions that are unstable and therefore unphysical despite providing good fits to the RV data. The large uncertainty of eccentricities should warn us that there is not enough RV data yet to be able to constrain orbital parameters in a Keplerian model. Forced circular orbits in the RV fits might be too strong a restriction though.

Long-term numerical integrations show that the all-circular four-planet model of Vogt et al. (2012) remains dynamically stable for 1 Myr, even for an inclination i = 5°, i.e. with ∼12.5 times higher planetary masses than the minimum mass. In this model, the innermost planet e's orbit remains stable during the integration time; however, the stable region on the a-e plane around it is rather small.

A fifth planetary body in the four-planet system could have a stable orbit between the two super-Earth planets c and d, and there is a stable region for an additional planet beyond the orbit of planet d. The five-planet model of Vogt et al. (2012), which includes the disputed planet g, is also dynamically stable on a longer time-scale. The existence of the low-mass planet g in the habitable zone is supported from the dynamical point of view and it can be considered stable for a wider range of orbital configurations.

While the models of Vogt et al. (2012) with initial circular orbits are dynamically stable over a longer time period, the true orbits might be elliptic as a consequence of dynamical interactions between the planets. Our simulations put dynamical constraints on the eccentricity of the planets, and these limits can in turn be used to help find the best stable model to fit the RV measurements. The Neptune-sized planet, Gliese 581 b given its high mass greatly influences the stability of the adjacent planets e and c, and the eccentricity of these planets is therefore dependent on planet b's eccentricity. In case of e_b ≃ 0, dynamical stability of the inner planets is only ensured, if planet e's eccentricity stays below e_e ≤ 0.18, and planet c's eccentricity below e_c ≤ 0.32. In case planet b's orbit is slightly elliptic, the orbit of planet e and c has to be elliptic as well for the inner planets to remain stable. The outermost Gliese 581 d is far away from the neighbouring planet c to considerably influence its stability; however, the existence of a planet between the two super-Earth planets puts constraint on the eccentricity of planet d. If a low-mass planet like the 32-d planet g orbits in the habitable zone, then planet d's eccentricity cannot be larger than ∼0.28. On the other hand, the stability of a planet like g is strongly dependent on the eccentricity of the outermost planet d. This is also true in general for an additional planet orbiting in the habitable zone between planet c and d, which would likely only be stable on circular or near-circular orbit. (Note that eccentricity limits are valid for an inclination of i = 90°. If i < 90°, the limits of the eccentricities decrease as well.)

We greatly appreciate the help of William Brocas with the manuscript and we thank the anonymous reviewers for providing constructive comments.

REFERENCES

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article:

Figure 1. (a)–(l) Maximum eccentricity of a fifth massless test planet in the Gliese 581 exoplanetary system as a function of its semi-major axis and eccentricity. The integrations ran for 5000 P_test for inclinations i = 30, 50, 70° – in the range where the inclination of the debris disc around the star lies; while the argument of pericentre was varied between 0 ≤ w ≤ 360° with steps of ▵ω = 90°. Neither of these parameters changed significantly the possible stable regions in the system for a fifth planet, the a–e stability maps resemble very closely that in Fig. 2(b) in each case (Supplementary Data).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.