© M. Damasso et al. 2022

1 Introduction

Thanks to the most rewarding detection techniques, nowadays the search for and characterisation of low-mass planets orbiting within the habitable zone (HZ) of M dwarfs are popular topics, and particularly promising prospects for further investigation are provided by exoplanets detected around the nearest and brightest stars. Some of the earlier radial velocity (RV) surveys focused on low-mass stars (e.g. Endl et al. 2003; Zechmeister et al. 2009; Bonfils et al. 2013), and these are still among the main targets of new spectrographs developed and commissioned in the last decade (e.g. Mahadevan et al. 2012; Artigau et al. 2014; Affer et al. 2016; Quirrenbach et al. 2018; Pepe et al. 2021). Following the milestone Kepler/K2 mission, in the realm of transit photometry, the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2016) is presently the major provider of planet candidates around the brightest nearby stars, and M dwarfs are the main targets to detect small-sized planets in the HZ. Mostly RV and transit surveys have been used to make the first assessments of the occurrence rate of HZ planets around M dwarfs (e.g. Bonfils et al. 2013; Dressing & Charbonneau 2015; Hsu et al. 2020; Sabotta et al. 2021; Pinamonti et al. 2022), and stunning detections have been made in the last few years, especially around mid-to-late-type M dwarfs (e.g. Anglada-Escudé et al. 2016; Astudillo-Defru et al. 2017b; Cloutier et al. 2017; Gillon et al. 2017; Bonfils et al. 2018; Zechmeister et al. 2019).

So far, the detection and characterisation of HZ planets has proven to be more complicated for early-type M dwarfs (M0V–M2V), as testified by the lower number of discoveries compared to mid-to-late-type M dwarfs. We used the NASA exoplanet archive¹, complemented with the archive of confirmed potentially habitable low-mass planets maintained by the Planetary Habitability Laboratory (PHL)², and with the most recent papers on specific planetary systems for which no data are reported in either catalogue to compile a list of planets with radius R < 2.5 R_⊕ and/or mass m < 10 M_⊕, and an insolation flux S in the range 0.25–1.65 S_⊕. We specifically chose targets that have been detected around M-dwarfs with T_eff in the range 3650–3900 K. We found three transiting planets detected by Kepler, namely Kepler-186f (Quintana et al. 2014), Kepler-705b, and Kepler-1229b (Morton et al. 2016) (even the bona-fide KOI-4427.01 planet could be added to this group; see Torres et al. 2015). Because of the faintness of their host stars, masses have not been measured using RVs for any of these planets. The list of transiting planets also includes K2-3 d (Crossfield et al. 2015), a 1.6 R_⊕ planet orbiting near the inner edge of the HZ (P_d = 44.57 d; S_d = 1.6S_⊕) and detected by Kepler/K2. Its brightness (V = 12.2) compared to that of the previous systems makes K2-3 an excellent target for determining the planetary masses via high-precision Doppler spectroscopy. This task proved to be much more challenging than expected due to the difficulty in filtering out the signal induced by stellar activity. Only a mass upper limit has been determined for K2-3 d, despite the precise orbital ephemeris provided by transits and more than 400 RVs available (Damasso et al. 2018; Kosiarek et al. 2019). In the list compiled for M0V–M2V dwarfs, as described above, we found only one non-transiting planet detected using the RV technique, the ~7 M_⊕ super-Earth Gl 229 A c (Feng et al. 2020) in a red-brown dwarf binary system at a distance of 5.75 pc.

In the present paper, we report the detection of a super-Earth that moves on an eccentric orbit through the HZ of Gl 514 (BD+11 2576), an M0.5V–M1.0V star (V = 9 mag) located at a distance of 7.6 pc. The detection is based on nearly 25 yr of RV monitoring with the HIRES, HARPS, and CARMENES spectrographs. The manuscript is organised as follows. In Sect. 2, we summarise the main astrophysical properties of the host Gl 514. In Sect. 3, we provide general information about the spectroscopic and photometric dataset analysed in this work. In Sect. 4, we discuss the properties of the stellar magnetic activity of Gl 514 as derived from the analysis of spectroscopic activity diagnostics, and use the results to interpret the frequency content observed in the RV dataset (explored in Sect. 5). In Sect. 6, we present and discuss the detection of the super-Earth Gl 514 b based on an extensive analysis of 540 RVs. In Sect. 7, we present the analysis of TESS photometric data, and in Sect. 8, we discuss the astrometric sensitivity to wide-separation companions. Conclusions and future perspectives are discussed in Sect. 9. Additional results supporting our findings are shown in the appendix.

2 Stellar fundamental parameters

Gl 514 is a M0.5-M1 dwarf located at a distance of 7.6 pc from Earth. Its main astrophysical properties are summarised in Table 1. An earlier estimate of the stellar rotation period was provided by Suárez-Mascareno et al. (2015) from the analysis of spectroscopic activity indexes based on part of the same HARPS spectra analysed in this work. A following and independent result by Fuhrmeister et al. (2019), and based on CARMENES spectra (also part of our dataset), shows that P_{*, rot} = 30.8 ± 0.3 days and 30.3 ± 0.2 days, which were derived from the analysis of the Hα and CaII infrared triplet lines, respectively. As we discuss in Sects. 4 and 6.1, our derived values for P_*,rot are in agreement with previous determinations.

The Galactic space velocities U, V, and W of Gl 514 were derived using the Gaia coordinates and proper motions. We also employed the Gaia radial velocity to calculate the U, V, and W heliocentric velocity components in the directions of the Galactic centre, Galactic rotation, and north Galactic pole, respectively, with the formulation developed by Johnson & Soderblom (1987). We note that the right-handed system is used and that we did not subtract the solar motion from our calculations. The uncertainties associated with each space velocity component were obtained from the observational quantities and their error bars after the prescription of Johnson & Soderblom (1987). According to the resulting space velocities, Gl 514 has kinematics that deviate from that of known young stellar moving groups in the UV versus WV planes (e.g. Gagné & Faherty 2018). Therefore, the kinematical age of Gl 514 is likely >0.8 Gyr. Given its positive systemic RV, Gl 514 is moving away from us, and based on data from Gaia DR2, Bailer-Jones et al. (2018a) calculated that the closest approach to the Sun occurred nearly 30 Myr ago at a distance of 7.4 pc.

3 Overview of the spectroscopic and photometric dataset

Gl 514 was monitored with the HIRES, HARPS, and CARMENES spectrographs, collecting a total of 540 (104+162+274) RV measurements over 24 yr. HIRES observations coverthe time spanbetween 7 April 1997 and14 December 2014. We used the publicly available data derived by Tal-Or et al. (2019)³.

With HARPS, mounted at the ESO La Silla (Chile) 3.6 m telescope, the star was observed between 27 May 2004 and 26 April 2016. The spectra are publicly available through the ESO Archive⁴. The majority of the spectra (142 out of 162) were collected before the fibre upgrade intervention occurred in May 2015 (we label these as the pre-2015 dataset), and have a mean signal-to-noise ratio (S/N) of 82 measured at a reference wavelength of 5500 Å. The remaining 20 spectra collected after May 2015 (the post-2015 dataset) have mean S/N = 94. We extracted the HARPS RVs using different pipelines, all based on template-matching. One is the Template-Enhanced Radial velocity Re-analysis Application (TERRA V9.0; Anglada-Escudé & Butler 2012), an algorithm which has been proven to be particularly efficient with M dwarfs (e.g. Perger et al. 2017). We used the RVs extracted from orders ≥22 (corresponding to wavelengths λ ≥ 4400 Å, according to the notation adopted by TERRA), excluding the bluest and lower S/N orders from the RV computation (in this case, we did not use the spectral region between ~3800–4000 Å in the computation). This selection minimises the RMS and uncertainties of the RVs. Pre- and post-2015 spectra were treated as coming from different instruments, and two template spectra were calculated for each group. To account for the reported offset introduced by the fibre upgrade, when modelling the RVs we considered a zero-point (γ) and an uncorrelated jitter term (σ_jit) as free parameters for each dataset separately. We also exploited two other RV datasets extracted with alternative template-matching pipelines as a sanity cross check. One dataset was extracted with the NAIRA algorithm (Astudillo-Defru et al. 2015). The NAIRA recipe is based on a stellar template built from the median of the whole set of spectra, with uncertainties σ_RV computed following Bouchy et al. (2001), and the perspective acceleration subtracted using proper motion and parallax from Gaia DR2 (Gaia Collaboration 2016, 2018). The third dataset is represented by the nightly zero-point(NZP)-corrected RVs published by Trifonov et al. (2020)⁵, which were derived using the template-matching code SERVAL (Zechmeister et al. 2018). The RVs extracted with TERRA and NAIRA are listed in Tables A.1 and A.2.

The CARMENES spectrograph at the 3.5 m telescope of the Calar Alto Observatory in Spain (Quirrenbach et al. 2018) observed Gl 514 from 10 January 2016 to 2 April 2021, with only a small overlap with the last HARPS observations. In this work, we use the RVs extracted from the VIS channel, covering the wavelength range 5200–9600 Å. The data of the near-infrared (NIR) channel are characterised by a far lower S/N, and a mean σ_RV uncertainty of 8.7 m s⁻¹, and are therefore not particularly useful for detecting and characterising RV signals with a semiamplitude of few m s⁻¹, which are those we are interested in. Originally, CARMENES collected 321 spectra of Gl514, which were reduced using the CARACAL pipeline to obtain calibrated 1D spectra (Zechmeister et al. 2014; Caballero et al. 2016), from which we extracted the RVs using the template-matching code S ERVAL (Zechmeister et al. 2018). During a few nights, CARMENES observed the star more than once. We binned the corresponding RVs on a nightly basis⁶, and then discarded 14 measurements with σ_RV > 3 m s⁻¹ (nearly twice the mean value of σ_RV), resulting in a total of 274 RVs, which are listed in Table A.4.

A summary of the main properties of each RV dataset is provided in Table 2. The whole RV time-series is shown in Fig. 1 (for clarity, only the HARPS measurements extracted with TERRA are shown). HIRES RVs have the longest time-span and a less dense sampling. Due to its properties, we only use the HIRES dataset to test the presence of long-term signals, because it is not suitable to search for low-amplitude and shorter period signals compared to the other two dataset. CARMENES RVs are characterised by a typical internal error σ_RV of 1.6 m s⁻¹, equal to that of HIRES data, but the time-span is nearly three times shorter and the sampling is much denser. The larger RV uncertainties compared with those of the HARPS data are mostly due to the lower average exposure time of 440 ± 180 s, while for HARPS this is 900 ± 40 s.

