© ESO, 2013

1. Introduction

M dwarfs are the predominant stellar population of our Galaxy (e.g. Chabrier & Baraffe 2000). Compared to our Sun, they are cooler, smaller, and lower-mass stars. These characteristics ease the detection of planets for many techniques and M dwarfs were therefore included from an early stage in planet-search samples. While the first claimed detections (e.g. around Barnard’s star – van de Kamp 1963) were later found to be incorrect (Gatewood 1995; Gatewood & Eichhorn 1973), targeting M dwarfs has since proven to be more successful.

At the forefront of planet discoveries, the radial-velocity (RV) technique was first to unveil a candidate giant planet orbiting an M dwarf. Three years after the discovery of 51Peg b (Mayor & Queloz 1995), the detection of a giant planet orbiting the M dwarf GJ 876 (Delfosse et al. 1998; Marcy et al. 1998) proved that planets could also formed around M dwarfs. GJ 876 was actually one of a few tens of M stars monitored by radial-velocity surveys, and its detection provided the early impression that giant planets could be common around late-type stars. Today, only 6 M dwarfs are known to host a planet with a minimum mass >0.5 M_Jup (see Table 1) and it has become progressively evident that there is a lower rate of occurrence for giant planets, than for Sun-like stars (Bonfils et al. 2006; Butler et al. 2004; Cumming et al. 2008; Endl et al. 2006; Johnson et al. 2007).

Improvements in the RV technique have led to the discovery of lower-mass planets down to masses of msini ≃ 1.9 M_⊕ (Mayor et al. 2009). Below 25 M_⊕, there are 8 known M-dwarf hosts which altogether host 12 such low-mass exoplanets. Hence, despite the greater difficulty in their detection, planets of low mass appear to orbit M dwarfs more frequently than giant planets (Bonfils et al. 2007). Among the detected low-mass planets, GJ 581d and GJ 667Cc are noticeably interesting because they have msini < 10 M_⊕ and receive closely the same amount of light received by Earth in our solar system (Mayor et al. 2009; Udry & Santos 2007; Delfosse et al., in prep.). Depending on their atmosphere (thickness, albedo, and chemistry), liquid water may flow on their surface – the standard criterium to define a habitable planet (Kasting et al. 1993; Selsis et al. 2007).

The transit technique has also been successful in detecting two planets transiting an M dwarf. One is GJ 436 b, a Neptune-mass planet initially detected by means of Doppler measurements (Butler et al. 2006; Maness et al. 2007) and subsequently seen in transit (Gillon et al. 2007b). Finding that GJ 436 b undergoes transits has enabled a wealth of detailed studies, such as the determinations of the planet’s true mass and radius and measurements of its effective temperature and orbital eccentricity (Deming et al. 2007; Demory et al. 2007; Gillon et al. 2007a). Most recently, the Mearth project, a search for transiting planets dedicated to late M dwarfs (Nutzman & Charbonneau 2008), has unveiled a ~6 M_⊕ planet transiting the nearby M4.5 dwarf GJ 1214 (Charbonneau et al. 2009). Like GJ 436b, it has a planetary to stellar radius ratio that is well-suited to in-depth characterizations with current observatories. Both planets are considered Rosetta stones to the physics of low-mass planets.

Anomalies in gravitational microlensing light curves can reveal planetary systems kiloparsecs away from our Sun. Most frequently, the lenses are low-mass stars of masses ≲0.6 M_⊙ and of spectral types M and K. Up to now, this technique has found 12 planets in 11 planetary systems. Among those, 7 are giant planets and 5 fall in the domain of Neptunes and super-Earths (Table 1). The technique is mostly sensitive to planets a few AUs away from their host, which, for M dwarfs, is far beyond the stellar habitable zone. The microlensing technique probes a mass-separation domain that is complementary to those studied by the RV and transit techniques and has shown evidence that, at large separations, low-mass planets outnumber giant planets (Gould et al. 2006; Sumi et al. 2010).

Ground-based astrometry applied to planet searches has been cursed by false positives, of which van de Kamp’s attempts around Barnard’s star are probably the most famous examples (Gatewood & Eichhorn 1973; van de Kamp 1963). Fifty years ago, van de Kamp first claimed that a 1.6 M_Jup planet orbits Barnard’s star every 24 years. Over the subsequent decades, he continued to argue that a planetary system orbited around the star (van de Kamp 1982), despite growing evidence of systematics in the data (e.g. Gatewood & Eichhorn 1973; Hershey 1973). Radial-velocity and astrometric data have now completely excluded the van de Kamp planets (Benedict et al. 2002; Gatewood 1995; Kürster et al. 2003), but Barnard’s star has been far from the only target with false astrometric detections.

Nevertheless, astrometry has proven useful in confirming the planetary nature of a few radial-velocity detections, and removing from the planet sample the low-mass stars seen with an unfavorable inclination (e.g. Halbwachs et al. 2000). Moreover, thanks to HST/FGS astrometric observations, GJ 876b was the second exoplanet for which a true mass determination was obtained (Benedict et al. 2002), soon after the detection of the transiting planet HD 209458 b (Charbonneau et al. 2000; Henry et al. 2000; Mazeh et al. 2000).

To complete the view of planetary-mass objects formed at the lower end of the main sequence, we highlight the cases of a ~5 M_Jup companion detected with RV measurements around the young brown dwarf Cha Hα 8 (Joergens & Müller 2007) and the ~5 M_Jup companion imaged around another young brown dwarf 2M1207 (Chauvin et al. 2004). The protoplanetary disks of both brown dwarfs were most likely not massive enough to form such massive objects, since observations show that protoplanetary disk masses scale at most linearly with the mass of the star. Both 2M1207b and Cha Hα 8b therefore probably formed in a similar way to stars rather than in protoplanetary disks.

Table 1 lists the known M-dwarf hosts and their procession of planets. For each planet, it gives the basic characteristics and a reference to the discovery papers. In total, there are 35 planets in 28 planetary systems.

Planets orbiting M dwarfs formed in a different environment from those around solar-type stars, and therefore reflect a different outcome of the planetary formation mechanism. The mass of the proto-planetary disk, its temperature and density profiles, gravity, the gas-dissipation timescale, etc. all change with stellar mass (e.g. Ida & Lin 2005). Following our construction of the HARPS spectrograph for ESO, our consortium has been granted 500 observing nights with the instrument spread over six years. We chose to dedicate 10% of this guaranteed time to characterize the planetary population of stars with masses <0.6 M_⊙.

This paper reports on our six-year RV search for planets around M dwarfs and the outline is as follows. We first describe our sample (Sect. 2) and present the RV dataset collected (Sect. 3). We next perform a systematic analysis to search for variability, long-term trends, and periodic signals (Sect. 4). We close that section with an automated Keplerian multi-planet search. For all signals detected with sufficient confidence, we apply a suite of diagnostics to disentangle Doppler shifts caused by bona fide planets from Doppler shifts caused by stellar surface inhomogeneities (Sect. 5). Section 6 presents the detection limits for individual stars in our sample and Sect. 7 collates these limits together to estimate both the survey sensitivity and the frequency of planets orbiting M dwarfs. Section 8 summarizes our results and presents our conclusions.

2. Sample

