© ESO, 2016

1. Introduction

More than 600 exoplanets were discovered and many others confirmed using radial-velocity (RV) observations (www.exoplanet.eu, Schneider et al. 2011). The method is more sensitive to massive planets orbiting close to their stars, since Neptune- and Earth-mass planets induce very small amplitude radial-velocity signals, often at the level of the current observational uncertainties of 1 ms^-1 and weaker. In addition, stellar activity can induce false-positive signals that mimic and mask the radial velocity signature of a low-mass planet (e.g. Bonfils et al. 2007; Robertson & Mahadevan 2014; Santos et al. 2014). Only with very high precision Doppler spectroscopy and dense sampling is one able to detect these planets (e.g. Dumusque et al. 2012; Hatzes 2014) and to push the detection limits in the direction of a possibly habitable Earth-like planet, the same stars must be followed for a very long time.

The increasing number of detected planets provides constraints for models of planet formation and evolution (e.g. Udry & Santos 2007; Ida & Lin 2004; Mordasini et al. 2009). In particular, the metallicity and, in general, the chemical abundances of the host star, are now known to be a key ingredient in planet formation (Ida & Lin 2004; Mordasini et al. 2015). According to the core-accretion theory, proto-planetary disks with higher metallicity are able to form rocky or icy cores in time for runaway accretion to lead to the formation of a giant planet before disk dissipation occurs. In lower metallicity disks, the cores do not grow fast enough to accrete gas in large quantities before disk dissipation occurs, implying a lower fraction of giant planets.

Gonzalez (1997) showed the first observational hints for a correlation between the presence of giant planets and metallicity. As more and more planets were discovered, this correlation was solidly confirmed: it is more likely to detect a giant planet orbiting a metal-rich star (Santos et al. 2004; Fischer & Valenti 2005; Sousa et al. 2011b). This result was also confirmed in data from transit surveys (e.g. Buchhave et al. 2012). Furthermore, it is also now known that the metallicity, or even specific chemical abundance ratios, can have a strong impact on the planet formation efficiency, composition, and architecture (Guillot et al. 2006; Adibekyan et al. 2013, 2015; Dawson et al. 2015; Santos et al. 2015).

Because giant planets are more frequent around metal-rich stars, several RV surveys became biased towards metal-rich samples (Tinney et al. 2002; Fischer & Valenti 2005; Melo et al. 2007; Jenkins et al. 2013a). Nevertheless, some programmes focused on the metal-poor regime and tried to determine the frequency of giant planets orbiting metal-poor stars and the metallicity limit below which no giant planets can be observed (Sozzetti et al. 2009; Santos et al. 2011; Mortier et al. 2012, 2013). Until recently, however, these programmes focused on the search for massive, giant planets.

Neptune-mass planets, in contrast, have been found to have a relatively flat metallicity distribution (Udry et al. 2006; Sousa et al. 2008, 2011b; Neves et al. 2013). In systems with only hot Neptunes, the metallicity distribution actually becomes slightly metal-poor (Mayor et al. 2011; Buchhave et al. 2012). These results are supported by theoretical models that predict lower mass planets are common around stars with a wide range of metallicities. One should note here the works of Wang & Fischer (2014) and Buchhave & Latham (2015), who find conflicting results about the existence of a universal planet-metallicity correlation extending to the terrestrial planets, and Adibekyan et al. (2012), who show that, in the metal-poor regime, planets are more prevalent on α-enhanced stars (i.e. of higher global metallicity).

The number of detected low-mass planets orbiting metal-poor stars is nevertheless still small. In an effort to detect these planets and explore the low-metallicity regime, our team started a Large Programme with HARPS to follow a sample of about 100 metal-poor stars. The programme and sample were presented in Santos et al. (2014), to which we point the reader for more details. The main goal of the survey is to derive observational constraints on the frequency of Neptunes and super-Earths in the metal-poor regime. These estimates will be compared with the ones from the HARPS-GTO programme (Mayor et al. 2011), which searched for very low-mass planets orbiting a sample of solar-neighbourhood stars (thus with solar metallicities). When combined, the two surveys will set important constraints for planet formation and evolution models and help in providing an estimate of the frequency of planets in our Galaxy.

This paper presents the analysis of the stars in the metal-poor sample that were observed on more than 75 nights. Some of these stars now have precise radial velocity observations covering a baseline of over ten years. In Sect. 2 we present the stellar parameters for these stars and their radial-velocity observations. The method used to search for Keplerian signals is outlined in Sect. 3 and in Sect. 4 the results of its application. Section 5 presents a simple comparison of the model selection criteria we consider, and detection limits for individual stars are derived in Sect. 6. A first constraint on the occurrence of planets in this sample and our conclusions are presented in Sect. 7.

2. Sample and observations

Our complete sample of metal-poor stars contains 109 targets, chosen from a survey for giant planets orbiting metal-poor stars (Santos et al. 2007), and a programme to search for giant planets orbiting a volume-limited sample of FGK dwarfs (Naef et al. 2007). The criteria for defining this sample are detailed in Santos et al. (2014). Stellar parameters were derived for the complete sample from a set of high-resolution HARPS spectra (Sousa et al. 2011a; Santos et al. 2014).

In this work we analyse those stars that, up to December 2014, were observed on more than 75 nights. Though not strictly motivated, this number of nights means that in the worst case scenario, we have six data points per parameter for a model with two planets. A total of 15 stars meet this criterion. In Table 1 we present the stellar parameters for this subsample. The table lists the stellar mass, effective temperature, surface gravity, metallicity, B − V colour, activity level, and estimated rotation period.

