# The HARPS search for southern extra-solar planets XXXVIII. Bayesian re-analysis of three systems. New super-Earths, unconfirmed signals, and magnetic cycles

R. F. Diaz, D. Segransan, S. Udry, C. Lovis, F. Pepe, X. Dumusque, M.

Marmier, R. Alonso, W. Benz, F. Bouchy, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

R. F. Diaz, D. Segransan, S. Udry, C. Lovis, F. Pepe, et al.. The HARPS search for southern extra-solar planets XXXVIII. Bayesian re-analysis of three systems. New super-Earths, unconfirmed signals, and magnetic cycles. Astronomy and Astrophysics - A&A, 2016, 585, 10.1051/0004-6361/201526729. hal-01440711

HAL Id: hal-01440711

https://hal.science/hal-01440711

Submitted on 14 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# The HARPS search for southern extra-solar planets

## XXXVIII. Bayesian re-analysis of three systems. New super-Earths, unconfirmed signals, and magnetic cycles\*,\*\*

R. F. Díaz<fnref n="1" />, D. Ségransan<fnref n="1" />, S. Udry<fnref n="1" />, C. Lovis<fnref n="1" />, F. Pepe<fnref n="1" />, X. Dumusque<fnref n="2" />,<fnref n="1" />, M. Marmier<fnref n="1" />, R. Alonso<fnref n="3" />,<fnref n="4" />, W. Benz<fnref n="5" />, F. Bouchy<fnref n="1" />,<fnref n="6" />, A. Coffinet<fnref n="1" />, A. Collier Cameron<fnref n="7" />, M. Deleuil<fnref n="6" />, P. Figueira<fnref n="8" />, M. Gillon<fnref n="9" />, G. Lo Curto<fnref n="10" />, M. Mayor<fnref n="1" />, C. Mordasini<fnref n="5" />, F. Motalebi<fnref n="1" />, C. Moutou<fnref n="6" />,<fnref n="11" />, D. Pollacco<fnref n="12" />, E. Pompei<fnref n="10" />, D. Queloz<fnref n="1" />,<fnref n="13" />, N. Santos<fnref n="8" />,<fnref n="14" />, and A. Wyttenbach<fnref n="1" />

Received 12 June 2015 / Accepted 9 October 2015

## ABSTRACT

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 long-period and low-mass 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 long-term 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 super-Earth 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 habitable-zone candidate HD 40307 g. We also fully characterise the Neptune-mass planet that orbits HD 204313 in 34.9 days.

**Key words.** techniques: radial velocities – methods: data analysis – methods: statistical – planetary systems

## 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 long-period orbits and low-mass 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 multi-planet 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 in-depth exploration of the planetary systems around them. However, even for these very weakly active

<footnote n="1">¹ Observatoire astronomique de l'Université de Genève, 51 Ch. des Maillettes, 1290 Versoix, Switzerland
e-mail: Rodrigo.Diaz@unige.ch</footnote>

<footnote n="2">² Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="3">³ Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="4">⁴ Dpto. de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="5">⁵ Physikalisches Institut, Universität Bern, Silderstrasse 5, 3012 Bern, Switzerland</footnote>

<footnote n="6">⁶ Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, 13388 Marseille, France</footnote>

<footnote n="7">⁷ School of Physics and Astronomy, University of St. Andrews, North Haugh, St. Andrews, Fife KY16 9SS, UK</footnote>

<footnote n="8">⁸ Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="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</footnote>

<footnote n="10">¹⁰ European Southern Observatory, 3107 Alonso de Cordova, Vitacura, Casilla 19001, Santiago 19, Chile</footnote>

<footnote n="11">¹¹ Canada France Hawaii Telescope Corporation, Kamuela, 96743, USA</footnote>

<footnote n="12">¹² Department of Physics, University of Warwick, Coventry, CV4 7AL, UK</footnote>

<footnote n="13">¹³ Cavendish Laboratory, J J Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="14">¹⁴ Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal</footnote>

<footnote>* Based on observations made with the HARPS instrument on the ESO 3.6m telescope at La Silla Observatory under the GTO programme ID 072.C-0488, and its continuation programmes ID 183.C-0972, 091.C-0936, and 192.C-0852.</footnote>

