# The HARPS search for southern extra-solar planets XXXIX. HD175607, the most metal-poor G dwarf with an orbiting sub-Neptune

A. Mortier, J. P. Faria, N. C. Santos, V. Rajpaul, P. Figueira, I. Boisse, A.

Collier Cameron, X. Dumusque, G. Lo Curto, C. Lovis, et al.

► To cite this version:

A. Mortier, J. P. Faria, N. C. Santos, V. Rajpaul, P. Figueira, et al.. The HARPS search for southern extra-solar planets XXXIX. HD175607, the most metal-poor G dwarf with an orbiting sub-Neptune. Astronomy and Astrophysics - A&A, 2016, 585, 10.1051/0004-6361/201526905. hal-01440134

HAL Id: hal-01440134

https://hal.science/hal-01440134

Submitted on 8 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The HARPS search for southern extra-solar planets*,**

## XXXIX. HD 175607, the most metal-poor G dwarf with an orbiting sub-Neptune

A. Mortier<fnref n="1" />, J. P. Faria<fnref n="2" />,<fnref n="3" />, N.C. Santos<fnref n="2" />,<fnref n="3" />, V. Rajpaul<fnref n="4" />, P. Figueira<fnref n="2" />, I. Boisse<fnref n="5" />, A. Collier Cameron<fnref n="1" />, X. Dumusque<fnref n="6" />, G. Lo Curto<fnref n="7" />, C. Lovis<fnref n="8" />, M. Mayor<fnref n="8" />, C. Melo<fnref n="7" />, F. Pepe<fnref n="8" />, D. Queloz<fnref n="8" />,<fnref n="9" />, A. Santerne<fnref n="2" />, D. Ségransan<fnref n="8" />, S. G. Sousa<fnref n="2" />, A. Sozzetti<fnref n="10" />, and S. Udry<fnref n="8" />

Received 6 July 2015 / Accepted 2 November 2015

## ABSTRACT

**Context.** The presence of a small-mass planet ($M_p < 0.1 M_{Jup}$) seems, to date, not to depend on metallicity, however, theoretical simulations have shown that stars with subsolar metallicities may be favoured for harbouring smaller planets. A large, dedicated survey of metal-poor stars with the HARPS spectrograph has thus been carried out to search for Neptunes and super-Earths.

**Aims.** In this paper, we present the analysis of HD 175607, an old G6 star with metallicity [Fe/H] = -0.62. We gathered 119 radial velocity measurements in 110 nights over a time span of more than nine years.

**Methods.** The radial velocities were analysed using Lomb-Scargle periodograms, a genetic algorithm, a Markov chain Monte Carlo analysis, and a Gaussian processes analysis. The spectra were also used to derive stellar properties. Several activity indicators were analysed to study the effect of stellar activity on the radial velocities.

**Results.** We find evidence for the presence of a small Neptune-mass planet ($M_p \sin i = 8.98 \pm 1.10 M_\odot$) orbiting this star with an orbital period $P = 29.01 \pm 0.02$ days in a slightly eccentric orbit ($e = 0.11 \pm 0.08$). The period of this Neptune is close to the estimated rotational period of the star. However, from a detailed analysis of the radial velocities together with the stellar activity, we conclude that the best explanation of the signal is indeed the presence of a planetary companion rather than stellar related. An additional longer period signal ($P \sim 1400$ d) is present in the data, for which more measurements are needed to constrain its nature and its properties.

**Conclusions.** HD 175607 is the most metal-poor FGK dwarf with a detected low-mass planet amongst the currently known planet hosts. This discovery may thus have important consequences for planet formation and evolution theories.

**Key words.** planetary systems – stars: individual: HD 175607 – techniques: radial velocities – stars: solar-type – stars: activity – stars: abundances

## 1. Introduction

Very early after the first exoplanets were discovered, it was suggested that stars with a higher metallicity have a higher probability of hosting a Jupiter-like planet than stars with lower metallicity (Gonzalez 1997). This result was confirmed in a number of subsequent studies (e.g. Santos et al. 2001; Fischer & Valenti 2005; Johnson et al. 2010; Mortier et al. 2013a). Taken at face value, it favours planet formation theories based on the core-accretion model (e.g. Pollack et al. 1996; Mordasini et al. 2009, 2012). According to this model, dust and grains coagulate to

form planetesimals and combine to make larger cores and thus planets. Metal-rich stars and disks can form these cores more quickly, so they have time to accrete gas before the disk dissolves resulting in more gas giants around metal-rich stars.

For lower-mass planets ($M_p < 0.1 M_{Jup}$), such as Neptunes and (super-)Earths, the same correlation is not observed and the planet occurrence rate even appears to be independent of the host-star metallicity (e.g. Udry & Santos 2007; Sousa et al. 2011b; Buchhave & Latham 2015). This is also in agreement with core-accretion theories; see, however, Adibekyan et al. (2012b) or Wang & Fischer (2015). Planet synthesis simulations based on the theories of core-accretion and planet migration showed that the correlation may even be reversed in the case of Earth-sized planets where stars with subsolar metallicities are favoured for harbouring an Earth-sized planet (Mordasini et al. 2012).

For these reasons, a sample of 109 metal-poor stars was chosen for an extensive radial velocity survey with the HARPS

<footnote n="1">¹ SUPA, School of Physics and Astronomy, University of St. Andrews, St. Andrews KY16 9SS, UK
e-mail: am352@st-andrews.ac.uk</footnote>

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

<footnote n="3">³ Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, 4099-002 Porto, Portugal</footnote>

<footnote n="4">⁴ Sub-department of Astrophysics, Department of Physics, University of Oxford, Oxford OX1 3RH, UK</footnote>

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

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

<footnote n="7">⁷ European Southern Observatory, 19001 Casilla, Santiago, Chile</footnote>

<footnote n="8">⁸ Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, 1290 Sauverny, Switzerland</footnote>

<footnote n="9">⁹ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK</footnote>

<footnote n="10">¹⁰ INAF-Osservatorio Astrofisico di Torino, via Osservatorio 20, 10025 Pino Torinese, Italy</footnote>

<footnote>* Based on observations taken with the HARPS spectrograph (ESO 3.6-m telescope at La Silla) under programmes 072.C-0488(E), 082.C-0212(B), 085.C-0063(A), 086.C-0284(A), and 190.C-0027(A).</footnote>

<footnote>** Radial velocity and stellar activity data are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via</footnote>

<footnote>http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/585/A135</footnote>
spectrograph (Mayor et al. 2003) to search for Neptunes and (super-)Earths (Santos et al. 2014). The targets in this survey are bright, chromospherically quiet FGK dwarfs with metallicities between −2.0 and −0.4 dex. More details about this programme can be found in Santos et al. (2014).

To this date, no low-mass planets have been detected in this metal-poor sample, although there is a debate over one star, HD 41248, that shows clear signs of radial velocity variability. Jenkins et al. (2013) reported on the existence of two planets orbiting this star, close to the 7:5 mean motion resonance. However, using the extended dataset coming from our large programme, these planets could not be confirmed (Santos et al. 2014). One of the signals can clearly be seen in the activity indicators and is thought to be due to the stellar rotation and stellar spots on the surface of the star. The other signal could not be detected any more in an extended dataset and may have shown up as a result of the time sampling of the data or as a signature of differential rotation (though Jenkins & Tuomi 2014, reported that the signals are coherent over time).