The activity level of the stars, denoted here by the weighted average of the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ values (Noyes et al. 1984b), was derived from the analysis of the CaII H and K lines in the HARPS spectra (e.g. Santos et al. 2000; Lovis et al. 2011). Estimates for the rotation period of each star were derived from the activity-rotation calibrations of Noyes et al. (1984b) and Mamajek & Hillenbrand (2008). The typical error bar on these estimates, computed from the dispersion in the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ time series, is close to the difference between the estimates from the two calibrations. For HD 31128, the rotation period derived in this way is around 2.8 d. The very low metallicity of this star, however, places it outside the calibration range of the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ flux calibration and of the adopted activity-rotation relations. This renders the rotation period estimate uncertain, so we chose to omit it from the table.

Table 2 shows the number of nights each star was observed, the mean error bar $\overline{{\mathit{\sigma }}_{\mathit{i}}}$, and the weighted standard deviation of the radial velocities s_RV. The ratio of the last two quantities gives an indication of whether the radial-velocities vary in excess of the internal errors and by how much. This table also shows the baseline of observations.

The distributions of metallicity, effective temperature, and surface gravity of the full sample and of the 15 stars studied here are shown in Fig. 1. By selecting stars based on the number of observations, one can introduce biases that depend on the criteria used for scheduling observations. From the distributions in Fig. 1, we do not see a clear bias towards warmer or more evolved stars, which may show more significant (stellar) variability. The subsample can nevertheless be biased in favour of stars hosting low-mass planets, since these stars will also show more significant variability. This can have an impact on our determination of the planet frequency (Sect. 7), although it is difficult to determine to what extent. Regarding metallicity, the subsample is representative of the overall distribution but does not cover the complete metallicity range.

For these 15 stars, we analyse a total of 1540 radial velocity measurements obtained using the HARPS spectrograph at the 3.6 m ESO telescope (La Silla-Paranal Observatory). All observations were reduced with the HARPS pipeline (version 3.5), where the radial velocities were obtained by performing a cross-correlation of the observed spectrum with a numerical mask (Baranne et al. 1996; Pepe et al. 2002).

For the first measurements (obtained before the present Large Programme and still within the HARPS GTO programme), the exposure times were not long enough to average out the noise coming from stellar oscillations, and the radial velocities had a limiting precision of around 2 ms^-1. The observing strategy was revised later by setting exposure times to 900 s, resulting in a precision of the order of 1 ms^-1. With the start of the Large Programme, in October 2012, we aimed at obtaining more than one spectrum of the star in a given night, separated by several hours.

We use nightly binned data in our analysis, thus focusing on signals with periods longer than one night. The objective of this strategy is to minimize the impact of granulation phenomena on the radial velocities and improve the planetary detection limits (see Dumusque et al. 2011c). We note, however, that this observing strategy has some caveats since it has been designed for and tested in stars with solar metallicities, and its optimality is not guaranteed for metal-poor stars. Also, considering the binned data means we cannot assess the presence of shorter period planets (e.g. Dawson & Fabrycky 2010; Hatzes 2014; Grunblatt et al. 2015).

3. Methodology

As a first step in the analysis of each star, and as is standard when searching for planets in RV data, we computed the generalised Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009) of the radial velocities and searched for significant peaks. The significance of a peak is evaluated by the false alarm probability (FAP) estimated with a bootstrap permutation procedure, first devised by Murdoch et al. (1993; see also Endl et al. (2001), Mortier et al. (2012). The RV measurements and errorbars are (together) randomly redistributed with repetition, keeping the time stamps of observations fixed. For each randomization, a new GLS periodogram is computed, and its maximum power is compared with the original periodogram, thus obtaining the FAP levels. After repeating this process 1000 times, we show the estimated 1% and 0.1% FAP levels.

If any periodogram peak is deemed significant, we proceed with the analysis of the RV data together with other activity proxies to try to determine the nature of the signal. We focus on the FWHM of the HARPS cross-correlation function (CCF) and the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ proxy as the most reliable indicators of activity, even though we have also analysed other indicators. It may also be that a signal is present in the data that was created by a planet with a moderately high eccentricity. The standard GLS periodogram (which searches for sinusoidal signals) is then not the best tool for recovering it. Therefore, for the stars in which the periodogram does not show significant peaks, we still run a planet detection algorithm that compares the evidence for models with 0, 1, and 2 Keplerian signals (see below). For some stars that show evidence of a Keplerian signal, we try to assert its nature further.

As is now ubiquitous in the literature, our model for the radial-velocity observations v_i is described by a Gaussian sampling distribution $\begin{array}{ccc}{\mathit{v}}_{\mathit{i}}\mathrm{|}{\mathit{M}}_{\mathit{n}}\mathrm{~}\mathrm{𝒩}\mathrm{(}{\mathit{V}}_{\mathrm{0}}\mathrm{+}{\mathit{V}}_{\mathit{n}}\hspace{0.17em}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathit{,\theta }\mathrm{\right)}\mathit{,}\hspace{0.17em}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}\mathrm{+}{\mathit{s}}^{\mathrm{2}}\mathrm{)}\mathit{,}& & \end{array}$(1)where M_n represents a model with n planets, V₀ is a constant offset, V_n (t_i,θ) is the radial-velocity shift caused by n orbiting planets, which depends on time t_i and the planets’ orbital parameters θ. The parameter s is an extra “jitter” term, added in quadrature to the total measurement uncertainties σ_i, which allows for additional white noise to be taken into account.

We assume each planet induces a Keplerian signal k (t), such that V_n (t) = ∑ _nk (t). This signal depends on the orbital period P, semi-amplitude K, eccentricity e, argument of the periastron ω, and time of periastron passage T₀ (e.g. Perryman 2014, Chap. 2). When n = 0, we have V_n = 0 by definition.

