Issue 
A&A
Volume 585, January 2016



Article Number  A134  
Number of page(s)  24  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201526729  
Published online  11 January 2016 
The HARPS search for southern extrasolar planets
XXXVIII. Bayesian reanalysis of three systems. New superEarths, unconfirmed signals, and magnetic cycles ^{⋆,}^{⋆⋆}
^{1} Observatoire astronomique de l’Université de Genève, 51 Ch. des Maillettes, 1290 Versoix, Switzerland
email: Rodrigo.Diaz@unige.ch
^{2} HarvardSmithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA
^{3} Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife, Spain
^{4} Dpto. de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain
^{5} Physikalisches Institut, Universitat Bern, Silderstrasse 5, 3012 Bern, Switzerland
^{6} AixMarseille Université, CNRS, LAM (Laboratoire d’Astrophysique de Marseille) UMR 7326, 13388 Marseille, France
^{7} School of Physics and Astronomy, University of St. Andrews, North Haugh, St. Andrews, Fife KY16 9SS, UK
^{8} Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150762 Porto, Portugal
^{9} Institut d’Astrophysique et de Géophysique, Université de Liège, Allée du 6 Août 17, Bât. B5C, 4000 Liège, Belgium
^{10} European Southern Observatory, 3107 Alonso de Cordova, Vitacura, Casilla 19001, Santiago 19, Chile
^{11} Canada France Hawaii Telescope Corporation, Kamuela, 96743, USA
^{12} Department of Physics, University of Warwick, Coventry, CV4 7AL, UK
^{13} Cavendish Laboratory, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{14} Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169007 Porto, Portugal
Received: 12 June 2015
Accepted: 9 October 2015
We present the analysis of the entire HARPS observations of three stars that host planetary systems: HD 1461, HD 40307, and HD 204313. The data set spans eight years and contains more than 200 nightly averaged velocity measurements for each star. This means that it is sensitive to both longperiod and lowmass planets and also to the effects induced by stellar activity cycles. We modelled the data using Keplerian functions that correspond to planetary candidates and included the short and longterm effects of magnetic activity. A Bayesian approach was taken both for the data modelling, which allowed us to include information from activity proxies such as log R'_{HK} in the velocity modelling, and for the model selection, which permitted determining the number of significant signals in the system. The Bayesian model comparison overcomes the limitations inherent to the traditional periodogram analysis. We report an additional superEarth planet in the HD 1461 system. Four out of the six planets previously reported for HD 40307 are confirmed and characterised. We discuss the remaining two proposed signals. In particular, we show that when the systematic uncertainty associated with the techniques for estimating model probabilities are taken into account, the current data are not conclusive concerning the existence of the habitablezone candidate HD 40307 g. We also fully characterise the Neptunemass planet that orbits HD 204313 in 34.9 days.
Key words: techniques: radial velocities / methods: data analysis / methods: statistical / planetary systems
Based on observations made with the HARPS instrument on the ESO 3.6 m telescope at La Silla Observatory under the GTO programme ID 072.C0488, and its continuation programmes ID 183.C0972, 091.C0936, and 192.C0852.
Full Tables 3, 6, and 10 are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/585/A134
© ESO, 2016
1. Introduction
The detailed architecture of multiplanet systems is a key observable to constrain formation and evolution theories. Observations made over many years with very stable instruments are necessary to fully unveil the structure of planetary systems. Detecting companions on longperiod orbits and lowmass planets at shorter periods usually requires many dozens of radial velocity measurements (see e.g. Mayor et al. 2009; Pepe et al. 2011), especially when the candidates are found in multiplanet systems, as is common (Mayor et al. 2011; Lissauer et al. 2012, 2014). These data sets tend to span many years.
The HARPS search for extrasolar planets in the southern hemisphere (Mayor et al. 2003) has recently celebrated its tenth anniversary. The most inactive stars in the solar neighbourhood have been monitored continuously for over a decade, producing data sets with over two hundred measurements. These are expected to permit an indepth exploration of the planetary systems around them. However, even for these very weakly active stars, the presence of magnetic cycles complicates the detection of lowmass companions (Santos et al. 2010; Lovis et al. 2011a; Dumusque et al. 2011a). For this reason, additional observables are routinely obtained from the HARPS spectra: activity proxies based on the line fluxes (mainly the proxies based on the Ca II H&K lines, but recently also Hα and Na I D), the mean line bisector span, its fullwidth at halfmaximum (FWHM), etc. All of these can help identify activity cycles and ultimately correct for their effect on the radial velocity time series.
Which analysis method is best applied on these data sets has been the subject of some debate. While the classical frequentist approach of studying the velocity periodograms is known to have drawbacks (e.g. Sellke et al. 2001; Lovis et al. 2011b; Tuomi 2012), the alternative Bayesian model comparison has led to different numbers of signals reported for the same system. For the star Gl667C, for example, for which the periodogram analysis revealed two planetary companions (Delfosse et al. 2013), different groups have reported a number of planets ranging from two (Feroz & Hobson 2014) to six or even seven (Gregory 2012; AngladaEscudé et al. 2013) when they used Bayesian methods. Moreover, stellar activity can mimic planetary signals, and the nature of the detected signals is frequently difficult to identify (see the case of Gl581: Udry et al. 2007; Robertson et al. 2014) and not always agreed upon (AngladaEscudé et al. 2014, 2015; Robertson et al. 2015).
The radial velocity signal produced by stellar activity can be separated into two types of effects: the shortterm effect produced by active regions (spots and plages) that rotate in and out of view as the star revolves, and the longterm effect associated with changes in the global activity level of cyclic stars (e.g. Baliunas et al. 1995; Hall et al. 2007, 2009; Santos et al. 2010; Isaacson & Fischer 2010; Gomes da Silva et al. 2012; Robertson et al. 2013). The shortterm modulation is produced by the difference in flux and convective blueshift of active regions with respect to the surrounding photosphere (e.g. Saar & Donahue 1997; Hatzes 2002; Desort et al. 2007, see also Boisse et al. 2012; Dumusque et al. 2014). This creates a radial velocity signal whose frequency power is concentrated on the stellar rotational period and its harmonics (Boisse et al. 2011). This shortterm signal strongly depends on the detailed configuration of the active regions and is therefore difficult to model precisely because no clear correlation with activity proxies is systematically seen. On the other hand, longterm activity variations are related to global changes in the convective pattern of the star (Lindegren & Dravins 2003; Meunier et al. 2010; Meunier & Lagrange 2013) that are produced by a change in the typical number of active regions. The effect is therefore less sensitive to the details of individual active regions, and a clear correlation is seen with activity proxies such as (Noyes et al. 1984) and the width and asymmetry of the mean spectral line (Dravins 1982; Santos et al. 2010; Lovis et al. 2011a; Dumusque et al. 2011a).
Here we analyse the radial velocity data of HD 1461, HD 40307, and HD 204313; the data include more than ten years of HARPS observations. All three stars have been reported to host at least one planet. They all exhibit longterm variability in their activity levels, as would be produced by magnetic cycles, and their effect on the radial velocity data are clearly detected. However, a full cycle is only observed for HD 1461. We employed the traditional periodogram analysis to identify potential periodic signals and used a Bayesian model comparison to asses their significance. Taking the discussion from the previous paragraph into consideration, we take a different approach to model each of the activity effects: the shortterm variability is described using a nondeterministic model that does not require a detailed description of the activity signal, while the longterm variability is modelled assuming a simple relation between the activity effect on the radial velocities and the activity proxy . Effects produced by stellar pulsations and granulation have timescales of minutes and can be efficiently averaged out (Dumusque et al. 2011b) or treated as additional Gaussian noise.
The paper is organised as follows: Sect. 2 briefly presents the data used in the analysis, Sect. 3 reports the results from the spectroscopic analysis of the three target stars and their main characteristics, Sect. 4 describes the models employed, including the model for the short and longterm activity effects. In this section we also detail the technique used to compare and select the models, we discuss the algorithm for sampling from the model posterior distribution, and we present the choice of parameter priors. In the following three sections we describe the results for each system. Finally, we discuss the results and present our conclusions in Sect. 8.
2. Observations and data reduction
All three targets were observed as part of the Guaranteed Time Observations programme to search for southern extrasolar planets and its continuation highprecision HARPS programmes. The observations were reduced using the HARPS pipeline (version 3.5), and the stellar radial velocities were obtained through a weighted crosscorrelation with a numerical mask (Baranne et al. 1996; Pepe et al. 2002). The FWHM and bisector span of the peak in the crosscorrelation function (CCF) were also measured for each spectrum, as well as the activity proxy based on the Ca II H and K lines, , calibrated as described in Lovis et al. (2011a).
The number of measurements and basic characteristics of the observations studied here are presented in Table 1, and the nightly averaged radial velocity measurements are given in Tables 3, 6, and 10, available online. The HARPS observations of HD 204313 started around three years later than for the other two stars because this target was regularly monitored by the CORALIE instrument on the 1.2 m Swiss telescope at La Silla. For the analysis of HD 204313 we also included 104 RV measurements by CORALIE, 56 of which were obtained after the instrument upgrade performed in 2007 (Ségransan et al. 2010).
3. Stellar characteristics
Basic characteristics of the HARPS observations of the three target stars.
Observed and inferred stellar parameters.
The atmospheric parameters for the three stars studied here have been obtained by Sousa et al. (2008) based on HARPS spectra. Although in some cases now a larger number of spectra are available, the gain in signaltonoise ratio (S/N) and therefore in the precision of the parameters, is surely limited. We therefore decided to use the parameters as reported in Sousa et al. (2008), which are listed in Table 2. We note that the reported uncertainties for the atmospheric parameters do not consider potential systematic errors and may therefore be underestimated.
The stellar mass is given without uncertainty because the statistical error bar is certainly plagued with systematic errors (the choice of the stellar tracks, the physics used to compute the track, etc.). To compute the masses and semimajor axes of the detected companions, we conservatively fixed the uncertainty of the stellar mass to 10%. For all companions reported in this article, except for HD 1461 c and HD 40307 f, the contribution of the uncertainty in the stellar mass to the companion mass is larger than the contribution of the uncertainty in the orbital parameters. This illustrates the importance of improving our knowledge of the fundamental stellar characteristics, and the relevance of space missions such as Gaia and PLATO.
All three stars are magnetically quiet, with a mean below –4.9, but they all exhibit magnetic variability on timescales similar to the duration of the HARPS observations. These variations are reminiscent of the solar activity cycle and are described in more detail in the following sections. Despite their relative brightness and solarlike characteristics, these stars have not been systematically included in the southern surveys of stellar activity targeting solarlike stars (e.g. Henry et al. 1996; Cincunegui et al. 2007; Mauas et al. 2012). To the best of our knowledge, the only mention of one of the target stars in southern surveys appears in Arriagada (2011). Based on eight observations of HD 204313 from the Magellan planet search programme, Arriagada (2011) computed = – 5.0. Longterm observations of the magnetic activity level of HD 1461 exist from northern surveys (Hall et al. 2009), however. These observations show that the low level of activity of HD 1461 was maintained for at least 15 yr and seem to confirm the cyclic behaviour detected in the HARPS data (see Sect. 5).
The mean activity level is used to estimate the rotational period of the targets using the relation reported by Mamajek & Hillenbrand (2008). This estimate is also reported in Table 2.
4. Data analysis
4.1. Description of the models
The stellar radial velocity (RV) variations are described using a physical model M_{n} consisting of n Keplerian curves representing potential planetary companions and activity signals with timescales of the order of the rotational period, plus an additional signal produced by longterm stellar activity effects, a(t), which could take different forms depending on the knowledge we have on the stellar activity cycle (its period, etc.). The Keplerian functions plus the longterm activity signal constitute the deterministic part of the model (m). To this we add a statistical noise component ϵ that represents the stellar activity “jitter” – that is, the shortterm activityinduced variability that is not correctly modelled by a Keplerian curve – and all remaining systematic errors not considered in the reported uncertainties. As mentioned in the introduction, the shortterm activity signal depends on the details of the active regions visible at a given time on the stellar disk, their positions, sizes, and evolutions. This signal is therefore very hard to model using a single deterministic model, which is why we decided to add to it a statistical (nondeterministic)component that we call the stellar jitter. Therefore, the RV prediction for model M_{n} at time t_{i} can be written as (1)where k_{j} is the Keplerian curve of companion j. The Keplerian curves were parametrised using their period, amplitude, eccentricity, argument of periastron ω, and mean longitude at epoch L_{0}. We describe the statistical model for the stellar jitter in detail below, but we note that we explicitly added the subindex i to the noise component of the model ϵ_{i} to indicate that it can potentially depend on time. Additionally, the data errors are assumed to be uncorrelated and normally distributed.
Furthermore, assuming that the error term ϵ is distributed as a zerocentred normal, with standard deviation σ_{J}(t_{i}) = σ_{Ji} (potentially a function of time), the likelihood function of our model takes the form (Gregory 2005, Sect. 4.8) (2)where (v_{i},σ_{i}) is the velocity measurement and its associated uncertainty at time t_{i}.
Stellar jitter
The stellar jitter is included in our model as an additional, statistical error in the model prediction. We note that this is different from the approach taken by other authors (e.g. Tuomi et al. 2013a,b), who constructed a deterministic model of the stellar activity at short timescales. In this analysis, we make the strong assumption that the additional noise is uncorrelated and normally distributed.We note, however, that this term appears in addition to any potential rotational signal modelled by a Keplerian curve. It aims at accounting for the parts of the rotational activity signal that are not represented by the deterministic model. In that sense, the whitenoise assumption is probably less dramatic than if we were to use it to model the entire rotational activity signal. Two simple models of the amplitude of the stellar jitter were explored.
In the first one, the added jitter term has a constant amplitude σ_{J} for all times. In this case, the global value of σ_{J} is the sole parameter of the statistical model. However, we know beforehand that this model does not correctly describe the data because it is known that the dispersion in RV time series is larger for more active stars. To illustrate this, the constant jitter model was fitted separately to the RV measurements of HD 40307 (see Sect. 6) obtained before and after BJD = 2 454 800. As these data sets have different stellar activity levels, it is no surprise to find a clear difference in the distribution of the parameter σ_{J} (Fig. 1). Therefore, we decided to use a second model of the stellar activity jitter, in which the standard deviation of the noise component ϵ increases linearly with . The dependency on is motivated by the fact that the scatter in the RV measurements increases when the activity proxy does, but the linearity is an additional assumption of the model that needs to be tested. The second model was parametrised using the jitter level when =−5.0 and the slope of the dependence of the jitter amplitude with . The jitter level when =−5.0 (the baselevel jitter) represents any RV effect that might exist for solartype stars with such a low level of magnetic activity, such as the granulation noise (e.g. Dumusque et al. 2011b) and the undetected instrumental systematics changing from one night to another, which do not appear anywhere else in our model. The model requires an extra parameter and therefore suffers the Occam penalty described in Sect. 4.4. However, it is preferred by the available data for all three systems studied here, and we therefore only consider the evolving jitter model in the analyses presented below.
Stellar activity cycles
All three stars analysed here exhibit longterm activity variations reminiscent of the solar activity cycle. We claim that the effect on the RV time series is clearly detected and include it in the model in the form of the term a(t). The functional form of a(t) is not fixed a priori, but is instead taken from a fit to the time series. We tried a number of different models to fit the time series (see Sect. 5), but in all cases, the “shape” of the time series was used to model the longterm variations seen in the RV data.
To transform the variations in into variations in RV, a scaling constant α is included in the model. Previous studies exploring the effect of magnetic cycles on RV data have parametrised an equivalent scaling constant as a function of the effective temperature and the metallicity of the star (Lovis et al. 2011a). However, the dispersion around the fit is considerable, and we therefore decided not to include a prior for parameter α. We instead compared in each case the obtained result with the expected value based on the T_{eff}− [Fe/H] parametrisation.
Fig. 1 Stellar jitter for the constantjitter model, fitted to the RV data of HD 40307 before and after BJD = 2 454 800, which have different mean levels of activity (see Sect. 6 for details). 
We note that by modelling the effect of the activity cycle in this way, we are assuming a linear relation between the variations in the proxy and the longterm activityinduced RV. This is a different assumption from the one used for the jitter model, which states that the amplitude of the additional Gaussian noise scales linearly with .
4.2. Posterior sampling
To estimate the model parameter credible regions, we obtained samples from the posterior distribution using the Markov chain Monte Carlo (MCMC) algorithm described in Díaz et al. (2014) with normal proposal distributions for all parameters. The algorithm uses an adaptive principal component analysis to efficiently sample densities with nonlinear correlations.
To increase the efficiency of the MCMC algorithm, the starting point for the chain was chosen using the genetic algorithm (GA) implemented in the yorbit package (Ségransan et al. 2011). This drastically reduces the burnin period and guarantees that the entire parameter space has been explored. We nevertheless launched a number of independent chains to explore the possibility of multimodal posterior distributions. The chains were combined after thinning using their autocorrelation length.
The posterior of the stellar mass, which is not constrained in our model, was assumed to be a normal distribution centred at the value reported in Table 2, with a width equivalent to 10% of this value. A randomly drawn sample from this distribution was coupled to the MCMC sample of the remaining model parameters to obtain the posteriors of model parameters such as the planet minimum masses or the semimajor axes.
4.3. Choice of priors
The priors of the model parameters are presented in detail in the Appendix for each system. In general, the only parameters with informative priors are the orbital eccentricity and the parameters of the activity signal a(t). For the eccentricity we chose a Beta distribution, as advocated by Kipping (2013), who derived the shape parameters that best match a sample of around 400 RVmeasured orbital eccentricities (a = 0.867, b = 3.03). The priors for the longterm activity signal were chosen as normal distributions around the leastsquares fit to the time series, neglecting the covariances between the fit parameters (see Tables A.1–A.3). This is the practical way in which we incorporated the information present in the time series to our model.
For the remaining parameters we used uninformative priors (i.e. uniform or Jeffreys). The limits chosen for each parameter are shown in the tables in the Appendix.
4.4. Bayesian model comparison
One of the aims of our analysis is to establish the number of periodic signals present in a given RV data set, independently of their nature. Traditionally, this is addressed by studying the periodogram of the RV time series and by estimating the significance of the highest peak found. To do this, a series of synthetic datasets are obtained by reshuffling or permuting the original data points. The periodogram is computed on each newly created data set and the power of the highest peak is recorded. The histogram of the maximum peak powers is used to estimate the pvalue as a function of power level. This pvalue is estimated under the null hypothesis – in this case no (further) signal – since no real signal is expected in the reshuffled data sets. If the pvalue of the highest peak in the original histogram is lower than a predefined threshold, the bestfit Keplerian signal at the peak frequency is subtracted from the data and the periodogram analysis is again performed on the velocity residuals. This process is repeated until no further peaks appear above the threshold. Finally, a global fit including all detected frequencies is performed.
This technique has the advantage of being computationally inexpensive and is expected to produce the correct number of significant signals ifthe threshold pvalue is chosen to be low enough and provided the removed signals are well constrained. However, it hastwo main limitations: a) the interpretation of the pvalue as a falsealarm probability is in general incorrect and leads to an overestimation of the evidence against the null hypothesis (Sellke et al. 2001); and b) the uncertainties of the signals subtracted from the data are not taken into account when computing the statistical significance of any potential remaining signal (see e.g. Lovis et al. 2011b; Tuomi 2012; Hatzes 2013; Baluev 2013). Therefore, when dealing with signals with amplitudes below ~1 or 2 m s^{1}, which are similar to the activity signals and to the uncertainty of the individual observations, it is not advisable to conclude on the significance of the signals based on the pvalues obtained from the periodogram. In these cases we resorted to the more rigorous technique of Bayesian model comparison. We note, however, that the periodogram was used throughout the analysis to identify possible periodicities in the data, and when the associated pvalue was low enough (typically below 0.1%), more sophisticated analyses were not deemed necessary to declare the signal significant.
Bayesian statistics permits, unlike the frequentist approach, computing the probability (p) of any logical proposition, where the probability is understood to be a degree of plausibility for that proposition. In this framework, comparing two models (M_{1} and M_{2}) in the face of a given data set D and some information I can be made rigorously by computing the ratio of their posterior probabilities, known as the odds ratio: (3)where the first term on the righthand side is called the prior odds and is independent of the data, and the second term is the Bayes factor and encodes all the support the data give to one model over the other.
The Bayesian approach to model comparison treats models with different numbers of parameters and nonnested models. The Bayes factor has a builtin mechanism that penalises models according to the number of free parameters they have (known as Occam’s factor, see Gregory 2005, Sect. 3.5). We note that when there is no prior preference for any model (p(M_{1}  I) /p(M_{2}  I) = 1), the Bayes factor is directly the odds ratio. To compute the Bayes factor, the evidence or marginal likelihood of each model are needed, defined as the weighted average of the model likelihood (p(D  θ_{i},M_{i},I) = ℒ(θ_{i})) over the prior parameter space^{1}: (4)where θ_{i} denotes the parameter vector of model i, and π(θ_{i}  M_{i},I) is the parameter prior distribution. In multidimensional parameter spaces, such as those associated with models of multiplanet systems, the integral of Eq. (4) is often intractable and has to be estimated numerically. Moreover, the basic Monte Carlo integration estimate consisting of obtaining the mean value of ℒ(θ_{i}) over a sample from the prior density is expected to fail for highdimensional problems if the likelihood is concentrated relative to the prior distribution because most elements from the prior sample will have very low likelihood values.
A considerable number of methods for estimating the evidence exist in the literature (see Friel & Wyse2012, for a recent review and Kass & Raftery1995). Estimating the evidence is difficult in multidimensional spaces, and different techniques can lead to very different results (see e.g. Gregory 2007). We therefore decided to use three different methods and compare their results.All of them rely on posterior distribution samplesand are therefore relatively fast to compute because they use the sample obtained with the MCMC algorithm described above.In some cases, further samples are needed from known distributions from which they can be drawn in a straightforward manner.

