© ESO, 2017

1. Introduction

M dwarfs lie at the bottom of the main sequence, and consequently are small, cool, and intrinsically faint. These physical characteristics provide many advantages when looking for smaller, cooler, and fainter planets. For orientation, an Earth-mass planet in the habitable zone (HZ) of a 0.2 M_⊙ M dwarf, produces a radial velocity (RV) wobble that is over an order of magnitude larger than that of the Earth generates in the Sun. When caught in transit, this Earth-size planet decreases the flux of the M dwarf (with a 0.2 R_⊙ radius) by 25 times as much as the Earth does when it crosses the Sun. This makes planets around M dwarfs easier to detect, and even more importantly easier to characterize. M dwarfs with transiting planets – like GJ 436 (Butler et al. 2004; Gillon et al. 2007), GJ 1214 (Charbonneau et al. 2009), GJ 3470 (Bonfils et al. 2012), GJ 1132 (Berta-Thompson et al. 2015), and K2-3 (Crossfield et al. 2015; Almenara et al. 2015) – are therefore very popular targets in the exoplanet characterization community. Our HARPS program has made significant contributions to the detection, confirmation, and/or characterization of all of them. Some of the more numerous HARPS detections that do not transit may become amenable to characterization once future instruments are able to resolve them angularly, for example the Mid-infrared E-ELT Imager and Spectrograph (METIS) at the E-ELT(Snellen et al. 2015).

In this paper, we report a total of twelve planet detections in five planetary systems. The five stars, GJ 3138, GJ 3323, GJ 273, GJ 628, and GJ 3293, were observed with the HARPS spectrograph as part of our long-term effort to search for planets around bright, nearby M dwarfs. Section 2 summarizes the stellar properties of each host star, while Sect. 3 briefly describes our data set and how both radial velocities and spectroscopic indices are measured from the HARPS spectra. In Sect. 4 we identify periodic signals and ascribe them to either planets or activity-induced variations. Finally, Sect. 5 presents our conclusions.

2. Stellar properties

GJ 3138, GJ 3323, GJ 273, GJ 628, and GJ 3293 are part of our HARPS sample of M dwarfs because their distance is d< 11 pc, their V band magnitude V< 14 mag, and their rotational velocity v sini< 6.5 km s^-1 (Bonfils et al. 2013). We retrieved their spectral types, photometry, coordinates, proper motions, parallaxes, and multiplicity from the literature; the full references are listed in Table 1. The secular acceleration (dv/dr) is calculated following Eq. (2) from Zechmeister et al. (2009). We computed their luminosities using the Leggett et al. (2001) bolometric corrections. We estimated their physical parameters using calibrations by Delfosse et al. (2000) for stellar mass, Boyajian et al. (2012) for stellar radius and effective temperature, and Neves et al. (2013) for metallicity. We computed their UVW galactic velocities using the Johnson & Soderblom (1987) conventions and assigned them to dynamical population following the Leggett (1992) prescription. We computed the limit of the habitable zone of each star with the Kopparapu et al. (2013, 2014) online calculator¹, and adopted the new inner limit of the HZ from Kopparapu et al. (2016). This new inner limit incorporates 3D global climate models (GCM) in tidally locked exoplanets around cool stars. We quantified their magnetic activity by measuring their Mount Wilson S index (Wilson 1968; Vaughan et al. 1978) on their average HARPS spectra. We converted the S index into ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$ values following the procedure described in Astudillo-Defru et al. (2017). With $\log \mathrm{\left(}{\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}\mathrm{\right)}$ between −4.78 and −5.60, all five stars are magnetically quiet or, at most, very moderately active. While GJ 628 has been catalogued as a BY-Dra variable (Nakajima & Morino 2012), its photometric variability is only moderate since the dispersion of its V band flux is just 16 ± 7 mmag (Hosey et al. 2015). As discussed below, we find its variations consistent with a magnetic cycle. Table 1 summarizes these stellar properties.

3. Data

We obtained spectra of the five targets with HARPS (Mayor et al. 2003), a fiber-fed cross-dispersed echelle spectrograph at the 3.6 m telescope at ESO/La Silla Observatory (Chile). The instrument has a resolving power of R ~ 115 000 and almost full wavelength coverage between 380 nm and 690 nm. It is enclosed in a vacuum vessel and temperature-stabilized so as to reach a long-term radial-velocity precision <1 m s^-1. Calibrations are performed daily with a ThAr lamp (Lovis & Pepe 2007) and, since 2010, additionally with a white-light-illuminated Fabry-Pérot étalon (Wildi et al. 2010). Using both of these calibrations, the sub-m/s drift of the instrument can be measured by illuminating a second fiber with a calibration lamp, while the first fiber collects the light of the star. The offset between the simultaneous and daily calibrations can be subtracted from the radial velocity measurement of the star to improve its precision to a few tens of cm/s (see Sect. 2.6 in Fischer et al. 2016). We did not initially rely on this simultaneous calibration for our program, because some very strong ThAr lines spill from the calibration fiber into the stellar spectra to degrade our measurements of some important spectroscopic diagnostics in the blue-part of M dwarf spectra. The Fabry-Pérot calibration source has much more uniform line intensities, and therefore does not significantly contaminate the stellar spectrum. Starting in 2012 we thus changed our observing strategy and now use simultaneous calibration with the Fabry-Pérot source for stars brighter than V = 11 mag. For fainter targets, we use the second fiber to monitor the sky background.

Octagonal fibers are used in the entrance of a spectrograph to stabilize the illumination and to decrease the effect of imperfect guiding or seeing variations (Perruchot et al. 2011). The HARPS fiber link was upgraded with such a fiber on May 28, 2015 (Lo Curto et al. 2015), stabilizing its line-spread function but introducing a zero point in our measurement series.

3.1. Radial velocities