This paper reports on the presence of at least one Neptune around one of the stars of the metal-poor HARPS survey, HD 175607. In Sect. 2 we describe the observations made. Section 3 presents the stellar properties. We analyse the stellar activity in Sect. 4 and the radial velocities in Sect. 5. We discuss our findings in Sect. 6.

## 2. Observations

HD 175607 was observed with the HARPS spectrograph on the 3.6-m telescope at La Silla Observatory. A total of 119 spectra over 110 nights were taken between July 2004 and October 2013 under different observing programmes<fnref n="1" />. Most spectra were observed with an exposure time of 15 min. This is done to average out noise (signals) coming from short-term stellar oscillations (e.g. Santos et al. 2004). When the large programme started in October 2012, if possible, we tried to obtain two spectra separated by several hours in one given night to reduce granulation effects, following the optimised observational strategies from Dumusque et al. (2011). Since the signals analysed in this work are on much longer timescales, we then averaged over these two measurements per night. The spectra have a mean signal-to-noise ratio of 104 around 6200 Å.

Radial velocities (RVs) were homogeneously derived using the HARPS Data Reduction Software (DRS). This pipeline cross-correlates the observed spectra with a mask representing a G8 dwarf (the spectral type of HD 175607 is G6V). By fitting a Gaussian to the cross-correlation function (CCF), the value and uncertainty of the RV is determined (e.g. Baranne et al. 1996; Pepe et al. 2002). We end up with 110 precise RV measurements with a mean error bar of 0.95 m s⁻<fnref n="1" />, including photon, calibration, and instrumental noise. This mean error bar is slightly lower than the average error bar of all the stars in our sample. The data are taken over a time span of 3390 days (i.e. 9 years and 3 months).

From the DRS, we also get measurements for different stellar activity indicators: full width at half maximum (FWHM) of the CCF, line bisector inverse slope (BIS), contrast from the CCF, chromospheric activity indicator log R'HK from the CaII H&K lines, Hα index<fnref n="2" />. Error bars for the FWHM, BIS,

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Top to bottom: time series of the radial velocities, FWHM, BIS, contrast, log R'HK, and Hα index (the mean value is subtracted for the RVs and FWHM).</figcaption></figure>

and contrast were scaled from the radial velocity error, following Santerne et al. (2015). Figure 1 shows the radial velocity time series, together with the time series of all these indicators.

## 3. Stellar properties

HD 175607 is a bright dwarf star of spectral type G6. It is located at a distance of 45.27 pc from the Sun, according to the new HIPPARCOS reduction (van Leeuwen 2007). All relevant stellar parameters can be found in Table 1.

The stellar atmospheric parameters, effective temperature, surface gravity, and metallicity have been derived by a spectroscopic line analysis on a spectrum resulting from the sum of five individual HARPS spectra, with a total signal-to-noise ratio of 246.40 (Sousa et al. 2011a). Equivalent widths of iron lines (Fe I and Fe II) were automatically determined. These were then used, along with a grid of ATLAS plane-parallel model atmospheres (Kurucz 1993), to determine the atmospheric parameters, assuming local thermodynamic equilibrium in the MOOG code<fnref n="3" /> (Sneden 1973). More details on the method are found in Sousa et al. (2011a) and references therein.

They found a temperature of 5392 ± 17 K. Casagrande et al. (2011) used photometry to derive stellar parameters and obtained a slightly hotter temperature of 5521 K. Given the known issues with the spectroscopic derivation of the surface gravity (e.g. Torres et al. 2012; Mortier et al. 2013b), we corrected the

<footnote n="1">¹ It was first part of a GTO run, then part of three smaller, metal-poor programmes and eventually part of the large programme.</footnote>

<footnote n="2">² The FWHM and contrast were corrected with a second-degree polynomial to account for the telescope losing focus over time.</footnote>

<footnote n="3">³ http://www.as.utexas.edu/~chris/moog.html</footnote>
<table><caption>**Table 1.** Stellar parameters for HD 175607.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Note</th></tr></thead><tbody><tr><td>RA [h m s]</td><td>19 01 05.49</td><td>(1)</td></tr><tr><td>Dec [d m s]</td><td>-66 11 33.65</td><td>(1)</td></tr><tr><td>Spectral type</td><td>G6V</td><td></td></tr><tr><td>m<sub>v</sub></td><td>8.61</td><td></td></tr><tr><td>B-V</td><td>0.70</td><td></td></tr><tr><td>Parallax [mas]</td><td>22.09 ± 1.01</td><td>(1)</td></tr><tr><td>Distance [pc]</td><td>45.27 ± 2.07</td><td></td></tr><tr><td>T<sub>eff</sub> [K]</td><td>5392 ± 17</td><td>(2)</td></tr><tr><td>log g</td><td>4.64 ± 0.03</td><td>(2)</td></tr><tr><td>[Fe/H]</td><td>-0.62 ± 0.01</td><td>(2)</td></tr><tr><td>[α/Fe]</td><td>0.26</td><td>(3)</td></tr><tr><td>Mass [M<sub>☉</sub>]</td><td>0.74 ± 0.05</td><td>(4)</td></tr><tr><td>Radius [R<sub>☉</sub>]</td><td>0.71 ± 0.03</td><td>(4)</td></tr><tr><td>Mass [M<sub>☉</sub>]</td><td>0.71 ± 0.01</td><td>(5)</td></tr><tr><td>Radius [R<sub>☉</sub>]</td><td>0.70 ± 0.01</td><td>(5)</td></tr><tr><td>Age [Gyr]</td><td>10.32 ± 1.58</td><td>(5)</td></tr><tr><td>⟨log R'<sub>HK</sub>⟩</td><td>-4.92</td><td></td></tr><tr><td>P<sub>Rot</sub> [days]</td><td>28.95 ± 0.33</td><td>(6)</td></tr><tr><td>P<sub>Rot</sub> [days]</td><td>29.68 ± 0.47</td><td>(7)</td></tr><tr><td>v sin i [km s<sup>-1</sup>]</td><td>0.9</td><td>(8)</td></tr><tr><td>v sin i [km s<sup>-1</sup>]</td><td>1.31</td><td>(9)</td></tr></tbody></table>

**References.** (1) van Leeuwen (2007); (2) Sousa et al. (2011a), with the surface gravity corrected following Mortier et al. (2014); (3) Adibekyan et al. (2012c); (4) using the Torres et al. (2010) calibration; (5) Bayesian estimation (da Silva et al. 2006) using the PARSEC isochrones (Bressan et al. 2012); (6) using the empirical relationships from Noyes et al. (1984, their Eqs. (3) and (4)); (7) using the empirical relationship from Mamajek & Hillenbrand (2008, their Eq. (5)); (8) Glebocki & Gnacinski (2005); (9) using the recipe of Santos et al. (2002), adapted to the HARPS CCF.

surface gravity from Sousa et al. (2011a) to a more accurate value with the formula provided in Mortier et al. (2014). The spectroscopic metallicity of $-0.62 \pm 0.01$ shows that this star is indeed metal poor, although within the metal-poor survey, it belongs to the more metal-rich half of the sample. The presented errors are precision errors, intrinsic to the spectroscopic method we used, and are very small. A discussion on the systematic errors of our method can be found in Sousa et al. (2011a), their Sect. 3.1. For effective temperature, a systematic error of 60 K is quoted while for metallicity, they quote a systematic error of 0.04 dex.

