© ESO 2018

1 Introduction

In the last decade, the search for potentially habitable exoplanets has focussed particularly on M-dwarf stars. In this context, we regard a potentially habitable exoplanet as, in the broadest sense, one that orbits within or close to the habitable zone (HZ) of the parent star. Using the two main indirect detection methods, photometrictransits and radial velocity (RV), potentially habitable exoplanets can be detected more easily around M dwarfs than around earlier type stars.

However, even M dwarfs present some challenges because of their faintness, their magnetic activity, and the difficulty in measuring accurate stellar parameters, which are essential if we wish to accurately determine the planetary parameters. Several exoplanet surveys have been devised to target M dwarfs specifically, such as the current ground-based photometric experiments MEarth (Irwin et al. 2015), APACHE (Sozzetti et al. 2013), and TRAPPIST (Gillon et al. 2017); the upcoming SPECULOOS¹ and ExTrA (Bonfils et al. 2015) projects; and surveys exploiting high-resolution and high-stability spectrographs (e.g. Bonfils et al. 2013; Delfosse et al. 2013b; Affer et al. 2016 and the HADES paper series; Quirrenbach et al. 2016).

Results from the Kepler and K2 missions have also been used to provide estimates for the occurrence of planets in the HZ of M dwarfs (e.g. Dressing & Charbonneau 2013, 2015; Howell et al. 2014; Dressing et al. 2017a, b). These analyses suggest that small, low-mass planets are abundant around M dwarfs, and a significant percentage of these planets may have properties suitable for the emergence of life (Dressing & Charbonneau 2015; Tuomi et al. 2014). Space-based missions planned for the very near future, such as TESS (Ricker et al. 2014), CHEOPS (Fortier et al. 2014), and PLATO (Rauer et al. 2014) will target bright and nearby M dwarfs to detect potentially habitable planets, while JWST (Beichman et al. 2014) and ground-based 25–40 m class telescopes will characterize the atmospheres of these planets.

Some of the most intriguing exoplanet discoveries of the last few years are temperate planets around nearby M dwarfs, which orbit close to, or within, the predicted circumstellar HZ. For example, temperate, low-mass planets have been discovered around a number ofnearby M-dwarf stars including Proxima Centauri (Anglada-Escudé et al. 2016), TRAPPIST-1 (Gillon et al. 2017), LHS 1140 (Dittmann et al. 2017), GJ 273 (Astudillo-Defru et al. 2017), K2-18 (Cloutier et al. 2017), GJ 667 C (Anglada-Escudé et al. 2013; Delfosse et al. 2013a; Feroz & Hobson 2014), and Ross 128 (Bonfils et al. 2018). These new discoveries have stimulated many theoretical studies addressing the habitability of planets orbiting M dwarfs (see e.g. Shields et al. 2016 for a review on this topic). Open questions regarding the true habitability of these temperate planets include the influence of the spectral energy distribution and the activity of M dwarfs on planetary atmospheres, as well as the effect of tidal locking. Potentially habitable planets around M dwarfs are likely subject to physical conditions that are very different from those experienced on Earth, and their properties may also depend strongly on the M-dwarf spectral subtype.

Indeed, M dwarfs that host temperate rocky planets exhibit a wide range of physical properties that might influence habitability. GJ 273, Ross 128, LHS 1140, Proxima Centauri, and TRAPPIST-1 are mid- to late-M dwarfs (M3.5V, M4V, M4.5V, M5.5V, and M8V, respectively). GJ 273 and Ross 128 have weak magnetic activity, while Proxima Centauri and TRAPPIST-1 have high activity levels, which probably has a large impact on the potential habitability of their planets. K2-18 is, instead, an earlier type M2.5 dwarf with low chromospheric activity, and GJ 667 C is a quiet M1.5 dwarf. The fact that these temperate planets are located at different distances from host stars spanning different spectral subtypes makes their comparative characterization particularly interesting. Unfortunately, many of these planets either do not transit or their parent star is too faint to permit detailed follow-up, making it difficult to robustly characterize these planets and hence study their interior structures and compositions.

Within this context, the planetary system around K2-3 (EPIC 201367065), a nearby (45 pc) M0 dwarf (V = 12 mag; J = 9.4 mag), presents an interesting opportunity for follow-up studies. Observations from the K2 mission revealed that K2-3 hosts at least three transiting small planets (Crossfield et al. 2015): K2-3 b (R_p = 2 R_⊕, P_orb = 10 days), K2-3 c (R_p = 1.7 R_⊕, P_orb = 24.6 days), and K2-3 d (R_p = 1.6 R_⊕, P_orb = 45.5 days). According to the optimistic HZ boundaries derived by Kopparapu et al. (2013, 2014), planet K2-3 d orbits close to the inner edge of the HZ of its host star. A particularly intriguing property of the K2-3 system is that the host star is bright enough to estimate the masses of these planets using existing high-resolution stabilized spectrographs.

Measuring the masses of the K2-3 planets would be interesting for several reasons. First, by determining the mass and bulk density of the temperate planet K2-3d, we can extend the study of planets orbiting close or within the HZ to earlier type host stars than those discussed earlier. Moreover, because of their measured sizes, planets K2-3 c and K2-3 d are optimal targets to test the results of Rogers (2015), who found that the majority of the observed planets with radius R_p ≥ 1.6 R_⊕ have densities too low to have a bulk rocky composition. Fulton et al. (2017) cast light on this result by finding evidence for a bimodal distribution of small planet sizes orbiting stars with T_eff > 4500 K with periods less than 100 days. Planets are preferentially found with radii around ~1.3 R_⊕ and ~2.4 R_⊕, with a region between 1.5 and 2.0 R_⊕ where planet occurrence is rare. This bi-modality likely reflects a separation between purely rocky planets and rocky cores surrounded by varying amounts of lower density volatile material. K2-3 c and K2-3 d fall within the gap between the two peaks in the planet radius distribution and orbit a star cooler than 4000 K, making them interesting test cases for understanding the composition of planets in this size regime at lower stellar irradiation than most of the planets analysed by Fulton et al. (2017).

Previously, other groups have recognized the appeal of the K2-3 system and have begun conducting RV follow-ups (Almenara et al. 2015; Dai et al. 2016). Almenara et al. (2015) found that stellar activity has a major impact on RV observations of K2-3 and, based on their single semester of monitoring, were unable to measure the masses of K2-3 c and K2-3 d robustly. Accurate determination of the masses of K2-3 c and K2-3 d evidently requires a larger dataset and denser sampling to trace out the activity signal.

We present the results of an intense RV follow-up of K2-3 conducted over three seasons with the HARPS-N and HARPS spectrographs. Despite having data from several independent teams, and having precise transit ephemeris from high-precision photometry (Beichman et al. 2016; Fukui et al. 2016), measuring the mass of K2-3 c and K2-3 d has been challenging.