The HARPS Data Reduction Software (DRS) provides a quasi-real-time estimate of the radial velocity through cross-correlation (CCF) of the stellar spectrum with a software mask (Baranne et al. 1979; Pepe et al. 2002). This measurement, although a very good first estimate, uses a fixed template across M subtypes and discards the information contained in the many blended and weak lines were rejected at its construction. It therefore has a suboptimal signal-to-noise ratio, and we recompute improved velocities off-line through a few steps (Astudillo-Defru et al. 2015). We first align all spectra of a given star to a common frame using RV_CCF as a first guess, and co-add them into a high signal-to-noise stellar template. We then produce a high signal-to-noise template of the telluric absorption, subtracting the average stellar template contribution from each observed spectrum and co-adding the residuals in the laboratory frame. The stellar spectrum moves on the CCD according to the projection of the radial velocity of the Earth toward the star, but the telluric lines are fixed on the detector. The stellar and telluric spectra can thus be disentangled. Armed with these two templates, we measure a radial velocity from each observed spectrum as the Doppler shift that maximizes the likelihood between this spectrum and the high signal-to-noise template, using the telluric template to mask out the pixels most affected by telluric lines. Extracting the RV through maximum-likelihood optimization against a high signal-to-noise template was proposed long ago (e.g., Howarth et al. 1997; Zucker & Mazeh 2006), and independent implementations were applied to HARPS data by Chelli (2000), Galland et al. (2005, in Fourier space) and by Anglada-Escudé & Butler (2012).

Compared to the online pipeline, this reprocessing reduces the error bars by typically 20 to 30%. Since the HARPS upgrade with an introduction of a section of octagonal fiber altered the line spread function, we opted to process the pre- and post-upgrade epochs independently. For each star we thus produce pre- and post-upgrade templates, and we adjust a zero point offset when analyzing the time series.

3.2. Activity proxies

Because stellar activity can confuse planet searches, we need to monitor its strength and its variation on timescales that range from the stellar rotation period (a few to 100 days) to the length of a magnetic cycle (a few hundred days to a dozen years). We therefore monitor the shape of the spectral lines and measure spectroscopic indices in our HARPS spectra, and we also use publicly available photometry.

3.2.1. Spectroscopic indices

The HARPS DRS provides the full width at half maximum (FWHM) and the bisector span of the cross-correlation function (Queloz et al. 2001), which are both useful probes of stellar activity. We also measure from the spectra the chromospheric emission in the Ca ii H&K lines, and in the Hα, Hβ, Hγ, and Na D lines (Wilson 1968; Giampapa et al. 1978; Linsky et al. 1982). Our measurements of the Ca ii H&K follow Vaughan et al. (1978, S and those of HαGomes da Silva et al. (2011). For the other lines that are sensitive to stellar activity, we compute the ratio of the integrated flux in the line (C) over that in two control bands (V, R): $\mathrm{H}\mathit{\beta }\mathrm{=}\frac{{\mathit{C}}_{\mathit{\beta }}}{{\mathit{V}}_{\mathit{\beta }}\mathrm{+}{\mathit{R}}_{\mathit{\beta }}}\mathrm{;H}\mathit{\gamma }\mathrm{=}\frac{{\mathit{C}}_{\mathit{\gamma }}}{{\mathit{V}}_{\mathit{\gamma }}\mathrm{+}{\mathit{R}}_{\mathit{\gamma }}}\mathrm{;}\mathrm{Na}\mathrm{D}\mathrm{=}\frac{\mathit{C}{\mathrm{1}}_{\mathrm{Na}}\mathrm{+}\mathit{C}{\mathrm{2}}_{\mathrm{Na}}}{{\mathit{V}}_{\mathrm{Na}}\mathrm{+}{\mathit{R}}_{\mathrm{Na}}}\mathit{,}$(1)where the wavelength passbands are listed in Table 2.

3.2.2. Photometry

Since the contrast between spots and the rest of the photosphere increases at bluer wavelength, we use photometry in the optical rather than near-infrared band, and retrieved photometric time series of the five targets from the All Sky Automated Survey (ASAS², Pojmanski 1997). Since M dwarfs are relatively faint objects in this band, the ASAS photometry through the smallest syntheticapertures provide the highest quality measurements by minimizingthe sky background.

We systematically searched for spot induced periodicities for the whole photometric time series and for each seasonal observations. The seasonal analysis helps in the detection of the stellar rotation period assuming that active regions evolve and last for an unknown period of time.

4. Data analysis

4.1. Methods

We applied the same method for the five stars, searching the time series for variability and periodicities, and modeling the periodic signals as Keplerian orbits. In an effort to disentangle activity signals from planets, we compare the detected periods with the stellar rotation period and we search all RV signals for a counterpart in the time series of the activity indicators.

Using the CCF computed by the HARPS DRS, we start by building a representative measurement uncertainty for each star that accounts for the photon noise (Bouchy et al. 2001) and the 0.60 ms^-1 calibration error (see Fig. 6 in Fischer et al. 2016), but not for the stellar jitter. We quadratically add the two noise components for each epoch, and compute the typical expected error σ_i as a weighted average of these epoch uncertainties. Prior to the RV analysis, we subtract the secular acceleration (see Table 1). The observed dispersions of the pipeline RV measurements of GJ 3138, GJ 3323, GJ 273, GJ 628, and GJ 3293 stand above this expected error by factors of respectively 1.5, 1.5, 3.2, 3.7, and 3.1.

This motivated the extraction of the offline RVs described above, which we searched for periodicities by computing a generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009). Our GLS normalization choice is such that a GLS power of 1.0 means that a sine wave is a perfect description of the data, while a power of 0.0 means that a constant model is an equally good fit. To translate GLS power values to false alarm probabilities (FAP), we use bootstrap Monte Carlo simulations, creating 10 000 synthetic time series by shuffling the original velocity values between the observation dates. We compute a periodogram for each bootstrapped data set, and record the distribution of the maximum value of the bootstrapped periodograms. We obtain the FAP for a given GLS power as the fraction of the bootstrapped periodograms that had a higher maximum value. The window function of the data sets is analyzed in order to identify aliasing or periodicities due to our temporal sampling.

