Bayesian search for low-mass planets around nearby M dwarfs. Estimates for occurrence rate based on global detectability statistics

Due to their higher planet-star mass-ratios, M dwarfs are the easiest targets for detection of low-mass planets orbiting nearby stars using Doppler spectroscopy. Furthermore, because of their low masses and luminosities, Doppler measurements enable the detection of low-mass planets in their habitable zones that correspond to closer orbits than for Solar-type stars. We re-analyse literature UVES radial velocities of 41 nearby M dwarfs in a combination with new velocities obtained from publicly available spectra from the HARPS-ESO spectrograph of these stars in an attempt to constrain any low-amplitude Keplerian signals. We apply Bayesian signal detection criteria, together with posterior sampling techniques, in combination with noise models that take into account correlations in the data and obtain estimates for the number of planet candidates in the sample. More generally, we use the estimated detection probability function to calculate the occurrence rate of low-mass planets around nearby M dwarfs. We report eight new planet candidates in the sample (orbiting GJ 27.1, GJ 160.2, GJ 180, GJ 229, GJ 422, and GJ 682), including two new multiplanet systems, and confirm two previously known candidates in the GJ 433 system based on detections of Keplerian signals in the combined UVES and HARPS radial velocity data that cannot be explained by periodic and/or quasiperiodic phenomena related to stellar activities. Finally, we use the estimated detection probability function to calculate the occurrence rate of low-mass planets around nearby M dwarfs. According to our results, M dwarfs are hosts to an abundance of low-mass planets and the occurrence rate of planets less massive than 10 M$_{\oplus}$ is of the order of one planet per star, possibly even greater. ...