<footnote>** Full Tables 3, 6, and 10 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/585/A134</footnote>
stars, the presence of magnetic cycles complicates the detection of low-mass 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 full-width at half-maximum (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; Anglada-Escudé 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 (Anglada-Escudé 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 short-term effect produced by active regions (spots and plages) that rotate in and out of view as the star revolves, and the long-term 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 short-term 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 short-term 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, long-term 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 log(R'_{HK}) (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 long-term 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 assess their significance. Taking the discussion from the previous paragraph into consideration, we take a different approach to model each of the activity effects: the short-term variability

is described using a non-deterministic model that does not require a detailed description of the activity signal, while the long-term variability is modelled assuming a simple relation between the activity effect on the radial velocities and the activity proxy log (R'_{HK}). 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 long-term 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 high-precision HARPS programmes. The observations were reduced using the HARPS pipeline (version 3.5), and the stellar radial velocities were obtained through a weighted cross-correlation with a numerical mask (Baranne et al. 1996; Pepe et al. 2002). The FWHM and bisector span of the peak in the cross-correlation function (CCF) were also measured for each spectrum, as well as the activity proxy based on the Ca II H and K lines, log (R'_{HK}), 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

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 signal-to-noise 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 semi-major 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
<table><caption>**Table 1.** Basic characteristics of the HARPS observations of the three target stars.</caption><thead><tr><th rowspan="2">Target</th><th rowspan="2">N</th><th rowspan="2">N<sub>nights</sub></th><th colspan="3">Dates</th><th rowspan="2">time span [yr]</th><th rowspan="2">⟨S/N⟩<sub>550 nm</sub></th></tr><tr><th>start</th><th>end</th><th>end</th></tr></thead><tbody><tr><td>HD 1461</td><td>448</td><td>249</td><td>2003-09-16</td><td>2013-11-28</td><td>10.2</td><td>277</td></tr><tr><td>HD 40307</td><td>441</td><td>226</td><td>2003-10-29</td><td>2014-04-05</td><td>10.4</td><td>246</td></tr><tr><td>HD 204313</td><td>96</td><td>95</td><td>2006-05-05</td><td>2014-10-17</td><td>8.5</td><td>151</td></tr></tbody></table>

**Notes.** N is the total number of spectra, and $N_{\text{nights}}$ is the total number of nights on which the target was observed. The average signal-to-noise ratio (⟨S/N⟩) is computed over the nightly average spectra at 550 nm.

<table><caption>**Table 2.** Observed and inferred stellar parameters.</caption><thead><tr><th>Parameters</th><th></th><th>HD 1461</th><th>HD 40307</th><th>HD 204313</th></tr></thead><tbody><tr><td>Sp. T.<sup>(1)</sup></td><td></td><td>G0V</td><td>K3V</td><td>G5V</td></tr><tr><td>V<sup>(1)</sup></td><td></td><td>6.47</td><td>7.17</td><td>7.99</td></tr><tr><td>B - V<sup>(1)</sup></td><td></td><td>0.674</td><td>0.935</td><td>0.697</td></tr><tr><td>π<sup>(2)</sup></td><td>[mas]</td><td>43.02 ± 0.51</td><td>76.95 ± 0.37</td><td>21.11 ± 0.62</td></tr><tr><td>T<sup>(3)</sup><sub>eff</sub></td><td>[K]</td><td>5765 ± 18</td><td>4977 ± 59</td><td>5776 ± 22</td></tr><tr><td>[Fe/H]<sup>(3)</sup></td><td>[dex]</td><td>+0.19 ± 0.01</td><td>-0.31 ± 0.03</td><td>0.18 ± 0.02</td></tr><tr><td>log (g)<sup>(3)</sup></td><td>[cgs]</td><td>4.38 ± 0.03</td><td>4.47 ± 0.16</td><td>4.38 ± 0.02</td></tr><tr><td>M<sup>(3)</sup><sub>*</sub></td><td>[M<sub>☉</sub>]</td><td>1.02</td><td>0.77</td><td>1.02</td></tr><tr><td>log (R'<sub>HK</sub>)<sup>(4)</sup></td><td></td><td>-5.021 ± 0.013</td><td>-4.940 ± 0.058</td><td>-5.024 ± 0.019</td></tr><tr><td>P<sup>(5)</sup><sub>rot</sub></td><td>[days]</td><td>30.2 ± 3.5</td><td>47.9 ± 6.4</td><td>34.1 ± 3.7</td></tr><tr><td>P<sup>(6)</sup><sub>rot</sub></td><td>[days]</td><td>–</td><td>37.4</td><td>–</td></tr><tr><td>v sin (i<sub>*</sub>)</td><td>[km s<sup>-1</sup>]</td><td>–</td><td>1.61</td><td>2.4</td></tr></tbody></table>

**Notes.** (1) As listed in Perryman et al. (1997); (2) F. van Leeuwen (2007); (3) Sousa et al. (2008); (4) this work: mean and standard deviation; (5) Estimated from log(R'_{HK}) using the relationship by Mamajek & Hillenbrand (2008). The reported uncertainty is the quadratic addition of the individual errors propagated from the log(R'_{HK}) uncertainty and the deviation in the values produced by the change in the activity level; (6) Measured in the log(R'_{HK}) FWHM, and/or bisector time series.

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 log (R'_{HK}) 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 solar-like characteristics, these stars have not been systematically included in the southern surveys of stellar activity targeting solar-like 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 log (R'_{HK}) = -5.0. Long-term 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 long-term 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 long-term activity signal constitute the deterministic part of the model (m). To this we add a statistical noise component $\epsilon$ that represents the stellar activity “jitter” – that is, the short-term activity-induced 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 short-term 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 (non-deterministic) component that we call the stellar jitter. Therefore, the RV prediction for model $M_n$ at time $t_i$ can be written as

$$m_i + \epsilon_i = \sum_{j=1}^{n} k_j(t_i) + a(t_i) + \epsilon_i, \quad (1)$$

where $k_j$ is the Keplerian curve of companion j. The Keplerian curves were parametrised using their period, amplitude, eccentricity, argument of periastron $\omega$, 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 $\epsilon_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 $\epsilon$ is distributed as a zero-centred normal, with standard deviation $\sigma_J(t_i) = \sigma_{Ji}$ (potentially a function of time), the likelihood function of our model takes the form (Gregory 2005, Sect. 4.8)

$$ \mathcal{L} = \prod_i \frac{1}{\sqrt{2\pi} \sqrt{\sigma_i^2 + \sigma_{Ji}^2}} \exp\left[-\frac{(v_i - m_i)^2}{2(\sigma_i^2 + \sigma_{Ji}^2)}\right], \quad (2) $$

where $(v_i, \sigma_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 white-noise 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 $\sigma_J$ for all times. In this case, the global value of $\sigma_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 $\sigma_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 $\epsilon$ increases linearly with $\log(R'_{HK})$. The dependency on $\log(R'_{HK})$ is motivated by the fact that the scatter in the RV measurements increases when the $\log(R'_{HK})$ 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 $\log(R'_{HK}) = -5.0$ and the slope of the dependence of the jitter amplitude with $\log(R'_{HK})$. The jitter level when $\log(R'_{HK}) = -5.0$ (the base-level jitter) represents any RV effect that might exist for solar-type 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 long-term 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

<figure><img src="image_3.png" /><figcaption>Fig. 1. Stellar jitter for the constant-jitter 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).</figcaption></figure>

of $a(t)$ is not fixed a priori, but is instead taken from a fit to the log ($R'_{HK}$) time series. We tried a number of different models to fit the log ($R'_{HK}$) time series (see Sect. 5), but in all cases, the “shape” of the log ($R'_{HK}$) time series was used to model the long-term variations seen in the RV data.

To transform the variations in $\log(R'_{HK})$ into variations in RV, a scaling constant $\alpha$ 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 $\alpha$. We instead compared in each case the obtained result with the expected value based on the $T_{eff} - [Fe/H]$ parametrisation.

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 $\log(R'_{HK})$ proxy and the long-term activity-induced 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 $\log(R'_{HK})$.

## 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 non-linear 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 burn-in period and guarantees that the entire parameter space has been explored. We nevertheless launched a number of independent chains to explore
the possibility of multi-modal 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 semi-major 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 RV-measured orbital eccentricities ($a = 0.867, b = 3.03$). The priors for the long-term activity signal were chosen as normal distributions around the least-squares fit to the $\log(R'_{HK})$ 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 $\log(R'_{HK})$ 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 $p$-value as a function of power level. This $p$-value 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 $p$-value of the highest peak in the original histogram is lower than a predefined threshold, the best-fit 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 if the threshold $p$-value is chosen to be low enough and provided the removed signals are well constrained. However, it has two main limitations: a) the interpretation of the $p$-value as a false-alarm 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⁻<fnref n="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 $p$-values 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 $p$-value 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:

$$O_{1,2} = \frac{p(M_1|D, I)}{p(M_2|D, I)} = \frac{p(M_1|I)}{p(M_2|I)} \cdot \frac{p(D|M_1, I)}{p(D|M_2, I)}, \quad (3)$$

where the first term on the right-hand 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 non-nested models. The Bayes factor has a built-in 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|\theta_i, M_i, I) = L(\theta_i)$) over the prior parameter space<fnref n="1" />:

$$\mathcal{E}_i = p(D|M_i, I) = \int \pi(\theta_i|M_i, I) \cdot L(\theta_i) \cdot d\theta_i, \quad (4)$$

where $\theta_i$ denotes the parameter vector of model $i$, and $\pi(\theta_i|M_i, I)$ is the parameter prior distribution. In multi-dimensional parameter spaces, such as those associated with models of multi-planet 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 $L(\theta_i)$ over a sample from the prior density is expected to fail for high-dimensional 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 & Wyse 2012, for a recent review and Kass & Raftery 1995). Estimating the evidence is difficult in multi-dimensional 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 samples and 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 based on the fact that the marginal likelihood is the normalising constant of the posterior density. The method requires

<footnote n="1">¹ The prior distribution was assumed to be normalised to unity.</footnote>
estimating the posterior density at a single point in parameter space $\theta^*$. 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 multi-modal 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.

- 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

$$ \hat{\mathcal{E}}_i = N^{-1} \sum_{n=1}^{N} \frac{\mathcal{L}(\theta^{(n)})\pi(\theta^{(n)}|M_i, I)}{\prod_{j=1}^{q_i} p_j(\theta^{(n)}|D, I)}, \quad (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 $\theta^{(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 time-consuming 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 non-parametric 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 $\mathcal{I}$ is a mixture of posterior distribution samples at different stages of a MCMC:

$$ \mathcal{I} \propto (1 - \lambda) \mathcal{L}(\theta_i)\pi(\theta_i|I) + \lambda \mathcal{L}(\theta_{i-h})\pi(\theta_{i-h}|I). \quad (6) $$

The level of mixture ($\lambda$) and the lag between samples ($h$) are two parameters of the method that the authors explored. The result is called a truncated posterior-mixture (TPM) estimate. This estimator is designed to solve the well-known 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 $\lambda$ tends to zero, and therefore its variance also tends to infinity<fnref n="2" />. However, when $\lambda$ 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 prior-insensitivity 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 $\lambda = 0$, the TPM estimate is equivalent to the problematic HME, and if $\lambda > 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 super-Earth on a 5.77-day 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 \text{ 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 one-year oscillation with an amplitude of ~1.4 m s⁻¹. This one-year 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,

<footnote n="2">² TPM converges in probability to the HME, which implies convergence in distribution (E. Cameron, priv. comm.).</footnote>
<table><caption>**Table 3.** HARPS measurements of HD 1461.</caption><thead><tr><th>BJD<br>-2 450 000</th><th>RV<br>(km s<sup>-1</sup>)</th><th>σ<sub>RV</sub><br>(km s<sup>-1</sup>)</th><th>BIS<br>(m s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>log (R'<sub>HK</sub>)</th><th>σ<sub>log (R'<sub>HK</sub>)</sub></th></tr></thead><tbody><tr><td>52 899.7648</td><td>-10.0528</td><td>0.0011</td><td>-11.32</td><td>7.1703</td><td>-4.9972</td><td>0.0040</td></tr><tr><td>52 900.7412</td><td>-10.0549</td><td>0.0010</td><td>-11.55</td><td>7.1782</td><td>-4.9941</td><td>0.0028</td></tr><tr><td>52 901.7596</td><td>-10.0560</td><td>0.0009</td><td>-12.80</td><td>7.1792</td><td>-4.9998</td><td>0.0021</td></tr><tr><td>52 902.7191</td><td>-10.0585</td><td>0.0010</td><td>-13.28</td><td>7.1767</td><td>-4.9999</td><td>0.0026</td></tr><tr><td>52 903.7775</td><td>-10.0566</td><td>0.0010</td><td>-14.92</td><td>7.1736</td><td>-5.0105</td><td>0.0035</td></tr></tbody></table>

**Notes.** The full version is available at the CDS.

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 log(R'_{HK}) activity proxy and two spectral line measurements (FWHM and bisector velocity span) that can also be affected by activity. A similar long-term 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 Lomb-Scargle 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⁻¹ also agrees with Rivera et al. (2010) and corresponds to a minimum mass of around 6.4 M<sub>⊕</sub>. 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 three-Keplerian model. The model probabilities are plotted in Fig. 4.

## 5.1. Two-Keplerian model. A new super-Earth candidate

The residuals around the one-Keplerian model (Fig. 3) reveal a new significant signal with *P* = 13.5 days and an amplitude of 1.5 m s⁻¹, 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 long-term 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.5-day peak is an alias of the longer activity-induced 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 activity-induced 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 log(R'_{HK}), 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⁻¹ and 1.16 m s⁻¹, respectively, before and after the degree-three polynomial is used to correct for the effect of the activity cycle

<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

(see below). The GLS of the bisector exhibits significant power at 29.2 days with an amplitude of around 60 cm s⁻¹, most likely caused by the stellar rotational modulation (see Table 2). The FWHM time series does not present any significant periodicity when the long-term trend is removed and exhibits a dispersion of only 3 m s⁻¹ over more than ten years. The time series of the log (R'_{HK}) activity proxy, after correction for the long-term 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.
<figure><img src="image_5.png" /><figcaption>**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% *p*-value levels.</figcaption></figure>

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_\odot^3$. The parameters of the new companion are listed in Table 5.

## 5.2. Activity cycle and search for additional signals

A common long-term 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 solar-like 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 long-period signal seen in the RV time series is produced by a magnetic activity cycle, with a period of $P_{cycle} = 9.64 \pm 0.21$ yr, as measured by

a Keplerian fit to the log ($R'_{HK}$) time series. Hall et al. (2009) presented seasonally averaged log ($R'_{HK}$) 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 \pm 0.4$ yr, an amplitude above $3 \ \text{m s}^{-1}$, and a significant eccentricity of $e = 0.43 \pm 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 long-term variation with real periodic signals. For example, the one-year 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 best-fit Keplerian curve to the log ($R'_{HK}$) has a slightly different period and eccentricity ($e = 0.17 \pm 0.07$ for log ($R'_{HK}$)) 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.

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 log ($R'_{HK}$) variations are modelled using a sinusoidal function, and the best-fit parameters are used as Gaussian priors for a fit to the RV time series, with the exception of the 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 low-pass filter (cutoff at 100 days) to the log ($R'_{HK}$) 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 log ($R'_{HK}$) 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 log ($R'_{HK}$).

We note that all the methods used to account for the activity cycle assume a linear relation between the variations observed in the log ($R'_{HK}$) 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 log ($R'_{HK}$) 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 = 2454850, which correspond to the last five observing seasons and to the period of lesser magnetic activity, according to the log ($R'_{HK}$) 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

<footnote n="3">³ This companion was previously announced by Mayor et al. (2011).</footnote>
<table><caption>**Table 4.** Model probabilities for HD 1461.</caption><thead><tr><th rowspan="2">n</th><th rowspan="2">Periods [d]</th><th colspan="4">log p(D|M<sub>n</sub>, I) - 1100</th><th rowspan="2">p(M<sub>n</sub>|D, I)/p(M<sub>3</sub>|D, I)</th><th rowspan="2">σ<sub>J<sub>n</sub></sub>|log R'<sub>HK</sub>|<sub>=5</sub> [m s<sup>-1</sup>]</th><th rowspan="2">σ<sub>O-C</sub> [m s<sup>-1</sup>]</th></tr><tr><th>CJ01</th><th>P14</th><th>TPM</th><th>BIC</th></tr></thead><tbody><tr><td>3</td><td>{5.77, 13.5, 3500}</td><td>80.44 ± 0.06</td><td>81.15 ± 0.05</td><td>132.51</td><td>93.76</td><td>1.0</td><td>1.0</td><td>2.20 ± 0.19</td><td>2.257 ± 0.071</td></tr><tr><td>4</td><td>{5.77, 13.5, 22.9, 3500}</td><td>76.31 ± 0.11</td><td>79.84 ± 0.05</td><td>143.52</td><td>93.18</td><td>0.016 ± 0.002</td><td>0.274 ± 0.016</td><td>2.15 ± 0.18</td><td>2.242 ± 0.075</td></tr><tr><td>4</td><td>{5.77, 13.5, 620, 3500}</td><td>84.57 ± 0.04</td><td>84.16 ± 0.11</td><td>145.74</td><td>97.07</td><td>60.8 ± 4.0</td><td>19.6 ± 2.4</td><td>1.87 ± 0.19</td><td>1.965 ± 0.092</td></tr><tr><td>5</td><td>{5.77, 13.5, 22.9, 620, 3500}</td><td>77.24 ± 0.17</td><td>76.92 ± 0.14</td><td>153.97</td><td>93.94</td><td>(2.25 ± 0.30) × 10<sup>-4</sup></td><td>0.015 ± 0.003</td><td>1.81 ± 0.19</td><td>1.946 ± 0.098</td></tr></tbody></table>

**Notes.** Estimate of the evidence (marginal likelihood) for models with *n* = 3, ..., 5 Keplerians. The periods of the included signals are listed in the second column. The estimates based on the method of Chib & Jeliazkov (2001, CJ01) and Perrakis et al. (2014, P14) are given in Cols. 3 and 4. For comparison, we also report the estimate obtained based on the Tuomi & Jones (2012) TPM estimator and on the Bayesian information criterion (BIC) in Cols. 5 and 6. For clarity, we have subtracted 1100 from each estimation, which corresponds to a change of units in the observed radial velocities. Columns 7 and 8 list the odds ratio between model *i* with respect to the three-Keplerian model. Additionally, the posterior estimates of the amplitude of the additional base-level noise and the dispersion of the residuals are given in Cols. 9 and 10.

<table><caption>**Table 5.** Parameter posteriors for the HD 1461 system.</caption><thead><tr><th>Orbital parameters</th><th></th><th>Planet b</th><th>Planet c</th><th>Magnetic cycle</th></tr></thead><tbody><tr><td>Orbital period, P<sup>•</sup></td><td>[days]</td><td>5.77152 ± 0.00045</td><td>13.5052 ± 0.0029</td><td>3503 ± 80</td></tr><tr><td>RV amplitude, K<sup>•</sup></td><td>[m s<sup>-1</sup>]</td><td>2.28 ± 0.15</td><td>1.49 ± 0.17</td><td>1.51 ± 0.26</td></tr><tr><td>Eccentricity, e<sup>•</sup></td><td></td><td>&lt;0.131; &lt;0.172<sup>†</sup></td><td>&lt;0.228; &lt;0.305<sup>†</sup></td><td>0.103<sup>+0.063</sup><sub>-0.030</sub></td></tr><tr><td>Argument of periastron, ω<sup>•</sup></td><td>[deg]</td><td>−<sup>†</sup></td><td>−<sup>†</sup></td><td>294 ± 15</td></tr><tr><td>e<sup>1/2</sup> cos(ω)</td><td></td><td></td><td></td><td>0.131 ± 0.09<sup>•</sup></td></tr><tr><td>e<sup>1/2</sup> sin(ω)</td><td></td><td></td><td></td><td>-0.298 ± 0.068<sup>•</sup></td></tr><tr><td>Mean longitude at epoch, L<sub>0</sub><sup>•</sup></td><td>[deg]</td><td>271.6 ± 4.1</td><td>317.9 ± 6.6</td><td>148.6 ± 6.7</td></tr><tr><td>Systemic velocity, V<sub>0</sub><sup>•</sup></td><td>[km s<sup>-1</sup>]</td><td></td><td>-10.05960 ± 1.7 × 10<sup>-4</sup></td><td></td></tr><tr><td>Semi-major axis of relative orbit, a</td><td>[AU]</td><td>0.0634 ± 0.0022</td><td>0.1117 ± 0.0039</td><td>—</td></tr><tr><td>Minimum mass, M sin i</td><td>[M<sub>⊕</sub>]</td><td>6.44 ± 0.61</td><td>5.59 ± 0.73</td><td>—</td></tr><tr><td>Noise model<sup>‡</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>Additional noise at log (R'<sub>HK</sub>) = -5, σ<sub>J</sub>|<sub>-5.0</sub><sup>•</sup></td><td>[m s<sup>-1</sup>]</td><td></td><td>2.13<sup>+0.25</sup><sub>-0.11</sub></td><td></td></tr><tr><td>Slope, α<sub>J</sub><sup>•</sup></td><td>[m s<sup>-1</sup>/dex]</td><td></td><td>25.5 ± 6.4</td><td></td></tr><tr><td>Additional noise at log (R'<sub>HK</sub>)</td><td>[m s<sup>-1</sup>]</td><td></td><td>1.67 ± 0.09</td><td></td></tr><tr><td>rms(O - C)</td><td>[m s<sup>-1</sup>]</td><td></td><td>2.260 ± 0.067</td><td></td></tr></tbody></table>

**Notes.** The epoch is BJD = 2455155.3854 for planets b and c and BJD = 2455195.8367 for the magnetic cycle. (•) MCMC jump parameter. (†) Eccentricity does not differ significantly from zero; the 95% and 99% upper limits are reported. The argument of periastron ω is therefore unconstrained. (‡) The additional (stellar) noise for measurement *i* is σ<sub>J<sub>i</sub></sub> = σ<sub>J</sub>|<sub>-5.0</sub> + α<sub>J</sub> · (log (R'<sub>HK</sub>)<sub>i</sub> + 5.0).

Keplerian signal to the model with two planets and the magnetic cycle. We initialised the MCMC algorithm using the best-fit solution of the two-Keplerian model for the parameters of the two super-Earths 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(Lπ) in each maximum, where L 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<sup>-1</sup>, and

another one around 620 days with an amplitude of 1.2 m s<sup>-1</sup>. We note that their frequencies are not significantly present in the periodogram of the residuals of the three-Keplerian 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 log (R'<sub>HK</sub>) time series, as explained above.

### 5.3. Four-Keplerian model I. The 22.9-day signal

The parameters of the 22.9-day 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.9-day signal. The Bayesian information criterion (BIC)<sup>4</sup> is inconclusive in this respect. We therefore discarded the possibility that only the 22.9-day signal is present. We tested, on

<footnote n="4"><sup>4</sup> 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</footnote>
<figure><img src="image_6.png" /><figcaption>**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*<sub>*n*+1,*n*</sub> = 3) and strong (*O*<sub>*n*+1,*n*</sub> = 20) and their inverses are shown as dashed and dotted lines, respectively. The model with four signals contains the 620-day Keplerian.</figcaption></figure>

the other hand, the inclusion of both 22.9-day and 620-day sig-
nals, 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 pref-
ference for more complex models, as is clearly seen in Table 4.

5.4. Four-Keplerian model II. The 620-day period signal

Both the CJ01 and P14 techniques favour the inclusion of a sig-
nal 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 activ-
ity 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<sub>d</sub> ~ 14.5 ± 1.3 M<sub>⊕</sub>. No sig-
nificant power is present at similar periods in the time series of
the log (R'<sub>HK</sub>), the bisector velocity span, or the FWHM, even af-
ter subtracting the long-term trend associated with the magnetic
cycle.

If only the inactive data set (BJD > 2 454 850) is consid-
ered, the GLS periodogram of the residuals of a two-Keplerian
model exhibits significant (*p*-value <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 am-
plitude of the peak is strongly reduced (Fig. 6), indicating that
the periodicity may be an alias of the long-term trend and not
a real signal. The 22.9-day signal exhibits the same behaviour.
This would explain why the signal is recovered for all the correc-
tion methods of the activity cycle, as well as the slight change of
the period under different corrections. Since a long-period 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
long-term signal, still not fully sampled, is present in the data.
We conclude that although the periodicity at 620-day period is

evidence of the model as the size of the data set increases. However, the
authors warn that that the relative error is *O*(1), meaning that even for
large samples the correct value is not achieved.

<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

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 ad-
ditional Keplerian curve to model the activity cycle. The planet
signals are independent of the method used to model the activity
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** HD 1461. GLS periodogram of the RV residuals of the two-Keplerian model (top) and two-Keplerian + 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 long-term trend caused by the activity cycle is removed, indicating that the observed periodicities are aliases of a long-period signal present in the data.</figcaption></figure>

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 semi-amplitude 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 semi-amplitudes and the eccentricities is much 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.5-day period has an amplitude of $1.49 \pm 0.17$ m s$^{-1}$, which implies a minimum mass of $5.59 \pm 0.73$ M$_\oplus$, in agreement with the values reported by Mayor et al. (2011). The RV amplitude associated with the activity cycle is $1.51 \pm 0.26$ m s$^{-1}$, which means that the scaling constant $\alpha$ between $\log(R'_{HK})$ and RV is $74.2 \pm 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 best-fit 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 semi-major axis change is smaller than $10^{-6}$ for the outer companion and around $8 \times 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 \pm 0.1$ m s$^{-1}$ for the mean $\log(R'_{HK})$ 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.

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 super-Earth-type 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 \sim 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 cross-correlation – 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 moving-average model to take into account the correlation between individual observations taken during a single night, a deterministic model of the short-term 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.

The radial velocity data set is plotted in Fig. 9 together with the time series of the log ($R'_{HK}$), the bisector velocity span, and the FWHM of the CCF. A low-frequency signal is clearly visible in all the observables. As for HD 1461, it seems reasonable to assume that this long-term trend is linked to a stellar magnetic cycle. This signal is currently stronger than the reflex motion produced by the known companion at a four-day period, which illustrates the hindering effect magnetic cycles have on the detection of low-mass planets. The period of the signal is largely unconstrained, and we therefore decided to model it with a third-degree polynomial instead of a periodic function.
<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Radial velocity data phase-folded to the best-fit 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).</figcaption></figure>

<table><caption>**Table 6.** HARPS measurements of HD 40307.</caption><thead><tr><th>BJD<br>-2 450 000</th><th>RV<br>(km s<sup>-1</sup>)</th><th>&sigma;<sub>RV</sub><br>(km s<sup>-1</sup>)</th><th>BIS<br>(m s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>log (R'<sub>HK</sub>)</th><th>&sigma;<sub>log (R'<sub>HK</sub>)</sub></th></tr></thead><tbody><tr><td>52 942.8215</td><td>31.3372</td><td>0.0009</td><td>8.98</td><td>5.9072</td><td>-4.9460</td><td>0.0029</td></tr><tr><td>52 999.7639</td><td>31.3340</td><td>0.0011</td><td>8.84</td><td>5.9011</td><td>-4.9971</td><td>0.0045</td></tr><tr><td>53 000.7606</td><td>31.3375</td><td>0.0012</td><td>11.18</td><td>5.9002</td><td>-4.9849</td><td>0.0054</td></tr><tr><td>53 001.6684</td><td>31.3397</td><td>0.0012</td><td>8.41</td><td>5.8989</td><td>-4.9859</td><td>0.0054</td></tr><tr><td>53 002.6686</td><td>31.3393</td><td>0.0011</td><td>8.40</td><td>5.8993</td><td>-4.9689</td><td>0.0034</td></tr></tbody></table>

**Notes.** The full version is available at the CDS.

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 $\langle \log(R'_{HK}) \rangle = -4.99$, and a typical dispersion of 0.014. In the active data set we find $\langle \log(R'_{HK}) \rangle = -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⁻¹ to 1.7 m s⁻¹ between the inactive and active data sets and justifies the use of the varying-jitter model.

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 long-term trend is also present. The remaining panels in Fig. 10 present the GLS of the residuals to fits with three, four, and five Keplerian signals, plus an additional third-degree polynomial to account for the activity cycle. The corresponding model evidence estimates 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 \sim 20.42$ d instead of $P \sim 20.46$ d). This is
<figure><img src="image_11.png" /></figure>

**Fig. 9.** Upper panel: HARPS time series of HD 40307. For $log(R'_{HK})$, the empty circles are the data included in M09, and have $log(R'_{HK}) = -4.99$. The filled circles are the new data presented here, with $log(R'_{HK}) = -4.87$. The vertical dashed line separates the active (BJD < 2454 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.

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 long-term trend is corrected. However, when a least-squares fit is performed on each observing season individually, the amplitude of the bisector signal is seen to anti-correlate with the one in the RVs. The bisector amplitude varies from below $50 \text{ cm s}^{-1}$ during the first three seasons to around $2.5 \text{ 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 anti-correlation 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 base-level additional noise is below $1 \text{ 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 six-signal models, the noise level is around $60 \text{ cm s}^{-1}$. On the other hand, the models with a weaker base-level jitter have a higher sensitivity to $log(R'_{HK})$, that is, a larger slope parameter.

<figure><img src="image_12.png" /></figure>

**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 long-term 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.

## 6.1. Four-Keplerian model. A super-Earth companion on a 51.6-day period orbit

A 51.6-day period signal appears as significant when the model with three Keplerian signals and a degree-three polynomial is subtracted (Fig. 10). Given that a long-term 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.6-day signal has an amplitude of $75 \text{ cm s}^{-1}$. All models with at least four Keplerians converge to a period of $P \sim 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 \text{ cm s}^{-1}$, while in more complex models the amplitude is closer to $85–90 \text{ cm s}^{-1}$.
<table><caption>**Table 7.** Model probabilities for HD 40307.</caption><thead><tr><th rowspan="2">n</th><th rowspan="2">Periods [d]</th><th colspan="4">log p(D|M<sub>n</sub>, I) - 1000</th><th colspan="2">p(M<sub>n</sub>|D, I)/p(M<sub>3</sub>|D, I)</th><th rowspan="2">σ<sub>J<sub>n</sub></sub>|log R'<sub>HK</sub><sup>-5</sup> [m s<sup>-1</sup>]</th><th rowspan="2">σ<sub>O-C</sub><sup>(inactive)</sup> [m s<sup>-1</sup>]</th></tr><tr><th>CJ01</th><th>P14</th><th>TPM</th><th>BIC</th><th>CJ01</th><th>P14</th></tr></thead><tbody><tr><td>2</td><td>{9.6, 20.4}</td><td>1.95 ± 0.09</td><td>0.57 ± 0.04</td><td>52.41</td><td>16.60</td><td>(1.93 ± 0.23) × 10<sup>-21</sup></td><td>(1.65 ± 0.12) × 10<sup>-21</sup></td><td>1.70 ± 0.15</td><td>1.81 ± 0.04</td></tr><tr><td>3</td><td>{4.3, 9.6, 20.4}</td><td>49.65 ± 0.10</td><td>48.43 ± 0.07</td><td>123.90</td><td>76.73</td><td>1.0</td><td>1.0</td><td>0.972 ± 0.099</td><td>1.19 ± 0.03</td></tr><tr><td>4</td><td>{4.3, 9.6, 20.4, 51.6}</td><td>51.88 ± 0.16</td><td>51.54 ± 0.12</td><td>143.41</td><td>84.76</td><td>9.3 ± 1.7</td><td>22.6 ± 2.8</td><td>0.818 ± 0.090</td><td></td></tr><tr><td>5</td><td>{4.3, 9.6, 20.4, 51.6, 204.9}</td><td>46.46<sup>+0.42</sup><sub>-0.72</sub></td><td>54.30<sup>+0.19</sup><sub>-0.28</sub></td><td>164.15</td><td>99.91</td><td>0.041 ± 0.026</td><td>356 ± 11</td><td>0.630 ± 0.079</td><td>0.92 ± 0.03</td></tr><tr><td>6</td><td>{4.3, 9.6, 20.4, 51.6, ...<sup>†</sup>}</td><td>8.7 ± 7.0</td><td>47.7 ± 2.3</td><td>170.43</td><td>93.44</td><td>1.64 × 10<sup>-18</sup></td><td>0.46<sup>+4.2</sup><sub>-0.41</sub></td><td>0.574 ± 0.084</td><td>0.85 ± 0.04</td></tr></tbody></table>

**Notes.** Estimate of the evidence (marginal likelihood) for models with *n* = 2, ..., 5 Keplerians and an additional long-term drift, common to all of them. Columns are the same as in Table 4. (†) the orbital periods of the fifth and sixth Keplerian curves are not well constrained in this model. See text for details.

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** GLS periodogram of the bisector velocity span of HD 40307. A peak at the period of planet d is clearly detected.</figcaption></figure>

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 high-dimensional spaces (Gregory 2007).

We note that this signal is not far from the rotational period estimated based on the log ($R'_{HK}$) level (Table 2). Indeed, the active period of the bisector velocity span exhibits a significant peak ($p$-value < 0.01) at $P = 51.5$ days. However, the peak power is reduced to below the level of the $p$-value = 0.1 when the data are detrended to account for the long-term evolution. Additionally, no equivalent peak is seen in the inactive period. Although this may cause concern at first sight, the 51.6-day signal is present in the RV data set even after the long-term effect has been removed. The fourth Keplerian signal is therefore probably not attributable to stellar activity. The parameters are presented in Table 9.

## 6.2. Five-Keplerian model. A doubtful 205-day period signal

The periodogram of the residuals around the four-Keplerian model exhibits peaks at ~200 days and 28.6 days with power above the $p$-value = 0.01 level (Fig. 10). The signal at 28.6 days is most probably an alias introduced by an incomplete correction of the long-term variability. The most prominent peak of the window function after the one-year peak is at 27.8 days, probably

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** HD 40307. Posterior distribution of the amplitude the additional white noise for an inactive star log ($R'_{HK}$) = -5.0. The dashed vertical lines represent the mean of each distribution.</figcaption></figure>

linked to the Moon orbital period<fnref n="5" />. Including a fifth Keplerian in the MCMC model leads to a period of $P \sim 205$ d for the fifth signal, with an amplitude of 80 cm s⁻¹ 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 five-Keplerian model is around 15.8 times more probable than the model with only four signals, the CJ01 estimate indicates that the four-Keplerian model is 2.6 times more probable than the five-Keplerian 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

<footnote n="5">⁵ Observations are usually scarcer during full Moon.</footnote>
<figure><img src="image_15.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**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*<sub>n+1,n</sub> = 3) and strong (*O*<sub>n+1,n</sub> = 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.</figcaption></figure>

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 four-Keplerian 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

<table><caption>**Table 8.** Summary of the Bayesian evidence computation for each model.</caption><thead><tr><th>n</th><th>log p(D|M<sub>n</sub>, I) - 1000</th><th>log<sub>10</sub> p(M<sub>n</sub>|D, I)/p(M<sub>3</sub>|D, I)</th></tr></thead><tbody><tr><td>2</td><td>1.26 ± 0.70</td><td>-20.75 ± 0.61</td></tr><tr><td>3</td><td>49.04 ± 0.62</td><td>0.0</td></tr><tr><td>4</td><td>51.71 ± 0.23</td><td>1.16 ± 0.29</td></tr><tr><td>5</td><td>50.38 ± 3.98</td><td>0.6 ± 1.8</td></tr><tr><td>6</td><td>28.2 ± 20.7</td><td>-9.0 ± 9.0</td></tr></tbody></table>

**Notes.** The geometric mean between the estimates by CJ01 and P14 is computed. The uncertainties are the quadratic mean between the statistical errors reported in Table 7 and the difference between the two estimates.

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)<fnref n="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 205-day variability.

### 6.3. Six-Keplerian 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 six-signal 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 four-Keplerian one, but the conclusion holds even if the systematic error between the estimators is taken into account (Table 8).

### 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

<footnote n="6">⁶ There is significant power at 180 days and 240 days, the latter probably related to the duration of the observing season.</footnote>
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Periodogram of the active FWHM time series of HD 40307 after detrending using a third-degree polynomial. The single peak reaching the *p*-value = 0.01 level (indicated by the dot-dashed lines) has *P* = 37.4 days. We interpret it as a signature of the stellar rotation period.</figcaption></figure>

is 1.25 m s⁻¹ over more than ten years and only 0.86 m s⁻¹ for the inactive data (defined in Sect. 6; around three years of data). This shows the exquisite precision and long-term 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 long-term evolution using a second-degree polynomial fit (Fig. 15). A peak at the same period appears in the GLS periodogram of the active log (*R*′<sub>HK</sub>) time series, albeit with *p*-value > 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 log (*R*′<sub>HK</sub>) is considered.

Other signals are present in the log (R'_{HK}) 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 long-term 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 log (R'_{HK}) and FWHM time series as well. Most, if not all of them, are probably caused by an imperfect detrending of the long-term evolution, which introduces periodicities at frequencies associated with the spectral window function.

## 6.5. The planetary system around HD 40307

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 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 super-Earths. 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 phase-folded velocity variation for each planet candidate after subtracting the effect of the remaining ones and the long-term 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 semi-major 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 –although the 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 six-Keplerian 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 log (*R*′<sub>HK</sub>), bi-sector, 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 ~236-day peak is much more prominent than in our data set, probably because of the longer observation time-span (Fig. 18). The 236-day alias of *P* = 41.9 days (the rotational period present in the inactive log (*R*′<sub>HK</sub>) 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 long-term drift attributed to stellar activity, the posterior distribution of the polynomial coefficients resemble closely the priors imposed from the analysis of the log (R'_{HK}) time series (Table A.1). The conversion constant between log (R'_{HK}) and RV is found to be (22.2 ± 2.2) m s⁻¹/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⁻¹ per unit of log (R'_{HK}), but where the error in the calibration coefficients has not been considered. This could indicate the presence of long-period 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.

## 7. HD 204313

Ségransan et al. (2010) discovered a companion to HD 204313 with $m \sin i = 4 M_⊕$ on a 5-yr 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 Neptune-mass planet on an inner orbit of 34.9 days. More recently, Robertson et al. (2012)
<table><caption>**Table 9.** Parameter posteriors for the HD 40307 system.</caption><thead><tr><th colspan="2">Orbital parameters</th><th>Planet b</th><th>Planet c</th></tr></thead><tbody><tr><td colspan="2">Orbital period, P<sup>•</sup></td><td>[days]</td><td>4.31150 ± 0.00027</td><td>9.62070 ± 0.0012</td></tr><tr><td colspan="2">RV amplitude, K<sup>•</sup></td><td>[m s<sup>-1</sup>]</td><td>1.79 ± 0.13</td><td>2.31 ± 0.13</td></tr><tr><td colspan="2">Eccentricity, e<sup>•</sup></td><td></td><td>&lt;0.161; &lt;0.168<sup>†</sup></td><td>&lt;0.098; &lt;0.103<sup>†</sup></td></tr><tr><td colspan="2">Mean longitude at epoch, L<sub>0</sub><sup>•</sup></td><td>[deg]</td><td>163.5 ± 4.0</td><td>321.5 ± 3.2</td></tr><tr><td colspan="2">Semi-major axis of relative orbit, a</td><td>[AU]</td><td>0.0475 ± 0.0011</td><td>0.0812 ± 0.0018</td></tr><tr><td colspan="2">Minimum mass, M sin i</td><td>[M<sub>⊕</sub>]</td><td>3.81 ± 0.30</td><td>6.43 ± 0.44</td></tr><tr><th colspan="2">Orbital parameters</th><th>Planet d</th><th>Planet f</th></tr><tr><td colspan="2">Orbital period, P<sup>•</sup></td><td>[days]</td><td>20.4185 ± 0.0052</td><td>51.56 ± 0.14</td></tr><tr><td colspan="2">RV amplitude, K<sup>•</sup></td><td>[m s<sup>-1</sup>]</td><td>2.44 ± 0.13</td><td>0.75 ± 0.13</td></tr><tr><td colspan="2">Eccentricity, e<sup>•</sup></td><td></td><td>&lt;0.117; &lt;0.122<sup>†</sup></td><td>&lt;0.335; &lt;0.352<sup>†</sup></td></tr><tr><td colspan="2">Mean longitude at epoch, L<sub>0</sub><sup>•</sup></td><td>[deg]</td><td>185.4 ± 3.1</td><td>142 ± 11</td></tr><tr><td colspan="2">Semi-major axis of relative orbit, a</td><td>[AU]</td><td>0.1340 ± 0.0029</td><td>0.2485 ± 0.0054</td></tr><tr><td colspan="2">Minimum mass, M sin i</td><td>[M<sub>⊕</sub>]</td><td>8.74 ± 0.58</td><td>3.63 ± 0.60</td></tr><tr><th colspan="2">Velocity drift<sup>•</sup></th><td></td><td></td><td></td></tr><tr><td colspan="2">Systemic velocity, V<sub>0</sub><sup>•</sup></td><td>[km s<sup>-1</sup>]</td><td>31.334376 ± 9.9 × 10<sup>-5</sup></td><td></td></tr><tr><td colspan="2">Scaling constant, α<sup>•</sup></td><td>[m s<sup>-1</sup>/dex]</td><td>22.1 ± 2.1</td><td></td></tr><tr><td colspan="2">Linear<sup>•</sup></td><td>[dex yr<sup>-1</sup>]</td><td>0.03761 ± 8.9 × 10<sup>-4</sup></td><td></td></tr><tr><td colspan="2">Quadratic<sup>•</sup></td><td>[dex yr<sup>-2</sup>]</td><td>0.00488 ± 4.2 × 10<sup>-4</sup></td><td></td></tr><tr><td colspan="2">Cubic<sup>•</sup></td><td>[dex yr<sup>-3</sup>]</td><td>-0.001617 ± 8.5 × 10<sup>-4</sup></td><td></td></tr><tr><th colspan="2">Noise model<sup>‡</sup></th><td></td><td></td><td></td></tr><tr><td colspan="2">Additional noise at log (R'<sub>HK</sub>)=-5, σ<sub>J</sub>|<sub>-5.0</sub><sup>•</sup></td><td>[cm s<sup>-1</sup>]</td><td>81.8 ± 8.5</td><td></td></tr><tr><td colspan="2">Slope, α<sub>J</sub><sup>•</sup></td><td>[m s<sup>-1</sup>/dex]</td><td>8.1 ± 1.6</td><td></td></tr><tr><td colspan="2">rms(O - C)</td><td>[m s<sup>-1</sup>]</td><td>1.733 ± 0.022</td><td></td></tr></tbody></table>

Notes. The epoch is BJD = 2 454 521.6791. The argument of periastron $\omega$ is unconstrained for these nearly-circular orbits. (*) MCMC jump parameter. (†) Eccentricity does not differ significantly from zero; the 95% and 99% upper limits are reported. (•) See Sect. 4.3. (‡) The additional (stellar) noise for measurement *i* is $\sigma_{Ji} = \sigma_{J|_{-5.0}} + \alpha_J \cdot (\log(R'_{HK})_i + 5.0)$.

<figure><img src="image_18.png" /><figcaption>**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.</figcaption></figure>

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.7-MJ 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⁻¹. Additionally, the hardware changes (see Ségransan et al. 2010 for details) introduced a zero-point offset of around 20 m s⁻¹. For this reason, the data taken before and
<figure><img src="image_19.png" /></figure>

**Fig. 17.** Radial velocity data phase-folded to the best-fit 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 long-term 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).

<figure><img src="image_20.png" /></figure>

**Fig. 18.** HD 40307. Window function around the one-year 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.

after the upgrade were modelled as if they had been taken by different instruments<fnref n="7" />. The data span 13 yr and contain 215 points. The weighted mean error of the measurements is 63 cm s⁻¹. Two HARPS data points (BJD = 2 453 951 and BJD = 2 456 468) were discarded because they had an unusually low signal-to-noise ratio or exhibited an anomalous CCF (high contrast), probably linked to an incorrect background correction. The data are presented in Table 10.

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 2050-day 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. Two-Keplerian model. A Neptune-mass object on a 35-day period orbit

The residuals of the one-Keplerian model show a sharp significant peak at 34.9-days and an additional trend with unconstrained period (Fig. 19). The signal at 34.9 days corresponds to the Neptune-mass 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 one-Keplerian model does not exhibit any significant peak (Fig. 19)<fnref n="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 log ($R'_{HK}$) time series does not exhibit any remaining significant peak after the long-term 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⁻¹. 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.9-day signal detected in the RV lead us to conclude that it is most likely planetary in origin.

The long-term trend is still seen in the residuals of the two-Keplerian 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 log ($R'_{HK}$), 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 third-degree polynomial. We note that the RV variation in the past ten years is about 10 m s⁻¹, which is similar to the dispersion around the one-Keplerian 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

<footnote n="7">⁷ The instrument is referred to as CORALIE98 before the upgrade and CORALIE07 after the upgrade.</footnote>

<footnote n="8">⁸ When a two-Keplerian model is fitted to the data from instruments other than HARPS, an RV amplitude significantly different from zero is found for the signal at 34.9 days. This does not mean, however, that the detection is significant. See discussion below.</footnote>
<table><caption>**Table 10.** HARPS and CORALIE measurements of HD 204313.</caption><thead><tr><th>BJD<br>-2 450 000</th><th>RV<br>(km s<sup>-1</sup>)</th><th>σ<sub>RV</sub><br>(km s<sup>-1</sup>)</th><th>BIS<br>(m s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>log (R'<sub>HK</sub>)</th><th>σ<sub>log (R'<sub>HK</sub>)</sub></th><th>Instrument</th></tr></thead><tbody><tr><td>53861.9016</td><td>-9.7468</td><td>0.0004</td><td>-26.22</td><td>7.0415</td><td>-5.0043</td><td>0.0026</td><td>HARPS</td></tr><tr><td>53862.8961</td><td>-9.7465</td><td>0.0004</td><td>-28.17</td><td>7.0409</td><td>-5.0026</td><td>0.0028</td><td>HARPS</td></tr><tr><td>53863.9005</td><td>-9.7430</td><td>0.0004</td><td>-29.17</td><td>7.0377</td><td>-5.0020</td><td>0.0028</td><td>HARPS</td></tr><tr><td>53864.8695</td><td>-9.7434</td><td>0.0004</td><td>-27.87</td><td>7.0383</td><td>-4.9913</td><td>0.0029</td><td>HARPS</td></tr><tr><td>53865.9012</td><td>-9.7446</td><td>0.0004</td><td>-27.85</td><td>7.0390</td><td>-5.0062</td><td>0.0029</td><td>HARPS</td></tr></tbody></table>

**Notes.** The full version is available at the CDS.

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 two-Keplerian model with an additional third-degree polynomial, which was employed to account for the long-term 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 least-squares fit to the log ($R'_{HK}$) time series (Table A.3).

As the log ($R'_{HK}$) is not available for the other instruments, the long-term drift cannot be monitored. Instead, we decided to perform the analysis in two steps equivalent to using a model where the long-term drift is only present for HARPS data: first, we modelled data from all instruments except for HARPS using a two-Keplerian model, without long-term drift, and with a constant-jitter 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 long-term third-degree polynomial and the varying-jitter model. In other words, we updated (in the Bayesian sense) the probability distributions for the parameters present in the two-Keplerian model using HARPS data, and sampled from the posterior of the new parameters associated with the long-term drift. The posterior distribution of the model without long-term trend was approximated by an uncorrelated multi-normal distribution for all parameters, except for the eccentricity of the 39.4-day 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 two-step procedure described above is based on two assumptions: 1) the long-term drift has a negligible effect on data from all instruments other than HARPS; and 2) the posterior distributions of the two-Keplerian 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 goodness-of-fit found and on the fact that no strong correlations are seen in the posterior distributions of the two-Keplerian model parameters.

We note that even without the HARPS data, both the period and RV semi-amplitude of the 34.9-day 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 iswise 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

<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Top: radial velocities of HD 204313 with the corresponding one-Keplerian 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 Lomb-Scargle periodogram of the residuals to the model plotted in the upper panel, exhibiting excess power at $P = 34.9$ days and a low-frequency trend. The horizontal dotted and dashed lines represent the 10% and 1% p-value levels, respectively. Bottom: GLS periodogram of the RV residuals to the one-Keplerian model without the HARPS data. The peak at 34.9 days is no longer significant, and the long-term trend has been replaced by a definite period at $P \sim 3500$ days with insignificant power.</figcaption></figure>

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 one-Keplerian
<figure><img src="image_22.png" /><figcaption>**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 two-Keplerian model and the red curve in the *lower panel* is the corresponding periodogram.</figcaption></figure>

model is orders of magnitude higher than the corresponding two-Keplerian model probability, confirming that the 34.9-day signal is not detected without the HARPS data.

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 phase-folded orbits and the model confidence intervals are shown in Fig. 23. The companion on a 34.9-day orbit is a Neptune-like 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 long-term stability of the system. We nevertheless performed numerical integrations using Mercury (Chambers 1999) for over $7 \times 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 semi-major 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 long-period companion announced by Robertson et al. (2012). Indeed, the long-term trend

seen in the residuals of the one-Keplerian 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 long-term 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 long-term drift exhibit a peak at 4.7 days with an amplitude of $70 \ \text{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 \pm 37$ using the P14 estimator, and $2180 \pm 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

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 long-period
<table><caption>**Table 11.** Parameter posteriors for the HD 204313 system.</caption><thead><tr><th>Orbital parameters</th><th></th><th>Planet c</th><th>Planet b</th></tr></thead><tbody><tr><td>Orbital period, P<sup>•</sup></td><td>[days]</td><td>34.905 ± 0.012</td><td>2024.1 ± 3.1</td></tr><tr><td>RV amplitude, K<sup>•</sup></td><td>[m s<sup>-1</sup>]</td><td>3.42 ± 0.22</td><td>68.45 ± 0.30</td></tr><tr><td>Eccentricity, e<sup>•</sup></td><td></td><td>0.155 ± 0.071; &lt; 0.310<sup>†</sup></td><td>0.0946 ± 0.0032</td></tr><tr><td>Argument of periastron, ω<sup>•</sup></td><td>[deg]</td><td>238 ± 22</td><td>292.5 ± 2.7</td></tr><tr><td>Mean longitude at epoch, L<sub>0</sub><sup>•</sup></td><td>[deg]</td><td>46.5 ± 3.5</td><td>107.55 ± 0.24</td></tr><tr><td>Semi-major axis of relative orbit, a</td><td>[AU]</td><td>0.2099 ± 0.0071</td><td>3.167 ± 0.12</td></tr><tr><td>Minimum mass, M sin i</td><td>[M<sub>⊕</sub>]</td><td>17.6 ± 1.7</td><td>1360 ± 94</td></tr><tr><td>Minimum mass, M sin i</td><td>[M<sub>J</sub>]</td><td>0.0553 ± 0.0053</td><td>4.28 ± 0.30</td></tr><tr><td>Velocity drift<sup>*</sup></td><td></td><td></td><td></td></tr><tr><td>Systemic velocity, V<sub>0</sub><sup>•</sup></td><td>[km s<sup>-1</sup>]</td><td>-9.73922 ± 3.3 × 10<sup>-4</sup></td><td></td></tr><tr><td>Scaling constant, α<sup>•</sup></td><td>[m s<sup>-1</sup>/dex]</td><td>140 ± 35</td><td></td></tr><tr><td>Linear<sup>•</sup></td><td>[dex yr<sup>-1</sup>]</td><td>-0.0100 ± 0.0013</td><td></td></tr><tr><td>Quadratic<sup>•</sup></td><td>[dex yr<sup>-2</sup>]</td><td>(-8.1 ± 2.4) × 10<sup>-4</sup></td><td></td></tr><tr><td>Cubic<sup>•</sup></td><td>[dex yr<sup>-3</sup>]</td><td>(2.97 ± 0.68) × 10<sup>-4</sup></td><td></td></tr><tr><td>Noise model<sup>‡</sup></td><td></td><td></td><td></td></tr><tr><td>Additional noise at log (R'<sub>HK</sub>) = -5, σ<sub>J</sub>|<sub>-5.0</sub><sup>•</sup></td><td>[m s<sup>-1</sup>]</td><td>1.28 ± 0.21</td><td></td></tr><tr><td>Slope, α<sub>J</sub><sup>•</sup></td><td>[m s<sup>-1</sup>/dex]</td><td>9.6 ± 6.3</td><td></td></tr><tr><td>rms(O - C)</td><td>[m s<sup>-1</sup>]</td><td>1.32 ± 0.11</td><td></td></tr><tr><td>Other parameters<sup>◇</sup></td><td></td><td></td><td></td></tr><tr><td>CORALIE98 systemic velocity</td><td>[km s<sup>-1</sup>]</td><td>9.7888 ± 0.0016</td><td></td></tr><tr><td>CORALIE07 systemic velocity</td><td>[km s<sup>-1</sup>]</td><td>9.7754 ± 0.0015</td><td></td></tr><tr><td>McDonald systemic velocity</td><td>[km s<sup>-1</sup>]</td><td>9.6878 ± 0.0011</td><td></td></tr></tbody></table>

**Notes.** The epoch is BJD = 2454993.8485. (*) MCMC jump parameter. (†) Eccentricity does not differ significantly from zero at the 99%-level. Upper limit is reported. (**) See Sect. 4.3. (‡) The additional (stellar) noise for measurement *i* is $\sigma_{Ji} = \sigma_J|_{-5.0} + \alpha_J \cdot (\log(R'_{HK})_i + 5.0)$. (◦) Parameters from the preliminary analysis without HARPS data (see text for details), listed here for convenience.

<figure><img src="image_23.png" /><figcaption>**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 than the corresponding distribution of planet c, the vertical scale was set to 20 times larger than the scale shown for planet c.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 23.** Radial velocity data phase-folded to the best-fit period of each of the two planetary candidates in the HD 204313 system after subtracting the effect of the other signal and the long-term drift. The blue lines represent the 95%- highest density interval (HDI).</figcaption></figure>

companion. The data sets span more than ten years and reveal long-term 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 short-term activity effect (the so-called jitter of the star), which is much harder to describe using a deterministic model. The long-term activity signal has been shown to correlate with the evolution of the activity proxy log (R'_{HK}) (Lovis et al. 2011a), and we therefore used the HARPS

measurement of this activity proxy in time to provide priors for the long-term activity effect. Some assumptions are key to our model: 1) the RV signal produced by the long-term activity variability scales linearly with log (R'_{HK}); and 2) the stellar jitter is
adequately described as an additional Gaussian noise with an
amplitude that scales linearly with log (R'_{HK}). We have tested an
alternative model with a fixed amplitude for the additional noise,
but the data mildly preferred the varying-jitter model.

Potential periodicities were identified using a generalised
Lomb-Scargle periodogram and by resorting to a genetic algo-
rithm. 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 com-
parison. 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 peri-
odogram analysis that is usually employed by the HARPS team
to detect weak signals. One of the reasons involves the uncertain-
ties 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 *p*-value obtained for the periodogram
peak power as a false-alarm 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
well-established method for computing the probability of all in-
volved 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 signal-to-noise
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% *p*-value 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 super-Earth planet candidates, with orbital periods $P = 5.77$ days and $P = 13.5$ days. The 13.5-day 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 $\log(R'_{HK})$. The minimum masses are $M \sin i = 6.4 M_\oplus$ and $M \sin i = 5.6 M_\oplus$ for the 5.77-day and 13.5-day signals, respectively. The long-term activity signal has an amplitude of $1.5 \ \text{m s}^{-1}$ and a period of $9.6$ yr.

- HD 40307 hosts four certain planetary companions at peri-
  ods of between 4.3 and 51.6 days and with masses of be-
  tween $M \sin i = 3.6 M_\oplus$ and $M \sin i = 8.7 M_\oplus$. We find
  inconclusive evidence for an additional companion at $P \sim$
  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 Neptune-mass planet on
  a 34.9-day orbit and a 4.3-$M_J$ candidate on an outer orbit
  ($P = 2024$ days).

Only six systems with a massive ($M > 4 M_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.8-$M_J$ candidate on a 14.3-day orbit,
and a 12.3-$M_J$ object on an 2140-day orbit. The mass ratio of
this system is around five times lower than for HD 204313, how-
ever. The system around HD 74156 (Naef et al. 2004) consists of
two companions in orbits alike to those of the HD 204313 sys-
tem (a 1.8 $M_J$ object at 51.6 days and an 8.2 $M_J$ companion on
a ~2500-day 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 fu-
ture permit unveiling further planetary signals and improving the
modelling of the activity effect.

HD 1461 and HD 40307 are two of only six solar-type plan-
etary systems with at least two super-Earth companions whose
mass is known to better than 50%, and they will become prime
targets for the follow-up 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 Kepler-102 (Marcy
et al. 2014). Except for Kepler-102, which was first discov-
ered 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, ac-
tivity cycles are omnipresent. Correctly modelling their effect
on RV data is therefore necessary to unveil the full system archi-
tecture, from the short-period companions out to the habitable-
zone planets. We have shown evidence that periodicities arise at
the aliasing frequencies when an incomplete correction of the
long-term 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 high-activity periods can save hours
of telescope time.

**Acknowledgements.** 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/2007-2013) 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

Anglada-Escudé, G., Tuomi, M., Gerlach, E., et al. 2013, A&A, 556, A126

Anglada-Escudé, G., Arriagada, P., Tuomi, M., et al. 2014, MNRAS, 443, L89

Anglada-Escudé, G., Tuomi, M., Arriagada, P., et al. 2015, ApJ, submitted [arXiv:1506.09072]

Arriagada, P. 2011, ApJ, 734, 70

Baliunas, S. L., Donahue, R. A., Soon, W. H., et al. 1995, ApJ, 438, 269

Baluev, R. V. 2013, MNRAS, 436, 807

Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373

Boisse, I., Bouchy, F., Hébrard, G., et al. 2011, A&A, 528, A4

Boisse, I., Bonfils, X., & Santos, N. C. 2012, A&A, 545, A109

Chambers, J. E. 1999, MNRAS, 304, 793

Chib, S., & Jeliazkov, I. 2001, J. Amer. Statis. Assoc., 96, 270

Cincunegui, C., Díaz, R. F., & Mauas, P. J. D. 2007, A&A, 469, 309

Dawson, R. I., & Fabrycky, D. C. 2010, ApJ, 722, 937

Delfosse, X., Bonfils, X., Forveille, T., et al. 2013, A&A, 553, A8

Desort, M., Lagrange, A.-M., Galland, F., Udry, S., & Mayor, M. 2007, A&A, 473, 983

Díaz, R. F., Almenara, J. M., Santerne, A., et al. 2014, MNRAS, 441, 983

Dravins, D. 1982, ARA&A, 20, 61

Dumusque, X., Lovis, C., Ségransan, D., et al. 2011a, A&A, 535, A55

Dumusque, X., Udry, S., Lovis, C., Santos, N. C., & Monteiro, M. J. P. F. G. 2011b, A&A, 525, A140

Dumusque, X., Pepe, F., Lovis, C., et al. 2012, Nature, 491, 207

Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132

Dumusque, X., Pepe, F., Lovis, C., & Latham, D. W. 2015, ApJ, 808, 171

van Leeuwen, F. 2007, Hipparcos, the New Reduction of the Raw Data, Astrophys. Space Sci. Lib., 350

Feroz, F., & Hobson, M. P. 2014, MNRAS, 437, 3540

Figueira, P., Marmier, M., Boué, G., et al. 2012, A&A, 541, A139

Fischer, D. A., Marcy, G. W., Butler, R. P., et al. 2001, ApJ, 551, 1107

Fischer, D. A., Marcy, G. W., Butler, R. P., et al. 2003, ApJ, 586, 1394

Friel, N., & Wyse, J. 2012, Statistica Neerlandica, 66, 288

Fulton, B. J., Weiss, L. M., Sinukoff, E., et al. 2015, ApJ, 805, 175

Geweke, J. 1989, Econometrica, 57, 1317

Gomes da Silva, J., Santos, N. C., Bonfils, X., et al. 2012, A&A, 541, A9

Gray, D. F., & Baliunas, S. L. 1995, ApJ, 441, 436

Gray, D. F., Baliunas, S. L., Lockwood, G. W., & Skiff, B. A. 1996, ApJ, 465, 945

Gregory, P. C. 2005, Bayesian Logical Data Analysis for the Physical Sciences: A Comparative Approach with Mathematica Support (Cambridge, UK: Cambridge University Press)

Gregory, P. C. 2007, MNRAS, 381, 1607

Gregory, P. C. 2011, MNRAS, 415, 2523

Gregory, P. C. 2012, ArXiv e-prints [arXiv:1212.4058]

Hall, J. C., Lockwood, G. W., & Skiff, B. A. 2007, AJ, 133, 862

Hall, J. C., Henry, G. W., Lockwood, G. W., Skiff, B. A., & Saar, S. H. 2009, AJ, 138, 312

Hatzes, A. P. 2002, Astron. Nachr., 323, 392

Hatzes, A. P. 2013, ApJ, 770, 133

Henry, T. J., Soderblom, D. R., Donahue, R. A., & Baliunas, S. L. 1996, AJ, 111, 439

Isaacson, H., & Fischer, D. 2010, ApJ, 725, 875

Kass, R. E., & Raftery, A. E. 1995, J. Am. Statist. Assn., 90, 773

Kipping, D. M. 2013, MNRAS, 434, L51

Lindegren, L., & Dravins, D. 2003, A&A, 401, 1185

Lissauer, J. J., Marcy, G. W., Rowe, J. F., et al. 2012, ApJ, 750, 112

Lissauer, J. J., Marcy, G. W., Bryson, S. T., et al. 2014, ApJ, 784, 44

Lovis, C., Dumusque, X., Santos, N. C., et al. 2011a, A&A, unpublished [arXiv:1107.5325]

Lovis, C., Ségransan, D., Mayor, M., et al. 2011b, A&A, 528, A112

Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264

Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20

Mauas, P. J. D., Buccino, A., Díaz, R., et al. 2012, in IAU Symp. 286, eds. C. H. Mandrini, & D. F. Webb, 317

Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20

Mayor, M., Udry, S., Lovis, C., et al. 2009, A&A, 493, 639

Mayor, M., Marmier, M., Lovis, C., et al. 2011, A&A, submitted [arXiv:1109.2497]

Meunier, N., & Lagrange, A.-M. 2013, A&A, 551, A101

Meunier, N., Desort, M., & Lagrange, A.-M. 2010, A&A, 512, A39

Motalebi, F., Udry, S., Gillon, M., et al. 2015, A&A, 584, A72

Naef, D., Mayor, M., Beuzit, J. L., et al. 2004, A&A, 414, 351

Newton, M. A., & Raftery, A. E. 1994, J. Roy. Statist. Soc., Ser. B (Methodological), 3

Noyes, R. W., Hartmann, L. W., Baliunas, S. L., Duncan, D. K., & Vaughan, A. H. 1984, ApJ, 279, 763

Pepe, F., Mayor, M., Galland, F., et al. 2002, A&A, 388, 632

Pepe, F., Lovis, C., Ségransan, D., et al. 2011, A&A, 534, A58

Perrakis, K., Ntzoufras, I., & Tsionas, E. G. 2014, Comp. Stat. Data Analysis, 77, 54

Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49

Rivera, E. J., Butler, R. P., Vogt, S. S., et al. 2010, ApJ, 708, 1492

Robert, C. P., & Wraith, D. 2009, in AIP Conf. Ser. 1193, eds. P. M. Goggans, & C.-Y. Chan, 251

Robertson, P., Horner, J., Wittenmyer, R. A., et al. 2012, ApJ, 754, 50

Robertson, P., Endl, M., Cochran, W. D., & Dodson-Robinson, S. E. 2013, ApJ, 764, 3

Robertson, P., Mahadevan, S., Endl, M., & Roy, A. 2014, Science, 345, 440

Robertson, P., Roy, A., & Mahadevan, S. 2015, ApJ, 805, L22

Saar, S. H., & Donahue, R. A. 1997, ApJ, 485, 319

Santos, N. C., Gomes da Silva, J., Lovis, C., & Melo, C. 2010, A&A, 511, A54

Schwarz, G. et al. 1978, The annals of statistics, 6, 461

Ségransan, D., Udry, S., Mayor, M., et al. 2010, A&A, 511, A45

Ségransan, D., Mayor, M., Udry, S., et al. 2011, A&A, 535, A54

Sellke, T., Bayarri, M., & Berger, J. O. 2001, Am. Statist., 55, 62

Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373

Tuomi, M. 2012, A&A, 543, A52

Tuomi, M., & Jones, H. R. A. 2012, A&A, 544, A116

Tuomi, M., Anglada-Escudé, G., Gerlach, E., et al. 2013a, A&A, 549, A48

Tuomi, M., Jones, H. R. A., Jenkins, J. S., et al. 2013b, A&A, 551, A79

Udry, S., Bonfils, X., Delfosse, X., et al. 2007, A&A, 469, L43

Wright, J. T., Fischer, D. A., Ford, E. B., et al. 2009, ApJ, 699, L97

Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577
## Appendix A: Parameter priors

<table><caption>**Table A.1.** Parameter prior distributions for the HD 1461 system.</caption><thead><tr><th rowspan="2">Orbital parameters</th><th colspan="2">Prior distribution</th></tr><tr><th>Planet b / c</th><th>Magnetic cycle<sup>*</sup></th></tr></thead><tbody><tr><td>Orbital period, <i>P</i></td><td>[days]</td><td><i>J</i>(1.0, 10<sup>4</sup>)</td><td><i>N</i>(3522.4, 75.4)</td></tr><tr><td>RV amplitude, <i>K</i></td><td>[m s<sup>-1</sup>]</td><td><i>U</i>(0.0, 200)</td><td><i>U</i>(0.0, 200)</td></tr><tr><td>Eccentricity, <i>e</i></td><td></td><td><i>B</i>(0.867, 3.03)</td><td>—</td></tr><tr><td>Argument of periastron, <i>ω</i></td><td>[deg]</td><td><i>U</i>(0.0, 360.0)</td><td>—</td></tr><tr><td><i>e</i><sup>1/2</sup> sin(<i>ω</i>)</td><td></td><td>—</td><td><i>N</i>(-0.388, 0.078)</td></tr><tr><td><i>e</i><sup>1/2</sup> cos(<i>ω</i>)</td><td></td><td>—</td><td><i>N</i>(0.129, 0.096)</td></tr><tr><td>Mean longitude at epoch, <i>L</i><sub>0</sub></td><td>[deg]</td><td><i>U</i>(0.0, 360.0)</td><td><i>N</i>(124.0, 21.0)</td></tr><tr><td>Systemic velocity, <i>V</i><sub>0</sub></td><td>[km s<sup>-1</sup>]</td><td><i>U</i>(-10.061, -10.055)</td><td></td></tr><tr><td colspan="4"><strong>Noise model<sup>‡</sup></strong></td></tr><tr><td>Additional noise at log (<i>R</i>'<sub>HK</sub>) = -5, σ<sub>J</sub>|<sub>-5.0</sub></td><td>[m s<sup>-1</sup>]</td><td><i>U</i>(0, 50)</td><td></td></tr><tr><td>Slope, α<sub>J</sub></td><td>[m s<sup>-1</sup>/dex]</td><td><i>U</i>(-200, 200)</td><td></td></tr></tbody></table>

**Notes.** The epoch is BJD = 2 455 155.3854 for planets *b* and *c* and BJD = 2 455 195.8367 for the magnetic cycle. *U*(*x*<sub>min</sub>; *x*<sub>max</sub>): uniform distribution between *x*<sub>min</sub> and *x*<sub>max</sub>. *J*(*x*<sub>min</sub>; *x*<sub>max</sub>): Jeffreys (log-flat) distribution between *x*<sub>min</sub> and *x*<sub>max</sub>. *N*(*μ*; *σ*): normal distribution with mean *μ* and standard deviation *σ*. *B*(*a*, *b*): beta distribution. (*) See Sect. 4.3. (‡) The additional (stellar) noise for measurement *i* is σ<sub>Ji</sub> = σ<sub>J</sub>|<sub>-5.0</sub> + α<sub>J</sub> · (log (<i>R</i>'<sub>HK</sub>)<sub>i</sub> + 5.0).

<table><caption>**Table A.2.** Prior distributions for the HD 40307 system.</caption><thead><tr><th>Orbital parameters</th><th>Prior distribution</th></tr></thead><tbody><tr><td>Orbital period, <i>P</i></td><td>[days]</td><td><i>J</i>(1.0, 10<sup>4</sup>)</td></tr><tr><td>RV amplitude, <i>K</i></td><td>[m s<sup>-1</sup>]</td><td><i>U</i>(0.0, 200)</td></tr><tr><td>Eccentricity, <i>e</i></td><td></td><td><i>B</i>(0.867, 3.03)</td></tr><tr><td>Argument of periastron, <i>ω</i></td><td>[deg]</td><td><i>U</i>(0, 360)</td></tr><tr><td>Mean longitude at epoch, <i>L</i><sub>0</sub></td><td>[deg]</td><td><i>U</i>(0, 360)</td></tr><tr><td>Systemic velocity, <i>V</i><sub>0</sub></td><td>[km s<sup>-1</sup>]</td><td><i>U</i>(28.996, 33.668)</td></tr><tr><td colspan="3"><strong>Velocity drift (long-term activity effect)<sup>*</sup></strong></td></tr><tr><td>Scaling constant, α</td><td>[m s<sup>-1</sup>/dex]</td><td><i>U</i>(0, 100)</td></tr><tr><td>Linear</td><td>[10<sup>-4</sup> dex yr<sup>-1</sup>]</td><td><i>N</i>(373.8, 8.6)</td></tr><tr><td>Quadratic</td><td>[10<sup>-4</sup> dex yr<sup>-2</sup>]</td><td><i>N</i>(46.3, 4.6)</td></tr><tr><td>Cubic</td><td>[10<sup>-4</sup> dex yr<sup>-3</sup>]</td><td><i>N</i>(-17.3, 1.0)</td></tr><tr><td colspan="3"><strong>Noise model<sup>‡</sup></strong></td></tr><tr><td>Additional noise at log (<i>R</i>'<sub>HK</sub>) = -5, σ<sub>J</sub>|<sub>-5.0</sub></td><td>[m s<sup>-1</sup>]</td><td><i>U</i>(0, 50)</td></tr><tr><td>Slope, α<sub>J</sub></td><td>[m s<sup>-1</sup>/dex]</td><td><i>U</i>(0, 50)</td></tr></tbody></table>

**Notes.** Priors are identical for all signals in the model. The epoch is BJD = 2 454 521.6791. The argument of periastron *ω* is unconstrained for these nearly-circular orbits. *U*(*x*<sub>min</sub>; *x*<sub>max</sub>): uniform distribution between *x*<sub>min</sub> and *x*<sub>max</sub>. *J*(*x*<sub>min</sub>; *x*<sub>max</sub>): Jeffreys (log-flat) distribution between *x*<sub>min</sub> and *x*<sub>max</sub>. *N*(*μ*; *σ*): normal distribution with mean *μ* and standard deviation *σ*. *B*(*a*, *b*): beta distribution. (*) See Sect. 4.3. (‡) The additional (stellar) noise for measurement *i* is σ<sub>Ji</sub> = σ<sub>J</sub>|<sub>-5.0</sub> + α<sub>J</sub> · (log (<i>R</i>'<sub>HK</sub>)<sub>i</sub> + 5.0).

<table><caption>**Table A.3.** Parameter prior distributions for the HD 204313 system.</caption><thead><tr><th rowspan="2">Orbital parameters</th><th colspan="2">Prior distribution</th></tr><tr><th>Planet b</th><th>Planet c</th></tr></thead><tbody><tr><td>Orbital period, <i>P</i></td><td>[days]</td><td><i>N</i>(2046.3, 9.1)</td><td><i>N</i>(34.989, 0.033)</td></tr><tr><td>RV amplitude, <i>K</i></td><td>[m s<sup>-1</sup>]</td><td><i>N</i>(66.9, 1.3)</td><td><i>N</i>(4.9, 1.0)</td></tr><tr><td>Eccentricity, <i>e</i></td><td></td><td><i>N</i>(0.125, 0.017)</td><td><i>B</i>(0.93, 5.5)</td></tr><tr><td>Argument of periastron, <i>ω</i></td><td>[deg]</td><td><i>N</i>(303.0, 8.3)</td><td><i>U</i>(0, 360)</td></tr><tr><td>Mean longitude at epoch, <i>L</i><sub>0</sub></td><td>[deg]</td><td><i>N</i>(109.10, 0.91)</td><td><i>N</i>(35.6, 15.0)</td></tr><tr><td>Systemic velocity, <i>V</i><sub>0</sub></td><td>[km s<sup>-1</sup>]</td><td><i>U</i>(-9, 762, -9.697)</td><td></td></tr><tr><td colspan="4"><strong>Velocity drift (long-term activity effect)<sup>*</sup></strong></td></tr><tr><td>Scaling constant, α</td><td>[m s<sup>-1</sup>/dex]</td><td><i>U</i>(0, 300)</td><td></td></tr><tr><td>Linear</td><td>[10<sup>-2</sup> dex yr<sup>-1</sup>]</td><td><i>N</i>(-1.02, 0.13)</td><td></td></tr><tr><td>Quadratic</td><td>[10<sup>-4</sup> dex yr<sup>-2</sup>]</td><td><i>N</i>(-5.5, 2.5)</td><td></td></tr><tr><td>Cubic</td><td>[10<sup>-4</sup> dex yr<sup>-3</sup>]</td><td><i>N</i>(3.68, 0.92)</td><td></td></tr><tr><td colspan="4"><strong>Noise model<sup>‡</sup></strong></td></tr><tr><td>Additional noise at log (<i>R</i>'<sub>HK</sub>) = -5, σ<sub>J</sub>|<sub>-5.0</sub></td><td>[m s<sup>-1</sup>]</td><td><i>U</i>(0, 50)</td><td></td></tr><tr><td>Slope, α<sub>J</sub></td><td>[m s<sup>-1</sup>/dex]</td><td><i>U</i>(-200, 200)</td><td></td></tr></tbody></table>

**Notes.** The epoch is BJD = 2 454 993.84858. *U*(*x*<sub>min</sub>; *x*<sub>max</sub>): uniform distribution between *x*<sub>min</sub> and *x*<sub>max</sub>. *N*(*μ*; *σ*): normal distribution with mean *μ* and standard deviation *σ*. *B*(*a*, *b*): beta distribution. (*) See Sect. 4.3. (‡) The additional (stellar) noise for measurement *i* is σ<sub>Ji</sub> = σ<sub>J</sub>|<sub>-5.0</sub> + α<sub>J</sub> · (log (<i>R</i>'<sub>HK</sub>)<sub>i</sub> + 5.0).