Main final dvi

• ǫ

i,lNS

are the residuals to a fit. We assume that each

ǫ

i,lNS

is a Gaussian random variable

with a zero mean and a variance of

σ

2

i

+ σ

2

INS

, where

σ

2

i

is the reported uncertainty of the

i-th measurement and σ

2

INS

is the jitter parameter and represents the excess white noise not

included in

σ

2

i

.

• γ

INS

is the zero-point velocity of each instrument. Each INS can have a different zero-point

depending on how the radial velocities are measured and how the wavelengths are calibrated.

• ˙γ is a linear trend parameter caused by a long term acceleration.

• The term κ(∆t

i

) is the superposition of k Keplerian signals evaluated at ∆t

i

. Each Keplerian

signal depends on five parameters: the orbital period

P

p

, semi-amplitude of the signal

K

p

,

mean anomaly

M

0

,p

, which represents the phase of the orbit with respect to the periastron of

the orbit at

t

0

, orbital eccentricity

e

p

that goes from

0 (circular orbit) to 1 (unbound parabolic

orbit), and the argument of periastron

ω

p

, which is the angle on the orbital plane with respect

to the plane of the sky at which the star goes through the periastron of its orbit (the planet’s

periastron is at

ω

p

+180 deg). Detailed definitions of the parameters can be found elsewhere.

37

• The Moving Average term

MA

i,INS

= φ

INS

exp

t

i−1

− t

i

τ

INS

ǫ

i−1,INS

(5)

is a simple parameterization of possible correlated noise that depends on the residual of the

previous measurement

ǫ

i−1,INS

. As for the other parameters related to noise in our model,

we assume that the parameters of the MA function depend on the instrument; for example the

different wavelength ranges used will cause different properties of the instrumental systematic

noise. Keplerian and other physical processes also introduce correlations into the data, there-

fore some degree of degeneracy between the MA terms and the signals of interest is expected.

As a result, including a MA term always produces more conservative significance estimates

than a model with uncorrelated random noise only. The MA model is implemented through a

coefficient

φ

INS

and a time-scale

τ

INS

.

φ

INS

quantifies the strength of the correlation between

the

i and i − 1 measurements. It is bound between −1 and 1 to guarantee that the process

is stationary (i.e. the contribution of the MA term does not arbitrarily grow over time). The

exponential smoothing is used to decrease the strength of the correlation exponentially as the

difference

t

i

− t

i−1

increases.

38

• Linear correlations with activity indices can also be included in the model in the following

manner,

A

i,INS

=

ξ

C

ξ,INS

ξ

i,INS

(6)

where

ξ runs over all the activity indices used to model each INS dataset (e.g. m

2

,

m

3

,

S-index, etc. whose description is provided below). To avoid any confusion with other dis-

cussions about correlations, we call these

C

ξ,INS

activity coefficients. Note that each activity

12

---

coefficient

C

ξ,INS

is associated to one activity index (

ξ

i

) obtained simultaneously with the i-th

radial velocity measurement (e.g. chromospheric emission from the H

α

line, second moment

of the mean-line profile, interpolated photometric flux, etc.). When fitting a model to the data,

an activity coefficient significantly different from

0 indicates evidence of Doppler variability

correlated with the corresponding activity index. Formally speaking, these

C

ξ,INS

correspond

to the coefficient of the first order Taylor expansion of a physical model for the apparent radial

velocities as a function of the activity indices and other physical properties of the star.

A simplified version of the same likelihood model is used when analyzing time-series of activity

indices. That is, when searching for periodicities in series other than Doppler measurements, the

model will consist of the

γ

INS

zero-points, a linear trend term

˙γ∆t

i

, and a sum of

n sinusoids

ˆ

κ(t

i

, θ)

=

n

k

A

k

sin

2π∆t

i

P

k

+ B

k

cos

2π∆t

i

P

k

(7)

where each

k-th sinusoid has three parameters A

k

,

B

k

, and

P

k

instead of the five Keplerian ones.

Except for the period parameters and the jitter terms, this model is linear with all the other parame-

ters, which allows a relatively quick computation of the likelihood-ratio periodograms.

1.3

Bayesian prior choices.

As in any Bayesian analysis, the prior densities of the model parameters have to be selected in a

suitable manner (for example see

39

). We used uniform and uninformative distributions for most

of the parameters apart from a few, possibly significant, exceptions. First, as we used a parameter

l = ln P in the MCMC samplings instead of the period P directly, the uniform prior density π(l) = c

for all

l ∈ [ln T

0

, ln T

max

], where T

0

and

T

max

are some minimum and maximum periods, does not

correspond to a uniform prior in

P . Instead, this prior corresponds to a period prior such that

π(P ) ∝ P

−1

.

40

We made this choice because the period can be considered a “scale parameter” for

which an uninformative prior is one that is uniform in

ln P .

41

We selected the parameter space of

the period such that

T

0

= 1 day and T

max

= T

obs

, where

T

obs

is the baseline of the combined data.

For the semi amplitude parameter

K, we used a π(K) = c for all K ∈ [0, K

max

], where K

max

was selected as

K

max

= 10 ms

−1

because the RMSs of the Doppler series did not exceed 3 ms

−1

in any of the sets. Following previous works,

40, 42

we chose the prior for the orbital eccentricities as

π(e) ∝ N (0, Σ

2

e

), where e is bound between zero (circular orbit) and 1. We set this Σ

2

e

= 0.3 to

penalize high eccentricities while keeping the option of high

e if the data strongly favours it.

We also used an informative prior for the excess white noise parameter of

σ

IN S

for each instru-

ment. Based on analyses of a sample of M dwarfs,

15

this “stellar jitter” is typically very close to

a value of 1 ms

−1

. Thus, we used a prior such that

π(σ

l

) ∝ N (µ

σ

, σ

2

σ

) such that the parameters

were selected as

µ

σ

= σ

σ

= 1 ms

−1

. Uniform priors were used in all the activity coefficients

C

ξ

∈ [−C

ξ,max

, C

ξ,max

]. For practical purposes, the time-series of all activity indices were mean

subtracted and normalized to their RMS. This choice allows us to select the bounds of the activity

coefficients for the renormalized time-series as ˆ

C

ξ,max

= 3 ms

−1

, so that adding correlation terms

13

---