Concerning the photometric observations, TESS monitored Gl 514 from 19 March to 15 April 2020 (sector 23). We analysed the long-cadence light curve extracted from the Full Frame Images (FFIs) using the PATHOS pipeline (Nardiello et al. 2019, 2020). Before the measurement of the flux of the target star, the light-curve extractor IMG2LC (Nardiello et al. 2015, 2016) subtracts all its neighbouring stars from each FFI by adopting empirical point spread functions (PSFs) and information from the Gaia DR2 catalogue (Gaia Collaboration 2018). The flux of the target star is measured with different photometric apertures (with radii from 1 to 4 pixels), fitting the empirical PSF. The systematic effects that change the raw light curve are corrected using the cotrending basis vectors (CBVs) extracted and applied as in Nardiello et al. (2020). We also analysed the short-cadence light curve released by the TESS team. We did not use the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) flux because of some systematic effects due to over-corrections and/or injection of spurious signals. We corrected the Simple Aperture Photometry (SAP) flux by applying CBVs obtained using the SAP light curves of the stars in the same Camera/CCD in which Gl 514 falls, and following the same procedure used for extracting the CBVs for the long-cadence light curves (Nardiello et al. 2020). Gl 514 falls on a CCD heavily affected by stray-light contamination, which in turn badly affects the overall quality of the data.

Fig. 1 Time-series of the HIRES (orange dots), HARPSTERRA (red and black dots), and CARMENES-VIS (green dots) RVs analysed in this work. The average has been subtracted from each dataset.

4 Stellar magnetic activity characterisation

We used spectroscopic diagnostics to characterise the time variability induced by the magnetic activity of Gl 514. This is a fundamental analysis to correctly interpret the nature of significant periodic signals possibly found in the RV time-series. We analysed the frequency content of the data for each individual instrument using the Generalised Lomb-Scargle algorithm (GLS; Zechmeister et al. 2009). The main peaks identified in the periodogram are summarised in Table 3 together with their false alarm probabilities (FAPs) calculated through a bootstrap with replacement analysis (i.e. an element of the dataset may be drawn multiple times from the original sample).

4.1 HIRES

We analysed the chromospheric S -index reported by Tal-Or et al. (2019). The time-series and GLS periodogram are shown in Fig. 2. The main peak occurs at 734 days, with almost the same power as its alias at 629 days (the alias frequency in the window function is 0.000229 day⁻¹, corresponding to P = 4361 days). Due to the low number of points sparsely sampled over a large time-span, we can only conclude that this period suggests the presence of a long-term modulation whose accurate properties cannot be determined. The other two significant peaks (FAP < 1%) occur at 30.7 and 14.9 days in the original data, and correspond to the stellar rotation period P_{*, rot} and its first harmonic (see Table 1).

4.2 HARPS

We derived activity diagnostics based on the CaII H&K, Hα, and NaI spectral lines using the code ACTIN V1.3.6 (Gomes da Silva et al. 2018). We also use the FWHM calculated from the cross-correlation function (CCF). We did not use the FWHM calculated by the standard data reduction system (DRS) of HARPS, because, as opposed to FGK-type stars, it does not perform a colour-correction in the case of M dwarfs, which stabilises the CCF against airmass and seeing variations. We found that, for M dwarfs, the colour-corrected FWHM time-series show a reduced RMS which allows better characterisation of the variability caused by stellar activity. Based on this improvement, we corrected the CCFs of Gl 514 by re-weighting them against a fixed flux distribution using files already available within the HARPS DRS. The time-series are shown in Fig. 3 together with the GLS periodograms, and listed in Table A3. For the FWHM, we show the periodogram after removing a long-term trend. Long-term variability (greater than a few hundred days) is detected for all the indices, although it is not possible to identify any periodicity (if present) unambiguously. Signals associated to P_{*, rot} appear significant (FAP < 0.1%) for the FWHM, CaII H and K, and NaI diagnostics. No significant signals are detected at the first harmonic of P_{*, rot}.

4.3 CARMENES

We derived spectroscopic activity diagnostics from the CARMENES-VIS spectra using the code SERVAL. These include the strength of emission lines of molecules sensitive to the chromospheric magnetic fields, such as the calcium infrared triplet (Ca-IRT) and the sodium doublet (NaD). We further derived the time-series of the differential line width of the spectral lines (dLW), and the chromatic index (CRX) by measuring the wavelength dependency of the RVs in the different orders of the spectra (see Zechmeister et al. 2018 for details). The time-series of these activity indicators are shown in the left panels of Fig. 4, and the data are listed in Table A5. All the time-series show a clear long-term trend, especially the Ca-IRT, NaD, and dLW indexes. For the Ca-IRT, NaD, and dLW indexes, we corrected this trend using a least-square linear fit. The right panels of Fig. 4 show the GLS periodograms for each activity diagnostic (Ca-IRT, NaD, and dLW: pre-whitened residuals; CRX: original data), together with the FAP levels determined through a bootstrap (with replacement) analysis. Evidence for signals related to the stellar rotation period, or its first harmonic, are present in all the periodograms, except for CRX. The peaks with the highest power occur in the range 410-636 days, suggestive of a mid-term modulation with no immediate interpretation.

5 Frequency content analysis of the RVs

We investigate the frequency content of the different RV dataset and their combinations by calculating the maximum-likelihood periodograms (MLPs; Zechmeister et al. 2019) using the implementation of the MLP code included in the EXO-STRIKER package⁷. Differently from the commonly used generalised Lomb-Scargle periodogram (GLS; Zechmeister et al. 2009), for each tested frequency, the MLP shows the difference d(lnL) between the logarithm of the likelihood function corresponding to the best-fit sine function and that of a constant function. The MLP algorithm includes RV zero points as free parameters in the cases where datasets come from different instruments, as well as instrumental uncorrelated jitter terms. In the following, first we report the results for the HIRES, HARPS, and CARMENES data separately. We then inspect the periodograms for the combined HARPS+CARMENES, and HIRES+HARPS+CARMENES time-series. The FAPs indicated in each case are analytical and are calculated by Exo-Striker.

Fig. 2 Time series and GLS periodogram of the chromospheric activity S-index derived from HIRES spectra (first and second row). The levels of FAPs are indicated as horizontal dashed lines (green: FAP 0.1%; yellow: 1%; red: 10%), and have been determined through a bootstrap analysis. Third row: Window function of the data.

Fig. 3 Time-series and GLS periodograms for activity diagnostics derived from HARPS spectra. Rows 1–4: time series (left panels) and periodograms (right panels) of the CCF FWHM, and activity indexes derived from the spectral lines CaII H&K, Hα and NaI. For the FWHM, we show the periodogram of the residuals, after removing the long-term trend clearly seen in the time series. The levels of FAPs are indicated as horizontal dashed lines (green: FAP 0.1%; yellow: 1%; red: 10%), and are determined through a bootstrap analysis. Last row: Window function of the data.