The paper is organized as follows. We first introduce the RV datasets and discuss the significant signals present in the data through a frequency analysis (Sect. 2). In Sect. 3 we present updated stellar parameters and analyse the photometric light curves of K2-3 and the Hα line activity indicator time series. The analysis of the RVs with Gaussian processes (GP) is described in Sect. 4. The significance of our best-fit solution for the masses of K2-3 c and K2-3 d is investigated through Monte Carlo simulations and is discussed in Sect. 5. Finally, we present the mass-radius diagram for the K2-3 planets and discuss the implications of our findings for their bulk composition.

2 Description and first look on the HARPS(-N) radial velocity datasets

Northern and southern spectroscopic observations of K2-3 were carried out between January 22, 2015 and July 7, 2017, producing a total of 211 HARPS-N and 138 HARPS spectra. Of the 349 spectra, 283 are unpublished observations while 66 HARPS observations were previously published by Almenara et al. (2015). The HARPS-N spectra come from two independent programmes: the HARPS-N Collaboration Guaranteed Time Observations (GTO)² and the Global Architecture of Planetary Systems programmes (GAPS; Benatti et al. 2016). The two collaborations shared observing time on this target to maximize the number of RV measurements and to optimize the observing strategy. The spectra were reduced with the version 3.7 of the HARPS-N Data Reduction Software (DRS; Cosentino et al. 2014). The extraction and wavelength calibration of the HARPS spectra were performed using the on-line pipeline (Lovis & Pepe 2007). The exposure time was fixed at 1800 s for both instruments, producing a typical signal-to-noise ratio S∕N = 30 at a wavelength λ ~ 550 nm.

2.1 Definition of the final dataset

We excluded some observations from our analysis for two reasons. We did not consider spectra with a S∕N ≤ 11 as measured at λ = 550 nm (echelle orders 46 and 49 for HARPS-N and HARPS, respectively). We then identified RV measurements potentially contaminated by scattered moonlight. This is particularly important for faint targets such as K2-3, which lie close to the ecliptic plane. Using the procedure described in Sect. 2.1 of Malavolta et al. (2017), we identified five HARPS-N spectra that are probably contaminated by scattered moonlight, and we discarded them from the dataset³. We could not perform the same analysis on the full HARPS dataset because nearly half of those spectra were acquired with simultaneous Fabry-Perot and not with fibre B on sky. Nonetheless, we adopted a heuristic approach⁴ to identify four potentially contaminated measurements, which we removed from the final dataset⁵. After these cuts, our final dataset includes 197 HARPS-N spectra and 132 HARPS spectra.

2.2 Radial velocity extraction

In this study we primarily use RVs extracted with the TERRA pipeline (Anglada-Escudé & Butler 2012). The TERRA pipeline is commonly used to extract RVs of M dwarfs because it typically provides measurements with better precision and lower scatter than those extracted for low-mass stars with the CCF recipe of the on-line DRS pipeline (Perger et al. 2017). The HARPS fibre link was upgraded with octagonal fibres on May 28, 2015 (Lo Curto et al. 2015), introducing an RV offset between data acquired before and after the upgrade. We accounted for this offset by producing pre- and post-upgrade spectral templates to extract the RVs, and by introducing a velocity zero point offset for each dataset as free parameter when fitting the time series.

The TERRA RV time series are listed in the on-line Tables A.1 and A.2 for HARPS-N and HARPS, respectively. The median values of the internal errors are 1.88 m s⁻¹ for HARPS-N and 2.04 m s⁻¹ for HARPS.

2.3 Preliminary RV frequency analysis

We conducted a preliminary analysis of the RV datasets using the generalized Lomb-Scargle (GLS) algorithm (Zechmeister & Kürster 2009) to identify significant signals. We performed the GLS frequency analysis of the northern and southern datasets separately along with the combined RV dataset, removing RV offsets as determined by the analysis described inSect. 4. We show the resulting periodograms in Figs. 1–3.

The HARPS-N periodogram is dominated by a signal at 37.2 ± 0.1 days⁶, which we identify as the rotation period of the star (see Sect. 3). This signal has a p-value = 0.01%, as estimated by boot strapping the data, i.e. by randomly drawing the RV measurements (with replacement) and generating 10 000 mock datasets. Its semi-amplitude is 2.9 ± 0.3 m s⁻¹, as estimated with GLS, slightly lower than the RMS of the data (4.1 m s⁻¹). A significant peak with a slightly lower power appears at the orbital period of K2-3 b, while the peak corresponding to the orbital period of K2-3 c appears with less significance (p-value ~1%). The orbital period of K2-3 d is undetected in the HARPS-N periodogram.

The HARPS periodogram is dominated by the orbital period of K2-3 b (Fig. 2; p-value = 0.01%). The signal produced by the stellar rotation has a p-value around 1%, therefore it is much less significant than in the HARPS-N data. The window function is responsible for a pattern of alias frequencies around these signals. In addition to the one-year aliases, one-month aliases (synodic month frequency f_s.m. = 0.03386 c/d) are also probably present, which is expected for a star near to the ecliptic. For example, a one-month alias of the 37-day signal occurs at 0.061 c/d (=1/37+f_s.m. c/d), while a one-monthalias of P_b occurs at 0.065 c/d. Both aliases are visible in Fig. 2.

The signals that most clearly emerge in the periodogram of the combined dataset are those of the planet K2-3 b and the stellar rotation period, which are highly significant and have almost the same power (Fig. 3). The signal of planet K2-3 c is more weakly present with a p-value ~1%.

Fig. 1 Upper plot: time series of the K2-3 RVs extracted with the pipeline TERRA using the HARPS-N spectra. Middle plot: GLS periodogram of the RV time series (red line) is shown. Three levels of p-values are indicated by the horizontal dashed lines. The green dashed line overplotted on the RV periodogram represents the window function of the measurements (bottom plot) shifted in frequency to be superimposed on the strongest peak of the RV periodogram. This helps to identify the alias frequencies of the most relevant peaks. The stellar rotation and planetary orbital frequencies are indicated by vertical lines and corresponding labels.

Fig. 2 As in Fig. 1 but for TERRA RVs extracted from the HARPS spectra. Pre- and post-upgrade RV offsets have been applied, as derived from our analysis described in Sect. 4.

Fig. 3 As in Fig. 1 but for all TERRA HARPS-N and HARPS RVs. We applied an offset to each separate dataset, as derived from our analysis described in Sect. 4.

3 Stellar parameters and activity

The main properties of the star K2-3 are summarized in Table 1. While K2-3 has previously been characterized by Crossfield et al. (2015) and Almenara et al. (2015), we independently determined the stellar parameters from our HARPS-N spectra using the method developed by Maldonado et al. (2015)⁷. We estimated the effective temperature and iron abundance of the star using measurements of the pseudo-equivalent widths of the spectral features. We then determined stellar mass, radius, and surface gravity using empirical relationships. Our results are consistent with, but slightly more precise than, previous estimates. We adopted these new values for the analysis in the rest of the paper.

3.1 Photometry

