A&A 479, 277-282 (2008)
DOI: 10.1051/0004-6361:20078794

Dynamical evolution of the Gliese 581 planetary system

H. Beust¹ - X. Bonfils² - X. Delfosse¹ - S. Udry³

1 - Laboratoire d'Astrophysique de Grenoble, UMR 5571 CNRS, Université J. Fourier, BP 53, 38041 Grenoble Cedex 9, France
2 - Centro de Astronomia e Astrofísica da Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal
3 - Observatoire de Genève, 51 ch. des Maillettes, 1290 Sauverny, Switzerland

Received 4 October 2007 / Accepted 2 December 2007

Abstract
Context. The M dwarf Gliese 581 has recently been found to harbour two super Earths in addition to an already known close-in Neptune-mass planet. Interestingly, these two planets are considered as potentially habitable, and recent theoretical works give further credit to this hypothesis, in particular for the outermost planet (Gl 581 d).
Aims. In this paper, we address the issue of the dynamical stability evolution of this planetary system. This is important because the basic stability ensures that a 3-planet model is a physically adequate description of the radial-velocity (RV) data. It is also crucial when considering the planets' habitability because the secular evolution of the orbits may regulate their climate, even in the case where the system is stable.
Methods. We have numerically integrated the planetary system over 10⁸ yr, starting from the present fitted solution. We also performed additional simulations where i) we vary the inclination of the system relative to the line of sight, ii) assume eccentricities at the upper limit of the error bars in the radial velocity fit and where iii) we consider additional (yet undetected) outer planets. We also compute Lyapunov exponents to quantify the level of dynamical chaos in the system.
Results. In all cases, the system appears dynamically stable, even in close to pole-on configurations. The system is actually chaotic, but nevertheless stable. The semi-major axes of the planets are extremely stable, and their eccentricities undergo small amplitude variations. The addition of potential outer planets does not affect this result.
Conclusions. Consequently, from the dynamical point-of-view, a 3-planet model is an adequate description of the present RV-data set. Only a limited range of inclinations can be excluded for coplanar orbits ( $i<10^\circ$). The climate on the planets is expected to be secularly stable, thus not precluding the development of life. Gl 581 remains the best candidate for a planetary system with planets that potentially bear primitive forms of life.

Key words: planetary systems - methods: N-body simulations - celestial mechanics - stars: individual: Gliese 581 - astrobiology - stars: low-mass, brown dwarfs

1 Introduction

The M dwarf Gliese 581 has been the subject of a recent investigation with the identification of its 3-planet system. One of the planets (Gl 581 b), a Neptune-mass object orbiting the star on a 5.4-day orbit, has already been known for two years (Bonfils et al. 2005). Recently, Udry et al. (2007) have reported the discovery of two additional super Earths (Gl 581 c and d), revolving around the star in 12.9 and 83 days (see details in Table 1). Considering Gl581's luminosity, Udry et al. (2007) inferred the equilibrium temperature of both planets and conclude they may lie within the habitable zone of the star.

Table 1: Orbital parameters of the Gl 581 planetary system, as derived from the fit of Udry et al. (2007).

Detailed further modeling by von Bloh et al. (2007) and Selsis et al. (2007) addresses the habitability of the planets. Likewise, von Bloh et al. (2007) find that the greenhouse effect increases Gl 581 c temperature so much that they no longer consider the planet to be in the habitable zone. For similar reasons, Selsis et al. (2007) find that Gl 581 c's surface temperature is very likely higher than the equilibrium temperature calculated by Udry et al. (2007). However, they do not rule out habitability for this planet, as a large cloud coverage (>$75\%$) would cool down the planet enough. Conversely, both studies agree that the outermost planet (Gl 581 d) is a good candidate for habitability (although close to the outer edge of the habitable zone) and actually consider it as the better of the two candidates.

An important and unsettled issue about this system concerns its dynamical behaviour. It is first important to know whether the planetary system is dynamically stable and for which range of orbital inclinations. If verified, the basic stability of the system ensures that the model used (3 planets) is a physically adequate description of the observations (the radial-velocity measurements). If not verified, the planet detections are not necessarily invalidated (as RV periodogrammes clearly show that three coherent signals sum up at specific periods). It would instead mean that either not enough data were collected to converge toward the ``true'' parameters of the system or that the model (3 planets) is not complex enough to describe the data. Further varying the orbital inclination, we expect to find that below a given value - or, equivalently, above given masses for the planets - the system becomes unstable. Not valid physically, this range of inclinations should be rejected among the possible solutions. This partially constrains the $\sin i$ degeneracy inherent to radial-velocity detections.