The Chib & Jeliazkov (2001, hereafter CJ01) estimator is basedon the fact that the marginal likelihood is the normalising constantof the posterior density. The method requires estimating theposterior density at a single point in parameter spaceθ^{⋆}. To do this, a sample from the posterior density and from the proposal density used to produce the trial steps in the MCMC algorithm are needed. The method is straightforward and relatively fast, but can run into problems for multimodal posterior distributions (Friel & Wyse 2012). In this study, the sample from the proposal distribution was obtained by approximating the proposal density by a multivariate normal with covariance equal to the covariance of the posterior sample. The uncertainty was estimated by repeatedly sampling from the proposal density and using different subsets of the posterior sample.A weakness of this method as implemented here is the approximation of the proposal density. Moreover, computing the likelihood on the sample from this distribution is the most computationally expensive step in the process, which limits the sample size that can be drawn. Additionally, some of the draws from the proposal distribution fall outside the prior domain, reducing the effective sample size further.
Table 3HARPS measurements of HD 1461.

The Perrakis et al. (2014, hereafter P14) estimator is based on the importance sampling technique. Importance sampling improves the efficiency of Monte Carlo integration of a function over a given distribution using samples from a different distribution, known as the importance sampling density (see e.g. Geweke 1989; Kass & Raftery 1995).This technique can readily be employed to estimate the integral in Eq. (4). P14 proposed using the product of the marginal posterior densities of the model parameters as importance sampling density, which yields the estimator (5)where the p_{j}, with j = 1,...,q_{i}, are the marginal posterior densities of the model parameters, q_{i} is the number of parameters in model i, and θ^{(n)}, with n = 1,...,N, are the parameter vectors sampled from the marginal posterior densities. We here obtained the sample from the marginal posterior densities by reshuffling the N elements from the joint posterior sample obtained with the MCMC algorithm, so that correlations between parameters are lost, as suggested by P14. We note that if the marginal posterior sample is obtained in this way, no further draws are necessary from the posterior distribution, although computing the likelihood in the reshuffled sample is still necessary and is the most timeconsuming step in the estimation. The technique also requires evaluating the marginal posterior probabilities that appear in the denominator of Eq. (5). We estimated these densities using the normalised histogram of the MCMC sample for each parameter. The error produced by this estimation of the marginal posterior distributions is a weak point of our implementation because it increases with the number of parameters as a result of the product in the denominator. We are currently studying more sophisticated techniques such as the nonparametric kernel density estimation. The uncertainty was estimated by repeatedly reshuffling the joint posterior sample to produce new samples from the marginal posterior distributions.