The K2 light curve of K2-3, with all the transit signals removed, is shown in Fig. 4. We used the K2SFF⁸ light curve processed asdescribed by Vanderburg & Johnson (2014) and Vanderburg et al. (2016a). The light curve shows a quasi-periodic modulation with a flux semi-amplitude of ~0.1%. There is also clear evidence for changes from one rotation to the next, likely due to the evolution of active regions. In the bottom panel of Fig. 4 we show the GLS periodogram of the binned K2 light curve (one point per day), which shows a strong peak at P_rot = 38.3 ± 0.7 days⁹. Even though the K2 light curve was obtained about six months before the first spectroscopic observations and only covers about two rotationperiods, the high S/N photometry is useful for constraining the stellar rotation period. It is interesting to note that, while the amplitude of the photometric rotational variability of the star is fairly low, the stellar rotation frequency is nonetheless the strongest signal in the RV time series of the star.

We also obtained time-series photometry of K2-3 from the MEarth survey (Irwin et al. 2015). We monitored K2-3 with one of the 40 cmtelescopes of the Southern MEarth array from January 21, 2015 to April 4, 2016. With a total of 8669 data points, the light curve overlaps with the first two seasons of the spectroscopic observations. The fact that the MEarth data are contemporaneous with many of our RV observations and their high photometric precision (due in part to the excellent observing conditions at Cerro Tololo), make the MEarth light curve potentially useful for characterizing the stellar activity. We analysed the MEarth data (in nightly bins) by calculating a GLS periodogram. The highest peak in the periodogram is at 31.8 ± 0.2 days with a semi-amplitude of ~1 mmag (Fig. 5). This period does not correspond to the more robust stellar rotation period derived from the uninterrupted observations of K2 or using the H_α spectroscopic indicator (Sect. 3.2); nonetheless the result is worth mentioning because the data have been collected with unevensampling by a ground-based small-aperture telescope and they show very low amplitude modulation compared to the typical photometric error. The periodogram also shows an additional peak at a period of about ~230 days, without any counterpart in the window function. While the signal could be due to systematics present in the MEarth data, its nature is unclear and an astrophysical origin cannot be ruled out (see Sect. 3.2).

Fig. 4 Top plot: K2 light curve of K2-3 with the planetary transits removed. Bottom plot: GLS periodogram of the binned light curve (one point per day), showing a peak at Prot = 38.3 days.

Fig. 5 Top plot: light curve of K2-3 from the MEarth-south survey, folded at the best period P = 31.8 days found using GLS. Middle plot: GLS periodogram of the MEarth light curve. The green line corresponds to the window function of the measurements (bottom plot), which is shifted in frequency so that the peak is superimposed on the major peak of the RV periodogram to directly identify alias frequencies.

3.2 Spectroscopic activity indexes

We studied the activity level of K2-3 during the time period of our observations using the Hα line as a spectroscopic activity indicator, which we extracted from our spectra using the method described by Gomes da Silva et al. (2011) 10. The time series of the Hα activity indicator is shown in Fig. 6 and the data are listed in Table A.5. We analysed the time series by calculating a GLS periodogram to identify periodicities related to a possible activity cycle and stellar rotation. The GLS periodogram is shown in the second and third panels of Fig. 6.

The highest peak in the periodogram occurs at P = 211 ± 3 days, and has a power comparable to that of its one-year alias frequency at ~ 450 days. The existence of a long-term modulation is especially clear after looking at the data of the first two seasons. A bootstrap (with replacement) Monte Carlo analysis based on 10 000 mock datasets reveals that these peaks are statistically significant and have false alarm probabilities lower than 0.1%. We are, however, unable to ascertain whether the 211-day period or its alias at 450 days is the true underlying period. The origin of this long-period signal is unclear. If the signal is astrophysical, one possible explanation is an intermediate-duration activity cycle. Some tentative evidence exists for such cycles in low-mass stars (Savanov 2012; Robertson et al. 2013)and could represent sub-cycles superimposed on longer duration activity cycles, as observed for the Sun (so-called Rieger cycles). We note that the 211-day period is close to the 230-day signal observed in the MEarth photometry.

The second highest peak in the periodogram of the Hα indicator is close to the expected stellar rotation period and is highly significant (P = 40.3 ± 0.1 days, p-value < 0.1%). This signal is particularly strong in the last season of observations; during this time, a clear modulation related to P_rot is visible in the Hα time series, which covers nearly four stellar rotations. A GLS analysis of only the last season of observations identifies a periodicity of 43.5 ± 1.0 days. Folding the data at this period reveals that the modulation does not have a simple sinusoidal shape (Fig. 7).

Fig. 6 Top plot: time series of the activity indicator based on the Hα line extracted from the HARPS and HARPS-N spectra. Second and third plots: GLS periodogram of the dataset is shown. The dashed horizontal lines indicate the p-value levels as derived from a bootstrap analysis. The greencurve corresponds to the window function of the measurements, which is shifted in frequency so that the peak is superimposed on the major peak of the RV periodogram (second plot), to directly identify alias frequencies. In the third plot the window function is shifted to be superimposed on the P ~ 40 days rotational signal. Bottom: time series phase-folded at the period P = 211 days. Red points represent the average of the data within 20 bins in the phase range [0,1].

Fig. 7 Top plot: time series of the activity indicator based on the Hα line extracted from the HARPS and HARPS-N spectra. Here, only the dataset of the third season is shown. Middle plot: GLS periodogram of the dataset is shown. The green line corresponds to the window function of the measurements, which is shifted in frequency so that the peak superimposes to the major peak of the RV periodogram, to directly identify possible alias frequencies. Bottom plot: time series phase-folded at the period P = 43.5 days.

4 Gaussian process regression analysis of radial velocities

We used thestellar activity information derived in Sect. 3 to perform a detailed analysis of the combined HARPS-N and HARPS RV datasets within a Bayesian framework based on GP regression. It has now become standard in RV analysis to use GPs to model the stellar contributionto the RV variations jointly with a number of Keplerian functions describing the planetary orbital motion (see Haywood et al. 2014 for the first application of this technique). This approach has proven to be a powerful to retrieve planetary masses from high-precision RV measurements when a signal closely related to the stellar rotation period, or its harmonics, is present in the RV time series (Dumusque et al. 2017).

For a general description of the GP method, and its performance when applied to RV time series, we refer, amongst others, to the recent works by López-Morales et al. (2016), Cloutier et al. (2017), Damasso & Del Sordo (2017), and Dittmann et al. (2017). The GP regression is based on the choice of a specific kernel; i.e. a covariance matrix describing the correlation between measurements taken at two different epochs. For the case of K2-3, the so-called quasi-periodic (q-p) kernel is particularly useful because a signal very likely related to the stellar rotation period P_rot dominates the RV time series (Sect. 2), and there is evidence for an evolutionary timescale of the active regions close to P_rot (Sects. 3.1 and 3.2). Each element of the covariance matrix has the form $$\begin{array}{lll}K(t,t^{\prime })\hfill & =\hfill & h^2\cdot \exp [-\frac{{(t-t^{\prime })}^2}{2{\lambda }^2}-\frac{si{n}^2(\frac{\pi (t-t^{\prime })}{\theta })}{2w^2}]\hfill \\ \hfill & \hfill & +[{\sigma }_{\text{RV},\text{instr}}^2(t)+{\sigma }_{\text{jit},\text{instr}}^2]\cdot {\delta }_{t,t^{\prime }},\hfill \end{array}$$(1)