Fig. 4 Time-series of activity diagnostics extracted from CARMENES-VIS spectra of Gl 514 (left panels) and their GLS periodograms (right panels). Except for the CRX index, the other periodograms were calculated after a pre-whitening by removing the clear long-term trend visible in the time-series. The levels of FAPs, determined through a bootstrap analysis, are indicated as horizontal dashed lines (green: FAP 0.1%; yellow: 1%; red: 10%. Last panel: Window function of the CARMENES-VIS data.

5.1 HIRES RVs

The MLP of the HIRES RVs (Fig. 5) is dominated by a significant peak at 14.97 days. We note that the peak observed at 15.02 days is likely an alias due to the sampling. We recall that a quite significant peak at 14.9 days is present in the periodogram of the S-index extracted from HIRES spectra, and therefore the peaks observed in the RVs are likely related to the first harmonic of P_{*, rot}, and should be attributed to stellar activity. In addition to a less significant peak at ~9.4 days, which could be related to the second harmonic of P_{*, rot}, there is one peak at lower frequency (P = 2746.76 days) with low significance, which is likely related to the long-term modulation that can be guessed by eye looking at the time-series in Fig. 1.

Fig. 5 Maximum-likelihood periodogram of the HIRES RVs. The periods corresponding to the main peaks are indicated, as are the analytical FAP levels.

5.2 HARPS RVs

Figure 6 shows the MLP of the full HARPS RV dataset for TERRA, NAIRA, and Trifonov et al. (2020) datasets⁸. Signals in three different frequency ranges appear significant. As in the case of HIRES data, the highest peak occurs at ∼15 days. Peaks with FAP = 0.1% (TERRA) occur at ∼33 and 36 days, and are likely related to P_{⋆, rot}. For the NAIRA dataset, their significance is slightly lower. The MLP of TERRA and Trifonov et al. (2020) data shows another peak with FAP < 0.1% at a period of 138 days, which has a counterpart at ∼147 days (FAP ∼ 1%) in the MLP of the NAIRA RVs. Due to its high significance in two datasets over three, the nature of this signal needs to be examined in more detail using more sophisticated modelling than a simple and inaccurate pre-whitening. A more thorough investigation of this signal is the main focus of our study. We note that there is no evidence for the peak at ∼9 days observed in the MLP of HIRES data, and that there are no significant peaks at low frequencies.

Fig. 6 Maximum-likelihood periodogram of the HARPS TERRA (first panel), NAIRA (second panel), and Trifonov et al. (2020; third panel) RVs (pre- and post-2015 combined dataset). The periods corresponding to the four highest peaks are indicated, as are the analytical FAP levels.

5.3 CARMENES RVs

The MLP of the CARMENES-VIS RVs is shown in the upper plot of Fig. 7. The RVs show very significant signals at 354 and 478 days and at nearly half the rotational period (15.8 days). We also note the peak at 180 days, with lower FAP. By removing the 354-day signal with a pre-whitening, the 478-day and 180-day signals both disappear, meaning that they are related to each other, while the signal related to the stellar rotation becomes the most significant (second panel of Fig. 7). The pre-whitening increases the significance of the 140-day signal. However, its significance in the CARMENES MLP is low, but the MLP itself is dominated by a signal at 15.8 days. We emphasise that the signal at 354 days is detected only in the CARMENES RVs, and that, despite its high significance, we cannot attribute it to a companion of Gl 514; otherwise it would have been detected even in the periodogram of the HARPS RVs. This signal is likely spurious, and might be due to micro-telluric spectral lines which are not masked by the standard RV-extraction pipeline SERVAL, as a few preliminary tests suggest. Given its high significance, we took this signal into account when modelling the RVs, and treated it as a sinusoid to fit only the CARMENES RVs. This choice appears reasonable because the MLP provides evidence for the presence of a very significant sinusoidal modulation.

Fig. 7 Upper panel: maximum-likelihood periodogram of the original CARMENES-VIS RVs. The periods corresponding to the highest peaks are indicated, as are the analytical FAP levels. Lower panel: maximum-likelihood periodogram of the pre-whithened RVs after removing the 354-day signal.

5.4 HARPS+CARMENES RVs

The MLPs of the combined HARPS and CARMENES-VIS dataset are shown in Fig. 8 for all the different HARPS RV dataset, together with the window function. The dominant peak is located at 15.82 days in all cases, with a best-fit semi-amplitude of the model sinusoid of 1.2–1.3 m s⁻¹. The peak at ~140 days reaches a FAP of 0.01% with the HARPSterra RVs, therefore becoming more significant after combining the HARPS and CARMENES RVs together. This period becomes more clear and its significance increases in the case of the HARPS_NAIRA RVs, while the FAP increases when using the RVs extracted by Trifonov et al. (2020). In all cases, the peaks with FAP ~ 1% delimit the range of periods where P_{*, rot} is expected to be located.

5.5 HIRES + HARPS + CARMENES RVs

The MLP of the joint HIRES, HARPS_TERRA, and CARMENES-VIS dataset is shown in Fig. 9. It is dominated by activity-related signals, and the inclusion of the HIRES data does not increase the significance of the ∼140-day signal, which has a higher FAP close to 0.1%.

In summary, the most remarkable results from the frequency content analysis are that the RVs derived from each instrument are dominated by signals related to stellar activity (corresponding to the first harmonic of the 30-day stellar rotation period), and that a peak at ∼140 days appears in the periodograms of the HARPS and CARMENES dataset separately, with different levels of significance. This peak does not have a counterpart in the periodograms of the activity diagnostics, and becomes significant (FAP ~ 0.01%) when the HARPS and CARMENES data are combined together thanks to the increased time-span of the whole dataset. In light of this, we use more sophisticated models to fit the RVs which take into account the presence of activity-induced variability, and examine the possible planetary nature of the 140-day signal in the combined HARPS+CARMENES dataset.

Fig. 8 Maximum-likelihood periodograms of the joint HARPS (preand post-2015) and CARMENES-VIS RVs (panel 1: HARPSTERRA; panel 2: HARPSNAIRA; panel 3: HARPSTrifonov et al.). The periods corresponding to the highest peaks are indicated, as are the analytical FAP levels. Offsets have been subtracted from the HARPS and CARMENES RVs, as derived by the MLP procedure. Lower panel: window function of the combined dataset.

Fig. 9 Upper panel: maximum-likelihood periodogram of the joint HIRES, HARPS (TERRA; pre- and post-2015), and CARMENES-VIS RVs. The periods corresponding to the highest peaks are indicated, as are the analytical FAP levels. Lower panel: Window function of the combined dataset.

6 Radial velocity analysis

We identified three different GP kernels, and based on the results in the previous section, we deem the use of all of them for testing on the RVs of Gl 514 to be physically justified. They are the standard quasi-periodic (QP) kernel (e.g. Damasso et al. 2018); the recently proposed quasi-periodic with cosine (QPC) kernel (Perger et al. 2021); and the so called ‘rotational’, or double simple harmonic oscillator (dSHO) kernel. In this study, we use all of them to fit the correlated stellar activity signal, and verify whether or not the detection of the planetary candidate signal depends on the choice of the kernel. A significant detection with all three kernels would indeed strongly support the case for the existence of the 140-day signal.

Application of the rotational kernel to RVs has not been discussed many times in the literature (see, e.g. Benatti et al. 2021 for an application to a young and very active star), therefore little is known about the performance of this kernel in recovering planetary Doppler signals compared to that of the quasi-periodic family. The properties and hyper-parameters of each GP model are described in detail in Appendix B. The QPC and rotational kernels contain terms that explicitly depend on both P_{*, rot} and its first harmonic, and therefore they appear particularly suitable for modelling the activity term in the RVs of Gl 514, assuming that the dominant signal at ~15 days can be attributed to stellar activity.

We modelled the stellar activity term using the GP regression package george (Ambikasaran et al. 2015). We tested models with or without planetary signals (testing both circular and eccentric orbits), and explored the full (hyper-)parameter space with the publicly available Monte Carlo (MC) nested sampler and Bayesian inference tool MultiNest v3.10 (e.g. Feroz et al. 2019), through the pyMultiNest wrapper (Buchner et al. 2014). The priors used throughout the analysis described below are summarised in Table 4. Hereafter, subscripts b, c, and d refer to the signal with a period of ~ 140 days, to a possible innermost planet, and to a possible outermost planetary signal, respectively. The upper limit of the prior on the orbital period P_b was set to 200 days based on the results of the MLP analysis, and in order to guarantee an unbiased analysis. The MC sampler was setup to run with 500 live points and a sampling efficiency of 0.5 in all the cases considered in our study. We performed a model comparison analysis by calculating the difference ∆ ln Z between the natural logarithm of the Bayesian evidences Z determined by MultiNest for each tested model. The a-priori probability is assumed to be the same for each model, and we follow the scale in Feroz et al. (2011) to assess their statistical significance. Hereafter, first we present the results of the analysis of the combined HARPS+CARMENES dataset, which is the most suitable – due to the overall high number of measurements, dense sampling, and precision - for testing the presence of a low-mass companion in the HZ. We then performed an analysis including the HIRES RV dataset in addition, which, because of its large time-span, is suitable for testing the existence of a long-term signal.

6.1 HARPS+CARMENES dataset

We tested the three GP kernels on the 378 RVs collected with HARPS and CARMENES-VIS. In addition to the candidate planetary signal, we included a sinusoid to model the one-year signal only in the CARMENES-VIS data (see Sect. 5.3). Concerning the choice of the prior for the period P_365–day of this signal, initially we considered a uniform prior in the range 300–600 days in order to also sample the ~480-day signal seen in the MLP of Fig. 8. The aim was to check which period would have been selected by the MC sampler. The resulting posterior is not multi-modal, and it appears symmetric around 355 days (±8 days), as expected given that the highest peak in the MLP occurs at this period. This test allowed us to safely restrict the range of the priors for P_365–day and T_{0, 365–day}, with a considerable decrease in the computing time. We tested both models with the eccentricity e_b of the planet candidate either fixed to zero or treated as a free parameter together with the argument of periastron ω_*,b, using the parametrization and . Taking into account the different wavelength range covered by HARPS and CARMENES, the spectrographs could be sensitive to stellar activity at a different level. Therefore, even at the cost of increasing the number of free parameters, we adopted different GP amplitudes for each instrument to get a more reliable fit (the same approach is used when including the HIRES RVs, as described below).

As stated, we repeated the same analyses using HARPS RVs extracted with alternative pipelines. We disclose that we obtained best-fit values of the planetary parameters which are all in agreement within the errors for all three RV-template-matching extraction methods. This is indeed an important outcome which supports the detection of Gl 514b. For the sake of simplicity, and with no loss of information, hereafter we discuss only the results obtained using TERRA. A summary of the results obtained for the main planetary parameters using NAIRA and Trifonov et al. (2020) RVs is provided in Appendix D. The main outcomes of the GP analysis are as follows. Independently from the GP kernel, the model including a Keplerian for a candidate planet is the most significant (strong evidence in all the cases, i.e. ln Z_1p-ln Z_0p > 5); the planetary-like signal is detected in all the three cases with K_b ~ 1.2m s⁻¹(5.5–6σ significance), an orbital period of P_b ~ 140 days, determined with an error bar lower than 1 day, and a 3σ significant eccentricity e_b ~ 0.4–0.5; and we get moderate-to-strong evidence (∆ ln Z = 2.5–2.8) for the model with a Keplerian over a model with a planet on a circular orbit in two out of three cases⁹. Overall, the model with the lowest Bayesian evidence is that using the QP kernel, while the highest statistical evidence is found for the model with the QPC kernel. We therefore elect the latter as our reference model. We emphasise that the parameters of the planetary signal are all in agreement for each tested GP kernel.

The best-fit values for the parameters of this model are reported in Table 5, while we summarise the results for the QP and dSHO kernel in Tables C.1 and C.2. Figure 10 shows the spectroscopic orbit of the candidate planet Gl 514 b based on the best-fit model, and in Fig. C.3 we show the QPC component of the RV time-series corresponding to stellar activity. For comparison, we show in Fig. C.2 the spectroscopic orbit of Gl 514 b for the circular case (QPC kernel).

Concerning the results for the hyper-parameters of the correlated activity term QPC, we note that the stellar rotation period θ is retrieved with high-precision and is well in agreement with the literature values, and the characteristic evolutionary time-scale λ_QPC is nearly four times larger than the value of θ. It is not unusual that, for early-type M dwarfs with rotation periods similar to that of Gl 514, the value of λ_QPC obtained from the fit of the RVs is of the same order of magnitude as the value of P_{*, rot}, as found for instance in Gl 686 ( and λ_QPC ~ 47 days; Affer et al. 2019), K2-3 ( and λ_QPC ~ 80 days; Damasso et al. 2018), GJ 3998 ( and λ_QPC ~ 69 days; Affer et al. 2016), or Gl 15A ( and λ_QPC ~ 72 days; Pinamonti et al. 2018), where the RVs were modelled with a quasi-periodic GP kernel, and we used the relation λ_QPC = 2 · λ_QP for deriving the timescales reported in parentheses.

6.1.1 Stability of the 140-day signal over time

If the 140-day signal is due to a companion of Gl 514, then its properties, such as semi-amplitude and period, must tend to values that are in agreement within the error bars as a function of the progressive increase in the number of RVs. To verify this behaviour, we repeated the GP QPC fit described above on the HARPSTERRA+CARMENES RVs (keeping all the priors unchanged) after removing the last 100, 150, and 200 CARMENES measurements. This corresponds to a reduction of 23%, 34.4%, and 45.9% of the total HARPSTERRA+CARMENES data, respectively, and to a decrease in the total time baseline (equal to 6154 days) of 675, 797, and 1059 days. The posterior distributions of the semi-amplitude K_b, orbital period P_b, and eccentricity e_b of the candidate planet Gl 514 b derived for each cropped dataset are shown in Fig. 11, and are compared with the posteriors obtained for the whole sample of HARPSTERRA+CARMENES RVs. The posteriors of K_b are all in agreement within the uncertainties, and those of P_b become narrower and move closer to that of the whole dataset with the increasing number of RVs. The posteriors for e_b are all in agreement within the uncertainties, and the eccentricity moves to higher values and becomes more significant with the increasing number of data. The model including a Keplerian for planet b becomes more significant over the model without the planetary signal with the increasing number of RVs, as indicated by the values of the Bayesian evidence differences ∆ ln Z in the plot legend. This check demonstrates the persistence of the 140-day signal over time. In Appendix E we present an independent cross-check analysis to test the nature of this signal. The results of this latter analysis also support the planetary hypothesis.

Fig. 10 Spectroscopic orbit of Gl 514 b (blue solid line) based on the best-fit solution shown in Table 5 (case eb ≠ 0). The error bars of the measurements include uncorrelated jitter terms added in quadrature to the formal RV uncertainties. Cyan dots correspond to RV data averaged over 15 phase bins.

6.1.2 Testing the two-planet model

Given the large number of HARPSTERRA+CARMENES RVs, and their dense sampling, we investigated the existence of an additional planetary signal, focusing on orbits internal to that of Gl 514 b. The RMS of the RV residuals of the one-planet GP models is in the range 1.2–1.4m s⁻¹, and Fig. C.4 shows their GLS periodograms calculated for each GP model. The periodograms in blue, which correspond to residuals with the activity signal not removed from the original dataset, show power at periods longer than P_*,rot. As signals with periods longer than P_*,rot are usually strongly suppressed in the residuals of a GP model, the periodograms shown in red, corresponding to residuals with the activity signal removed from the original dataset, are not particularly informative when searching for signals with P > P_*,rot. From this analysis, we expect to get at best hints for a possible sub m s⁻¹ signal worthy of future follow-up with extreme-precision RVs. The setup is the same as that adopted for the previous analysis, and we tested all the GP kernels used for the case with one Keplerian. To keep the analysis unbiased, we adopted uniform and large priors to sample the orbital periods, setting them to U(0,100) and U(100,200) days. In the following, the second planetary signal is identified with the subscript c, while the subscript b still refers to the 140-day planet.

Fig. 11 Posterior distributions of the semi-amplitude Kb (top panel), orbital period Pb (middle panel), and eccentricity eb (bottom panel) of the candidate planet Gl514 b derived from a GP QPC analysis of the HARPSterra+CARMENES VIS RVs after removing an increasing number of CARMENES measurements, beginning from the last epoch. Vertical lines indicate the 16th, 84th (dashed), and 50th (solid) percentiles for the whole HARPSTERRA+CARMENES dataset. Here, ∆ ln Z denotes the Bayesian evidence difference between a model with and without a Keplerian included.

Test 1: QP kernel

For planet Gl 514 b we find , days, and , and those of the possible innermost companion are K_c = 0.74 ± 0.31 m s⁻¹, days, and . The full set of posteriors are shown in Fig. F.1. The Bayesian evidence is ln Z = −960.4, which is equal to the Bayesian evidence for the QP model that includes only one Keplerian.

Test 2: QPC kernel

The retrieved parameters for planet b are K_b = 1.21 ± 0.17 m s⁻¹, days, and , while those of the possible innermost companion are m s⁻¹, days, and . The Bayesian evidence is ln Z = −953.9, which is slightly higher than the Bayesian evidence for the model that includes only one Keplerian (∆ ln Z = +1.3). This is not enough to claim that this model is significantly favoured over the one with only the Keplerian for planet b included, but it is suggestive that a second Keplerian signal at a shorter period could be present. The posteriors for this model are shown in the second panel of Fig. F.1.

Test 3: dSHO kernel

For the candidate planet b, we find days, and , while for the possible innermost companion we find K_c = 0.88 ± 0.30 m s⁻¹, days, and . We note that P_c is well constrained, and in agreement with the results of the QP and QPC kernel. However, this model is not statistically favoured over the simpler one-planet model. The Bayesian evidence is ln Z = −956.2, which is only slightly higher than the Bayesian evidence for the model that includes only one Keplerian (∆ ln Z = +0.5). The posteriors are shown in the third panel of Fig. F.1.

Table 6 summarises the Bayesian evidence for the two-planet models. We did not find strong evidence for an additional companion orbiting at a closer distance from Gl 514 than planet b, whose main parameters remain unchanged with respect to the model with one Keplerian. However, two of three models are characterised by a slightly higher Bayesian evidence and are not disfavoured, at least suggesting the presence of a signal with a period of close to 64 days and a sub m s⁻¹ semi-amplitude that is significant to ~2.5–3σ. This signal comes with an eccentricity of around 0.5⁻0.6, and this naturally raises concerns about the dynamical stability of such a two-planet system, especially against orbit crossing. To assess whether there are stable orbital configurations that are compatible with our solutions, we repeated the analysis for the QPC and dSHO cases using a more complete dynamical model (Almenara et al., in prep.; Rein & Liu 2012; Rein & Tamayo 2015) with an unchanged setup, but including the orbit inclination angles as free parameters. In our models, we investigated the scenarios corresponding to co-planar and non-co-planar orbits, to assess the dynamical effects linked to a non-zero relative orbital inclination angle. After sampling from the posterior (Foreman-Mackey et al. 2013), we got 100 000 posterior samples (each with a different orbital configuration) for each of the four scenarios. To these, we applied a ‘filter’ in order to select only the configurations that guarantee the dynamical stability of the system over 10⁵ orbits of Gl 514 b. We adopted the following stability criteria: (i) avoiding orbit crossing, and (ii) ensuring a MEGNO chaos indicator (Cincotta & Simó 2000; Cincotta et al. 2003) in the range between 1.99 and 2.01. In the cases of both co-planar and non-co-planar orbits, we find that there are thousands of possible stable configurations that survive our stability filtering, with 11.46%/3.76% and 3.6%/1.63% of the whole posterior samples for the QPC and the dSHO kernels, respectively. The outcome of this analysis is that, in principle, a model with two Keplerians cannot be ruled out with our data based on dynamical stability criteria, despite it is weakly favoured at best (i.e. in the case of the QPC model). Given the complexity of the problem, different analysis techniques and approaches could be used to further investigate the significance of a two-planet model, and collecting additional RVs with higher precision could be worthwhile given that the signal we find would be compatible with a second planet moving through the HZ.

6.2 HIRES+HARPS+CARMENES dataset: Testing the existence of a longer period companion

Up to this point in our investigation, the existence of Gl 514 b is the most likely hypothesis according to our analysis. It is interesting to search for an external companion that could be responsible for the observed high eccentricity. Here we examine this possibility by first including the HIRES measurements to the RV dataset, then searching for astrometric anomalies in Gaia and Hipparcos data (Sect. 8). As the HIRES data are not sensitive to the 140-day signal, to speed up the analysis without affecting it, we examined the RV residuals obtained by removing our adopted best-fit solution for the Keplerian of planet b from all the datasets. The MLP of these data is shown in Fig. 12. The periodogram does not show significant signals at long periods, as already visible in Fig. 9. Based on these premises, we considered it sufficient to test only the GP model with the dSHO kernel which, as we demonstrated, provides a reliable modelling of the correlated activity signal and is computationally less demanding. For the Keplerian parameters, we obtain m s⁻¹, days (with an over-density region around 1300 days), and an unconstrained eccentricity. The Bayesian evidence for this model (ln Z = −1214.8) is lower than that for the model without the Keplerian included (∆ ln Z = −2.1), and therefore it is not statistically significant over a pure correlated noise model.

To explore the presence of a longer period companion, we fitted the RV residuals including an acceleration term in place of a Keplerian. We obtain m s⁻¹ day⁻¹, which differs from zero only within 1.5σ. The lower Bayesian evidence (ln Z = −1222) confirms that this model is statistically less significant than that including a Keplerian or the reference model (only GP). In conclusion, based on our RV data and analysis framework, we do not find statistical evidence for the presence of an external companion to Gl 514b.

7 TESS photometry analysis

We searched for transit-like signals in the short-cadence TESS light curve from sector 23, following the procedure described in Nardiello (2020). To remove the imprint of stellar variability and other trends, we modelled the light curve with a fifth-order spline defined over a grid of N_K knots with a space of 24h from one to another; we also removed all the points with the TESS quality parameter DQUALITY>0, high values of the background (>5σ_sky above the mean sky value), and all the points with flux more than 4σ above the mean of the flattened light curve. The light curve is shown in Fig. 13. We calculated the box least squares (BLS) periodograms (Kovács et al. 2002) of the flattened light curve searching for transit-like signals with a period in the range between 1 day and the time-span of the light curve. We found a peak in the BLS periodogram at ~2.5 days (Fig. 13) associated to a signal detection efficiency (SDE) of ~3.9, and a S/N of ~27. The depth of the transit model fitted to this signal is found to be ~ 73 ppm, which would correspond to a planet with radius R ~ 0.5 R_⊕. Assuming an empirical threshold SDE ≥9 to claim a significant detection, we conclude that no transits are detected in the light curve of sector 23.

In our work, we are also interested in searching for a single transit event of Gl 514 b, the occurrence of which is associated to a geometrical transit probability of 0.5%. An analysis based on BLS can be useful for detecting events even when only one transit falls within the time-span of the observations but, in such a case, a detection threshold defined using SDE appears meaningless, as we show below. We investigated the likelihood that the signal we detected in the TESS data corresponds to a single transit produced by a planet with the orbital period of Gl 514 b. To this purpose, we devised injection and retrieval simulations. Each simulated dataset is built by injecting into the original light curve a transit signal produced by a planet with an orbital period of 140 days and a radius selected from a grid of values (0.5, 0.75, 0.85, 1.0, and 2.0 R_⊕). We constrained the phase of the orbit in such a way that the transit falls within the time-span covered by the light curve. For each radius, we performed 10 simulations, moving the time of central transit T_conj recursively two days ahead of the T_conj generated in the previous simulation. In this way, we checked how the photometric systematic errors that locally characterise the light curve affect the detection efficiency of the transit. Each simulated dataset was analysed with BLS in the same way as the original light curve. We flagged as recovered the transits for which T_{conj, injected} − T_{conj,recovered} < 0.1 d. We recovered all the injected transits with R ≥ 1.0 R_⊕. All of these have a transit depth in agreement with that of the injected model: for R = 1R_⊕ we recover depths in the range 310 ± 40 ppm (the expected value is 350 ppm), and for R = 2R_⊕ we obtain depths in the interval 1450 ± 170 ppm (the expected value is 1350 ppm). For all the recovered transits, the main peak in the BLS periodogram has S/N > 100, and the SDE values are very low (between 1 and 4), likely due to the presence of only a single transit. For this reason, we conclude that, in the case of Gl 514 data, the SDE is not a useful figure of merit to discriminate between recovered and missed single transits. For R = 0.85 R_⊕ we recovered 85% of the injected transits with a S/N > 85 and transit depths in the range 250 ± 35 ppm (the expected value is 240 ppm). The detection efficiency falls down for R = 0.75 R_⊕, and we are able to recover only 30 % of the injected transits with S/N > 80 and depths in the interval 210 ± 50 ppm (the expected value is 190 ppm). Finally, we could not detect transits for R = 0.50 R_⊕.

The results of the simulations allow us to conclude about the presence in the TESS light curve of a transit signal ascribable to Gl 514 b. Assuming an inclination angle of 90° for the orbital plane (i.e. forcing the planet to transit), the measured mass of 5.6 ± 0.9 M_⊕ would correspond to a radius , following the probabilistic mass-radius relation of (Chen & Kipping 2017). We show that we are able to detect a single transit due to a planet with a radius as large as that predicted for Gl 514 b within 1σ, but we find no evidence for such a transit in the TESS light curve of sector 23. The signal found by BLS corresponds to a radius for which the simulations provided a null result. The non-detection may be due to one of three factors: (i) the planet does not transit because of its geometrical configuration; (ii) the planet transits, but the planet radius is smaller than 1 R_⊕, and we are not able to detect it with high confidence; or (iii) the planet transits, but the transit does not fall within the time-span covered by the light curve of sector 23. Looking at Fig. 14, the latter option appears possible. The figure shows that there is a small likelihood that the time of inferior conjunction T_{conj, b} falls within the TESS observing window during sector 23. On the contrary, the same plot shows that there is a high probability that the time of centre transit will occur during sector 50, scheduled between 26 March and 22 April 2022. In the even more favourable case that the light curve from sector 50 has better quality than sector 23, we consider the detection of the transit to be a realistic possibility.

Fig. 12 Maximum likelihood periodogram of the HIRES, HARPS, and CARMENES RV residuals after removing only the best-fit signal of the candidate planet with an orbital period of ∼140 days from the original data (second column of Table 5). The best-fit sinusoid with period P = 357.76 days has been subtracted from the CARMENES-VIS data. Stellar activity has not been filtered out, because the best-fit model for the activity signal is linked only on RVs measured with HARPS and CARMENES.

Fig. 13 TESS short-cadence light curve of Gl 514. Panels a and b show the normalised light curves before and after the subtraction of the model, respectively; the model represented in red in panel a was calculated on a grid of knots spaced every 24h. Time is expressed as the TESS Barycenter-corrected Julian Day (BTJD = BJD - 2457 000). Panel c shows the BLS periodogram. The red triangle indicates the period associated with the highest S/N.

Fig. 14 Posterior distributions for the conjuction times Tconj,b of Gl 514 b that are close to TESS Sectors 23 (March-April 2020) and 50 (March-April 2022). The posteriors are derived from that obtained for our best-fit QPC model. The vertical lines mark the time-span of TESS observations.

8 Astrometric sensitivity to wide-separation companions

Gl 514 is astrometrically “quiet”. Its entry in the Gaia EDR3 archive reports values of astrometric excess noise and reduced unit weight error (RUWE) of 0.143 mas and 1.09, respectively. These numbers are typical of sources whose motion is well described by the standard five-parameter astrometric model (positions, proper motions, and parallax) in Gaia astrometry¹⁰. No statistically significant difference in the stellar proper motion at the Hipparcos and Gaia mean epochs is reported in the Kervella et al. (2022) and Brandt (2021) Hipparcos-Gaia catalogues of astrometric accelerations. Given that Gl 514 is really in the Sun’s backyard (its distance is 7.6 pc), the proper motion anomaly technique allows us to place rather interesting limits on the presence of wide-separation companions in the planetary and substellar mass regime. Using the formalism presented in Kervella et al. (2019) (Eqs. (13)-(15)), we show in Fig. 15 the mass-orbital separation sensitivity diagram from the Hipparcos-Gaia absolute astrometry. In the approximate range of 3–10 au, an object with ~0.2 M_J can be ruled out at the 1-σ level. A companion with the same value of m_p sin i_p at 6.5 au (the approximate limit in orbital period sampled by the combined RV dataset) would induce an RV modulation in Gl 514 with a semi-amplitude of ~3 m s⁻¹, which would have likely been picked up in the RV analysis. However, the absolute astrometry limits apply to companions with any orbital inclination, implying that the presence of such a low-mass gas giant can be ruled out even for quasi-face-on configurations (effectively producing overly small RV signals). The proper motion anomaly sensitivity diagram in Fig. 15 also indicates that massive giant planets and brown dwarf companions out to several tens of astronomical units (au) would have been detected if present. Jupiter-mass companions at Neptune-like separations would remain undetected. At 0.42 au, the upper limits on the true mass of the RV-detected planet in this work are not particularly strong (≲1.5 M_jup, with a constraint on the true inclination angle ≳1°).

Fig. 15 Detectability curve for companions of given mass and orbital distance around Gl 514 based on the proper motion anomaly technique described in Kervella et al. (2019).

9 Conclusions and perspectives

We presented an analysis of nearly 25 yr of RVs of the nearby M dwarf Gl 514 collected with the HIRES, HARPS, and CARMENES spectrographs. The data appear to be dominated by signals related to stellar magnetic activity, which we corrected testing three different GP kernels of proven efficacy on quasi-periodic modulations. In all cases, we found strong evidence for a signal which we attribute to the presence of a super-Earth with minimum mass m_b sin i_b = 5.2 ± 0.9 M_⊕ (with an 87% probability that the true mass is within a factor of two larger) moving on an eccentric orbit () with orbital period P_b = 140.43 ± 0.41 days and semi-major axis au. The parameters of our adopted best-fit model are highlighted in bold face in Table 5. This result was corroborated using RVs extracted from the HARPS spectra with three different pipelines, all based on a template matching technique. We also investigated the possibility that a signal induced by an additional companion to Gl 514b is present in the RVs. Exploring orbital periods P < 100 days, we did not find evidence in favour of a two-planet model, even though our solution suggests the existence of a sub m s⁻¹ signal around 64 days which we deem worthy of future follow-up using very high-precision RVs. We find no evidence for a longer period planet (P > 200 days). Available astrometric data rule out a ~0.2 M_jup planet at a distance of ~3–10au, and massive giant planets or brown dwarfs out to several tens of au, while Jupiter-mass companions at Neptune-like separations would remain undetected.

One point emerging from our analysis that deserves further attention is the 3σ significant orbital eccentricity of Gl 514 b. In the first place, the eccentricity is a relevant parameter in that it can provide crucial information about the evolutionary history of a planetary system, and it is therefore important to get an accurate and precise measurement, when possible. However, in our case it is necessary to proceed with caution, because Gl 514b could be an ‘eccentric impostor’. This is a well-known problem that has been extensively discussed in literature (e.g. Rodigas & Hinz 2009; Trifonov et al. 2017; Wittenmyer et al. 2013, 2019; Boisvert et al. 2018). Works such as those of Anglada-Escudé et al. (2010) and Kürster et al. (2015) have discussed in detail the issue that a system of two planets on 2:1 resonant and nearly circular orbits can be confused with a single planet on an eccentric orbit. However, their results apply to systems with very different properties from those of Gl 514: generally, they have sparse data sampling, and the ‘eccentric’ planet has a much greater semi-amplitude than that of Gl 514 b, and significantly greater than the RMS of the residuals of the one-Keplerian model. For Gl 514, the RMS of the residuals (1.35 m s⁻¹) is greater than the semi-amplitude of the Keplerian signal, and it does not appear to be feasible to distinguish between the two models. Nonetheless, we have a large number of RVs with a dense sampling, and we used GPs to efficiently model the stellar activity contribution. Therefore, we deemed it interesting to perform an analysis of the HARPSterra and CARMENES data anyway and check whether or not Gl 514 b could be an eccentric impostor composed of two planets on 2:1 resonant and nearly circular orbits. Our analysis shows that, independently from the choice of the GP kernel, the two-planet solution is statistically never favoured over that with a single eccentric planet. This is in agreement with the more general results we present in Sect. 6.1.2, where we show that Gl 514 b is still fitted as an eccentric planet when using a model with two Keplerians. On the same subject, more recent works (e.g. Hara et al. 2019 and Faria et al. 2022) discuss the issue where orbital eccentricities fitted using only RVs can be spurious, as a consequence of an inappropriate modelling or suboptimal data quality relative to the low semi-amplitude of the signal. Hara et al. (2019) showed how one can get an incorrect inference of the eccentricity when an uncorrelated jitter term is not included in the model and, more importantly, when correlated signals such as stellar activity are not properly modelled. Faria et al. (2022) showed how the value of the eccentricity of Proxima d depends on the method used to extract the RVs from ESPRESSO spectra, with a better constraint obtained using a technique based on template matching. In our analysis, we included uncorrelated jitter terms, we showed that the result is independent of the model used to fit the correlated stellar activity signal, and we tested our finding against different RV-extraction methods applied to HARPS spectra. We conclude that the eccentric solution for the orbit of Gl 514 b is likely non-spurious¹¹, although the RVs analysed in this work, and the analysis framework we have selected, do not allow us to obtain strong statistical evidence in favour of this model. Given the low ~ 1 m s⁻¹ semi-amplitude of the planetary signal, we acknowledge that further investigation of the signal and additional data are necessary to get a more accurate and precise measurement of e_b. Any future follow-up of Gl 514 with high-resolution spectrographs such as ESPRESSO, given a sufficiently dense sampling, could help in this regard. With this caveat in mind, for the rest of our discussion we assume that our result concerning e_b is realistic.

Fig. 16 Schematic representation of a few possible orbits for planet Gl 514 b (blue ellipses) randomly selected taking into account the parameter values and their errors of our adopted solution (Table 5, second column). The best-fit orbit is shown with a dashed black line. The red star identifies one of the foci of the ellipses occupied by the host star Gl 514 b, and the green annulus corresponds to the conservative habitable zone.

9.1 Gl 514 b in the context of planets orbiting in the habitable zone of M dwarfs

According to the theoretical calculations by Kopparapu et al. (2013, 2014), the conservative HZ for Gl 514 and for a planet with mass m_p = 5 M_⊕ extends between 0.207 and 0.411 au. Gl 514 b spends nearly 34% of its orbital period within the conservative HZ of its host star, as shown in the sketch of Fig. 16. Among the known low-mass exoplanets orbiting in the HZ of nearby M dwarfs¹², none have an eccentricity as high and as significant as Gl 514 b.

Although we have not yet been able to measure the radius of Gl 514 b or constrain its average composition and physical structure, it would be interesting to use Gl 514 b as a case study to investigate the habitability of an eccentric super-Earth orbiting a low-luminosity star using climate models. The question of the habitability of planets that experience insolation variations along their orbits and spend a considerable fraction of time outside the HZ is a complex problem (see, e.g. Williams & Pollard 2002; Kane & Gelino 2012; Bolmont et al. 2016; Méndez & Rivera-Valentín 2017). While Williams & Pollard (2002) concluded that long-term climate stability depends primarily on the average stellar flux received over an entire orbit, Bolmont et al. (2016) found that, for water worlds of higher eccentricities (or host stars with higher luminosity), the mean flux approximation becomes less reliable for assessing their ability to sustain a liquid water ocean at their surface. When discussing the equilibrium temperature of an exoplanet with orbital properties similar to those of Gl 514 b (with implications about climate and habitability), one fundamental parameter to be taken into account is the thermal timescale, which is defined as the timescale on which the planetary temperature fluctuations adjusts around the flux-averaged equilibrium value to the changing stellar irradiation (Quirrenbach 2022). Any temperature fluctuation is expected to be damped following an exponential decay, and the thermal timescale depends on the heat capacity per unit surface area of the planet. General considerations suggest that substantial liquid surface water reservoirs or atmospheres of several tens of bars can work as an efficient climate buffering, avoiding temperature fluctuations on short timescales. If these properties characterise Gl 514 b, its surface equilibrium temperature should be damped around the flux-averaged value, which to a first approximation we estimate to be 202 ± 11 K (assuming zero albedo). It is clear that, lacking any detailed knowledge of the physical and chemical properties, any consideration about the habitability of Gl 514 b is presently only speculative, and further discussion of this topic is out of the scope of this paper. Nonetheless, we support Gl 514 b as a benchmark system for investigating the habitability of a super-Earth using sophisticated climate models, and tools such as VPLANET (Barnes et al. 2020).

Using the same list of potentially habitable planets, compiled as described above, we put Gl 514 b in the context of other low-mass planets detected in the HZ of M dwarfs (Fig. 17). Planetary masses are represented versus the effective stellar temperature (panel (a)), insolation¹³ (panel (b)), and distance (panel (b)). We note that in panel (a), only two super-Earths are close to Gl 514 b, namely Gl 229 A c and K2-3 d, and therefore Gl 514 b enters the very small group of low-mass planets moving within the HZ of nearby stars with spectral type earlier than M2V. When looking at panel (b), Gl 514 b is located next to Gl 682 b, but this planet has a much shorter period (~17.5 d), and its minimum mass is far less precise. Gl 514 b has a minimum mass that is compatible with those of LHS 1140 b (real mass) and Gl 357 d (minimum mass) within the uncertainties, and an insolation that is ~25% lower. The potential habitability of Gl 357 d was discussed in detail by Kaltenegger et al. (2019). Gl 357 d is considered a prime target for observations with Extremely Large telescopes as well as future space missions, and the same expectations are even more valid for Gl 514 b, which is closer than Gl 357 d and has a brighter host star. Currently, we can only speculate about realistic outcomes of high-contrast imaging of Gl 514 b with the Extremely Large Telescope (ELT), considering the expected performance of the Planetary Camera and Spectrograph (PCS) that will be dedicated to detecting and characterising exoplanets with sizes from subNeptune to Earth-sized in the solar neighbourhood (Kasper et al. 2021). Assuming a maximum angular star-planet separation of ~80 mas (~55 mas for a circular orbit), PCS could in principle detect the planet if the planet-to-star I-band flux ratio is roughly greater than ~2 ∙ 10^-9, a threshold that the Gl 514 system could realistically exceed (see Fig. 1 in Kasper et al. 2021).

Fig. 17 Gl 514 b in the context of planets with mass or minimum mass <10 M⊕ that have been detected to date within the HZ of M dwarfs. The insolation flux of Gl 514 b is the temporal average orbit value.

9.2 Radio emission from subAlfvénic star-planet interaction

Very recent studies based on radio-frequency observations proposed the star-planet interaction as a possible mechanism to explain the detection of radio emission coming from nearby M dwarfs (Turnpenney et al. 2018; Vedantham et al. 2020; Callingham et al. 2021; Pérez-Torres et al. 2021). This interaction is expected to yield auroral radio emission from stars and planets alike, because of the electron cyclotron maser (ECM) instability (Melrose & Dulk 1982), whereby plasma processes within the star (or planet) magnetosphere generate a population of unstable electrons that amplifies the emission. In some favourable cases, the emission could be detectable. The characteristic frequency of the ECM emission is given by the electron gyrofrequency, v_G = 2.8 B MHz, where B is the local magnetic field in the source region in Gauss. ECM emission is a coherent mechanism that yields broadband (∆v ~ v_G/2), highly polarized (sometimes reaching 100%), amplified non-thermal radiation.

If the velocity v_rel of the plasma relative to the planetary body is less than the Alfvén speed, v_A, that is M_A = v_rel/v_A < 1, where M_A is the Alfvén Mach number, then energy and momentum can be transported upstream of the flow along Alfvén wings. The interaction between Jupiter and its Galilean satellites is a well-known example of subAlfvénic interaction, producing detectable auroral radio emission (Zarka 2007). In the case of star-planet interaction, the radio emission arises from the magnetosphere of the host star, induced by the exoplanet crossing the star magnetosphere, and the relevant magnetic field is that of the star, B_*, and not that of the exoplanet. As M-dwarf stars have magnetic fields ranging from about 100 G and up to above 2–3 kG, their auroral emission falls in the range from a few hundred MHz up to a few GHz, with flux densities that could be well detected by present and future radio telescope arrays.

Given the growing interest around the topic of star-planet interaction at low-frequencies, we believe it would be interesting to investigate this possibility for the case of Gl 514 b. We followed the prescriptions in Appendix B of Pérez-Torres et al. (2021) to estimate the flux density expected to arise from star-planet interaction at the frequency of 1.4 GHz, which corresponds to the cyclotron frequency of the local magnetic field of 500 G that we assume, which is a reasonable value for a star with a rotation period of 30 days, as indicated by the data and Fig. 5 in Shulyak et al. (2019). We computed the radio emission arising from star-planet interaction for two different magnetic field geometries: a closed dipolar geometry, and an open Parker spiral geometry. For the dipolar case, the motion of the plasma relative to Gl 514 b happens in the supra-Alfvénic regime. Therefore, no energy or momentum can be transferred to the star through Alfvén waves. However, in the open Parker spiral case, the plasma motion proceeds in the subAlfvénic regime. We show in Fig. 18 the predicted flux density as a function of orbital distance arising from the interaction of a magnetized exoplanet (1 G) with its host star. The yellow and blue shaded areas correspond to the predictions of two models, and encompass the range of values from 0.01 to 0.1 for the efficiency factor, є, in converting Poynting flux into ECM radio emission. We note that the Zarka/Lanza model (blue; Zarka 2007; Lanza 2009) predicts flux densities from a few hundred µJy up to a few mJy at the orbital distance of Gl 514 (dashed line). Those values are almost two orders of magnitude larger than expected in the Saur/Turnpenney model (yellow; Saur et al. 2013; Turnpenney et al. 2018), which predicts less than about 10 µJy in the most favourable case. If the magnetic field responsible for this putative cyclotron radio emission is of 500 G, as we assume here, then observations at frequencies in the 1.0–1.8 GHz range could be used to rule out one of the models above. However, we caution that the magnetic field of Gl 514 is not known, and our assumption is no more than an educated guess, and is likely to be uncertain to within a factor of two. As a consequence, observations at multiple frequencies starting from about 150 MHz and up to about 2–3 GHz would be advisable in order to better constrain the radio emission in this system.

Fig. 18 Expected flux density for auroral radio emission arising from star-planet interaction in the system Gl 514 as a function of orbital distance. The interaction is expected to be in the subAlfvénic regime (i.e. MA = vrel/vAlfv ≤ 1; top panel) at distances equal to the semi-major axis and periapsis (vertical dashed line).

Acknowledgements

We thank the anonymous referee for her/his useful comments. This work is partly based on observations collected with the CARMENES spectrograph, which is an instrument at the Centra Astronómico Hispano-Alemán (CAHA) at Calar Alto (Almería, Spain), operated jointly by the Junta de Andalucía and the Instituto de Astrofísica de Andalucía (CSIC). CARMENES was funded by the Max-Planck-Gesellschaft (MPG), the Consejo Superior de Investigaciones Científicas (CSIC), the Ministério de Economía y Competitividad (MINECO) and the European Regional Development Fund (ERDF) through projects FICTS-2011-02, ICTS-2017-07-CAHA-4, and CAHA16-CE-3978, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l’Espai, Institut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the MINECO, the Deutsche Forschungsgemeinschaft through the Major Research Instrumentation Programme and Research Unit FOR2544 “Blue Planets around Red Stars”, the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía. We acknowledge financial support from the Agencia Estatal de Investigación of the Mnisterio de Ciencia, Innovación y Universidades through project PID2019-109522GBC5[1:4]. The authors acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación and the ERDF “A way of making Europe” through projects PID2020-120375GB-I00, PID2019-109522GB-C5[l:4], and PGC2018-098153-B-C33, and the Centre of Excellence “Severo Ochoa” and “Maria de Maeztu” awards to the Instituto de Astrofísica de Canarias (CEX2019-000920-S), Instituto de Astrofísica de Andalucía (SEV-2017-0709), Centro de Astrobiologia (MDM-2017-0737), Institut de Ciències de l’Espai (CEX2020-001058-M), and the Generalitat de Catalunya/CERCA programme. They acknowledge financial contribution from the agreement ASI-INAF n.2018-16-HH.0. M. Damasso acknowledges financial support from the FP7-SPACE Project ETAEARTH (GA no. 313014). M. Perger and I. Ribas acknowledge the support by Spanish grant PGC2018-098153-B-C33 funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, by the programme Unidad de Excelencia Maria de Maeztu CEX2020-001058-M, and by the Generalitat de Catalunya/CERCA programme. M. Perger also acknowledges support from Spanish grant PID2020-120375GB-I00 funded by MCIN/AEI. D. Nardiello acknowledges the support from the French Centre National d’Etudes Spatiales (CNES). N. Astudillo-Defru acknowledges the support of FONDECYT project 3180063. M. Pérez-Torres acknowledges financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award to the Instituto de Astrofísica de Andalucía (SEV-2017-0709) and through the grant PID2020-117404GB-C21 (MCI/AEI/FEDER, UE). A. Suárez Mascareno acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) under 2018 Juan de la Cierva program IJC2018-035229-I. A.S.M. acknowledges financial support from the MICINN project PID2020-117493GB-I00 and from the Government of the Canary Islands project ProID2020010129. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work is dedicated to the memory of M.T.