Our search for planets orbiting M dwarfs originates from RV observations that started in 1995 with the 1.93 m/ELODIE spectrograph (CNRS/OHP, France). This former program was designed to determine the stellar multiplicity of very-low-mass stars (Delfosse et al. 1999), as well as to detect their most massive planets (Delfosse et al. 1998). In 1998, we initiated a similar program in the southern hemisphere, with the 1.52 m/FEROS spectrograph (La Silla, Chile). With FEROS’s early settings, we were unsuccessful in improving on its nominal precision of ~50 m/s. Nevertheless, we used our experience gained in these observations to perform a superior sample selection for observations with HARPS, from which we removed spectroscopic binaries and rapid rotators for which precision RV measurement is more difficult.

Our HARPS sample corresponds to a volume-limited list of M dwarfs closer than 11 pc, with declinations δ < + 20°, magnitudes brighter than V = 14 mag and projected rotational velocities of vsini ≲ 6.5  km s^-1. We removed known spectroscopic binaries and visual pairs with separation <5′′ (to avoid light contamination from the unwanted component). We however encountered a few spectroscopic binaries and rapid rotators that have been undiscovered before our observations. We list them in Table 2 and discard them from the sample presented here. We note that we also dismiss GJ 1001, Gl 452.1, and LHS 3836. The first two stars were initially included in our volume-limited sample (with π = 103 and 96 mas, respectively – Gliese & Jahreiß 1991) and have had their parallax revised since, now placing them beyond 11 pc (with π = 76.86 ± 3.97 and 88.3 ± 3.7 mas – Henry et al. 2006; Smart et al. 2010). LHS 3836 was initially included based on its V magnitude in Gliese & Jahreiß (1991)’s catalog but our first measurements were indicative of a much fainter brightness.

Table 3 lists the 102 stars selected for the sample. Their coordinates, proper motions, and parallaxes are primarily retrieved from the revised Hipparcos catalog (van Leeuwen 2007). A fraction of the parallaxes, which are unavailable in the Hipparcos database, were obtained from van Altena et al. (1995), Henry et al. (2006), Reid et al. (1995), and the 4th Catalog of Nearby Stars (CNS4 – Jahreiss, priv. comm.). V-band magnitudes were taken from Simbad and infrared J- and K-band magnitudes from 2MASS (Cutri et al. 2003). We used the empirical Delfosse et al. (2000)’s mass-luminosity relationship together with parallaxes and K-band photometry to compute the mass of each star. Infrared K-band photometry and (J − K) colors were converted to luminosities with Leggett et al. (2001)’s bolometric correction.

We also indicate in Table 3 the inner and outer limits to the distance of the habitable zone using the recent-Venus and early-Mars criterions, respectively, and Eqs. (2) and (3) of Selsis et al. (2007). The boundaries of the Habitable Zone are uncertain and depend on the planet’s atmospheric composition. Extra-solar planets found close to these edges should therefore meet more stringent conditions to be inhabitable. For more detailed considerations, we refer the reader to more comprehensive models (e.g. Selsis et al. 2007).

Our sample is composed of the closest neighbors to the Sun. Nearby stars tend to have large proper motions and the projection of their velocity vector may change over time, by up to few m/s/yr (Kürster et al. 2003; Schlesinger 1917). We therefore report the value of their secular acceleration in Table 3.

To portray our sample, we show its V-mag and mass distributions in Fig. 1. For both distributions, the average (resp. median) value is plotted with a vertical straight (resp. dashed) line. The magnitudes and masses of planet hosts are also marked with vertical ticks on top of the histograms. The target brightness spans V = 7.3 to 14 mag with a mean (resp. median) value of 11.25 mag (resp. 11.43 mag). The stellar mass ranges from 0.09 to 0.60 M_⊙ with an average (resp. median) value of 0.30 M_⊙ (resp. 0.27 M_⊙). The lower count seen in the 0.35–0.40 M_⊙ bin remains unexplained but may be due to statistical fluctuations. Interestingly, we note that our sample covers a factor of about six in stellar mass, while the mass step between our typical M dwarf (~0.27 M_⊙) and the typical Sun-like star (~1 M_⊙) corresponds to a factor of less than four. This means that planetary formation processes depending on stellar mass could lead to larger observable differences across our sample than between our M-dwarf sample and Sun-like stars.

There are overlaps between our sample and others that similarly target M dwarfs to search for planets. Among the data for these other samples, we found published RV time-series for Gl 1, Gl 176, Gl 229, Gl 357, Gl 551, Gl 682, Gl 699, Gl 846, and Gl 849 in Endl et al. (2006, hereafter E06), for Gl 1, Gl 229, Gl357, Gl 433, Gl 551, Gl 682, Gl 699, Gl 846, and Gl 849 in Zechmeister et al. (2009, hereafter Z09), and, a few others in detection papers, such as Gl 176 (Endl et al. 2008), Gl 832 (Bailey et al. 2009), and Gl 849 (Butler et al. 2006). When possible, we compare our results to these time series and, for completeness, provide an additional comparison in Appendix A.

3. Observations

To gather RV observations for the sample described above, we used the HARPS instrument (Mayor et al. 2003; Pepe et al. 2004), a spectrograph fiber fed with the ESO/3.6-m telescope (La Silla, Chile). This instrument covers the 3800−6800 Å wavelength domain, has a resolution R ~ 115   000, and an overall throughput (instrument+telescope) greater than 6%. It is enclosed in a vacuum vessel, is pressure controlled to ± 0.01 mbar and thermally controlled to ±0.01 K. This ensures that there is a minimum instrumental shift of the position of the spectral image on the CCD or, practically, a RV drift ≲0.5 m/s/night. To determine the origin of the instrumental RV drift, the fiber can be illuminated with a thorium-argon (ThAr) lamp at any time. HARPS also offers a second fiber that can be illuminated with ThAr light simultaneously while the scientific fiber receives starlight. This mode avoids the need to record frequent calibrations between scientific exposures and can correct the small instrumental drift that occurs while the stellar spectrum is being recorded. Since the instrumental drift during the night is small, this mode is only used when a sub-m/s precision is required. Our observational strategy for M dwarfs aims to achieve a precision of ~1 m/s per exposure for the brightest targets. We therefore chose not to use the second fiber and relied instead on a single calibration done before the beginning of the night. As the science and calibration spectral orders are interlaced, avoiding the second fiber eludes light cross-contamination between the science and calibration spectra. This can be a source of noise for the blue-most spectral orders, where we can practically reach only low signal-to-noise ratios for M dwarfs. We were particularly interested in clean Ca ii H&K lines because they are useful diagnostics of stellar activity (see Sect. 5).

From the first light on HARPS on February 11, 2003 to the end of our guaranteed time program on April, 1 2009, we recorded 1965 spectra for the M-dwarf sample, for a total integration time of 460 h.

We computed RVs by cross-correlating the spectra with a numerical weighted mask following Baranne et al. (1996) and Pepe et al. (2002). Our numerical mask was generated by adding several exposures taken on Gl 877, a bright M2 star of our sample. Co-addition of spectra requires knowledge of their relative Doppler shifts. We computed RVs for Gl 877 spectra after the first iteration of the template and then re-built the template more precisely. We achieved convergence after only a few iterations. The numerical mask counts almost 10 000 lines and most of the Doppler information in the spectrum. No binning is done.

The RV uncertainties were inferred from the Doppler content of individual spectra, using the linear approximation of a Doppler shift (Eq. (12) – Bouchy et al. 2001). This formula gives more weight to spectral elements with higher derivatives because these are more sensitive to phase shifts and provide a greater contribution to the total Doppler content. We note that we do not sum the Doppler content of individual spectral elements over the whole spectrum. The derivative of the spectrum has a higher variability against noise than the spectrum itself and, for low signal-to-noise ratio, doing so would over-estimate the RV precision. To compute more realistic uncertainties, we instead applied the formula directly to the cross-correlation profile, which has a ~30 times higher S/N than the individual spectral lines (see Appendix A in Boisse et al. 2011). To account for the imperfect guiding (~30 cm/s) and wavelength calibration (~50 cm/s), we quadratically added 60 cm/s to the Doppler uncertainty.