where t and t′ represent two different epochs. The first term represents the quasi-periodic kernel, which is composed of a periodic term and an exponential decay term. This functional form is suitable for modelling a recurrent signal linked to stellar rotation and takes into account the finite lifetime of the active regions. Eq. (1) contains four covariance matrix hyperparameters: h represents the amplitude of the correlations; θ represents the rotation period of the star; w is the length scale of the periodic component, linked to the size evolution of the active regions; and λ is the correlation decay timescale, which can be physically related to the active regions lifetime. The remaining parameters in Eq. (1) are σ_{RV, instr}(t), which is the RV internal error at time t for each spectrograph (or independent dataset); σ_{jit, instr}, which are additional uncorrelated jitter terms, one for each instrument (or independent dataset), which we add in quadrature to the internal errors to account for additional instrumental effects and noise sources neither included in σ_{RV, instr}(t) nor modelled by the q-p kernel; and ${\delta }_{t,t^{\prime }}$, which is the Kronecker delta function.

Our GP analysis was based on a Markov chain Monte Carlo (MCMC) algorithm. The model, algorithms, and statistical framework used in this work are the same as described in Damasso & Del Sordo (2017), and we refer to that work for a detailed description. In this work, we assumed circular orbits for all the planets. The best-fit values and uncertainties for each jump parameter were calculated as the median of the marginal posterior distributions and the 16% and 84% quantiles.

4.1 Choice of the priors

The priors adopted in our analysis are listed in Table 2. In this subsection, we describe and justify some of our choices for the priors.

In our fits, we allowed the semi-amplitude h of the correlated stellar signal to vary up to a value of 5 m s⁻¹, which represents nearly twice the value estimated by GLS at P = 37 days. The prior range used for the stellar rotation hyperparameter θ is defined based on our rotation period estimates from the K2 photometry and Hα activity indicator. This also takes into account the results of a trial MCMC analysis that was run after the conclusion of the second observing season, which adopted a larger prior range and indicated that the posterior distribution was well constrained within the range [35, 43] days.

For the active regions evolutionary timescale λ we adopted a uniform prior between 20 and 60 days, corresponding nearly to P_rot/2 and 1.5 ⋅ P_rot. This choice was first motivated by the marginal posterior obtained with the trial MCMC analysis of the dataset for the first and second seasons, which is symmetric around λ = 36 days (σ ~ 8 days). An evolutionary timescale of the order of the stellar rotation period could also be guessed directly by looking at the K2 light curve, which shows changes in its pattern from one rotation to the next¹¹. We also used the results of Giles et al. (2017; Eq. (8)) to get a loose estimate of the timescale based on the RMS of the K2 data and the stellar effective temperature. The equation for determining the timescale was derived by Giles et al. (2017) using a calibration sample not biased by spectral type and, although faster rotators than K2-3 were analysed in the Giles et al. work, this result could be tentatively used for a star with P_rot = 40 days. With an RMS = 7.6 × 10⁻⁴ mag and T_eff = 3835 K, we get a timescale of 38${}_{-13}^{+19}$ days.

For the planetary orbital parameters, we fixed the upper limit of the Keplerian semi-amplitude to 5 m s⁻¹ for K2-3 b and K2-3 c, and to 3 m s⁻¹ for planet K2-3 d. For K2-3 b, K = 5 m s⁻¹ is more than twice the value estimated by GLS after removing the stellar rotation signal, and represents a conservative upper limit also for the planet K2-3 c. For the outermost planet, K2-3 d, the absence of a signal in the GLS periodogram suggests that K_d should be significantly lower. The priors on the orbital period and time of transit were taken from Beichman et al. (2016), who determined the ephemerides by combining K2 observations with additional transits observed with the Spitzer space telescope.

4.2 Analysis of the combined radial velocity dataset

We analysed the full RV dataset without binning the data when more than one observation is available during the same night. To guarantee a wide exploration of the parameter space, we adopted 150 independent chains properly initialized to start from well separated locations. We discarded the first 3 000 steps of each chain by resetting the sampler. The MCMC chains reached the convergence according to the Gelman-Rubin statistics after 18 000 steps, and a further burn-in (0.75% of the total steps) was applied to calculate the best-fit values of the parameters (see Eastman et al. 2013 and references therein).

We show the best-fit results for all the free parameters and derived quantities in Table 3. The mass of the super-Earth K2-3 b is robustly determined with a S/N of ~7σ, while the mass of K2-3 c is detected with a lower significance of 2.6σ. The mass of K2-3 d remains undetermined, with a 1σ upper limit of about ~2M_⊕. This suggests that the mass of K2-3 d is too low to be detected, but it is also plausible that the low amplitude of the Doppler signal, compared to that of the stellar activity component, and its proximity to the stellar rotation period allow the signal to be absorbed into our GP activity model (Vanderburg et al. 2016b). We investigate the robustness of this upper limit in Sects. 5.1 and 5.2.

We show the RV curves in Fig. 8, with the stellar activity contribution subtracted, and folded at the best-fit orbital periods. The stellar contribution to the RV time series resulting from our fit is shown in Fig. 9. The stellar activity component of the model is reliably described by our best-fit results. The stellar rotation period and the evolutionary time scale of the active regions appear to be well characterized for both datasets. As we have found from the analysis of the RV datasets of other targets, the GP quasi-periodic regression tends to suppress low frequencies in the residuals. This means that if an additional long-term modulation is actually present in the original data, and it is not explicitly included in the fitted model, it would be absorbed by the stellar activity component and disappear in the residuals, after removing the planetary solutions. However, we note that the data in Fig. 9 do not show any long-term trend residual suggestive of an activity cycle or a longer period companion to K2-3 not included in our global model. We also note that the stellar rotation period, as retrieved by the GP regression, differs from that derived with GLS by ~3 days. This difference is because the functional form of the quasi-periodic kernel used to model the stellar contributionis different from the simple sinusoid used by the GLS periodogram. We calculated the GLS periodogram of the GP quasi-periodic stellar activity timeseries and found the strongest peak at P = 37 days, reproducing the period found in our GLS analysis of the RV dataset. The RMS of the residuals is 2.6 m s⁻¹, only slightlyhigher than the median of the internal errors.

Fig. 8 First row: TERRA RV residuals, after removing our best-fit stellar component, phase-folded to the three planetary solutions found for the quasi-periodic GP model (represented by a red curve). Second row: histograms of the number of data divided in bins of phase are represented, showing a fairly uniform coverage for each planet.

Fig. 9 Stellar signal contribution to the TERRA RV times series, for each observing season, as fitted with our GP quasi-periodic model (Table 3). The blue line indicates the best-fit curve, while the shaded grey area represents the ± 1σ confidence interval.

Fig. 10 TERRA RV residuals, after removing the planetary and stellar activity signals, as fitted within a GP quasi-periodic framework (Table 3).

4.3 Cross-check with the alternative RV extraction method