When we identify a periodic signal with a FAP under 1%, we adjust a Keplerian to the velocities and look for additional periodic signals in the fit residuals, iterating until no significant signal is detected.

To protect against overfitting the data with a more complex model than they warrant, we compute the Bayesian information criterion (BIC: minus twice the Schwarz criterion, Schwarz 1978) for models with increasing numbers of Keplerian components. The BIC is widely used for model comparison because it is easily computed³; however, the BIC’s reliability as an approximation of the logarithm of the (more rigorous but harder to derive) Bayes factor depends on the sample size (Kass & Raftery 1995) and quantifies the evidence in favor (or not) of a model over the null hypothesis model H₀ (generally the simplest model).

We adjusted the (multi-)Keplerian models using the YORBIT (Ségransan et al. 2011) code, which efficiently explores the model parameter space through a genetic algorithm. After model convergence, we estimated the parameter uncertainties with a Markov chain Monte Carlo algorithm. As noted above, all models include a free zero-point offset between the velocities acquired before and after the HARPS fibers upgrade (Lo Curto et al. 2015).

4.2. GJ 3138

4.2.1. Periodicity analysis

We acquired 199 radial velocities of GJ 3138, of which 156 were obtained prior to the mid-2015 HARPS fibers upgrade and 43 afterwards. We discarded one measurement (epoch BJD = 2 457 258.81) which is a strong outlier. The root mean square dispersion of these velocities is σ_{(O − C)} = 2.87 ms^-1, well above their expected uncertainty of σ_i = 1.82 ms^-1.

Our iterative periodogram analysis detects four significant periodicities (Fig. 1). The periodogram of the radial velocities shows a strong peak around a period of 6.0 days, with a power far higher than any found in the 10 000 reshuffled data sets. Its FAP is thus well below 0.01%. The residuals of a fit of the corresponding Keplerian model (model k1) have periodogram peaks at periods of 48.2, 1.2, and 257 days, with FAPs of respectively 0.46%, 0.49%, and 4.72%. In the next section, we show that the peak at 48.2 days is likely due to stellar activity since after subtraction of the other periodic radial velocity signals we observe a strong correlation with the Hα flux. Instead of modeling this 48.2-day periodicity as a Keplerian signal, we use the RV-Hα correlation derived in the next section to describe it. After subtraction of both that activity signal and the first planet, the strongest periodogram peak is at 1.2 days and has a 0.53% FAP. We adjust a two-Keplerian model (model k2) to the activity-corrected radial velocity time series (a RV time series now corrected from the RV-Hα correlation). In a fourth iteration, the periodogram of the residuals exhibits a significant power excess at a period of 257 days, with a FAP = 0.29%. A fifth iteration finds two peaks at 1.05 and 20.4 days, which are 1-day^-1 aliases of each other. The strongest of the two has a 0.9% FAP that it just at our detection threshold of 1% FAP . However, the BIC values of the models with three or four Keplerians (ΔBIC_K3,K4 = 4.4) do not favor this more complex model. In addition, a periodicity of 20.4 days is half the rotation of GJ 3138. We therefore interpret that this low RV signature is possibly provoked by stellar activity.

4.2.2. Validating the planet interpretation

The 48-day signal as stellar rotation.

The stellar rotation period estimated from the average ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$ is 34 days (Table 1). This is safely remote from the 1.2-day, 6.0-day, and 257-day periods, but within its likely error bar of the 48-day period. We therefore search the time series of the activity indicators for periodicities (Fig. 2) and for correlation with the measured radial velocity (Fig. 3, bottom). Since the fiber upgrade introduced an offset in both the FWHM_CCF and BIS_CCF time series, we restrict this analysis to the pre-upgrade measurements as they represent the bulk of the data set. The ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$, Hα, and FWHM_CCF periodograms show significant peaks at periods of respectively 41.8, 44.6, and 41.2 days, while neither the periodogram of the bisector nor that of the ASAS photometry (497 points over 9.0 yr) show significant periodicities. The periodograms for the seasonal photometry do not show evidence of stellar rotation. The three identified signals have similar periods, which are also close to the 34-day stellar rotation period estimate from ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$ and to the 48.1-day periodicity in the RV time series. We surmise that these periods all reflect modulation of photospheric inhomogeneities by the stellar rotation; their dispersion is due to some combination of differential rotation and measurement noise.

We turn to subsamples of the time series to confirm this suspicion. We subtract the 1.2- and 6.0-day signals from the RV series to focus on the 48-day signal, and divide the activity diagnostic and RV time series into three BJD−2 450 000 intervals: 6500−6700, 6700−7200, and 7200−7500, with 57, 77, and 41 measurements. Figure 3 shows the corresponding periodograms. The periodograms of the RV residuals and Hα for the second epoch ([6700−7200]) are strikingly similar; both have significant power excess at ~43 days. The RV residuals clearly anti-correlate with Hα for that epoch (Fig. 3, bottom), with a Pearson coefficient of − 0.42. The same figure also shows Hα phased to a 48.2-day period, with points for phases between 0.2–0.5 mostly in the upper left corner and points for phases between 0.7–0.9 mostly in the lower right part of the plot. We fit the anti-correlation and obtain $\mathrm{RV}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1.953}\mathrm{\pm }\mathrm{0.262}\mathrm{)}\mathrm{\times }\mathrm{H}\mathit{\alpha }\mathrm{+}\mathrm{(}\mathrm{0.106}\mathrm{\pm }\mathrm{0.001}\mathrm{)}\mathit{.}$(2)We use this relation to correct the activity effect over the central interval. Without evidence for a similar relation over the other two intervals, we do not correct the rest of the RV time series. This partial correction reduces the strength of the 48-day signal very significantly, as demonstrated by a comparison of the second and third periodograms in Fig. 1. This adds to the evidence for stellar activity, and we classify the 48-day periodicity as an activity signal rather than a planet.

