A&A 462, 769-776 (2007)
DOI: 10.1051/0004-6361:20066194

The HARPS search for southern extra-solar planets

VIII. $\mu $ Arae, a system with four planets[*],[*]

F. Pepe1 - A. C. M. Correia2 - M. Mayor1 - O. Tamuz1 - J. Couetdic3 - W. Benz4 - J.-L. Bertaux5 - F. Bouchy6 - J. Laskar3 - C. Lovis1 - D. Naef7,1 - D. Queloz1 - N. C. Santos8,1,9 - J.-P. Sivan10 - D. Sosnowska1 - S. Udry1

1 - Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, 1290 Sauverny, Switzerland
2 - Departamento de Física da Universidade de Aveiro, Campus Universitário de Santiago, 3810-193 Aveiro, Portugal
3 - Astronomie et Systèmes Dynamiques, IMCCE-CNRS UMR 8028, 77 avenue Denfert-Rochereau, 75014 Paris, France
4 - Physikalisches Institut Universität Bern, Sidlerstrasse 5, 3012 Bern, Switzerland
5 - Service d'Aéronomie, Route des Gatines, 91371 Verrières-le-Buisson Cedex, France
6 - Institut d'Astrophysique de Paris, 98bis Bd Arago, 75014 Paris, France
7 - ESO La Silla Observatory, Alonso de Cordova 3107, Vitacura Casilla 19001, Santiago 19, Chile
8 - Observatório Astronómico de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal
9 - Centro de Geofisica de Évora, Rua Romão Ramalho 59, 7002-554 Évora, Portugal
10 - Laboratoire d'Astrophysique de Marseille, Traverse du Siphon BP8, 13376 Marseille Cedex 12, France

Received 5 August 2006 / Accepted 5 September 2006

Context. The $\mu $ Arae planetary system is fairly complex, because it contains two already known planets, $\mu $ Arae b with P=640 days and $\mu $ Arae c with P=9.64 days , and a third companion on a wide, but still poorly defined, orbit.
Aims. Even with three planets in the system, the data points keep anomalously high dispersion around the fitted solution. The high residuals are only partially due to the strong p-mode oscillations of the host star. We therefore studied the possible presence of a fourth planet in this system.
Methods. During the past years we carried out additional and extremely precise radial-velocity measurements with the HARPS spectrograph. These data turned out to be highly important for constraining the many free parameters in a four-planet orbital fit. Nevertheless, the search for the best solution remains difficult in this complex and multi-dimensional parameter space. The new Stakanof software, that uses an optimized genetic algorithm, helped us considerably in this task and made our search extremely efficient and successful.
Results. We provide a full orbital solution of the planetary system around $\mu $ Arae. It turns out to be the second system known to harbor 4 planetary companions. Before this study, $\mu $ Arae b was already well known and characterized. Thanks to the new data points acquired with HARPS we can confirm the presence of $\mu $ Arae c at P=9.64 days, which produces a coherent RV signal over more than two years. The new orbital fit sets the mass of $\mu $ Arae c to 10.5  $M_{\oplus}$. Furthermore, we present the discovery of $\mu $ Arae d, a new planet on an almost circular 310 day-period and with a mass of 0.52  $M_{\rm Jup}$. Finally, we give completely new orbital parameters for the longest-period planet, $\mu $ Arae e. It is the first time that this companion has been constrained by radial-velocity data into a dynamical stable orbit, which leaves no doubt about its planetary nature. We take this opportunity to discuss naming conventions for poorly characterized planets.

Key words: instrumentation: spectrographs - techniques: radial velocities - stars: individual: $\mu $ Arae (HD 160691) - stars: planetary systems

1 Introduction

The discovery of the first planet around $\mu $ Arae (HD 160691) was reported by Butler et al. (2001). At that time, the authors announced a 1.97  $M_{\rm Jup}$-mass planet orbiting its host star in P=743 days on a very eccentric orbit (e=0.62). The acquisition of new data points during the following year called, however, for a significant correction of the preliminary parameters of $\mu $ Arae b: Jones et al. (2002) indicated a new orbital solution with P=638 days and much lower eccentricity e=0.31. From a clear trend in the radial-velocity data, the same authors also deduced the presence of a second companion, for which the orbital parameters were, however, unconstrained. Their data only set a lower limit for the period and the minimum mass, still leaving the possibility for a stellar companion open. The $\mu $ Arae system clearly deserved further investigation.

Indeed, asteroseismology observations, which have aimed at testing whether the metallicity of $\mu $ Arae was of a primordial nature or rather the product of planet engulfment (Bouchy et al. 2005), unveiled the presence of $\mu $ Arae c, a very low-mass planet of only 14  $M_{\oplus}$ on a 9-day orbit (Santos et al. 2004). The low radial-velocity wobble amplitude of 4 m s-1 was detected thanks to the extreme precision of the HARPS instrument (Mayor et al. 2003) and to the high measurement density. In fact, the stellar pulsation modes, which may attain amplitudes up to 6 to 9 m s-1 p-p on $\mu $ Arae, must be averaged out by integrating over times scales considerably longer than the pulsation period. Under these conditions, the HARPS asteroseismology run provided data points with night-to-night dispersion below 0.5 m s-1. The Neptune-mass planet $\mu $ Arae c was also confirmed by a continuous follow-up with the HARPS instrument during subsequent months.