Beyond the basic stability, the secular evolution of the orbits may play an important role regarding planets' habitability. All climate calculations (Selsis et al. 2007; von Bloh et al. 2007) have been done with the currently determined orbits. The secular evolution of the orbits has the potential of affecting the climate on the planets. A given planet may lie within the habitable zone but, if subject to significant eccentricity changes, it can undergo strong climate variations that eventually preclude life development. The presently determined eccentricities (Table 1) are small enough to ensure climate stability. But one needs to know which maximum values they reach due to secular perturbations.

In the present paper, we numerically investigate the secular evolution of the Gl 581 system, starting from the solution of Table 1. In Sect. 2, we study this solution (that we subsequently refer to as the nominal case). In Sect. 3, we perform other integrations, assuming different inclinations from edge-on and letting the initial eccentricities of the planets reach their maximum values within their error bars. In order to quantify the dynamical chaos in this system, we compute Lyapunov exponents for all these solutions in Sect. 4. In Sect. 5, we investigate the perturbing action of potentially additional outer planets that have not been detected yet, provided their contribution to the radial velocity signal is small enough compared to the residuals of the 3-planets fit. Our conclusions are presented in Sect. 6.

2 The nominal case

The best 3-planet orbital fit for Gl 581 is explained in Table 1. This solution with the assumptions of coplanarity and $\sin i = 1$ ( $i=90\hbox{$^\circ$ }$, an edge-one system) will constitute our nominal case. We numerically integrate this system taking $0.31~M_\odot$ for the mass of Gl 581. The integration is performed using the symplectic N-body code SyMBA (Duncan et al. 1998), which handles close encounters. The initial timestep is fixed to $2\times10^{-4}~\mbox{yr}=0.18~\mbox{day}$, i.e. 1/30 of the smallest orbital period. Symplectic integration schemes usually ensure energy conservation with 10^-6 relative accuracy as soon as the timestep is taken to $\sim$1/20 of the smallest orbital period (Levison & Duncan 1994). The integration is carried out over 10⁸ yr. On more limited timespans, we checked that taking a significantly smaller timestep does not change the result. In Fig. 1, we display the fractional errors on the total energy and angular momentum over the 10⁸ yr integration. The energy is preserved to less than 10^-7 relative accuracy. Hence we are confident in our integration. Figure 2 shows the first 10⁴ years of the integration. We see that the secular variations of the 3 planetary orbits are very regular. The eccentricities undergo quasi-periodic modulations, while the longitudes of periastra precess regularly. This solution is in fact very close to the one we can compute with a linear secular theory (Laplace - Lagrange), such as described the one by Bretagnon (1974,1990). In the linear approximation, the secular evolution of the eccentricity vectors of the 3 planets is a combination of sine and cosine terms with 3 characteristic frequencies (g_i, i=1,2,3) that are listed in Table 2. These frequencies are obtained by solving the linear secular equations for their eigenvalues. We obviously see these characteristic frequencies in Fig. 2. An interesting outcome is that these precession frequencies are much higher than in the Solar System, which do not exceed $25\hbox{$^{\prime\prime}$ }/\mbox{yr}$, and g₁ is basically the main precession frequency of the periastron of Gl 581 b and Gl 581 c. This is obviously due to the much smaller size of the system. Given the error bar on the fits of the arguments of periastra $(\omega_1, \omega_2)$, the secular motion should be detectable within $\sim$30 years, and probably less if the orbital fits get more constrained in the near future thanks to further monitoring.

	Figure 1: Fractional errors on the total energy E with respect to the initial one E0 (black curve), and on the total angular momentum H with respect to the initial one H0 (grey curve), as a function of time over the 108 yr integration.
Open with DEXTER

	Figure 2: The 104 first years of the integration of the nominal solution. The upper plots show the temporal evolution of the eccentricity of the three planets Gl 581  b), c), and  d), from left to right, respectively; the lower plots are the same for the longitude of periastron $\varpi$.