As a trade-off between exposure time and precision, we chose to fix the integration time to 900 s for all observations. We obtained a precision¹σ_i ~ 80 cm/s from V^mag = 7 − 10 stars and σ_i ~ 2.5^{(V − 10) / 2} m/s for V^mag = 10 − 14. Our internal errors σ_i (composed of photon noise + instrumental errors) are shown in Fig. 2, where we report, for all stars with more than six measurements, the mean σ_i (blue filled circle) as a function of the star’s magnitude, with error bars corresponding to σ_i’s dispersion. For comparison, Fig. 2 shows the observed dispersions (or external error) σ_e for all stars with more than six measurements (black squares, changed to triangles for clipped values). The σ_e are the observed weighted rms and are related to the χ² value (1)We discuss the difference between the internal and external errors in Sect. 4.

Our RV time series are given in the solar system barycentric reference frame and can be retrieved online from CDS. Prior to their analysis, the time series are corrected from the secular RV changes reported in Table 3. We also show the time series in Fig. 3 (only available on-line) after subtraction of half the min+max value for a better readability.

4. Data analysis

Several planet-search programs have presented a statistical analysis of their survey (e.g. Cumming et al. 1999, 2008; Endl et al. 2002; Murdoch et al. 1993; Walker et al. 1995; Zechmeister et al. 2009). Statistical tests are often applied to the time series in order to appraise the significance of trends or variability. The time series are then searched for periodicities and, if a significant periodicity is found, the corresponding period is used as a starting point for a Keplerian fit. Statistical tests are again applied to decide whether a sinusoidal or a Keplerian model is a good description of the time series. In this section, we follow the same strategy and add a heuristic method based on genetic algorithms to the systematic search for Keplerian signals.

4.1. Excess variability

After computing the RVs, the Doppler uncertainties (σ_i), and the RV dispersion (σ_e) for each star, the first step is to test our time series for variability. In Fig. 2, we have reported both σ_i and σ_e as a function of stellar magnitudes for all stars with more than six measurements. Apart from one star (Gl 447, which has only six RVs), we have σ_e ≳ σ_i. The σ_e have a lower envelope that matches the σ_i’s envelope for V^mag = 10 − 14 but is slightly higher (~2 m/s) in the brighter range V^mag = 7 − 10.

To test whether the observed RVs vary in excess of our internal errors, we first compare the expected variance (the mean internal error) to the measured variance by applying a F-test, which evaluates a probability P(F) using the F-value (e.g. Zechmeister et al. 2009). As another test of variability, we also compute the for the constant model and , the probability of having given the σ_i. For both the F-test and the χ²-test, a low probability means that the photon noise, the calibration uncertainty, and the guiding errors are insufficient to explain the observed variability. One has then to invoke an additional source of variability, from unaccounted noise to planetary companions.

We report σ_i, σ_e, P(F), χ², and P(χ²) in Table 4 and change the P(F) and the P(χ²) values to boldface when lower than 1%, i.e. when they indicate a confidence level for variability higher than 99%. Using this criterion, the F probabilities (resp. the χ² probabilities) indicate that 45% (resp. 63%) of our sample displays an excess variability. When focusing on stars with V^mag = 7 − 10, all stars but two are found to be more variable than expected according to P(F), and all stars according to P(χ²). The reason is that our external error never reaches the ~70 cm/s threshold estimated for the brighter range of our sample (but rather 1.5−2 m/s), and is dominated by photon noise for a third of the sample in the V^mag = 10 − 14 range.

4.2. Trends

We now examine the time series to help us discern any possible trends which may correspond to incomplete orbits and betray the presence of long-period companions. For each star, we fit a slope α (RV = αt + β) to the RV data and evaluate the value of that model.

To know determine the slope is a closer description of the data than the no slope model, we use two statistical tests. First, we use the F-test to gauge whether a lower than is a sufficient improvement to justify an additional free parameter (two for the slope model against one for the constant model).

In addition, because the F-test statistics are ill-behaved for non-normally distributed uncertainties, we use a less-sensitive test based on bootstrap randomization. We generate virtual RV time series by shuffling the original data, i.e., we keep the same observing dates and attribute to each date a measurement randomly chosen among the other dates, without repeating a given measurement. On each virtual time series, we adjust a slope and compute its value. The fraction of virtual data sets with then gives us the false-alarm probability (FAP) for the slope model. For all stars, apart from those with six or fewer measurements, we generate 1000 virtual time series. Since this method only probes FAPs greater than O(1 / N !), we limit the number of trials to N ! for stars with fewer measurements.

Table 4 gives the slope coefficient α as well as the P(F) and FAP values for all time series. The P(F) and the FAP values are reported in boldface when below a threshold of 1%, to indicate the high confidence level of the slope model.

Among our sample and according to the FAP values, we find that 15 stars have time series that are more accurately described by a slope than a simple constant. They are Gl 1, LP 771-95A, Gl 205, Gl 341, Gl 382, Gl413.1, Gl 618A, Gl 667C, Gl 680, Gl 699, Gl 701, Gl 752A, Gl 832, Gl 849, and Gl 880. We also see that, while LP 771-95A, Gl 367, Gl618A, Gl 680, and Gl 880 display smooth RV drifts, the other 10 stars seems to obey a more complex variability. According to P(F) values, we find that the same 15 stars plus 8 more have a significant chi-squared improvements when we fit a slope. They are Gl54.1, Gl 250B, Gl 273, Gl 367, Gl 433, Gl 551, Gl 674, and Gl 887.

4.3. Periodicity

In our search for planets, variability selects the stars to focus on and trends reveal the still uncomplete orbits. Our next step in the search for planetary candidates is to look for periodic signals. The classical diagnostic of periodic coherent signals in unevenly spaced time series is the Lomb-Scargle periodogram (Lomb 1976; Scargle 1982) or, in order to account for the unknown system’s mean velocity, its generalized version: the floating-mean periodogram (Cumming et al. 1999).

We therefore compute generalized Lomb-Scargle periodograms for all our time series with at least 6 measurements. We follow the formalism developed in Zechmeister & Kürster (2009), and choose a power normalization, where 1.0 means that a sinusoidal fit is perfect (χ² = 0) and 0.0 means that a sinusoidal fit does not improve χ² over a constant model. We calculate the FAPs as we did in our trend analysis, with a bootstrap randomization (which is superior to a F-test). We create 1000 virtual time series by shuffling the original data set. For each individual data set, we compute a periodogram and locate the highest power. Considering all periodograms from all virtual data sets, we compute the distribution of the power maxima. The power value that is above 99% of all trial powers is then assumed to represent the 1% FAP. More generally, we attribute a FAP to the maximum power value found in the original data set by counting the fraction of the simulated power maxima that have a larger value.

We show all periodograms in Fig. 16 (only available in electronic format). The periods corresponding to the highest power, the corresponding signal semi-amplitude, the χ², as well as the associated FAPs are reported in Table 5. As previously, we boldface the significant signal detections, i.e. the periods with a power excess that has a FAP < 1%.