We inquire if the RV signals for the two shorter periodicities (1.22 and 6 days) are present along the entire data set – favoring the scenario where RV variations are due to the presence of an exoplanet – or concentrate in a particular subsample – challenging the planetary interpretation. We thus divide in three our RV time series, with subsamples that satisfy BJD-2 450 000 in the ranges 6500–6700, 6700–7200, and 7200–7500 and compute their periodograms. Figure 4 shows those periodograms specifically around 1.22 and 6 days, where both the 1.22- and 6-day signals seen as power excess in most epochs.

In parallel, we need to evaluate if the power excess are actually expected given the sampling of each one of the three subsamples. To that intend, we make additional synthetic time series (100) by drawing random RV points with a normal distribution centered around the best Keplerian fit and with a 2 m s^-1 standard deviation. For each one of the synthetic time series we make a new periodogram, and the gray areas in Fig. 4 represent their ±1 sigma power distribution.

Around both 1.22 and 6 days, the periodogram of the original time series (black lines in Fig. 4) is never found below the ±1 sigma region (gray areas). This means both signals are indeed detected every time the sampling allows so, which in turn means the signals seem coherent rather than transitory.

Overall, we conclude on the planetary nature of both the 1.22- and 6-day signals.

4.2.3. Orbital parameters

Considering the above analysis, we detect 4 periodicities and conclude to the detection of 3 planets and the stellar rotation. Our Keplerian solution is illustrated in Fig. 5 and presented in Table 3, which reports separate noise estimates and rms residuals around the solution for the data acquired before and after the HARPS fibers upgrade. For the stellar mass reported in Table 1, the semi-amplitudes convert to minimum masses of 1.9, 4.4 and 10.9 M_⊕, in the regime of super-Earths and mini-Neptunes. Their irradiance are 116, 14, and 0.09 times the terrestrial one, which for extreme Bond albedos of 0.75 and 0.0 correspond to equilibrium temperatures ranges of respectively 640–900 K, 375–530 K, 107–150 K.

4.3. GJ 3323

4.3.1. Periodicity analysis

We acquired 157 radial velocities on GJ 3323, including 132 before and 25 after the HARPS upgrade. We rejected three points because two have very low SNR (<3 at 550 nm; BJD = 2 456 630.65, 2 456 926.86) and one is a >5σ outlier (BJD = 2 457 258.92). The 154 remaining RVs show a standard deviation σ_{(O − C)} = 2.85 m s^-1, in excess of estimated uncertainties σ_i = 2.18 m s^-1. The periodogram of raw RVs (Fig. 6) shows a powerful peak at around 5.4 days with FAP ≪ 0.01%. We fit a Keplerian to that first signal and compute the periodogram of the residuals. We found a second periodicity at ~40 days with, again, a FAP ≪ 0.01%. We made a third iteration and found no more significant periodicity.

4.3.2. Challenging the planet interpretation

GJ 3323 b.

GJ 3323 b is found with a period of 5.4 days, which is significantly shorter than the 33-day period estimated for the stellar rotation (Table 1). We know that in some configurations, a spot coming in and out of view could produce a Doppler signal with a significant power excess at one-half or one-third of the stellar rotation period (Boisse et al. 2011). It is thus reassuring to know the 5.4-day period is also significantly different than any harmonic of the rotation.

As we did for the GJ 3138 RV analysis, we split the RV time series into three time intervals. They are chosen such that BJD-2 450 000 is in the ranges 6500–6800, 6800–7200, and 7200–7500, with 56, 67, and 24 points, respectively. Figure 8 shows the periodograms with, on the left-hand side, a focus on a period of around 5 days. The signal is recovered at all epochs.

We can confidently accept that the first signal is a planet.

GJ 3323 (c).

The periodicity of around 40 days is close to our estimate of the stellar rotation. It therefore requires more attention to the different activity proxies.

Periodograms of ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$, FWHM_CCF, and Hα present significant power excess in the period range of ~160 days (Fig. 7). The FWHM_CCF also shows additional power excess around 40 days. Conversely, in the same figure, the BIS_CCF index does not show any significant periodicity and the periodogram for the whole ASAS photometry time series (654 points spanning 9.0 yr) shows power excess at about 665 (2% FAP) and 88 days (10% FAP). From the seasonal analysis, a periodicity close to 88 days is detected for the ASAS seasons between BJD 2800−3200 (0.24% FAP) and 3500−4000 (0.42% FAP). This value is compatible with the 88.5 days reported by Kiraga (2012) for the rotational period. However, the reliability of the 88.5-day stellar rotation is questionable because only 2/10 of ASAS seasons shows the ~90-day periodicity even though spots on M dwars are usually more stable over time (e.g., GJ 674, GJ 581, GJ 174, GJ 551) and the stellar rotation period estimated from the $\log \mathrm{\left(}{\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}\mathrm{\right)}$ is 33 days. Except for the smaller power excess at 40 days in FWHM_CCF, the variability is thus seen in a period range much different than that of the 33-day period estimated for the stellar rotation by the ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}\mathrm{-}{\mathit{P}}_{\mathrm{rot}}$ relation (Table 1). It is possible that these longer timescale changes are related to magnetic cycles instead (Gomes da Silva et al. 2011).