Adibekyan et al. (2012c) calculated the chemical abundances of this star and found that it is alpha-enhanced ($[\alpha/Fe] = 0.26$). Kinematically this star would belong to the thin disk, or transitioning between the thin and thick disk (Adibekyan et al. 2012c). The alpha-enhancement could hint that this star is more likely to be a planet host since Adibekyan et al. (2012a) found in the HARPS GTO and Kepler samples that iron-poor planet hosts (in all mass regimes) are alpha-enhanced, while single iron-poor stars show no enhancement in other metals.

Stellar masses and radii were derived using two methods. First, to maintain homogeneity with the online catalogue for stellar parameters of planet hosts (SWEET-Cat<fnref n="4" /> – Santos et al. 2013), we used the corrected calibration formulae of Torres et al. (2010)<fnref n="5" />. This gives us a stellar mass of $0.74 \pm 0.05 M_{\odot}$ and a stellar radius of $0.71 \pm 0.03 R_{\odot}$. Second, we also used a Bayesian estimation of stellar parameters (da Silva et al. 2006) through

their web interface<fnref n="6" />. For this, we used the apparent V magnitude, the HIPPARCOS parallax, the effective temperature and metallicity from the spectroscopic analysis, and the PARSEC isochrones (Bressan et al. 2012). From the models, we obtain a stellar mass of $0.71 \pm 0.01 M_{\odot}$ and a stellar radius of $0.70 \pm 0.01 R_{\odot}$ which are comparable with the results from the calibration formulae. Using the same input and through the same web interface for the Bayesian isochrone fitting (da Silva et al. 2006; Bressan et al. 2012), we also get an estimate for the stellar age (10.32 Gyr) that makes it a fairly old star. It also returns a value for the surface gravity, $4.57 \pm 0.01$, which is close to the corrected spectroscopic value. Since the isochronal stellar mass value is more precise, we use that value for the duration of this paper.

HD 175607 is a slowly rotating star. Glebocki & Gnacinski (2005) report a value for the projected rotational velocity $v \sin i = 0.9$ km s⁻¹. Following a similar recipe in Santos et al. (2002), we used the B-V colour and the mean FWHM of all 110 measurements to obtain an estimate of $v \sin i = 1.31$ km s⁻¹. We get an estimate for the rotational period with the empirical relationships of Noyes et al. (1984, their Eqs. (3) and (4)) or Mamajek & Hillenbrand (2008, their Eq. (5)) via the chromospheric activity indicator $log R'_{HK}$. The weighted mean value of $log R'_{HK}$ is -4.92 over all 110 measurements. Combining this with the B-V, we obtain an estimated rotational period of about 29 days. This is just an estimate resulting from calibrations and the true rotational period is not known. All stellar parameters are in Table 1.

## 4. Activity analysis

Even in relatively inactive stars, radial velocity variations can be induced by stellar mechanisms other than orbiting planets, such as intrinsic stellar variations coming from stellar spots and/or faculae on the surface of the star (e.g. Boisse et al. 2011; Haywood et al. 2014; Santos et al. 2014; Robertson et al. 2015a). It is thus important that we study the stellar activity to be able to distinguish between RV signals coming from a planet and those from the star itself. As mentioned in Sect. 2, we have measurements of different activity indicators. If periodic variations in the RV signal were also present in one or more of these activity indicators, that could mean that the RV variation is activity induced rather than planet induced.

Figure 2 shows the General Lomb-Scargle (GLS) periodograms (Zechmeister & Kürster 2009) from the RV and the four main activity indicators provided by the HARPS DRS pipeline: FWHM, BIS, contrast, and $log R'_{HK}$. A bootstrapping method is used to determine the 1% false alarm probability (FAP, for details see Mortier et al. 2012). There are four significant peaks in the RV periodogram (see more in Sect. 5). The most significant peak is seen around 29 days, which is the same as the estimated rotational period from the activity level (see previous section). Studying the activity indicators as proxies of stellar activity is thus even more important in this specific case.

When we look at the GLS periodograms of the CCF parameters (FWHM, BIS, contrast), none of the peaks seen in the RV periodogram are observed. In fact, none of these indicators show strong periodical patterns. There is some short-term (3–5 days), non-significant variation in the BIS, but none of these signals could be found in the RV periodogram. In fact, the estimated rotational period is not clear from these indicators. The periodogram of the $H\alpha$ index shows significant peaks at 24.5 and

<footnote n="4">⁴ https://www.astro.up.pt/resources/sweet-cat/</footnote>

<footnote n="5">⁵ See Santos et al. (2013) for details on the correction.</footnote>

<footnote n="6">⁶ http://stev.oapd.inaf.it/cgi-bin/param</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 4.** Pre-whitening the radial velocities using GLS periodograms. Top panel: raw RVs. Middle panel: residual RVs after subtracting the best-fitted signal at 29.03 days. Bottom panel: residual RVs after subtracting the best-fitted signals at 29.03 and 713.65 days. The horizontal black dashed lines represent the 1% FAP.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 5.** BGLS periodogram of the raw RVs. The highest peak has been normalised to 100% probability. This shows that the period at 29 days is ~10¹⁰ times more probable than the period at 713 days.</figcaption></figure>

Figure 4 shows the results of this data pre-whitening. In the GLS periodogram of the original RV data, the strongest peak can be seen at 29.03 days. After removing this period from the data, we find that the peak at around 18 days also disappeared. This hints at the fact that this period could be associated with the monthly alias of the 29-day period. The long-term periods are still significant, the highest of which is at 713.65 days. After subtracting this solution from the data, the other long-term period peak, at around 1400 days, also vanished. In the residual periodogram, the highest peak is now around 21 days, but this is not significant and at the level of the noise. We thus find two significant periodicities in the data: one at 29 or 18 days and one at 713 or 1400 days.

To assess the relative probability of the peaks in the periodograms, we used the Bayesian Generalized Lomb-Scargle Periodogram (BGLS) as described in Mortier et al. (2015)<fnref n="8" />. Figure 5 shows this BGLS where the probability of the highest peak (at 29 days) is normalised at 100%. This analysis shows that the period at 29 days is ~10¹⁰ times more probable than the period at 713 days. The periods have a median relative probability of $\langle P \rangle \sim 10^{-55}\%$, so it is highly probable that the observed periodicities are associated with real periodic signals in the data.

A multi-frequency periodogram (e.g. Baluev 2013) can also be used to detect multiple periodicities in the data and assess their significance. We used FREDEC (for details see Baluev 2013). We looked for all tuples of significant periodicities in the data with periods between 2 and 10 000 days. We find several significant possibilities for a two-period solution. The strongest solution, with a tuple FAP of 1.66 × 10⁻⁷% (and the lowest χ²-value), is found for the combination of periods at 29 and 706 days. All combinations are made up of a short period (29 or 18 days) and a longer period (700 or 1400 days).

## 5.2. Statistical analysis

Periodograms are tools to check which sinusoidal periodicities are present in a dataset. They are important for a first interpretation of the data, but to get a more robust fit of the data and to assess error bars on the parameters, other methods should be employed. We used a genetic algorithm, an MCMC algorithm, and a Gaussian processes (GP) analysis.

### 5.2.1. Genetic algorithm