In ideal cases, long-term variability only affects the long-period range of the periodogram. However, the sparse sampling of our time series can cause a power excess caused by a long-period signal to leak to shorter periods. Removing the long-term significant trends cleans the periodogram and may reveal periodic signal at shorter periods. Thus, we also compute the periodogram for the time series after removing the slope adjusted in the previous section, as well as the corresponding FAPs (labeled FAP₂ in Table 5). We record a noticeable change for only Gl 618A and Gl 680. No power excess remains in their periodograms after subtraction of the slope, meaning that their RV variation are mostly linear drifts.

Our periodicity analysis identifies 14 stars that have power excesses with FAP₂ < 1%, namely Gl 176, G 205, Gl 273, Gl 358, Gl 388, Gl 433, Gl 479, Gl 581, Gl 667C, Gl 674, Gl 832, Gl 846, Gl 849, and Gl 876. The star Gl 887 additionally has a FAP approaching our 1% threshold.

After identifying possible periodicities, we use the candidate periods as starting values for Keplerian fits. We search the residuals again for periodic signals, computing periodograms (see Fig. 17 in the online material) and locating power excesses with FAP < 1%. We find a probable second periodic signal for Gl 176, Gl 205, Gl 581, Gl 674, and Gl 876 (Table 6). We note that Gl 667C, Gl 832, and Gl 846 have a FAP approaching our 1% threshold. Although not shown here, a third iteration only reveals coherent signals for Gl 581 and, possibly, Gl667C. A fourth iteration found a signal for Gl 581 only.

4.4. Keplerian analysis

Even in its generalized form, the Lomb-Scargle periodogram is optimal for sinusoidal signals only. Eccentric orbits spread the spectral power over various harmonics, decreasing the power measured at the fundamental period, and fitting a Keplerian function to provide an obvious improvement. Compared to a periodogram search, it can detect planets with high eccentricities in the presence of a higher level of noise. It is often not used because it is non-linear in some of its orbital parameters. Traditional non-linear minimizations can only converge to a solution close to a starting guess, which is possibly far from the global optimum. Applying non-linear minimization from many starting guesses becomes quickly impracticable when the number of planets increases. Keplerian functions are indeed transcendent and evaluating the radial velocity at a given time therefore requires an integration that is computationally expensive. Finally, the sequential approach outlined above requires data of higher S/N for systems with several planets on commensurable orbits (i.e. with RV pulls of similar amplitudes).

To work around these shortcomings, we use a hybrid method based on both a fast non-linear algorithm (Levenberg-Marquardt) and genetic operators (breeding, mutations, and cross-over). The algorithm was developed by one of us (D.S.) in an orbit analysis software named YORBIT. We give a brief overview of this software here, but defer its detailed description to a future paper (Ségransan et al., in prep.). While a population of typically 1000 solutions evolves, the top layer of YORBIT evaluates the performances of the different minimization methods and tunes their time allocation in real time. Genetic algorithms efficiently explore the parameter space on large scales. They thus score well and are given more CPU time during the early phase of the minimization. Once promising solutions are found, non-linear methods converge much more efficiently toward the local optima. Hence, when new solutions arise outside local minima, the non-linear methods are given more CPU time. This heuristic approach has proven to be very efficient and the solution to multi-planet systems can be found in only a few minutes, on common desktop computers.

We search for planets using YORBIT on all stars with more than 12 measurements, neglecting planet-planet interactions at this point. We scale the complexity of the tested models with the number of measurements. Although in principle 5N+1 RVs are enough to obtain a Keplerian fit to an N-planet system, we wish to minimize the number of spurious solutions and arbitrary require at least 12 RVs per planet in the model. Hence, we use a one-planet model for 12 RVs, both a one-planet and a two-planet model for stars with more than 24 RVs, and four different models (1−4 planets) for stars with more than 48 RVs. To complement those models, we also use the same suite of models with the addition of a linear drift. We allow YORBIT to run for 150 s per planet in the model.

As n our periodicity analysis, evaluating the credence of the model is essential. The χ² of solutions necessarily improves for more complex models as the number of degrees of freedom increases, and we wish to evaluate whether this improvement is either statistically significant or occurs by chance. As previously, we generate virtual datasets by shuffling the original data and retaining the dates. We create 1000 virtual datasets and refit all tested models 1000 times with YORBIT. For each star and model, we obtain 1000 χ² values and count the number that are below the χ² measured for the original data. This gives the FAP of a particular model, relative to the no-planet hypothesis. A model is considered to improve χ² significantly when fewer than 1% of the virtual trials give as low a χ².

When we find a significant model, more complex models are then compared to that model, and no longer to any simpler models. We consider that signals in this model are detected (i.e., for instance, we assume the system is composed of two planets if that model is a two-planet model). To generate the virtual datasets, we use the residuals around the best-fit solution for the new reference model (i.e., in our example, the residuals around the two-planet model). Shuffling the residuals (and retaining the dates), we again create 1000 virtual datasets and fit the more complex models using YORBIT. The number of χ² in these virtual trials that are lower than the best-fit χ² obtained on the actual RVs then indicates the FAP for the complex models, which can be compared to that of the simpler model.

We report the parameters of the best-fit solutions for each star and model in Table 7. In Tables 8 and 9, we compare the FAP of the models to those of selected simpler models. Except for the column with χ², all numbers are FAP, given in percent. Each row lists a model, which is compared with simpler models given in the different columns. The FAP are in boldface when the more complex model is found to be a statistically significant improvement (FAP < 1%) over the simpler model. Of the 43 stars with more than 12 measurements, 19 have RVs that can be well-modeled by one or more planets. Among them, we recover all stars with a probable RV periodicity except Gl 680 which has fewer than 12 measurements and was not tested. In the following section, we discuss the interpretation of these candidate signals.

5. Interpretation

The above analysis reveal the Keplerian signals in our time series, but Doppler shifts do not always correspond to planets. To vet a RV variability against stellar activity, we make use of several diagnostics. The shape of the cross-correlation function (CCF) in particular is often informative. While its barycenter measures the radial-velocity, its bisector inverse slope (BIS – Queloz et al. 2001) and its full width at half maximum (FWHM) can diagnose inhomogeneities at the stellar surface. Alternatively, spectral indices based on Ca ii H&K or H_α lines can also help to identify inhomogeneities in the stellar chromosphere and photosphere, respectively (Bonfils et al. 2007). Finally, we obtain photometric observations of a few stars to check whether plages or spots could produce the observed Doppler changes.

5.1. Planetary systems

Among the stars with clear periodic/Keplerian signals, we recover of course several planets that were previously known. In total, 14 planets are known to orbit 8 M dwarfs in our sample. Nine were found by this program (Gl 176 b, Gl 433 b, Gl 581 b, c, d & e, Gl 667C b, and Gl 674 b), one more was found in 1998 by our former program on ELODIE (Gl 876 b), two were detected by concurrent programs and confirmed by this program (Gl 876 c & d), and two were detected by concurrent programs and confirmed in this paper (Gl 832 b and Gl 849 b).

• Gl 176:

from HET radial-velocity data, Endl et al. (2008) proposed that it hosts a msini = 24 M_⊕ planet in a 10.2-d orbit. However, we found that our HARPS data are incompatible with the existence of such a planet (Forveille et al. 2009), finding evidence instead of a lower-mass planet with a shorter period (msini = 8 M_⊕; P = 8 d). As in the case of Gl 674 (see below), Gl 176 is a moderately active M dwarf. We also observe a second periodic signal (P ~ 40 d), which has marginally higher power than the 8-d signal in our periodogram. Thanks to photometric observations and spectroscopic indices measured in the same spectra, we identified the 40-d signal with a spot rather than a second planet (Forveille et al. 2009). We note that our systematic Keplerian search for planets converges to solutions with different periods and very high eccentricities. This is mostly because the signal associated with the 40-d period might be poorly described by a Keplerian function with a large eccentricity. The method converges to the same periods as the periodogram analysis when we restrict the range of eccentricities to <0.6.