In order to test the robustness of the results using RVs obtained with the TERRA dataset, we also extracted the RVs using an independent pipeline described by Astudillo-Defru et al. (2015, 2017). In brief, this pipeline works by aligning all the spectra to a common reference frame by removing the Earth’s barycentric RV and by using the RVs of K2-3measured by the DRS as an initial guess. Then a median template is computed from the aligned spectra and regions of the spectra contaminated by tellurics are rejected. The template is used to calculate a chi-square profile as a function of RVs for each individual spectrum, whose minimum corresponds to the stellar RV.

The RVs from this alternative pipeline are listed in Tables A.3, A.4. We performed the same analysis as described in Sect. 4.2 with the data from this alternative pipeline, and found that the results of this analysis (shown in Table 3) are in very good agreement with those obtained with TERRA. In particular, the planetary parameters are all in agreement within 1σ, strengthening our findings.

5 Assessing the reliability of the derived planetary masses

The results of the analysis presented in Table 3 show that while the Doppler signal of the innermost planet is retrieved with a significance of ~7σ, the Doppler signal of K2-3 c is retrieved with a significance lower than 3σ, and K2-3 d is undetected. With the present dataset and the adopted model, we cannot deter- mine a bulk density of K2-3 c accurately enough to put reliable constraints on its possible composition, and we can only place an upper limit on the bulk density of K2-3 d. This impacts our understanding of the formation and evolutionary scenarios, which led to the observed system architecture, and makes inferences about the habitability of K2-3 d especially challenging.

Why is the characterization of this system so challenging, despite the large dataset we have collected with two world-class spectrographs? Can we rule out a rocky composition for K2-3 d? To provide answers to these questions we performed simulations to investigate the impact of the observing sampling, stellar activity, and internal RV uncertainties on our ability to retrieve masses for the K2-3planets.

5.1 Simulations using real epochs

We ran a first set of simulations based on the real epochs of our observations with the objective to improve the estimate of K2-3 d’s mass, for which we could provide only an upper limit, and to assess the robustness of our result for K2-3 c’s mass, despite its low significance. First, we adopted our best-fit values of the GP hyperparameters (TERRA dataset) to generate the stellar activity RV signal. We then injected planetary signals (circular orbits) with semi-amplitudes K_b = 2.7 m s⁻¹, K_c = 0.95 m s⁻¹, and K_d = 1 m s⁻¹, and the known periods and times of transit (Table 3) into the simulated activity signal. We refer to the dataset built in this way as the exact solution. While the simulated Doppler semi-amplitudes for K2-3 b and K2-3 c are those we derived from real data, for K2-3 d we assumed a value corresponding approximately to a purely rocky composition, given the measured radius¹².

We then created N = 156 mock RV time series¹³. Each dataset was obtained by randomly drawing RV values normally distributed around the exact solution, where ${\sigma }_{\text{RV}}^2(t)+{\sigma }_{\text{jit},\text{instr}}^2$ was used as the variance of the normal distributions at each epoch. We analysed each synthetic dataset following the same procedure as for the real data (Sect. 4.2) and recorded the best-fit values of the free parameters once the MCMC chains reached convergence.

Results from these simulations depend on the actual properties of the stellar activity observed during our campaign. More complex simulations to explore in detail the effects of the stellar activity could also be carried out, where each hyperparameter is randomly drawn within the uncertainties while keeping the others fixed. However, such simulations, which would require a large number of mock datasets (i.e. thousands) and a correspondingly huge amount of computational time, are beyond the scope of this paper. Thus, our simulations do not explore the possibility that the non-detection is due to the proximity of the stellar rotation period θ to the orbital period of K2-3 d, which could be a limiting factor (Vanderburg et al. 2016b). We note that the amplitude h of the stellar activity term is precisely known (~ 10σ, Table 3), so we expect that drawing values from the posterior distribution of this hyperparameter would not significantly change the results of our simulations.

Analysis framework. For each ith simulated dataset we derive a posterior distribution for the Doppler semi-amplitudes of planets c and d that we call K_{p, i}, where p = (c, d). Each semi-amplitude is characterized by a median value $K_{\text{p},\text{i}}^{\text{med}}$ and upper and lower uncertainties ${\sigma }_{\text{p},\text{i}}^{\text{+}}$ and ${\sigma }_{\text{p},\text{i}}^{\text{-}}$ (as derived from the 16th and 84th percentiles). We then calculate the median recovered semi-amplitude $K_{\text{p},\text{N}}^{\text{med}}$ of all the $K_{\text{p},\text{i}}^{\text{med}}$. We compare the median recovered semi-amplitude $K_{\text{p},\text{N}}^{\text{med}}$ with the injected value K_p,inj to draw conclusions about the results obtained for the real RV dataset.

For K2-3 d, we define the ratio r_d,i = (K_d,inj-$K_{\text{d},\text{i}}^{\text{med}}$)/${\sigma }_{\text{d},\text{i}}^{\text{+}}$ to measure the discrepancy between the best-fit estimate and the injected value in units of ${\sigma }_{\text{d},\text{i}}^{\text{+}}$. The term (K_d,inj-$K_{\text{d},\text{i}}^{\text{med}}$) is weighted by ${\sigma }_{\text{d},\text{i}}^{\text{+}}$ to take the skewness of each posterior distribution into account. By averaging r_d,i over the number, N, of simulated datasets we get the metric ${\overline{r}}_{\text{d}}$, that we propose as a way to correct the measured semi-amplitude K_d,meas using the equation $$K_{\text{d},\text{real}}=K_{\text{d},\text{meas}}+{\overline{r}}_{\text{d}}\cdot {\sigma }_{\text{d},\text{meas}}^{\text{+}},$$(2)

where K_d,meas and ${\sigma }_{\text{d},\text{meas}}^{\text{+}}$ come from our best-fit solution (Table 3).

As an alternative approach, for each marginal distribution K_d,i we calculate the percentile corresponding to the position of K_d,inj. We use th e median over N of these percentiles to derive an estimate for K_d,real from the posterior distribution obtained for the real dataset.

Results for K2-3c. The distribution of the median values $K_{\text{c},\text{i}}^{\text{med}}$ is shown in Fig. 11. The median of this distribution is $K_{\text{c},\text{N}}^{\text{med}}$ = 0.96${}_{-0.22}^{+0.27}$ m s⁻¹, where the uncertainties represent the 16^th and 84^th percentiles. Looking at the values for ${\sigma }_{\text{c},\text{i}}^{\text{+}}$ and ${\sigma }_{\text{c},\text{i}}^{\text{-}}$ over all the marginal posteriors K_{c, i}, we note that ⟨(${\sigma }_{\text{c},\text{i}}^{\text{+}}$-${\sigma }_{\text{c},\text{i}}^{\text{-}}$)⟩ = 0.01 m s⁻¹, indicating that the distributions are normal-shaped and their average is $\langle {\sigma }_{\text{c},\text{i}}^{\text{+}}\rangle$ = $\langle {\sigma }_{\text{c},\text{i}}^{\text{-}}\rangle$ = 0.33 m s⁻¹. These results show that the injected signal K_c,inj = 0.95 m s⁻¹ is well recovered, and indicate that our estimate of K_c from the real dataset is reliable, despite our detection having a significance <3σ.