As previously, we subdivide the time series of both RVs and activity proxies into the same three time intervals. We subtract the best Keplerian fit to RVs to focus on the 40-day periodicity (right-hand side of Fig. 8). We note that only the first two intervals have a meaningful sampling able to probe the 40-day signal. The third interval only has 24 points spread on just one period and a half. A power excess appears at about 40 days for both the RVs and FWHM_CCF periodograms in the 6500−6800 and 6800−7200 subsamples.

In Fig. 9, we further inspect RVs against FWHM_CCF for two time intervals and it is hard to exclude possible resemblances, in particular if we allow for a shift in phase. However, no other activity indicator shows evidence of stellar activity at about 40 days (Fig. 7) and currently we have no example of FWHM_CCF variation without Ca ii H&K or Hα variation; additionally, the FWHM_CCF periodograms for the 6500−6800 and 6800−7200 subsamples show most of the power centered on periods 80−250. We therefore search for other possible origins of the power excess seen at about 40 days in the FWHM_CCF periodogram.

We show in Fig. 10 the time series for the FWHM_CCF subsample satisfying BJD-2 450 000 ϵ [6800−7200]. It is evident that this activity tracer is higher at about BJD-2 450 000 = 7050. In the same figure, we note that our sampling generates a peak at ~60 days in the window function. This means that the power excess at about 40 days is an alias of the strongest peak located at about 100 days (1/37.5 = 1/100+1/60); when subtracting a 100-day periodicity, power excess is no longer present at about 40 days. However, this peak at ~100 days and its ~40-day alias remains unexplained because the other activity indicators show a power excess concentrated at ~150 days for the ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$ and Hα, and at 665 days for the photometry (Fig. 7), while a stellar rotation period is expected to be about 33 days (Table 1). Taking these data all together, and before more data could be acquired, we think GJ 3323c should be considered a planet candidate.

4.3.3. Keplerian analysis

We use YORBIT to fit a model with two Keplerians with semi-amplitudes of 2.55 ± 0.32 and 1.49 ± 0.32 m s^-1. We estimated a stellar mass of 0.164 M_⊙ (Table 1). This converts the semi-amplitudes to 2.02 ± 0.25 and 2.31 ± 0.50M_⊕ planets and makes GJ 3323 one of the lowest-mass stars (0.16 M_⊙) known with planets discovered by radial velocity. GJ 3323b and GJ 3323(c) receive 2.51 and 0.17 times the Earth irradiance and have an equilibrium temperature in the range of 420–595 K and 214–303 K, respectively, if bond albedos between 0 and 0.75 are assumed. Details of the derived parameters are given in Table 4, and Fig. 11 shows the phase folded RVs.

4.4. GJ 273

4.4.1. Periodicity analysis

We intensively monitored GJ 273, recording 280 spectra between December 2, 2003, and September 30, 2016. From this 12.8 yr of RV monitoring, 43 of them were acquired after the HARPS fibers upgrade. The RVs have a standard deviation of σ_{(O − C)} = 2.75 m s^-1, a value which is well above the estimated uncertainties of σ_i = 0.94 m s^-1.

Figure 12 shows the periodograms, first for the raw time series, and then for the residuals of each iteration. Most of the power in the raw RV periodogram is concentrated at low frequencies with FAP ≪ 0.01, and YORBIT converges to a solution with a periodicity of ~420 days. The periodogram of the residues shows a powerful peak largely above our 1% FAP threshold at about 20 days, modeled by a second Keplerian. After subtraction of this model, the two highest peaks in the periodogram of the residues are located at about 700 days and 5 days (both with FAP ≪ 0.01). When adding a third Keplerian to the model, the favored solution contains the two previous periodicities plus the one at 700 days. Now the 5-day peak and its 1-day alias dominate the periodogram of the residues. Although the residues still have power excesses at about 100 days (0.5% FAP), our final model consists of four Keplerian because this periodicity is explained well by a signal generated by the stellar rotation (99 days; Table 2). As we describe below, we suspect that both the two longest periodicities (with P ~ 420 and ~700 days) also have a stellar origin.

4.4.2. Challenging the planet interpretation

GJ 273 b, c.

The two RV signals with the shortest periodicities (4.7 and 18.6 days) have periods that are very different from that estimated for the stellar rotation (P_Rot ~ 99 days; see Table 1), very different from that measured for the two long-period RV signals (which we suspect are caused by stellar activity), and very different from the harmonics of those periodicities (~P/2, ~P/3, ...). They thus already appear as robust planet detections.

We wanted to check whether the signals are coherent across different epochs. We divided our RV time series into four parts, with subsamples satisfying BJD-2 450 000 within the ranges 6200−6600, 6600−6800, 6800−7200, and 7200−7500. We computed the periodograms for each subsample. As shown in Fig. 14, power excess is seen in almost all epochs.

As we did for GJ 3138, we needed to understand whether the signal would be detected at each epoch. We proceeded by making a 100 synthetic RV time series with the same sampling as our observation, normally distributed (σ = 2 m s^-1) and centered on the best Keplerian fit. The periodogram for each of the 100 time series was computed, from which we derived the ±1 power distribution. Because the periodogram of the original subsamples fell inside the ±1 region (gray zone in Fig. 16), we conclude that there is no evidence that either the 5- or the 18-day signal is transitory.

Overall, we accept both the 5- and 18-day signals as bona fide planet detections. Figure 16 shows the phase folded RVs.

Magnetic cycle and rotation.

In Fig. 12 we note that the raw RV periodogram exhibits power excess for periodicities from 300–4000 days, with peaks at 400–500 and 2000 days. Similarly, the periodograms of activity indicators show significant power excess in the same range of periodicities (Fig. 13). Time series for radial velocities and activity indicators are shown in Fig. 15. The correlation and p-value are 0.18 and 0.002 respectively. However, the maximum RV corrections from this correlation are of a few cm/s and are too low for a noticeable improvement.