• Gl 433:

this nearby M2V dwarf is rather massive (M_⋆ = 0.48 M_⊙) for our sample. The periodogram of our HARPS RVs shows a clear power excess at a 7.2-d period. We failed to find a counterpart to that signal in our activity indicators. In addition, on the basis of the intensity rather than the variability of the Hα and Ca ii lines, the star seems to have a weak magnetic activity, and most probably a low rotational velocity. It is therefore likely that a planet revolves around Gl 433 every 7.2 days. A χ² minimization of HARPS RVs lead to a minimum mass of msini = 6 M_⊕ for that planet. We note that Zechmeister et al. (2009), who measured UVES radial velocities for Gl 433, found no significant periodicity. The semi-amplitude of our solution 3.5 ± 0.4 m/s translates to an rms of 5.0 ± 0.6 m/s, which is nonetheless compatible with the 4.4 m/s rms reported in Z09. In addition, we found of Z09 that the RVs are compatible with our data, provided that we use a model composed of a Keplerian plus a low-order polynomial to fit the merged data sets. We refer the reader to Delfosse et al. (in prep.) for a detailed description.

• Gl 581:

this system contains at least four planets, the discoveries of which were reported in the three papers: Bonfils et al. (2005), Udry & Santos (2007), and Mayor et al. (2009). We also performed a stability analysis for the system in Beust et al. (2008), which was updated in Mayor et al. (2009). Composed of one Neptune-mass planet and three super-Earths, the system is remarkable because it includes both the first two possibly habitable exoplanets (Gl 581c & d – Selsis et al. 2007; von Bloh et al. 2007) and the lowest-mass planet known to date (Gl 581e – msini = 1.9 M_⊕ ). In addition, we found that stability constrains its configurations. In the coplanar cases, inclinations lower than ~40° induce too strong interactions between “b” and “e”, and “e” is ejected after a few hundred thousands years. A lower bound to the inclinations translates to upper bounds for planetary masses. For instance, Gl 581e, would not be more massive than ~3 M_⊕, if the system were coplanar. In 2010, Vogt et al. proposed that two additional planets orbit Gl 581, one having msini = 3.1 M_⊕ and being in the middle of the habitable zone, between Gl581 c and d. However, we demonstrate elsewhere that these planets do not exist (Forveille et al. 2011b).

• Gl 667C:

this is a M2V dwarf that we have intensively observed. We find several coherent signals in RV data and estimate the rotation period using FWHM measurements of the cross-correlation function. Fitting a one-planet model plus a ~1.8 m   s^-1   yr^-1 linear drift to account for the A+B stellar-binary companion to Gl 667C, converges toward a minimum mass of a super-Earth (msini = 5.9 M_⊕) on a short-period orbit (7.2 day). Adding one more planet to the model causes the fit to converge toward a ~180 day period and a very eccentric solution, while the power excess identified in the periodogram of the first model residuals is located around 90 days. A finer analysis actually interprets that second signal as a possible harmonic of a (half-)yearly systematic affecting a few data points (Delfosse et al., in prep.). Filtering the signals of two planets plus a linear drift reveals another candidate planet (P₃ = 28 d; msini = 3.9 M_⊕). This candidate receives about 90% of the amount of light received by Earth in our solar system and we speculate that the planet is habitable (see Delfosse et al., in prep., for a detailed description).

• Gl 674:

only 4.5 pc away from our Sun, this M2.5 dwarfs hosts at least one low-mass planet (msini = 11 M_⊕; P = 4.7 d – Bonfils et al. 2007). Although a second periodic signal exists for Gl 674 (P₂ ~ 35 d), analysis of spectroscopic indices and photometric observations shows that this additional signal originated from stellar surface inhomogeneities. Today and with additional measurements, after subtracting the 4.7-d periodic signal and computing a periodogram of the residuals we were able to identify a power excess at a period of ~25 d instead of 35 d. If this excess were due to a planet, this would be a super-Earth in Gl 674’s habitable zone. However, the semi-amplitude K of the Keplerian orbit of that second planet is ~3.8 m/s, which is significantly above the residuals around the 2007 combined fit (rms ~ 80 cm/s). There are two different interpretations of this apparent inconsistency: either the 2007 solution excludes the present solution and today’s 25-d periodicity is spurious or, the 2007 fit absorbed both the 35- and the 25-d signals, simultaneously.

The Keplerian analysis presented in Sect. 4.4 measured a lower significance of only ~94.6% for the second signal. Nevertheless we which to test this result using additional diagnostics. Restricting the data set to the 2007’ RVs, we attempt to fit a 1 planet+sine wave model instead of two Keplerians and found almost equally low residuals (rms ~ 1.1 m/s). This is strongly inconsistent with the presence of an additional planet with a 3.8 m/s semi-amplitude and either a short or moderate orbital period. On the other hand, periodograms of the Hα and Ca ii H+K indices continue to peak at a ~35-d period (Fig. 4), indicating that the signal remains coherent for these indicators. A decorrelation between spectral indices and RVs is then hard to explain. We conclude that the case for an additional planet is not strong enough with the present data set, and that we will have to analyze data gathered after Apr. 1 2009, before drawing any stronger conclusions.

• Gl 832:

a decade-long RV campaign with the Anglo Australian Telescope (AAT) has revealed that Gl 832 hosts a long-period companion with a period that is almost as long as the survey itself (~9.5 yr – Bailey et al. 2009). The best-fit to AAT data implies a minimum mass msini = 0.64 M_Jup. Our HARPS data do not span as long a time period and while they do confirm with a high confidence level the long-period RV variation, they cannot confirm the planetary nature of Gl 832b in themselves. Together with the AAT data², our HARPS RVs refine the orbit of Gl 832b. A Keplerian fit using YORBIT converges to msini = 0.62 ± 0.05 M_Jup and P = 3507 ± 181 d (Fig. 5). Our Keplerian analysis (Sect. 4.4) identifies a possible second signal with a 35-d period, which reaches a 99% significance level in neither our periodicity nor our Keplerian analysis. Some power excess is seen around ~40 d in the BIS periodogram, which is uncomfortably close to the possible 35-d periodicity. We continue to monitor Gl 832 to clarify whether a second periodic signal is present and assess its true nature.

• Gl 849:

it was discovered in 2006 that this M3V dwarf hosts a Jupiter-mass companion (Butler et al. 2006). The RV variation is clearly seen in our HARPS observations, and has no counterpart in our activity indicators (based on the shape of the cross-correlation function or Hα and Ca ii spectral indices). Our observations confirm that Gl 849 hosts a Jupiter-mass companion.