Tuomi & Jones (2012) also used importance sampling to estimate the marginal likelihood. The importance sampling function ℐ is a mixture of posterior distribution samples at different stages of a MCMC: (6)The level of mixture (λ) and the lag between samples (h) are two parameters of the method that the authors explored. The result is called a truncated posteriormixture (TPM) estimate. This estimator is designed to solve the wellknown stability problem of the harmonic mean estimator (HME; Newton & Raftery 1994; Kass & Raftery 1995), which uses the posterior density as importance sampling density. The HME converges very slowly to the evidence (Friel & Wyse 2012; Robert & Wraith 2009) and usually produces an estimator with infinite variance (Robert & Wraith 2009). In addition, as the HME is based solely in samples from the posterior, which is typically much more peaked than the prior distribution, it will generally not be very sensitive to changes in the prior. This is documented in Friel & Wyse (2012) and is a clear drawback of the HME because the evidence is known to be extremely sensitive to prior choice. The TPM estimate aims at solving the stability problem by using a mixture for the importance sampling density.This estimator converges to the HME of Newton & Raftery (1994) as λ tends to zero, and therefore its variance also tends to infinity^{2}. However, when λ is different from zero, the TPM estimate is inconsistent, that is, it does not converge to the evidence as the sample size increases. In addition, the TPM estimate has the very important drawback of inheriting the priorinsensitivity of the HME. It is therefore unable to correctly reproduce the effect of Occam’s penalisation found in the Bayes factor. This is documented in the Tuomi & Jones (2012) article where TPM is introduced, but is presented as an advantage of the estimator. In summary, for λ = 0, the TPM estimate is equivalent to the problematic HME, and if λ> 0 the estimator is inconsistent. Therefore we do not expect this estimator to produce reliable results, but we included it for comparison.
5. HD 1461
HD1461 hosts a superEarth on a 5.77day period orbit. Rivera et al. (2010) reported its discovery based on 167 radial velocity measurements taken with HIRES on the Keck telescope over 12.8 yr. The presence of two additional companions in longer period orbits (P = 446.1 days and P = 5017 days) is also discussed by the authors. We analysed 249 nightly averaged HARPS measurements spanning more than ten years with a mean internal uncertainty of 49 cm s^{1}, which include photon noise and the error in the wavelength calibration.
A preliminary analysis of the HARPS radial velocities produced by the instrument pipeline revealed a periodic oneyear oscillation with an amplitude of ~1.4 m s^{1}. This oneyear signal has previously been identified as a systematic effect in HARPS data (Dumusque et al. 2015). Its origin is the manufacturing of the E2V CCD by stitching together (512 × 1024)pixel blocks to reach the total detector size (4096 × 2048 pixels). The spacing between these blocks is not as regular as the spacing between the columns within a block. Such discontinuities are at the moment not taken into account in the HARPS wavelength calibration. Despite the great stability of HARPS, the position of the stellar spectral lines on the detector varies throughout a year due to the changes in the Earth orbital velocity. Depending on the content of the spectrum and the systemic velocity of the star some spectral lines may go through these stitches and produce the observed yearly oscillation. This is the case for HD 1461, and we have corrected for this effect by removing the responsible lines from the spectral correlation mask. When this is done, the signal at one year disappears. The average uncertainty in the velocity increases around 13% due to the smaller number of lines used for the correlation. The velocities reported in Table 3 and plotted in Fig. 2 are the corrected version.
In Fig. 2 we plot time series of the RV, the activity proxy and two spectral line measurements (FWHM and bisector velocity span) that can also be affected by activity. A similar longterm evolution of all four observables is clearly visible, indicating the presence of a magnetic activity cycle (Sect. 5.2). The top panel of Fig. 3 presents the generalized LombScargle periodograms (GLS; Zechmeister & Kürster 2009) of the RV data. The periodogram is dominated by a signal with a period of 5.77 days, compatible to the planet candidate reported by Rivera et al. (2010). The amplitude of 2.37 ± 0.20 m s^{1} also agrees with Rivera et al. (2010) and corresponds to a minimum mass of around 6.4 M_{⊕}. In the remaining panels of Fig. 3 the GLS of the residuals around models with increasing number of Keplerian components are shown. Table 4 presents the Bayesian evidence of models with at least three signals (including the activity cycle; see below) and the associated Bayes factors with respect to the threeKeplerian model.The model probabilities are plotted in Fig. 4.
Fig. 2 Top panel: HARPS time series of HD 1461. The vertical dashed line separates the active (BJD < 2 454 850) from the inactive data set. Lower panel: GLS at periods over 400 days for the four time series plotted in the top panel. The vertical dotted lines represent the time span of observations and twice this value. 
5.1. TwoKeplerian model. A new superEarth candidate
The residuals around the oneKeplerian model (Fig. 3) reveal a new significant signal with P = 13.5 days and an amplitude of 1.5 m s^{1}, accompanied by its yearly and seasonal aliases. We employed the technique described by Dawson & Fabrycky (2010) to identify the peak corresponding to the real signal, but the data were not sufficient to obtain a definitive answer. A longterm signal associated with the magnetic cycle (see below) is also significant. There is no peak in the spectral window function that might indicate that the 13.5day peak is an alias of the longer activityinduced signal.On the other hand, the signal period is close to the first harmonic of the estimated rotational period. Boisse et al. (2011), among others, showed that activityinduced signals preferentially appear at the rotational period and its two first harmonics. However, the signal is recovered with the same period and amplitude if only the last five seasons of observations (BJD > 2 454 850) are considered, when the activity level of HD 1461 was at a minimum. This indicates that the signal is coherent over many years, which is not expected from a signal induced by stellar magnetic activity.
Furthermore, none of the , bisector or FWHM time series show any significant power at this period. The bisector velocity span time series exhibits a dispersion of 1.24 m s^{1} and 1.16 m s^{1}, respectively, before and after the degreethree polynomial is used to correct for the effect of the activity cycle (see below). The GLS of the bisector exhibits significant power at 29.2 days with an amplitude of around 60 cm s^{1}, most likely caused by the stellar rotational modulation (see Table 2). The FWHM time series does not present any significant periodicity when the longterm trend is removed and exhibits a dispersion of only 3 m s^{1} over more than ten years. The time series of the activity proxy, after correction for the longterm variation interpreted as the activity cycle, still exhibits power at periods ~500 days. As discussed below, this is surely due to an incorrect modelling of the activity cycle, which introduces aliasing frequencies in the periodogram of the corrected time series.
We conclude that the signal is best explained as produced by an additional planetary companion to HD 1461, with a minimum mass of around 6 M_{⊕}^{3}. The parameters of the new companion are listed in Table 5.
Fig. 3 Periodograms of the RV data of HD 1461 (top panel) and residuals around models with 1, 2, and 3 Keplerians. The horizontal lines are the 10% and 1% pvalue levels. 
5.2. Activity cycle and search for additional signals
A common longterm evolution is conspicuous in the time series plotted in the upper panel of Fig. 2. The periodograms in the lower panel of the same figure show peaks at periods of around 3000 days, close to the time span of the observations (Table 1). Changes in the bisector span throughout a solarlike magnetic activity cycle are expected from changes in the convective blueshift pattern (see for example Gray & Baliunas 1995; Gray et al. 1996; Dumusque et al. 2011a; Lovis et al. 2011a). These variations have a slightly longer period that the variation observed in the RV time series. The period of the FWHM is not yet constrained. We conclude that the longperiod signal seen in the RV time series is produced by a magnetic activity cycle, with a period of P_{cycle} = 9.64 ± 0.21 yr, as measured by a Keplerian fit to the time series. Hall et al. (2009) presented seasonally averaged measurements between late 1998 and late 2007. These data agree well with the trend observed in the HARPS time series and seem to confirm the amplitude of activity variations. On the other hand, the combined data set hints at a longer period and at a shorter active season around the year 2007. A fit of the combined data set gives a period between 16 yr and 18 yr.
The RV signature of the activity cycle has a period P = 9.1 ± 0.4 yr, an amplitude above 3 m s^{1}, and a significant eccentricity of e = 0.43 ± 0.07. It is the dominant feature in the residuals of the model, including the planets at 5.77 and 13.5 days. This signal has to be corrected for to continue searching for signals in the RV time series and to avoid mistaking an alias of this longterm variation with real periodic signals. For example, the oneyear aliases of the signal with the period of the activity cycle are located at 330 and 407 days. The GLS periodogram of the RV residuals shows significant power at these frequencies. On the other hand, the bestfit Keplerian curve to the has a slightly different period and eccentricity (e = 0.17 ± 0.07 for ) than the one for the RVs, as also seen in the periodograms of Fig. 2. An incorrect correction for the effect of activity can introduce spurious signals in the data.
Model probabilities for HD 1461.
Parameter posteriors for the HD 1461 system.
We therefore decided to study different functional forms for the activity function a(t) included in our model (Eq. (1)) and to compare the signals obtained under each method. A signal independent of the correction method intuitively has more support than a signal that is only found for one particular correction method. The activity cycle was included in the model of the RV data in two different manners:

a)
the variations are modelled using a sinusoidalfunction, and the bestfit parameters are used as Gaussian priorsfor a fit to the RV time series, with the exception ofthe sine amplitude, which is free to vary, and

b)
same as (a), but using a Keplerian function instead of a sinusoid.
Additionally, we tested other methods of removing the activity signal from the RV time series a priori:

c)
applying a lowpass filter (cutoff at 100 days) to the time series and using the filtered time series to detrend the RV data (see Dumusque et al. 2012), and

