© ESO 2019

1 Introduction

The radial velocity (RV) planet search programs with the HARPS spectrograph on the ESO 3.6-m telescope (Pepe et al. 2000; Mayor et al. 2003) have contributed in a tremendous way to our knowledge of the population of small-mass planets around solar-type stars. The HARPS planet-search program on guaranteed time observations (GTO, PI: M. Mayor) was on-going for 6 yr between autumn 2003 and spring 2009. The high-precision part of this HARPS GTO survey aimed at the detection of very low-mass planets in a sample of quiet solar-type stars already screened for giant planets at a lower precision with the CORALIE Echelle spectrograph mounted on the 1.2-m Swiss telescope on the same site (Udry et al. 2000). The GTO was then continued within the ESO Large Programs 183.C-0972, 183.C-1005 and 192.C-0852 (PI: S. Udry), from 2009 to 2016.

Within these programs, HARPS has allowed for the detection (or has contributed to the detection) of more than 100 extra-solar planet candidates (see detections in Díaz et al. 2016; Moutou et al. 2015, 2011; Lo Curto et al. 2013; Dumusque et al. 2011a; Mayor et al. 2011; Pepe et al. 2011; Lovis et al. 2011a). In particular, HARPS has unveiled the existence of a large population of low-mass planets including super-Earths and hot Neptunes previous to the launch of the Kepler satellite which provided us with an overwhelming sample of thousands of small-size transiting candidates (Coughlin et al. 2016; Mullally et al. 2015; Borucki et al. 2011; Batalha et al. 2011). A preliminary analysis of the HARPS data showed at that time that at least 30% of solar-type stars were hosting low-mass planets on short-period orbits with less than 50 days (Lovis et al. 2009). A comprehensive analysis of our high-precision sample combined with 18 yr of data from CORALIE allowed us later to precise this occurrence rate: about 50% of the stars surveyed have planets with masses below 50 M_⊕ on short to moderate period orbits (Mayor et al. 2011). Furthermore, a large fraction of those planets are in multi-planetary systems. This preliminary statistics of hot super-Earth and Neptune frequency is now beautifully confirmed by the impressive results of the Kepler mission.

With the RV technique, the variation of the velocity of the central star due to the perturbing effect of small-mass planets becomes very small, of the order or even smaller than the uncertainties of individual measurements. The problems to solve and characterize the individual systems are then multi-fold, requiring to disentangle the planetary from the stellar, instrumental, and statistical noise effects. Efficient statistical techniques, mainly based on a Bayesian approach, have been developed to optimise the process and thus the outcome of ongoing RV surveys (see e.g., Dumusque et al. 2017; Díaz et al. 2016, and references therein). A large number of observations is however paramount for a complete probe of the planetary content of the system, and to take full advantage of the developed technique of analysis. In this context, focusing on the most observed, closest and brightest stars, an on-going HARPS LP (198.C-0836) is continuing the original observing efforts and providing us with an unprecedented sample of well observed stars. In this paper we describe 8 planetary systems hosting 20 planets. The detection of these planets has been announced in Mayor et al. (2011)¹ studying statistical properties of the systems discovered with HARPS. Most of them are super-Earths and Neptunian planets on relatively short periods, and member of a multi-planet system. We present here the orbital solutions for one 4-planet system around HD 20781, three 3-planet systems around HD 20003, HD 31527 and HD 136352, and four 2-planet systems around HD 21693, HD 45184, HD 51608, and HD 134060. In addition, we give updated parameters for the very eccentric planet orbiting HD 20782 (Jones et al. 2006), derived by combining HARPS, CORALIE and the published UCLES data. The paper is organised as follows. In Sect. 2, we discuss the primary star properties. Radial-velocity measurements and orbital solutions of each system are presented in Sects. 3 and 5, while Sect. 4 describes the framework used for analysing the data. We provide concluding remarks in Sect. 6.

2 Stellar characteristics

This section provides basic information about the stars hosting the planets presented in this paper. Effective temperature, gravity and metallicity are derived from the spectroscopic analysis of HARPS spectra and provided in Sousa et al. (2008). We used the improved HIPPARCOS astrometric parallaxes re-derived by van Leeuwen (2007) to determine the absolute V -band magnitude using the apparent visual magnitude from HIPPARCOS (ESA 1997). Metallicities, together with the effective temperatures and M_V are then used to estimate basic stellar parameters (ages, masses) using theoretical isochrones from the grid of Geneva stellar evolution models, including a Bayesian estimation method (Mowlavi et al. 2012).

Individual spectra were also used to derive measurements of the chromospheric activity S-index, log (R${}_{HK}^{\prime }$) and Hα, following a similar approach as used by Santos et al. (2000) and Gomes da Silva et al. (2011). Using Mamajek & Hillenbrand (2008), we estimate rotational periods for our stars from the empirical correlation of stellar rotation and chromospheric activity index (Noyes et al. 1984). We also derived the vsin(i) from a calibration of the Full Width at Half Maximum (FWHM) of the HARPS Cross-Correlation Function (CCF) following a standard approach (e.g. Santos et al. 2002). All those extra indicators are used to disentangle small-amplitude planetary signals from stellar “noise” (see e.g. the cases of CoRoT-7 and α Centauri B; Queloz et al. 2009; Dumusque et al. 2012). The derived stellar parameters are summarized in Table 1.

3 HARPS RV measurements

Radial velocities presented here have been obtained with the HARPS high-resolution spectrograph installed on the 3.6m ESO telescope at La Silla Observatory (Mayor et al. 2003). The long-term m s⁻¹ RV precision is ensured by nightly ThAr calibrations (Lovis et al. 2008). On the short timescale of a night, the high precision is obtained using simultaneous ThAr (from 2003 to 2013) or Fabry-Perot étalon (since 2013) reference calibrations. The data reduction was performed with the latest version of the HARPS pipeline. In addition to the barycentric radial velocities with internal error bars, the reduction provides the Bisector Inverse Slope (BIS) of the HARPS CCF (Pepe et al. 2002; Baranne et al. 1996), as defined by Queloz et al. (2000), the FWHM and the Contrast of the CCF, as well as the Ca II activity indexes S and log (R${}_{HK}^{\prime }$).

By construction of the HARPS high-precision sample from low-activity stars in the CORALIE volume-limited planet-search sample, the eight stars discussed here present low activity levels. The average values of the log (R${}_{HK}^{\prime }$) activity index, estimated from the chromospheric re-emission in the Ca II H and K lines at λ = 3933.66 and 3968.47 Å, are low, ranging from −5.01 to −4.87. Activity-induced RV jitter on stellar rotation time scales, due to spots and faculae on the stellar surface, is thus expected to remain at a low level. The potential influence of stellar activity on RVs is nevertheless scrutinised closely when long-term variation of activity indexes are observed (magnetic cycle) or when the planetary and stellar rotation periods (or harmonics) are of similar values.

Low-mass planets are very often found in multi-planet systems (e.g. Lovis et al. 2006, 2011a; Udry & Santos 2007; Mayor et al. 2009; Vogt et al. 2010; Díaz et al. 2016; Motalebi et al. 2015; Latham et al. 2011; Lissauer et al. 2011, 2014; Fabrycky et al. 2014, for examples of RV results and Kepler findings). The several components in the system give rise to often complex, low-amplitude RV signals, not easy to solve for. The optimal observing strategy for a given star is a priori unknown, the relevant planetary periods possibly extending over three orders of magnitude. An adequate and efficient observing strategy has been developed coping with the need to accumulate a large number of measurements and to probe different timescales of variation. We typically follow our targets every night during a few initial observing runs, and after gathering a couple of dozens of observations we compare the observed dispersion with the expected jitter for the considered stellar spectral type. For significant variations we then continue monitoring at the same cadence when high-frequency variations are seen, or we adapt the frequency to possible variations at longer timescales.

The majority of the stars in our HARPS high-precision program have been followed since 2003, gathering observations spanning more than 11 yr. We are not considering here observations obtained after May 2015 when a major upgrade of the instrument was implemented (change of the optical fibres)² . We also removed measurements obtained on the nights JD = 2455115 and 2455399 as an unexplained instrumental systematic produced RV measurements that were off by more than 10 sigma on several stars observed those nights. Each RV measurement corresponds normally to a 15-min HARPS exposure. This long exposure time allows mitigating the short-timescale variations induced by stellar oscillations and therefore to improve RV precision (Santos et al. 2004; Dumusque et al. 2011b). For bright stars (V ≤ 6.5), the exposures are split in several sub-exposures within the 15 min in order to avoid CCD saturation while keeping the same strategy to mitigate stellar oscillations. The different time series presented in this paper are composed of binned points calculated through a weighted average of all points taken within an hour, so that all the observations taken within 15 min are binned together. The stars presented here have between 178 and 245 observations typically spread over ~4000 days and with a sampling allowing the detection of planets with periods from below 1 day to the full span of the measurements. The obtained signal-to-noise ratio (S/N) at 550 nm typically ranges from 100 to 250 (up to 400 for HD 134060), depending on the star andweather conditions. The corresponding quantified uncertainties on the RVs range then from 0.33 to 0.76 m s⁻¹, including photon noise, calibration errors and instrumental drift uncertainty (see Table 2). This does not include other instrumental systematics like telescope centering and guiding errors, which are expected to be small but difficult to estimate. Possible additional uncertainties not included in this estimate might originate from the RV intrinsic variability of the star (jitter) due to stellar oscillations, granulation (Dravins 1982; Dumusque et al. 2011b) and magnetic activity (Dumusque et al. 2014a, 2011c; Meunier et al. 2010; Desort et al. 2007; Saar & Donahue 1997). Those extra errors induced by instrumental systematics or stellar signals will be taken into account when modeling the RVs (see Sect. 4).