At about the same time, McCarthy et al. (2004) published a possible orbital solution for the long-period companion, although, in their analysis, they did not take the presence of the innermost planet into account. The single Keplerian with a linear trend no longer fitted the new data points correctly. They showed that the residuals were significantly reduced when instead assuming a two-Keplerian model with the second companion of 3.1  $M_{\rm Jup}$ orbiting the parent star in 8.2 years on a eccentric orbit with e=0.57. Since the long-period planet orbit was only covered partially, large uncertainties still persisted for the orbital parameters of the outer planet. A recent dynamical analysis of the $\mu $ Arae system performed by Gozdziewski et al. (2005) shows that the published orbit was dynamically unstable and confirms that the system is still unconstrained, leaving room for longer periods and more massive companions. By investigating possible stable orbital solutions within a 3$\sigma$ range from the best fit, these authors indicate the range for the planet's orbital period at [4, 6) AU, therefore setting a lower limit at 4 AU, but no clear upper limit, and guessing the most likely value to lie around 5.5 AU.

In order to definitively constrain the orbital solutions of this complex system, a long-term follow-up was required. The many free parameters also called for a high measurement density. Last but not least, small signatures, like that of the innermost planet, required high instrumental precision and needed the effect of stellar pulsations to remain below the 1 m s-1. We decided therefore to continue a precise radial velocity follow-up of this star with HARPS and adopted a strategy for minimizing the pulsation signal. Finally, we combined the obtained HARPS data with CORALIE and UCLES data to obtain a time coverage that is as complete as possible[*].

2 Description of the data set

More than 25% of the known extrasolar planets populate systems with at least two planets. In most of these systems, the various planets were discovered sequentially, starting with the planet inducing the highest radial-velocity amplitude at short or intermediate period. The subsequently discovered planets unveiled themselves when, after a while, the residuals became high compared to the expected measurement precision and showed some structure as a function of time.

Clearly $\mu $ Arae falls under this category of systems. Both increased data amount and the enlarged time span have led to the discovery of two additional companions. While the orbit of $\mu $ Arae c is well constrained by the HARPS data, some uncertainties remain on the parameters of the third companion. With an increasing number of planets, and hence free fit parameters, the amount of possible solutions increases drastically, and the only way to constrain the orbits is to acquire many new and precise data points.

2.1 High-precision radial-velocity follow up of $\mu $ Arae with HARPS