Appendix A Dataset extracted from the HARPS and CARMENES spectra

We report in Tables A1 and A2 the RVs extracted from the HARPS spectra. The spectroscopic activity diagnostics discussed in this work are listed in Tables A3 and A5. The HARPS RVs extracted by Trifonov et al. (2020) are also available at https://www2.mpia-hd.mpg.de/homes/trifonov/HARPS_RVBank.html. Radial velocities derived from HIRES spectra are available at https://cdsarc.unistra.fr/viz-bin/cat/VMNRAS/484/L8.

Appendix B Description of the GP kernels tested in this study

Appendix B.1 Quasi-periodic (QP)

The QP covariance matrix (e.g. Haywood et al. 2014) has now become a standard tool to model the activity term in RV time-series showing variability modulated over the stellar rotation period. In our work, an element of the matrix is defined as follows: (B.1)

Here, t and t′ represent two different epochs of observations, σ_RV is the RV uncertainty, and δ_t/t′ is the Kronecker delta. Our analysis takes into account other sources of uncorrelated noise - instrumental and/or astrophysical – by including a constant jitter term σ_jit which is added in quadrature to the formal uncertainties σ_RV. The GP hyper-parameters are h, which denotes the scale amplitude of the correlated signal; θ, which represents the periodic timescale of the correlated signal, and corresponds to the stellar rotation period; w, which describes the ‘weight’ of the rotation period harmonic content within a complete stellar rotation (i.e. a low value of w indicates that the periodic variations contain a significant contribution from the harmonics of the rotation periods); and λ_QP, which represents the decay timescale of the correlations, and is related to the temporal evolution of the magnetically active regions responsible for the correlated signal observed in the RVs.