The analysis of long series of HARPS radial velocities have recently revealed that, at the level of precision of HARPS, stitching of the detector (i.e. difference between the inter-pixel size and the gap between imprinted patches of pixels, see Sect. 4.3) may introduce a residual signal in the radial velocities at periods close to 6 months or 1 yr when a strong stellar line crosses the gaps between pixel patches due to the yearly motion of the Earth around the Sun (Dumusque et al. 2015). To avoid this effect, new sets of radial velocities for the stars presented in this paper have been obtained by removing the corresponding zone of the spectra in the template used for the cross correlation.

The final one-hour binned HARPS rsadial velocities and log (R${}_{HK}^{\prime }$) measurements are displayed in the left column of Fig. 1 for the 8 planet-host stars discussed in the paper. These velocities and the parameters inferred from the spectra and CCFs are provided in electronic form at CDS. A sample of these data is provided in Table 3. The statistics of the RV series are listed in Table 2.

4 General approach of the data analysis

The data analysis presented in this paper is performed using a set of online tools available from the DACE platform³ . The RV toolon this platform allows to upload any RV measurements and then perform a multiple Keplerian adjustment to the data using the approach described in Delisle et al. (2016). Once preliminary Keplerian plus drift parameters are found using this iterative approach, it is possible to run a full Bayesian Markov chain Monte Carlo (MCMC) analysis using an efficient algorithm, allowing various choices of models for signals (and noise) of different origins (Díaz et al. 2016, 2014).

For each system, we performed a full MCMC analysis, probing the following set of variables for planetary signals: log P, log K, $\sqrt{e}\cos \text{ ω}$, $\sqrt{e}\sin \text{ ω}$ and λ₀ , each one corresponding to the period, the RV semi-amplitude, the eccentricity, the argument of periastron and the mean longitude at a given reference epoch. We used $\sqrt{e}\cos \text{ ω}$ and$\sqrt{e}\sin \text{ ω}$ as free parameters rather than the eccentricity and the argument of periastron because they translate into a uniform prior in eccentricity (Anderson et al. 2011). The mean longitude (λ₀ ) is also preferred as a free parameter (instead of the mean anomaly or the date of passage through periastron) since this quantity is not degenerated at low eccentricities.

When no hint of magnetic activity can be seen in the different activity observables (S-index, log (R${}_{HK}^{\prime }$), Hα-index), we fitted a RV model composed of Keplerians plus a polynomial up to the second order. The MCMC analysis is performed with uniform priors for all variables, with the exception of the stellar mass for which a Gaussian prior is used based on the information given in Table 1. We chose an error in stellar mass of 0.1 M_⊙ for all systems to propagate the stellar mass error to the estimation of the planet masses. To take into account uncertainty due to instrumental systematics and/or stellar signals not estimated by the reduction pipeline, we included in the MCMC analysis a white-noise jitter parameter, σ_JIT, that is quadratically added to the individual RV error bars.

When a magnetic cycle is detected in the different activity observables (S-index, log (R${}_{HK}^{\prime }$), Hα-index), we decided, in addition to the model described above, to include two extra components in our RV model to account for the RV variation induced by this magnetic cycle. As explained in Meunier et al. (2017) and Dumusque et al. (2011a), the variation of the total filling factor of spots and faculae along a magnetic cycle change the total amount of stellar convective blueshift, as it is reduced inside spots and faculae due to strong magnetic fields, which thus change the absolute RV of the star. A positive correlation between the different activity observables and the RV is thus expected. In this paper, we decided to use the method proposed by Meunier & Lagrange (2013) to mitigate the impact of magnetic cycles, i.e. to adjust a linear correlation between RV and one of the activity index. Here we chose to use the log (R${}_{HK}^{\prime }$) (Vaughan et al. 1978; Wilson 1968; Noyes et al. 1984). A larger number of spots and faculae are present on the stellar surface during high-activity phases of the magnetic cycle, which implies a stronger stellar jitter due to those surface features coming in and out of view. To account for this stellar jitter that changes in amplitude along the magnetic cycle, we included in the MCMC analysis two extra white-noise jitter parameters, σ_{JIT LOW} and σ_{JIT HIGH}, that correspond to the jitter during the minimum and the maximum phases of the magnetic cycle. We therefore replace σ_JIT described in the precedent paragraph by σ_{JIT LOW} + (σ_{JIT HIGH} − σ_{JIT LOW}).norm( log (R${}_{HK}^{\prime }$)), where norm(log (R${}_{HK}^{\prime }$)) corresponds to log (R${}_{HK}^{\prime }$) normalized from 0 to 1. For each RV measurement, a new jitter parameter is derived according to the activity level and is quadratically added to the corresponding RV error bar. This is similar to the approach adopted in Díaz et al. (2016). The list of all the parameters probed by our MCMC is given in Table 4.

Before running a full MCMC analysis to obtain reliable posteriors for the orbital parameters of the planets present in the RV measurements, we first iteratively look for significant signals in the data using the approach of Delisle et al. (2016), which gives us a good approximation of the orbital parameters that are used as initial conditions for the MCMC analysis. The iterative approach follows several independent but complementary steps.

Fig. 1 Left column, from top to bottom panels: HARPS RV measurements as a function of barycentric Julian date obtained for HD 20003, HD 20781, HD 21693, HD 31527, HD 45184, HD 51608, HD 134060 and HD 136352. Error bars only include photon and calibration noise. Middle column: log (R${}_{HK}^{\prime }$) activity indicator as a function of time for the same stars. Right column: GLS periodogram of the corresponding RV measurements.

4.1 Removing long-term trends

Long-term trends in the data, due to either long-period companions (stellar or planetary) or magnetic cycles (Dumusque et al. 2011a; Lovis et al. 2011b that provides a broad view of magnetic cycle behavior for the entire HARPS high-precision RV survey), might perturb the detection of planets on shorter period orbits due to the spread of the energy of the signal at high frequencies in the Generalized Lomb-Scargle periodogram (GLS, Zechmeister & Kürster 2009; Scargle 1982) and aliases of these signals appearing at shorter periods. We fit a polynomial up to the second order to account for a non-resolved long-period companion. In addition, if a significant long-period signal is seen in the GLS periodogram of the calcium activity index log (R${}_{HK}^{\prime }$), we fit the RVs with a polynomial up to the second order, plus a Keplerian model reproducing the log (R${}_{HK}^{\prime }$) variation, leaving only the amplitude as a free parameter. Note that when correcting magnetic cycle effect in RVs when adopting a step-by-step approach, we do not consider a simple linear correlation between RV and log (R${}_{HK}^{\prime }$) as explained above because significant planetary signals not yet removed from the data might destroy any existing correlation. There is no a priori reason why the magnetic cycle should look like a Keplerian, however a Keplerien has more degrees of freedom than a simple sinusoid and can therefore better estimate the long-term variation seen in log (R${}_{HK}^{\prime }$). Regarding signal significance, a signal is considered worth looking into when its p-value, which gives the probability that the signal appears just by chance, is smaller or equal to a chosen threshold. Experience has shown us that a p-value of 1% provides a good empirical threshold. Note that in this paper, p-values are estimated in a Monte Carlo approach by bootstrapping 10 000 times over the dates of the observations.

4.2 Detection of periodic signals in the data

Planetary signals are searched for in the GLS periodogram of the RV residuals after correcting long-term trends. As for magnetic cycles, we consider a signal as significant if its p-value is smaller than 1%⁴ . If the p-values of some signals in the GLS periodogram are smaller than 1%, we adjust a Keplerian orbital solution to the signal presenting the smallest p-value following the approach of Delisle et al. (2016) with the detected period as a guess value.

Additional planets in the systems are then considered if other significant signals are present in the GLS periodogram of the RV residuals. Note that each time a Keplerian signal is added to the model, a global fit including all the previously detected signals is performed. We stopped when no more signals in the residuals present a p-value smaller than 1% and we finally visually inspected the solutions.

4.3 Origin of the different signals found

After the main periodic signals have been recognized in the data, it is important to identify the ones that are very likely not from planetary origin. They may be of different natures: stellar, instrumental or observational. Here is a brief description of some cases encountered.

Stellar origin

The signal is due to stellar intrinsic phenomena. The most common one is the variation of the measured RVs due to spectral line asymmetries induced by spots and faculae coming in and out of view on the surface of the star, modulating the signal over a star rotational period (Haywood et al. 2016, 2014; Dumusque et al. 2014b; Meunier et al. 2010; Desort et al. 2007; Santos et al. 2000; Saar & Donahue 1997). In order to recognize RV variations of intrinsic stellar origin, we compare the periods of the derived orbital solutions with the rotational periods of the stars and their harmonics estimated using the log (R${}_{HK}^{\prime }$) average activity level (see Table 1, Mamajek & Hillenbrand 2008; Noyes et al. 1984). We also compare the orbital periods with the variation timescales of spectroscopic activity indicators derived from the CCF, such as the BIS and the FWHM,and derived from spectral lines sensitive to activity, such as the Ca II H and K lines (S-index and log (R${}_{HK}^{\prime }$)) or the Hα line (Hα-index).