Results for K2-3d. A sample of the N posterior distributions K_d,i is shown in Fig. 12. The distribution of $K_{\text{d},\text{i}}^{\text{med}}$ is shown in Fig. 13. The median of this distribution is $K_{\text{d},\text{N}}^{\text{med}}$ = 0.54${}_{-0.14}^{+0.18}$ m s⁻¹, while we get $\langle {\sigma }_{\text{d},\text{i}}^{\text{+}}\rangle$ = 0.47 m s⁻¹ and $\langle {\sigma }_{\text{d},\text{i}}^{\text{-}}\rangle$ = 0.35 m s⁻¹ for the 68.3% confidence interval of each distribution. Therefore, the semi-amplitude K_d is generally underestimated with respect to K_d,inj. This result is suggestive when compared to what we get for K_c. As shown in Fig. 8, the phase coverage of our data is uniform for planet K2-3 c, and fairly uniform for planet K2-3 d. Then, we would expect two signals with an equal semi-amplitude (here K_c,inj ≃ K_d,inj) to be recovered with similar significance in absence of other hampering factors. For the metric ${\overline{r}}_{\text{d}}$ we find ${\overline{r}}_{\text{d}}$ = 0.99 ± 0.04, where the error is calculated as RMS[r_d,i]/$\sqrt[]{N}$. We use this result and our best-fit value for K_d,meas = 0.29${}_{-0.18}^{+0.34}$ to draw N = 10 000 random values for ${\overline{r}}_{\text{d}}$ and K_d,meas, obtaining a distribution of N samples for K_d,real from Eq. (2). By taking the median of this distribution and the 68.3% confidence interval we get K_d,real = 0.63${}_{-0.18}^{+0.32}$ m s⁻¹. This corresponds to M_d,real = 2.7${}_{\text{-0}.\text{8}}^{\text{+1}.\text{2}}$ M_⊕ and ρ_d,real = 3.1${}_{\text{-1}.\text{2}}^{\text{+2}.\text{0}}$ g cm^-3 for the planet mass and density. Using the second approach mentioned above, we find that K_d,real corresponds on average at the 84^th percentile of the K_d posterior distribution for the real dataset, that is K_d,real = 0.62 m s⁻¹, which is a result equal to that obtained with the first method.

Our statistical analysis shows that K_d,real is <1 m s⁻¹ with 1σ confidence. This in turn suggests that the planets in the K2-3 system may have very similar bulk densities and thus share a similar composition, although this outcome should be taken with caution because of the large uncertainties in the masses and densities for K2-3 c and K2-3 d.

Fig. 11 Distribution of N = 156 median values for the Doppler semi-amplitude Kc of planet K2-3 c, normalized to the injected value Kc, inj, as derived from posterior distributions of the N mock RV datasets. This result refers to the case of real observing epochs.

Fig. 12 Sample of posterior distributions (grey histograms) for the Doppler semi-amplitude Kd of planet K2-3 d obtained from the GP regression analysis of the mock RV datasets described in Sect. 5.1. Each plot shows the marginal distribution for a single mock dataset, and this is compared to the posterior distribution obtained forthe real TERRA RV dataset, represented by the yellow histogram. The vertical dashed lines indicate the 50th (green) and 68.3th (black) percentiles for the posterior distributions of the mock datasets.

Fig. 13 Distributions of the 50th percentiles for all the posterior distributions of the semi-amplitude Kd obtained from 156 mock datasets, normalized to the injected value Kd, inj. This result refers to the case of real observing epochs.

5.2 Role of the observing sampling

We devised new simulations to investigate how the detection significance of the signals induced by the planets K2-3 c and K2-3 d would change by increasing the number of the RV data collected with HARPS-N and HARPS, still assuming K_d = 1 m s⁻¹. We simulated an intensive observing strategy conducted during the third season, which has the lowest amount of real data, in a similar way as carried out for the high-cadence campaign devised to detect Proxima b with HARPS (Anglada-Escudé et al. 2016). In order to keep our simulation realistic, we did not include mock epochs later than the 2017 observing season because our representation of the RV stellar signal cannot be considered predictive in the far future.

We created the mock datasets as described in Sect. 5.1. We generated all the epochs suitable for observationsduring the 2017 observing season in addition to the real epochs. We simulated only one measurement per night avoiding superposition with the epochs corresponding to the real observations. Every random epoch was selected by placing constraints on the Moon phase illumination and distance from the target (they have to be <90% and >45°, respectively), and on the altitude of the star above the horizon (airmass <1.7). Using these criteria we got 112 and100 additional new epochs for HARPS-N and HARPS, respectively. We randomly removed 10% of these epochs at each simulation run to account for bad weather, thus simulating an optimistic scenario for a feasible follow-up¹⁴. The uncertainties σ_RV(t) of the mock RV data were randomly drawn from normal distributions with mean and σ equal to the average and RMS values of the HARPS-N and HARPS post-upgrade internal errors derived with TERRA. The final mock dataset is obtained by randomly shifting each data point of the exact solution within the error bars by a quantity ΔRV(t) drawn from a normal distribution with mean zero and σ equal to $\sqrt[]{{\sigma }_{\text{RV}}^{\text{2}}(t)+{\sigma }_{\text{jit}}^{\text{2}}}$. An example mock dataset is show in Fig. 14.

Our final sample is composed of N = 97 mock datasets. Also in this case, we analysed each simulated dataset within the same GP quasi-periodic framework applied to the real dataset, except for the σ_jit terms that were not included as free parameters, and we analysed the outcomes as described in Sect. 5.1.

Results for K2-3c. The median and 68.3% confidence interval of the distribution for the K_c,i semi-amplitudes are K_c = 0.96${}_{\text{-0}.\text{26}}^{\text{+0}.\text{27}}$ m s⁻¹. For the upper and lower uncertainties we get ⟨(${\sigma }_{\text{c},\text{i}}^{\text{+}}$-${\sigma }_{\text{c},\text{i}}^{\text{-}}$)⟩ = 0.006 m s⁻¹ over all the posterior distributions, indicating that they are generally normal-shaped. In addition, $\langle {\sigma }_{\text{c},\text{i}}^{\text{+}}\rangle$ = $\langle {\sigma }_{\text{c},\text{i}}^{\text{-}}\rangle$ = 0.26 m s⁻¹. This result not only confirms that the estimate obtained from real data is robust, but also represents an improvement in the significance of the detection that is now increased to 3.7σ.

Results for K2-3d. The median of the K_d,i best-fit values is now K_d = 0.58${}_{\text{-0}.\text{15}}^{\text{+0}.\text{18}}$ m s⁻¹, while $\langle {\sigma }_{\text{d},\text{i}}^{\text{+}}\rangle$ = 0.40 m s⁻¹ and $\langle {\sigma }_{\text{d},\text{i}}^{\text{-}}\rangle$ = –0.34 m s⁻¹. This result is not very different from that presented in Sect. 5.1, and shows that, despite the 190 additional data to the real dataset, the Doppler signal of K2-3 d is still underestimated and not significantly detected.

Fig. 14 Example of a simulated RV dataset used to explore the effects on the characterization of the K2-3 planets of additional measurements taken over the 2017 season. The upper plot shows the complete mock dataset, while only the third seasonis shown in the second plot, to better appreciate the intensive simulated sampling.