d)
running a principal component analysis on the combined and RV time series; the corrected RV are constructed by using only the second principal component, which is orthogonal to the direction of the joint variation of RV and .
We note that all the methods used to account for the activity cycle assume a linear relation between the variations observed in the proxy and those in the RV data. The alternative of fitting an additional Keplerian curve to the RV data without any prior information on the variations was discarded because it does not fully consider all available information. The search for additional signals was also performed on the RV data obtained after JD = 2 454 850, which correspond to the last five observing seasons and to the period of lesser magnetic activity, according to the proxy. These inactive data set contains 191 nightly averaged observations spanning 4.5 yr. The activity cycle is less prominent in these data and appears as a weak drift in the radial velocities.
Additional signals were searched for in the RV data using each of the models of the activity cycle and adding a further Keplerian signal to the model with two planets and the magnetic cycle. We initialised the MCMC algorithm using the bestfit solution of the twoKeplerian model for the parameters of the two superEarths and randomly drawing parameters from the prior joint density for the third potential planet. We note that although the two known planets were started at fixed points, no informative priors were used for their parameters, and they therefore were able to change freely if the data required it in the model with three planets. To thoroughly explore the parameter space, we launched 75 chains thus initialised. We list the priors used for each parameter in Table A.1.
As expected, the chains became trapped in the numerous local likelihood maxima associated with different values of the period of the putative additional planet. By comparing the value of log (ℒπ) in each maximum, where ℒ is the likelihood function and π is the prior probability density, those that clearly produced a much poorer fit were discarded.
Two signals were found irrespective of the method employed to model the activity cycle or correct its effect: a signal at 22.9 days with an amplitude of around 75 cm s^{1}, and another one around 620 days with an amplitude of 1.2 m s^{1}. We note that their frequencies are not significantly present in the periodogram of the residuals of the threeKeplerian model (Fig. 3). Table 4 presents the results of the Bayesian model comparison between models with and without these two additional periodicities, using a Keplerian model for the activity cycle with priors based on the time series, as explained above.
5.3. FourKeplerian model I. The 22.9day signal
The parameters of the 22.9day period Keplerian are approximately the same for different methods and for the inactive data set. In Table 4 both the CJ01 and P14 estimators indicate that the improvement in the data fit is not enough to justify the inclusion of the 22.9day signal. The Bayesian information criterion (BIC)^{4} is inconclusive in this respect. We therefore discarded the possibility that only the 22.9day signal is present. We tested, on the other hand, the inclusion of both 22.9day and 620day signals, but this model is not favoured by the data, probably due to the larger number of parameters. As discussed above, the TPM estimator overestimates the evidence for all cases and does not incorporate the Occam penalisation correctly, leading to a preference for more complex models, as is clearly seen in Table 4.
Fig. 4 HD 1461. Odds ratio for models with n Keplerian curves with respect to models with n−1 Keplerian curves as a function of model complexity n, assuming equal unity prior odds in all cases. The estimates using different techniques are shown and the customary limits for positive (O_{n + 1,n} = 3) and strong (O_{n + 1,n} = 20) and their inverses are shown as dashed and dotted lines, respectively. The model with four signals contains the 620day Keplerian. 
Fig. 5 Marginal posterior distributions of the orbital period, eccentricity, phase, and signal amplitude for the signal at around 620 days for each method used to account for the RV effect of the activity cycle. Also included are the posterior distributions using only the RVs obtained during the period of lower activity. 
5.4. FourKeplerian model II. The 620day period signal
Fig. 6 HD 1461. GLS periodogram of the RV residuals of the twoKeplerian model (top) and twoKeplerian + linear drift (bottom) for data taken after JD = 2 454 850 (the inactive data set). The two peaks standing out as significant signals in the top panel have periods of 22.9 days and around 650 days. Note that the significance is reduced drastically when the longterm trend caused by the activity cycle is removed, indicating that the observed periodicities are aliases of a longperiod signal present in the data. 
Both the CJ01 and P14 techniques favour the inclusion of a signal at 620 days, with a Bayes factor of 60 and 20, respectively, which is considered as strong evidence (Kass & Raftery 1995). The estimation based on the BIC leads to a similar conclusion.
For all models used to describe the effect of the activity cycle of HD 1461 on the RV measurements, a signal at around 630 days is seen, albeit its period changes slightly with the method employed. Methods (a) and (b) produce a signal closer to 615 days, while for methods (c) and (d), the signal is found closer to 640 days (Fig. 5). For all methods, the amplitude is compatible. If this signal is of planetary origin, the minimum mass of the companion would be M_{d} ~ 14.5 ± 1.3M_{⊕}. No significant power is present at similar periods in the time series of the , the bisector velocity span, or the FWHM, even after subtracting the longterm trend associated with the magnetic cycle.
If only the inactive data set (BJD > 2 454 850) is considered, the GLS periodogram of the residuals of a twoKeplerian model exhibits significant (pvalue <0.01) power at the period of the signal (Fig. 6). However, when a linear drift is added to the model to account for the effect of the activity cycle, the amplitude of the peak is strongly reduced (Fig. 6), indicating that the periodicity may be an alias of the longterm trend and not a real signal. The 22.9day signal exhibits the same behaviour. This would explain why the signal is recovered for all the correction methods of the activity cycle, as well as the slight change of the period under different corrections. Since a longperiod signal must remain in the data for the alias frequencies to be present, this either means that the correction of the activity cycle is not fully satisfactory with any of the methods or that an additional longterm signal, still not fully sampled, is present in the data. We conclude that although the periodicity at 620day period is significantly present in the data, its nature is still uncertain and might originate in an incomplete correction of the activity cycle.
5.5. The planetary system around HD 1461
Our final model of the RV series includes two Keplerian curves for the known planet candidates at 5.77 and 13.5 days and an additional Keplerian curve to model the activity cycle. The planet signals are independent of the method used to model the activity cycle. For simplicity, we chose the Keplerian model, which also allows us, unlike the filtering and principal components method, to include the uncertainties in the parameters in the error budget of the planet signals.
The resulting posterior distributions for the semiamplitude and the orbital eccentricity are shown in Fig. 7, which clearly shows that the three signals have amplitudes significantly different from zero. The covariance between the three semiamplitudes and the eccentricities ismuch smaller than the variance of each parameter. In Table 5 the mode and 68.3% credible intervals are listed for all MCMC parameters and for a series of derived parameters. The reflex motion induced by the new companion at the 13.5day period has an amplitude of 1.49 ± 0.17 m s^{1}, which implies a minimum mass of 5.59 ± 0.73M_{⊕}, in agreement with the values reported by Mayor et al. (2011). The RV amplitude associated with the activity cycle is 1.51 ± 0.26 m s^{1}, which means that the scaling constant α between and RV is 74.2 ± 12.8 m s^{1}/dex, in good agreement with the value calibrated as a function of effective temperature and metallicity by Lovis et al. (2011a), which gives 74.5 m s^{1}/dex.
The RV data folded to the bestfit period of each signal are shown in Fig. 8 after subtracting the effect of the remaining signals. This correction was made by computing the model corresponding to each Keplerian curve for each step of the MCMC, sampled at the data times. The mean value of these models in each data time was subtracted from the observed data. The blue curves in Fig. 8 are the 95% highest density interval (HDI) of the curve sampled in 300 phase points. We computed this in a similar way as Gregory (2011): the period of each signal was sampled at 300 points, and the corresponding RV model was computed for each posterior sample element obtained with the MCMC algorithm.
To study the stability of the system, we performed a numerical integration of the system over half a million years using the Mercury code (Chambers 1999). Two simulations were run: the first using the minimum masses as the true masses of the companions and coplanar orbits, and the second one increasing the masses by a factor two and including a mutual inclination of ten degrees. The initial eccentricities were set to the 95% upper confidence level. In both cases the system was stable over the integrated time scale. Additionally, the eccentricities did not increase beyond 0.24 and 0.23 for planet b and c, respectively. The fractional semimajor axis change is smaller than 10^{6} for the outer companion and around 8 × 10^{5} for the inner one, which is similar to the precision of the integrator.
The residuals of the model with three Keplerians still show significant scatter (2.3 m s^{1}), which forces the additional noise component of the model to be 1.7 ± 0.1 m s^{1} for the mean value. This is caused partially by the remaining signal originated in an incomplete cycle correction and by other effects that were not taken into account in our model, such as rotational modulation of the RV data due to stellar spots. It could also be indicative of additional planets in the system. Further observations of this system are needed to fully characterise it.
Fig. 7 Posterior distributions of the amplitude (top row) and eccentricity (bottom row) of the three Keplerian curves used to model the HARPS radial velocities of HD 1461. The grey dotted curves represent the eccentricity prior for the planetary signals. To facilitate comparison, the axis scales are the same for the three signals. 
Fig. 8 Radial velocity data phasefolded to the bestfit period of each of the three Keplerian curves used in the final modelling of HD 1461, after subtracting the effect of the remaining signals. The error bars include the additional noise term. The blue lines represent the 95% highest density interval (HDI). 
Rivera et al. (2010) reported two potential signals with periods around 450 and 5000 days in their HIRES data set. In the light of the present analysis, their detection might be related to an incompletely sampled magnetic cycle, although more data are needed to reach a firm conclusion on the nature of those suggested periodicities.
6. HD 40307
HD 40307 was reported to host three superEarthtype planets with orbital periods P_{b} = 4.311 d, P_{c} = 9.6 d, and P_{d} = 20.5 d by Mayor et al. (2009, hereafter M09 based on 2.4 yr of HARPS data. More recently, Tuomi et al. (2013a, hereafter T13 analysed the publicly available HARPS data, which included the M09 data and observations taken on three additional nights, and claimed the presence of three additional planets in the system, with orbital periods P_{e} = 34.62 d, P_{f} = 51.76 d, and P_{g} = 197.8 d. Planet g would be in the habitable zone of the star. They also detected a periodic signal with P ~ 320 days that they attributed to magnetic activity effects because its amplitude changes depending on the fraction of the spectrum used to compute the radial velocities. The analysis by T13 differs from the one by M09 mainly in the way the radial velocities are obtained – by template matching instead of mask crosscorrelation – and in that they used the complete HARPS data set instead of the nightly binned velocities, including seven points taken during the commissioning of the instrument. Additionally, T13 used a movingaverage model to take into account the correlation between individual observations taken during a single night, a deterministic model of the shortterm activity signal. The analysis presented here includes four additional years of data, for a total of 226 nightly averaged radial velocity measurements taken over eight years. This represents around 70% more data points than used by M09. We list the data in Table 6.
HARPS measurements of HD 40307.
Fig. 9 Upper panel: HARPS time series of HD 40307. For , the empty circles are the data included in M09, and have ⟨⟩ =−4.99. The filled circles are the new data presented here, with ⟨⟩ =−4.87. The vertical dashed line separates the active (BJD < 2 454 800) from the inactive data set. Lower panel: corresponding GLS periodograms for periods longer than 400 days. The vertical dotted lines represent the time span of observations and twice this value. 
The radial velocity data set is plotted in Fig. 9 together with the time series of the , the bisector velocity span, and the FWHM of the CCF. A lowfrequency signal is clearly visible in all the observables. As for HD 1461, it seems reasonable to assume that this longterm trend is linked to a stellar magnetic cycle. This signal is currently stronger than the reflex motion produced by the known companion at a fourday period, which illustrates the hindering effect magnetic cycles have on the detection of lowmass planets. The period of the signal is largely unconstrained, and we therefore decided to model it with a thirddegree polynomial instead of a periodic function.
It is interesting to compare the evolution of the time series as the level of activity changes. Throughout this section we consider the inactive and active data sets, corresponding to the observations obtained before and after BJD ~ 2 454 800, respectively (Fig. 9). The inactive data set includes only one additional observing night with respect to the data presented by M09, and two less than the data analysed by T13. The inactive period has ⟨⟩ =−4.99, and a typical dispersion of 0.014. In the active data set we find ⟨⟩ =−4.87 and rms = 0.024 dex over 90 observations, that is, an increased level of activity and significantly larger dispersion. We expect these differences to reflect on the radial velocities. As illustrated in Fig. 1, the additional noise increases from 1.0 m s^{1} to 1.7 m s^{1} between the inactive and active data sets and justifies the use of the varyingjitter model.
Fig. 10 Periodogram of the RV data of HD 40307 (top panel), and of the residuals around models with three (second from top), four (second from bottom), and five (bottom) Keplerian signals in addition to a cubic function to take into account the longterm trend produced by the magnetic cycle of the star. The horizontal dotted and dashed lines represent the 10% and 1% false alarm probability levels, respectively. 
In the top panel of Fig. 10 we present the periodogram of the HARPS RV data. The periods of the three planets reported by M09 (4.3, 9.6, and 20.5 days) are seen as narrow spikes in the periodogram. The longterm trend is also present. The remaning panels in Fig. 10 present the GLS of the residuals to fits with three, four, and five Keplerian signals, plus an additional thirddegree polynomial to account for the activity cycle. The corresponding model evidenceestimates are listed in Table 7 and plotted in Fig. 14. The methods of CJ01 and P14 agree remarkably well for models with up to four planets, with differences smaller than 1.4 in log p(D  M_{n},I). As expected, the TPM estimator of Tuomi & Jones (2012) largely overestimates the evidences.
In all models with at least three planets, the system announced by M09 is recovered, albeit with a slightly shorter period for planet d (P ~ 20.42 d instead of P ~ 20.46 d). This is probably due to the effect of activity at a similar period. Indeed, a significant peak appears at P = 21.4 days in the bisector time series (Fig. 11) when the longterm trend is corrected. However, when a leastsquares fit is performed on each observing season individually, the amplitude of the bisector signal is seen to anticorrelate with the one in the RVs. The bisector amplitude varies from below 50 cm s^{1} during the first three seasons to around 2.5 m s^{1} when the activity increases. If the signal in the RV data were produced by magnetic activity, we would expect a correlation to exist between its amplitude and that of the bisector signal. The fact that an anticorrelation is seen indicates that the activity signal is scrambling the signal seen in the RV, but does not cast doubt on its interpretation as a planetary companion. Otherwise, the amplitude and eccentricity distributions of the three companions are compatible in all models. The baselevel additional noise is below 1 m s^{1} for all models with at least three signals (Fig. 12), illustrating the high precision of HARPS. As the level of complexity of the model increases, the needed additional noise level decreases. For the five and sixsignal models, the noise level is around 60 cm s^{1}. On the other hand, the models with a weaker baselevel jitter have a higher sensitivity to , that is, a larger slope parameter.
Model probabilities for HD 40307.
Fig. 11 GLS periodogram of the bisector velocity span of HD 40307. A peak at the period of planet d is clearly detected. 
6.1. FourKeplerian model. A superEarth companion on a 51.6day period orbit
A 51.6day period signal appears as significant when the model with three Keplerian signals and a degreethree polynomial is subtracted (Fig. 10). Given that a longterm signal was subtracted, particular care should be given to the spectral window: as seen for HD 1461, if the signal is not correctly corrected for and some power remains at very long periods, peaks will appear at the frequencies present in the spectral window function. In this case, no peak is present at frequencies corresponding to ~51 days in the window function of the HD 40307 data. The 51.6day signal has an amplitude of 75 cm s^{1}. All models with at least four Keplerians converge to a period of P ~ 51.6 days, with an eccentricity distribution that in all cases is compatible with a circular orbit. Its amplitude, however, depends mildly on the model (Fig. 13). In the model with four signals, the amplitude is around 75 cm s^{1}, while in more complex models the amplitude is closer to 85–90 cm s^{1}.
When this fourth Keplerian is included, the model probability increases by a factor of 9.3 ± 1.7 or 22.6 ± 2.8 using the CJ01 and P14 estimates, respectively. This corresponds to positive and strong evidence in favour of the fourth signal, according to the scale presented by Kass & Raftery (1995). The evidence estimates based on these two techniques agree within 30%, which is remarkable given the difficulties associated with estimating the evidence in highdimensional spaces (Gregory 2007).
We note that this signal is not far from the rotational period estimated based on the level (Table 2). Indeed, the active period of the bisector velocity span exhibits a significant peak (pvalue < 0.01) at P = 51.5 days. However, the peak power is reduced to below the level of the pvalue = 0.1 when the data are detrended to account for the longterm evolution. Additionally, no equivalent peak is seen in the inactive period. Although this may cause concern at first sight, the 51.6day signal is present in the RV data set even after the longterm effect has been removed. The fourth Keplerian signal is therefore probably not attributable to stellar activity. The parameters are presented in Table 9.
Fig. 12 HD 40307. Posterior distribution of the amplitude the additional white noise for an inactive star = – 5.0. The dashed vertical lines represent the mean of each distribution. 
Fig. 13 HD 40307. Posterior distribution of the amplitude of the fourth signal. The dashed vertical lines represent the mean of each distribution. The model with four Keplerians produces a slightly weaker signal than the more complex models. 
Fig. 14 HD 40307. Odds ratio for models with n Keplerian curves with respect to models with n−1 Keplerian curves as a function of model complexity n, assuming equal unity prior odds in all cases. The estimates using different techniques are shown and the customary limits for positive (O_{n + 1,n} = 3) and strong (O_{n + 1,n} = 20) and their inverses are shown as dashed and dotted lines, respectively. The odds ratio in favour of the model with three Keplerians is out of scale. 
6.2. FiveKeplerian model. A doubtful 205day period signal
The periodogram of the residuals around the fourKeplerian model exhibits peaks at ~200 days and 28.6 days with power above the pvalue = 0.01 level (Fig. 10). The signal at 28.6 days is most probably an alias introduced by an incomplete correction of the longterm variability. The most prominent peak of the window function after the oneyear peak is at 27.8 days, probably linked to the Moon orbital period^{5}. Including a fifth Keplerian in the MCMC model leads to a period of P ~ 205 d for the fifth signal, with an amplitude of 80 cm s^{1} and an eccentricity posterior distribution peaked at 0.32 but compatible with zero at the 90%level.
Up to this point, the estimates of the model evidence provided by the CJ01 and P14 techniques have agreed, at least qualitatively. However, they disagree strongly concerning the significance of the fifth Keplerian signal. While based on the P14 estimator the fiveKeplerian model is around 15.8 times more probable than the model with only four signals, the CJ01 estimate indicates that the fourKeplerian model is 2.6 times more probable than the fiveKeplerian model.
The difference is certainly due to the approximations and limitations of each estimator and highlights the difficulty of obtaining a robust estimate of the integral in Eq. (4) for highly dimensional models. This problem was also reported by Gregory (2007) in his study of the radial velocity data of HD 11964. The author compared three different techniques for computing the evidence and found large differences (factors of around 20) for models with more than two planets. In our case, the estimates for models with n = 2,3 differ by a factor of a few (3–4), and even less for the model with four planets. On the other hand, a difference of three orders of magnitudes exists for the model with five Keplerians.
As final value for the evidence of the fourKeplerian model, we adopted the geometric mean between both estimates, with an uncertainty that contains both estimates (Table 8). We did this to account for the systematic errors associated with the estimation technique. The values in Table 8 indicate that it is not possible to decide between models with four and five signals with the current data. However, it can still be claimed that the model with four signals is more probable than the model with three planets (between 7.4 and 28 times).
On the other hand, the spectral window function of the analysed data exhibit peaks near this period (Fig. 18)^{6}. Even if present in the RV data, the observed signal at 205 days could also be an alias introduced by an incomplete correction of the activity cycle. Moreover, with a period of around 200 days, the putative planet would be well inside the habitable zone of HD 40307 (T13). For such planet it is particularly important to have strong and robust evidence for this signal and its nature. We therefore decided to remain cautious and to conservatively retain the model with four signals over the one including the 205day variability.
6.3. SixKeplerian model
The residuals of the model with five Keplerians still show power at 28.6 days, but with a much smaller amplitude (Fig. 10, lower panel). No additional signal is significant according to the periodogram analysis. However, as discussed above, the periodogram can produce an incorrect estimate of the significance of a signal, and we therefore estimated the Bayesian evidence for a model with six signals.
When adopting this model for the MCMC algorithm the fifth Keplerian does not converge to a clear maximum. There is posterior mass at periods of ~200 days, but also at around 28 days and 330 days (see Fig. 10). The sixth Keplerian also exhibits a multimodal behaviour, with power at 200 days, 34.6 days, and 20.5 days. This shows that the fifth and sixth Keplerians, if present in the data, do not have a clearly defined period in this model. We note that the signal at 34.6 days and the one at around 200 days have been announced by T13 as planets e and g, respectively. We discuss this further in Sect. 6.5.
The Bayesian model comparison rules out the sixsignal model as less probable than models with three or four Keplerians (Table 7). The discrepancy between the two estimates becomes much stronger for this model than for the fourKeplerian one, but the conclusion holds even if the systematic error between the estimators is taken into account (Table 8).
Summary of the Bayesian evidence computation for each model.
6.4. Signals in other observables. The stellar rotation period
The remaining observables obtained routinely from the HARPS spectra are used to validate the nature of the planetary signals. They often exhibit a high stability. In the case of HD 40307, for example, the dispersion of the bisector measurements is 1.25 m s^{1} over more than ten years and only 0.86 m s^{1} for the inactive data (defined in Sect. 6; around three years of data). This shows the exquisite precision and longterm stability of the HARPS spectrograph.
A significant period of around 37.4 days is found in the active time series of the FWHM, after correcting for the longterm evolution using a seconddegree polynomial fit (Fig. 15). A peak at the same period appears in the GLS periodogram of the active time series, albeit with pvalue > 0.01 (see, however, Tuomi 2012; Lovis et al. 2011b), and is also seen in the BIS time series at a slightly different period (39.1 days). We interpret this as the signature of the rotational period of the star, which is estimated to be around P = 47.9 ± 6.4 days (Table 2) taking the entire data set, and P = 41.6 ± 1.7 days (Mamajek & Hillenbrand 2008) if only the active part of the is considered.
Other signals are present in the time series. Considering the entire data set, a significant periodicity at 27.2 days appears, which is probably related to the first harmonic of the rotational period. The inactive time series exhibits power at 340 days and 41.2 days. Both periods have also been reported in the analysis of the S index by T13. The former is probably due to the longterm evolution of the stellar activity, while the latter is related to the rotational period as well. It is tempting to hypothesise that the difference with the period found in the inactive data set (P = 39.4 days) is due to the migration of active regions towards lower (faster) latitudes as the cycle progresses.
Signals at periods of between 1000–2000 days are seen in the and FWHM time series as well. Most, if not all of them, are probably caused by an imperfect detrending of the longterm evolution, which introduces periodicities at frequencies associated with the spectral window function.
Fig. 15 Periodogram of the active FWHM time series of HD 40307 after detrending using a thirddegree polynomial. The single peak reaching thepvalue = 0.01 level (indicated by the dotdashed lines) has P = 37.4 days. We interpret it as a signature of the stellar rotation period. 
6.5. The planetary system around HD 40307
Parameter posteriors for the HD 40307 system.
HD 40307 hosts four planetary candidates. We confirm the candidate announced by T13 with a period of 51.6 days, although the amplitude is slightly smaller in our analysis, implying a decrease of the minimum mass from 5.2 M_{⊕} to 3.6 M_{⊕}. The posterior distributions of the model parameters are presented in Table 9. The histogram of the MCMC sample for the signal amplitude and eccentricity are plotted in Fig. 16. All planets are on orbits that are compatible with circular. Indeed, the circular model explains the data equally well, and because it has fewer parameters, it is about 40 times more probable than the full Keplerian solution. We have, however, retained the eccentric model to provide upper limits to the eccentricity. In Fig. 16 the posterior distributions are more concentrated than the corresponding prior distributions, indicating that that the current data constrain the planet eccentricities beyond the prior level, albeit mildly for companion f. All the planets in the system have minimum masses below 10 M_{⊕} and can therefore be classified as superEarths. The least massive planet of the system is planet f, with a mass similar to that of the innermost planet (around 3.7 M_{⊕}). In Fig. 17 we plot the phasefolded velocity variation for each planet candidate after subtracting the effect of the remaining ones and the longterm drift. We chose different symbols for the active and inactive data set to show the increased dispersion during the active period (see also Fig. 1).
Dynamical integrations of the system were preformed over half a million years similar as for HD 1461. The mutual inclination of the planetary orbits were randomly drawn from a uniform distribution extending between –5 and +5 deg, which is a conservative assumption in the light of the observed distribution of mutual inclinations (Figueira et al. 2012). The true masses were set to twice the minimum mass. The system is stable over the explored timescale, the eccentricities remain lower than 0.34 – for the outer planet –, and the semimajor axes do not evolve significantly: the fractional changes are around 10^{6}, except for the inner planet, which exhibits a fractional change of ~10^{4}, related with the precision in the energy conservation of the system.
Two other signals were reported in this system by T13. Their companion g at P = 197.8 days could correspond roughly to the signal detected at 205 days –althoughthe upper limit of their 99% HDI is 203.5 days – but as discussed above, its significance and interpretation are doubtful. The periodicity of candidate e is found in the sixKeplerian model, but this model is clearly disfavoured by the data. Additionally, signals with periods P ~ 41–42 days and P ~ 37–39 days were found in the , bisector, and FWHM time series, indicating that these frequencies are associated with activity phenomena; they are close to the period of candidate e. More importantly, in the spectral window function of the data set employed by T13 the ~236day peak is much more prominent than in our data set, probably because of the longer observation timespan (Fig. 18). The 236day alias of P = 41.9 days (the rotational period present in the inactive time series) is P = 35.6 days, not far from the frequency of candidate e. Depending on the detailed structure of the sampling window, the aliased period may have more power than the real signal (Dawson & Fabrycky 2010).
Concerning the longterm drift attributed to stellar activity, the posterior distribution of the polynomial coefficients resemble closely the priors imposed from the analysis of the time series (Table A.1). The conversion constant between and RV is found to be (22.2 ± 2.2) m s^{1}/dex, which is higher than the value expected for this type of star according to Lovis et al. (2011a), which is around 12.1 m s^{1} per unit of , but where the error in the calibration coefficients has not been considered. This could indicate the presence of longperiod planets in the system. However, the dispersion around the fit by Lovis et al. (2011a) is large, and it would be premature to conclude based on this discrepancy.
Fig. 16 Posterior distributions of the amplitude (top row) and eccentricity (bottom row) of the four Keplerian curves used to model the HARPS radial velocities of HD 40307. The grey dotted curves represent the eccentricity prior. To facilitate comparison, the axis scales are the same for all signals. 
Fig. 17 Radial velocity data phasefolded to the bestfit period of each of the four Keplerian curves used in the final modelling of HD 40307 after subtracting the effect of the remaining signals and the longterm drift. The error bars include the additional noise term (see text). The inactive data set is plotted using filled red circles, while for the active data set we chose lighter empty circles. This promptly shows that the dispersion around the model is largely caused by the active data set. The blue lines represent the 95% highest density interval (HDI). 
Fig. 18 HD 40307. Window function around the oneyear alias (0.00273 cycles/day) for the data set used by T13 (solid curve), consisting of the inactive data set plus a few points taken during commissioning. The corresponding spectral window function for the data set analysed here is shown as a dotted curve. The frequency corresponding to 236 days is indicated by the vertical dashed line. 
7. HD 204313
Ségransan et al. (2010) discovered a companion to HD 204313 with msini = 4M_{J} on a 5yr period orbit based on CORALIE measurements. HARPS data were obtained on this system starting in 2006, and they permitted Mayor et al. (2011) to announce the presence of an additional Neptunemass planet on an inner orbit of 34.9 days. More recently, Robertson et al. (2012) combined the Ségransan et al. (2010) data set with 36 RV measurements obtained at McDonald observatory that span almost eight years. They reported an additional 1.7M_{J} companion on an outer orbit with P ~ 2800 days. The authors stated that the combined baseline of both data sets is responsible for the detection of this outer planet in 3:2 resonance with the first detected companion. Here we revisited this system using 28 new CORALIE data points taken between August 2009 and October 2014 and 95 HARPS data points obtained between May 2006 and October 2014, together with the already published data from CORALIE and McDonald data. Of the 104 CORALIE data points, 56 were obtained after the instrument upgrade in 2007, which mainly increased the efficiency of the spectrograph. This led to a reduced photon noise uncertainty on a given target and exposure time, but did not improve the overall precision of around 5 m s^{1}. Additionally, the hardware changes (see Ségransan et al. 2010 for details) introduced a zeropoint offset of around 20 m s^{1}. For this reason, the data taken before and after the upgrade were modelled as if they had been taken by different instruments^{7}. The data span 13 yr and contain 215 points. The weighted mean error of the measurements is 63 cm s^{1}. Two HARPS data points (BJD = 2 453 951 and BJD = 2 456 468) were discarded because they had an unusually low signaltonoise ratio or exhibited an anomalous CCF (high contrast), probably linked to an incorrect background correction. The data are presented in Table 10.
HARPS and CORALIE measurements of HD 204313.
The RV time series are plotted in Fig. 19. The variability produced by the planet first reported by Ségransan et al. (2010) is clearly seen by eye. When the data are blindly searched for signals using the GA, a 2050day nearly circular orbit is found, that is, a period slightly longer than originally reported, but in agreement with the findings of Robertson et al. (2012) and Mayor et al. (2011).
7.1. TwoKeplerian model. A Neptunemass object on a 35day period orbit
The residuals of the oneKeplerian model show a sharp significant peak at 34.9days and an additional trend with unconstrained period (Fig. 19). The signal at 34.9 days corresponds to the Neptunemass planet announced by Mayor et al. (2011), whose full discovery report is given here. The detection of this signal is solely due to the HARPS data. Indeed, without the 93 HARPS measurements, the GLS periodogram of the residuals to the oneKeplerian model does not exhibit any significant peak (Fig. 19)^{8}. The other instruments are useful for constraining the period of the massive outer planet, but contribute only negligibly to the identification of this new planet candidate.
The estimates of the rotational period based on Noyes et al. (1984) and Mamajek & Hillenbrand (2008) are around 32 days, which is close to the frequency of the signal at 34.9 days. However, the time series does not exhibit any remaining significant peak after the longterm activity evolution is corrected for (Fig. 21, top panel). The same is true for the bisector velocity span (Fig. 21, middle panel), whose dispersion after detrending is 1.5 m s^{1}. The GLS periodogram of the FWHM does show a peak that could correspond to the rotational period of the star (at period P = 29.5 days; Fig. 21, lower panel). However, its significance is below the 1%level. These facts, together with the general inactive state of the star and the relatively large amplitude of the 34.9day signal detected in the RV lead us to conclude that it is most likely planetary in origin.
The longterm trend is still seen in the residuals of the twoKeplerian model. It is also detected exclusively in the HARPS data. As for HD 1461 and HD 40307, a similar trend is also observable in the HARPS , line bisector, and CCF FWHM (Fig. 20). Once again we are led to conclude that the trend is produced by a change in the activity level of the star in the past ten years. No cyclic behaviour is detected so far, and we therefore decided, as for HD 40307, to model this effect using a thirddegree polynomial. We note that the RV variation in the past ten years is about 10 m s^{1}, which is similar to the dispersion around the oneKeplerian model for the CORALIE and McDonald data. This explains why the activity trend is detected solely by HARPS. Unlike HD 40307, the active and inactive periods of the target are sampled very differently, with only 28 measurements after 2011, the period that could be considered as inactive. This impedes a systematic separate analysis of the inactive and active periods.
All available data (CORALIE98, CORALIE07, McDonald, and HARPS) were employed to constrain the parameters of a twoKeplerian model with an additional thirddegree polynomial, which was employed to account for the longterm effect of the magnetic activity of the star seen in the HARPS data. As for HD 40307, the priors for the coefficients of the polynomial were taken from the leastsquares fit to the time series (Table A.3).
Fig. 19 Top: radial velocities of HD 204313 with the corresponding oneKeplerian model. Red points are HARPS data; blue and orange points represent CORALIE07 and CORALIE98 data, respectively; magenta points are the Robertson et al. (2012) McDonald radial velocities. Middle: generalised LombScargle periodogram of the residuals to the model plotted in the upper panel, exhibiting excess power at P = 34.9 days and a lowfrequency trend. The horizontal dotted and dashed lines represent the 10% and 1% pvalue levels, respectively. Bottom: GLS periodogram of the RV residuals to the oneKeplerian model without the HARPS data. The peak at 34.9 days is no longer significant, and the longterm trend has been replaced by a definite period at P ~ 3500 days with insignificant power. 
Fig. 20 Same as Figs. 2 and 9 for HD 204313, except that the RV data in the upper panel are the residuals to the twoKeplerian model and the red curve in the lower panel is the corresponding periodogram. 
Fig. 21 GLS of the (from top to bottom) , bisector, and FWHM time series of HD 204313 obtained with HARPS after subtracting a thirddegree polynomial to account for the longterm activity evolution. The horizontal linescorrespond to pvalue = 0.1 and 0.01. 
As the is not available for the other instruments, the longterm drift cannot be monitored. Instead, we decided to perform the analysis in two steps equivalent to using a model where the longterm drift is only present for HARPS data: first, we modelled data from all instruments except for HARPS using a twoKeplerian model, without longterm drift, and with a constantjitter model. We used wide, uninformative priors and sampled the parameter posterior distribution with the MCMC algorithm. The MCMC sample was used to estimate the posterior density of the system. The result from this modelling is used as the prior distribution for the analysis of the HARPS data, this time including the longterm thirddegree polynomial and the varyingjitter model. In other words, we updated (in the Bayesian sense) the probability distributions for the parameters present in the twoKeplerian model using HARPS data, and sampled from the posterior of the new parameters associated with the longterm drift. The posterior distribution of the model without longterm trend was approximated by an uncorrelated multinormal distribution for all parameters, except for the eccentricity of the 39.4day candidate, which was described using a Beta distribution with parameters a = 0.93;b = 5.50. The modelled distributions are listed in Table A.3.
The validity of the twostep procedure described above is based on two assumptions: 1) the longterm drift has a negligible effect on data from all instruments other than HARPS; and 2) the posterior distributions of the twoKeplerian model are correctly sampled and modelled. Concerning 1), the fact that no such trend is visible when the HARPS data are left out (see above) partly justifies the assumption. Concerning point 2), we rely on the excellent goodnessoffit found and on the fact that no strong correlations are seen in the posterior distributions of the twoKeplerian model parameters.
We note that even without the HARPS data, both the period and RV semiamplitude of the 34.9day signal are well constrained and the amplitude is significantly different from zero. This seems to contradict our previous statement that the signal would not be detected if HARPS data were not included. However, as discussed by Gregory (2005, Sect. 3.9), it is unwise to conclude on the significance of a signal based on the posterior probability of a given parameter instead of on the computation of the Bayes factor. Indeed, the Bayes factor includes an additional Occam penalisation for the more complex model that is not accounted for otherwise. In this case, when HARPS data are left out, the posterior probability for the oneKeplerian model is orders of magnitude higher than the corresponding twoKeplerian model probability, confirming that the 34.9day signal is not detected without the HARPS data.
Parameter posteriors for the HD 204313 system.
The mean and 68.3% credible intervals of the model parameters are listed in Table 11. The posterior distributions of the amplitudes and eccentricities are presented in Fig. 22 and the phasefolded orbits and the model confidence intervals are shown in Fig. 23. The companion on a 34.9day orbit is a Neptunelike candidate and lies in a nearly circular orbit (e< 0.3 at the 99% level). The parameters agree with those reported previously (Mayor et al. 2011).
The large separation between the two planetary candidates leaves little doubt as to the longterm stability of the system. We nevertheless performed numerical integrations using Mercury (Chambers 1999) for over 7 × 10^{5} yr. We proceeded as for HD 1461 and HD 40307. No secular evolution of the orbital parameters was observed, except for planet c, whose semimajor axis exhibits a fractional change of 19 parts per million, which is of the same order as the energy conservation during the integration.
7.2. Search for the additional signals
Our analysis of the combined CORALIE, HARPS, and McDonald data does not detect the longperiod companion announced by Robertson et al. (2012). Indeed, the longterm trend seen in the residuals of the oneKeplerian model is not constrained from above, unlike the peak seen in the data analysed by Robertson et al. (2012, see their Fig. 2, which shows a peak slowly decreasing for periods longer than around 10 000 days). If we leave out the HARPS data set, which is the only one exhibiting the longterm drift, as mentioned above, a peak appears at a period of around 3500 days, but with insignificant power.
The GLS of the residuals of a model with two Keplerian and a longterm drift exhibit a peak at 4.7 days with an amplitude of 70 cm s^{1}. Although the peak is well below the 0.01 level, we decided to perform a Bayesian comparison between a model with and without a Keplerian at this period. The two estimators of the evidence agree that the proposed signal is not significantly detected: the Bayes factor in favour of the simpler model is 108 ± 37 using the P14 estimator, and 2180 ± 190 for the CJ01 technique. Once again, we see a large difference between the two estimators, and as for HD 40307, the CJ01 estimate seems to punish more complex models more severely. We therefore conclude that the only two significant signals in HD 204313 are those at 2020 and at 34.9 days.
8. Summary and discussion
Fig. 22 Posterior distributions of the amplitude (top row) and eccentricity (bottom row) of the two Keplerian curves used to model the HARPS radial velocities of HD 204313. Because the eccentricity distribution of planet b is much more concentrated that the corresponding distribution of planet c, the vertical scale was set to 20 times larger than the scale shown for planet c. 
Fig. 23 Radial velocity data phasefolded to the bestfit period of each of the two planetary candidates in the HD 204313 system after subtracting the effect of the other signal and the longterm drift. The blue lines represent the 95% highest density interval (HDI). 
We analysed the entire data set produced by the HARPS southern extrasolar planet search programme on HD 1461, HD 40307, and HD 204313, three systems that each contain at least one known planet (Rivera et al. 2010; Mayor et al. 2009; Tuomi et al. 2013a; Ségransan et al. 2010; Robertson et al. 2012). For HD 204313, we also employed CORALIE and McDonald observations to constrain the mean motion of a longperiod companion. The data sets spanmore than ten years and reveal longterm variability associated with changes in the mean activity level of the star on timescales of around a decade.
The model employed contains a term that represents the effect of the evolving activity level of the target stars, as well as a statistical term to account for the shortterm activity effect (the socalled jitter of the star), which is much harder to describe using a deterministic model. The longterm activity signal has been shown to correlate with the evolution of the activity proxy (Lovis et al. 2011a), and we therefore used the HARPS measurement of this activity proxy in time to provide priors for the longterm activity effect. Some assumptions are key to our model: 1) the RV signal produced by the longterm activity variability scales linearly with ; and 2) the stellar jitter is adequately described as an additional Gaussian noise with an amplitude that scales linearly with . We have tested an alternative model with a fixed amplitude for the additional noise, but the data mildly preferred the varyingjitter model.
Potential periodicities were identified using a generalised LombScargle periodogram and by resorting to a genetic algorithm. After we chose the candidate signals, we determined their significance and the total number of Keplerian signals present in each data set using the full machinery of Bayesian model comparison. Estimating the marginal likelihood is problematic (see, for example Gregory 2007). Using different techniques to obtain this estimate (e.g. Chib & Jeliazkov 2001; Perrakis et al. 2014) has permitted us to compare them, and to identify the cases needing more data or a more sophisticated model. The Bayesian model comparison is supposed to perform better than the periodogram analysis that is usually employed by the HARPS team to detect weak signals. One of the reasons involves the uncertainties of the subtracted signals and is discussed, for example, by Lovis et al. (2011b) and Tuomi (2012). A more general reason is discussed, for example, by Sellke et al. (2001) and involves the interpretation of the pvalue obtained for the periodogram peak power as a falsealarm probability. The periodogram is still an essential tool for identifing periodicities in time series, but fails at providing robust estimates of the significance of each signal. In this sense, Bayesian model comparison constitutes a wellestablished method for computing the probability of all involved models, albeit not without technical difficulties, such as estimating the Bayesian evidence. The periodogram analysis and the Bayesian model comparison are expected to produce equally good results for detections with a strong and high signaltonoise ratio, provided that the threshold for the peak power is chosen correctly. A detailed study of the limits of each method is needed to fully understand the techniques used to mine signals in radial velocity data. In the cases analysed here, all signals supported by the Bayesian approach would have also been obtained with the periodogram analysis using the 1% pvalue criterion.
After we established the number of signals, we resorted to other observables obtained from the HARPS spectra (CCF asymmetry and width, and an activity proxy) to conclude on the nature of the RV variations. In general, if the frequency of a RV variation is seen in any of the other observables, it casts doubt on its planetary origin.
Finally, we conclude the following.