Fig. 2 GLS Periodogram of the RV residuals of HD 20003 at each step after removing, from top-left to bottom-right, the magnetic cycle effect plus a second-order polynomial drift, then the two planets at 11 and 33 days one after the other, and finally the signal at 184 days. The GLS periodogram of the raw RVs is shown in Fig. 1.

Instrumental origin

As already mentioned earlier, spurious signals of small amplitudes with periods close to one year or one of its harmonics (half a year, a thirdof a year) can be created by a discontinuity in the wavelength calibration introduced by the stitching of the detector. The 4k × 4k HARPS CCD is composed of 32 blocks of 512 × 1024 pixels. When each block were imprinted to form the detector mosaic, the technology at the time was not precise enough to ensure that pixels between blocks had the same size than intra-pixels. Therefore, every 512 pixels in the spectral direction, the CCD presents pixels that differs in size (Wilken et al. 2010). Block stitching may introduce a residual signal in the RVs at periods close to 6 months or 1 yr when strong stellar lines crosse block stitchings due to the yearly motion of the Earth around the Sun. As mentioned in Sect. 3, in order to mitigate this effect new sets of RVs for the stars presented in this paper have been obtained by removing from the correlation masks used the spectral lines potentially affected. For each star, the systemic velocity will shift the stellar spectrum on the CCD, therefore an optimization of the correlation mask is done on a star-by-star basis (for more information, see Dumusque et al. 2015).

Observational limitations

Aliases have to be taken into account. The ones due to one year or one day sampling effects are well known (e.g., Dawson & Fabrycky 2010). They apply on all signals and not only on the planetary ones.

5 Description and analysis of individual systems

In this section, we present eight planetary systems including Neptune-mass and super-Earth planets. Presence of planets around these stars and preliminary system characterizations have been announced at the “Extreme Solar System II” conference held at Grand Teton USA in September 2011 and published in Mayor et al. (2011). We present here the detailed analysis of each system with updated data, describing and discussing the derived planetary system characteristics, as well as other (i.e. non-planetary) signals in the data.

5.1 HD 20003: Eccentric, close-in Neptune-mass planets close to the 3:1 commensurability

5.1.1 Radial velocity analysis

An ensemble of 184 high S/N observations (⟨ S/N⟩ of 110 at 550 nm) of HD 20003 have been gathered covering about 11 yr (4063 days). The typical photon-noise plus calibration uncertainty of the observations is 0.74 ms⁻¹, well below the observed dispersion of the RVs at 5.35 m s⁻¹. The RVs with their GLS periodogram and the log (R${}_{HK}^{\prime }$) time-series are displayed in Fig. 1. A long-period variation is observed in the RV data, as well as in log (R${}_{HK}^{\prime }$), indicative of a magnetic cycle effect on the velocities. To correct the velocities for this effect, we modeled the long-term variation by a Keplerian with parameters fixed and determined from the log (R${}_{HK}^{\prime }$), except for the amplitude that was left free to vary.

After correction for the long-period variation due to the magnetic cycle plus fitting a second-order polynomial drift, two peaks very clearly emerge in the periodogram well above the 0.1% p-value limit, at periods around 11.9 and 33.9 days (Fig. 2). Once those two signals are fitted for along with the magnetic cycle and the second-order polynomial drift, a significant signal at 184 days appears in the GLS periodogram of the RV residuals, with strong aliases at 127 and 359 days. A global fit including a second-order polynomial drift, a Keplerian for the magnetic cycle and three Keplerians for signals at 11.9, 33.9 and 184 days allows us to model all the variations seen in the RVs. No significant signal (with p-values smaller than 5%) are then present in the RV residuals (Fig. 2, bottom-right plot).

After this preliminary phase looking for significant signals in the data, we searched for the best-fit parameters with an MCMC, using a model composed of a second-order polynomial drift, a linear correlation of the RVs with log (R${}_{HK}^{\prime }$) to adjust the magnetic cycle effect, three Keplerians to fit for the signals at 11.9, 33.9 and 184 days, and two jitters that correspond to the instrumental plus stellar noise at the minimum and maximum of the magnetic cycle (see Sect. 4). This model converged to a stable solution. However, a significant signal in the residuals was still present near 3000 days (Fig. 3), a period similar to the one of the magnetic cycle, an indication that our modeling of the long period variation was not satisfactory, very probably because the magnetic cycle of HD 20003 is not fully covered, as seen in Fig. 1. We finally decided then to adjust a model composed of 3 Keplerians at 11.9, 33.9 and 184 days, an extra Keplerian atlong-period to absorb the effect of the magnetic signal plus any long-term trend that could be present in the data, and one jitter corresponding to the instrumental plus stellar noise. In this case, the GLS periodogram of the residuals is extremely similar to the last plot in Fig. 2 and therefore no more significant signal is present in those RV residuals. Note that the period of the signal adjusting the magnetic cycle changed however from 4038 days (log (R${}_{HK}^{\prime }$)) to 3298 days (RVs). Although we do not understand fully the long-term signal present in the RVs of HD 20003, adjusting it with different models gives similar orbital parameters for the signals at 11.9, 33.9 and 184 days. We adopted as final solution the model with the magnetic cycle modeled by a long-period Keplerian with free parameters.

The best-fit for each inner planet, the signal at 184 days, and the RV residuals are displayed in Fig. 4. The best-fit parameters can be found in Table 5. A careful look at the activity indicators, log (R${}_{HK}^{\prime }$), BIS SPAN and FWHM, after removing the long-period signal induced by the stellar magnetic cycle, reveals no significant peaks that match the planetary signals found in our analysis (see Fig. A.1). We are therefore confident that those signals are not induced by stellar activity.

The correlation between the activity index log (R${}_{HK}^{\prime }$) and the RV residuals when removing all the detected signals except the magnetic cycle effect can be seen in the left panel of Fig. 5. Note that this dependence in activity is included in our model when adjusting the long-period Keplerian.

Fig. 3 GLS Periodogram of the RV residuals of HD 20003 obtained from the MCMC solution including a linear correlation between RV and log (R${}_{HK}^{\prime }$) to account for the magnetic cycle, a second-order polynomial drift and three Keplerian.

5.1.2 A 3rd planet with a period of 184 days?

Although the signals at 11.9 and 33.9 days are clearly induced by planets orbiting around HD 20003, the signal at 184 days is more difficult to interpret as it is very close to half a year, with strong aliases at 127 days and 359 days (close to 1/3 and 1 yr). As explained in Sect. 4.3, signals at a year or harmonics of it can be induced by a discontinuity in the wavelength calibration introduced by tiny gaps between the different quadrants of the detector. However, the RV data that we analyze here have been corrected for this effect following the work of Dumusque et al. (2015). This correction has been already applied to the RVs of several stars, and has always been successful in removing the spurious signal. The fact that the RV residuals after removing the two planets, and the long-period effect of the magnetic cycle do not correlate with the barycentric Earth velocity (BERV, see right panel of Fig. 5) disfavors the hypothesis that the 184-day signal is due to a discontinuity in the wavelength calibration introduced by tiny gaps between the different quadrants of the detector. In conclusion, we would be inclined to claim that this signal at 184 days is a real planetary signal. However, although intensively tested, there is always the possibility that the data reduction system does not correct for all the instrumental effects in this particular case. We leave the signal as a potential one, and encourage other teams to confirm it using other pipelines or other facilities than HARPS.

5.1.3 Discussion of the planetary system

A surprising aspect of this system is the high eccentricity of the inner planet, e = 0.38, while the planet at 33.9 days has an eccentricity close to 0. From a planet formation or a dynamical point of view, this isnot straightforward to explain. There might be actually several possibilities to explain the observed configuration. It is out of the scope if this paper to check all of them in details, we just mention here the basic ideas.

1) The first and simplest possibility is that the high eccentricity of the planet at 11.9 days actually hides a planet at half the orbital period (Anglada-Escudé et al. 2010), in (or close to) a 2:1 resonance configuration. We therefore tried to fit a model with an extra planet with initial period at 5.95 days, setting the eccentricities of the 5.95 and 11.9-day planets to zero or letting them free tovary. Model comparison shows that, in both cases, the more complex solution is disfavoured with a ΔBIC of 8.3 and 24.3, respectively. It seems therefore that the significant eccentricity of planet b is real.

2) The two planets are close to the 3:1 resonance. An increase of the eccentricity could have happened if, in the past, the two planets havebeen migrating trapped in the resonance. Then a specific event made the two planets go slightly out of resonance on the nowadays observed orbits. Possible scenarios for this event could be instability when the gas disappeared or the presence of an additional planet (e.g. resonant chain) that was ejected when the eccentricity increased (scattering).

3) Invoking additional bodies in the system raises the possibility of an external companion that could induce a Kozai-lidov effect on theinner planets in case of a large mutual inclination (Kozai 1962; Lidov 1961). Such a companion could actually be hinted by our difficulty of modelling the long-period variation observed in the radial velocities. It is however not very clear to understand the differential effect on the two inner planets.

4) Finally, a very appealing potential explanation might rely on the spin-orbit evolution of planets feeling the tidal effect of the central star and the gravitational perturbation of another planetary companion. Correia et al. (2012) have shown that “under some particular initial conditions, orbital and spin evolution cannot be dissociated and counter-intuitive behaviors can be observed, such as the secular increase of the eccentricity. This effect can last over long timescales and may explain the high eccentricities observed for moderate close-in planets”.