Open with DEXTER

Table 2: Precession frequencies for the nominal solution, as computed from the linear secular theory.

In Table 3, we list the maximum evolution ranges for the orbital elements of the three planets. The semi-major axes are extremely stable, revealing a regular dynamics out of any mean-motion resonance configuration. The evolution ranges of the eccentricities are narrow, so that we may claim than the system is stable with a high level of confidence. While the time span of the integration is 10⁸ yr, most of the characteristic features of the secular evolution of the orbital parameters occur on a 10⁴ yr-time scale. Therefore, even if the star is believed to be older than $2\times10^9~$yr, the current integration clearly explores all the dynamical possible outcomes of the system. Actually, due to the short orbital periods of the planets (and to the high precession frequencies), integrating the Gl 581 over 10⁸ yr is basically equivalent to integrating the Solar System over $\sim$100 Gyr!

Interestingly, the present-day eccentricity of Gl 581 c roughly corresponds to its maximum values along its secular evolution, and the eccentricity of Gl 581 d only has small variations. Hence we expect the climate of both outer planets to be secularly stable.

3 Other solutions

The nominal solution corresponds to an inclination $i=90\hbox{$^\circ$ }$ (so the lowest possible planetary masses) and to the orbital parameters of the discovery paper (Table 1). Lower inclinations and/or parameter's values slightly outside the best solution may lead to different dynamical behaviours that are worth investigating.

Table 3: Variation ranges for some orbital parameters fo the 3 planets over the 10⁸ yr integration.

In a first set of additional simulations, we assume various inclinations ranging from 0 (pole on) to  $90\hbox{$^\circ$ }$ (edge on), but still holding the initial eccentricities to their nominal values. The mass of each planet is augmented by a factor $1/\sin i$ with respect to the nominal case. In a second set of simulations we assume different inclinations and, moreover investigate the impact of eccentricities larger than in the nominal case (as less stability is only expected if the eccentricity is larger). For that set, we take the initial eccentricities for the three planets at the upper limit of their error bars (we add $1\sigma$ to the eccentricities) For both sets of integrations, we plot the width of the evolution ranges obtained over the 10⁵ yr integration for both the semi-major axis and the eccentricities of the three planets (Fig. 3).

As can be seen from Fig. 3, when the inclination decreases, the dynamical interactions increase accordingly and we expect the system to become unstable below a given inclination. As for the nominal case, the integrations are carried out over 10⁵ yr. They naturally show that both the semi-major axis and the eccentricity take a wider range of values than in the nominal case with decreasing inclinations. In the first set of integrations (nominal initial eccentricities), the system nonetheless remains stable down to $i=10\hbox{$^\circ$ }$. Almost pole-on configurations ( $i<10\hbox{$^\circ$ }$) are unstable and should be rejected from possible solutions. Although, such low inclinations are very improbable, and from the statistical point-of-view, the actual masses of the planets are probably close to those listed in Table 1.

In the second set of simulations ($1\sigma$ augmented initial eccentricities), the dynamical interactions are slightly enhanced and the semi-major axis and the eccentricity take a wider range of values than for the first set of additional simulations. The system is therefore found unstable below larger inclinations (< $20\hbox{$^\circ$ }$).

In all cases, the instability appears very unlikely. If we assume that the rotation axis of the system is randomly distributed in space, $i>20\hbox{$^\circ$ }$ occurs with a probability of 0.94. In conclusion, irrespective of its actual inclination, the Gl 581 planetary system is very probably stable.

Figure 3: Stability of the three-planet system in various configurations. The maximum variation range for the semi-major axes ( upper plots) and for the eccentricities ( lower plots) is displayed as a function of the assumed viewing inclination of the system with respect to pole-on. Each cross corresponds to a single simulation. The left plots correspond to simulations with the nominal eccentricities as initial conditions, and the right plots to simulations with eccentricities increased by $1\sigma$ relative to the error bars given in Table 1.

Open with DEXTER

4 Lyapunov exponents