Appendix B.2 Quasi-periodic with cosine (QPC)

The QPC kernel was recently introduced by Perger et al. (2021), who found that in general it guarantees a better performance over of the QP kernel. The covariance matrix element of the QPC kernel, as implemented in our work, is (B.2)

Again, t and t′ represent two different epochs of observations; h₁ and h₂ are scale amplitudes; θ still represents the periodic timescale of the modelled signal, and corresponds to the stellar rotation period; w still describes the weight of the rotation period harmonic content within a complete stellar rotation; λ_QPC is defined as 2 · λ_QP, better representing the average lifetime of the activity-related features responsible for the stellar correlated signal in the RVs (Perger et al. 2021); σ_RV and σ_jit are the RV uncertainty and instrument-dependent jitter, respectively, and δ_t/t′ is the Kronecker delta.

Appendix B.3 Rotational (dSHO)

The ‘rotational’ kernel is defined by a mixture of two stochastically driven, damped simple harmonic oscillators (SHOs) with undamped periods of P_*,rot and P_*,rot/2. This can be obtained by combining two SHOTerm kernels included in the package celerite (Foreman-Mackey et al. 2017)¹⁴. The power spectral density corresponding to this kernel is (B.3)

where (B.4) (B.5) (B.6) (B.7) (B.8) (B.9)