5.2 HD 20781: A packed system with 2 super-Earths and 2 Neptune-mass planets

HD 20781 is part of a visual binary system including another star, HD 20782, known to host a planet on a 595-day very eccentric orbit (Jones et al. 2006). We take the opportunity of the discovery of a compact system of small planets around HD 20781 to provide here as well an updated solution for the planet around HD 20782 using HARPS data.

5.2.1 Analysis of the HD 20781 system

HD 20781 was part of the original high-precision HARPS GTO survey and the star has been then followed for more than 11 yr (4093 days). Over this time span, we gathered a total of 226 high S/N spectra (⟨S/N⟩ of 112 at 550 nm) corresponding in the end to 216 RV measurements binned over 1 h. As reported in Table 2, the typical precision of individual measurements is 0.76 m s⁻¹ including photon noise and calibration uncertainties. The raw RV rms is significantly higher, at 3.41 m s⁻¹, pointing towards additional variations in the data of potentially stellar or planetary origin, assuming the instrumental effects are kept below the photon-noise level of the observations.

As a first approach we looked at the RV and activity index time series shown in Fig. 1. No significant long-term variation is observed in log (R${}_{HK}^{\prime }$) data and no long-term variation is visible in the GLS periodogram of the velocity time series. We conclude that there is no noticeable sign of a magnetic activity cycle for this star. The average value of log (R${}_{HK}^{\prime }$) at − 5.03 ± 0.01 is also very low with a very small dispersion, similar to the Sun at minimum activity.

Due to the small activity level and the large number of observations, the GLS periodogram of the velocity series is very clean, with significant peaks clearly coming out of the noise background. So, even if the final characterization of the planetary system parameters is performed through a Bayesian-based MCMC approach, a step-by-step analysis of the system, characterizing and then removing one planet after the other from the data, will provide an excellent illustration of the significance of the planet detection in this system. The most prominent peak in the GLS periodogram is at a period around 20 days. Deriving a Keplerian solution for this planet and then removing the corresponding signal from the raw RVs makes a clear signal at 86 days appear in the GLS periodogram of the residuals (Fig. 6). Keeping on with the same approach, we can clearly identify, sequentially, significant signals first at 13.9 and then at 5.3 days. It has to be noted here that none of the significant periods in the data are close to the stellar rotational period, 46.8 days, estimated from the activity index (Table 1).

As described in Sect. 4, the final determination of the planetary system parameters is performed using a MCMC sampler and a model composed of four Keplerians representing the planetary signals, and an extra white-noise jitter to consider potential stellar or instrumental noise not included in the RV error bars (the σ_JIT parameter, see Table 4). Phase-folded planetary solutions, as well as the RV residuals around the best solution, are displayed in Fig. 7. The best-fit parameters are reported in Table 6.

Looking at the periodograms of the log (R${}_{HK}^{\prime }$), the BIS SPAN and the FWHM in Fig. A.2, to check if any announced planet matches any signal in the activity indicators, we see that the log (R${}_{HK}^{\prime }$) times series presents signals at 115, 81 and 68 days, and the FWHM time series at 380 days. This later signal is probably due to interaction with the window functions, creating a signal near a year. The former signals are more difficult to interpret as they are not compatible with the estimated rotational period of the star, 46.8 days (Table 1). In particular, we can note that the 86-day Neptune-mass planet has a period close to the detected peaks in the log (R${}_{HK}^{\prime }$) time-series at 81 days, and furthermore is close to the 1-yr alias of the 115-d signal (days). The amplitude of these signals is however very low, at the level of ~0.01 dex, 20 times smaller than the solar magnetic cycle variation. Such signals can therefore probably not be responsible for the 2.6 m s⁻¹ periodic variation detected in the RVs at 86 days.

Our conclusion is thus that the system HD 20781 hosts two inner super-Earths with periods of 5.3 and 13.9 days and two outer Neptune-mass planets with periods of 29 and 86 days. With small eccentricities, stability is not a concern for this system of small-mass planets in a moderately compact configuration.

Fig. 4 Phase-folded RV measurements of HD 20003 with the best Keplerian solution represented as a black curve for each of the signals in the data. From top left to bottom right panels: planet b, planet c and the signal at 184 days. The residuals around the solution are displayed in the lower right panel. Contrarily to Fig. 1 error bars include now the instrumental + stellar jitter from the MCMC modeling. Corresponding planetary orbital elements are listed in Table 5.

Fig. 5 Left panel: activity index log (R${}_{HK}^{\prime }$) as a function of the RV residuals when removing all the detected signals except the magnetic cycle effect for HD 20003. The observed correlation indicates that most of the RV residual variation is due to activity-related effects. Right panel: barycentric Earth velocity projected on the line of sight (BERV) as a function of the RV residuals around the best derived solution without considering the 184-day signal. The fact that no correlation can be observed disfavours the hypothesis that the 184-day signal is due to a discontinuity in the wavelength calibration introduced by tiny gaps between the different quadrants of the detector.

Fig. 6 GLS periodogram of the residuals at each step, after removing one planet after the other in the analysis of HD 20781 (from left toright and top to bottom). The GLS periodogram of the raw RVs is shown in Fig. 1.

Fig. 7 Phase-folded RV measurements of HD 20781 with the best-fit solution represented as a black curve for each of the signal foundin the data (from left to right and top to bottom: planets b–e). Error bars include photon and calibration noise, as well as instrumental + stellar jitters derived from the MCMC modeling. The residuals around the solution are displayed in the lower panel. Corresponding orbital elements are listed in Table 6.

5.2.2 HD 20782: More planets in this visual binary system

The star HD 20782 is the brightest companion of the HD 20781-HD 20782 binary system. It is known to harbor a 595-day very eccentric Jupiter planet (Jones et al. 2006). A total of 71 high S/N spectra (⟨S/N⟩ of 181 at 550 nm) were obtained with HARPS on this target, which translates into 68 RV measurements after binning the data over 1 h. Because of the very large semi-amplitude and eccentricity of the RV signal induced by the planet, we decided to include in addition to the HARPS very precise measurements the lower precision RV data obtained with UCLES (published in Jones et al. 2006) and CORALIE.

We searched for the best-fit parameters using a MCMC sampler. The solution converges to an extremely eccentric Jupiter-mass planet with a period of 597 days. The phase-folded planetary solution, as well as the RV residuals around the best solution and their GLS periodogram, are displayed in Fig. 8. The best-fit parameters are reported in Table 6. Although the eccentricity, the amplitude and the argument of periastron are compatible within one sigma with the values presented in Jones et al. (2006), the period and argument of periastron are not compatible, even if different by less than 2%. This can be explained by the addition of the CORALIE and HARPS data that allows to sample much better the periastron passage.

The high eccentricity of HD 20782 b has commonly been explained via the Kozai-Lidov mechanism (Kozai 1962; Lidov 1961) induced by HD 20781.

5.3 HD 21693: A system of 2 Neptune-mass planets close to a 5:2 resonance

Over a time span of 11 yr (4106 days), 212 high S/N spectra (⟨ S/N⟩ of 141 at 550 nm) of HD 21693 were gathered, resulting in a total of 210 RV measurements when binning the data over 1 h. The typical photon-noise and calibration uncertainty is 0.60 m s⁻¹, which is significantly below the 4.72 m s⁻¹ observed dispersion of the RVs, pointing towards the existence of additional signals in the data. The raw RVs, their GLS periodogram and the calcium activity index of HD 21693 are shown in Fig. 1. The variation of the log (R${}_{HK}^{\prime }$) ranging from −5.02 to −4.83 highlights a significant magnetic cycle. This magnetic cycle is very similar in magnitude to the one of the Sun, but slightly shorter with a period of 10 yr. When fitting a Keplerian signal to the log (R${}_{HK}^{\prime }$), we are left with significant signals in the residuals at 740 and 33.5 days (see Fig. A.3). Those signals are also present in the FWHM and the bisector span of the CCF, although less significant. The signal at ~740 days, close to two years is probably due to the sampling of the data, and the 33.5-day signal is likely the stellar rotation period. This value is compatible with the rotation period derived from the mean activity log (R${}_{HK}^{\prime }$) level, i.e. 36 days (Table 1).

Looking at the raw RVs and their GLS periodogram in Fig. 1, it is clear that the observed magnetic cycle has an impact on the measured RV measurements. To remove the RV contribution of the magnetic cycle, we remove from the RVs a Keplerian model that has the same parameters as the Keplerian fitted to log (R${}_{HK}^{\prime }$) but with an amplitude free to vary. The GLS periodogram of the RV residuals after correcting for the magnetic cycle effect are displayed in the top panel of Fig. 9. A highly significant signal at 54 days is present in the data. When removing this signal by fitting a Keplerian, an extra signal at 23 days is seen in the residuals (middle panel of Fig. 9). No signal with p-value smaller than 1% appears in the residuals of a two-Keplerian model; we therefore stop looking for extra signals. Note however that the most significant signal left in the GLS periodogram corresponds to a period of 16 days, likely the first harmonic of the stellar rotation period, which is expected from stellar activity (Boisse et al. 2011).