Fitting a Keplerian orbit to HARPS observations alone enables us to determine a minimum mass of msini = 1.17 ± 0.06 M_Jup and a period P = 2165 ± 132 d. In addition to the Keck RVs however, one planet is insufficient to explain all the RV motion. As already suspected from Keck data, a long-term change is superimposed on the first periodic signal. We therefore fit the merged data set with a 1 planet+drift, a 2 planets, and a 2 planets+drift model and calculate their respective FAPs. We find that a model more complex than the 1 planet+a drift model is not justified. For that model, our best-fit solution () corresponds to a Jupiter-mass planet (msini = 0.99 ± 0.02 M_Jup; P = 1852 ± 19 d; e = 0.04 ± 0.02) plus, a RV drift with a slope of − 3.84 m/s/yr (Fig. 6; Table 10). Since Gl 849 is a nearby star (d = 8.77 ± 0.16 pc), the long-period massive companion makes it an excellent target for astrometric observations and direct imaging with high angular resolution.

• Gl 876:

this system was known to harbor planets before our observations started, when it was the only planetary system centered on an M dwarf. The first giant planet found to orbit GJ 876 was detected simultaneously by members of our team using the ELODIE and CORALIE spectrographs (Delfosse et al. 1998) and by the HIRES/Keck search for exoplanets (Marcy et al. 1998). The system was later found to host a second giant planet in a 2:1 resonance with GJ 876b (Marcy et al. 2001). The third planet detected around GJ 876 was the first known super-Earth, Gl 876d (Rivera et al. 2005). Because the 2:1 configuration of the two giant planets leads to strong interactions, the orbits differ significantly from Keplerian motions. To model the radial velocities, one has to integrate the planet movement with a N-body code. This lifts the sini degeneracy and measures the masses of the giant planets. A full N-body analysis was first performed for GJ 876 by Laughlin & Chambers (2001) and Rivera & Lissauer (2001). From an updated set of Keck RVs, Rivera et al. (2005) found the third planet only because planet-planet interactions were properly accounted for in the fitting procedure. Those authors still had to assume coplanar orbits to assign a mass to each planet. Bean et al. (2009) then combined the Keck RVs with HST astrometry to both measure the masses in the coplanar case and to measure the relative inclination between planets “b” and “c”. Most recently, in Correia et al. (2010), we used our HARPS data and the published Keck measurements to model the system and measure the relative inclination of both giant planets (<1°), relying on RVs only. We also analyzed the dynamical stability and showed that the libration amplitudes are smaller than 2° thanks to a damping process acting during the planet formation.

5.2. Activity-dominated variations

We now gather the “active” cases, which are stars that tested positively for periodicity and/or a Keplerian signal, and their measurement variability correlates with an activity indicator. We do not show all diagnostics for each star, but instead select the most illustrative. A statistical discussion of all activity indicators will be presented in a separate paper (Bonfils et al. 2010, in prep.).

• Gl 205:

a periodogram of the velocities identifies excess power around 32.8 d. Our Keplerian search with YORBIT, for a either single planet or a 1 planet+drift model, converges toward either a similar period or 0.970 d, an alias with the 1 day sampling. We find indications that the variation is intrinsic to the star from both spectral indices and photometric observations. Considering the whole dataset, we identify excess power around 33 d in periodograms of Hα and Ca ii H+K indices, though those peaks are not the highest. Restricting the dataset instead to one observational season, a strong power excess around the period 33 d dominates the periodogram (Fig. 7). We also note a strong correlation between the Hα and Ca ii H+K indices (their Pearson correlation coefficient is 0.97), suggesting that both variations originate from the same surface homogeneity. in addition, photometric monitoring of Gl 205 reports a similar period for the stellar rotation (33.61 d – Kiraga & Stepien 2007). The observed RV modulation is most probably due to surface inhomogeneities, which remain stable over one season but not over several.

• Gl 358:

we gathered 28 measurements for Gl 358, which show evidence of significant variability with a periodicity of ~26 d. In the RV time series periodogram, we also identify power excess at the first harmonic of that period (~13 d). The RV modulation is well described by a Keplerian orbit. However, we also observe similar variability in the FWHM of the CCF as well as a possible anti-correlation between the RV and spectral line asymmetry (see Fig. 8), as measured by the CCF bisector span (Queloz et al. 2001). The Pearson’s correlation for BIS and RV is −0.40, and rises to −0.67 for the 2007 measurement subset. The photometric monitoring of Kiraga & Stepien (2007) found a rotational period of 25.26 d for Gl 358, which is consistent with a scenario where a stellar surface inhomogeneity such as a spot or a plage produces the RV change, rather than a planet.

• Gl 388 (AD Leo):

the periodogram of its RV time series displays significant power excess at short periods, with two prominent peaks at 1.8 d and 2.22 d, which are consistent with the 2.24-d rotational period reported by Morin et al. (2008). These are 1-d aliases of each other, the latter being slightly stronger. Here, the bisector span demonstrates that stellar activity is responsible for the variation. Its periodogram shows a broad power excess at short periods, and it is strongly anti-correlated to RV (with a Pearson’s correlation coefficient of − 0.81 – see Fig. 9). Correcting for the BIS-RV correlation by subtracting a linear fit does decrease the rms from 24 m/s to 14 m/s, but leaves some power excess around ~2 d.

• Gl 479:

we observe significant power excesses in the RV time series at two periods, ~11.3 d and 23–24 d, with the shorter period being roughly half the longer one. The RVs vary with an amplitude of ~27 m/s and an rms of 4.13 m/s. Modeling that RV variability with Keplerians converges toward two planets with very similar periods (23.23 and 23.40 d), which would clearly be an unstable system. Gl 479 shares its M3 spectral type with Gl 674 and Gl581, which we use as benchmarks for their moderate and weak magnetic activity, respectively. From a spectral index based on the Ca ii H&K lines, we find that Gl 479 has a magnetic activity that is in-between that of Gl 674 and Gl 581. Neither the bisector nor the spectral indices show any significant periodicity or correlation with RVs. However, we complemented our diagnostic for stellar activity with a photometric campaign with the Euler Telescope (La Silla). The periodogram of the photometry shows a maximum power excess for the period 23.75 d, which is similar to the RV periodicity. The photometry phase folded to the 23.75-d period varies with a peak-to-peak amplitude of 5% and complex patterns. We cannot ascertain the rotational period of Gl 479, but a 5% variability can explain the observed RV variations including the very slow rotation. We are unable to correlate the photometry with a RV variation because they were not taken at the same time but, phasing both to a 23 d period, we found a 0.24-phase shift that is consistent with a spot. The observed RV variability is therefore probably due to magnetic activity on Gl 479 rather than to planets.

• Gl 526:

we observe a RV periodicity with a power excess close to 50 d and a FAP approaching 1%. As expected for either a moderate or long period, we find corresponding changes in neither the BIS nor FWHM. As for Gl 674, spectral indices or photometry are then more informative. For Gl 526, we find that RV is weakly correlated to Hα and that the Ca II H+K index varies with a clear 50-d period (Fig. 11). Since the observed period is similar for the calcium index and the RV shift, we infer that the RV changes as due to magnetic activity rather than a planet.

• Gl 846:

we observe a RV variability with significant power excesses in the periodogram at several periods (7.4, 7.9, and 10.6 d), plus their aliases with the 1-d sampling, near 1-d. The BIS is well correlated with the RV (Fig. 12). As for both Gl 358 and Gl 388, Gl 846 also clearly displays stellar intrinsic variability rather than a planetary-companion Doppler shift.

5.3. Unclear cases

• Gl 273:

this M3.5 dwarf shows RV variability and has a RV drift according to the χ²-probability test (and approaches FAP = 1% for the permutation test). The RV periodogram indicates that there is a significant power excess at a period of ~434 d, and YORBIT find a good solution for a planet with P ~ 440 d, with or without a supplementary drift. However, good solutions are found in cases of an uncomfortably high eccentricity and, most importantly, a poor phase coverage. Among the activity indicators, only the Hα index has a power excess for long periods but with a different period (~500–600 d). After a linear drift has been subtracted, 2 RV points clearly deviate from the general trend (with BJD = 2 454 775 and BJD = 2 454 779). They have a value 8–10 m/s lower than the residual mean. If we fit a drift again to the original data, considering all but these 2 points, the 440-d power excess disappears from the periodogram of the residuals. This suggests that, if these two particular points were outliers, the 440-d period could be an alias between a long-term RV drift and the one-year sampling. However, a direct inspection of Gl 273 spectra, cross-correlation functions, and spectral indices give no reason to exclude those values. We conclude that a firm conclusion on Gl 273 would be premature and will require more data.

• Gl 887:

formally significant models are found for one and three planets but all converge to solutions with very high or unrealistic eccentricities. Most probably, the RV variability of Gl 887 cannot be described by a Keplerian motion, and our automatic search became confused. In addition, we do not find any significant periodicity or correlation with RVs in our diagnostics.

6. Detection limits

While the previous sections focused on signal detection, the present one attempts to determine upper limits to the signals we were unable to detect. For individual stars, the upper limit translates into which planet, as a function of its mass and period, can be ruled out given our observations. For the sample, all upper limits taken together convert into a survey efficiency and can be used to measure statistical properties (Sect. 7).

To derive a period-mass limit above which we can rule out the presence of a planet, given our observations, we start with the periodogram analysis presented in Sect. 4.3. We hypothesize that the time series consists only of noise. As for the FAP calculations, we evaluate the noise in the periodogram by generating virtual data sets. The virtual data are created by shuffling the time series while retaining the observing dates (i.e. by bootstrap randomization). For each trial, we repeat our computation of the periodogram. This time however, we do not look for the most powerful peak. We instead keep all periodograms and build a distribution of powers, at each period. For a given period then, the power distribution tells us the range of power values compatible with no planet, i.e. compatible with our hypothesis that the time series consists only of noise. In the same manner as for the FAP computation, we can also evaluate the probability that a given power value occurred by chance just by counting the fraction of the power distribution with lower values.

Once we know that the power distribution is compatible with no planet, we inject trial planetary orbits into the data. In this paper, we restrict our analysis to circular orbits and therefore fix the eccentricity and argument of periastron to zero. Nevertheless we note that eccentricities as high as 0.5 do not have a strong effect on much the upper limit to planet detection (Cumming & Dragomir 2010; Endl et al. 2002). We thus add sine waves, choosing a period P, a semi-amplitude K, and a phase T. We explore all periods computed in the periodograms, from 1.5 days to 10 000 days, with a linear sampling in frequency of step 1/20 000 day^-1, and for 12 equi-spaced phases. At each trial period, we start with the semi-amplitude of the best sine fit to the original time series K_obs, and compute the new periodogram power p_sim for that period. We next increase the semi-amplitude K until p_sim reaches a value with a probability as low as or lower than 1%, if the data were noise only. On the one hand, we impose a power threshold on all our 12 trial phases and obtain a conservative detection limit. On the other hand, we average this power threshold over our 12 trial phases and obtain a statistical (or phase-averaged) detection limit. Eventually, for both detection limits, we convert the K semi-amplitude to planetary mass³, and orbital period to orbital semi-major axis⁴, using the stellar mass listed in Table 3.

The method described above was previously applied by Cumming et al. (1999, 2008) and Zechmeister et al. (2009). We note a small difference between Cumming et al.’s and Zechmeister et al.’s approaches. On the one hand, Cumming and collaborators add the trial sine wave to normally distributed noise and choose as a variance the rms around the best sine wave to the observed data. On the other hand, Zechmeister and collaborators do not choose to make trial versions for the noise and consider the observed data as the noise itself (to which they add the trial orbit). We choose Zechmeister et al.’s approach because we believe it is more conservative when, in some cases, the noise departs from a normal distribution.

We report both conservative and phase-averaged upper limits to the 98 time series with more than 4 measurements (an example is shown in Figs. 13 and 14 for Gl 581 and we group the figures of all stars in Figs. 18 and 19, which is only available online). When a periodic variation is attributed to stellar magnetic activity, we know the main variability is not due to a planet. We therefore apply a first order correction to the RV data by subtracting the best-fit sine function. We choose a simple sine wave rather than a more complex function (such as a Keplerian) because the fundamental Fourier term is the least informative choice, and hence the most conservative. When instead, we identified that the RV variability is due to one or more planets, we are interested in the upper-limit imposed by the residuals around the solution. We therefore subtract the best-fit Keplerian solution to the time series before computing the upper limit. For the multi-planetary systems Gl 581, Gl 667C, and Gl 876, we compute the periodogram once with raw time series and once with the residuals around the full orbital solution (with all detected planets). Since the giant planets in Gl 876 system undergo strong mutual interactions, we use a N-body integration to compute the residuals (Correia et al. 2010). For Gl 674 and Gl 176 (which appear to be affected by both planet- and activity-induced variations), we use a 1 Keplerian+sine model to fit the RVs. Finally, in cases when we observed no variability, variability without periodicity, or periodic variability without a well-identified cause (planetary versus magnetic activity), we use the raw time series to compute periodograms and upper limits.

7. Planet occurrence

Interpreted together, the phase-averaged detection limits calculated for individual stars give the survey efficiency, which is eventually used to correct for the detection incompleteness and derive the precise frequency of occurrence of planets. Although the statistical analysis of our survey is the purpose of a companion paper (Bonfils et al., in prep.), we provide here a first aperçu.

In an msini-period diagram, we collate the phase-averaged detection limits computed in Sect. 6 and, for each period, count the number of excluded planets more massive than a given mass (with our 99% criterion). This synthesis is shown in Fig. 15 together with the iso-contours for 1, 10, 20, 30, 40, 50, 60, 70, 80, and 90 stars. We also overlay the planet detections including all planets of multi-planetary systems. This diagram is especially useful for comparing the planetary occurrence for different domains of masses and periods. For instance, for periods P < 100 days, our sample counts 1 host star (Gl 876) with 2 giant planets (msini = 0.5−10 M_Jup) but 7 super-Earths (msini = 1−10 M_⊕). Although our survey is sensitive to Gl 876b-like planets for 92 stars, short-period super Earths could be detected for only ~5–20 stars. It is therefore obvious that super-Earths are much more frequent than Gl 876b-like giants.

For more precise estimate, we delineate regions in the mass-period diagram and approximate the frequency of planets by the ratio f = N_d / N_{⋆, eff}, where N_d is the number of planets detected in that region and N_{⋆, eff} is the number of stars whose detection limits confidently exclude planets of a similar mass and period. We evaluate N_{⋆, eff} with a Monte Carlo sampling: we draw a random mass and period within the region delineated (assuming a log-uniform probability for both quantities), use the msini-period diagram of Fig. 15 to give a local estimate of N_⋆, eff, and, with many trials, compute an averaged N_⋆, eff value. Table 11 reports the number of detections, which are the average values for N_{⋆, eff} and the corresponding occurrence of planets for different regions chosen in Fig. 15.

The numbers of stars that we presented in Table 11 confirm that planets are increasingly abundant toward lower-mass and longer-period planets.