The parameters in [B.4–B.9], where the subscripts 1 and 2 refer to the primary (P_*,rot) and secondary (P_*,rot/2) modes, represent the inputs to the SHOTerm kernels. However, instead of using them directly, we adopt a different parametrisation using the following variables as free hyper-parameters in the MC analysis, from which S_i, Q_i, and ω_i are derived through Eq. [B.4-B.9]: the variability amplitude A, the stellar rotation period P_*,rot, the quality factor Q₀, the difference ∆Q between the quality factors of the first and second terms, and the fractional amplitude f of the secondary mode relative to the primary.

Appendix C Additional material relative to one-planet models fitted to the HARPS_TERRA+CARMENES RV time-series

In Fig. C.1 we show the posterior distributions, and their two-by-two correlations (corner plot), for all the free (hyper-)parameters of our assumed best-fit solution for the Gl 514 system based on the analysis of the combined HARPS_terra+CARMENeS RVs.

We show in Fig. C.2 the spectroscopic orbit of Gl 514 b for the circular case, obtained using the QPC kernel to model the activity component. In Fig. C.3 we show two subsets of the HARPS_TERRA+CARMENES RV residuals (i.e. after subtracting the Keplerian for planet b, and the additional sinusoid for the CARMENES data only) together with the GP QPC best-fit solution for the correlated term. These data correspond to the stellar activity component present in the RV time-series, possibly with other not modelled small-amplitude signals included (astrophysical or instrumental), according to our adopted best model for the Gl 514 system (see the second column of Table 5 for the list of all the free parameters).