Looking for variation ranges for the orbital elements is a basic tool for investigating the stability of a planetary system. A more sophisticated way for quantifying chaos is to compute Lyapunov exponents. For all simulations decribed above, we compute the maximum Lyapunov exponents (MLE) for the three planets, following the technique by Benettin et al. (1978; see also Morbidelli 2002). The exponents are computed by integrating fictitious bodies having initial conditions that are very close to those of the planets and estimating their exponential diverge rate. We coupled this algorithm with the SyMBA integrator.

When we start integrating the bodies with initial coordinate vector $\vec{p_0}$ (holding for the positions and velocities of all the bodies), we also integrate another system of bodies with identical masses, but with initial coordinate vector $\vec{p_0}+\delta\vec{p_0}$, where $\vert\vert\delta\vec{p_0}\vert\vert\ll\vert\vert\vec{p_0}\vert\vert$. After a fixed normalization time $t_{\rm norm}$, we compute the error vector $\delta\vec{p}(t_{\rm norm})$ as the difference at $t=t_{\rm norm}$ (after integration) between the coordinate vector of the fictitious bodies and of the regular bodies. We then compute

$$s_1=\frac{\vert\vert\delta\vec{p}(t_{\rm norm})\vert\vert}{\v... ...elta{\vec{p_1}}= \frac{\delta\vec{p}(t_{\rm norm})}{s_1} \cdot$$ (1)

We then use $\delta\vec{p_1}$ as a new initial error vector for the fictitious bodies relative to the coordinates vector of the regular bodies at $t_{\rm norm}$, and we iterate the above process. Each $t_{\rm norm}$, the error vector is renormalized this way, and we obtain a sequence $s_1,s_2,\ldots$ of renormalization factors. Benettin et al. (1978) proved that the MLE $\mathcal{L}$ can be computed as

$$\mathcal{L}=\lim_{n\rightarrow+\infty}\frac{1}{n~t_{\rm norm}} \sum_{i=1}^{n}s_i\equiv\lim_{n\rightarrow+\infty}L_n.$$ (2)

The result is independent of the choice of $t_{\rm norm}$, provided it is chosen small enough to avoid too large an exponential divergence. During the integration, we compute $\log L_n$ for every $t_{\rm norm}$, and we try to derive an asymptotic behaviour. Two cases can occur: i) $\log L_n$ converges towards a finite limit. Then the system is chaotic and we have reached $\mathcal{L}>0$ within our integration time; ii) or $\log L_n$ keeps decreasing monotonically at the end of the integration. This means that the orbit is very probably regular or, more precisely, less chaotic than a given level that depends on the integration time.

In practice we compute MLEs for all the individual bodies in the integration: for each orbit we compute the associated (partial) error vector, and perform individual renormalizations. Each body has its own sequence of s_i's. Note that more than the absolute values of the MLE derived, the comparison between the values derived for the individual bodies and between different integrations is more relevant. This shows clearly which orbits are the most chaotic.

In our example applications, we integrate over 10⁴ yr to compute the MLEs. $t_{\rm norm}$ has been fixed to 0.02 yr, i.e., 100 times the timestep. The progress of the computation of the MLEs in the nominal case for the three planets is plotted as a function of the integration time in Fig. 4. We see that they have all converged towards finite limits at t=10⁴ yr, showing that the orbits are actually chaotic.

The global result of MLE calculation is shown in Fig. 5, where we have computed the MLEs for the three planets for all the simulations described in Fig. 3 (stopping at t=10⁴ yr). In all cases, we obtain non-zero exponents, showing that the system is actually chaotic.

We see that the MLEs slowly increase with decreasing inclinations, showing as expected that solutions at smaller inclinations are more chaotic, due to higher planetary masses. We nevertheless note that the variation is small except for $i<20\hbox{$^\circ$ }$. The system is not much more chaotic at $i=20\hbox{$^\circ$ }$ than at $i=90\hbox{$^\circ$ }$. This confirms that there is no real dynamical constraint on the inclination. We also note that solutions with $1\sigma$ increased eccentricities are not more chaotic than those with nominal eccentricities. As a result the dynamical stability does not put any additional constraint on the planet eccentricities other than those derived from the radial velocity analysis.