We have pursued this objective by observing $\mu $ Arae regularly with HARPS. In the time span between May 2004 (asteroseismology run) and July 2006, we measured $\mu $ Arae on 85 different nights. One additional measurement was recorded in August 2003, before HARPS entered into operation. Apart from this older measurement, each data point consists of a short series of measurements covering a total integration time of 12 min. A typical observation scheme foresees 4 exposures of 3 min. The SNR reached in one single frame is on the order of 250 per extracted pixel, corresponding to a photon-noise limited accuracy of about 0.25 m s-1. Therefore, the photon-noise contribution of a complete observation to the total radial-velocity error lies between 0.15 m s-1 and 0.20 m s-1. The radial-velocity errors given in the data of Table 3 were, however, assumed more conservatively: To the photon-noise error, we added quadratically a constant term of 0.8 m s-1, which includes calibration errors, pulsation noise, stellar jitter, and other possible short and long-term effects. This error estimation is well-confirmed by typical dispersion values obtained on stable stars (see Lovis et al. (2005); Pepe et al. (2005), although the recent discovery of three planets around HD 69830 (Lovis et al. 2006) indicate that HARPS can provide long-term accuracy on the order of 0.3 m s-1 on very stable stars. The nightly averages of the HARPS radial-velocity data and the estimated errors are shown in Table 3. Note that we adopted an error value of 0.5 m s-1 for the asteroseismology nights. This choice was motivated by the high number of data points that averages out completely any possible oscillation, guiding and atmospheric noise. This value is also motivated by the fact that Santos et al. (2004) obtain a dispersion of about 0.45 m s-1 on the asteroseismology-night data alone.

2.2 Long-term CORALIE data

The HARPS data are of excellent quality but cover only the last two observation years. Consequently, they would not constrain possible planets with longer orbital periods. We therefore used our CORALIE data, which provide measurements as old as 8 years. The 39 CORALIE measurements acquired in 15-minute exposures show a typical photon error of 1-2 m s-1, to which we quadratically added a calibration noise of 3 m s-1. In contrast to HARPS, stellar pulsation and jitters do not significantly add to the error so can be neglected. It should be noted here that the CORALIE data were reduced with the HARPS data-reduction software which was recently adapted for CORALIE. However, some of the very first calibration frames did not pass the severe quality control, mainly because of the calibration exposures not yet being optimized in the first months of the CORALIE operation. We therefore chose to remove these measurements from our data set, these data are not presented here and were not used for the orbital fit.

Although the CORALIE data points are of lower quality compared with HARPS, new measurements were carried out recently, in order to ensure a good overlap with the HARPS data - and thus a good RV-offset determination. The CORALIE radial-velocity data and the estimated errors are shown in Table 4.

2.3 UCLES data

As mentioned above, the total amount of data points is critical for constraining the orbital parameters of a multi-planetary system unambiguously. In order to complete our data set in terms of time coverage and sampling, we integrated in this study the UCLES data on $\mu $ Arae published by McCarthy et al. (2004), to which we refer the reader for more details.

3 The $\mu $ Arae multi-planetary system

The complete data set covers about 8 years of observations and contains 171 nightly-averaged measurements of quite different precision depending on the respective instrument. The various instruments cover different periods and show poor overlap in time. One difficulty consists in assigning the correct error bars to the measurements, in particular, of one instrument relative to the other. In fact, an incorrect error estimation, and thus weights, could lead to emphasizing one time domain more than another. Even worse, an instrument that covered a short time span with high density will probe short-period planets, while it would be more sensitive to long-period planets if its observation were to cover a longer time span.

The computation of the HARPS and CORALIE measurement errors was performed as described in the previous section. For the UCLES data, we used the error values indicated in McCarthy et al. (2004). Effects from jitter or the p-mode pulsation of $\mu $ Arae, as well as long-term instrumental accuracy, cannot be taken into account precisely, especially at the level of 1 m s-1 precision. We therefore consider the used error bars for the various instruments as estimates that are realistic although not fully reliable.

The mean radial velocity provided by the three instruments is expected to be quite different between the various instruments, especially between the UCLES instrument, on one hand, and CORALIE and HARPS, on the other. The reason is that the UCLES data were recorded using the iodine-cell technique, which does not provide the stellar velocity offset, since it measures radial-velocity changes compared to some reference spectrum. Instead, the thorium technique used on HARPS and CORALIE measures the stellar radial velocity relative to the laboratory rest frame. Although these two instruments are operated in the same mode and use the same pipeline, we cannot exclude a priori some small instrumental offset due, for instance, to the different optical configuration, resolution, cross-dispersion, etc. We therefore decided to leave the offsets of the data of each instrument as additional free parameters. The two additional degrees of freedom are fully acceptable, if we consider the large number of data points.

Despite the large amount of measurements, the search for the best orbital solution remains challenging, because the irregular time sampling and the time span, in combination with the possibly large number of free fit parameters, produce a complex and rapidly varying $ \chi ^2 $ function. The major risk is to remain trapped in a local minimum that does not correspond to the solution with an absolute minimum. On the other hand, we might come up with different types of solutions showing insignificantly different $ \chi ^2 $.

In order to explore the entire parameter range and to find the best solution in terms of residuals minimization, we used the Stakanof software package recently developed at the Geneva Observatory (Tamuz 2006). Stakanof combines genetic algorithms (Goldberg 1989) for efficiently exploring the entire parameter space with algorithms allowing for fast convergence into local minima, similarly to the approach used earlier by other authors in this domain (see e.g. Laughlin 2002; Gozdziewski & Konacki 2004; Stepinski et al. 2000). Stakanof delivers a full set of solutions (same $ \chi ^2 $ within the uncertainty range defined by the measurement errors) in less than 20 min on an ordinary personal computer. However, Stakanof only fits a sum of any number of Keplerians but does not take planet-planet interaction into account.

3.1 From a two-planet to a four-planet system

In the present case we started our exploration by assuming only 2 planets in the system: the well known 1.3  $M_{\rm Jup}$ planet orbiting $\mu $ Arae in about 640 days and a planet with a long period. This was motivated by the observed trends. For simplicity, and in order to disentangle various questions, we did not fit for the short-period 14  $M_{\oplus}$ planet on the 9.6-day orbit. The orbit-search program Stakanof immediately finds the well-known planet on an eccentric orbit (e=0.24) of 643 days. A second planet is found on an orbit of 3973 days and shows eccentricity of e=0.15. As we will see later, the $\chi^2=1740$ and the weighted dispersion of the residuals of 3.5 m s-1 is quite high. In particular, the residuals to the data show a time-correlated structure, as can been seen from the upper panel of Fig. 1. We think that this solution is not satisfactory and that a fourth planet, in addition to the short-period low-mass planet, must be present in the system. Therefore, we explored a solution that foresees the presence of four planets. Besides the 640 day-period and the short-period planets and a planet with a period longer than 2000 days, we guessed the presence of an additional planet with an intermediate period of about 300 days. Despite any possible a priori knowledge we might have introduced, we decided to leave all orbital parameters for all the planets completely free and to let Stakanof search for a solution with four planets in the full parameter space.
\end{figure} Figure 1: Phase-folded radial-velocity measurements (in km s-1) and best fit for the two planets in the two-planet solution. For every planet the contribution by the other planets has been subtracted. In the upper panel showing the long-period planet, a time-correlated structure becomes evident. HARPS data are plotted in red, UCLES data in blue, and CORALIE data in green.

In contrast to our expectations, the program converges rapidly towards a class of very good solutions with weighted dispersion values around 1.75 m s-1 and $\chi^2=442$. The system contains the known planet $\mu $ Arae b at P=643.3 days, with e=0.128, K=37.8 m s-1, and M=1.68  $M_{\rm Jup}$, as well as the Neptune-mass planet ($\mu $ Arae c, with 10.5  $M_{\oplus}$ in this fit) at 9.64 days, with e=0.171 and K=3.1 m s-1. In addition, the fit confirms our guess by finding a new 0.522  $M_{\rm Jup}$ planet at 310.6 days. The eccentricity of this planet is low e=0.067 and the amplitude is K=14.9 m s-1. Finally, a clear signal is also found at a longer period. The best fit delivers a planet with M=1.81  $M_{\rm Jup}$, P=4205.8 days, e=0.098, and K=21.8 m s-1.

3.2 Exploring the parameter space

The Stakanof software also finds other solutions with an $ \chi ^2 $ only slightly higher (within 10%), but with parameters that are quite different for the outer planet. In fact the orbital period can easily vary from 3200 to 5000 days and its mass by a similar fraction. While the orbital parameters of the three inner planets essentially remain constant, the orbit of the outermost planet can vary considerably, mainly because the data points have not yet covered a full orbital cycle. We thus conclude that, although the existence of the outermost planet is not in doubt, its period and mass are fairly uncertain. Furthermore, the orbital parameters of this planet slightly influence the eccentricity and the mass of the third planet, which may in turn change the dynamic of the system. An investigation of the system stability, as in the next section, may therefore contribute to the characterization of the system or, better, help us exclude unstable systems from our solution space.

Still using Stakanof we also found a completely different solution with an $ \chi ^2 $ identical to the previously presented best fit. This system also has four planets with the inner planet at 9.64 days and the outer one at about 2741 days, but, in contrast to the previous solution, the intermediate pair of planets are trapped in a 1:1 mean motion resonance around 600 days; a pair of Trojan planets. As shown in the next section, we can exclude this solution, however, from the set of the possible ones, because it corresponds to a highly unstable configuration.

3.3 Characteristics of the $\mu $ Arae four-planet system

In Fig. 2 we present our best 4-planet Keplerian fit. The plot shows the phase-folded data and the best fit for the four planets individually. In each of the plots, the contribution by the remaining three planets was subtracted. For the detailed description of the orbital parameters of the planets, we used an N-body fit that takes the planet-planet interaction into account (Table 1). Since the observational data cover only a short period of time (about 8.5 years), the orbital elements obtained with an N-body fit do not present significant differences with respect to the one obtained with Stakanof. Figure 3 shows the global solution, i.e. the combination of the motion of the four planets as a function of time and over-plotted with the actual measurements. Figure 4 represents a zoom on the time span covered by HARPS data, where only the HARPS data are shown to emphasize how well the curve is constrained by these observations. Even by covering less than two cycles of the P=640-day planet, it becomes clear that two other planets (in addition to the very short-period planet) must be present. In particular, the flattening of the curve after the first relative maximum is a clear indication of the presence of the previously unknown companion with P=310 days.

\end{figure} Figure 2: Phase-folded radial-velocity measurements (in km s-1) and best fit for any of the four planets. For every planet the contribution by the remaining three planets has been subtracted. HARPS data are plotted in red, UCLES data in blue, and CORALIE data in green. Note that on the third panel showing the longest-period planet the HARPS data clearly indicate a flattening of the radial velocity curve. The recent HARPS data consequently constrain the orbital period of this planet.

The Neptune-mass planet $\mu $ Arae c is hardly recognized in the UCLES and CORALIE data. The semi-amplitude of 3.1 m s-1 is only constrained by the HARPS measurements. We therefore plotted only the phase-folded HARPS data in Fig. 5, after subtracting the contribution of the three other planets, in order to make the low-mass planet amplitude variation appear clearly. The scatter of the residuals is 1.41 m s-1 rms, which must be compared to UCLES and CORALIE data, with their 3.34 m s-1 rms and 6.09 m s-1 rms, respectively (Table 2). It is important to note that the HARPS data cover a time span of 530 days, even without accounting for the very first measurements recorded one year before the asteroseismology series. This demonstrates that the oscillations remained fully in phase during the covered period and that it is therefore very unlikely that they are due to spots and the activity of the star, as the absence of correlation between radial velocities and line bisector also proves. It is worth noticing that the amplitude of $\mu $ Arae c is smaller than reported in the discovery paper of K=4.1 m s-1. This difference cannot be explained only by the error bar. The main difference comes from the fact that, at the time of its discovery, the longest-period planet was poorly characterized and the P=310.6-day planet was even not known, both leading to a wrong estimation of the local "slope'' in the radial-velocity data.

We estimated that planet-planet interaction cannot be neglected in the present system, since the second and the third planets are rather close to each other. We therefore re-fitted the radial-velocity data by including interaction in our model and starting from the values computed from the Keplerian fit. The obtained orbital elements (Table 1) are very similar to the Keplerian solution. Since they are more general, we adopt these values and use them in the following as the baseline for our study of the system dynamics. For this fit we got a weighted dispersion of 1.75 m s-1 and $\chi^2=442$. The reduced $ \chi^2_{\rm red}$ of this best fit is 1.729. Note that the errors on the orbital parameters have been derived from the covariance matrix and are likely to be underestimated.

Table 1: Orbital parameters and inferred planetary characteristics derived from the best-fit orbital solution with planetary perturbations for the four planets around $\mu $ Arae at JD = 2 453 000 and keeping the inclination fixed at $90^\circ $.
$a$ &[AU] & \multicolumn{2}{c}{1.497}\\
...$a$ &[AU] & \multicolumn{2}{c}{0.09094}\\
$a$ &[AU] & \multicolumn{2}{c}{0.9210}\\
$a$ &[AU] & \multicolumn{2}{c}{5.235}\\

\end{figure} Figure 3: Global solution containing 4 planets and covering the full time range of observations. The lower plot shows the residuals of the fitted function to the actual measurements. The fit takes planet-planet interaction into account. The weighted dispersion of the residuals is 1.75 m s-1.

\end{figure} Figure 4: The figure shows a zoom on the time span covered by the HARPS data. It must be note that HARPS data constrain in an exceptional way the three shorter-period orbit. It particular, these complex curve would be incompatible even with a two-planet solution.

\end{figure} Figure 5: The upper panel shows the phase-folded HARPS data and the best Keplerian best fit to the P=9.64 days-planet after subtracting the contribution of the other planets. The lower panels shows the residuals of the data to the best fit. The dispersion of the HARPS residuals is 1.41 m s-1 rms.

4 Stability of the system

Running a simulation over one million years using the initial conditions from Table 1, we observed that the orbits of the two intermediate-period planets, $\mu $ Arae b and $\mu $ Arae d, can come as close together as 0.62 AU (Fig. 6). Consequently, the inner region of this system is subject to important gravitational interactions that may destabilize the orbits. Indeed, tracking the dynamical evolution for longer periods of time, the system is destroyed after 76 million years.

The estimated age of the $\mu $ Arae system is about 6.4 Gyr (Saffe et al. 2005), indicating that the previous orbits were not very well-determined. One reason is that the fitted parameters still present some uncertainties around the best-fitted value. This is particularly true for the outermost planet, whose orbit is not yet completely constrained because the orbital period is longer than the present 8.5 years of observational data. Moreover, additional planets may exist in the system that are disturbing the present solution, but that could not be detected yet. We should thus consider that the set of parameters given in Table 1 constitutes the best determination one can make so far in terms of minimum $ \chi ^2 $, and we will therefore search for more stable solutions in its vicinity.

Since the orbits of the innermost planet $\mu $ Arae c and the third planet $\mu $ Arae b are well-established, with small standard errors, we kept the parameters of these two planets constant. We also did not change the inclination of the orbital planes, keeping both at $90^\circ $. For the second planet, $\mu $ Arae d, we let a, $\lambda$, e and $\omega$ vary. Typically, as shown in Fig.  7, we fixed $\lambda$ and $\omega$ to specific values, and spanned the (a,e) plane for initial conditions with a step size of 0.005 AU for a and 0.05 for e, in a similar way as for the HD 202206 system (Correia et al. 2005). For each initial condition, the orbits of the three outer planets were integrated over 2000 years with the symplectic integrator SABA4 of Laskar & Robutel (2001), using a step size of 0.02 year. The innermost planet has a very small influence on the stability of the rest of the system, so removing it allows us to use a step size as big as 0.02 year. The stability of the orbit is then measured by frequency analysis (Laskar 1993,1990). Practically speaking, a refined determination of the mean motion nd, nd' of the second planet is obtained over two consecutive time intervals of T=1000 years, and the measure of the difference D = | nd-nd'|/T (in deg/yr2 in Fig.  7) is a measure of the chaotic diffusion of the trajectory. It should be close to zero for a regular solution, and high values will correspond to strong chaotic motion. In the present case a regular motion will require D < 10-6. The contour plot in Fig. 7 gives the value of $ \chi ^2 $ obtained for each choice of parameters.

Table 2: Table of the instrument-specific characteristics.
...{2}{c}{3.34} & \multicolumn{2}{c}{6.01}\\

\end{figure} Figure 6: Orbital evolution of the four planets over one million years starting with the orbital solution from Table 1 and an inclination of $90^\circ $. The panel shows a face-on view of the system; x and y are spatial coordinates in a frame centered on the star. Present orbital solutions are traced with solid lines and each dot corresponds to the position of the planet every 50 years. The semi-major axes are constant, and the eccentricities undergo small variations ( 0.09 < eb < 0.13; 0.16 < ec < 0.21; 0 < ed < 0.19; 0.08 < e(e) < 0.11).

\end{figure} Figure 7: Global view of the dynamics of the $\mu $ Arae system for variations in the eccentricity with the semi-major axis of $\mu $ Arae d. The color scale is the stability index (D) obtained through a frequency analysis of the longitude of the outer planet over two consecutive time intervals of 1000 yr. Red areas correspond to high orbital diffusion (instability) and the blue ones to low diffusion (stable orbits). The labeled contour plot shows the value of $ \chi ^2 $ obtained for each choice of parameters; the line without label indicates $ \chi ^2 = 450 $. There are two different types of blue zones where the system can be stabilized, one of them corresponding to the 2:1 mean motion resonance (for ad > 0.915 AU and ed > 0.15).

We find that the vicinity of $\mu $ Arae d (P=310.6 days) is somewhat chaotic (red regions in Fig.  7) and almost no regular motion is found. Because of the proximity of $\mu $ Arae b, the chaotic behavior was expected. Studying the (a,e) plane of initial conditions (Fig.  7), we find three main regions with qualitatively different dynamical behavior. A first region is delimited by

\begin{displaymath}e_d < e_d^{\rm (lim)} \approx 0.2 \times \left( \frac{0.937-a_d}{0.937-0.89} \right) \cdot
\end{displaymath} (1)

The fit from Table 1 lies at the edge of this region, as shown by the $ \chi ^2 $ level curves. Orbits in this particular region are chaotic, and the second and third planets experience close encounters within a few tens of million years. More regular orbits can be found in the small blue stripes with ad roughly between 0.916 and 0.92, where the system can be stabilized. Above this region, and for ad<0.915 AU, orbits are highly unstable, with most of them lost in less than 2000 years. Finally, for ad > 0.915 AU, we find the region that corresponds to the 2:1 mean motion resonance island; in Fig.  7 we actually see only part of the island, as it expands up to $a_d\approx0.975$ AU and to very high eccentricities. However, most of it is unstable. The thin blue stripes are much more regular, although still chaotic: the frequency of the main resonant argument is not very well-defined, which means that those orbits lie near the separatrix. The main resonant argument is given by:

\begin{displaymath}\theta = \lambda_d - 2 \lambda_b + \omega_d ,
\end{displaymath} (2)

with a libration period of about 30 years. Moreover, due to a localized chaotic behavior, they cross the separatrix. We should not choose such orbits for this system, but they indicate that one might be able to find stable resonant orbits nearby by changing the right angle variable to move away from the separatrix, inside the resonance.

Beside the solution listed in Table 1, we have found two more solutions with an identical $ \chi ^2 $ and residuals (see previous section). One of these solutions is very similar to the solution from Table 1 except for the outer planet (less constrained) for which we find a(e) = 5.65 AU, e(e) = 0.138, $ \lambda_{0{(e)}} = 55.85^\circ $, $ \omega_{(e)} = 39.71^\circ$, and K(e) = 23.25 m/s, corresponding to a minimum mass of 2  $M_{\rm Jup}$. The dynamical behavior of this system is very similar to the one listed in Table 1. Indeed, the modifications in the orbit of the outer planet do not change the global picture of the inner planetary system. A numerical long-term simultation is stable over about 14 million years (after which the system is destroyed), but we can slightly modify the parameters of the $\mu $ Arae d in order to stabilize the system for longer time scales.

The other solution provided by the genetic algorithm was a system of four planets, where $\mu $ Arae b and d were replaced by a pair of Trojan planets orbiting with the same period as the former $\mu $ Arae b in a 1:1 mean motion resonance. The innermost planet has the same orbital parameters as the previous two solutions, but the outermost one is closer to the inner system, orbiting just at 3.93 AU with e(e) = 0.35 and a minimum mass of about 1.3  $M_{\rm Jup}$ (see legend of Fig. 8). These parameters and the quality of the fit remain essentially the same even when we consider the planet-planet interactions. However, due to the strong interactions, if we follow this system for a few decades, we see that the orbits of the Trojan planets undergo large variations, and the system is destroyed in less than a century due to a collision between the two Trojan planets (Fig. 8). We also looked for stability in the vicinity of the Trojan planets, in a similiar study as for the other solutions (Fig. 7), but no stable configurations were found when changing ad, ed, and $\omega_d$.

During our search with Stakanof, we encountered various other solutions similar to the best fit. The main difference between them was the period of the outer most planet, ranging from 3000 to 5000 days. Given the large distance of this planet from the inner planets, its exact period will have no direct impact on the system stability. However, when fitting the radial-velocity data and minimizing $ \chi ^2 $, different orbital elements for the outermost planet will also lead to slightly different orbital elements for the inner planets. This may in turn change the system dynamics, such that it would be necessary to investigate the parameter space in the vicinity of each of these solutions. We restrict ourselves in this paper to the statement that several stable solutions exists that are compatible with the measured radial-velocity data. The presented solution (Table 1) definitively belongs to the group of the most probable ones (best fitted and more stable).

\end{figure} Figure 8: Orbital evolution of the four planets with a 1:1 resonance: inclination is equal to $90^\circ $ and ab = 1.47 AU, eb = 0.128, $ \lambda _{0b} = 36.59^\circ $, $ \omega _b = 255.68^\circ $, Kb = 62.00 m/s; ac = 0.09 AU, ec = 0.175, $ \lambda _{0c} = 176.60^\circ $, $ \omega _c = 206.04^\circ $, Kc = 3.20 m/s; ad = 1.42 AU, ed = 0.122, $ \lambda _{0d} = 235.13^\circ $, $ \omega _d = 358.80^\circ $, Kd = 27.94 m/s; and a(e) = 3.93 AU, e(e) = 0.350, $ \lambda _{0{(e)}} = 89.79^\circ $, $ \omega _{(e)} = 164.00^\circ $, K(e) = 18.81 m/s. The panel shows a face-on view of the system; x and y are spatial coordinates in a frame centered on the star. Due to the strong interactions, the orbits of the Trojan planets undergo large variations and the system is destroyed in only 92 years. A dynamical analysis of the stability around the Trojan planets also provides no alternative stable possibility.

In summary, we have found an orbital solution (Table 1), which is quite stable, but not completely locked over the time scales of the planetary system's age. By modifying the initial orbital elements, one can find stable solutions in its vicinity, leading to only slightly increased $ \chi ^2 $. There are two different zones where the system can be stabilized (Fig. 7). In one of them, the two planets $\mu $ Arae b and $\mu $ Arae d are close to the 2:1 resonance. The second one corresponds to the 2:1 resonance, although those solutions seem less likely for the moment. This result proves the existence of a solution that is compatible with the observed radial velocities, on one hand, and with the system age, on the other.

5 About the naming of planetary companions

Our results unambiguously demonstrate the presence of 3 planetary companions around $\mu $ Arae: $\mu $ Arae b, the long-known P=643 days planet; $\mu $ Arae c, the Neptune-mass planet on the short P=9.64 days discovered during the HARPS asteroseismology run; and the new planet $\mu $ Arae d, in an orbit of P=310.6 days. These planets have been fully characterized and only little uncertainty remains on their orbital parameters. We therefore named these three planets according to the time sequence of their discovery and characterization.

The orbital elements of the outermost planet as indicated by Jones et al. (2002) first, and revised by McCarthy later, are no longer consistent with the present (extended) data set. The recently acquired HARPS data have strongly contributed to constraining the system parameters. In contrast to the three well-characterized inner planets, the period remains however somewhat uncertain. Fitting the data still yields a wide range of possible solutions with very similar $ \chi ^2 $ values. The best-fit is obtained with a period of 4206 days, which is slightly longer than the time span of the observations. Without degrading the residuals significantly, it is possible to find a solution within about 20% of the period of the best solution. Nevertheless, for the first time it is possible to exclude stellar mass, since the uncertainty on the mass given by the radial-velocity data is small. In addition, the data fit provides us with a dynamically stable solution. We have therefore chosen to name this planet $\mu $ Arae e in this paper.

We propose that a similar naming strategy be adopted for all new discoveries. It is not sufficient to indicate trends in the RV-data or loose orbital elements to assess the planetary nature of the companion. There is also the risk of giving the wrong information to the theoreticians, who have to rely on assessed results (within the measurement uncertainty) to produce models that can be compared with the observations. It is worth noting that the data we had acquired until January 2006 led us first to a 4-planet solution, which turned out to be dynamically unstable. Although the stability of the system is not sufficient for proving the correctness of the solution, it can be considered as a necessary condition for making it plausible. This led us to the decision to wait for additional HARPS data, in order to provide results at a higher confidence level.

6 Conclusions

In this paper we have characterized a system of planets around $\mu $ Arae in great detail. In addition to the already well-known planet $\mu $ Arae b on a 643-day orbit and with a minimum mass of 1.68  $M_{\rm Jup}$, we confirm the presence of the Neptune-mass at P=9.64 days. The new data points constrain the mass of this planet to only 10.5  $M_{\oplus}$. We have also presented the discovery and characterized a new, to date unknown planet on a P=310.6-day orbit, $\mu $ Arae d. This planet alone accounts for the previous, anomalously high dispersion of the data points around the fitted orbit. The minimum mass of this planet is well-constrained to 0.52  $M_{\rm Jup}$. Finally we confirm the presence of a fourth companion, $\mu $ Arae e, with a likely period P=4206 days and probable mass of 1.81  $M_{\rm Jup}$. The orbital period has, however, not yet been fully constrained.

It is worth mentioning that these results could be obtained in particular thanks to the superb quality of the HARPS data and to the long time coverage delivered by the UCLES and CORALIE data, which helped for constraining the many free parameters of the orbital solution. The weighted rms of 1.75 m s,-1 obtained on 171 data points over a time span of 8.5 years and on a strongly-pulsating star, is exceptionally low and, for the moment, leaves only little space for a fifth companion.

The $\mu $ Arae system is the second known to harbor 4 planetary companions. According to new theoretical models and the experience acquired from the latest discoveries, we can expect to find many more of these systems in future. Precise radial-velocity machines like HARPS and the increasing time coverage of the various ongoing follow-up programs may provide us with many new and amazing results.

We are grateful to all technical and scientific collaborators of the HARPS Consortium, ESO Head Quarter, and ESO La Silla, who have all contributed with their extra-ordinary passion and valuable work to the success of the HARPS project. We would like to thank the European RTN "The Origin of Planetary Systems (PLANETS, contract number HPRN-CT-2002-00308) for supporting this project. This research made use of the SIMBAD database, operated at the CDS, Strasbourg, France. Finally, we acknowledge the Swiss National Science Foundation for its continuous support in our research.



7 Online Material


Table 3: HARPS radial-velocity data of $\mu $ Arae. Indicated values are weighted means of various exposures during a single night.
Julian date RV RV error
[T - 2 400 000] [ km s-1] [ km s-1]
52 906.51936 -9.2909 0.00084
53 160.72599 -9.3398 0.00050
53 161.72780 -9.3428 0.00050
53 162.72601 -9.3448 0.00050
53 163.72588 -9.3477 0.00050
53 164.72576 -9.3482 0.00050
53 165.68275 -9.3455 0.00050
53 166.78196 -9.3427 0.00050
53 167.72693 -9.3421 0.00050
53 201.61987 -9.3611 0.00085
53 202.64137 -9.3598 0.00085
53 203.61075 -9.3619 0.00085
53 204.63545 -9.3563 0.00085
53 205.56147 -9.3542 0.00085
53 206.63707 -9.3556 0.00084
53 207.66322 -9.3542 0.00085
53 216.79388 -9.3598 0.00085
53 217.60908 -9.3630 0.00085
53 218.56732 -9.3639 0.00084
53 219.67984 -9.3647 0.00091
53 229.49489 -9.3683 0.00084
53 230.57363 -9.3738 0.00083
53 232.49288 -9.3707 0.00086
53 234.66632 -9.3676 0.00089
53 237.49386 -9.3688 0.00086
53 264.52049 -9.3826 0.00088
53 265.50482 -9.3811 0.00088
53 266.49340 -9.3801 0.00089
53 267.50962 -9.3844 0.00088
53 268.51030 -9.3877 0.00088
53 269.51171 -9.3853 0.00090
53 271.48362 -9.3862 0.00089
53 272.48518 -9.3855 0.00088
53 273.48732 -9.3852 0.00088
53 274.52460 -9.3838 0.00084
53 297.52223 -9.3926 0.00086
53 298.50560 -9.3935 0.00085
53 307.50815 -9.3976 0.00084
53 308.49826 -9.3989 0.00084
53 309.50581 -9.3969 0.00085
53 310.52537 -9.3941 0.00085
53 311.51144 -9.3951 0.00084
53 312.52384 -9.3920 0.00085
53 315.50462 -9.3931 0.00085
53 400.88939 -9.3651 0.00100
53 404.89478 -9.3676 0.00092
53 410.89394 -9.3584 0.00090
53 428.86141 -9.3487 0.00084
53 466.78683 -9.3363 0.00087
53 467.84365 -9.3367 0.00085
53 470.82430 -9.3368 0.00085
53 489.88910 -9.3301 0.00093
53 490.85565 -9.3325 0.00092
53 491.84209 -9.3339 0.00089
53 492.82606 -9.3304 0.00089
53 511.85243 -9.3281 0.00085
53 515.78177 -9.3265 0.00090
53 517.79817 -9.3261 0.00089
53 520.80548 -9.3265 0.00092
53 542.68983 -9.3226 0.00086
53 543.70021 -9.3213 0.00089
53 544.76392 -9.3228 0.00088
53 545.77313 -9.3231 0.00101
53 546.76437 -9.3220 0.00088
53 547.67799 -9.3270 0.00089
53 550.67900 -9.3271 0.00091
53 551.67212 -9.3223 0.00096
53 573.65657 -9.3203 0.00091
53 576.61243 -9.3244 0.00091
53 606.63942 -9.3282 0.00088
53 608.57512 -9.3304 0.00087
53 668.48502 -9.3309 0.00086
53 670.54485 -9.3297 0.00087
53 671.52852 -9.3332 0.00087
53 672.53005 -9.3353 0.00086
53 673.54838 -9.3385 0.00087
53 809.91489 -9.3676 0.00087
53 813.91743 -9.3662 0.00091
53 829.88048 -9.3787 0.00084
53 831.91068 -9.3730 0.00084
53 834.90788 -9.3747 0.00084
53 836.90174 -9.3809 0.00084
53 862.80026 -9.3876 0.00086
53 865.79772 -9.3936 0.00086
53 869.81768 -9.3928 0.00085
53 886.76366 -9.4053 0.00087


Table 4: CORALIE radial-velocity data of $\mu $ Arae, with errors significantly higher than on HARPS data and dominated by photon noise.
Julian date RV RV error
[T - 2 400 000] [ km s-1] [ km s-1]
51 298.862816 -9.4013 0.0036
51 299.887396 -9.3949 0.0038
51 755.598346 -9.3344 0.0113
51 973.901056 -9.3838 0.0036
51 975.880104 -9.3819 0.0037
52 136.669644 -9.3727 0.0037
52 140.637318 -9.3658 0.0038
52 142.662767 -9.3802 0.0037
52 143.625355 -9.3660 0.0037
52 144.641898 -9.3768 0.0037
52 145.571119 -9.3738 0.0036
52 150.613056 -9.3620 0.0036
52 164.524007 -9.3603 0.0036
52 168.534181 -9.3541 0.0033
52 169.522324 -9.3646 0.0037
52 171.538214 -9.3557 0.0039
52 172.559328 -9.3573 0.0035
52 173.532866 -9.3560 0.0033
52 177.550376 -9.3597 0.0037
52 178.552546 -9.3544 0.0042
52 180.553288 -9.3491 0.0033
52 436.705549 -9.3390 0.0035
52 715.915800 -9.4003 0.0036
53 137.825826 -9.3568 0.0035
53 172.732248 -9.3804 0.0035
53 237.734580 -9.4093 0.0057
53 238.544750 -9.4093 0.0036
53 239.574868 -9.4124 0.0031
53 239.584890 -9.4144 0.0031
53 447.912310 -9.3706 0.0031
53 630.526063 -9.3627 0.0032
53 632.492165 -9.3770 0.0033
53 633.480736 -9.3638 0.0034
53 635.509067 -9.3725 0.0046
53 638.610152 -9.3568 0.0034
53 639.509958 -9.3639 0.0034
53 641.512632 -9.3552 0.0036
53 643.536254 -9.3674 0.0042
53 644.515441 -9.3670 0.0033
53 781.889598 -9.3823 0.0033

Copyright ESO 2007