In Tables C.1 and C.2 we summarise the best-fit values of the parameters for the models including the QP and dSHO kernels, which are characterised by a Bayesian evidence lower than that for the QPC model. In Fig. C.4 we compare the GLS periodogram of the HARPS+CARMENES RV residuals for all the GP kernels, after removing the best-fit models shown in the second columns of Tables C.1–C.2. We show the cases both with and without the correlated GP-fitted signal included in the residuals. We note that, when the activity signal is not removed from the dataset, all the periodograms peak at 54 days, with semi-amplitude of the best-fit sinsuoid equal to 1.2 ± 0.2 m s⁻¹. Nothing significant is seen at frequencies larger than 0.1 d⁻¹. When the activity signal is removed, none of the periodograms show significant peaks, especially at low frequencies. This does not come as a surprise, because the suppression of low frequencies has been commonly observed in residuals of GP-filtered RVs, especially at frequencies larger than the stellar rotational frequency. Nonetheless, we note that, despite their low power, peaks at around periods of 50–60 days are not completely suppressed. Their nature is discussed in Sect. 6.1.2.

Fig. C.1 Posterior distributions of the free (hyper)parameters of our assumed best-fit model, including a Keplerian for planet b and a GP QPC correlated activity signal fitted to the RVs obtained from HARPS (TERRA dataset) and CARMENES VIS spectra.

Fig. C.2 Phase-folded RVs of HARPSterra and CARMENES-VIS showing the spectroscopic orbit of Gl 514 b for the circular case. The blue solid line is the best-fit solution using a QPC kernel to model the activity. The error bars include the uncorrelated jitter terms added in quadrature to the formal RV uncertainties. Cyan dots correspond to RV data averaged over 15 phase bins.