To infer the number of Keplerians supported by a dataset d = { t_i,v_i,σ_i,i = 1,...,N } , we can calculate the marginal likelihood, or evidence E_n, for each of the models M_n: $\begin{array}{ccc}{\mathit{E}}_{\mathit{n}}\mathrm{=}\mathit{p}\left(\mathit{d}\mathrm{|}{\mathit{M}}_{\mathit{n}}\right)\mathrm{=}{\mathrm{\int }}_{\mathrm{\Theta }}\mathit{p}\left(\mathit{d}\mathrm{|}\mathrm{\Theta }\mathit{,}{\mathit{M}}_{\mathit{n}}\right)\mathit{p}\left(\mathrm{\Theta }\mathrm{|}{\mathit{M}}_{\mathit{n}}\right)\mathrm{d\Theta }\mathit{,}& & \end{array}$(2)where p(d | Θ,M_n) is the likelihood function and p(Θ | M_n) the prior distribution for the parameters of the model. Here, Θ has different dimensions, depending on the value of n and corresponds to the full set of parameters to be marginalised over, Θ = { V₀,s,θ }. The likelihood is given by the product of the terms in Eq. (1)for each data point, $\begin{array}{ccc}\mathit{p}\left(\mathit{d}\mathrm{|}\mathrm{\Theta }\mathit{,}{\mathit{M}}_{\mathit{n}}\right)& \mathrm{=}& {\mathrm{(}\mathrm{2}{\mathit{\pi }}^{\mathrm{)}}}^{\mathrm{-}\mathit{N}\mathit{/}\mathrm{2}}\left[\underset{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}{\mathrm{(}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}\mathrm{+}{\mathit{s}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\right]\\ & & \mathrm{\times }\exp \begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\mathrm{-}\sum _{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathit{N}}\frac{{\left[{\mathit{v}}_{\mathit{i}}\mathrm{-}{\mathit{V}}_{\mathrm{0}}\mathrm{-}{\mathit{V}}_{\mathit{n}}\hspace{0.17em}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathit{,\theta }\mathrm{\right)}\right]}^{\mathrm{2}}}{\mathrm{2}\mathrm{(}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}\mathrm{+}{\mathit{s}}^{\mathrm{2}}\mathrm{)}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathrm{·}\end{array}$(3)The comparison between any pair of models is made by evaluating the associated odds ratio (e.g. Gregory 2010) $\begin{array}{ccc}{\mathrm{𝒪}}_{\mathit{ij}}\mathrm{=}\frac{\mathit{p}\left({\mathit{M}}_{\mathit{i}}\mathrm{|}\mathit{d}\right)}{\mathit{p}\hspace{0.17em}\mathrm{\left(}{\mathit{M}}_{\mathit{j}}\mathrm{\right|}\mathit{d}\mathrm{)}}\mathrm{=}\frac{\mathit{p}\mathrm{(}{\mathit{M}}_{\mathit{i}}\mathrm{)}}{\mathit{p}\hspace{0.17em}\mathrm{\left(}{\mathit{M}}_{\mathit{j}}\mathrm{\right)}}\mathrm{·}\frac{\mathit{p}\left(\mathit{d}\mathrm{|}{\mathit{M}}_{\mathit{i}}\right)}{\mathit{p}\hspace{0.17em}\mathrm{\left(}\mathit{d}\mathrm{\right|}{\mathit{M}}_{\mathit{j}}\mathrm{)}}\mathit{,}& & \end{array}$(4)which simplifies to the ratio of the evidence when there is no prior preference for any model, that is, when $\frac{\mathit{p}\hspace{0.17em}\mathrm{\left(}{\mathit{M}}_{\mathit{i}}\mathrm{\right)}}{\mathit{p}\hspace{0.17em}\mathrm{\left(}{\mathit{M}}_{\mathit{j}}\mathrm{\right)}}\mathrm{=}\mathrm{1}$. If there is a preference for a given model, it can be included directly in Eq. (4). In this work we chose to consider all models equally probable a priori. Kass & Raftery (1995) provide a qualitative scale for interpretating evidence ratio values.

To evaluate the difficult (multidimensional) integral in Eq. (2), we use the nested sampling algorithm (Skilling 2004) implemented in MultiNest (Feroz et al. 2009, 2013). The key idea in nested sampling is to update a set of points, originally sampled from the prior, with new samples that are subject to a hard likelihood constraint. At each iteration, the algorithm progressively moves towards regions of higher likelihood. Each time a new sample is obtained, it is assigned a value X ∈ [ 0,1 ] representing the amount of prior mass estimated to lie at a higher likelihood than that of the discarded point¹. These X values introduce a mapping from the parameter space to the interval [ 0,1 ] or, in other words, they divide the prior volume into a large number of points with equal prior mass. In the new space, the prior becomes a uniform distribution, and the likelihood is a decreasing function of X. Then, the evidence can be computed by simple quadrature. For further details, the reader is referred to Skilling (2004), Mukherjee et al. (2006), Sivia & Skilling (2006), and Feroz et al. (2009).

Nested sampling also provides posterior samples for all parameters as a by-product. The necessary descriptive statistics can then be calculated from these samples. Estimating all the posterior distributions for a particular model is what we mean by fitting a Keplerian in the remainder of the paper. This algorithm has been tested and validated in the analysis of radial-velocity data by Feroz et al. (2011a,b) and Feroz & Hobson (2014).

Another way to approach the problem is to treat n (the number of Keplerians) itself as a free parameter of the model (so that it is part of Θ) and sample from its own posterior distribution, p(n | d), to infer the number of Keplerians present in the data. We use the recent algorithm proposed by Brewer & Donovan (2015) to sample from this distribution. The method uses trans-dimensional birth-death moves, within the context of diffusive nested sampling (Brewer et al. 2011), and it was found to give results comparable to MultiNest, often with a smaller number of likelihood evaluations and less computation time. The same likelihood function and priors were used for both methods.

We furthermore need to specify the priors for each of the parameters. It is through the priors that the Bayesian analysis can take the principle of parsimony and Occam’s razor into account. What is important is not (only) the number of parameters in a model but also the amount of support in their priors. This is the reason to choose uninformative and physically motivated priors.

The form and limits of the priors we use are shown in Table 3. A few of our choices merit further discussion. For the orbital periods, we use a Jeffrey’s prior between 1d and the timespan of each dataset. Apart from HD 56274 and HD 175607, which are analysed separately, none of the stars show clear long-term trends that could suggest the presence of long-period planets. We thus limit our search for periods that are shorter than the timespans, which are in fact very long already. The lower limit is due to our use of nightly binned observations. For the semi-amplitudes, a modified Jeffrey’s prior is used, spanning the range 0−10 ms^-1 and with a “knee” at the mean error bar of each dataset. The choice of upper limit comes from the cuts used when defining the complete metal-poor sample (see Santos et al. 2014), which was stripped of stars showing large RV dispersions². The same arguments apply to the jitter prior. The prior for the orbital eccentricities was suggested by Kipping (2013) after an analysis of the population of known exoplanets. Priors for the other parameters are uniform between what we consider to be sensible limits.

Model selection criteria other than the evidence can be more easily calculated using just the maximum-likelihood value $\hat{\mathit{L}}$ and the total number of parameters in the model, N_par. When presenting our results in Table 4, we include the values of the sample-corrected Akaike information criteria (AIC_C) and the Bayesian information criteria (BIC), calculated as described in Burnham & Anderson (2010): $\begin{array}{ccc}& & {\mathrm{AIC}}_{\mathrm{C}}\mathrm{=}\mathrm{-}\mathrm{2}\ln \hat{\mathit{L}}\mathrm{+}\mathrm{2}\hspace{0.17em}{\mathit{N}}_{\mathrm{par}}\mathrm{+}\frac{\mathrm{2}\mathit{K}\mathrm{(}{\mathit{N}}_{\mathrm{par}}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathit{n}\mathrm{-}{\mathit{N}}_{\mathrm{par}}\mathrm{-}\mathrm{1}}\\ & & \mathrm{BIC}\mathrm{=}\mathrm{-}\mathrm{2}\ln \hat{\mathit{L}}\mathrm{+}{\mathit{N}}_{\mathrm{par}}\ln \mathit{N}\mathit{.}\end{array}$The preferred model when using these criteria is the one that gives the lowest AIC_C and BIC values. For the choice of one model to be significant one can require differences ΔAIC_C and ΔBIC around 10 (Burnham & Anderson 2010).

4. Case-by-case results

In this section we present the individual analysis for each star. Plots with the radial-velocity data and GLS periodograms are shown in Figs. A.1–A.15, and a summary of the model-selection results for the stars without significant periodogram peaks is shown in Table 4.

4.1. Candidate planetary signals

HD 87838.

This star was observed 104 times on 87 different nights. The radial velocities show a dispersion of 2.4 ms^-1 and a mean uncertainty of 1.3 ms^-1. The GLS periodogram shows no significant peaks (Fig. A.1), but the Bayesian analysis finds slight evidence of one Keplerian signal (Table 4).

The strongest peaks in the periodogram are at 988 d, 68 d, and 334 d. The most convincing (highest likelihood) Keplerian fit is a signal with P = 68.34 ± 0.46 d, a semi-amplitude K = 1.9 ± 0.38 ms^-1, and an eccentricity of 0.46 ± 0.16. If caused by a planet, this signal would correspond to a planet with a minimum mass of 11.14 ± 2.47 M_⊕ at an orbital distance of 0.322 ± 0.001 AU. The phase-folded radial-velocities for this solution are shown in Fig. 2.

Analysing the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ and FWHM activity indicators, we do not find strong periodicities close to 68 d (Fig. 3), so it is unlikely that this signal is caused by stellar activity. (Indeed, the periodogram of the FWHM shows a non-significant peak close to the estimated rotation period.) Since the Bayesian analysis only gives marginal evidence of a Keplerian signal (a Bayes factor of 1.75 is “not worth more than a bare mention” in the scale of Kass & Raftery 1995), we consider that this detection is marginal and needs further confirmation. More data would certainly shed light on its true origin, and we will continue to observe this star in the future.

HD 175607.

This star’s periodogram (Fig. A.2) shows a very significant peak at 29 d, coinciding almost exactly with the estimated rotation period (see Table 1). Some of the longer periods found in the RVs are also present in the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ indicator, but no indication of the 29 d period is found in any of the activity proxies. A full analysis of these data is presented in Mortier et al. (2016) where the authors find evidence for (at least) one Neptune-mass planetary companion. In Sect. 6, we consider the detection limits after having removed this orbital solution.

4.2. Signals caused by stellar activity

The effects of stellar activity on radial-velocity observations can be broadly divided into short-term effects, which are modulated by the stellar rotation, and long-term effects caused by global activity cycles (Baliunas et al. 1995). The activity-induced signals can be diagnosed using activity proxies such as line-profile indicators (Queloz et al. 2001; Dumusque et al. 2011a; Boisse et al. 2011) or the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ values (e.g. Bonfils et al. 2007). One way to correct the radial velocities for the effects of stellar activity is to assume a linear correlation between, for example, the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ activity index and the activity-related radial-velocity variations and to remove this correlation from the RVs (e.g. Melo et al. 2007). Assuming an estimate of the rotation period of the star, one can also subtract sinusoidal functions from the RVs at the rotation period and its harmonics (e.g. Boisse et al. 2011). The following four stars show evidence of these activity-induced signals.

HD21132.

A total of 85 (107 before binning) radial-velocity measurements were obtained for this star, spanning a baseline of over ten years. The average error bar and dispersion are relatively high, probably because of the high effective temperature (implying a higher oscillation and granulation noise, Dumusque et al. 2011b). All the peaks in the periodogram are below the 1% FAP (Fig. A.3). The model selection analysis nevertheless finds positive evidence of one Keplerian signal, with the best solution having a period of 3712 d and a 4.9 ms^-1 amplitude.

The very long period of this solution (the timespan of observations is 3935 d) leads us to believe that it is not of planetary origin. In fact, some of peaks at long periods, seen in the periodogram of the RVs, are also found in the periodogram of the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$. When the RVs are corrected with a linear variation on this indicator, these peaks show a clear decrease in power (Fig. 4). No additional signals were found after applying this correction. We thus conclude that the marginally significant signal found in the data of HD 21132 is probably best explained by a long-term magnetic cycle with a period around ten years.

HD 41248.

This star was the subject of a recent debate regarding the presence of a planetary system composed of two super-Earths. From an analysis of the first ~60 measurements of this dataset, Jenkins et al. (2013b) found evidence of two planets at orbital periods of 18.357 d and 25.648 d. Adding about 160 new radial-velocities, Santos et al. (2014) attributed one of the signals (25 d) to activity and could not recover convincing evidence for the second planet. More recently, Jenkins & Tuomi (2014) re-analysed the new data and again found the two-planet model to be the most probable. No new measurements were made for this star following the last two works, but we include it here for completeness.

The period at 25 d is found to be significant both in the periodogram and from the model selection analysis (with an odds ratio³ of ~6 × 10⁴). This would be very strong evidence in favour of (at least one) planetary companion. However, both the FWHM of the HARPS CCF and the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ show significant peaks at the same period (Fig. 5). In light of this evidence we must consider this signal to be caused by stellar activity (see also Suárez Mascareño et al. 2015).

The 18 d peak is not recovered clearly in the full dataset. In addition, a convincing detection of this planet presupposes a correction for the activity-induced signal. When attempting to perform this correction by removing sinusoids at 25 d and its harmonics, we do not find convincing evidence for the one-planet model. The orbital period of the reported planet, as detected by Jenkins & Tuomi (2014), is also very close to our estimate of the rotation period. It is clear that the currently available data do not allow a firm conclusion about the presence of a planetary system around HD 41248.

HD 56274.

This star was followed closely for over 11 years with 228 observations in 161 nights. The RVs show a clear long-term modulation (Fig. 6), which is also evident from the periodogram (Fig. A.5). This modulation is correlated with the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ indicator (Fig. 7), which shows a clear long-term magnetic cycle. The periodogram of the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ values suggests a period of over ten years for this activity cycle.

From the relations of Noyes et al. (1984a) and Mamajek & Hillenbrand (2008), we estimate the rotation period of HD 56274 to be around 13 d (see Table 1). A peak at this period appears in the periodogram of the RV data, slightly above the 1% FAP line, but not in the observed values of the activity indicators. Interestingly, when removing a cubic polynomial from the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ values (as a correction for the magnetic cycle), a peak at 14.5 d appears in the periodogram of the residuals, although not at a significant level. This suggests that the rotation period of the star is indeed around 14 d, such that this signal is present in both RVs and the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$.

To try to correct the RVs for the long-term activity-induced signal, we removed a linear correlation both with $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ and with FWHM (which also correlates with the RVs). The periodograms of the corrected RVs are shown in Fig. 8. Correcting with $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ removes most of the long-period peaks but leaves the peaks associated with one year aliases in the periodogram. The cause for these peaks is probably the same instrumental effect identified in the case of HD 88725 (see Sect. 4.3). No other significant signals can be found in the residual data. The correction with the FWHM reveals a significant peak at 14 d (close to the estimated rotation period) and two peaks at 62 d and 74 d. With this second correction, however, the long-term periodicities are not completely removed. We conclude that with the available data, and after correcting for the long-term magnetic cycle, it is not possible to confirm the presence of planets orbiting this star.

HD 114076.

The rotation period for this star is estimated to be between 41 and 43 d, depending on the calibration (Table 1). The periodogram (Fig. 9 top panel, and Fig. A.6) shows a cluster of peaks around 40 d and also around 80 d. The strongest peak is at 41.3 d and very close to the 1% FAP.

To correct the RVs from what we believe is an activity-induced signal, a sine function was fitted and subtracted from the data. The best fit period was 41.27 d, and the periodogram of the residuals, shown in Fig. 9 (bottom panel), does not present additional significant signals.

To confirm that this signal is caused by the rotation of the star, we analysed the activity indicators. In particular, the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ shows a clear long-term trend, which might be caused by a magnetic cycle. We tried to remove the long-term variations by fitting both quadratic and cubic polynomials to the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ observations. The residuals of both fits show a clear periodic signal at 44.4 d, which helps to corroborate the activity-induced nature of the 41 d signal in the RVs.

4.3. Sampling and instrumental effects

One additional obstacle in the detection of periodic RV variations due to planets is the presence of spurious signals and periodicities caused by the discrete time sampling of the observations or by instrumental effects. These can stem, amongst other effects, from Earth’s rotation and orbital motion. The following stars show evidence of contaminations of this kind.

HD 22879.

Although the periodogram of this star (Fig. A.7) does not show any significant peaks, the Bayesian analysis finds evidence of one Keplerian signal with a period of ~770 d and an eccentricity of 0.7. This corresponds to the second highest peak in the periodogram.

We find that this signal can, however, be caused by the time sampling of the observations. In the top and middle panels of Fig. 10, we show the periodogram together with the window function of the observations. The window function is calculated as the Fourier transform of the observation times (Roberts et al. 1987; Dawson & Fabrycky 2010), and it provides information on the power that is introduced in the periodogram due to the sampling.

The clean algorithm, introduced for frequency analysis by Roberts et al. (1987), is a well-known method of deconvolving the observed spectrum and the window function, thereby reducing the artefacts introduced by the sampling. Applying the clean algorithm to the RV measurements of HD 22879 results in the spectrum shown in the bottom panel of Fig. 10. The highest peak is at 14 d, which is close to the estimated stellar rotation period. It is not clear if this power is activity-induced, since the activity indicators do not show any clear signal at this period (or any long-term drifts).

We conclude that the Keplerian signal found at 770 d is best explained as originating in the time sampling. It is worth mentioning here that the Bayesian analysis is vulnerable to this kind of signal because it does not consider the information present in the window function at any stage.

HD 79601.

The periodogram of the radial velocities of this star shows one significant peak at ~550 d (Fig. A.8). A Keplerian fit at this period results in a solution with P = 532.85 ± 9.63 d, a semi-amplitude of 1.92 ± 0.28 ms^-1, and an eccentricity of 0.55 ± 0.18. A planet with these parameters would have a mass of 20.4 ± 4.16 M_⊕. The phase-folded radial velocities for this solution (Fig. 11) show that the phase coverage is not ideal. The residuals from this fit do not show any significant periodogram peaks or evidence of any additional signals.

The lack of a good phase coverage of the orbital solution puts the planetary hypothesis in question. From Fig. A.8, it is clear that the measurements were obtained in several series, separated by large time gaps. Looking only at the observations after BJD = 2 454 500, the data show hints of a linear drift. When a linear drift is removed from this subset of the data (and also when a quadratic drift is removed from the full dataset), the peak at 550 d vanishes, and no other significant peaks remain in the periodogram of the residuals. At this point, our data do not allow us to reach any firm conclusions about the nature of the signal observed for HD 79601.

HD 88725.

The radial velocities of this star show a significant periodicity at 368.3 d (Fig. A.9). Because this period is close to one year, we measured the correlation between the radial velocity and the barycentric Earth radial velocity (BERV). The posterior distribution for Spearman’s rank correlation coefficient (which was sampled using Markov chain Monte Carlo (MCMC); see Figueira et al. 2016) has a mean of –0.29, and a 95% credible interval (the highest posterior density interval) is [−0.46,−0.13 ]. This anti-correlation suggests that the signal detected in the RVs might be induced by the orbit of the Earth around the Sun.

Dumusque et al. (2015) recently described the cause for signals of this type to show up in HARPS radial velocities and suggested a method of correcting for this effect: the spectral lines that cross the so-called block stitchings in the HARPS CCD can be identified in the correlation mask used to derive the RVs. We created two different correlation masks, one with only the lines falling close to the block stitchings and one without those lines.

For the case of HD 88725, the periodograms of the RVs derived with the original mask and both custom masks are shown in Fig. 12. The spectral lines that cross the block stitching are the origin of the peak at one year (middle panel), and when those lines are removed from the correlation mask, the peak vanishes from the periodogram (bottom panel). After the RVs are corrected, the periodogram does not show significant peaks, and we do not find evidence of any Keplerian signals.

HD 111777.

In Fig. 13, the periodogram of the 79 measurements of this star is plotted together with the spectral window function (see also Fig. A.10). Although there are significant peaks in the periodogram, it becomes impossible to identify clear peaks produced by physical signals, because of the confusion introduced by the sampling. This does not mean that a long-period signal is not present, only that the highest peak in the periodogram can show up in a different location owing to the sampling.

Since the window function shows only excess power for long periods, and not distinct peaks, the CLEAN algorithm does not provide a good correction in this case. The large gap without observations is what causes these issues. After considering only the better-sampled subset of the data (after BJD = 2 456 000), we do not find any significant signals. Also, after removing a linear trend⁴ from the full dataset, the periodogram of the residuals does not show any significant peaks.

The periodograms of the activity indicators do not help in corroborating the presence of a planetary signal: the FWHM shows a non-significant periodicity close to 730 days, but the highest peak in the periodogram of the $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ is at 16 days. For the case of HD 111777, we conclude that we are currently not able to assess the presence of a signal that might originate in a long-period planetary companion.

HD 224817.

The periodogram of the RVs of this star (Fig. A.11) does not show peaks above the 1% FAP line, but the model-selection analysis finds the model with one Keplerian to be the most probable. On the scale of Kass & Raftery (1995), an odds ratio of 1.25 (Table 4) is “not worth more than a bare mention”.

The highest peak in the periodogram is at ~245 d – the best-fit Keplerian also converges to this solution – but this period is probably caused by the deformation of the spectral lines that cross the block stitchings of the HARPS CCD, as for the case of HD 88725. We noticed this when fitting and removing a sinusoid with a period of one year from the RVs decreased the power of the 245 d peak.

To go on to check that the signal is indeed caused by the instrument, we created different correlation masks both with and without the lines falling close to the block stitchings (Dumusque et al. 2015). The periodograms for the original RVs and those from the two custom masks are shown in Fig. 14. When the spectral lines that cross the block stitchings are removed from the correlation mask, the 245 d peak vanishes from the periodogram (bottom panel). The periodogram of the corrected RVs does not show significant peaks, and we do not find evidence of any Keplerian signals.

The second highest peak in the periodogram of the original radial velocities, around 10 d, is not far from the estimated rotation period of the star (Table 1). A clear (though not significant) peak at this period is present in the periodogram of the FWHM indicator. We conclude that the data available for HD 224817 does not suggest the presence of any planetary companions.

4.4. Stars without identified signals

HD 31128.

This star was followed very closely with a total of 215 radial-velocity measurements on 169 different nights. These numbers already take an outlier observation at BJD = 2 455 304.48 into account, which was removed prior to the analysis. This observation has a very low CCF contrast value, making us suspect some instrumental or observational error at this date.

The RV data show high dispersion though the average uncertainty is also above 3 ms^-1. These can be explained by the high effective temperature and very low metallicity of this star, which hinder the precise calculation of the radial-velocity with the CCF technique. Low-metallicity stars within this temperature range have shallower lines, making the determination of precise RVs more difficult, while the higher temperature introduces higher oscillation and granulation noise levels; Dumusque et al. 2011b.

Despite the large number of measurements, none of the periodogram peaks are significant (Fig. A.12). According to the model selection analysis, there is no evidence of any planetary companions orbiting this star. The constant model is preferred relative to the one-planet model according to the value of the evidence (with an odds ratio of 2.5) and the BIC (with ΔBIC of 5.4). However, the AIC_C selects instead the two-planet model, but in light of our general results (see Sect. 5), we tend to trust the Bayesian criteria.

HD 119173.

This star was observed on 98 different nights (totalling 128 measurements), which cover a total timespan of 8.1 years. No significant periodogram peaks were found (Fig. A.13), and in the model selection analysis, the constant model is preferred over the one-planet model (with odds ratio and ΔBIC = 8.7). We conclude that there is no evidence for planetary signals in these data.

HD 119949.

The periodogram of the 81 measurements, shown in Fig. A.14, does not present any significant peaks. The Bayesian search for Keplerian signals reaches the same conclusions (odds ratio ). The AIC_C again selects a more complex one-planet model with a non-significant ΔAIC_C. We report no clear evidence for planetary companions.

HD 126793.

A total of 96 measurements on 79 different nights were obtained. The periodogram (Fig. A.15) does not show any significant peaks, and in this case as well, the model comparison does not find evidence for Keplerian signals in the data. For these data, all model selection criteria select the constant model as the most probable with , ΔAIC_C= 1.1, ΔBIC = 11.5.

5. Comparison of model selection criteria

For the maximum-likelihood solutions, the AIC_C and BIC can be calculated using the number of free parameters in a given model. They are sometimes used as model-selection criteria and are much easier to calculate than the value of the evidence. In Table 4 we show the AIC_C and BIC values for each star and model, using the maximum-likelihood obtained from all MultiNest samples. We highlight in bold the best model according to each criteria.

As general trends in our results we can note that the AIC_C tends to choose a model with a higher number of Keplerian signals. The BIC is found to be more conservative and tends to agree with the model with the highest evidence. These results suggest that the BIC can be used as a viable approximation to the evidence, given that the (global) maximum of the likelihood function can be found.

6. Detection limits

Many of the stars we analysed in Sect. 4 do not contain significant periodic signals or show signals stemming from activity and sampling contaminations. In this section, we derive upper limits to the planetary signals that can be present and yet still undetectable in these data. For any given star, we determine which planet, as a function of its mass and orbital period, can already be ruled out with our observations.

The detection limits are calculated with a standard procedure: injecting trial circular orbits into the observed data (Cumming et al. 1999; Endl et al. 2001; Zechmeister et al. 2009; Dumusque et al. 2011c; Mayor et al. 2011; Mortier et al. 2012). We explore all periods in the range from 1 d to 500 d and semi-amplitudes from 0 to 10 ms^-1 (using a binary search). For ten linearly spaced phases, we compare the periodogram power of the injected period with the FAP level of that period in the original dataset. If, for all phases, the former is higher we consider the planet to be detected. We convert the semi-amplitude to planetary mass using the stellar masses listed in Table 1.

It is important to note that this method assumes that the original dataset (in which mock planets are injected) only contains uncorrelated noise. For the stars analysed in Sect. 4.4, this assumption is valid. For the stars analysed in Sect. 4.1, we substracted the putative orbital solution from the data before calculating the detection limits. Since correcting for activity-induced signals and artefacts originating in the sampling is much more prone to error, for the stars discussed in Sects. 4.2 and 4.3, we assume that the observed radial-velocities contain only noise. Therefore, for this last group of stars, the detection limits estimates can be taken as conservative.

The 99% detection limits for all 15 stars are shown in Fig. 15 where we highlight the regions corresponding to planets with P< 50 d, and masses in the ranges M = 10−30 M_⊕ and M = 30−100 M_⊕. These can be compared to the same regions considered in (Mayor et al. 2011) for the calculation of the planet ocurrence rate.

Almost uniformly in this sample, we are sensitive to planets that are more massive than 10 M_⊕ up to orbital periods of 50 days. For shorter period planets, with P ~ 10 days and less, the current data would already allow for detecting planets with just a few Earth masses.

7. Discussion and conclusions

From over 100 targets in our metal-poor large programme, we selected those that were observed on more than 75 nights. A homogeneous analysis of the radial-velocity observations was carried out by comparing models with a varying number of Keplerian signals. We used the nested sampling algorithm to estimate the evidence of a given model and compared the results with the posterior distribution of the number of Keplerians n. We also calculated the AIC_C and BIC and analysed the periodograms of the radial velocities and activity indicators.

In some cases,we find a disagreement between the evidence results and the presence of significant periodogram peaks. The Bayesian analysis was able to identify a few Keplerian signals that are not significant in the periodograms. Nevertheless, the values of the odds ratios indicate that these detections are not very significant. The results obtained with MultiNest almost always⁵ agree with the posterior distribution p(n | d), which is reassuring but not surprising given that the same priors and likelihood were used in both algorithms.

Signals induced by stellar activity were detected for HD 21132, HD 56274, and HD 114076, besides the known case of HD 41248. After attempting to correct for these signals, we were not able to recover any additional signs of planetary companions.

For HD 79601, HD 22879, and HD 111777, the time sampling of the observations induces spurious signals in the radial velocities, hindering the detection of planetary signals. We also identified a clear one-year periodicity and its harmonics, caused by instrumental effects, on the radial velocities of HD 88725 and HD 224817. These instrument-induced RV variations are not expected to be present in all stars since each star may have different values of the BERV and different spectral lines crossing the CCD block stitchings (see Dumusque et al. 2015).

HD 87838 shows an interesting signal around 68 d, which can be fitted with a slightly eccentric Keplerian. The model selection does not provide strong evidence for the presence of this planet, but we hope to obtain more data to assert its nature. It is important to note that the value of the evidence (and our estimate) is sensitive to the priors on the parameters of a given model. This is not necessarily a limitation, meaning only that care should be taken when choosing the priors and that they should be stated clearly in any analysis.

As an example of these points, we analysed the data of HD 87838 again with n = 1, considering the prior for the orbital period to be uniform between 65 d and 75 d⁶. The resulting value for the evidence was found to be −209.5, in contrast to the value determined previously, of −213.9 (Table 4). This corresponds to an odds ratio of ~156 against the constant model, which would imply very strong evidence for a Keplerian signal and a confident planet detection.

We performed a sensitivity analysis of our results, namely by considering different priors for the orbital periods and semi-amplitudes. Even if, in some cases, a different model was preferred (that is, yielded a stronger proof), the absolute value of the Bayes factors did not change considerably. This means that our level of confidence in the presence or not of a planetary signal does not strictly depend on the parameter priors.

Even though our working sample is small, it already allows for an analysis of the frequency of planets around these stars. We do not attempt to provide constraints on the planet frequency as a function of stellar metallicity, since our sample does not span the metallicity space completely (Fig. 1).

The probability of obtaining k detections in a sample of size N is given by the binomial distribution, assuming a planet frequency f_p$\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathit{k}\hspace{0.17em}\mathrm{\right|}\hspace{0.17em}\mathit{N,}{\mathit{f}}_{\mathrm{p}}\mathrm{)}\mathrm{=}\frac{\mathit{N}\mathrm{!}}{\mathit{k}\mathrm{!}\mathrm{(}\mathit{N}\mathrm{-}\mathit{k}\mathrm{)}\mathrm{!}}\hspace{0.17em}{\mathit{f}}_{\mathrm{p}}^{\hspace{0.17em}\mathit{k}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{p}}{\mathrm{)}}^{\mathit{N}\mathrm{-}\mathit{k}}\mathit{.}& & \end{array}$(5)This equation can be thought of as a function of f_p for given values of k and N, therefore representing a (un-normalised) posterior distribution for f_p. The number of planets detected in a given sample places a constraint on the probable values of the planet frequency, which can be expressed by the mode and the range that covers 68% of the distribution, for example.