Initially, we ran a genetic algorithm using yorbit (Ségransan et al. 2011). This algorithm uses a population of 4800 genomes where each genome (defined by frequency, phase, and eccentricity) corresponds to a planetary system. We ran the genetic algorithm twice, once assuming one planet and once assuming two planets. No conditions were set on any of the parameters. A restriction on the eccentricity is automatically set to avoid the planet colliding with the star. Initial starting positions are chosen based on the peaks in the periodogram. The evolution ended when more than 95% of the population converged within 3 sigma of the best solution.

The one planet model ended with a population of planets with periods $P = 29.022 \pm 0.014$ days and eccentricities $e = 0.148 \pm 0.084$. For the two planet model, we again find this planet around 29 days ($P = 29.007 \pm 0.014$ and $e = 0.091 \pm 0.037$). The second planet, however, is not that well constrained. Similar to the frequency analysis carried out in Sect. 5.1, the algorithm finds two types of solutions with periods equally distributed around 700 or 1400 days. The longer period would also be slightly more eccentric, but all solutions have an eccentricity lower than 0.6.

### 5.2.2. Markov chain Monte Carlo

The solutions explored by the genetic algorithm do not provide a reliable statistical population from which to perform inference. It only provides a small parameter space that could be a good starting point for more robust fitting methods such as sampling from the posterior probability through MCMC. This alternative method allows the posterior distribution of each parameter to be inferred. We employ the following model for the RVs:

$$RV(t) = \gamma + \sum_{i} K_{i}[\cos(\omega_{i} + v(t, e_{i}, T_{0,i}, P_{i})) + e_{i} \cos \omega_{i}], \quad (1)$$

where $\gamma$ is the constant systemic velocity, $K$ the RV amplitude, $e$ the eccentricity, $\omega$ the argument of periapse, and $v(t)$ the true anomaly. A sum is taken over all possible Keplerian signals. The true anomaly is a function of time, eccentricity, the period $P$, and the time of periastron passage $T_0$. It is defined as

$$\tan \frac{v}{2} = \sqrt{\frac{1+e}{1-e} \tan \frac{E}{2}}, \quad (2)$$

<footnote n="8">⁸ https://www.astro.up.pt/exoearths/tools.html</footnote>
<table><caption>**Table 2.** Priors for the MCMC procedure.</caption><thead><tr><th>Parameter</th><th>Prior</th><th>Limits</th></tr></thead><tbody><tr><td>γ [m/s]</td><td>Uniform</td><td>[-91 906.42, -91 871.82]</td></tr><tr><td>jitter</td><td>Mod. Jeffreys*</td><td>[0.0, 5.0]</td></tr><tr><td>K<sub>1</sub> [m/s]</td><td>Mod. Jeffreys*</td><td>[0.0, 10.0]</td></tr><tr><td>P<sub>1</sub> [d]</td><td>Jeffreys</td><td>[27.0, 32.0]</td></tr><tr><td>e<sub>1</sub></td><td>Uniform</td><td>[0, 1]</td></tr><tr><td>ω<sub>1</sub></td><td>Uniform</td><td></td></tr><tr><td>T<sub>0</sub>, 1 [JDB]</td><td>Uniform</td><td>[2 455 500.0, 2 455 560.0]</td></tr><tr><td>K<sub>2</sub> [m/s]</td><td>Mod. Jeffreys</td><td>[0.0, 10.0]</td></tr><tr><td>P<sub>2</sub> [d]</td><td>Jeffreys</td><td>[200.0, 2000.0]</td></tr><tr><td>e<sub>2</sub></td><td>Uniform</td><td>[0, 1]</td></tr><tr><td>ω<sub>2</sub></td><td>Uniform</td><td></td></tr><tr><td>T<sub>0</sub>, 2 [JDB]</td><td>Uniform</td><td>[2 454 300.0, 2 456 300.0]</td></tr></tbody></table>

**Notes.** (*) Knee for the modified Jeffreys prior is taken to be the mean error bar $\bar{\sigma}_i$.

with *E* the eccentric anomaly, which in turn can be found by solving Kepler's equation

$$E - e \sin E = 2\pi \frac{t - T_0}{P}. \qquad (3)$$

An additional jitter term is quadratically added to the error bars to incorporate the underestimation of these RV error bars and account for any additional noise present in the data. The final Gaussian likelihood function is

$$p(D|\theta) = \prod_{i=1}^{N} \left[ \frac{1}{\sqrt{2\pi(\sigma_i^2 + \text{jitter}^2)}} \exp\left(-\frac{[y_i - \text{RV}(t_i)]^2}{\sigma_i^2 + \text{jitter}^2}\right) \right], \quad (4)$$

where *N* is the number of datapoints, *θ* the set of parameters in the RV model, and *D* the data. This data consists of the times of observation *t*ᵢ, the measured radial velocities *y*ᵢ, and the estimated error bars $\bar{\sigma}_i$.

The parameter set *θ* has a prior distribution $p(\theta)$. We assume that all parameters are independent so that the total prior distribution can be expressed as the product of the prior distributions of each parameter. We take uniform priors for *γ*, *T*₀, *e*, and *ω*, a Jeffreys prior for the period *P*, and a modified Jeffreys prior for the amplitude *K* and the jitter term (as in Gregory 2005). The knee for this modified Jeffreys prior is taken to be the mean error bar $\bar{\sigma}_i$. All priors used for the MCMC are listed in Table 2.

Using Bayes' theorem, the posterior density is then expressed as

$$p(\theta|D) = \frac{p(\theta)p(D|\theta)}{p(D)}. \qquad (5)$$

Herein, the data probability $p(D)$ is seen as a normalisation constant and is kept at 1 for the MCMC procedure. We calculate $p(D)$ later to compare the different models.

In the MCMC routine, we calculate the natural logarithm of the posterior probability density. Furthermore, we perform a coordinate transformation and use $\sqrt{e} \cos(\omega)$ and $\sqrt{e} \sin(\omega)$ instead of *e* and *ω* (see e.g. Ford 2006). This can be done easily since the Jacobian factor for this transformation is 1. To run the MCMC, we use *emcee* (Foreman-Mackey et al. 2013), a Python code that implements an affine invariant MCMC ensemble sampler (Goodman & Weare 2010). An initial guess for the walkers is randomly chosen inside the final population of the genetic algorithm. We used 700 walkers with 2000 steps. We allow for

<figure><img src="image_5.png" /><figcaption>**Fig. 6.** Full orbit, using the MCMC results of a two Keplerian model. Top panel: relative RVs versus time; bottom panel: residuals.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 7.** Phased orbits, using the MCMC results of a two Keplerian model. Top panel: 29 d signal; bottom panel: 1400d signal.</figcaption></figure>

a burn-in period, which is chosen to be ten times the maximum autocorrelation time of the resulting walkers. Afterwards, we additionally perform a declustering method to remove the walkers with significantly lower posterior probabilities (as in Hou et al. 2012). This removes the walkers that got stuck inside local maxima.

Results for the one- and two-Keplerian models are listed in Table 3. The best fit for the two Keplerian model is shown in Figs. 6 and 7. A periodogram of the residuals reveals just noise, so it was chosen not to run a model with three Keplerians.

In order to compare the two models statistically, one would want to assess the Bayes factor, i.e. the ratio of the model evidence. In the case of two models $M_1$ and $M_2$, each with the parameter set $\theta_1$ and $\theta_2$, the Bayes factor to assess model two over model one is expressed as:

$$B_{21} = \frac{P(D|M_2)}{P(D|M_1)} = \frac{\int P(\theta_2|M_2)P(D|\theta_2, M_2)d\theta_2}{\int P(\theta_1|M_1)P(D|\theta_1, M_1)d\theta_1}. \qquad (6)$$

Calculating these integrals over the complete parameter space is tricky. However, there are ways to solve it. The *emcee* package provides a parallel-tempering ensemble sampler that can be used to estimate this integral. It makes use of thermodynamic integration as described in Goggans & Chi (2004). For a more detailed calculation, see Appendix A. We applied this formalism, using 20 different temperatures (each one increasing with $\sqrt{2}$) with 200 walkers each. As a burn-in, we used 1000 steps and then an
<table><caption>**Table 3.** Planetary parameters from the MCMC and GP fitting procedures.</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="3">MCMC - 1 planet</th><th colspan="3">MCMC - 2 planets</th><th colspan="2">GP - 1 planet</th></tr><tr><th>median</th><th>+σ</th><th>-σ</th><th>median</th><th>+σ</th><th>-σ</th><th>MAP value</th><th>±σ</th></tr></thead><tbody><tr><td>γ [m/s]</td><td>-91 889.69</td><td>0.22</td><td>0.22</td><td>-91 890.41</td><td>0.29</td><td>0.28</td><td>-918 89</td><td>1</td></tr><tr><td>K<sub>1</sub> [m/s]</td><td>2.21</td><td>0.33</td><td>0.33</td><td>2.37</td><td>0.29</td><td>0.30</td><td>1.8</td><td>0.4</td></tr><tr><td>P<sub>1</sub> [d]</td><td>29.03</td><td>0.03</td><td>0.03</td><td>29.01</td><td>0.02</td><td>0.02</td><td>29.0</td><td>0.2</td></tr><tr><td>m<sub>1</sub> sin i [M<sub>⊕</sub>]</td><td>8.26</td><td>1.25</td><td>1.25</td><td>8.98</td><td>1.10</td><td>1.10</td><td>6.7</td><td>1.5</td></tr><tr><td>e<sub>1</sub></td><td>0.16</td><td>0.14</td><td>0.11</td><td>0.11</td><td>0.09</td><td>0.07</td><td>0.17</td><td>0.10</td></tr><tr><td>ω<sub>1</sub></td><td>0.55π</td><td>0.36π</td><td>0.34π</td><td>0.79π</td><td>0.29π</td><td>0.29π</td><td>1.0π</td><td>0.40π</td></tr><tr><td>T<sub>0, 1</sub> [BJD]</td><td>2 455 528.01</td><td>4.71</td><td>5.16</td><td>2 455 532.17</td><td>4.15</td><td>4.18</td><td>2 453 219</td><td>6</td></tr><tr><td>K<sub>2</sub> [m/s]</td><td>—</td><td>—</td><td>—</td><td>2.86</td><td>0.51</td><td>0.51</td><td>—</td><td>—</td></tr><tr><td>P<sub>2</sub> [d]</td><td>—</td><td>—</td><td>—</td><td>1336.61</td><td>103.27</td><td>45.50</td><td>—</td><td>—</td></tr><tr><td>m<sub>2</sub> sin i [M<sub>⊕</sub>]</td><td>—</td><td>—</td><td>—</td><td>34.97</td><td>6.93</td><td>—</td><td>—</td><td>—</td></tr><tr><td>e<sub>2</sub></td><td>—</td><td>—</td><td>—</td><td>0.42</td><td>0.15</td><td>0.14</td><td>—</td><td>—</td></tr><tr><td>ω<sub>2</sub></td><td>—</td><td>—</td><td>—</td><td>0.08π</td><td>0.10π</td><td>0.09π</td><td>—</td><td>—</td></tr><tr><td>T<sub>0, 2</sub> [BJD]</td><td>—</td><td>—</td><td>—</td><td>2 455 244.26</td><td>63.39</td><td>73.95</td><td>—</td><td>—</td></tr><tr><td>jitter</td><td>2.01</td><td>0.17</td><td>0.19</td><td>1.40</td><td>0.16</td><td>0.17</td><td>—</td><td>—</td></tr><tr><td>P<sub>gp</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>29.9</td><td>0.2</td></tr><tr><td>λ<sub>p</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>0.16</td><td>0.02</td></tr><tr><td>τ [d]</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>67</td><td>11</td></tr></tbody></table>

**Notes.** Errors are the 1σ uncertainties taken from the posterior distributions.

additional 2000 steps for the integral calculation. We find that $B_{21} \sim \exp(15)$, supporting the model with two Keplerians with very strong evidence (e.g. Kass & Raftery 1995).

We emphasize that this evidence is dependent on the chosen priors. Specifically, the prior on the period of the inner planet may be seen as too narrow. We ran tests where the prior on this period is 1 to 100 days. We get comparable results as with the more narrow prior, although the time of periastron (whose prior is then also widened) is less constrained because it is cyclic. Thermodynamic integration with these wider priors gives us a Bayes factor $B_{21} \sim \exp(19)$, even higher than before. We can thus be confident that the strong evidence is not due to our choice of priors.

### 5.2.3. Gaussian processes

Gaussian processes provide a mathematically-tractable and flexible framework for performing Bayesian inference about functions. They are particularly suitable for the joint modelling of deterministic processes (such as signals induced by planets) with stochastic processes of unknown functional forms such as activity signals (Aigrain et al. 2012; Haywood et al. 2014). Despite not knowing the functional form of these stochastic processes, we usually know some of its properties.

Rajpaul et al. (2015), hereafter R15, developed this kind of framework to model RV time series jointly with one or more ancillary activity indicators. This allows the activity component of the RV time series to be constrained and disentangled from planetary components. Their framework treats the underlying stochastic process, giving rise to activity signals in all available observables (RVs and ancillary time series) as being described by a GP, with a suitably-chosen covariance function. They then use physically-motivated and empirical models to link this GP to the observables; with the addition of noise and deterministic components (e.g. dynamical effects for the RVs), all observables can be modelled jointly as GPs, with the ancillary time series thus serving to constrain the activity component of the RVs. They showed their framework can be used to disentangle activity

and planetary signals. This is the found even when the planetary signal is much weaker than the activity signal ($\Delta RV \lesssim 0.5$ m/s) and has a period identical to the activity signal. Since the period of the first signal in the data for HD 175607 is very close to the estimated rotational period of the star, we performed a fit for this signal with the GP framework as described in R15.

The marginal likelihood $\mathcal{L}(\theta, \phi)$ for the data, given a GP model, can be expressed as

$$ \ln[\mathcal{L}(\theta, \phi)] = -\frac{1}{2}r^T K^{-1}r - \frac{1}{2}\ln(\det K) - \frac{N}{2}\ln(2\pi), \quad (7) $$

where $r(t, \theta) = y - m(t, \theta)$ is the vector of residuals of the data after the mean function $m$ has been subtracted and $N$ is the number of datapoints. The free hyper-parameters $\theta$ and $\phi$ can then be varied to maximise $\mathcal{L}$; this process is known as Type-II maximum likelihood, or marginal likelihood maximisation (Gibson et al. 2012). In so doing, we refine vague distributions over many, very different functions, the forms of which are controlled by $\theta$ and $\phi$, to more precise distributions that are focused on functions that best explain our observed data.