Fig. C.3 Sections of the residuals of the whole set of HARPSterra+CARMENES RVs (corresponding to the model shown in the second column of Table 5, showing the QPC correlated signal mostly related to variations in stellar activity. The error bars include the uncorrelated jitter terms added in quadrature to the formal RV uncertainties.

Fig. C.4 GLS periodograms of the HARPSterra+CARMENES RV residuals, after removing the best-fit models shown in the second columns of Table C.1–C.2. For each kernel, periodograms of data with and without the correlated activity signal included are shown. The periodograms were calculated with the uncorrelated jitter terms added in quadrature to the formal RV uncertainties.

Appendix D Results using alternative RV extraction pipelines applied to HARPS spectra

As discussed in Sect. 3, in this study we also examined two alternative RV dataset extracted from HARPS spectra, namely the NAIRA and Trifonov et al. (2020) nightly zero-point-corrected (NZP) dataset, all coming from techniques based on template matching as for the case of TERRA. We repeated the analysis described in Sect. 6.1 as a sanity cross-check for confirming the existence of the candidate planetary signal at ~140 days. The analysis were performed using the same MC sampler, setup, and priors as those adopted for the RVs extracted with TERRA. We summarise the results in Table D.1, by reporting (i) the best-fit values of some of the signal parameters (semi-amplitude K_b, period P_b, and eccentricity e_b) obtained for each tested GP kernel, (ii) the corresponding logarithm of the Bayesian evidence, ln Z, and (iii) the difference ∆ ln Z with respect to the model with no planetary signal included. Independently from the algorithm used to extract the RVs, the model with the eccentricity e_b treated as a free parameter is generally strongly favoured (∆ ln Z > 5 in five cases over six) over the model with no planet included. It is generally moderately-to-strongly favoured over the model with e_b fixed to zero ( in five cases over six). The eccentricity is retrieved with significance within ~3–4σ. These results are in agreement with those obtained for the HARPS TERRA dataset, and strengthen our conclusion that a planetary-like signal is present in the data. Moreover, the parameters K_b, P_b, and e_b are all consistent with those derived using the TERRA RVs. Based on the results shown in Table D.1, we conclude that the GP QPC plus a Keplerian is the best-fit model to the RVs of the two alternative pipelines, as in the case of TERRA, because the corresponding ln Z is significantly the highest among all the values of the Bayesian evidence. We note that the dSHO model provides the highest significance in favour of the ~140-d signal over the model with no planet included, as we found for the case of the TERRA dataset.

Appendix E Additional tests to cross-check the existence and temporal stability of the 140-d signal

To independently cross-check the results of Section 6, we analysed HARPSterra and CARMENES-VIS RVs with other methods. We first use the ℓ₁ periodogram (Hara et al. 2017), a tool based on a sparse recovery technique called the basis pursuit algorithm (Chen et al. 1998). The ℓ₁ periodogram takes in a frequency grid and a covariance matrix to model noise sources as inputs. It aims to find a representation of the RV time-series as a sum of a small number of sinusoids whose frequencies are in the input grid. It outputs a figure which has a similar aspect as a GLS periodogram, but with fewer peaks due to aliasing. FAPs can be calculated for each peak, and their interpretation is equivalent to that of common periodograms.

To determine the influence of the noise model, we followed (Hara et al. 2020) and considered a grid of covariance models and rank the alternatives with cross validation. We defined the covariance matrix V so that its element at index (k, l) is (E.1)

where σ_RV is the nominal measurement uncertainty; σ_W is an additional white noise jitter term; σ_C is a calibration noise term; δ_k,l is equal to one if measurements k and l are taken within the same night, and zero otherwise; σ_R and τ_R define a correlated term to model contributions due to granulation (Cegla 2019) or instrumental effects, as defined in Hara et al. (2020); σ_act, τ_act, and P_act are the hyper-parameters of a quasi-periodic covariance term to model stellar activity.

For σ_W, σ_r, σ_qp, we use the grid of values 0.5, 1., 1.5 m/s. σ_C is fixed to 0.5 m/s. τ = 3 or 6 days, P_act is fixed to 31.4 days based on the analysis of activity indicators and τ_act = 60 days. We try every combination of these values and find that the model with highest cross-validation has σ_W = 1.5 m/s, σ_R = 0.5 m/s, σ_QP = 1 m/l, τ_R = 6 days. We use a free offset for each of the three datasets: pre- and post-2015 HARPS and CARMENES. The corresponding ℓ₁ periodogram is shown in Figure E.1. We find three peaks at 15.1, 140, and 15.8 days in decreasing order. Their corresponding FAP is 1.9 · 10⁻⁵, 2 · 10⁻³ and 0.97. The first two signals are statistically significant, and the second is the orbital period of the candidate planet Gl 514 b. It is interesting to note that the peak at P = 15.8 d has a very high FAP, contrary to the result of the MLP. The ℓ₁ method searches for several signals simultaneously and accounts for correlated noise at the stellar rotation period. As a result, signals that could appear dominant in a regular periodogram, either because of stellar activity or aliasing, appear damped in ℓ₁.

Fig. E.1 ℓ1 periodogram with free offsets and covariance models with the highest cross-validation score.

To further examine their properties, we apply the method of Hara et al. (2022), which is designed to test whether the signals have consistent amplitude, phase, and frequencies over time. We compute the apodized sine periodogram (ASP), which is defined as follows. We compare the χ² of a linear base model H, and the χ² of a model K(ω, t₀, τ) defined as the linear model of H plus an apodized sinusoid . (A cos ωt + B sin ωt). The ASP is defined as the χ² difference between the two models (E.2)

Like in the ℓ₁ periodogram, our base model includes the offsets of the three datasets. For the apodized sinusoid, we consider a grid of timescales τ = 10 · T_obs, T_obs/3, T_obs/9 and T_obs/27, where T_obs is the total observation timespan of the data. In Fig. E.2, we represent z(ω, t₀,τ) as defined in Eq. (E.2), maximised over t₀. The dominating peak has a period of 15.2 days. We want to test whether the signal is statistically compatible with a constant one (i.e. with τ = 10 · T_obs). Using t_(τ,ω) to denote the value of t₀ maximising the value of the periodogram (E.2) for a given ω and τ, we compute the distribution (E.3)

with the hypothesis that model K(ω, t_(τ,ω), τ, A*, B*) is correct, where the fitted cosine and sine amplitudes A*, B* are obtained by fitting model K to the data. D_z can easily be expressed as a generalised χ² distribution, with mean and variance given by an analytical expression (Hara et al. 2022). In Fig. E.2, horizontal dashed lines correspond to the values of z at the frequency where the maximum is attained. For each value of τ, we assume that the data contain a signal at frequency 15.2 days and with timescale τ. The points with error bars correspond to the expected value of the peaks in the ASP and their standard deviation, assuming that the timescale τ in abscissa is the correct one. These are computed as z(ω, t_(τ,ω),τ) − E(D_Z), where E(D_Z) is the expectancy of D_z as defined in Eq. (E.3), and as the variance of D_z. The more each expected value (with error bars) is closer to the level defined by the dashed horizontal line of the same colour, the more the value of τ in the abscissa is the best representation for the timescale associated to the signal of period P. From Fig. E.2, we see that assuming τ = 10 · T_obs as the correct timescale, the expected position of the ASP peaks of the other timescales is more than three sigma away from where they are observed. Their expected positions are in agreement with the observed levels for τ = T_obs/27, meaning that the period P = 15.2 days is not stable over time.

Fig. E.2 Left panel: ASP corresponding to different apodization timescales: τ = 10 · Tobs, Tobs/3, Tobs/9 and Tobs/27, respectively in blue, orange, green, and red. Middle panel: Zoom onto the maximum peak. Right panel: Statistical tests as defined in Hara et al. (2022).

Fig. E.3 As in Fig. E.2. The base model now includes the signal corresponding to the highest peak of the first ASP (Fig. E.2).

We pre-whitened the original dataset by including the 15.2-day signal in the null hypothesis/base model H, and we repeated the calculation of the ASP. We obtain the periodogram in Fig. E.3, which shows a peak occurring at a period of 140 days. The statistical test in the right panel of the figure shows that this signal is compatible with having a timescale τ = T_obs/3 or greater. This means that the signal can be credibly approximated as a sinusoid (we note that this method does not take the orbital eccentricity into account), meaning that it is compatible with a signal induced by a planet.

Appendix F Results relative to RV models including two Keplerians

In this section, we show plots related to the tests performed on the HARPS_TERRA+CARMENES RVs discussed in Sect. 6.1.2, where we examined the hypothesis of a two-planet system.

Fig. F.1 Corner plots showing the posterior distributions of all the free (hyper-)parameters of the two-planet models described in Sect. 6.1.2. Upper panel: QP kernel; Middle panel: QPC kernel; Lower panel: dSHO/rotational kernel.

References

All Tables

All Figures

	Fig. 1 Time-series of the HIRES (orange dots), HARPSTERRA (red and black dots), and CARMENES-VIS (green dots) RVs analysed in this work. The average has been subtracted from each dataset.
In the text

	Fig. 2 Time series and GLS periodogram of the chromospheric activity S-index derived from HIRES spectra (first and second row). The levels of FAPs are indicated as horizontal dashed lines (green: FAP 0.1%; yellow: 1%; red: 10%), and have been determined through a bootstrap analysis. Third row: Window function of the data.
In the text

Fig. 3 Time-series and GLS periodograms for activity diagnostics derived from HARPS spectra. Rows 1–4: time series (left panels) and periodograms (right panels) of the CCF FWHM, and activity indexes derived from the spectral lines CaII H&K, Hα and NaI. For the FWHM, we show the periodogram of the residuals, after removing the long-term trend clearly seen in the time series. The levels of FAPs are indicated as horizontal dashed lines (green: FAP 0.1%; yellow: 1%; red: 10%), and are determined through a bootstrap analysis. Last row: Window function of the data.

In the text

Fig. 4 Time-series of activity diagnostics extracted from CARMENES-VIS spectra of Gl 514 (left panels) and their GLS periodograms (right panels). Except for the CRX index, the other periodograms were calculated after a pre-whitening by removing the clear long-term trend visible in the time-series. The levels of FAPs, determined through a bootstrap analysis, are indicated as horizontal dashed lines (green: FAP 0.1%; yellow: 1%; red: 10%. Last panel: Window function of the CARMENES-VIS data.

In the text

	Fig. 5 Maximum-likelihood periodogram of the HIRES RVs. The periods corresponding to the main peaks are indicated, as are the analytical FAP levels.
In the text

	Fig. 6 Maximum-likelihood periodogram of the HARPS TERRA (first panel), NAIRA (second panel), and Trifonov et al. (2020; third panel) RVs (pre- and post-2015 combined dataset). The periods corresponding to the four highest peaks are indicated, as are the analytical FAP levels.
In the text

	Fig. 7 Upper panel: maximum-likelihood periodogram of the original CARMENES-VIS RVs. The periods corresponding to the highest peaks are indicated, as are the analytical FAP levels. Lower panel: maximum-likelihood periodogram of the pre-whithened RVs after removing the 354-day signal.
In the text

Fig. 8 Maximum-likelihood periodograms of the joint HARPS (preand post-2015) and CARMENES-VIS RVs (panel 1: HARPSTERRA; panel 2: HARPSNAIRA; panel 3: HARPSTrifonov et al.). The periods corresponding to the highest peaks are indicated, as are the analytical FAP levels. Offsets have been subtracted from the HARPS and CARMENES RVs, as derived by the MLP procedure. Lower panel: window function of the combined dataset.

In the text

	Fig. 9 Upper panel: maximum-likelihood periodogram of the joint HIRES, HARPS (TERRA; pre- and post-2015), and CARMENES-VIS RVs. The periods corresponding to the highest peaks are indicated, as are the analytical FAP levels. Lower panel: Window function of the combined dataset.
In the text

	Fig. 10 Spectroscopic orbit of Gl 514 b (blue solid line) based on the best-fit solution shown in Table 5 (case eb ≠ 0). The error bars of the measurements include uncorrelated jitter terms added in quadrature to the formal RV uncertainties. Cyan dots correspond to RV data averaged over 15 phase bins.
In the text

Fig. 11 Posterior distributions of the semi-amplitude Kb (top panel), orbital period Pb (middle panel), and eccentricity eb (bottom panel) of the candidate planet Gl514 b derived from a GP QPC analysis of the HARPSterra+CARMENES VIS RVs after removing an increasing number of CARMENES measurements, beginning from the last epoch. Vertical lines indicate the 16th, 84th (dashed), and 50th (solid) percentiles for the whole HARPSTERRA+CARMENES dataset. Here, ∆ ln Z denotes the Bayesian evidence difference between a model with and without a Keplerian included.

In the text

Fig. 12 Maximum likelihood periodogram of the HIRES, HARPS, and CARMENES RV residuals after removing only the best-fit signal of the candidate planet with an orbital period of ∼140 days from the original data (second column of Table 5). The best-fit sinusoid with period P = 357.76 days has been subtracted from the CARMENES-VIS data. Stellar activity has not been filtered out, because the best-fit model for the activity signal is linked only on RVs measured with HARPS and CARMENES.

In the text

Fig. 13 TESS short-cadence light curve of Gl 514. Panels a and b show the normalised light curves before and after the subtraction of the model, respectively; the model represented in red in panel a was calculated on a grid of knots spaced every 24h. Time is expressed as the TESS Barycenter-corrected Julian Day (BTJD = BJD - 2457 000). Panel c shows the BLS periodogram. The red triangle indicates the period associated with the highest S/N.

In the text

	Fig. 14 Posterior distributions for the conjuction times Tconj,b of Gl 514 b that are close to TESS Sectors 23 (March-April 2020) and 50 (March-April 2022). The posteriors are derived from that obtained for our best-fit QPC model. The vertical lines mark the time-span of TESS observations.
In the text

	Fig. 15 Detectability curve for companions of given mass and orbital distance around Gl 514 based on the proper motion anomaly technique described in Kervella et al. (2019).
In the text

Fig. 16 Schematic representation of a few possible orbits for planet Gl 514 b (blue ellipses) randomly selected taking into account the parameter values and their errors of our adopted solution (Table 5, second column). The best-fit orbit is shown with a dashed black line. The red star identifies one of the foci of the ellipses occupied by the host star Gl 514 b, and the green annulus corresponds to the conservative habitable zone.

In the text

	Fig. 17 Gl 514 b in the context of planets with mass or minimum mass <10 M⊕ that have been detected to date within the HZ of M dwarfs. The insolation flux of Gl 514 b is the temporal average orbit value.
In the text

	Fig. 18 Expected flux density for auroral radio emission arising from star-planet interaction in the system Gl 514 as a function of orbital distance. The interaction is expected to be in the subAlfvénic regime (i.e. MA = vrel/vAlfv ≤ 1; top panel) at distances equal to the semi-major axis and periapsis (vertical dashed line).
In the text

	Fig. C.1 Posterior distributions of the free (hyper)parameters of our assumed best-fit model, including a Keplerian for planet b and a GP QPC correlated activity signal fitted to the RVs obtained from HARPS (TERRA dataset) and CARMENES VIS spectra.
In the text

	Fig. C.2 Phase-folded RVs of HARPSterra and CARMENES-VIS showing the spectroscopic orbit of Gl 514 b for the circular case. The blue solid line is the best-fit solution using a QPC kernel to model the activity. The error bars include the uncorrelated jitter terms added in quadrature to the formal RV uncertainties. Cyan dots correspond to RV data averaged over 15 phase bins.
In the text

	Fig. C.3 Sections of the residuals of the whole set of HARPSterra+CARMENES RVs (corresponding to the model shown in the second column of Table 5, showing the QPC correlated signal mostly related to variations in stellar activity. The error bars include the uncorrelated jitter terms added in quadrature to the formal RV uncertainties.
In the text

	Fig. C.4 GLS periodograms of the HARPSterra+CARMENES RV residuals, after removing the best-fit models shown in the second columns of Table C.1–C.2. For each kernel, periodograms of data with and without the correlated activity signal included are shown. The periodograms were calculated with the uncorrelated jitter terms added in quadrature to the formal RV uncertainties.
In the text

	Fig. E.1 ℓ1 periodogram with free offsets and covariance models with the highest cross-validation score.
In the text

	Fig. E.2 Left panel: ASP corresponding to different apodization timescales: τ = 10 · Tobs, Tobs/3, Tobs/9 and Tobs/27, respectively in blue, orange, green, and red. Middle panel: Zoom onto the maximum peak. Right panel: Statistical tests as defined in Hara et al. (2022).
In the text

	Fig. E.3 As in Fig. E.2. The base model now includes the signal corresponding to the highest peak of the first ASP (Fig. E.2).
In the text

	Fig. F.1 Corner plots showing the posterior distributions of all the free (hyper-)parameters of the two-planet models described in Sect. 6.1.2. Upper panel: QP kernel; Middle panel: QPC kernel; Lower panel: dSHO/rotational kernel.
In the text