After this preliminary phase looking for significant signals in the data, we search for the best-fit parameters with an MCMC, using a model composed of a linear correlation of the RVs with log (R${}_{HK}^{\prime }$) to adjust the magnetic cycle effect, two Keplerian functions for the signals at 23 and 54 days, and two jitter values correspond to the instrumental plus stellar noise at the minimum and maximum of the magnetic cycle (see Sect. 4). The best-fit for each planet and the RV residuals are displayed in Fig. 10, and the best-fit parameters can be found in Table 7. None of the two detected planet corresponds to signals found in the activity indicators (see Fig. A.3).

When looking at the RV residuals, the scatter is still rather high even after removing all the significant signal detected in the data. Although in HIPPARCOS the star is catalogued as a G8 dwarf, the spectroscopic survey of nearby stars NSTAR finds that HD 21693 is a G9IV-V therefore a slightly evolved star (Gray et al. 2006). Evolved stars present higher photometric and RV jitter associated with more significant granulation, which might explain this significant residual jitter (Bastien et al. 2014, 2013; Dumusque et al. 2011b).

The correlation between the activity index log (R${}_{HK}^{\prime }$) and the RV residuals removing only the two-planet solution can be seen in Fig. 11. Note that this correlation is considered in the MCMC model we used to fit the RV data.

Our analysis of the HARPS RV measurements of HD 21693 provides strong evidence that 2 Neptune-mass planets orbit the star, with periods of 22.7 and 53.7 days. With such periods, this planetary system is close to a 5:2 resonance but with a period ratio smaller than the 5:2 commensurability. If the system formed through a convergent migration scenario trapping the planets in the 5:2 resonance, some process has to be invoked to push the system out of the resonance configuration after the disappearance of the protoplanetary disk.

Fig. 8 Best Keplerian solution for the eccentric planet orbiting HD 20782. Top panel: phase-folded RV measurements with the best solution represented as a black curve. Bottom panel: RV residuals around the best Keplerian solution.

Fig. 9 Top panel: GLS periodogram of the residual RVs of HD 21693 after removing the RV contribution of the magnetic cycle. Middle and bottom panels: GLS periodogram of the residuals at each step, after removing one planet after the other in the analysis. The GLS periodogram of the raw RVs is shown in Fig. 1.

5.4 HD 31527: A 3-Neptune system

In total, 257 high S/N spectra (⟨ S/N⟩ of 180 at 550nm) of HD 31527 were gathered over a time span of 11 yr (4135 days). This results in 245 observations of the star when data are binned over 1 h. The 3.19 m s⁻¹ observed dispersion of the RVs is much larger than the typical RV precision of 0.64 m s⁻¹, pointing towards the existence of extra signals in the data. The raw RVs, their GLS periodogram and the calcium activity index of HD 31527 are shown in Fig. 1. The calcium activity index log (R${}_{HK}^{\prime }$) does not show any significant variation as a function of time, therefore we do not expect the RVs to be affected by long-period signals generally induced by magnetic cycles. The mean activity level of the star is equal to log (R${}_{HK}^{\prime }$) = −4.96, close to solar minimum. The RVs should therefore be exempt as well of activity signal at the rotational period of the star and its harmonics due to active regions on the stellar surface (Boisse et al. 2011). This is confirmed when looking at the periodograms of the different activity indicators in Fig. A.4. Only a signal at 400 days in the BIS SPAN is significant. This signal, not too far from a year, might be due the interplay between the time series and the window function.

As we can see in the GLS periodogram of the raw RVs (Fig. 1), two extremely significant signals appear at 17 and 52 days. After fitting these signals with a two-Keplerian model, a third significant signal at 271 days can be seen in the RV residuals (middle panel of Fig. 12). Finally, the residuals of a three-Keplerian model do not show any signal with p-value smaller than 1% (bottom panel of Fig. 12). The fact that no signal is present at the estimated rotation period of the star (19 days, Table 1) or its harmonics confirms that the RVs of HD 31527 are not affected by significant activity signal.

After this first search for significant signals in the data, we fitted, using a MCMC sampler, a three-Keplerian model to the data including a white-noise jitter component to account for stellar and instrumental uncertainties not included in the RV error bars. The best-fit for each planet, as well as the RV residuals, can be seen in Fig. 13. We report in addition the best-fit parameters in Table 8.

None of the signals announced here matches signals in the different activity indicators (see Fig. A.4).We therefore conclude that HD 31527 hosts 3 Neptune-mass planets, with periods of 16.6, 51.2 and 272 days. The star is a G2 dwarf like the Sun, therefore the outer planet in this system lies on an orbit between the ones of Venus and the Earth, therefore in the habitable zone of its host star (Selsis et al. 2007). This planet, with a minimum mass of 13 Earth-masses, is however most likely composed ofa large gas envelope (Rogers 2015; Wolfgang & Lopez 2015; Weiss & Marcy 2014), except if it is similar to Kepler-10 c in composition (Dumusque et al. 2014a).

Fig. 10 Phase-folded RV measurements of HD 21693 with the best planet solution represented as a black curve for each of the signal in the data (top to bottom: planet b and planet c). The residuals around the solution are displayed in the lower panel.Corresponding orbital elements are listed in Table 7.

Fig. 11 RV residuals when removing all the detected signals except the magnetic cycle effect plotted as a function of the activity index log (R${}_{HK}^{\prime }$) for HD 21693. The observed correlation indicates that most of the RV residual variation is due to activity-related effects.

Fig. 12 Top panel: GLS periodogram of the residual RVs of HD 31527 after removing the best-fit Keplerian to account for the significant signal at 17 days seen in the raw RVs (bottom right panel of Fig. 1). Middle and bottom panels: GLSperiodogram of the residuals at each step, after removing the second and third planet from the RVs.

5.5 HD 45184: A system of two close-in Neptunes

We gathered a total of 309 high S/N spectra (⟨ S/N⟩ of 221 at 550 nm) of the G1.5 dwarf HD 45184 during a time span of 11 yr (4160 days). This results in 178 RV measurements, when the data are binned over 1 h. The average precision of the RVs is 0.41 m s⁻¹ considering only photon noise and calibration uncertainties. The raw RV rms is much higher, 4.72 m s⁻¹, which implies that significant signals are present in the data. In Fig. 1, we display the raw RVs and their GLS periodogram, and the log (R${}_{HK}^{\prime }$) time series. We see that a significant magnetic cycle affects log (R${}_{HK}^{\prime }$), with values ranging from −5.00 to −4.86 with a periodicity of 5 yr. This magnetic cycle is therefore smaller in amplitude and period than the solar cycle, despite the similarities of the derived stellar parameters with those of the Sun. To see if significant signals were present in the calcium activity index despite the long-period magnetic cycle, we fitted the log (R${}_{HK}^{\prime }$) time series with a Keplerian. In the residuals, a strong signal at 20 days is present, likely corresponding to the stellar rotation period (see Fig. A.5). This value is fully compatible with the rotation estimated using the log (R${}_{HK}^{\prime }$) average level (21.5 days, see Table 1; Mamajek & Hillenbrand 2008; Noyes et al. 1984).

Looking at the raw RVs of HD 45184 and their GLS periodogram in Fig. 1, we see that the magnetic cycle observed in log (R${}_{HK}^{\prime }$) has an influence on the RVs. To remove the RV contribution of the magnetic cycle, we fitted the log (R${}_{HK}^{\prime }$) with a Keplerian, and removed the same Keplerian from the RVs leaving the amplitude as a free parameter. In the residuals, displayed in the top panel of Fig. 14, we see a significant signal at 6 days. Once this signal is removed by fitting a Keplerian with the corresponding period, another signal at 13 days appears (see middle panel of Fig. 14). Finally after removing a two-Keplerian model to the RVs corrected for the magnetic cycle effect, no signal with p-value smaller than 10% remains. We therefore stop looking for extra signals in the data. Note that although not significant, the highest peak in the GLS periodogram of the RV residuals is at 18.6 days, likely the imprint of the stellar rotation period. The RVs are therefore slightly affected by stellar activity, however at a level that is not perturbing the detection and characterization of the two planets at 6 and 13 days.

After the preliminary phase of looking for significant signals, we fitted the RVs using a MCMC sampler and a model composed of a linear correlation of the RVs with log (R${}_{HK}^{\prime }$) to adjust the magnetic cycle effect, two Keplerians to fit for the signals at 5.9 and 13.1 days, and two jitters that correspond to the instrumental plus stellar noise at the minimum and maximum of the magnetic cycle (see Sect. 4). Each planet with its best fit can be seen in Fig. 15, as well as the RV residuals after the best solution has been removed. The best-fit parameters are reported in Table 9. None of the signals at 5.9 and 13.1 matched signals in the different activity indicators (see Fig. A.5), and therefore those signal are associated with bona-fide planets.

In Fig. 16, we show the RV residuals after removing the best-fit solution for planets b and c as a function of the log (R${}_{HK}^{\prime }$). The observed strong correlation indicates that most of the residuals are due to activity-related effects and motivates the use of our model that includes a correlation between log (R${}_{HK}^{\prime }$) and RVs to mitigate the effect of long-term activity.

Fig. 13 Phase-folded RV measurements of HD 31527 with the best Keplerian solution for planets b–e represented as black curves (fromleft to right and top to bottom). Error bars include photon and calibration noise, as well as a jitter estimated from the MCMC analysis. The residuals around the solution are displayed in the bottom right panel. Corresponding orbital elements are listed in Table 8.