6 Discussion and conclusions

This work was focussed on deriving the masses and bulk densities of the three planets transiting the nearby M dwarf K2-3, using 329 RV measurements collected with HARPS and HARPS-N over a period of 2.5 years. We found that stellar activity makes a significant contribution to the RV variations over the entire time period. We have also shown that, for the case of K2-3, this can be effectively mitigated using a GP regression with a quasi-periodic kernel. The results of our global model describe the stellar activity component in a plausible way, and this allowed us to derive a precise and accurate mass estimate for K2-3 b. We also derive a mass for K2-3 c with a significance of less than 3σ. However, using simulations, we demonstrate that our estimate is accurate. Conversely, we do not detect, in our data, the Doppler signal induced by the temperate planet K2-3 d.

Figure 15 shows a planetary mass-radius diagram that includes planets for which the mass and radius were both measured with a relative error better than 30%. Theoretical mass-radius curves for various chemical compositions (Zeng & Sasselov 2013; Zeng et al. 2016) are shown with solid lines. The precision of the radii of the K2-3 planets is mainly limited by that of the stellar radius (all the relative uncertainties are ~ 10%; see Table 1). One firm outcome of our analysis is that, for K2-3 b, an Earth-like composition (~ 33% of iron and ~67% of silicates) is rejected with high confidence; we note that in Fig. 15 the masses are represented on a logarithmic scale. Concerning K2-3 c, our mass determination excludes an Earth-like composition with a confidence level of ~ 4σ (assuming R_p = 1.77 R_⊕). The non-detection of K2-3 d was explored in detail through simulations showing that the real Doppler semi-amplitude K_d is likely less than 1 m s⁻¹ and its corresponding mass is $M_{\text{d}}={2.7}_{-0.8}^{+1.2}$ M_⊕ (Sect. 5). Looking at its position on the mass-radius diagram as derived from our simulations, the interior composition of K2-3 d would differ from that of the Earth with a confidence level greater than 2σ in mass and ~2σ in radius. We note that planets K2-3 c and K2-3 d occupy a region of the mass-radius diagram in which planet occurrence is rare, when only planets with mass and radius measured with a precision better than 30% are considered.

The corresponding bulk densities of all planets (ρ_p ~ 3 g cm^-3) show that they may have a very similar composition. If further measurements were to confirm our density estimates, excluding rocky compositions for K2-3 c and K2-3 d with higher significance, there are two scenarios that may explain the bulk properties of the K2-3 planets: water-poor planets with H/He envelopes or water worlds.

Several recent studies have investigated the water-poor hypothesis (Lopez 2017; Owen & Wu 2017; Jin & Mordasini 2018; Van Eylen et al. 2017). Fulton et al. (2017) analysed a sample of short-period planetary candidates detected by Kepler (P < 100 days) and demonstrated that the distribution of planetary radii is bi-modal: the planet candidates have radii that are predominantly either ~ 1.3 R_⊕ or ~ 2.4 R_⊕, and a gap is observed between 1.5 and 2 R_⊕. The evolutionary model of Owen & Wu (2017) reproduced the observed bi-modal radius distribution in terms of two populations of planets: those consisting of a bare core resulting from photoevaporation, and those with twice the core radius, where the size is doubled by a H/He envelopes. The “gap” detected by Fulton et al. (2017) is actually observed for planets orbiting FGK stars, while the M-dwarf regime was not explored. However, K2-3 is more similar to the FGK sample than to a late-M dwarf, and the results from Kepler could still be applicable¹⁵. Therefore, K2-3 b could have a significant volatile envelope, large enough to measurably change its radius with respect to that of a rocky core. The same may be the case for the other two planets in the K2-3 system, but neither of their radius estimates allow us to unambiguously associate them with one of the groups.

If all planets share the same composition, their densities can be explained by modest primordial hydrogen and helium envelopes atop Earth-like iron and silicate cores. Using the stellar and planet properties derived in this work and assuming a rocky Earth-like core and a solar composition H/He envelope, we find that K2-3 b, c, d are best fit with H/He envelopes comprising 0.7%, 0.3%, and 0.4% of their masses, respectively (Lopez & Fortney 2014). Moreover, using the planetary evolution models of Lopez (2017) we find that none of the planets in this system are vulnerable to losing significant mass through photo-evaporative atmospheric escape.

Even though the water-poor scenario has received a great deal of attention, alternatives should also be considered. Water worlds are planets having massive water envelopes comprising ≥ 50% of the planet total mass. Recently, bi-modal radius distributions have been derived for the complete Kepler sample of Q1–Q17 small exoplanet candidates with radii R_p < 4 R_⊕ (Zeng et al. 2017a,b). These distributions show that the limits and extent of the radius gap depends on the spectral type of the host star. One proposed explanation for the observed bi-modal distributions is the existence of two populations of planets: rocky worlds, with the lowest radii, and water worlds. The two populations likely share the same underlying rocky component by mass, but differ in the presence of a H₂ O-dominated mantle, which is similar to, or slightly more massive than, the rocky component. According to Zeng et al. (2017b), the radius gap for an M0 dwarf such as K2-3 is located at 1.6–1.7 R_⊕. Therefore K2-3 b has an observed mass and radius consistent with that expected for a water world. This could also be the case for K2-3 c and K2-3 d. However, the accuracy of their radii places these planets close to the transition limit that defines the gap, making their water-world membership not highly significant.

We therefore conclude that all three planets in the K2-3 system are likely sub-Neptunes, defined as small, non-rocky planets that have enough volatiles to change their bulk composition measurably. Both the H/He gas envelope and water-world scenarios are possible, particularly for K2-3 b. Based on our data, however, we cannot rule out that K2-3 c and K2-3 d have bare cores of purely rocky, Earth-like composition.

Within the H/He envelope scenario, planets likely formed within the first ~10 Myr before the dispersion of the gaseous proto-planetary disk. These water-poor planets could have formed in situ, and their similar masses would suggest similar formation histories (Lee & Chiang 2015, 2016). On the other hand, following Ginzburg et al. (2016) one may expect that the planets, during the cooling phase that follows their formation beyond the snow line (e.g. Selsis et al. 2007), were characterized by an intrinsic luminosity that could blow off, over a billion year timescale, any H/He envelope less than about ~ 5% by mass. Since a H₂O-layer has a much higher heat capacity than a H/He envelope, this evolutionary pathway could result in water worlds.

If K2-3 d is surrounded by a gaseous envelope, with the properties estimated here, this would likely result in surface pressure and temperatures that are too high to support a habitable planet scenario. A better characterization of K2-3 c and K2-3 d is left to the next generation of high-precision, high-stability spectrographs and to new photometric transit observations. We have shown that an intensive observing sampling over one season with two of the best spectrographs now available would still not have detected a signal with K = 1 m s⁻¹ and a period P = 44.5 days. Assuming that our GP result is a good representation of the stellar activity contribution, this means that the true mass of K2-3 d is expected to remain unmeasurable even with a dataset of more than 500 RVs, which is currently not possible for a single target and without a collaboration among various teams. Thus, detecting the real signal induced by K2-3 d is currently very challenging. K2-3 is, however, an ideal target for characterization studies with the VLT/ESPRESSO spectrograph (Pepe et al. 2014), or with near-infrared (NIR) spectrographs such as CARMENES (Quirrenbach et al. 2016), SPIRou (Artigau et al. 2014), and HPF (Mahadevan et al. 2014), provided that they reach their design RV precision and assuming, as expected, that RVs extracted from NIR spectra are less affected by stellar activity.