Another way to correct the RV-activity dependency is to proceed as in Dumusque et al. (2011). They noticed that both RVs and Ca ii H&K time series show remarkable resemblances. Dumusque et al. (2011) fit a Keplerian to the ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$ and used the resulting parameters to fix a Keplerian model to RVs, but leaving the amplitude free. Proceeding this way we fix a model with P ~ 2000 days, then the periodogram of RV residues shows subsequently the signals from GJ 273 b and c. The periodogram from a third iteration exhibits a clear peak at ~100 days, the expected stellar rotation period (Table 1). Evidence of the stellar rotation appears in the periodograms of the ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$ and the Na D-index measured between BJD 6200–6500 and 6500–6800, respectively.

The ~420-day peak in the RV periodogram could be explained as the 1-yr alias of the longer period of ~2000 days seen in RVs and activity tracers (Figs. 12, 13, 15). The peak arising at ~700 days in the RV residues is close to the second harmonic (P/3) of the magnetic cycle. Because of that, we suspect these two RV signals (P ~ 420, 700 days) are the product of a stellar activity cycle that has not been modeled well.

4.4.3. Keplerian analysis

GJ 273b is a super-Earth with a minimum mass of 2.89 ± 0.26M_⊕ with an orbital period of 18.650 ± 0.006 days. Orbiting at a distance of 0.09110 ± 0.00002 AU from its parent star, it is grazing the inner edge of the conservative habitable zone (Kopparapu et al. 2013), but with an incident flux of 1.06S_⊕ it is well within the HZ if one assumes the planet is surrounded by an atmosphere and accounts for GCM (Kopparapu et al. 2016). It has an equilibrium temperature in the range 206–293 K (using albedos of 0.75 and 0, respectively). GJ 273c is among the less massive planets detected by radial velocities, with a minimum mass of 1.18 ± 0.16M_⊕. It completes one orbit in 4.7234 ± 0.0004 days. At the orbital distance of 0.036467 AU it receives almost seven times the stellar flux received at Earth and has an equilibrium temperature within 327–462 K, and thus much too close to the star to be considered in the circumstellar habitable zone.

4.5. GJ 628

Between June 2, 2004, and September 28, 2015, we collected 190 spectra of GJ 628, 148 of which are publicly available and were previously analyzed by Wright et al. (2016). Three of the measurements show a difference between GJ 628 systemic and the Moon RVs < 1 km s^-1 and the angle between the objects <25°. We considered these measures (BJD = 2 454 661.6, 2 454 662.6, and 2 454 663.7) as contaminated by the Moon and they were rejected after the subsequent analysis. As we show below, we came to the same conclusions for the RV signals at shorter periodicity (4.9 and 17.9 days), but a different conclusion for the periodicity at 67 days. Here we detect a third RV signal at a longer periodicity and identify a stellar activity cycle.

GJ 628 RVs have a dispersion of σ_{(O − C)} = 2.72 ms^-1 and a mean internal error of σ_i = 1.02 ms^-1. The dispersion excess comes from several periodic variations as is shown in Fig. 17. The periodogram of raw RVs exhibits two narrow spikes at 4.9 and 17.9 days, with FAP < 10^-2% and 0.46%, respectively, and a third one at about 205 days with a 0.29% FAP. After subtracting a two-Keplerian model from the 4.9- and 17.9-day periods, the periodogram of resulting residues (third row in Fig. 17) shows several peaks located at 217, 48, 454, 115 days, detected with FAPs < 0.01%, 0.02%, 0.07%, and 0.17%, respectively. We added a third Keplerian to each period with power excess (one at a time) and found that the solution including the 217-day signal results in significantly lower ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}$ and BIC (see Table 8). We note that the 67-day RV signature discussed in Wright et al. is not significant here, as the weak peak located at this periodicity has a 2.5% FAP. After including the 217-day signal in our model the highest peak in the periodogram of RVs residues (P ~ 380 days) has a FAP of 6.3%, and therefore no more signatures are detected.

Stellar activity and orbital solution.

The Ca ii H&K chromospheric emission reveals that GJ 628 presents an activity cycle (Fig. 18). We combined the literature Ca ii H&K measurements (Wright et al. 2004; Isaacson & Fischer 2010) to better understand the cycle identified with the HARPS data alone. Although the data time span is increased by about 4 yr, the combined data does not cover a complete activity cycle, which is at least 19 yr. Figure 18 also shows one flare event detected in both ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$ and Hα (BJD = 2 456 481.55). Concerning the FWHM_CCF, in addition to the expected offset due to the HARPS fiber upgrade, this activity index – computed through the DRS – presents other unexplained offsets of about 1 km s^-1 (too large to be of stellar origin) and its analysis is omitted here.

The stellar rotation period of GJ 628 is estimated to be 93 days (Table 1) and the Hα periodogram (Fig. 18) presents a significant power excess at 64 days; however, we cannot assert that this periodicity is linked to a yearly alias of the stellar rotation. We do not identify any periodicity linked to the stellar rotation either from the seasonal Hα or from the seasonal ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$ measurements. The peak at about 64 days in the Hα periodogram and the 67 days RV signal reported in Wright et al. (2016) are worryingly close. The HARPS only ${\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}$ periodogram presents a power excess at 120 days that can be explained as the yearly alias of the stellar rotation period. Considering the possible aliases, and unlike Wright et al. (2016), we do not regard the 67-day signal as a planet detection. The ASAS photometry, which consists of 600 points spread over 8.7 yr, shows hints of the stellar rotation. The periodogram for the photometry satisfying BJD between 2 452 650 and 2 452 950 shows a strong peak at about 100 days, and the entire data set shows two peaks near this periodicity: one peak with 0.95% FAP at 129 days and another with a FAP slightly greater at about 86 days.