5.6 HD 51608: 2 Neptune-mass planets

Over a time span of 11 yr (4158 days), 218 high S/N spectra (⟨ S/N⟩ of 133 at 550 nm) of HD 51608 were gathered with HARPS, resulting in a total of 216 measurements binned over 1 h with a typical photon-noise and calibration uncertainty of 0.62 m s⁻¹. This value is significantly below the observed 4.07 m s⁻¹ dispersion of the RVs, pointing towards the existence of additional signals in the data. The raw RVs, their GLS periodogram and the calcium activity index of HD 51608 are shown in Fig. 1. A small, albeit significant, long-term variation can be seen in log (R${}_{HK}^{\prime }$), with values in log (R${}_{HK}^{\prime }$) ranging from −5.04 to −4.96 and with a period of 11 yr. Although the period of this magnetic cycle is very similar to that of the Sun, its amplitude is much lower. After fitting this long-period signal in log (R${}_{HK}^{\prime }$), a signal with a p-value of ~5% and a period of 37 days is detected in the log (R${}_{HK}^{\prime }$) residuals, and can alsobe seen in the BIS SPAN and the FWHM of the CCF (see Fig. A.6). This is likely a signature of stellar activity as the mean log (R${}_{HK}^{\prime }$) level gives an estimated rotation period of 40 ± 4 days (Table 1; Mamajek & Hillenbrand 2008; Noyes et al. 1984).

In the raw RVs, very strong signals at 14 and 96 days are present (see Fig. 1). Once fitting a two Keplerian model to account for those signals, a long-period signal with a p-value smaller than 0.1% appear in the GLS periodogram (see middle panel of Fig. 17). This signal is induced by the stellar magnetic cycle, and we remove it as in the precedent cases by fitting a Keplerian to the log (R${}_{HK}^{\prime }$), and removing the same Keplerian from the RVs leaving the amplitude as a free parameter. After fitting the effect of the two planets plus the magnetic cycle, no significant signal with p-value smaller than 5% is left in the residuals.

After this preliminary stage of checking significant signals in the data, we searched for the best-fit parameters using a MCMC sampler, and selecting a model composed of a linear correlation between the RVs and log (R${}_{HK}^{\prime }$) two Keplerians to fit for the planetary signals at 14.1 and 96.0 days, and two jitters that correspond to the instrumental plus stellar noise at the minimum and maximum of the magnetic cycle (see Sect. 4). The best-fit solution for the two planets are shown in Table 10 and illustrated in Fig. 18, along with the RV residuals. The two signals detected in RVs are not matching any significant signal in the different activity indicators (see Fig. A.6) and are therefore associated to bona-fide planets.

The correlation between the activity index log (R${}_{HK}^{\prime }$) and the RV residuals removing only the two-planet solution can be clearly seen in Fig. 19. Note that this correlation is considered in the model we used to fit the RV data.

With masses of 14.3 and 12.8 Earth-masses, the two planets orbiting HD 51608 are likely Neptune-type planets. There is however still the possibility to have (one of) them similar in composition to Kepler-10 c (Dumusque et al. 2014a).

Fig. 14 From top to bottom panels: GLS periodogram of the RVs after removing the effect induced by the stellar magnetic cycles and then one planet after the other in the analysis of HD 45184. The GLS periodogram of the raw RVs is shown in Fig. 1.

Fig. 15 Phase-folded RV measurements of HD 45184 with, from top to bottom, the best Keplerian solution for planets b and c, and the corresponding residuals. Error bars include photon and calibration noise, as well as a jitter effect (stellar + instrumental) determined in the MCMC analysis. Corresponding orbital elements are listed in Table 9.

Fig. 16 RV residuals when removing all the detected signals except the magnetic cycle effect plotted as a function of the activity index log (R${}_{HK}^{\prime }$) for HD 45184. The observed correlation indicates that most of the RV residual variation is due to activity-related effects.

Fig. 17 From top to bottom panels: GLS periodogram of the RV residuals of HD 51608 after removing the effect of the magnetic cycle and, at each step, after removing one planet after the other. The GLS periodogram of the raw RVs is shown in Fig. 1.

5.7 HD 134060: A short-period Neptune on an eccentric orbit with a long-period more massive companion

5.7.1 Data analysis

A total of 335 high S/N spectra (⟨ S/N⟩ of 199 at 550 nm) of HD 134060 have been gathered over a time span of 11 yr (4083 days). When binning the measurements over one hour, we are left with 155 RV measurements, with a typical photon-noise plus calibration uncertainty of 0.40 m s⁻¹. This is an order of magnitude below the observed dispersion of the RVs, 3.68 m s⁻¹. The raw RV with the corresponding GLS periodogram and the log (R${}_{HK}^{\prime }$) time series are displayed in Fig. 1. The log (R${}_{HK}^{\prime }$) time series presents only a very tiny long-term trend, with very low-level variations between −5.05 and −5. We therefore do not expect strong signals induced by stellar activity. This is confirmed by the fact that no significant signals appears in the periodogram of the different activity indicators in Fig. A.7.

One very significant signal at 3.3 days can be seen in the GLS periodogram of the raw RVs (see Fig. 1). Once this signal is fitted with a Keplerian, another significant peak appears at 1292 days, as can be seen in Fig. 20. After fitting simultaneously those two signals, nothing is left in the RV residuals with p-values smaller than 10 %.

To get the best possible orbital parameters for those two planets with reliable error bars, we perform an MCMC analysis with a model composed of two Keplerians plus a white-noise jitter to account for stellar and instrumental uncertainties not included in the RV error bars. The best-fit parameters can be found in Table 11. The planetary signals, folded in phase, can be seen in Fig. 21 along with the RV residuals shown in the bottom plot.

Fig. 18 Phase-folded RV measurements of HD 51608 with, from top to bottom the best Keplerian solution for planets b, c and the residuals around the best-fitted solution. Corresponding orbital elements are listed in Table 10.

5.7.2 Origin of the high eccentricity of the inner planet

The MCMC converges to a solution with an inner planet of minimum mass 10.1 M_⊕ and period 3.27 days, on a relatively high eccentricity orbit with e = 0.45. An inner eccentric planets in a system with an additional outer planet is reminiscent of the system HD 20003 above. The main differences are that in HD 134060 the outer companion is farther out and more massive, and the inner planet is experiencing stronger tides from the star making the circularization timescale shorter. Nevertheless some of the potential mechanisms proposed for explaining the high eccentricity might be considered as well.

1) The high eccentricity can hide a planetary system in 2:1 resonance, thus the existence of another planet at half its orbital period (Anglada-Escudé et al. 2010). We therefore tried to fit a model with an extra planet at 1.65 days, fixing the eccentricities at zero or leaving them free to vary. In both cases, the more complex solution is disfavoured by Bayesian model comparison with a ΔBIC of 3.1 and 23.5, respectively. We therefore keep the simplest solution with the relatively high eccentricity of the inner planet.

2) This planet has a long-period companion that has a minimum mass three times larger. In the case of the two planetary orbits being strongly non-coplanar it is likely that the long-period planet is perturbing its inner companion through a Kozai-Lidovmechanism causing libration of its orbit (Kozai 1962; Lidov 1961). During this process, the eccentricity of the inner planet can reach very high values. The inner planet starts then to interact with its host star during close approaches, inducing a circularization of the inner orbit. Because of the conservation of the total angular momentum, eccentricity can increase only if the inclination of the orbit changes. Inner planets under the influence of a Kozai-Lidovmechanism are therefore likely to be on inclined orbits relative to the stellar rotational plane. This can be measured using the Rossiter-McLaughlin effect if by chance the planet transits its host star. This has been checked for HD 134060 by Gillon et al. (2017) using the Spitzer space telescope with unfortunately a null result. Another difficulty here is that for a planet at 3.3 days, the circularization timescale is normally very short for typical values of the tidal dissipation parameters, therefore a priori preventing us from observing the system in such a configuration.

3) Scattering is always a valid possibility when high eccentricities are involved. As for the previous point, preventing the system to then circularize is also the challenge.

4) The explanation invoking spin-orbit coupling during the evolution of the planet orbit under the influence of tides from the central star and gravitational perturbation from the outer companion (Correia et al. 2012) might finally be the most appealing process. It has however still to be demonstrated that the process works for the present set of parameters. As often, the challenge is in the balance of timescales. To help, HD 134060 seems to be a relatively young sub-giant with an age estimate of 1.75 Gyr (Delgado Mena et al. 2015).

In any case, with the high eccentricity of its inner planet HD 134060 provides a very interesting system to study planetary formation and evolution.

Fig. 19 RV residuals around the best derived solution without considering the magnetic cycle effect plotted as a function of the activity index log (R${}_{HK}^{\prime }$) for HD 51608. The observed correlation indicates that most of the RV residual variation is due to activity-related effects.

Fig. 20 From top to bottom panels: GLS periodogram of the residuals at each step, after removing one planet after the other in the analysis of HD 134060. The GLS periodogram of the raw RVs is shown in Fig. 1.

5.8 HD 136352: A 3-planet system