Fig. 15 Mass-radius diagram for exoplanets for which the mass and radius have been both measured with a relative error better than 30%. The location of the K2-3 planets is emphasized. For K2-3 d we also plot the mass derived from the GP analysis (shaded point in violet), and the value corrected using the result of the simulations described in Sect. 5.1. The curve for bulk density ρ = 3 g cm-3 is shown in grey passing through the three positions occupied by the K2-3 planets. The planetary data are taken from the NASA exoplanet archive and updated to August 30, 2017.

Acknowledgements

MD acknowledges funding from INAF through the Progetti Premiali funding scheme of the Italian Ministry of Education, University, and Research. We used high performance computing resources made available by INAF through the pilot programme CHIPP (through the proposal Precise planetary mass determination in radial velocity data collected with the HARPS and HARPS-N spectrographs: facing the challenges posed by the time sampling and the presence of stellar noise). We especially thank F. Vitello (INAF-OACt) for his assistance within the CHIPP programme. This research has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant Agreement No. 313014 (ETAEARTH). Parts of this work have been supported by NASA under grants No. NNX15AC90G and NNX17AB59G issued through the Exoplanets Research Program. AVa acknowledges support from the California Institute of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. PF acknowledges support by Fundação para a Ciência e a Tecnologia (FCT) through Investigador FCT contract of reference IF/01037/2013/CP1191/CT0001, and POPH/FSE (EC) by FEDER funding through the programme “Programa Operacional de Factores de Competitividade - COMPETE”. NCS acknowledges support 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.

References

All Tables

All Figures

Fig. 1 Upper plot: time series of the K2-3 RVs extracted with the pipeline TERRA using the HARPS-N spectra. Middle plot: GLS periodogram of the RV time series (red line) is shown. Three levels of p-values are indicated by the horizontal dashed lines. The green dashed line overplotted on the RV periodogram represents the window function of the measurements (bottom plot) shifted in frequency to be superimposed on the strongest peak of the RV periodogram. This helps to identify the alias frequencies of the most relevant peaks. The stellar rotation and planetary orbital frequencies are indicated by vertical lines and corresponding labels.

In the text

	Fig. 2 As in Fig. 1 but for TERRA RVs extracted from the HARPS spectra. Pre- and post-upgrade RV offsets have been applied, as derived from our analysis described in Sect. 4.
In the text

	Fig. 3 As in Fig. 1 but for all TERRA HARPS-N and HARPS RVs. We applied an offset to each separate dataset, as derived from our analysis described in Sect. 4.
In the text

	Fig. 4 Top plot: K2 light curve of K2-3 with the planetary transits removed. Bottom plot: GLS periodogram of the binned light curve (one point per day), showing a peak at Prot = 38.3 days.
In the text

Fig. 5 Top plot: light curve of K2-3 from the MEarth-south survey, folded at the best period P = 31.8 days found using GLS. Middle plot: GLS periodogram of the MEarth light curve. The green line corresponds to the window function of the measurements (bottom plot), which is shifted in frequency so that the peak is superimposed on the major peak of the RV periodogram to directly identify alias frequencies.

In the text

Fig. 6 Top plot: time series of the activity indicator based on the Hα line extracted from the HARPS and HARPS-N spectra. Second and third plots: GLS periodogram of the dataset is shown. The dashed horizontal lines indicate the p-value levels as derived from a bootstrap analysis. The greencurve corresponds to the window function of the measurements, which is shifted in frequency so that the peak is superimposed on the major peak of the RV periodogram (second plot), to directly identify alias frequencies. In the third plot the window function is shifted to be superimposed on the P ~ 40 days rotational signal. Bottom: time series phase-folded at the period P = 211 days. Red points represent the average of the data within 20 bins in the phase range [0,1].

In the text

Fig. 7 Top plot: time series of the activity indicator based on the Hα line extracted from the HARPS and HARPS-N spectra. Here, only the dataset of the third season is shown. Middle plot: GLS periodogram of the dataset is shown. The green line corresponds to the window function of the measurements, which is shifted in frequency so that the peak superimposes to the major peak of the RV periodogram, to directly identify possible alias frequencies. Bottom plot: time series phase-folded at the period P = 43.5 days.

In the text

	Fig. 8 First row: TERRA RV residuals, after removing our best-fit stellar component, phase-folded to the three planetary solutions found for the quasi-periodic GP model (represented by a red curve). Second row: histograms of the number of data divided in bins of phase are represented, showing a fairly uniform coverage for each planet.
In the text

	Fig. 9 Stellar signal contribution to the TERRA RV times series, for each observing season, as fitted with our GP quasi-periodic model (Table 3). The blue line indicates the best-fit curve, while the shaded grey area represents the ± 1σ confidence interval.
In the text

	Fig. 10 TERRA RV residuals, after removing the planetary and stellar activity signals, as fitted within a GP quasi-periodic framework (Table 3).
In the text

	Fig. 11 Distribution of N = 156 median values for the Doppler semi-amplitude Kc of planet K2-3 c, normalized to the injected value Kc, inj, as derived from posterior distributions of the N mock RV datasets. This result refers to the case of real observing epochs.
In the text

Fig. 12 Sample of posterior distributions (grey histograms) for the Doppler semi-amplitude Kd of planet K2-3 d obtained from the GP regression analysis of the mock RV datasets described in Sect. 5.1. Each plot shows the marginal distribution for a single mock dataset, and this is compared to the posterior distribution obtained forthe real TERRA RV dataset, represented by the yellow histogram. The vertical dashed lines indicate the 50th (green) and 68.3th (black) percentiles for the posterior distributions of the mock datasets.

In the text

	Fig. 13 Distributions of the 50th percentiles for all the posterior distributions of the semi-amplitude Kd obtained from 156 mock datasets, normalized to the injected value Kd, inj. This result refers to the case of real observing epochs.
In the text

	Fig. 14 Example of a simulated RV dataset used to explore the effects on the characterization of the K2-3 planets of additional measurements taken over the 2017 season. The upper plot shows the complete mock dataset, while only the third seasonis shown in the second plot, to better appreciate the intensive simulated sampling.
In the text

Fig. 15 Mass-radius diagram for exoplanets for which the mass and radius have been both measured with a relative error better than 30%. The location of the K2-3 planets is emphasized. For K2-3 d we also plot the mass derived from the GP analysis (shaded point in violet), and the value corrected using the result of the simulations described in Sect. 5.1. The curve for bulk density ρ = 3 g cm-3 is shown in grey passing through the three positions occupied by the K2-3 planets. The planetary data are taken from the NASA exoplanet archive and updated to August 30, 2017.

In the text