INTRODUCTION
In recent years, planets have been discovered around the least massive stars, M dwarfs, in a diversity of different configurations with widely varying orbital properties ⋆ E-mail: mikko.tuomi@utu.fi; m.tuomi@herts.ac.uk and masses (e.g. Endl et al. 2006;Bonfils et al. 2013a, and references therein). For instance, there are several highmultiplicity systems around M dwarfs consisting of only low-mass planets that can be referred to as super-Earths or Neptunes, such as those orbiting GJ 581 (Bonfils et al. 2005;Udry et al. 2007;Mayor et al. 2009) 1 , GJ 667C Delfosse et al. 2013), and GJ 163 (Bonfils et al. 2013b;. Recent precision velocity surveys have also revealed the existence of more massive planetary companions orbiting nearby M dwarfs (e.g. Rivera et al. 2010;Anglada-Escudé & Tuomi 2012) showing that such companions do exist, but not in abundance (Bonfils et al. 2013a;Montet et al. 2014), and are less common than for K, G, and F stars (Endl et al. 2006). However, the most interesting planetary companions around these stars are the low-mass ones that orbit their hosts with such separations that, under certain assumptions regarding atmospheric properties, they can be estimated to enable the existence of water in its liquid form on the planetary surfaces (e.g. Selsis et al. 2007;Kopparapu et al. 2013). Planets of this type -sometimes called habitable-zone super-Earths -are easier to detect around M dwarfs than around more massive stars because the planet-star mass-ratios give rise to signals with sufficiently high amplitudes, and the shorter orbital periods allow for more orbital phases to be sampled in data covering a fixed length of time, to enable their detections (e.g. Mayor et al. 2009;Anglada-Escudé et al. 2013;. Recently, accurate estimates for the occurrence rate of planets in the Kepler's field have been reported in several studies (e.g. Howard et al. 2012;Dressing & Charbonneau 2013;Morton & Swift 2013). One of the most interesting features in the Kepler sample is that the occurrence rate of planets around stars appears to increase from roughly 0.05 planets per star around F2 stars to 0.3 per star around M0 dwarfs (Howard et al. 2012), although the functional form of this relation is far from certain. This increase applies to planets with orbital periods below 50 days because of the available baseline of the Kepler data. While Kepler will be able to provide occurrence rates for longer orbital periods, possibly up to 200-300 days, radial velocity surveys will be needed to probe the occurrence rate of planets on orbits longer than that. Moreover, unlike planets around more massive K, G, and F stars that have been targeted by the Kepler space-telescope in abundance, M dwarfs are not bright enough to be found in comparable numbers in the Kepler's field, which makes it difficult to estimate the occurrence rates and statistical properties of planets around such stars in detail.
According to Dressing & Charbonneau (2013), the Kepler's sample contains 3897 stars with estimated effective temperatures below 4000 K, out of which 64 are planet candidate host stars with a total of 95 candidate planets orbiting them. Dressing & Charbonneau (2013) concluded that with periods (P ) less than 50 days, the occurrence rate of planets with radii between 0.5 R⊕ < rp < 4R⊕ is 0.90 +0.04 −0.03 planets per star; with radii between 0.5 R⊕ < rp < 1.4 R⊕ is 0.51 +0.13 −0.06 planets per star, although this estimate might be underestimated as much as by a factor of two (Morton & Swift 2013); and that the occurrence rate 1 We note that the number of planets around GJ 581 is uncertain with different authors reporting different numbers from three to six (see Vogt et al. 2010Vogt et al. , 2012Gregory 2011;Tuomi 2011;Baluev 2013).
of planets with rp > 1.4 R⊕ decreases as a function of decreasing stellar temperature. Furthermore, the occurrence rate of planets appears to decrease heavily between 2 and 4 R⊕, which is indicative of overabundance of planets with low radii and therefore low masses (Morton & Swift 2013). These findings challenge the results obtained using radial velocity surveys that should be able to detect planets with similar statistics, although the comparison with Kepler's results is difficult due to the challenges in comparing populations described in terms of planetary radii and minimum masses in the absense of accurate population models for planetary compositions and therefore densities. The estimates based on transits detected by using the Kepler telescope might also be contaminated by a false positive rate of ∼ 10% due to astrophysical effects such as stellar binaries in the background (Morton & Johnson 2011;Fressing et al. 2013).
Far fewer planets around M dwarfs are known from radial velocity surveys of such stars (e.g. Bonfils et al. 2013a, who reported nine planet candidates in their sample). However, the ones that are known are among the richest and the most interesting extrasolar planetary systems in terms of numbers of planets, their orbital spacing and dynamical packedness, and their low masses (e.g. Mayor et al. 2009;Rivera et al. 2010;Anglada-Escudé & Tuomi 2012;Anglada-Escudé et al. 2013;. To a certain extent, this lack of known planets around M dwarfs is due to observational biases arising from the fact that early radial velocity surveys did not target low-mass stars because of the difficulties in obtaining sufficiently high signal-to-noise observations due to a lack of photons in the V band to enable high quality radial velocity measurements. Another reason was thatbased on a sample size of unity -Solar-type stars were considered more promising hosts to planetary systems. This observational bias is also likely caused by the fact thatin comparison with stars of the spectral classes F, G, and K -massive giant planets are not as abundant around M dwarfs (Bonfils et al. 2013a), and the planets that exist, if they indeed do exist, are likely so small that they induce radial velocity signals that have amplitudes comparable to the current high-precision measurement noise levels, which makes their detection difficult at best. Bonfils et al. (2013a) reported estimates for the occurrence rates of planetary companions orbiting M dwarfs based on radial velocity measurements obtained by using the High Accuracy Radial velocity Planet Searcher (HARPS) spectrograph. According to their results, super-Earths with minimum masses between 1 and 10 M⊕ are abundant around M dwarfs with an occurrence rate of 0.36 +0.25 −0.10 for periods between 1 and 10 days and 0.52 +0.50 −0.16 for periods between 10 and 100 days, respectively. Furhtermore, they reported an estimate for the occurrence rate of super-Earths in the habitable zones (HZs) of M dwarfs of 0.41 +0.54 −0.13 planets per star.
M dwarfs are the most abundant type of stars in the Solar neighbourhood. Therefore, the occurrence rate of planets around these stars will dominate any general estimates of the occurrence rate of planets. For this reason, we re-analyse the radial velocities obtained using the Ultraviolet and Visual Echelle Spectrograph (UVES) at VLT-UT2 of a sample of M dwarfs of Zechmeister et al. (2009) using posterior sampling techniques in our Bayesian search for planetary signals. We also extract HARPS radial velocities for these stars from the publicly available spectra in the European Southern Observatory (ESO) archive and analyse the combined UVES and HARPS velocities. The methods are presented in Section 2 in detail and we show the results based on combined HARPS and UVES data in Section 3. We present the statistics of the new planet candidates we detect and compare the obtained occurrence rates to other planet surveys targeting M dwarfs in Section 4, describe some of the interesting new planetary systems and the evidence in favour of their existence in greater detail in Section 5, and discuss the results in Section 6.

STATISTICAL METHODS AND BENCHMARK MODEL
We analyse the radial velocity data sets by using posterior sampling algorithms and estimations of Bayesian evidences for models with k = 0, 1, ... Keplerian signals. Throughout the analyses we apply a fully Bayesian data analysis framework as discussed and applied in several astronomy papers over the recent years (e.g. Ford 2005Ford , 2006Trotta 2007;Feroz et al. 2011;Gregory 2011;Loredo et al. 2012;Tuomi 2012;Tuomi et al. 2013a,b). In particular, we apply the adaptive Metropolis posterior sampling algorithm of Haario et al. (2001) that can be applied readily for analyses of radial velocity data (e.g. Tuomi 2012; Tuomi et al. 2013a,b). We estimate the Bayesian evidences in favour of a given number of signals (i.e. in favour of a model M k that contains k signals) by using the estimate of Newton & Raftery (1994) based on statistical samples drawn from both the prior and the posterior densities, and report parameter estimates using the maximum a posteriori (MAP) estimates and 99% Bayesian credibility intervals (BCS). We note that the acronym BCS stands for Bayesian credibility set, which is represented by an interval in a single dimension when the posterior density does not have multiple significant modes. This set with a probability threshold of δ is a set D δ defined for a posterior density π(θ|m) as where Ω is the parameter space of the parameter vector θ, π(θ|m) is a posterior probability density function given measurements m, ∂D δ represents the edge of the set D δ , and c is some positive constant. Formal definition and discussion can be found in textbooks of Bayesian statistics, e.g. Berger (1980) and Kaipio & Somersalo (2005). Our benchmark statistical model contains linear acceleration and correlation components. In the context of this model, we describe the radial velocity measurement (m i,l ) made at epoch ti with instrument l as in Tuomi et al. (2013a,b) and write it as where γ l andγ are free parameters of the model representing the reference velocity of the lth instrument and the linear acceleration, respectively, and ǫi,j is a Gaussian random variable with a zero mean and variance σ 2 i + σ 2 l describing the amount of Gaussian white noise in the ith measurement of the instrument l. Parameter(s) φ l represents the correlation between the deviations of the ith and i − 1th measurement from the mean, i.e. the dependence of the ith measurement on the deviation of the i − 1th measurement from the mean because there can be no causal relationship the other way around 2 . The parameter vector θ of a "baseline" model without Keplerian signals is then θ = (γ l ,γ, σ l , φ l ) for one instument. In our notation, function f k represents the superposition of k Keplerian signals.
In Eq.
(2), parameter α corresponds to the timescale of the exponential smoothing in the moving average (MA) component . We chose this parameter such that correlations on the time-scale of few days were taken into account because correlations in this timescale are known to occur in radial velocities (Baluev 2013;Tuomi et al. 2013b) but that measurements sufficiently far in time from one another, i.e. more than a dozen days or so, are unlikely to be correlated. Therefore, we chose α = 0.01 hours −1 -a value that was supported by the largest UVES datasets (GJ 551 and GJ 699). This value was also found to describe the data sets well because changing the value resulted in, at most, equally good performance of the statistical model for the data sets with large number of measurements. It was not necessary to treat this parameter as a free parameter of the statistical model because in most of the cases, it could not be constrained at all due to a small number of data points and in such cases the principle of parsimony should be applied to decrease the number of free parameters and e.g. fix the time-scale parameter as we have done in the current work. The value α = 0.01 hours −1 corresponds to a correlation time-scale of roughly 4 days. It is possible that this choice for the time-scale affects the analysis results of some data sets that in fact have correlations on a much shorter (longer) time-scale of few hours (dozen days) instead of days. However, we consider this to be rather unlikely because of the low number of measurements in most data sets. We also note that there are thus 3j + 5k + 1 free parameters in our model if j is the number of instruments that have been used to obtain the data. While our benchmark model in Eq.
(2) cannot be expected to be a perfect description of the radial velocities (e.g. Tuomi & Jenkins 2012;Tuomi et al. 2013b), we still estimate that it contains most of the important features of a good radial velocity model. Particularly, it contains the linear acceleration that could be present due to a previously unknown long-period substellar companion. Even without evidence of such a companion (i.e. of such a trend), we include this linear acceleration in the model to construct a standard model that can be applied to all the data sets we analyse. If there is no linear acceleration in a given data set, its inclusion in the model makes the model overparameterised because dropping the corresponding term from the formulation would provide a better model due to the principle of parsimony. However, we are willing to risk such overparameterisation to be able to use the same general bench-mark model for all the data sets and to obtain as trustworthy results as possible. The same principle applies to the parameter φ l that quantifies the amount of intrinsic correlation in the data obtained using the lth instrument (e.g. Tuomi & Jenkins 2012;Tuomi et al. 2013b). We use this MA component in the model even if we cannot show that it is significantly different from zero, which makes the model more complicated but enables us to see consistently whether the excess noise in the data has a significant redcolour component.
We note that even if there is no evidence for neither acceleration nor correlation in the sense that the parameterṡ γ and φ l have posterior densities that are consistent with zero, including the corresponding terms in the model makes the results more robust because the uncertainties of these terms can be taken into account directly. Furthermore, this ensures that none of the signals we detect are spuriously caused by the combination of these two factors.

Prior choice
Because prior densities of the model parameters are an integral part of Bayesian analyses, we discuss our prior choices briefly. Throughout the analyses, we use prior probability densities as described in Tuomi (2012) with one small but possibly significant exception. In comparison to the choice of Tuomi (2012), we follow  and use more restrictive eccentricity priors in our analyses by adopting π(e) ∝ N (0, 0.1 2 ). In addition to being a more approppriate functional form for an eccentricity prior than a uniform one (see also Kipping 2013, for comparison), we make this choice because the data sets we analyse have already been analysed by Zechmeister et al. (2009) who did not report any planetary signals in the UVES data sets. Therefore, if there are any significant signals in these data sets, they are likely to be deeply embedded in the noise in the sense that their amplitudes are comparable to the noise levels, which results in overestimated orbital eccentricities (Zakamska et al. 2011). Yet, because low-mass planets are mostly found on close-circular orbits , we prefer a slight underestimation of orbital eccentricities over their overestimation. We note that this prior choice is still much more conservative than the commonly made decision to fix eccentricities to zero a priori that correspond to choosing a delta-function prior such that π(e) = δ(e). For further justification of this prior choice, we refer to Appendix A in Anglada-Escudé et al. (2013).
There is also the possibility that our eccentricity prior decreases the significance of a signal that is actually caused by e.g. a planet on a moderately eccentric orbit. However, we consider that to be unlikely because in the current sample, any planetary signals we detect have amplitudes that are comparable to the measurement noise. This means that the corresponding eccentricities are ill-constrained from below and that using a prior density as described above, or indeed the choice proposed by Kipping (2013), is thus unikely to affect the results significantly even if some of the stars in the sample have low-mass planets on eccentric orbits.

Posterior samplings
While the adaptive Metropolis algorithm is an efficient tool in drawing statistically representative samples in cases of unimodal posteriors with little non-linear correlations between the parameters, it is not necessarily well-suited for such samplings of multimodal posteriors that are typical, especially with respect to the period parameters of the signals, when searching for periodic signals of planets. Therefore, when searching for a kth signal we simply divide the period space of this signal into parts that only contain one significant maximum that would, because of its significance, effectively prevent the chain from "jumping" between the different modes. This enables us to use the adaptive Metropolis algorithm. These different parts can then be sampled independently and treated as different (a priori ) models to assess the relative significances of the corresponding maxima. However, in practice, such cases were rare and we only applied such divisions of the period space to the data of two targets.
When estimating the significances of the solutions, i.e. that the Markov chains were sufficiently close to convergence, we used the Gelman-Rubin statistics (Gelman & Rubin 1992) as also described in Ford (2006). In particular, we required that the test statistics R(θ) that approaches unity from above as the Markov chains approach convergence was below 1.1 based on at least four chains to state that the chains were sufficiently close to convergence for all parameters θ. This corresponds to a situation where the variance of the parameters is lower between chains than within them indicating that all the chains have identified the same stationary distribution.

Signal detection criteria
Throughout the analyses, we used the signal detection criteria of Tuomi (2012). These criteria have recently been applied in e.g. Anglada-Escudé et al. (2013), Tuomi et al. (2013a), Tuomi et al. (2013b), , and they appear to be trustworthy in the sense that they are not particularly prone to false positives (Tuomi & Jenkins 2012), and enable the detections of signals that cannot be found by using more traditional detection criteria based on false alarm probabilities (FAPs) in the power spectrum of the model residuals (Anglada-Escudé & Tuomi 2012;Anglada-Escudé et al. 2013;Tuomi et al. 2013b). As an example, we refer to , who independently discovered the same three planet candidates around GJ 163 as Bonfils et al. (2013b) with only ∼ 35% of the data. This indicates that the Bayesian signal detection criteria indeed are very sensitive and robust ones in detecting weak signals of low-mass planets in radial velocity data.
The first criterion is that a model with k + 1 signals, denoted as M k+1 , has a posterior probability that is at least s times greater than the corresponding probability of a model with only k signals, i.e. P (M k+1 |m) sP (M k |m). If this criterion is not satisfied, the k + 1th signal has a considerable probability of being produced by noise instead of being a genuine periodicity in the data. When analysing the data with k = 0, 1, ... it might happen that P (M k |m) < sP (M k−1 |m) but that P (M k+1 |m) sP (M k |m). This means that it cannot be said there are significantly k signals, and if the model M k+1 was not tested, the conclusion would be that there are only k − 1 signals in the data. Therefore, if we observe putative probability maxima in the parameter space for model M k that do not satisfy the detection criterion, we typically analyse the data with models M k+1 and M k+2 as well to see if the superposition of more than k signals makes the detection of only k signals impossible. We note that the probability threshold s has to be chosen subjectively. A typical choice would be s = 150 because it has been interpreted as corresponding to strong evidence and recommended by e.g. Kass & Raftery (1995) based on the arguments of Jeffreys (1961), although it has been argued that a more conservative threshold of 1000 should be chosen (Evett 1991). This threshold of s = 150 has been applied in analyses of radial velocity data succesfully (e.g. Feroz et al. 2011;Gregory 2011;Tuomi 2012;Tuomi et al. 2013a,b). However, this choice might result in overinterpretation of the significance when the statistical model used to describe the data does not represent the data well. For the purpose of population studies such as the current work, one should be more cautious than usual in accepting new planet candidates to avoid the contamination of occurrence rate estimates by false positives. Therefore, we remain cautious about the interpretation of signals that are detected with low thresholds, and use a more conservative threshold of s = 10 4 as a requirement for a planet candidate.
The second criterion is related to the first one in the sense that it states that the radial velocity amplitude has to be statistically significantly (with e.g. 99% level) greater than zero. If this was not the case, there would be a considerable probability that the amplitude of the signal was actually negligible and that the signal did not thus exist in reality. The third criterion states that the period of a signal has to be well-constrained from below and above to consider it a genuine periodicity. Again, the justification of this criterion is simple because if the period parameter was not well constrained, it could not be claimed that the corresponding variations in the data were indeed of periodic nature. We note that in practice, based on the analyses in the present work, a signal whose significance exceeds s = 10 4 is always well constrained in the parameter space. Conversely, if a given signal is constrained in the period and amplitude space, its significance always exceeds some s > 1, though not necessarily s = 150 or s = 10 4 . This means that the detection criteria are complementary and robust in detecting low-amplitude signals.
To assess whether any given signal observed in the combined data set is supported by both data sets, we examine the residuals of each model using common weighted root-mean-square (RMS) statistics. If these statistics are decreased for both data sets when adding a signal to the model, we say that both data sets support the existence of the corresponding signal even when the signal cannot be detected in both data sets independently. This choice also serves as a test of whether a given weak signal could be a spurious one generated by noise and the correspondingly insufficiently accurate model. For a given signal to be a genuine one of stellar origin (as opposed to a spurious signal of instrumental origin) -possibly caused by the Doppler variations induced by an orbiting planet -it should result in a decrease in the RMS estimates of both data sets, although this is also subject to chance and does not necessarily hold when there are only few measurements.
Finally, to make the detections as robust as possible, we set the following criteria for signals to be able to call them planet candidates. First, we require that such signals are detected with a choice of s = 10 4 and can be concluded to be present in the corresponding combined data set(s) very significantly. If counterparts of these signals cannot be found in the activity indices calculated from the UVES and HARPS spectra, they are considered planet candidates and we refer to them according to the standard nomenclature by assigning letters b, c, ... to them. If, however, (1) the corresponding signals are detected only with a significance of s < 10 4 , (2) if there is data available from only one spectrograph such that the existence of a signal is not supported by an independent data set because such a set is not available, or (3) if a signal is present only in one of the data sets (in the sense that the RMS of the other does not decrease when adding the signal to the model) that dominates the posterior density given the combined data because of the low number or poor quality of the measurements in the other set, we do not call them candidate planets.
We note that if the posterior density consists of several other (local) maxima in the parameter space at several different periods with reasonably high probabilities (e.g. with n local maxima at θ i L , i = 1, ..., n, that have π(θ i L |m) 0.01π(θMAP|m) for all i, where θMAP is the MAP estimate), it is possible that there are also several significant solutions in the period space that satisfy the detection criteria. Such a situation is hard to interpret reliably, but can arise for two main reasons. First, all the significant periodicities can correspond to genuine periodic signals in the velocity dataseveral planetary signals in the data may have similar amplitudes and be detected roughly equally confidently (e.g. GJ 667C whose Doppler data was found to contain six signals of roughly equal amplitude; Anglada-Escudé et al. 2013). Alternatively, such a situation might be representative of poor noise modelling where periodic/quasiperiodic features and/or correlations in the noise are falsely interpreted as genuine signals. In such cases, we analyse the data by increasing k further but report the solution corresponding to the global maximum if there is no strong evidence in favour of additional signals in the sense of our detection criteria.

Search for significant periodicities
We performed the searches for significant periodicities in the data sets in several stages. First, we obtained a sample from the posterior density of the model without any Keplerian signals. In the second step, we sampled the posterior density of the model with k = 1 by using tempered chains such that the posterior density was raised to a power of β ∈ (0, 1). In particular, we used β = 0.3 -0.8 to find the most promising areas in the period space because such a choice of β enabled a rapid search of the whole period space [T0, T obs ], where T obs is the data baseline and we selected T0 = 1 day. However, if there were clear indications of probability maxima in excess of a period equal to T obs , we increased the upper limit to 2T obs . We did not use the adaptive Metropolis algorithm (Haario et al. 2001) but the standard Metropolis-Hastings version (Metropolis et al. 1953;Hastings 1970) in these periodicity searches to prevent the proposal density from adapting to the possibly very narrow global maximum in the posterior density because the purpose of these samplings was to identify the positions of the most significant probability maxima in the parameter (period) space and not to enable the chain to adapt to one of them.
For most data sets, we started by a test sampling with a "cold" chain, i.e. a chain with β = 1. When the number of measurements was lower than ∼ 15, these cold samplings were sufficient in exploring the whole period-space because typically there were no significant maxima in the parameter space that would have prevented the chain from exploring the whole period space rapidly. When such cold samplings were not possible due to an abundance of reasonably high maxima in the parameter space and thus poor mixing of the Markov chains in the period space, we used tempered chains by setting β = 0.3 to 0.5, depending on the number of measurements. The lowest values of 0.3 were used for the largest data sets of GJ 699 (NUVES = 226), GJ 551 (NUVES = 229), and GJ 433 (NUVES = 166). These samplings resulted in estimates for the positions of the most significant maxima in the period space. However, it must be noted that using tempered samplings corresponds actually to using different prior and likelihood models in the analysis such that π β (θ) and l β (m|θ) are used as a prior density and likelihood function instead of π(θ) and l(m|θ), respectively. Therefore, the shape of the posterior density, as estimated based on tempered samplings, does not necessarily correspond to the actual posterior density. However, the positions of the maxima are unchanged.
Given rough estimates for the positions of the probability maxima in the period space based on tempered samplings, we started several cold samplings with initial values in the vicinity of the observed maxima to enable fast convergence. If one (or some) of these maxima corresponded to a significant periodic signal(s) in the sense of the detection criteria we described above (Section 2.3), we continued by increasing k and by performing tempered samplings of the parameter space of a model with one more Keplerian signal. These samplings were performed by using cooler chains because after the strongest signal was accounted for by the model, samplings with β = 0.6 -0.8 were typically found to visit all areas of the period space of the additional signal with ∼ few 10 6 chain members.
Finally, we obtained samples from the posterior density of a model with k+1 Keplerian signals when the corresponding data was found to contain significant evidence in favour of only k signals, i.e. that only k signals satisfied the detection criteria. These samplings were performed to estimate which additional signals could be allowed by the data and which ones could thus be ruled out (e.g. Tuomi 2012).
We note that some data sets have only two UVES measurements and/or less than two HARPS measurements. Therefore, because we used a linear trend in the model, we did not analyse UVES data sets with less than three mesurements. Furthermore, we did not analyse the combined set when the number of measurements in either UVES or HARPS set is less than two. This means that we might analyse an individual data set with three measurements and a combined one with four (two UVES and two HARPS) measurements. Even in such extreme cases, when the number of free parameters in the statistical model exceeds the number of measurements, we expect to be able to obtain meaningful results because we use informative prior densities and because the Bayesian statistical techniques we use take the principle of parsimony, and thus the possible effects of such overparameterisation, into account and do not rely on assumptions regarding the number of parameters or measurements.

Example: GJ 229
To demonstrate the practicality of our periodicity search technique, we show the analysis results in detail for the UVES velocities of GJ 229 because it is a typical target in the sample with a reasonably large number of measurements (NUVES = 73) that show evidence in favour of periodic variations. We started by analysing the data with the benchmark model with k = 0.
We found that the UVES data of GJ 229 contained a significant amount of correlation with φUVES = 0.85 [0.45, 1.00]; evicence for a linear trend withγ = 1.26 [0.39, 2.17] ms −1 year −1 , small part of which is caused by secular acceleration 3 of 0.070 ms −1 year −1 (Zechmeister et al. 2009); and excess Gaussian white noise with σUVES = 2.51 [1.29, 3.74] ms −1 , which is rather low and implies that the star expresses very low levels of velocity variability and/or the UVES instrument uncertainties in the Zechmeister et al. (2009) velocities might be overestimated. Removing the MAP trend and correlations, we plotted the resulting residuals in Fig.  1. These residuals can indeed be described as "flat", which indicates that there cannot be periodic variations in these velocities with amplitudes greater than roughly 10 ms −1 . The RMS of the residuals is 3.61 ms −1 , which is considerably lower than that of the original velocities of 5.26 ms −1 .
We performed a tempered search for periodicities, and obtained a sample from the posterior density corresponding to a choice of β = 0.3. We plotted an example of a corresponding sampling in Fig. 2. As can be seen, the chain identifies a global maximum in the period space of the one-Keplerian model at a period of 2.8 days (red arrow). The chain also visits the whole period space between 1 day and Figure 2. Log-posterior density as a function of the log-period parameter of a Keplerian signal from a tempered sampling (β = 0.3) of the GJ 229 UVES data. The horizontal lines indicate the 10% (dotted), 1% (dashed), and 0.1% (solid) probability thresholds with respect to the MAP value that is denoted by using the red arrow.
T obs = 2325 days for GJ 229. The chain in Fig. 2 also shows that there are additional local maxima in the period space exceeding the 10%, 1% and/or 0.1% probability levels (dotted, dashed, and solid horisontal lines, respectively) with respect to the global maximum at periods of roughly 1.5, 10, 200, and 450 days. The existence of such multiple maxima in the scaled posterior suggests that none of them can be confidently considered a solution and thus a periodic signal that is reliably detected in the data. Our samplings also suggest that there are likely no other considerable probability maxima in the period space, or that if they exist, they are so narrow that the chains are unlikely to visit them. We performed such samplings several times and obtained consistent results with no indications of additional maxima with similar posterior values.
Because of the candidate periodicities, we started cold chains with initial periods in the vicinity of the highest maxima ( Fig. 2) and were able to identify two periodicities that satisfied the detection criteria discussed above. These periodicities were found at periods of 1.5 and 206 days with amplitudes of 5.3 and 4.3 ms −1 days, respectively, but it cannot be concluded that we have identified such signals confidently in the data because the posterior density has multiple maxima comparable to the posterior at the MAP estimate. Therefore, we can only conclude that if there are periodic signals in the UVES data, they are likely present at the periods corresponding to the maxima in the posterior density (Fig. 2).
We performed samplings of the UVES data with a model containing two Keplerian signals, but could not find significant two-Keplerian solutions to the data.
It should be noted that the samples from the (scaled) posterior densities in Fig. 2 can be roughly interpreted in a similar manner as e.g. common Lomb-Scargle (Lomb 1976;Scargle 1982) periodograms. This is because the plotted density represents the relative significances of the different periods and can thus be broadly interpreted according to the probabilities presented in the vertical axis showing the logposterior values. However, the relative probabilities of the corresponding periodicities, and indeed whether they are significant enough to satisfy the detection criteria, can only be assessed by additional samplings and by seeing whether the periods can be constrained from above and below.
One last test for a significance of a signal is whether it is present in two independent data sets or not. According to our results, the putative signals in the UVES data at periods of 1.5 and 206 days fail this test as the HARPS data effectively rules them out. This can be seen by looking at the corresponding log-posterior density as a function of the period of a signal in the combined data (Fig. 3, top panel). Accordingly, there is a strong and isolated global maximum at a period of roughly 470 days that satisfies all the detection criteria discussed above. Because this signal is clearly present as a local maximum in the UVES data alone, and because there are no comparable local maxima in the period space given the combined data, we conclude that there is a significant periodic signal in the combined HARPS and UVES velocities of GJ 229 at a period of 471 [459,493] days with an amplitude of 3.83 [2.15, 5.57] ms −1 that corresponds to a new and previously unknown planet candidate orbiting the star with a minimum mass of 32 [16,49] M⊕. This signal can be called a candidate planet because its existence is supported by both data sets according to our requirements and because it is detected very confidently in the combined data. Furthermore, as we discuss in Section 3.1, the signal does not correspond to variations in the activity data of GJ 229. To demonstrate its significance, we have plotted the phase-folded signal in Fig. 4. We note that there is no evidence in favour of a second signal in the combined data, as can also be seen in the bottom panel of Fig. 3 that does not have strong and isolated maxima in the period space.

Detection probability and planet occurrence
To obtain estimates of the underlying occurrence rate of planets in our sample, it is necessary to estimate the detection bias of the current data caused by the fact that radial velocity data sets are more sensitive to greater planetary masses and shorter orbital periods. It is also necessary to account for the quantity and quality of each data set, as well as their respective baselines. For this purpose, we calculated detection probabilities for the combined UVES and HARPS data sets of the sample by using samplings of the parameter space as shown in e.g. the bottom panel of Fig.  3.
First, we estimate that in any given data set (the ith data set), the observed number of signals in a given period interval ∆P and minimum mass interval ∆M , that we denote as ∆P,M for short, can be written as f obs,i (∆P,M) = focc,i(∆P,M)pi(∆P,M), where focc,i is the number of planets orbiting the ith star in the sample and pi is the detectability function that indicates whether a planet with parameters in the interval ∆P,M can be detected in the data set mi. However, while f obs,i is easy to obtain as it is directly the number of planets we detect in the respective mass-period interval in data set mi, pi is more difficult to calculate.
To estimate pi, we use the posterior samplings of a model with k + 1 Keplerian signals when there are only k signals in the data set. The sample drawn by using the posterior sampling enables us to reconstruct the areas of the parameter space where the parameters describing the hypothetical k + 1th signal visited. Because the signal cannot be detected in the sense that its amplitude and period could  be constrained, the chain visits all areas in the period space allowed by the chosen range of the period (from one day to data baseline). Moreover, the chain visits, given sufficiently good mixing properties such that the chain visits all relevant areas in the parameter space frequently, all amplitudes at all periods that are allowed by the data, i.e. amplitudes that are so low that they cannot be ruled out by the data because the likelihood function has reasonably high values. The chain does not visit amplitudes in excess of some limiting amplitude at each period because they would correspond to such low likelihoods that planetary signals with such amplitudes are effectively ruled out. Therefore, the sampling The white area corresponds to parameter values in the massperiod space where signals could be detected in the data (detection probability of unity) whereas the black area shows a corresponding mass-period space where signals cannot be detected (negligible detection probability). The candidate GJ 229 b is shown as a red circled dot that is in the former area because it was detected.
yields the areas of parameter space where there could be signals. The rest of the parameter space is thus the area where there are (very likely) no additional signals based on the data.
We demonstrate this by using GJ 229 as an example. It has one candidate and thus we use a model with two of them to estimate the detectability function. In Fig. 5, the black area shows the subset of the mass-period space where the mass and period parameters of the second signal visited during the Markov chain samplings. This area thus corresponds to signals that could exist in the data but that cannot be detected because we did not find a second signal in the GJ 229 data. The white area corresponds to massperiod space where the chains did not visit and thus rule out planetary signals because the likelihood function is so low in these areas that the chances of a planet existing in the white area can be approximated to be negligible. The planet candidate GJ 229 b is clearly above the threshold as it should because it was detected. We note that while this planet candidate has parameters close to the threshold, it is still several M⊕ heavier than the minimum-mass planet that could be detected at that period. In fact, if it was much more above the threshold, it is likely that its existence would already have been reported based on UVES or HARPS data sets alone.
We thus estimate that planet candidates could have been detected in the areas where the Markov chains did not visit (white areas in Fig. 5). This is the case because they correspond to the complement of the area where planets could not be detected. In reality, however, the threshold is not so strict because the detection probability is a (continuous) function of the parameters. We do not attempt to estimate this function more accurately in this work and use the "step-function" for each data set as shown in Fig. 5 for GJ 229 as a first-order approximation.
Assuming that there are some areas in the parameter space where the chains did not visit due to e.g. poor sampling even though they should have because these areas have sufficiently high likelihoods, we would effectively overestimate the detection probabilities and thus underestimate the occurrence rates. Thich means that our occurrence rate estimates are lower limits, although unlikely to be considerably lower to bias the results significantly given the rather large uncertainties due to the low number of planets in the sample. We approximate the detectability by setting the functionpi = 1−pi equal to unity in the areas where the Markov chains did visit, and equal to zero where they did not.
We express the detected frequency of planets for the whole sample of stars, with data sets m1, ..., mN , by summing the number of observed signals f obs,i of all the data sets and by assuming that the occurrence rate is common for all stars in the sample such that focc = focc,i for all i (in units of planets per star). Thus we obtain which implies a simple way of calculating the occurrence rate focc for the whole sample. In this equation, the term in square brackets on the right hand side is, when divided by N , the detection probability function of our sample that approximates the probability of being able to detect planets in the interval ∆P,M. Thanks to our samplings, we could estimate this function rather accurately -typically by using a 100×100 grid in the log-parameter space ranging from a minimum period of 1 day to a maximum of 10 4 days and from a minimum mass of 1 M⊕ to a maximum of 100 M⊕. Such a fine grid was not a practical choice for estimating the occurrence rates, because to obtain any meaningful estimates at all, it would be desirable to have at least one planetary signal in most grid points. For this reason, when calculating occurrence rates, we divided the interval into a 4×4 grid. When estimating the uncertainties for the occurrence rates in a given interval ∆P,M, we use rather conservatively the lowest and highest detection probabilities (obtained by using the finer grid) in this interval to calculate lower and upper limits, respectively.

THE SAMPLE AND KEPLERIAN SIGNALS
The sample of M dwarfs for which the UVES velocities were published in Zechmeister et al. (2009) contains a collection of nearby stars including the two nearest M dwarfs, namely GJ 551 and GJ 699 that are the nearest and fourth nearest stars to the Sun, respectively. We have listed the stars in this sample in Table 1 together with their estimated physical properties. While the parallaxes are obtained from Hipparcos (van Leeuwen 2007) and the mass estimates from Zechmeister et al. (2009), we estimated the effective temperatures and luminosities by using the empirical relations of Casagrande et al. (2008) and Boyajian et al. (2012). Based on the estimated luminosities and effective temperatures, we have also calculated the approximate inner and outer edges of the stellar habitable zones according to the equations of Kopparapu et al. (2013) and listed them in Table 1.
We note that the stars in the sample have also been extensively searched for brown dwarf companions. For instance, Dieterich et al. (2012) reported that 0.0 +3.5 −0.0 % occurrence rate of L and T companions in separations between 10-70 AU. Furthermore, the sample stars show no evidence for warm (e.g. Avenhaus et al. 2012) or cold (e.g. Lestrade et al. 2009) circumstellar material.
We obtained the HARPS-TERRA (Template Enhanced Radial velocity Re-analysis Application) velocities from the publicly available spectra 4 by using the data processing algorithms of Anglada-Escudé & Butler (2012) . While an abundance of such velocities could not be obtained for every star because a large fraction of the stars in our sample have not been primary targets of the HARPS-GTO survey (Bonfils et al. 2013a), there were still several stars for which the available HARPS data could be readily expected to provide better constraints for the possible planet candidates orbiting them due to its high precision, or help disputing the existence of putative signals in the UVES data as false positives. The numbers of HARPS measurements and the corresponding data baselines are listed in Table 2 and we have tabulated the corresponding HARPS-TERRA velocities in the Appendix B.
We present the results of our comparisons of models containing k = 0, 1, 2 Keplerian signals in Table 2 by presenting the numbers of favoured signals in the combined data sets together with the numbers of measurements and the baselines of the UVES and HARPS data sets. According to these results, there are 10 significant signals in the combined data that satisfy our detection criteria and are supported by both data sets according to our requirements.
When publishing the analysis results of the same UVES datasets, Zechmeister et al. (2009) did not report detections of planet candidates. According to our periodogram analyses with the standard Lomb-Scargle periodogram (Lomb 1976;Scargle 1982) of the UVES data, this conclusion is indeed justified because, excluding the massive stellar or substellar companions 5 around GJ 477, GJ 1046, GJ 3020, and GJ 3916; and excluding GJ 551 and GJ 699 for which Zechmeister et al. (2009) reported signals caused by data sampling and activity, none of the UVES data sets had significant powers in their Lomb-Scargle periodograms corresponding to low-mass planetary companions. This inability to detect the signals of the planet candidates we report is not surprising because all these signals have very low amplitudes whose detection is difficult without a large number of measurements and high precision.
Interestingly, unlike Zechmeister et al. (2009), we could not find the 44-day signal they reported in the velocities of the GJ 699. We believe the reason is that this signal is likely  caused by activity-related phenomena (Zechmeister et al. 2009) and is therefore quasiperiodic and/or time-dependent and explained rather well by correlations in the UVES measurements. Indeed, the estimate for the parameter φUVES was found to be 0.81 [0.62, 0.99], which implies a considerable amount of correlation in the UVES data that could -when coupled with data sampling -give rise to significant but spurious powers in periodogram when not accounted for.
In addition to the signals we find, five targets in our sample show linear acceleration that is not consistent with zero suggesting the existence of yet unknown long-period substellar companions (Table 2). Table 2. Properties of the UVES and HARPS data sets in terms of numbers of measurements (N UVES , N HARPS ) and data baselines (∆ UVES , ∆ HARPS ). Log-Bayesian evidence ratios, or Bayes factors (ln B i+1,i = ln P (m|M i+1 ) − ln P (m|M i )), are shown when there is evidence in favour of at least one Keplerian signal according to our detection criteria. GJ 190 and GJ 263 are not in the table because they only had two velocity measurements. The four stars with evidence of massive companions in the radial velocity data are also not shown.

Analysis of activity indicators
To determine whether the signals we observe in the UVES radial velocities are caused by Doppler fingerprints of planetary and/or substellar companions or periodic or aperiodic and/or quasiperiodic phenomena related to the stellar activities, we analysed the line bisector (BIS) spans as obtained from the UVES spectra available in the ESO archive. BIS values indicate the line asymmetry and can be used as a signature of variations caused by stellar activity (e.g. Queloz et al. 2001;Boisse et al. 2011). For most stars for which we observed radial velocity signals, the correlation coefficients between the velocities and BIS values were between -0.1 and 0.1, which implies no significant correlation. However, GJ 842 showed positive correlation coefficient of 0.43 between the radial velocities and BIS values that ap-pears significant although the corresponding data set only contained 16 data points 6 . We also obtained the BIS values from the HARPS spectra and tested whether there were correlations between these values and the radial velocities. None of the data sets showed significant correlations between the BIS values and velocities, which indicates that the variations in the velocity data are unlikely to have been induced by activity-related phenomena of the stellar surface such as starspots and/or active regions. The typical correlation coeficients were found to be between -0.1 and 0.1, which is consistent with the inter- 6 We could only obtain 16 spectra for GJ 842 because the ESO archive did not contain calibration frames although there are 17 radial velocities in the Zechmeister et al. (2009) data set (see also pretation that the signals we observe are genuine Doppler signatures of planetary nature.

UVES false positives
The UVES data contained some signals that were not supported by the HARPS velocities. This means that our benchmark model might not be optimal in the sense that features in the UVES measurement noise, or indeed biases, instrument-related variations, stability problems, etc. can produce variations that are interpreted as Keplerian signals. Furthermore, according to our analyses, even when we do detect a signal in the UVES data, the period space is typically littered by several other almost equally high local maxima that should be interpreted as alternative solutions before they can be ruled out. Such a multimodality occurs e.g. in the UVES data of GJ 160.2, GJ 229 (see Fig. 2), GJ 377, and GJ 27.1.
We identified only two false positives in the UVES data that were essentially ruled out by the HARPS velocities. The first one of these is a 15-day periodicity in the 14 UVES velocity measurements of GJ 377. This signal increases the model probability with respect to the zero-Keplerian model by a factor of 700, which does not exceed our detection threshold of 10 4 . In this case, the HARPS data did not confirm the existence of the signal but ruled it out rather confidently. This weak evidence for a signal in the UVES data was also suspect because the period parameter of the signal had local probability maxima around 2, 5, 50-200, and 700 days that exceeded the 0.1% or even 1.0% probability thresholds of the global maximum. We interpret this result as hints of a signal in the UVES data whose exact period or amplitude cannot be determined due to the low number of measurements. Instead, in the combined dataset the global maximum of the period parameter was found at a period of roughly 80 days but it did not satisfy the detection criteria.
The other example is provided by the solution to the UVES data of GJ 357 that consists of two periodicities at 5 and 26 days, respectively. The significances of these signals decrease below the detection threshold with only five HARPS velocities. This means that -even though there is not enough HARPS data to constrain any signals in the combined data set -some of the local maxima in the UVES data can be ruled out 7 , but others might actually correspond to signals whose existence could be verified when future highprecision data becomes available.
We note that the iodine cell used in the UVES to obtain precise reference lines to the spectra can give rise to differences with respect to the ThAr method used in the HARPS. This means that some of the spurious signals and/or spurious local maxima that are unrelated to sampling in the UVES data can arise from small-scale instability in the UVES instrument. The same argument applies to the HARPS velocities as well, although HARPS has much greater long-term stability and precision. However, at the moment it is not possible to distinguish whether the signals with low significances that are interpreted as false positives in the UVES data because they do not correspond to global maxima in the combined data are caused by such instability, periodic or quasiperiodic features in the noise, periodicities related to the stellar activity and magnetic phenomena, and/or genuine Doppler signatures caused by low-mass planets. Ideally, signals arising from instrumentation can be ruled out by obtaining support for them from two data sets, and those caused by stellar surface phenomena can in principle be ruled out by showing that they or their primary harmonics do not have any couterparts in the activity-indices of the two sets. However, any remaining signals can still be caused by other unknown sources of periodic variation.

PLANETS AROUND M DWARFS
We have plotted the known planets 8 (light blue dots), known planets around M dwarfs (blue circled dots), the new candidates we report (red circled dots) in Fig. 6. This figure also shows the estimated detection probability of the combined UVES and HARPS data sets as functions of minimum mass and orbital period for the whole sample (Section 2.6). The most remarkable feature in this plot appears to be that the new candidates are concentrated around minimum masses of 10 M⊕ and periods from few to few dozen days. Considering that such planets can be detected in this sample with rather low probabilities of roughly 10-30%, the abundance of such companions appears strikingly high. Yet, the results from analyses of Kepler's data appear to imply that this is a real feature in general (Howard et al. 2012), and applies to M dwarfs in particular, as shown in Dressing & Charbonneau (2013). They observed an abundance of transiting planets with the same period range and radii of up to 3 R⊕, which could correspond to the same population of planets of which we only observe the planets with the highest masses.
The detection probability of a given combination of minimum mass and period in Fig. 6 shows the different areas in the mass-period space where planets can be spotted easily and where it is very difficult (or unlikely) given the current sample. The highest gradient in this probability occurs along the line increasing from 4 M⊕ at a period of one day to 100 M⊕ at 1000 days. In our sample, only one candidate can be found above this line, even though planets in this region would be the easiest ones to detect. We note that the rather artificial-looking vertical line at roughly 3000 days in the top right corner of the Fig. 6 is in fact a threshold arising from the fact that we limited the period search to periods at most the baselines of the data sets. Furthermore, there is a weakly distinguishable feature that shows decreased detection probabilities for periods around 365-day period that demonstrates how the existence of planetary signals at one year period cannot be ruled out generally as easily as slightly shorter and longer ones due to poor phase-coverage caused by data samplings.
Based on the analyses of the combined UVES and HARPS radial velocities of the 41 nearby M dwarfs, these data sets contain the signals of eight new exoplanet candidates out of which seven can be classified as super-Earths Figure 6. Planet detection probability as determined by using Eq. (3) in the combined UVES and HARPS data set as functions of orbital period and minimum mass. The various dots represent the known planets orbiting all stars (light blue dots), known planets orbiting M dwarfs (circled blue dots), and planet candidates in our sample (circled red dots). The detection probabilities do not exceed 85% even at the high-mass short-period corner of the plot because there are six data sets where planetary signals could not be detected at all due a combination of low number of measurements and evidence of a massive companion that prevented detections of additional companions due to overparameterisation of the benchmark model. due to their minimum masses that are higher than that of the Earth but still lower, for most candidates considerably so, than one Neptune-mass. In addition to these super-Earths, we report a detection of a more massive sub-Saturnian companion orbiting GJ 229, and confirm the existence of the long-period planet orbiting GJ 443 detected by Delfosse et al. (2013). Eight of these candidate planets are thus new and previously unknown. We have listed the obtained orbital parameters together with the inferred minimum masses and semi-major axes of the new planet candidates we report in Table 3. The minimum masses and semimajor axes have been estimated by using the stellar masses as presented by Zechmeister et al. (2009) and by assuming conservatively that their standard deviations are 10% of the estimated values. The last column of Table 3 shows the estimated locations of the planets, i.e. whether they are in the habitable zone , in the cool zone outside the outer edge of the habitable zone, or in the hot zone inside the inner edge of the habitable zone. We have also plotted the estimated log-posterior densities as functions of signal periods, indicative of significant periodicities, in the Appendix (Figs. A1 -A7) together with phase-folded signals and probability distributions for the periods and signal amplitudes that demonstrate that the signals satisfy the detection criteria.
We also broadly verified the existence of the signals corresponding to the planet candidates in Table   3 by using an independent statistical method. We applied the log-likelihood periodograms of Baluev (2009) and Anglada-Escudé et al. (2013) and used the same statistical model as in the Bayesian analyses. The results of these periodogram analyses are shown in Table 4. Six of the signals are clear with FAPs below a threshold of 0.1%. However, GJ 180 c, GJ 422 b, and the two candidates of GJ 682 are only detected with FAPs greater than this threshold. While we suspect that the high FAP of GJ 422 b might be due to priors that were all assumed to be flat when calculating the log-likelihood periodograms and the assumptions behind the significance tests of these periodograms (Baluev 2009), it appears to be clear why GJ 180 c and the candidates GJ 682 b and c were not detected with the periodograms. This is because for GJ 180 c and GJ 682 c, the solution is actually obtained by assuming that there is only one signal when in reality there are two and the k = 1 model is therefore fitted to the superposition of the two signals that causes an apparent decrease to the significance of the second signal. Instead, the signal of GJ 682 b could not even be detected (it was also below the detection threshold with the Bayesian tools because the log-Bayesian evidence ratio corresponding to the detection threshold of 10 4 is 9.21, see Table 2) because the k = 1 model fitted the data much more poorly than a k = 2 model. With the posterior samplings, however, detecting these signals was possible (Figs. A3 and A7).  ); (C) cool planet, i.e. outside the outer edge of the HZ; (H) hot planet, i.e. inside the inner edge of the HZ. d The orbit cannot be constrained from above but there is little doubt about the existence of this signal in the combined data, and it has been detected by using a larger set of HARPS data (Delfosse et al. 2013), thus we list it as a candidate planet. a The first signal in the GJ 682 data could not be detected because a k = 1 model fits the combined data poorly due to a superposition of two signals with almost equal amplitudes.

Potential additional signals
The sample also contained ten additional signals that were well-constrained in period and amplitude but that did not exceed the detection threshold of s = 10 4 or were supported by data from only one instrument. We call these signals SRCs (signals requiring confirmation). They did, however, all exceed a less conservative threshold of s = 150 and it would thus be wrong to simply ignore the existence of these "emerging" signals. We took these signals into account when calculating the detection thresholds but we do not (yet) con-sider them to be candidate planets because we wish to avoid overestimation of the occurrence rates. However, we do calculate the occurrence rates under the assumption that these additional signals are planet candidates but remain cautious when interpreting the corresponding results as some of these signals might be false positives caused by noise, insufficient modelling, stellar activity, and/or instrumental artefacts. The ten additional signals are found in the combined data of GJ 433 (at a period of 36.0 days), GJ 551 (332 and 2200 days), GJ 821 (12.6 days), GJ 842 (190 days), GJ 855 (12.7 and 26.2 days), GJ 1009 (24.5 days), GJ 1100 (34.4 days) and in the UVES data of GJ 891 (30.6 days).

Occurrence rates
When calculating the planet occurrence rates as described in Section 2.6, we obtained some interesting estimates and show them in Table 5. We divided the period space into four bins: between 1 and 10 days, 10 and 100 days, 100 and 1000 days, and 1000 and 10 4 days and the minimum mass space into three bins between 3 and 10 M⊕, 10 and 30 M⊕, and 30 and 100 M⊕ (the last bin for masses above 100 M⊕ is omitted from the table). The resulting occurrence rates show several features that can be considered as representative of the underlying population of planets around M dwarfs. We also show the detection probabilities of each bin based on the whole sample in the bottom of Table 5 for comparison.
According to the estimated occurrence rates in Table 5 (top), the occurrence rate of super-Earths and more massive planets with mp sin i of up to 30 M⊕ increases dramatically for periods between 10 and 100 days. We find that the occurrence rate of planets with masses between 10 and 30 M⊕ is 0.06 +0.11 −0.03 for the sample in this period range. The most dramatic occurrence rate can be found for super-Earths (3 M⊕ < mp sin i 10 M⊕) on orbits between 10 and 100 Table 5. Expected numbers of planets per star based on the sample of stars studied in the current work (top), potential numbers assuming that the ten signals requiring confirmation (SRCs) correspond to planet candidates (middle), and detection probability (DP) of the whole sample in each mass-period bin (bottom).
Planet candidates 1 day < P 10 days 10 days < P 100 days 100 days < P 1000 days 1000 days < P days of 1.02 +1.48 −0.69 per star, which indicates that such planets are very common around M dwarfs in the Solar neighbourhood. Comparison with the results of Bonfils et al. (2013a), which is a similar survey of a larger sample of nearby M dwarfs in the southern sky, is difficult because they did not detect any candidates with masses between 10 and 100 M⊕ and periods in excess of 10 days. They did, however, detect two candidates with masses between 1 and 10 M⊕ and periods between 10 and 100 days, which implied an occurrence rate of such companions of 0.54 +0.50 −0.16 . These estimates are broadly consistent when bearing in mind the different sensitivities of the approaches and the fact that we were able to detect more such companions in the sample than Bonfils et al. (2005) could in their larger sample of 104 M dwarfs.
Second, the occurrence rate of planets more massive than 10 M⊕ on orbits with periods less than 10 days is very low, roughly 0.037 +0.011 −0.006 planets per star for planets less massive than 30 M⊕. In comparison to the results obtained by the HARPS M dwarf survey, Bonfils et al. (2013a) reported the the occurrence rate of planets with periods from 1 to 10 days and masses between 10 and 100 M⊕ to be 0.03 +0.04 −0.01 , which is consistent with our estimate of 0.037 +0.011 −0.006 . This result is also consistent with the observation of Dressing & Charbonneau (2013) that the occurrence rate of planets with radii in excess of 2.0 R⊕ on such orbital periods is roughly 0.05 and is therefore very likely a general feature for M dwarfs. However, comparing our results and the results of Bonfils et al. (2013a) with Kepler occurrence rates based on estimated planetary radii cannot be performed confidently without bias. We find only a slighly higher occurrence rate of 0.06 +0.11 −0.03 for planets with masses below 10 M⊕ because there was only one such candidate in our sample.
These results are very difficult to compare with the results of Dressing & Charbonneau (2013) and Morton & Swift (2013) because the mass-radius relation is far from a well-established one for a range of masses and radii in the super-Earth regime. Yet, Dressing & Charbonneau (2013) found the occurrence rate of planets to increase by a factor of ten when moving from radii interval of 2.8 -5.7 R⊕ down to 1.4 -2.8 R⊕, which is as dramatic increase as we find when moving from masses between 10 and 30 M⊕ down to the interval between 3 and 10 M⊕. While this does not mean that the results are consistent, it suggests that this large change in abundance may have the same origin. Obviously, this depends on whether the low-temperature Kepler sample and the M dwarf sample analysed in the current work are drawn from the same population. Unlike the Kepler sample which comprises of more massive and brighter M dwarfs further out from the galactic plane, samples analysed in the current work and in Bonfils et al. (2013a) are likely to be approximately volume-limited.
Third, our results suggest that planets with masses below 100 M⊕ and periods longer than 100 days might be abundant around nearby M dwarfs. However, conclusions are very hard to draw at the moment because of the low number of such candidates in the combined UVES and HARPS data. Bonfils et al. (2013a) did not observe any such planets in their sample even though the candidate orbiting GJ 433 could have been detected had they combined the HARPS data with the UVES velocities that were published in 2009. The same team of researchers reported this candidate in Delfosse et al. (2013), which means it could have been included in the results of Bonfils et al. (2013a).
Finally, out of the ten planet candidates in the sample, three appear to have orbits located within the respective stellar habitable zones according to the equations of Kopparapu et al. (2013) for the inner (moist greenhouse) and outer (maximum greenhouse) edge of the habitable zone. This enables us to state that 30% of M dwarf planets in the current sample are located within the HZs. However, we have to take into account the detection probability of the current data sets of habitable-zone planet candidates to estimate the occurrence rate of habitable-zone planets in a statistically representative, and thus meaningful, manner. We achieve this by calculating the detection probabilities in Fig. 6 as a function of semi-major axis instead of orbital period for each star and combine the resulting probabilities in the same way as in Section 2.6 but by limiting the analysis to the stellar habitable zones listed in Table 1. We find that the occurrence rate of planets with minimum masses between 3 and 10 M⊕ in the habitable zones of the sample stars is 0.21 +0.03 −0.05 planets per star and that of planets with minimum masses between 10 and 30 M⊕ is 0.035 +0.013 −0.007 planets per star.
These estimates can be compared with the results of Bonfils et al. (2013a) according to whom M dwarfs would have on average 0.41 +0.54 −0.13 super-Earths (1 M⊕ < mp sin i < 10 M⊕) per star in the habitable zones, which appears to be consistent with our estimate considering the slightly different mass range. These estimates also appear to be an order of magnitude greater than the occurrence rate estimates from the low-temperature sample of Kepler stars (Dressing & Charbonneau 2013). In particular, Dressing & Charbonneau (2013) estimated that the occurrence rate of Earth-size planets (0.5 -1.4 R⊕) is 0.06 +0.06 in the habitable zones of such stars. However, Kopparapu (2013) revised the estimates of Dressing & Charbonneau (2013) by calculating the habitable zones according to the modified equations of Kopparapu et al. (2013). As a result, there are 0.48 +0.12 −0.24 planets with 0.5 R⊕ < rp < 1.4 R⊕ per star in the habitable zones of the stars in the sample of Dressing & Charbonneau (2013), which increases to 0.51 +0.10 −0.20 when increasing the radius range to 2 R⊕. These estimates are broadly consistent with our results.
We also estimated the planetary mass function for lowmass companions with mp sin i < 100 M⊕. We used the same minimum mass bins as in Table 5 and plotted the resulting occurrence rate as a function of minimum mass in Fig. 7 (top panel) together with the mass distribution of the stars in the sample (bottom panel). This figure shows that the mass function increases dramatically with decreasing mass, similar to the increase found for planets orbiting more massive stars (Lopez & Jenkins 2012), which suggests that the high occurrence rate of super-Earths with mp sin i < 10 M⊕ corresponds to the rapid increase observed in the Kepler transit data for planets with radii less than 4 R⊕ (Morton & Swift 2013). However, the uncertainties are still too high to quantify this increase in a meaningful way and therefore we do not attempt to estimate the mass function quantitatively. Our occurrence rates denoted using the shaded histogram in Fig. 7  When calculating the occurrence rates under the assumption that the additional ten emerging signals, or SRCs, in the sample are also caused by planets, we obtained even higher occurrence rates for the low-mass planets around M dwarfs (Table 5, middle). These numbers are consistent with but slightly higher than those obtained for the ten candidate planets in the sample. However, we note that some of these signals could be false positives.

NEW PLANETARY SYSTEMS
We have listed the new planet candidates from our analyses of the sample velocities in Table 3. In this section, we discuss them briefly.

GJ 433
According to our results, we could also identify the two planetary signals that have been reported in the HARPS data of GJ 433 (Delfosse et al. 2013). We did not, however, detect any additional candidate planets satisfying our detection criteria in the combined HARPS and UVES data of the star.
We note that Bonfils et al. (2013a) also reported signals in the HARPS data around 30 days although there are significant problems with their tabulated solutions 9 . Nevertheless, the log-posterior of this target shows an emerging maximum at a period of 36 days for a three-Keplerian model together with a local maximum at 50 days, which suggests that a two-Keplerian model might not be a sufficiently accurate description of the combined UVES and HARPS velocities. If either one (or both) of these emerging signals is confirmed by future data, GJ 433 would become one of the highly populated planetary systems around M dwarfs together with the famous planet hosting stars GJ 581, GJ 667C, GJ 163, and GJ 676A.

GJ 682
In our sample, there are also two stars with two new candidate super-Earths orbiting them in the period space between 10 and 100 days (  17.478 [17.438, 17.540] and 57.32 [56.84, 57.77] days, respectively. The former candidate is located in the stellar habitable zone and can be classified as a habitable-zone super-Earth, although its minimum mass is consistent with (sub) Neptunian structure. An interesting detail in the model probabilities of the GJ 682 velocities is that the one-Keplerian model is only 9900 times more probable than the model without any signals, which implies that a one-Keplerian model is not sufficiently good description of the data to enable the detection of a planet candidate according to our criteria. However, a two-Keplerian model has a 1.8×10 7 times greater probability than the one-Keplerian model, which implies that there is very strong evidence in favour of two candidates orbiting the star. This means that when there are at least two signals of similar amplitude in the velocities, models with only one signal can be difficult to use to interpret the data because they are not good enough in describing the velocity variations.
GJ 682 is an inactive dwarf star with no signs of chromospheric activity (Walkowicz & Hawley 2009) and has a projected rotation period of 10.7 days (Reiners 2007). This suggests that the signals we have detected indeed are Doppler signatures of planets.
We note that possible additional signals in the GJ 682 data, and in particular their superpositions with the two detected ones, might cause biases to the obtained parameters of the two candidates listed in Table 3.

GJ 180
Another system with two candidate planets, GJ 180, corresponds to a remarkable configuration of super-Earths with an orbital period ratio of 7:5 -orbital periods of 17.380 [17.360, 17.398] and 24.329 [24.263, 24.381] days, respectively. This suggests (although does not imply) the existence of a stabilising 7:5 mean motion resonance (see also Jenkins et al. 2013a). The outer candidate is located in the stellar habitable zone, which makes this candidate interesting because of its reasonably low minimum mass of 6.4 [2.3, 10.1] M⊕ that enables its classification as a habitable-zone super-Earth. This system, its detection, dynamical stability, and formation history, are discussed in detail in a separate publication .

GJ 422
The candidate planet GJ 422 b has a minimum mass of 9.9 [5.9, 15.5] M⊕ and is thus a super-Earth or a sub-Neptunian planet, depending on the composition and atmospheric properties. Its orbit is likely located within the stellar habitable zone of GJ 422.

GJ 27.1
We also find a similar candidate orbiting GJ 27.1 with a minimum mass of 13.0 [6.4, 17.1] M⊕. This candidate can be classified as a sub-Neptunian planet candidate because it is likely too massive to be considered a super-Earth with rocky composition. With an orbital period of 15. 819 [15.793, 15.842] days this signal is unlikely to be caused by stellar activity coupled with the rotation period whose projected estimate is 11.9 days (Houdebine 2011).

GJ 160.2
In addition to GJ 433 b, GJ 160.2 b is the only candidate in our sample with an orbital period shorter than 10 days. Candidates of this kind are the easiest ones to observe in our sample (see Fig. 6) and the fact that we only found two of them demonstrates that such planets are not very common around M dwarfs as also quantified by the occurrence rates in Table 5. We note that GJ 160.2 might actually be a K dwarf (Koen et al. 2010), although Zechmeister et al. (2009) classified it as a M0 V star. The candidate GJ 160.2 b is a hot sub-Neptunian planet.
The star has a projected rotation period of 43.6 days (Houdebine 2011), which indicates that the signal in unlikely to be related to stellar rotation and thus activity.

GJ 229
Finally, we report a discovery of a planet candidate orbiting GJ 229 with an orbital period of 471 [459,493] days and a minimum mass of 32 [16,49] M⊕. This discovery makes the GJ 229 system one of the most diverse systems around M dwarfs because Nakajima et al. (1995) reported a brown dwarf companion to the star based on direct imaging.

CONCLUSIONS
We have presented our analysis of UVES velocities of a sample of 41 M dwarfs (Zechmeister et al. 2009) when combining the velocities with HARPS precision data as obtained from the spectra available in the ESO archive. As a result, we report the existence of eight new planet candidates around the sample stars (Tables 2 and 3) and confirm the existence of the two companions around GJ 433 (Delfosse et al. 2013) that exceed our conservative probabilistic detection threshold by making the statistical models more than 10 4 times more probable than models without the corresponding signals. Among the most interesting targets in our sample are GJ 433, GJ 180, and GJ 682, with at least two candidate planets each.
We have also presented estimates for the occurrence rate of low-mass planets around M dwarfs (Table 5) based on the current sample. We find that low-mass planets are very common around M dwarfs in the Solar neighbourhood and that the occurrence rate of planets with masses between 3 and 10 M⊕ is 1.08 +2.83 −0.72 per star. This estimate is likely consistent with that suggested based on the Kepler results for a sample of stars with T eff < 4000 K (Dressing & Charbonneau 2013;Morton & Swift 2013), although the comparisons are not easily performed because we could assess the occurrence rates of companions with periods up to the span of the radial velocity data of a few thousand days. On the other hand, we confirm the lack of planets with masses above 3 M⊕ on orbits with periods between 1-10 days. Such companions to low-mass stars have an occurrence rate of only 0.06 +0.11 −0.03 planets per star based on our sample.
There are nine targets in the sample that are also found in the sample of M dwarfs presented in Bonfils et al. (2013a): GJ 1, GJ 176, GJ 229, GJ 357, GJ 433, GJ 551, GJ 682, GJ 699, GJ 846, and GJ 849. Out of these nine stars, we found signals in the velocities of GJ 229, GJ 433, and GJ 682. Our results are essentially similar for GJ 433, for which Bonfils et al. (2013a) reported a signal at 7.4 days and the same group reported another long-period signal when analysing the HARPS data in combination with the UVES data analysed here (Delfosse et al. 2013 Bonfils et al. (2013a) did not report any such periodicities for these stars. We believe the reason is that we obtained HARPS-TERRA velocities from the HARPS spectra that are more precise for M dwarfs (Anglada-Escudé & Butler 2012), combined the HARPS velocities with the UVES ones which provides more information on the underlying periodic signals regardless of whether the signals can be detected in the two data sets independently or not, and accounted for correlations in the velocity data that could disable the detections of low-amplitude signals if not accounted for (Baluev 2013;Tuomi & Jenkins 2012;Tuomi et al. 2013b).
We have compared our results briefly with those obtained by using the Kepler space-telescope (e.g. Howard et al. 2012;Dressing & Charbonneau 2013) in Section 4. However, such a comparison is not necessarily reliable because the properties of Kepler's transiting planet candidates can only be discussed in terms of planetary radii and the radial velocity method can only be used to obtain minimum masses. Because of this, it is not surprising that there are remarkable differences that are unlikely to arise by chance alone. For instance, Dressing & Charbonneau (2013) estimated that there are roughly 0.15 +0.13 0.06 Earth-sized planets (radii between 0.5 and 1.4 R⊕) in the habitable zones of cool stars (with T eff < 4000 K) and that the nearest such planet could be expected to be found within 5 pc with 95% confidence. We calculated a similar estimate for candidates with masses between 3 and 10 M⊕ and obtained an occurrence rate estimate of 0.21 +0.03 −0.05 planets per star that appears to be higher than the estimate of Dressing & Charbonneau (2013) despite the fact that we cannot assess the occurrence rates of planets with masses below 3 M⊕ because we did not detect any such candidates orbiting the stars in the sample. However, these estimates can only be compared in detail with a range of robust planet composition and evolution models in hand, and is beyond the current work.
According to our results, M dwarfs have very high rates of hosting systems of low-mass planets around them and have a high probability of being hosts to super-Earths in their habitable zones. Together with the fact that radial velocity surveys can be used to obtain evidence for Earth-mass planets orbiting such stars, and the fact that M dwarfs are very abundant in the Solar neighbourhood, this makes them primary targets for searches of Earth-like planets, and possibly life, with current and future planet surveys. Figure A1. Solution for GJ 27.1 b. From left to right: estimated posterior density of the period from samplings with β < 1, where the horisontal lines denote probability thresholds of 10% (dotted), 1% (dashed), and 0.1% (solid) of the MAP estimate denoted by the red arrow; phase-folded MAP radial velocity curve, where UVES and HARPS measurements are denoted using red and blue circles, respectively; and the marginalised posterior densities of the period and velocity amplitude, where the solid curves denote Gaussian functions with the same mean and variance.      The HARPS-TERRA velocities obtained from the publicly available spectra in the ESO archive are presented in this section for all the targets for which at least two spectra were available. The secular acceleration has been subtracted from every HARPS-TERRA data set.