From Fig. 5, it also becomes clear that the two inner planets (Gl 581 b and c) are much more chaotic than the outer one (Gl 581 d) (the exponent is smaller). In fact, the two inner planets are significantly chaotic. This does not prevent them from being stable. Actually, chaos does not necessarily mean instability. The Solar System is known to be chaotic (but on longer timescales), yet it is nevertheless stable.

	Figure 4: Progress of the calculation of the partial Lyapunov exponents as a function of the integration time, in the nominal case for Gl 581 (zero inclination, nominal eccentricities). At t=104 yr, the three exponents have stabilised.
Open with DEXTER

	Figure 5: Lyapunov exponents as a function of the inclination i, computed for all the simulations described in Fig. 3. Left plot: simulations with nominal eccentricities; Right plot: simulations with $1\sigma$increased eccentricities.
Open with DEXTER

5 Other planets

Our simulations were made with the three known planets orbiting Gl 581. However, the system may harbour additional, unknown planets. The presence of these planets may affect the stability of the whole system. There are upper limits to the presence of additional (mainly outer) planets. The maximum amplitude of the residuals in the 3-planet fits of Udry et al. (2007) is $\pm$ $2.1~\mbox{m~s}^{-1}$. Any additional planet should not generate a radial velocity with a larger amplitude, otherwise it would have already been detected. Assuming $i=90\hbox{$^\circ$ }$ and a circular orbit, this puts severe constraints on the mass m and distance d of the unseen planet. We derive

$$\frac{m}{1~M_\oplus}\leq 13.227\times\sqrt{\frac{d}{1~\mbox{AU}}}\cdot$$ (3)

This constraint holds if the unseen planet generates full-amplitude variations within the timespan of the radial velocity data, i.e., $\sim$1000 days (Udry et al. 2007). This means that the orbital period of the unseen planet must not exceed $\sim$twice this time span to account for this constraint, i.e. an orbital distance $d\leq 5.5~$AU. For more distant planets, the constraint is much weaker. In fact, this upper limit is probably already too large. Dynamically speaking, we do no expect any hypothetical planet orbiting at 5.5 AU to significantly affect the dynamics of the inner system located inside 0.25 AU. We are thus confident in the conclusions we derive below, as all potentially destabilizing configurations have been explored.

Table 4: Semi-major axis and eccentricity variation ranges for the 3 known planet of Gl 581, plus additional outer planets (see text), computed over 10⁵ yr integrations.

We thus performed new simulations, each of them with the nominal conditions, but to which we add an additional planet orbiting the star on a circular orbit at an arbitrary distance d, and with the maximum mass allowed by Eq. (3). All the integrations were carried out over 10⁵ yr. We did 5 simulations with d=5, 4, 3, 2, and 1 AU. This gives masses of 29.6, 26.5, 22.9, 18.7, and $13.2~M_\oplus$, respectively. We also added a simulation with a 1 Jupiter mass ($M_{\rm J}$) planet orbiting the star at 5 AU, as the constraint (3) is less severe at this distance. Note that this case is by far the worst possible disturbing configuration that is still compatible with the constraints. More distant companions, even massive, are less destabilizing. In a first-order approximation, the perturbing effect of a distant planet of mass m orbiting at distance d on an inner planet orbiting at distance r scales as the tidal stripping effect on the orbit, i.e. $\propto$mr/d³. Hence a $1~M_{\rm J}$ planet at 5 AU is as disturbing as a $8~M_{\rm J}$ planet at 10 AU and a $32~M_{\rm J}$ brown dwarf at 20 AU. The AO surveys would likely have already detected such a massive companion.

In all cases, the whole system appears just as stable as without any additional planet. The result concerning the stability is summarised in Table 4 where we give the semi-major axis and eccentricity variation ranges for the three known planets. The results are very similar among the 6 different integrations, even in the case of a Jovian planet, showing that the additional planet has little influence on the stability of the inner system. Moreover, Table 4 is easily compared to Table 3. The varation ranges are very similar. Therefore, we may stress that any additional outer planet that fits into the constraint of the radial velocity residuals does not affect the stability of the 3-planet system. Note that the maximum eccentricity values in Table 4 are actually slightly lower than those in Table 3. This could appear surprising, since the integration in Table 3 is made without any additional perturber. Recall, however, that it extends over 10⁸ yr instead of 10⁵ for those in Table 4. This shows conversely that if we were entending the integrations of Table 4 up to 10⁸ yr, we should expect slightly wider variation ranges. The basic conclusion nevertheless remains: the stability of the system is not affected.