Additionally, the third strongest peak of the window function data set is located at the 66-day periodicity (bottom panel in Fig. 17), hence we cannot exclude that the 67-day signal seen in RVs in Wright et al. and Hα are generated by the data temporal sampling itself. In the same figure we note the strong peak located at one year, showing that one-year aliases are highly expected.

Stellar differential rotation could be inferred from the APT photometry (spanning 6 yr) presented in Kane et al. (2017) plus the ASAS photometry. Minimum and maximum periods from the photometry are P_min = 84.2 ± 3.0 and P_min = 106.3 ± 7.3. Well-determined periodicities from the APT photometry shows no evidence of P_Rot. ≥ 120 days, and we interpret the 120-day and 129-day peaks in Figs. 18 and 19, respectively, as aliases of the true stellar rotation period.

4.5.1. Keplerian analysis

When adjusting a Keplerian to the 48-day periodicity (second highest peak in the RV_omck2 periodogram) the resulting orbit eccentricity is relatively high (e ~ 0.7). Then, the periodogram of the residues shows the 217-day peak with high significance. Regarding the stability of the orbits, a solution including this periodicity is uncomfortable because it crosses the orbits with 17.9- and 217-day periodicities, if coplanar orbits are assumed.

The more conservative solution, and that preferred by both minimum χ² and BIC, consists in three Keplerians adjusting the 4.9-, 17.9-, and 217-day periodicities. The orbital parameters for such a solution are listed in Table 6, and Fig. 20 depicts the phase folded RVs with the model. Gl 628c falls outside the conservative HZ (Kopparapu et al. 2013); however, when considering the GCM in synchronous rotating planets (Kopparapu et al. 2016), the incident radiation means that the planet lies within the circumstellar HZ of its parent star.

4.6. GJ 3293

Astudillo-Defru et al. (2015) reported the detection of two Neptune-like companions orbiting GJ 3293 with periodicities of 30.6 and 124 days, and a possible super-Earth orbiting with a periodicity of 48.1 days. To shed light on the nature of the 48.1-day signal, we gathered 61 new RVs from October 13, 2014, to April 4, 2015. In this work we measure a RV dispersion of σ_{(O − C)} = 7.73 ms^-1 and a mean internal error of σ_i = 2.33 ms^-1.

The Astudillo-Defru et al.Hα analysis suggested that the stellar rotation period of GJ 3293 is 41 days, and from RVs they argued about an uncertain nature of the 48.1-day signal because 48.1 days is close to the stellar rotation and because of the lack of such a signal in the 2010-2011 data set. By using several chromospheric activity indicators ($\mathrm{H}\mathit{\beta ,}\mathrm{H}\mathit{\gamma ,}\mathrm{Na}\mathrm{D}\mathit{,}{\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime }}\mathit{,}\mathrm{FWH}{\mathrm{M}}_{\mathrm{CCF}}$) and improving the Hα detection, we confirm that the stellar rotation is 41.0 days. The ASAS photometry (599 points spanning 8.7 yr) does not show periodic variability compatible with stellar rotation. However, the Moon’s synodic period is evident in the photometry periodogram. The synodic periodicity appears with 2.3% FAP for the photometry acquired between BJD 1500 and 2020. Neither a 41.0-day signal nor its yearly aliases at 36.9 and 46.2 days appears to affect the 48.1-day RV signature in the periodogram. Currently, the 41-day period of the stellar rotation and the 48.1-day period seen in RVs seem distinct enough to interpret the 48.1-day signal as being due to a planetary companion. Moreover, after subtracting both the 30.6- and 124-day signals and after analyzing the RV residues (RV_omck2), we found that the 48.1-day periodicity is temporally stable when splitting the sample in a more optimal way (target seasonal visibility instead of calendar year). The four RV subsets used for the present analysis satisfies BJD between 4790 and 5750, 5750 and 6100, 6100 and 6350, and 6900 and 7150. Figure 21 shows such a stability of the 48.1-day signal and the periodograms of a set of activity indicators. This planet is within the habitable zone of its parent star, with an equilibrium temperature between 171 K and 241 K.

When modeling the 30.6-, 48.1-, and 124-day RV periodicities by three Keplerians, a fourth signal appears in the residues. The periodogram (Fig. 22) shows power excess at a periodicity of 13.3 days. Such a peak was also present in Astudillo-Defru et al. (2015, see their Fig. 4), but without enough significance (1σ). Here the corresponding peak is clearly visible despite its weakness, and the FAP decreases to 3.4% (2σ). Furthermore, a model consisting of four Keplerians is significantly favored over the three-Keplerian model (see the BIC difference in Table 8). Accordingly, our solution for GJ 3293 consists of two super-Earths (P = 13.3, 48.1 days) and two Neptune-mass planets (P = 30.6, 122.6 days) with minimum masses of 3.3 M_⊕, 7.6 M_⊕ and 23.5 M_⊕, 21.1 M_⊕. At their separation from GJ 3293, the Neptune-mass planet b and the super-Earth planet b are within the circumstellar habitable zone. Figure 23 shows the phase folded RVs with the four Keplerians, and Table 7 gives the orbital parameters.

In comparison to Astudillo-Defru et al. (2015), where a model with three Keplerians was used, GJ 3293 RVs are modeled here with four Keplerians. This may involve changes in some of the orbital parameters. The only significant change (2σ) comes from the eccentricity of the orbit with a 122.7-day periodicity, where we found e = 0.11 ± 0.09 instead of the 0.37 ± 0.06 reported in Astudillo-Defru et al. The other orbital parameters remain unchanged.