Considering the detection of the planet orbiting HD 175607 in a sample of 15 stars, this procedure gives a constraint of ${\mathit{f}}_{\mathrm{p}}\mathrm{=}\mathrm{6.67}{\mathrm{\%}}_{\mathrm{-}\mathrm{5.58}\mathrm{\%}}^{\mathrm{+}\mathrm{4.16}\mathrm{\%}}$. Taking also the unconfirmed planet around HD 87838 into account, one obtains ${\mathit{f}}_{\mathrm{p}}\mathrm{=}\mathrm{13.33}{\mathrm{\%}}_{\mathrm{-}\mathrm{8.84}\mathrm{\%}}^{\mathrm{+}\mathrm{5.80}\mathrm{\%}}$. This last result agrees with that of Mayor et al. (2011), who found 12.27 ± 2.45% for the combined occurrence rate of planets with orbital periods shorter than 50 days and masses in the range M = 10−100 M_⊕. We nevertheless stress that our constraints on the planet occurrence rate are preliminary, since they are based on a small sample and subject to selection biases.

Recent studies have identified correlations between the properties of the planetary system and the metallicity of the star. Lower metallicity stars have been found to host longer period planets (Beaugé & Nesvorný 2013; Adibekyan et al. 2013) or lower eccentricity orbits (Dawson & Murray-Clay 2013). These results, together with the detection limits calculated in Sect. 6, might explain why few planets have been detected so far in our sample. The search for low-mass planets around metal-poor stars requires not only high-precision radial velocities but also a very long baseline of observations.