HD 1461 hosts two superEarth planetcandidates, with orbital periods P = 5.77 days and P = 13.5 days. The 13.5day signal detected in the HD 1461 system is close to the first harmonic of the rotational period, which could mean that the signal is produced by activity (Boisse et al. 2011). We verified that the signal is coherent over a timescale of several years and that no sign of variability is seen at this frequency in the ancillary observables or in the activity indicator . The minimum masses are Msini = 6.4M_{⊕} and Msini = 5.6M_{⊕} for the 5.77day and 13.5day signals, respectively. The longterm activity signal has an amplitude of 1.5 m s^{1} and a period of 9.6 yr.

HD 40307 hosts four certain planetary companions at periods of between 4.3 and 51.6 days and with masses of between Msini = 3.6M_{⊕} and Msini = 8.7M_{⊕}. We find inconclusive evidence for an additionalcompanion at P ~ 200 days that it would be extremely interesting to confirm, as such a companion would orbit within the habitable zone of the star.

HD 204313 is a system with a Neptunemass planet on a 34.9day orbit and a 4.3M_{J} candidate on an outer orbit (P = 2024 days).
Only six systems with a massive (M> 4M_{J}) companion on an outer (P> 1000 days) orbit have been detected to date. Out of these, HD 204313 has the highest mass ratio between the orbiting companions. The system most closely resembling HD 204313 is arguably HD 38529 (Fischer et al. 2001, 2003; Wright et al. 2009), with a 0.8M_{J} candidate on a 14.3day orbit, and a 12.3M_{J} object on an 2140day orbit. The mass ratio of this system is around five times lower than for HD 204313, however. The system around HD 74156 (Naef et al. 2004) consists of two companions in orbits alike to those of the HD 204313 system (a 1.8 M_{J} object at 51.6 days and an 8.2 M_{J} companion on a ~2500day orbit), but with a much lower mass ratio.
In addition to HD 40307, only a handful of systems with more than three planets has been detected using RVs: μ Ara, 55 Cnc, HD 10180, and GJ876. The continuing monitoring of this target has permitted these detections, and may in the future permit unveiling further planetary signals and improving the modelling of the activity effect.
HD 1461 and HD 40307 are two of only six solartype planetary systems with at least two superEarth companions whose mass is known to better than 50%, and they will become prime targets for the followup mission CHEOPS. The other systems are HD 20794 (Pepe et al. 2011), HD 7924 (Fulton et al. 2015), HD 219134 (Motalebi et al. 2015), and Kepler102 (Marcy et al. 2014). Except for Kepler102, which was first discovered as a transiting candidate system, these systems required more than 100 RV measurements to be detected. This shows the difficulty associated with detecting this type of companions.
Additionally, as the required data sets span many years, activity cycles are omnipresent. Correctly modelling their effect on RV data is therefore necessary to unveil the full system architecture, from the shortperiod companions out to the habitablezone planets. We have shown evidence that periodicities arise at the aliasing frequencies when an incomplete correction of the longterm variability is performed. A fully satisfactory solution is not at hand. In the meantime, adapting the observing strategy based on the activity level of the star seems reasonable. We have shown here that the observations of HD 40307 obtained during the active period of the star contribute little to constraining the planetary system characteristics (see also Fig. 1). Reducing the observing cadence during highactivity periods can save hours of telescope time.
This companion was previously announced by Mayor et al. (2011).
The BIC is a very popular estimator based solely on the maximum likelihood of the model and the number of free parameters (Schwarz et al. 1978). It is therefore extremely simple to compute. According to Kass & Raftery (1995), minus half the BIC tends to the logarithmic evidence of the model as the size of the data set increases. However, the authors warn that that the relative error is , meaning that even for large samples the correct value is not achieved.
Acknowledgments
We acknowledge Ewan Cameron (Oxford) for bringing the Perrakis et al. (2014) estimator to our attention and for his help in understanding the limitations of the TPM estimator. This work has been carried out within the frame of the National Centre for Competence in Research “PlanetS” supported by the Swiss National Science Foundation (SNSF). It made use of the Data Analyses Center for Exoplanets – http://dace.unige.ch – a platform of PlanetS. R.F.D. acknowledges funding from the European Union Seventh Framework Programme (FP7/20072013) under Grant agreement number 313014 (ETAEARTH). S.U., C.L., D.S., F.P., C.M. and A.W. acknowledge the financial support of the SNSF. CM acknowledges the support of the Swiss National Science Foundation under grant BSSIO_155816 “PlanetsInTime”. We thank the University of Geneva for their continuous support to our planet search programmes. M. Gillon is Research Associate at the FNRS. P.F. and N.C.S. acknowledge support by Fundação para a Ciência e a Tecnologia (FCT) through Investigador FCT contracts of reference IF/01037/2013 and IF/00169/2012, respectively, and POPH/FSE (EC) by FEDER funding through the program “Programa Operacional de Factores de Competitividade – COMPETE”. P.F. further acknowledges support from Fundação para a Ciência e a Tecnologia (FCT) in the form of an exploratory project of reference IF/01037/2013CP1191/CT0001. This research has made use of the SIMBAD database and of the VizieR catalogue access tool operated at CDS, France.
References
 AngladaEscudé, G., Tuomi, M., Gerlach, E., et al. 2013, A&A, 556, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AngladaEscudé, G., Arriagada, P., Tuomi, M., et al. 2014, MNRAS, 443, L89 [NASA ADS] [CrossRef] [Google Scholar]
 AngladaEscudé, G., Tuomi, M., Arriagada, P., et al. 2015, ApJ, submitted [arXiv:1506.09072] [Google Scholar]
 Arriagada, P. 2011, ApJ, 734, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Baliunas, S. L., Donahue, R. A., Soon, W. H., et al. 1995, ApJ, 438, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Baluev, R. V. 2013, MNRAS, 436, 807 [NASA ADS] [CrossRef] [Google Scholar]
 Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boisse, I., Bouchy, F., Hébrard, G., et al. 2011, A&A, 528, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boisse, I., Bonfils, X., & Santos, N. C. 2012, A&A, 545, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chib, S., & Jeliazkov, I. 2001, J. Amer. Statis. Assoc., 96, 270 [Google Scholar]
 Cincunegui, C., Díaz, R. F., & Mauas, P. J. D. 2007, A&A, 469, 309 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dawson, R. I., & Fabrycky, D. C. 2010, ApJ, 722, 937 [NASA ADS] [CrossRef] [Google Scholar]
 Delfosse, X., Bonfils, X., Forveille, T., et al. 2013, A&A, 553, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desort, M., Lagrange, A.M., Galland, F., Udry, S., & Mayor, M. 2007, A&A, 473, 983 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Díaz, R. F., Almenara, J. M., Santerne, A., et al. 2014, MNRAS, 441, 983 [NASA ADS] [CrossRef] [Google Scholar]
 Dravins, D. 1982, ARA&A, 20, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Lovis, C., Ségransan, D., et al. 2011a, A&A, 535, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dumusque, X., Udry, S., Lovis, C., Santos, N. C., & Monteiro, M. J. P. F. G. 2011b, A&A, 525, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dumusque, X., Pepe, F., Lovis, C., et al. 2012, Nature, 491, 207 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Pepe, F., Lovis, C., & Latham, D. W. 2015, ApJ, 808, 171 [NASA ADS] [CrossRef] [Google Scholar]
 van Leeuwen, F. 2007, Hipparcos, the New Reduction of the Raw Data, Astrophys. Space Sci. Lib., 350 [Google Scholar]
 Feroz, F., & Hobson, M. P. 2014, MNRAS, 437, 3540 [NASA ADS] [CrossRef] [Google Scholar]
 Figueira, P., Marmier, M., Boué, G., et al. 2012, A&A, 541, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fischer, D. A., Marcy, G. W., Butler, R. P., et al. 2001, ApJ, 551, 1107 [NASA ADS] [CrossRef] [Google Scholar]
 Fischer, D. A., Marcy, G. W., Butler, R. P., et al. 2003, ApJ, 586, 1394 [NASA ADS] [CrossRef] [Google Scholar]
 Friel, N., & Wyse, J. 2012, Statistica Neerlandica, 66, 288 [CrossRef] [Google Scholar]
 Fulton, B. J., Weiss, L. M., Sinukoff, E., et al. 2015, ApJ, 805, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Geweke, J. 1989, Econometrica, 57, 1317 [CrossRef] [Google Scholar]
 Gomes da Silva, J., Santos, N. C., Bonfils, X., et al. 2012, A&A, 541, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gray, D. F., & Baliunas, S. L. 1995, ApJ, 441, 436 [NASA ADS] [CrossRef] [Google Scholar]
 Gray, D. F., Baliunas, S. L., Lockwood, G. W., & Skiff, B. A. 1996, ApJ, 465, 945 [NASA ADS] [CrossRef] [Google Scholar]
 Gregory, P. C. 2005, Bayesian Logical Data Analysis for the Physical Sciences: A Comparative Approach with Mathematica Support (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Gregory, P. C. 2007, MNRAS, 381, 1607 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Gregory, P. C. 2011, MNRAS, 415, 2523 [NASA ADS] [CrossRef] [Google Scholar]
 Gregory, P. C. 2012, ArXiv eprints [arXiv:1212.4058] [Google Scholar]
 Hall, J. C., Lockwood, G. W., & Skiff, B. A. 2007, AJ, 133, 862 [NASA ADS] [CrossRef] [Google Scholar]
 Hall, J. C., Henry, G. W., Lockwood, G. W., Skiff, B. A., & Saar, S. H. 2009, AJ, 138, 312 [NASA ADS] [CrossRef] [Google Scholar]
 Hatzes, A. P. 2002, Astron. Nachr., 323, 392 [NASA ADS] [CrossRef] [Google Scholar]
 Hatzes, A. P. 2013, ApJ, 770, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Henry, T. J., Soderblom, D. R., Donahue, R. A., & Baliunas, S. L. 1996, AJ, 111, 439 [NASA ADS] [CrossRef] [Google Scholar]
 Isaacson, H., & Fischer, D. 2010, ApJ, 725, 875 [NASA ADS] [CrossRef] [Google Scholar]
 Kass, R. E., & Raftery, A. E. 1995, J. Am. Statist. Assn., 90, 773 [Google Scholar]
 Kipping, D. M. 2013, MNRAS, 434, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Lindegren, L., & Dravins, D. 2003, A&A, 401, 1185 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lissauer, J. J., Marcy, G. W., Rowe, J. F., et al. 2012, ApJ, 750, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., Marcy, G. W., Bryson, S. T., et al. 2014, ApJ, 784, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Lovis, C., Dumusque, X., Santos, N. C., et al. 2011a, A&A, unpublished [arXiv:1107.5325] [Google Scholar]
 Lovis, C., Ségransan, D., Mayor, M., et al. 2011b, A&A, 528, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264 [NASA ADS] [CrossRef] [Google Scholar]
 Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Mauas, P. J. D., Buccino, A., Díaz, R., et al. 2012, in IAU Symp. 286, eds. C. H. Mandrini, & D. F. Webb, 317 [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 Mayor, M., Udry, S., Lovis, C., et al. 2009, A&A, 493, 639 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mayor, M., Marmier, M., Lovis, C., et al. 2011, A&A, submitted [arXiv:1109.2497] [Google Scholar]
 Meunier, N., & Lagrange, A.M. 2013, A&A, 551, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meunier, N., Desort, M., & Lagrange, A.M. 2010, A&A, 512, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Motalebi, F., Udry, S., Gillon, M., et al. 2015, A&A, 584, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Naef, D., Mayor, M., Beuzit, J. L., et al. 2004, A&A, 414, 351 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Newton, M. A., & Raftery, A. E. 1994, J. Roy. Statist. Soc., Ser. B (Methodological), 3 [Google Scholar]
 Noyes, R. W., Hartmann, L. W., Baliunas, S. L., Duncan, D. K., & Vaughan, A. H. 1984, ApJ, 279, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Pepe, F., Mayor, M., Galland, F., et al. 2002, A&A, 388, 632 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pepe, F., Lovis, C., Ségransan, D., et al. 2011, A&A, 534, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perrakis, K., Ntzoufras, I., & Tsionas, E. G. 2014, Comp. Stat. Data Analysis, 77, 54 [Google Scholar]
 Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49 [NASA ADS] [Google Scholar]
 Rivera, E. J., Butler, R. P., Vogt, S. S., et al. 2010, ApJ, 708, 1492 [NASA ADS] [CrossRef] [Google Scholar]
 Robert, C. P., & Wraith, D. 2009, in AIP Conf. Ser. 1193, eds. P. M. Goggans, & C.Y. Chan, 251 [Google Scholar]
 Robertson, P., Horner, J., Wittenmyer, R. A., et al. 2012, ApJ, 754, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Robertson, P., Endl, M., Cochran, W. D., & DodsonRobinson, S. E. 2013, ApJ, 764, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Robertson, P., Mahadevan, S., Endl, M., & Roy, A. 2014, Science, 345, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Robertson, P., Roy, A., & Mahadevan, S. 2015, ApJ, 805, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Saar, S. H., & Donahue, R. A. 1997, ApJ, 485, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Santos, N. C., Gomes da Silva, J., Lovis, C., & Melo, C. 2010, A&A, 511, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schwarz, G., et al. 1978, The annals of statistics, 6, 461 [Google Scholar]
 Ségransan, D., Udry, S., Mayor, M., et al. 2010, A&A, 511, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ségransan, D., Mayor, M., Udry, S., et al. 2011, A&A, 535, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sellke, T., Bayarri, M., & Berger, J. O. 2001, Am. Statist., 55, 62 [Google Scholar]
 Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Tuomi, M. 2012, A&A, 543, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tuomi, M., & Jones, H. R. A. 2012, A&A, 544, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tuomi, M., AngladaEscudé, G., Gerlach, E., et al. 2013a, A&A, 549, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tuomi, M., Jones, H. R. A., Jenkins, J. S., et al. 2013b, A&A, 551, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Udry, S., Bonfils, X., Delfosse, X., et al. 2007, A&A, 469, L43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wright, J. T., Fischer, D. A., Ford, E. B., et al. 2009, ApJ, 699, L97 [NASA ADS] [CrossRef] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Parameter priors
Parameter prior distributions for the HD 1461 system.
Prior distributions for the HD 40307 system.
Parameter prior distributions for the HD 204313 system.
All Tables
All Figures
Fig. 1 Stellar jitter for the constantjitter model, fitted to the RV data of HD 40307 before and after BJD = 2 454 800, which have different mean levels of activity (see Sect. 6 for details). 

In the text 
Fig. 2 Top panel: HARPS time series of HD 1461. The vertical dashed line separates the active (BJD < 2 454 850) from the inactive data set. Lower panel: GLS at periods over 400 days for the four time series plotted in the top panel. The vertical dotted lines represent the time span of observations and twice this value. 

In the text 
Fig. 3 Periodograms of the RV data of HD 1461 (top panel) and residuals around models with 1, 2, and 3 Keplerians. The horizontal lines are the 10% and 1% pvalue levels. 

In the text 
Fig. 4 HD 1461. Odds ratio for models with n Keplerian curves with respect to models with n−1 Keplerian curves as a function of model complexity n, assuming equal unity prior odds in all cases. The estimates using different techniques are shown and the customary limits for positive (O_{n + 1,n} = 3) and strong (O_{n + 1,n} = 20) and their inverses are shown as dashed and dotted lines, respectively. The model with four signals contains the 620day Keplerian. 

In the text 
Fig. 5 Marginal posterior distributions of the orbital period, eccentricity, phase, and signal amplitude for the signal at around 620 days for each method used to account for the RV effect of the activity cycle. Also included are the posterior distributions using only the RVs obtained during the period of lower activity. 

In the text 
Fig. 6 HD 1461. GLS periodogram of the RV residuals of the twoKeplerian model (top) and twoKeplerian + linear drift (bottom) for data taken after JD = 2 454 850 (the inactive data set). The two peaks standing out as significant signals in the top panel have periods of 22.9 days and around 650 days. Note that the significance is reduced drastically when the longterm trend caused by the activity cycle is removed, indicating that the observed periodicities are aliases of a longperiod signal present in the data. 

In the text 
Fig. 7 Posterior distributions of the amplitude (top row) and eccentricity (bottom row) of the three Keplerian curves used to model the HARPS radial velocities of HD 1461. The grey dotted curves represent the eccentricity prior for the planetary signals. To facilitate comparison, the axis scales are the same for the three signals. 

In the text 
Fig. 8 Radial velocity data phasefolded to the bestfit period of each of the three Keplerian curves used in the final modelling of HD 1461, after subtracting the effect of the remaining signals. The error bars include the additional noise term. The blue lines represent the 95% highest density interval (HDI). 

In the text 
Fig. 9 Upper panel: HARPS time series of HD 40307. For , the empty circles are the data included in M09, and have ⟨⟩ =−4.99. The filled circles are the new data presented here, with ⟨⟩ =−4.87. The vertical dashed line separates the active (BJD < 2 454 800) from the inactive data set. Lower panel: corresponding GLS periodograms for periods longer than 400 days. The vertical dotted lines represent the time span of observations and twice this value. 

In the text 
Fig. 10 Periodogram of the RV data of HD 40307 (top panel), and of the residuals around models with three (second from top), four (second from bottom), and five (bottom) Keplerian signals in addition to a cubic function to take into account the longterm trend produced by the magnetic cycle of the star. The horizontal dotted and dashed lines represent the 10% and 1% false alarm probability levels, respectively. 

In the text 
Fig. 11 GLS periodogram of the bisector velocity span of HD 40307. A peak at the period of planet d is clearly detected. 

In the text 
Fig. 12 HD 40307. Posterior distribution of the amplitude the additional white noise for an inactive star = – 5.0. The dashed vertical lines represent the mean of each distribution. 

In the text 
Fig. 13 HD 40307. Posterior distribution of the amplitude of the fourth signal. The dashed vertical lines represent the mean of each distribution. The model with four Keplerians produces a slightly weaker signal than the more complex models. 

In the text 
Fig. 14 HD 40307. Odds ratio for models with n Keplerian curves with respect to models with n−1 Keplerian curves as a function of model complexity n, assuming equal unity prior odds in all cases. The estimates using different techniques are shown and the customary limits for positive (O_{n + 1,n} = 3) and strong (O_{n + 1,n} = 20) and their inverses are shown as dashed and dotted lines, respectively. The odds ratio in favour of the model with three Keplerians is out of scale. 

In the text 
Fig. 15 Periodogram of the active FWHM time series of HD 40307 after detrending using a thirddegree polynomial. The single peak reaching thepvalue = 0.01 level (indicated by the dotdashed lines) has P = 37.4 days. We interpret it as a signature of the stellar rotation period. 

In the text 
Fig. 16 Posterior distributions of the amplitude (top row) and eccentricity (bottom row) of the four Keplerian curves used to model the HARPS radial velocities of HD 40307. The grey dotted curves represent the eccentricity prior. To facilitate comparison, the axis scales are the same for all signals. 

In the text 
Fig. 17 Radial velocity data phasefolded to the bestfit period of each of the four Keplerian curves used in the final modelling of HD 40307 after subtracting the effect of the remaining signals and the longterm drift. The error bars include the additional noise term (see text). The inactive data set is plotted using filled red circles, while for the active data set we chose lighter empty circles. This promptly shows that the dispersion around the model is largely caused by the active data set. The blue lines represent the 95% highest density interval (HDI). 

In the text 
Fig. 18 HD 40307. Window function around the oneyear alias (0.00273 cycles/day) for the data set used by T13 (solid curve), consisting of the inactive data set plus a few points taken during commissioning. The corresponding spectral window function for the data set analysed here is shown as a dotted curve. The frequency corresponding to 236 days is indicated by the vertical dashed line. 

In the text 
Fig. 19 Top: radial velocities of HD 204313 with the corresponding oneKeplerian model. Red points are HARPS data; blue and orange points represent CORALIE07 and CORALIE98 data, respectively; magenta points are the Robertson et al. (2012) McDonald radial velocities. Middle: generalised LombScargle periodogram of the residuals to the model plotted in the upper panel, exhibiting excess power at P = 34.9 days and a lowfrequency trend. The horizontal dotted and dashed lines represent the 10% and 1% pvalue levels, respectively. Bottom: GLS periodogram of the RV residuals to the oneKeplerian model without the HARPS data. The peak at 34.9 days is no longer significant, and the longterm trend has been replaced by a definite period at P ~ 3500 days with insignificant power. 

In the text 
Fig. 20 Same as Figs. 2 and 9 for HD 204313, except that the RV data in the upper panel are the residuals to the twoKeplerian model and the red curve in the lower panel is the corresponding periodogram. 

In the text 
Fig. 21 GLS of the (from top to bottom) , bisector, and FWHM time series of HD 204313 obtained with HARPS after subtracting a thirddegree polynomial to account for the longterm activity evolution. The horizontal linescorrespond to pvalue = 0.1 and 0.01. 

In the text 
Fig. 22 Posterior distributions of the amplitude (top row) and eccentricity (bottom row) of the two Keplerian curves used to model the HARPS radial velocities of HD 204313. Because the eccentricity distribution of planet b is much more concentrated that the corresponding distribution of planet c, the vertical scale was set to 20 times larger than the scale shown for planet c. 

In the text 
Fig. 23 Radial velocity data phasefolded to the bestfit period of each of the two planetary candidates in the HD 204313 system after subtracting the effect of the other signal and the longterm drift. The blue lines represent the 95% highest density interval (HDI). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.