6 Discussion

We have computed the secular evolution of the Gl 581 planetary system in various possible configurations. The main conclusion is that the system is almost always stable. It is stable for inclinations as low as $\sim$ $20\hbox{$^\circ$ }$ and even if the initial eccentricities are augmented by their 1-$\sigma$ error bars.

As expected for any planetary system with regular dynamics, the semi-major axes vary very little and the three planets are expected to remain at their current location with respect to the star. Meanwhile, the eccentricities of the two outer planets (both considered for habitability) reach values that are significantly above the Earth's value. Concerning Gl 581 c, the present-day eccentricity is close to its maximum value. This planet is not expected to get much farther away from its parent star and, to maintain a surface temperature cool enough to allow the presence of liquid water, a high water-cloud coverage ($\sim$75%) would be required at any time. Regarding Gl 581 d, the nominal eccentricity is non negligible ($\sim$0.12) and also found to be very stable. It is significantly above the maximum value reached by the Earth throughout its secular evolution ($\sim$0.06, see e.g. Laskar 1988) and corresponds to a 24% variation of the radiation flux received from the star between apoastron and periastron. The anomalistic season effects should therefore be strong, if not damped by the short orbital period (83 days). If we compare the periastron and apoastron values of Gl 581 d to the habitable zone calculations by Selsis et al. (2007) and von Bloh et al. (2007), we see that Gl 581 d is outside the habitable zone at apoastron but well inside at periastron. As pointed out by Selsis et al. (2007), the average stellar flux received by an eccentric orbit is enhanced by a factor $1/\sqrt{1-e^2}$ with respect to a circular orbit with the same semi-major axis. This can help maintaining Gl 581 d in the habitable zone. What we show here is that this effect is secularly permanent.

Now, if the obliquity of the rotation axis of this planet is non-zero, this should combine with the obliquity's seasonal effect and lead to climate differences between the hemispheres of this planet, much like Mars presently. The obliquity of Gl 581 d is of course unknown, but Selsis et al. (2007) and von Bloh et al. (2007) agree in claiming that, given the estimated age of the star (>2 Gyr), the rotation of Gl 581 d should already be tidally locked with the orbital motion. In that case, we would expect the obliquity to have been set to zero by tidal effects, and there should instead be climate differences between the night and day hemispheres. This could help in maintaining the day hemisphere habitable. Selsis et al. (2007) also show that tidal locking does not contradict the non-zero eccentricity of the orbit. Tides usually tend to both synchronize the rotation and circularize the orbit. The circularization time is almost always longer than the synchronization time (Hut 1981). For Gl 581 d, Selsis et al. (2007) estimate the synchronization time to 10 Myr and the circularization time to 10 Gyr, i.e. well above the present age of the system.

Acknowledgements

All the computations presented in this paper were performed at the Service Commun de Calcul Intensif de l'Observatoire de Grenoble (SCCI).

References

• Benettin, G., Galgani, L., Giorgilli, A., & Strelcyn, J. M. 1978, R. Acad. Sci. Paris, Ser. A, 286, 431 (In the text)
• Bonfils, X., Forveille, T., Delfosse, X., et al. 2005, A&A, 443, L15 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Bretagnon, P. 1974, A&A, 30, 141 [NASA ADS]
• Bretagnon, P. 1990, A&A, 231, 561 [NASA ADS]
• Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067 [NASA ADS] [CrossRef] (In the text)
• Hut, P. 1981, A&A, 99, 126 [NASA ADS] (In the text)
• Laskar J. 1988, A&A, 198, 341 [NASA ADS] (In the text)
• Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18 [NASA ADS] [CrossRef] (In the text)
• Morbidelli, A. 2002, in Modern Celestial Mechanics. Aspects of Solar System dynamics (London: Taylor & Francis) (In the text)
• Selsis, F., Kasting, J. F., Levrard, B., et al. A&A, 476, 1376 (In the text)
• Udry, S., Bonfils, X., Delfosse, X., et al. 2007, A&A, 469, L43 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• von Bloh, W., Bounama, C., Cuntz, M., & Franck, S. 2007, A&A, 476, 1365 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)