Our results show that the detection of small planets in these data is hindered by time sampling problems and activity-induced signals. Careful planning of the observational strategy and detailed analyses of the aliasing structures are needed to mitigate some of these problems, and we must seek ways to make use of all the information contained in the RVs and in all the activity indicators, in an optimal and physically meaningful way.

The analysis of the full sample of 109 metal-poor stars is ongoing, and we expect to provide constraints on the frequency of planets orbiting these stars when the observations of the Large Program are finished.

---

Acknowledgments

We would like to thank the referee, Dr. Silvano Desidera, for a careful reading and insightful comments on the paper that helped to improve its quality. This work was supported by Fundação para a Ciência e a Tecnologia (FCT) through the research grant UID/FIS/04434/2013. We also acknowledge the support from FCT in the form of grant reference PTDC/FIS-AST/1526/2014. J.P.F. acknowledges the support from FCT through the grant reference SFRH/BD/93848/2013 and thanks Brendon Brewer for his assistance with the diffusive nested sampling code. P.F. and N.C.S. acknowledge support by 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 FCT in the form of an exploratory project of reference IF/01037/2013CP1191/CT0001. A.M. received funding from the European Union Seventh Framework Program (FP7/2007-2013) under grant agreement number 313014 (ETAEARTH). A.S. is supported by the European Union under a Marie Curie Intra-European Fellowship for Career Development with reference FP7-PEOPLE-2013-IEF, number 627202.

References

Appendix A: RV time series and periodograms

This appendix contains figures with the radial-velocity time series and respective GLS periodogram for each star.

All Tables

All Figures