Finally, we estimate η_⊕, the frequency of habitable planets orbiting M dwarfs. For each star, we use the locations for both the inner (a_{HZ, in}) and outer (a_{HZ, out}) edges of the habitable zone computed in Sect. 2 and consider habitable planets to have msini between 1 M_⊕ and 10 M_⊕. To evaluate the sensitivity of our survey to habitable planets, we compute a new N_{⋆, eff} for the habitable zone with a Monte-Carlo approach again. We draw random masses from the range 1–10 M_⊕ and a random semi-major axis between a_{HZ, in} and a_{HZ, out}, choosing a log-uniform probability for both the mass and the semi-major axis. We screen the detection limits computed previously and increment N_{⋆, eff} when the trial falls above that threshold. After normalizing N_{⋆, eff} by the number of trial we found that N_{⋆, eff} = 4.84. As among our detections two planets (Gl 581d and Gl 667Cc) fall in the habitable zone, we have N_d = 2 and therefore .

Alternatively, we measure that 11 (resp. 3) stars of our sample have a time-series that is precise enough to detect planets with the same mass and period as Gl 581d (resp. Gl667Cc), which leads to a very similar estimate of η_⊕ (~42%).

8. Conclusions

We have reported on HARPS guaranteed-time observations for a volume-limited sample of nearby M dwarfs. The paper is based on a systematic analysis of the time series for 102 M dwarfs. It analyzes their variability and searches for possible trends, periodicities and Keplerian signals.

We have found significant periodic signals for 14 stars and linear trends for 15. We have recovered the signal for 14 known planets in 8 systems. In particular, we have confirmed the detection of 2 giant planets, Gl 849 b and Gl 832 b, and confirmed that an additional long-period companion is probably orbiting Gl 849. We have analyzed the RV periodicityand various stellar diagnostics of 8 other stars, and found evidences that the observed RV variation originate from stellar surface inhomogeneities for all but one (Gl 273). We have found a periodic RV variation in Gl 273 time series without any counterpart in activity indicators, though the phase coverage is too poor for a robust detection.

Our search for planets with HARPS has detected 9 planets in that sample alone, and a total of 11 planets when counting 2 M dwarfs from another complementary sample. Our detections include the lowest-mass planet known so far and the first prototype of habitable super-Earths. They are the fruit of slightly less than 500 h of observing time on a 3.6-m telescope, which nowadays is considered as a modestly sized telescope.

Beyond individual detections, we have also reported a first statistical analysis of our survey. We have derived the occurrence of M-dwarf planets for different regions of the msini-period diagram. In particular, we have found that giant planets (msini = 100−1000 M_⊕) have a low frequency (e.g. f < 1% for P = 1−10 d and for P = 10 − 100 d), whereas super-Earths (msini = 1 − 10 M_⊕) are likely very abundant ( for P = 1 − 10 d and for P = 10 − 100 d). We also found that the frequency of habitable planets orbiting M dwarfs . Considering that M dwarfs dominate the stellar count in the Galaxy, this estimate indicates that the frequency of habitable planets in our Galaxy is significantly high. In addition, for the first time, η_⊕ is a direct measure and not a number extrapolated from the statistics of more massive planets.

Many refinements are of course possible. For instance, back in Fig. 1 we indicated with vertical ticks (above the histograms) the V and mass values for the known planet-host stars included in our sample. It is striking that all planet-host stars are found in the brightest and more massive halves of the two distributions. This is reminiscent of what we observe between solar-type stars and M dwarfs. Solar-like stars have many more detected planets, but the observational advantages and disadvantages of targeting these stars compared to M dwarfs are difficult to weight. On the one hand, solar-type stars are brighter and have a lower jitter level (e.g. Hartman et al. 2011), but on the other hand they are more massive and the reflex motion induced by a given planet is weaker. To determine whether we face an observational bias or a true stellar mass dependance in the formation of planets, we need to evaluate the detection efficiency for all mass ranges, which is the purpose of a upcoming companion paper (Bonfils et al., in prep.).

Online material

Appendix A: Comparison with published time series

Appendix A.1: Variability

Compared to other published time-series, we measured lower dispersions for all M dwarfs apart from Gl 846 and known planet-host stars. Gl 1 is not found to be variable in E06 and Z09 but their dispersion is limited by a higher photon noise (~2.6 m/s, against σ_e = 1.9 m/s in our case). We report some variability of Gl 229 but at a level of <1.9 m/s, while the variability reported in E06 and Z09 implies a jitter of 3.9–4.7 m/s. The slightly smaller dispersion we observe for Gl 357 (σ_e = 3.2 m/s) compared to 3.7 m/s and 5.3 m/s for E06 and Z09, respectively, might not be significant given our small number of observations (6). For Gl 551, we measured only a dispersion slightly smaller (3.3 against 3.6 m/s). We observe significantly lower variability for Gl 682 (1.8 against 3.6 m/s) and Gl 699 (1.7 against 3.4 and 3.3 m/s), and a larger dispersion for Gl 846 (5.6 against 3.0 m/s). Although different time-spans, epochs of observations and activity levels at those epochs could explain the different dispersions for individual stars (as is certainly the case for Gl 846 – see Sect. 5.2), that we measure a smaller dispersion for most comparison stars most likely reflects the superior performance of the HARPS spectrograph.

Appendix A.2: Trends

As in this paper, Z09 reported non-significant slopes for Gl 357 and Gl 682 and significant slopes for Gl 1, Gl 551, and Gl 699 (although in our case Gl699 is attributed a significant trend only by the F-test). Nonetheless, the slopes that Z09 reported for Gl 551 and Gl 699 seem different and they moreover found a significant trend for Gl 229, whereas we do not. Time series have also been published for Gl 832 and Gl 849 as they were they were identified as likely hosts of orbiting planet (Bailey et al. 2009; Butler et al. 2006). For both stars, the planetary reflex motion clearly dominates the radial velocity signal so we discard them from any quantitative comparison. In Table A.1, we compare the slopes of linear fits to the time series in Z09 and to those of this paper. We note that the significant differences most often reflect a signal more complex than a simple linear drift.

Appendix A.3: Periodicity

Among stars with identified periodicity in RV data, Gl 832, Gl 849, and Gl 876 have time series published to report on detected planets. The periodicities we have found for those three stars are similar to their planetary orbital periods. Only Gl 876d is undetected with our automated procedure because one has to do a full N-body integration to subtract properly the signal induced by planets “b” and “c”. Besides known planet hosts, Z09 also report an absence of periodicities for Gl 229, Gl 357, Gl 433, and Gl 682, and significant periodicities for Gl 551 and Gl 699. Our results and Z09 are therefore in contradiction for three stars: Gl 433, Gl 551 and Gl 699. We noted in Sect. 5.1 that, for Gl 433, the RVs reported by Z09 and in this paper are not incompatible provided that the merged data set is fitted by a model composed of one planet plus a quadratic drift. The about one-year periodicity found for Gl 551 by Z09 and Endl & Kürster (2008) led these authors to attribute the signal to an alias of a low frequency signal with the typical one-year sampling. After Endl & Kürster (2008), the low frequency signal is believed to be caused by a clustering of points that are both blue-shifted and have a higher Hα index than other points in the time series. This putative activity signal might not be seen in our time series because it represented by only 24 measurements, against 229 in Z09. Finally, the periodicity found for Gl 699 is also attributed to activity, with a clear counterpart in Hα filling factor. In addition, if that activity signal is not seen in our time series, it is likely because it is represented by only 22 measurements, compared to 226 for Z09.

---

Acknowledgments

We express our gratitude to Martin Kuerster, who refered this paper. His comments were most useful and significantly improved the manuscript.

References

All Tables

All Figures