We performed an additional test on our GJ 3293 solution by running an N-body integrator. We ran GENGA (Grimm & Stadel 2014) through the Data and Analysis Center for Exoplanets (DACE⁴). The initial conditions in Table 7 are used and our simulation lasts for 10⁶ yr. Figure 24 shows the stability of the system over this time span, as well as the secular variation of the eccentricity.

5. Summary and conclusions

We have analyzed between 150 and 280 RVs of GJ 3138, GJ 3323, GJ 273, GJ 628, and GJ 3293. To fully exploit the Doppler information contained in M dwarf spectra, the RVs extracted from the HARPS data was done by maximizing the likelihood. For each star, a high signal-to-noise reference spectrum was built out of all the observed spectra. The RV analysis includes a careful monitoring of several spectroscopic activity indicators and photometry, which we used to measure stellar rotation periods for three stars, except for GJ 3323. For GJ 628 Ca II H&K chromospheric emission reveals a long-term variation compatible to a Sun-like magnetic cycle.

Conditioned by a knowledge of the orbital inclination, our analysis results in the detection of 12 planets (9 new), including Earth-mass planets and super-Earths. The GJ 3138 radial velocities reveal four signals: three of them are compatible with planets having minimum masses of 1.8, 4.2, and 10.5 Earth masses, orbiting with periods of 1.22, 5.97, and 258 days, respectively, while a signal with a periodicity of 48 days is more related to stellar activity as it shows dependence on Hα. We recall that the RV-Hα dependence was subtracted and that more RVs of this star are needed because the phase for the 258-day signal is not fully covered. These orbits are located outside the habitable zone of GJ 3138.

We interpret the RV variations of the fully convective star GJ 3323 as the signal coming from two low-mass planets. These planets have minimum masses of 2.0 and 2.3 Earth masses M_⊕ and are orbiting with periodicities of 5.36 and 40.5 days, respectively. The planets are orbiting outside the conservative boundaries of the HZ. GJ 3323b and GJ 3323c respectively receive 2.5 and 0.17 times the stellar flux compared to the Sun.

The M dwarf GJ 273 is orbited by a super-Earth and by an Earth-mass planet. They have minimum masses of 2.89 and 1.18 M_⊕ and orbital periods of 18.65 and 4.72 days, respectively. GJ 273b receives 1.06 times as much radiation as our Earth does from our Sun. It thus lies well within the habitable zone (Kopparapu et al. 2016) and water (if any) may flow on its surface. Located at only 3.8 parsec, GJ 273 is also the closest known planetary system with a planet in the HZ after Proxima Centauri (Anglada-Escudé et al. 2016). Compared to the flaring star Proxima Centauri, we highlight the quietness of GJ 273. This makes GJ 273b even more attractive because of its lower atmospheric erosion, which translates into a more favorable environment for habitability.

We found that GJ 628 is orbited by at least three planets. We found a Sun-like magnetic cycle and clues that the stellar rotation is about 128 days. From this, it follows that the 67-day RV signal is more likely produced by stellar activity, and not by a planetary companion as reported by Wright et al. Here we found that a model including three Keplerians with P = 4.89, 17.87, 217.2 days is favored. This solution corresponds to planets with msin(i) = 1.9, 3.4, 7.7 M_⊕. For GJ 628c we derive a mean orbital distance slightly larger than the Wright et al. value (using our uncertainty we differ by 1.67σ; their uncertainty for a is disturbingly low). At this distance and compared to the Sun-Earth system, GJ 628c receives 1.29 times the flux received at Earth.

New radial velocities of GJ 3293 allow us to confirm the detection of the planet previously suggested in Astudillo-Defru et al. (2015). Additionally, we report the detection of an extra RV signal. The planet we confirm has a minimum mass of 7.6 M_⊕ and orbital period of 48.1 days; it is therefore inside the habitable zone of GJ 3293. The additional planet is a super-Earth of msin(i) = 3.3M_⊕ and an orbit of P = 13.25 days. The previously detected planets are at P = 30.6, 122.6 days with msin(i) = 23.5, 21.1 M_⊕; here we report a slightly lower eccentricity for the farthest planet.

---

Acknowledgments

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. X.B., X.D., and T.F. acknowledge the support of the French Agence Nationale de la Recherche (ANR), under the program ANR-12-BS05-0012 Exo-atmos and of PNP (Programme national de planétologie). This work has been partially supported by the Labex OSUG@2020. X.B. acknowledges funding from the European Research Council under the ERC Grant Agreement n. 337591-ExTrA.

N.C.S. acknowledges the support by Fundação para a Ciência e a Tecnologia (FCT) (project ref. PTDC/FIS-AST/1526/2014) through national funds and by FEDER through COMPETE2020 (ref. POCI-01-0145-FEDER-016886), and through grant UID/FIS/04434/2013 (POCI-01-0145-FEDER-007672). N.C.S. was also supported by FCT through the Investigador FCT contract reference IF/00169/2012 and POPH/FSE (EC) by FEDER funding through the program “Programa Operacional de Factores de Competitividade – COMPETE”.

We are very thankful to Jean-Baptiste Delisle for his precious comments about the dynamics of planetary systems. We thank the anonymous referee for the comments that improved the manuscript.

This research made use of the databases at the Centre de Données astronomiques de Strasbourg (http://cds.u-strasbg.fr), NASA’s Astrophysics Data System Service (http://adsabs.harvard.edu/abstract_service.html), and the http://arXiv.org paper repositories. This work made use of SCIPY (Jones et al. 2001–), IPYTHON (Perez & Granger 2007), and MATPLOTLIB (Hunter 2007).

References

Appendix A: Radial velocities

We list here the radial velocities in the barycentric frame – without subtracting the secular acceleration – and activity indicators. The uncertainty on RVs accounts for an instrumental error of 0.60 ms^-1 (added quadratically). A maximum likelihood approach is used to compute RVs, using a high signal-to-noise stellar template.

All Tables

All Figures