HD 136352 was part of the original high-precision HARPS GTO survey and the star has been followed for nearly 11 yr (3993 days). Over this time span, we gathered a total of 649 high S/N spectra (⟨ S/N⟩ of 231 at 550 nm) corresponding in the end to 240 RV measurements binned over 1 h. As reported in Table 2, the typical precision of individual measurements is 0.33 m s⁻¹ including photon noise and calibration uncertainties, an order of magnitude smaller than the observed raw RV rms of 2.74 m s⁻¹. This suggests that significant signals, of stellar or planetary origin, are present in the data. As a first approach we looked at the log (R${}_{HK}^{\prime }$) activity index time series in Fig. 1. No significant long-term variation is observed in the log (R${}_{HK}^{\prime }$) data and no long-term variation is visible in the GLS periodogram of the velocity time series. We conclude that there is no noticeable sign of a magnetic activity cycle for this star. The average value of log (R${}_{HK}^{\prime }$) at − 4.95 is low with a small dispersion of ~0.01, close to the Sun at minimum activity. No significant effect of stellar activity in the RV measurements is thus expected for this star. This is confirmed by the fact that no significant signals appears in the periodogram of the different activity indicators in Fig. A.8.

Due to the small activity level and the large number of observations, the GLS periodogram of the velocity series is actually very clean, with peaks at 27.6, 11.6 and 108 days, in order of decreasing significance (see Fig. 22). After fitting those three signals with Keplerians, a study of the GLS periodogram of the RV residuals shows a peaks at 123 days with a p-value between 1 and 0.1 %, thus an interesting signal that we will consider in the MCMC analysis.

The best orbital parameters for the three planets orbiting HD 136352 are searched for using a MCMC sampler using a model composed of three Keplerians and an extra white-noise jitter to account for instrumental and stellar uncertainties not included in the RV error bars. Phase-folded planetary solutions are displayed in Fig. 23, as well as the RV residuals around the best solution. The best-fit parameters are reported in Table 12.

Our analysis of HD 136352 converges to the detection of three planet orbiting this G4V star. With minimum masses of 4.8, 8.6 and 10.8 M_⊕, HD 136352 host three super-Earth on orbits ranging from 11 to 108 days. The two inner planets, with periods of 11.6 and 27.6 days, are close to a 5:2 resonance and, contrary to the two planets in HD 21693, with a period ratio larger than 2.5.

As mentioned above, an interesting signal is present in the residuals at 123 days. Despite the period being suspicious (1/3 of a year), we ran another MCMC trial including this fourth signal as an additional planet (Keplerian). The fit converged towards a non-eccentric signal with an amplitude of 0.65 m s⁻¹ and a period of 122.6 days. However, when comparing the three- and four-planet model solutions, the case with three planets is strongly favored, with a ΔBIC of 39.6. We also tried to add to the three-planet model a polynomial of the first or the second order to check if this 123-day signal could be due to a very-long period companion, whose orbit is not covered by the data, creating an alias at a period of one year or one of its harmonics. This approach, however, did not reduce the signal at 123 days and was disfavoured by a model comparison using the BIC. We also looked for other possibilities of the three planet scenario, fitting the aliases of the 11 and 27.6-day signals, 0.91 and 0.96 day, respectively. These other possibilities are also ruled out with difference in BIC >17. Finally, we tested the sensitivity of the GLS periodogram to outliers. We found that by simply removing two of them, for example JD 2455411 and 24556168, the amplitude of the peak found at 123 days goes above a p-value of 1%. Therefore, several arguments points in the direction of this interesting, albeit not conclusively significant, signal at 123 days being more likely an artefact induced by noise in the data or interaction with the window function rather than a bona-fide planet. We therefore keep for HD 136352 the three planet solution, with periods of 11.6, 27.6 and 108 days.

Fig. 21 Phase-folded RV measurements of HD 134060 with the best fitted solution for planets b and c represented as a black curve. The RV residuals around the best-fit solution are displayed in the lower panel. Corresponding orbital elements arelisted in Table 11.

Fig. 22 From top to bottom panels: GLS periodogram of the residuals at each step, after removing one planet after the other in the analysis of HD 136352. The GLS periodogram of the raw RVs is shown in Fig. 1.

6 Conclusion

We have reported the discovery and detailed characterization of eight planetary systems hosting twenty low-mass planets, plus a possible additional planet on a 6-month orbit, discovered with the HARPS Echelle spectrograph mounted on the 3.6-m ESO telescope located at La Silla Observatory. We also improved the characterization of the extremely eccentric Jupiter orbiting HD20782. Figure 24 illustrates these detections in a diagram showing the planet distribution in periods, with indication of their masses and eccentricities. The planets in our sample can be divided in three categories:

• Very small-mass planets with minimum masses below 6 M_⊕, that are found orbiting HD 20781 and HD 136352, on short periods, less than 15 days.

• More massive planets in the super-Earth to Neptune transition regime, spanning a range in minimum masses from 8 to 17 M_⊕, and a range in periods from a few days to nearly a year. They represent the bulk of our detections.

• Massive planets found on long-period orbits.

The RV technique sensitivity goes down when moving towards small-mass and long-period planets. It is therefore not surprising that most of our detections are in the second mass category. The lack of massive objects on short-period orbits is already well established, and the detection of extremely small-mass planets, although numerous from Kepler statistics (e.g. Coughlin et al. 2016) and previous RV surveys (e.g. Butler et al. 2017; Mayor et al. 2011), is challenging for RV surveys because of the perturbing effect of additional stellar signals (Dumusque et al. 2017) and the preponderance of appearance of these planets in multiple systems.

The detection of these planets had been announced in Mayor et al. (2011) studying statistical properties of the systems discovered with HARPS. However, with more data in hand, it was possible to discover new bona-fide planets. This is the case for the two inner super-Earth orbiting HD20781 and the 13.1-day period super-Earth orbiting HD45184. There is also the 184-day signal found in the timeseries of HD20003 even if, as explained in Sect. 5.1, we cannot exclude an instrumental origin to this signal. These new detections show that gathering more data helps in detecting small-mass planets on short-period orbits as well as long-period signals. Although expensive, characterizing those signals is therefore a necessary task to get an as complete as possible view of planetary systems.

The most important targets in today’s exoplanet studies are the transiting small-mass planets orbiting bright stars. They allow for exquisite determination of their orbital and physical properties, bringing priceless constrains to our understanding of planet formation, of the establishment of planetary system architecture, and of the diversity of planet internal compositions. Those candidates would also be excellent candidates for further atmospheric characterization. A dedicated Spitzer survey looking for potential transit events induced by the shortest-period planets in our sample was conducted by Gillon et al. (2017). Unfortunately, after searching for transits of the innermost planets orbiting HD 20003, HD 20781, HD 31527, HD 45184, HD 51608 and HD 134060, no detection was reported in this study. The search will continue, in particular with the Swiss-ESA led CHEOPS satellite to be launched end of 2018 (Fortier et al. 2014).

Some systems presented here are interesting in terms of formation process and architecture. HD 20003 hosts two Neptune-mass planets with periods of 11.9 and 33.9 days, thus close to a 3:1 commensurability. This configuration might explain the relatively high eccentricity measured on the innermost planet in this system, especially if some resonance is involved in the process. We note also that the two Neptune-mass planets orbiting HD 21693 and the two innermost super-Earth found around HD 136352 are close to a 5:2 resonance. Finally, HD 134060 is also an interesting dynamical system as it harbors a small-mass planet on a fairly eccentric 3-day orbit, accompanied by a more massive long-period planet at ~1300 days. Several possibilities for the formation of such configurations have been mentioned in the text. The list is probably not exhaustive and further theoretical work is need to analyse these systems in a deeper way. Here we just checked that the stability of the proposed solutions is verified over a few tens of thousands of years, avoiding the systems to quickly fall apart.

Fig. 23 Phase-folded RV measurements of HD 136352 with, from left to right and top to bottom the best fitted solution for planets b–d. The residuals around the solution are displayed in the lower right panel. Corresponding orbital elements are listed in Table 12.

Fig. 24 Summary of the confirmed planetary systems with indication of their periods (x-axis), masses (symbol size) and eccentricities (colour code).

Acknowledgements