We implemented the GP framework exactly as described in R15. In particular, given that we have a physical reason to expect a degree of periodicity in the activity signals (as they are modulated by the periodic rotation of the star), we adopted the following quasi-periodic covariance function for the framework's underlying, activity-driving process. This covariance function was previously considered by Aigrain et al. (2012) to model observed variations in the Sun's total irradiance, and by Haywood et al. (2014) to model correlated noise in the CoRoT-7 data

$$ k(t, t') \propto \exp \left\{ -\frac{\sin^2 \left[ \frac{\pi(t-t')}{P_{gp}} \right]}{2\lambda_p^2} - \frac{(t-t')^2}{2\tau^2} \right\}, \quad (8) $$

where $P_{gp}$ and $\lambda_p$ correspond to the period and length scale of the periodic component of the variations and $\tau$ is an evolutionary timescale. While $\tau$ has units of time, $\lambda_p$ is dimensionless.

For HD 176507, we jointly modelled the $\Delta RV$ (after subtracting a polynomial to exclude longer period variations),
<figure><img src="image_7.png" /><figcaption>**Fig. 8.** GP model MAP fit to the HD 175607 data. The 110 observations in each time series were fit simultaneously, i.e. using a single set of (hyper)parameters. The dots indicate observed data, with estimated errors; solid lines are model posterior means; and shaded regions denote ±σ posterior uncertainty. Residuals are plotted below the observed data and fitted model, but for the sake of clarity, with an arbitrary vertical offset from the main time series.</figcaption></figure>

<table><caption>**Table 4.** Priors for the GP procedure.</caption><thead><tr><th>Parameter</th><th>Prior</th><th>Limits</th></tr></thead><tbody><tr><td>γ [m/s]</td><td>Uniform</td><td>[-91 906.42, -91 871.82]</td></tr><tr><td>K [m/s]</td><td>Mod. Jeffreys*</td><td>[0.0, 10.0]</td></tr><tr><td>P [d]</td><td>Jeffreys</td><td>[27.0, 32.0]</td></tr><tr><td>e</td><td>Uniform</td><td>[0, 1]</td></tr><tr><td>ω<sub>1</sub></td><td>Uniform</td><td>[0, 2π]</td></tr><tr><td>T<sub>0</sub> [JDB]</td><td>Uniform</td><td>[2 453 206.0, 2 453 206.0 + P]</td></tr><tr><td>P<sub>gp</sub></td><td>Uniform</td><td>[1, 100]</td></tr><tr><td>λ<sub>p</sub></td><td>Jeffreys</td><td>[0.01, 100]</td></tr><tr><td>τ</td><td>Jeffreys</td><td>[0.1, 1000]</td></tr></tbody></table>

**Notes.** (*) Knee for the modified Jeffreys prior is taken to be the mean error bar $\bar{\sigma}_i$.

$\log R'_{HK}$ and BIS time series as in R15. We chose not to include the FWHM since FWHM data are noisier than, but often very tightly correlated with $\log R'_{HK}$, and thus often do not contain useful extra information that the other indicators have not yet provided.

Non-informative priors (just as for the MCMC procedure) were placed on all Keplerian orbital parameters (incorporated into the GP's mean function, $m$). These priors and the priors on the hyper-parameters are listed in Table 4. Parameters for the Keplerian orbit are estimated using the MultiNest nested-sampling algorithm (Feroz & Hobson 2008; Feroz et al. 2009, 2013), with the GP hyper-parameters first fixed at their MAP values, as per the computational approximation motivated in Gibson et al. (2012).

Our findings were as follows. When not including a planetary component in the GP's mean function for the ΔRV time series, the MAP value of the hyper-parameter $P_{gp}$ ended up being $29.0 \pm 0.1$ d: because the 29.0-d signal was so significant in the ΔRV time series, the GP was forced to absorb this, whilst all but ignoring and thus failing to fit the ancillary time series.

On the other hand, when including a Keplerian component, the hyper-parameter $P_{gp}$ ended up being $29.9 \pm 0.4$ d, with the 29-d signal being absorbed entirely by the Keplerian component; under this model, the rms of the RV variations absorbed by the GP was reduced to the order of tens of centimetres per second. This is significant because whereas a GP can in principle model an arbitrarily-complex signal arbitrarily well (the key constraint in R15's framework, however, is that the same quasi-periodic GP basis functions must be used to model RV and ancillary time series simultaneously), a Keplerian function is far simpler, and is always be strictly periodic. Therefore, the fact that the simpler, less flexible Keplerian interpretation is favoured by the GP framework indicates that the 29 d signal must have a coherent phase over the entire dataset, strengthening the planetary interpretation of the 29-d signal. The planet parameters we inferred when using the GP framework are presented in Table 3. The evolution timescale for the activity signal is found to be 67 d, slightly more than two rotation periods, as would be expected for this type of star.

We used the sample size-adjusted Akaike information criterion (AIC; Burnham & Anderson 2002) to select between the one-planet vs. no-planet models. The AICc value for the no-planet model was -25.44, and the corresponding value for the one-planet model -33.06, indicating that the planetary explanation was favoured by about a factor of ten.

The MAP fit using the one-planet model is presented in Fig. 8 with a close-up in Fig. 9. After subtracting the one-planet GP model, the residual time series appeared white and normally-distributed, with no significant power on timescales smaller than one year, and with rms 0.95 m/s. This suggests that all of the RV variation (at least on timescales smaller than one year) can be explained fully with the planet + activity model. The $\log R'_{HK}$ and BIS residuals contained no significant periodicities on any timescales.
<figure><img src="image_8.png" /><figcaption>**Fig. 9.** GP model MAP fit to the HD 175607 data, with a close-up view of the region of densest time coverage (57 observations over the course of about seven months). The dots indicate observed data, with estimated errors; solid lines are model posterior means; and the shaded regions denote $\pm\sigma$ posterior uncertainty. Residuals are plotted below the observed data and fitted model, but for the sake of clarity, with an arbitrary vertical offset from the main time series.</figcaption></figure>

# 6. Discussion and conclusion

In this work we analysed the radial velocities of HD 175607, a metal-poor ($[Fe/H] = -0.62$) dwarf star. These radial velocities show a clear periodicity around 29 days and a significant longer period signal. The main question is whether these signals are caused by a planet or rather another phenomenon resulting from the star itself. We discuss each signal below.

## 6.1. Short period signal

The short period signal arises at 29 days. However, the rotational period is also estimated to be around 29 days and the Moon's orbital period is also close to 29 days, so caution is recommended. If this is due to a planet that would make the planet a small Neptune ($M_p \sin i = 8.98 \pm 1.10 M_\oplus$ if the two-planet model is assumed).

Radial velocities can be contaminated by scattered light from the Moon. Specifically, this contamination can produce an additional dip in the CCF. If the Moon's velocity is close to the mean stellar velocity, the two dips are blended, which affects the RV measurement of the star. In the case of HD 175607, the mean stellar velocity is about $-92 \text{ km s}^{-1}$. The Moon orbits the Earth at about $1 \text{ km s}^{-1}$, and the Earth orbits the Sun at about $30 \text{ km s}^{-1}$. Consequently, the additional CCF dip due to scattered moonlight contamination is always going to be equal or more than $60 \text{ km s}^{-1}$ redwards of the stellar CCF. This makes moonlight contamination in the RVs of this star impossible.

We emphasize that even if there would be contamination from the Moon in our RVs, Cunha et al. (2013) showed that for the spectral type and magnitude of HD 175607, the contamination would be around 10 cm/s, which is much lower than the signal seen here. We are thus confident that the 29 day signal is not due to the Moon.

Then remains the question of the rotational period. For several reasons listed below, we think that the signal is indeed best explained as being from a planet rather than activity-related:

- No significant correlations are found with any of the activity indicators provided by the HARPS DRS pipeline, nor with the extra activity indicators we calculated using the code in Figueira et al. (2013). If the signal were to be activity related, one would expect there to be some correlation with at least one of the activity indicators. The lack thereof suggests the signal is planet related.

- The $H\alpha$ index shows significant periodicities around 24 days. It could thus be that the estimated rotational period, coming from the $B-V$ colour and the mean log $R'_{HK}$, is not accurate and the rotational period is closer to 24 days. In this case, the RV signal would not be at the same period of the stellar rotation.

- We have data spanning over nine years with about 4.5 yr of intense datasets. This is of the order of 50 times the lifespan of a typical solar active region. Signals arising from activity are not expected to stay stable over this amount of time for this type of star. Since the period of the signal is still very well constrained, that hints that the signal is stable over time and thus not due to activity.

- In the GP analysis, the red noise is modelled separately from the Keplerian, though both are at similar periods. This analysis thus prefers the presence of a planet despite activity signals at similar periodicities. The planetary mass is lowest when using this model. We think this is because some of the signal's amplitude, swallowed by the GP, is treated as planetary in the other models.

- If a signal is not stable over time, such as one caused by activity, the peak in a periodogram would be variable, depending on the amount of activity on certain times. We tested this and the peak gets always stronger when adding more data.
- As a final test, we wanted to know what the expected periodogram power would be if we inject a noiseless Keplerian signal in the data with similar period, semi-amplitude, and eccentricity as the current signal. We thus injected a sinusoid with the same semi-amplitude and eccentricity but at a period of 21 days. As expected, we see a peak at 21 days. It has about the same power as the 29 d peak. Since we did not add any noise for the 21d signal, this again hints that the 29 d signal is of planetary nature.

There are other known cases where the orbital period is the same as the stellar rotation period, such as CoRoT-11b (Gandolfi et al. 2010) or XO-3b (Hébrard et al. 2008). However, these are all cases of close-in hot Jupiters around fast-rotating stars, where the synchronous planetary orbit may come from tidal locking with the host star (e.g. Lanza et al. 2011; Bolmont et al. 2012). The 29 d period of our mini-Neptune makes it plausible that the planet would have synchronised its host star since timescales for such a synchronisation scale with $(a/R_*)^5 \cdot 1/M_p$ (e.g. Dobbs-Dixon et al. 2004; Brown et al. 2011). The planet could be tidally locked to the star, but there is no way of verifying that without the planetary spin period. There are several discovered planets with periods between 10 and 40 days, which are the typical orbital periods for slowly rotating stars, making it not that unlikely that some of them have periods close to or similar to their estimated stellar rotational periods.

## 6.2. Long period signal

The long period signal is not as well constrained as the shorter period signal. From the MCMC, it was clear that the likelihood of a 1400d signal was much higher than the one from a 700d signal. The latter periods were sampled by the MCMC, but eventually removed in the declustering due to too low likelihood. If due to a planet, this planet would have a period of 1337 days and a minimum mass of about 35 Earth masses (i.e. 0.1 Jupiter masses), making it a large Neptune.

Though not statistically significant, similar long periodicities can be seen in the $\log R'_{HK}$ and contrast. However, after removing the inner planet, there is still no significant correlation between the residual RVs and these indicators, nor did it get stronger. If the longer period signal were due to activity, we would have expected the correlations to arise when removing the shorter period signal.

With the long data span, we cover about 2.5 orbits of ~1400 days. However, given the small number of datapoints in the first half of the dataset, we actually only span one full orbit. Furthermore, there are large gaps without data. We would need more data in order to confirm the nature of this signal and better constrain it in case of a planet. Follow up measurements are planned to resolve this.

## 6.3. Metal-poor survey

This detection is part of a large survey with the HARPS spectrograph for Neptunes around metal-poor FGK dwarfs. HD 175607b is the first Neptune-mass planet discovered in this survey. Despite the low metallicity of the host star ([Fe/H] = -0.62), it still belongs to the more metal-rich part of the sample. The metallicities for the entire sample range from -1.5 to -0.05 dex (Santos et al. 2014). In a forthcoming paper (Faria et al. 2015), the stars from this sample with more than 75 measurements, including HD 175607, are discussed. Neptune-mass

planets with periods lower than 50 days can be ruled out for these stars.

In the literature, there are only few examples of Neptunes or super-Earths orbiting such metal-poor stars. The planetary system around GJ 667C is one of them. It contains several super-Earths, while the star has a measured metallicity of -0.55 dex (Delfosse et al. 2013; Robertson & Mahadevan 2014). This star is an M-dwarf however and thus much cooler than HD 175607. Another Neptune system is claimed around Kapteyn's star (Anglada-Escudé et al. 2014; Bonfils et al. 2013; Robertson et al. 2015b), a very metal-poor ([Fe/H] = -0.86) halo star. This star is also an M-dwarf.

In this sense, HD 175607 would be the most metal-poor FGK dwarf to date with an orbiting Neptune. Giant planets are also rare around metal-poor stars and it has been proposed that a lower metallicity limit (~-0.7) could exist for the formation of giant planets (Mortier et al. 2012). Could the same be true for Neptunes or are we just still limited in the detection of lower-mass planets? This discovery may thus have important consequences for planet formation and evolution theories.

**Acknowledgements.** This work made use of the Simbad Database. The research leading to these results received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement number 313014 (ETAEARTH). JPF acknowledges support from FCT through grant reference SFRH/BD/93848/2013. This work was supported by Fundação para a Ciência e a Tecnologia (FCT) through the research grant UID/FIS/04434/2013. P.F., N.C.S., and S.G.S. also acknowledge the support from FCT through Investigador FCT contracts of reference IF/01037/2013, IF/00169/2012, and IF/00028/2014, respectively, and POPH/FSE (EC) by FEDER funding through the programme "Programa Operacional de Factores de Competitividade – COMPETE". This work results within the collaboration of the COST Action TD 1308. 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

Adibekyan, V. Z., Delgado Mena, E., Sousa, S. G., et al. 2012a, A&A, 547, A36

Adibekyan, V. Z., Santos, N. C., Sousa, S. G., et al. 2012b, A&A, 543, A89

Adibekyan, V. Z., Sousa, S. G., Santos, N. C., et al. 2012c, A&A, 545, A32

Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147

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

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

Bolmont, E., Raymond, S. N., Leconte, J., & Matt, S. P. 2012, A&A, 544, A124

Bonfils, X., Delfosse, X., Udry, S., et al. 2013, A&A, 549, A109

Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

Brown, D. J. A., Collier Cameron, A., Hall, C., Hebb, L., & Smalley, B. 2011, MNRAS, 415, 605

Buchhave, L. A., & Latham, D. W. 2015, ApJ, 808, 187

Burnham, K. P., & Anderson, D. R. 2002, Model selection and multimodel inference: a practical information-theoretic approach (Springer Science & Business Media)

Casagrande, L., Schönrich, R., Asplund, M., et al. 2011, A&A, 530, A138

Cunha, D., Figueira, P., Santos, N. C., Lovis, C., & Boué, G. 2013, A&A, 550, A75

da Silva, L., Girardi, L., Pasquini, L., et al. 2006, A&A, 458, 609

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

Dobbs-Dixon, I., Lin, D. N. C., & Mardling, R. A. 2004, ApJ, 610, 464

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

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

Faria, J. P., Santos, N. C., Figueira, P., et al. 2015

Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2013, ArXiv e-prints [arXiv:1306.2144]

Figueira, P., Santos, N. C., Pepe, F., Lovis, C., & Nardetto, N. 2013, A&A, 557, A93

Fischer, D. A., & Valenti, J. 2005, ApJ, 622, 1102
Ford, E. B. 2006, ApJ, 642, 505

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gandolfi, D., Hébrard, G., Alonso, R., et al. 2010, A&A, 524, A55

Gibson, N. P., Aigrain, S., Roberts, S., et al. 2012, MNRAS, 419, 2683

Glebocki, R., & Gnacinski, P. 2005, VizieR Online Data Catalog: III/244

Goggans, P. M., & Chi, Y. 2004, in Bayesian Inference and Maximum Entropy Methods in Science and Engineering, eds. G. J. Erickson, & Y. Zhai, AIP Conf. Ser., 707, 59

Gonzalez, G. 1997, MNRAS, 285, 403

Goodman, J., & Weare, J. 2010, Comm. App. Math. Comp. Sci., 5, 65

Gregory, P. C. 2005, ApJ, 631, 1198

Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517

Hébrard, G., Bouchy, F., Pont, F., et al. 2008, A&A, 488, 763

Hou, F., Goodman, J., Hogg, D. W., Weare, J., & Schwab, C. 2012, ApJ, 745, 198

Jenkins, J. S., & Tuomi, M. 2014, ApJ, 794, 110

Jenkins, J. S., Tuomi, M., Brasser, R., Ivanyuk, O., & Murgas, F. 2013, ApJ, 771, 41

Johnson, J. A., Aller, K. M., Howard, A. W., & Crepp, J. R. 2010, PASP, 122, 905

Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773

Kurucz, R. 1993, ATLAS9 Stellar Atmosphere Programs and 2 km s⁻¹ grid. Kurucz CD-ROM No. 13 (Cambridge, Mass.: Smithsonian Astrophysical Observatory), 13

Lanza, A. F., Damiani, C., & Gandolfi, D. 2011, A&A, 529, A50

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

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

Mordasini, C., Alibert, Y., Benz, W., & Naef, D. 2009, A&A, 501, 1161

Mordasini, C., Alibert, Y., Benz, W., Klahr, H., & Henning, T. 2012, A&A, 541, A97

Mortier, A., Santos, N. C., Sozzetti, A., et al. 2012, A&A, 543, A45

Mortier, A., Santos, N. C., Sousa, S., et al. 2013a, A&A, 551, A112

Mortier, A., Santos, N. C., Sousa, S. G., et al. 2013b, A&A, 558, A106

Mortier, A., Sousa, S. G., Adibekyan, V. Z., Brandão, I. M., & Santos, N. C. 2014, A&A, 572, A95

Mortier, A., Faria, J. P., Correia, C. M., Santerne, A., & Santos, N. C. 2015, A&A, 573, A101

Nardetto, N., Mourard, D., Kervella, P., et al. 2006, A&A, 453, 309

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

Pepe, F., Mayor, M., Rupprecht, G., et al. 2002, The Messenger, 110, 9

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62

Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, MNRAS, 452, 2269

Robertson, P., & Mahadevan, S. 2014, ApJ, 793, L24

Robertson, P., Endl, M., Henry, G. W., et al. 2015a, ApJ, 801, 79

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

Santerne, A., Díaz, R. F., Almenara, J.-M., et al. 2015, MNRAS, 451, 2337

Santos, N. C., Israeliian, G., & Mayor, M. 2001, A&A, 373, 1019

Santos, N. C., Mayor, M., Naef, D., et al. 2002, A&A, 392, 215

Santos, N. C., Bouchy, F., Mayor, M., et al. 2004, A&A, 426, L19

Santos, N. C., Sousa, S. G., Mortier, A., et al. 2013, A&A, 556, A150

Santos, N. C., Mortier, A., Faria, J. P., et al. 2014, A&A, 566, A35

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

Sneden, C. A. 1973, Ph.D. Thesis, The University of Texas

Sousa, S. G., Santos, N. C., Israeliian, G., et al. 2011a, A&A, 526, A99

Sousa, S. G., Santos, N. C., Israeliian, G., Mayor, M., & Udry, S. 2011b, A&A, 533, A141

Torres, G., Andersen, J., & Giménez, A. 2010, A&ARv, 18, 67

Torres, G., Fischer, D. A., Sozzetti, A., et al. 2012, ApJ, 757, 161

Udry, S., & Santos, N. C. 2007, ARA&A, 45, 397

van Leeuwen, F. 2007, A&A, 474, 653

Wang, J., & Fischer, D. A. 2015, AJ, 149, 14

Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577

## Appendix A: Model evidence from thermodynamic integration

In this section, we describe how an estimate of the model evidence can be determined using thermodynamic integration. We define the temperature-evidence function $E(\beta)$ as

$$E(\beta) = \int L^{\beta}(x)pr(x)dx, \qquad (A.1)$$

where $L(x)$ is the likelihood, $pr(x)$ the prior, and $\beta = 1/T$ with $T$ the temperature. The model evidence that we want to compute is equal to $E(1)$. Also, $E(0)$ is equal to the integrated prior. Since normalised priors are used in this work, this integrated prior $E(0)$ is equal to 1.

By using the formula for the differentiation of a natural logarithm,

$$\frac{\mathrm{d} \ln E}{\mathrm{d} \beta} = \frac{1}{E(\beta)} \frac{\mathrm{d} E(\beta)}{\mathrm{d} \beta}, \qquad (A.2)$$

and plugging in Eq. (A.1), we can write

$$\frac{\mathrm{d} \ln E}{\mathrm{d} \beta} = \frac{1}{E(\beta)} \int \ln L(x)L^{\beta}(x)pr(x)dx. \qquad (A.3)$$

The right-hand side of this equation is the average of the natural logarithm of the likelihood over the posterior at temperature $T = 1/\beta$. This is expressed as $\langle \ln L \rangle_{\beta}$.

$$\mathrm{d} \ln E = \langle \ln L \rangle_{\beta} \mathrm{d} \beta. \qquad (A.4)$$

If we now integrate both sides of this equation over the interval $[0, 1]$, we get

$$\ln E(1) = \int_{0}^{1} \mathrm{d} \ln E = \int_{0}^{1} \langle \ln L \rangle_{\beta} \mathrm{d} \beta. \qquad (A.5)$$

This integral can be estimated from the parallel-tempering ensemble sampler, embedded in *emcee*. For each temperature, the average logarithm of the likelihood is estimated from the chains. The integral can then be estimated using these values and applying a quadrature formula. From the estimation of the integral, we finally estimate the model evidence $E(1)$.