We thank the Swiss National Science Foundation (SNSF) and the Geneva University for their continuous support to our planet search programs. This work has been in particular carried out in the frame of the National Centre for Competence in Research “PlanetS” supported by SNSF. X.D. is grateful to the Society in Science–The Branco Weiss Fellowship for its financial support. N.C.S. was supported by Fundação para a Ciência e a Tecnologia (FCT, Portugal) through the research grant through national funds and by FEDER through COMPETE2020 by grants UID/FIS/04434/2013 & POCI-01-0145-FEDER-007672 and PTDC/FIS-AST/1526/2014 & POCI-01-0145-FEDER-016886, as well as through Investigador FCT contract nr. IF/00169/2012/CP0150/CT0002 funded by FCT and POPH/FSE (EC). P.F. acknowledges support by Fundação para a Ciência e a Tecnologia (FCT) through InvestigadorFCT contracts of reference IF/00169/2012/CP0150/CT0002, respectively, and POPH/FSE (EC) by FEDER funding through the program “Programa Operacional de Factores de Competitividade - COMPETE”. P.F. further acknowledges support from Fundação para a Ciência e a Tecnologia (FCT) in the form of an exploratory project of reference IF/01037/2013/CP1191/CT0001. A.C.C. acknowledges support from STFC consolidated grant ST/M001296/1. This research has made use of the SIMBAD database and of the VizieR catalogue access tool operated at CDS, France, and used the DACE platform developed in the frame of PlanetS (https://dace.unige.ch).

Appendix A: Periodogram of the activity and CCF indicators

Fig. A.1 From top to bottom panels: periodograms of the log (R${}_{HK}^{\prime }$), BIS SPAN, and FWHM residuals of HD20003 after fitting either a Keplerian to adjust at best the observed magnetic cycle or a second-order polynomial to take into account any drift that could be instrumental. Planetary signals announced in Sect. 5 are represented by dashed vertical red lines.

To check if any signal detected in RVs corresponds to signals measured in activity indicator, we show here the Lomb-Scargle periodograms of the log (R${}_{HK}^{\prime }$), the BIS SPAN, and the FWHM. Those indicators have been shown to be sensitive to activity, and therefore any signal appearing both in the RVs and at least one of those indicators might be induced by stellar activity. Because magnetic cycles will be seen as a long-period significant signals in all those indicators, we removed any long-term signal either by fitting a Keplerian to adjust at best the observed magnetic cycle, as in Fig. 1, or by adjusting a second-order polynomial.

Fig. A.2 Same as Fig. A.1 but for HD20781.

Fig. A.3 Same as Fig. A.1 but for HD21693.

Fig. A.4 Same as Fig. A.1 but for HD31527.

Fig. A.5 Same as Fig. A.1 but for HD45184.

Fig. A.6 Same as Fig. A.1 but for HD51608.

Fig. A.7 Same as Fig. A.1 but for HD134060.

Fig. A.8 Same as Fig. A.1 but for HD136352.

Appendix B: Parameters probed by MCMC

We present here the parameter estimates issued from the MCMC modeling of the planetary systems.

References

All Tables

All Figures

Fig. 1 Left column, from top to bottom panels: HARPS RV measurements as a function of barycentric Julian date obtained for HD 20003, HD 20781, HD 21693, HD 31527, HD 45184, HD 51608, HD 134060 and HD 136352. Error bars only include photon and calibration noise. Middle column: log (R${}_{HK}^{\prime }$) activity indicator as a function of time for the same stars. Right column: GLS periodogram of the corresponding RV measurements.

In the text

	Fig. 2 GLS Periodogram of the RV residuals of HD 20003 at each step after removing, from top-left to bottom-right, the magnetic cycle effect plus a second-order polynomial drift, then the two planets at 11 and 33 days one after the other, and finally the signal at 184 days. The GLS periodogram of the raw RVs is shown in Fig. 1.
In the text

	Fig. 3 GLS Periodogram of the RV residuals of HD 20003 obtained from the MCMC solution including a linear correlation between RV and log (R${}_{HK}^{\prime }$) to account for the magnetic cycle, a second-order polynomial drift and three Keplerian.
In the text

Fig. 4 Phase-folded RV measurements of HD 20003 with the best Keplerian solution represented as a black curve for each of the signals in the data. From top left to bottom right panels: planet b, planet c and the signal at 184 days. The residuals around the solution are displayed in the lower right panel. Contrarily to Fig. 1 error bars include now the instrumental + stellar jitter from the MCMC modeling. Corresponding planetary orbital elements are listed in Table 5.

In the text

Fig. 5 Left panel: activity index log (R${}_{HK}^{\prime }$) as a function of the RV residuals when removing all the detected signals except the magnetic cycle effect for HD 20003. The observed correlation indicates that most of the RV residual variation is due to activity-related effects. Right panel: barycentric Earth velocity projected on the line of sight (BERV) as a function of the RV residuals around the best derived solution without considering the 184-day signal. The fact that no correlation can be observed disfavours the hypothesis that the 184-day signal is due to a discontinuity in the wavelength calibration introduced by tiny gaps between the different quadrants of the detector.

In the text

	Fig. 6 GLS periodogram of the residuals at each step, after removing one planet after the other in the analysis of HD 20781 (from left toright and top to bottom). The GLS periodogram of the raw RVs is shown in Fig. 1.
In the text

Fig. 7 Phase-folded RV measurements of HD 20781 with the best-fit solution represented as a black curve for each of the signal foundin the data (from left to right and top to bottom: planets b–e). Error bars include photon and calibration noise, as well as instrumental + stellar jitters derived from the MCMC modeling. The residuals around the solution are displayed in the lower panel. Corresponding orbital elements are listed in Table 6.

In the text

	Fig. 8 Best Keplerian solution for the eccentric planet orbiting HD 20782. Top panel: phase-folded RV measurements with the best solution represented as a black curve. Bottom panel: RV residuals around the best Keplerian solution.
In the text

	Fig. 9 Top panel: GLS periodogram of the residual RVs of HD 21693 after removing the RV contribution of the magnetic cycle. Middle and bottom panels: GLS periodogram of the residuals at each step, after removing one planet after the other in the analysis. The GLS periodogram of the raw RVs is shown in Fig. 1.
In the text

	Fig. 10 Phase-folded RV measurements of HD 21693 with the best planet solution represented as a black curve for each of the signal in the data (top to bottom: planet b and planet c). The residuals around the solution are displayed in the lower panel.Corresponding orbital elements are listed in Table 7.
In the text

	Fig. 11 RV residuals when removing all the detected signals except the magnetic cycle effect plotted as a function of the activity index log (R${}_{HK}^{\prime }$) for HD 21693. The observed correlation indicates that most of the RV residual variation is due to activity-related effects.
In the text

	Fig. 12 Top panel: GLS periodogram of the residual RVs of HD 31527 after removing the best-fit Keplerian to account for the significant signal at 17 days seen in the raw RVs (bottom right panel of Fig. 1). Middle and bottom panels: GLSperiodogram of the residuals at each step, after removing the second and third planet from the RVs.
In the text

	Fig. 13 Phase-folded RV measurements of HD 31527 with the best Keplerian solution for planets b–e represented as black curves (fromleft to right and top to bottom). Error bars include photon and calibration noise, as well as a jitter estimated from the MCMC analysis. The residuals around the solution are displayed in the bottom right panel. Corresponding orbital elements are listed in Table 8.
In the text

	Fig. 14 From top to bottom panels: GLS periodogram of the RVs after removing the effect induced by the stellar magnetic cycles and then one planet after the other in the analysis of HD 45184. The GLS periodogram of the raw RVs is shown in Fig. 1.
In the text

	Fig. 15 Phase-folded RV measurements of HD 45184 with, from top to bottom, the best Keplerian solution for planets b and c, and the corresponding residuals. Error bars include photon and calibration noise, as well as a jitter effect (stellar + instrumental) determined in the MCMC analysis. Corresponding orbital elements are listed in Table 9.
In the text

	Fig. 16 RV residuals when removing all the detected signals except the magnetic cycle effect plotted as a function of the activity index log (R${}_{HK}^{\prime }$) for HD 45184. The observed correlation indicates that most of the RV residual variation is due to activity-related effects.
In the text

	Fig. 17 From top to bottom panels: GLS periodogram of the RV residuals of HD 51608 after removing the effect of the magnetic cycle and, at each step, after removing one planet after the other. The GLS periodogram of the raw RVs is shown in Fig. 1.
In the text

	Fig. 18 Phase-folded RV measurements of HD 51608 with, from top to bottom the best Keplerian solution for planets b, c and the residuals around the best-fitted solution. Corresponding orbital elements are listed in Table 10.
In the text

	Fig. 19 RV residuals around the best derived solution without considering the magnetic cycle effect plotted as a function of the activity index log (R${}_{HK}^{\prime }$) for HD 51608. The observed correlation indicates that most of the RV residual variation is due to activity-related effects.
In the text

	Fig. 20 From top to bottom panels: GLS periodogram of the residuals at each step, after removing one planet after the other in the analysis of HD 134060. The GLS periodogram of the raw RVs is shown in Fig. 1.
In the text

	Fig. 21 Phase-folded RV measurements of HD 134060 with the best fitted solution for planets b and c represented as a black curve. The RV residuals around the best-fit solution are displayed in the lower panel. Corresponding orbital elements arelisted in Table 11.
In the text

	Fig. 22 From top to bottom panels: GLS periodogram of the residuals at each step, after removing one planet after the other in the analysis of HD 136352. The GLS periodogram of the raw RVs is shown in Fig. 1.
In the text

	Fig. 23 Phase-folded RV measurements of HD 136352 with, from left to right and top to bottom the best fitted solution for planets b–d. The residuals around the solution are displayed in the lower right panel. Corresponding orbital elements are listed in Table 12.
In the text

	Fig. 24 Summary of the confirmed planetary systems with indication of their periods (x-axis), masses (symbol size) and eccentricities (colour code).
In the text

	Fig. A.1 From top to bottom panels: periodograms of the log (R${}_{HK}^{\prime }$), BIS SPAN, and FWHM residuals of HD20003 after fitting either a Keplerian to adjust at best the observed magnetic cycle or a second-order polynomial to take into account any drift that could be instrumental. Planetary signals announced in Sect. 5 are represented by dashed vertical red lines.
In the text

	Fig. A.2 Same as Fig. A.1 but for HD20781.
In the text

	Fig. A.3 Same as Fig. A.1 but for HD21693.
In the text

	Fig. A.4 Same as Fig. A.1 but for HD31527.
In the text

	Fig. A.5 Same as Fig. A.1 but for HD45184.
In the text

	Fig. A.6 Same as Fig. A.1 but for HD51608.
In the text

	Fig. A.7 Same as Fig. A.1 but for HD134060.
In the text

	Fig. A.8 Same as Fig. A.1 but for HD136352.
In the text