# A quarter century of spectroscopic monitoring of the nearby M dwarf Gl 514. A super-Earth on an eccentric orbit moving in and out of the habitable zone

M. Damasso, M. Perger, J. M. Almenara, D. Nardiello, M. Pérez-Torres, A. Sozzetti, N. C. Hara, A. Quirrenbach, X. Bonfils, M. R. Zapatero Osorio, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

M. Damasso, M. Perger, J. M. Almenara, D. Nardiello, M. Pérez-Torres, et al. A quarter century of spectroscopic monitoring of the nearby M dwarf Gl 514. A super-Earth on an eccentric orbit moving in and out of the habitable zone. Astronomy and Astrophysics - A&A, 2022, 666, 10.1051/0004-6361/202243522. insu-03860283

HAL Id: insu-03860283

https://insu.hal.science/insu-03860283

Submitted on 19 Nov 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# A quarter century of spectroscopic monitoring of the nearby M dwarf Gl 514

## A super-Earth on an eccentric orbit moving in and out of the habitable zone ⋆, ⋆⋆

M. Damasso<fnref n="1" /> Ⓓ, M. Perger²,³ Ⓓ, J. M. Almenara⁴ Ⓓ, D. Nardiello⁵,⁶ Ⓓ, M. Pérez-Torres⁷ Ⓓ, A. Sozzetti<fnref n="1" /> Ⓓ, N. C. Hara⁸ Ⓓ, A. Quirrenbach⁹, X. Bonfils⁴ Ⓓ, M. R. Zapatero Osorio¹⁰, N. Astudillo-Defru¹¹, J. I. González Hernández¹²,¹³ Ⓓ, A. Suárez Mascareño¹³,¹² Ⓓ, P. J. Amado⁷ Ⓓ, T. Forveille⁴, J. Lillo-Box¹⁴ Ⓓ, Y. Alibert¹⁵, J. A. Caballero¹⁴ Ⓓ, C. Cifuentes¹⁴ Ⓓ, X. Delfosse⁴, P. Figueira¹⁶,¹⁷ Ⓓ, D. Galadí-Enríquez¹⁸ Ⓓ, A. P. Hatzes¹⁹ Ⓓ, Th. Henning²⁰, A. Kaminski⁹ Ⓓ, M. Mayor⁸, F. Murgas¹³,¹² Ⓓ, D. Montes²¹, M. Pinamonti<fnref n="1" /> Ⓓ, A. Reiners²² Ⓓ, I. Ribas²,³ Ⓓ, V. J. S. Béjar¹³,¹², A. Schweitzer²³ Ⓓ, and M. Zechmeister²² Ⓓ

*(Affiliations can be found after the references)*

Received 10 March 2022 / Accepted 13 April 2022

## ABSTRACT

**Context.** Statistical analyses based on *Kepler* data show that most of the early-type M dwarfs host multi-planet systems consisting of Earth- to sub-Neptune-sized planets with orbital periods of up to ~250 days, and that at least one such planet is likely located within the habitable zone. M dwarfs are therefore primary targets to search for potentially habitable planets in the solar neighbourhood.

**Aims.** We investigated the presence of planetary companions around the nearby (7.6 pc) and bright (V = 9 mag) early-type M dwarf Gl 514, analysing 540 radial velocities collected over nearly 25 yr with the HIRES, HARPS, and CARMENES spectrographs.

**Methods.** The data are affected by time-correlated signals at the level of 2–3 m s⁻<fnref n="1" /> due to stellar activity, which we filtered out, testing three different models based on Gaussian process regression. As a sanity cross-check, we repeated the analyses using HARPS radial velocities extracted with three different algorithms. We used HIRES radial velocities and *Hipparcos-Gaia* astrometry to put constraints on the presence of long-period companions, and we analysed TESS photometric data.

**Results.** We find strong evidence that Gl 514 hosts a super-Earth on a likely eccentric orbit, residing in the conservative habitable zone for nearly 34% of its orbital period. The planet Gl 514 b has minimum mass $m_b \sin i_b = 5.2 \pm 0.9 M_⊕$, orbital period $P_b = 140.43 \pm 0.41$ days, and eccentricity $e_b = 0.45^{+0.15}_{-0.14}$. No evidence for transits is found in the TESS light curve. There is no evidence for a longer period companion in the radial velocities and, based on astrometry, we can rule out a ~0.2 $M_{Jup}$ planet at a distance of ~3–10 astronomical units, and massive giant planets and brown dwarfs out to several tens of astronomical units. We discuss the possible presence of a second low-mass companion at a shorter distance from the host than Gl 514 b.

**Conclusions.** Gl 514 b represents an interesting science case for studying the habitability of planets on eccentric orbits. We advocate for additional spectroscopic follow-up to get more accurate and precise planetary parameters. Further follow-up is also needed to investigate the presence of additional planetary signals of less than 1 m s⁻<fnref n="1" />.

**Key words.** techniques: radial velocities – techniques: photometric – planetary systems – stars: activity – stars: individual: Gl 514

## 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<fnref n="1" />, complemented with the archive of confirmed

<footnote>* Tables A.1–A.4 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/666/A187</footnote>

<footnote>** Based on HIRES, CARMENES and HARPS observations, the latter carried out during programs 072.C-0488, 183.C-0437, and 191.C-0873.</footnote>

<footnote n="1">¹ https://exoplanetarchive.ipac.caltech.edu/. Queried on 7 April 2022.</footnote>
potentially habitable low-mass planets maintained by the Planetary Habitability Laboratory (PHL)<fnref n="2" />, 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_{\text{eff}}$ in the range 3650–3900 K. We found three transiting planets detected by *Kepler*, namely *Kepler-186 f* (Quintana et al. 2014), *Kepler-705 b*, and *Kepler-1229 b* (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.6 S_⊕$) 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-Mascareño 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_{\star,\text{rot}} = 30.8 \pm 0.3$ days and $30.3 \pm 0.2$ days, which were derived from the analysis of the $H\alpha$ and CaII infrared triplet lines, respectively. As we discuss in

<table><caption>**Table 1.** Astrometric, photometric, and spectroscopically derived stellar properties of Gl 514.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Refs.</th></tr></thead><tbody><tr><td><strong>Astrometry</strong></td><td></td><td></td></tr><tr><td>α (J2000)</td><td>13h 29m 59.79s</td><td>[1,2,3]</td></tr><tr><td>δ (J2000)</td><td>+10° 22′37.8″</td><td>[1,2,3]</td></tr><tr><td>μ<sub>α</sub> · cos δ (mas yr<sup>-1</sup>)</td><td>1127.341 ± 0.027</td><td>[1,2,3]</td></tr><tr><td>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>-1073.888 ± 0.014</td><td>[1,2,3]</td></tr><tr><td>ω (mas)</td><td>131.101 ± 0.027</td><td>[1,2,3]</td></tr><tr><td>d (pc)</td><td>7.618 ± 0.003</td><td>[4]</td></tr><tr><td><strong>Photometry</strong></td><td></td><td></td></tr><tr><td>U</td><td>11.7</td><td>[5]</td></tr><tr><td>B</td><td>10.5</td><td>[5]</td></tr><tr><td>V</td><td>9.0</td><td>[5]</td></tr><tr><td>G</td><td>8.2053 ± 0.0027</td><td>[1,2,3]</td></tr><tr><td>J</td><td>5.902 ± 0.018</td><td>[6]</td></tr><tr><td>H</td><td>5.300 ± 0.033</td><td>[6]</td></tr><tr><td>K<sub>s</sub></td><td>5.036 ± 0.027</td><td>[6]</td></tr><tr><td><strong>Stellar parameters</strong></td><td></td><td></td></tr><tr><td>T<sub>eff</sub> (K)</td><td>3728 ± 68<br>3714 ± 51</td><td>[8]<br>[9]</td></tr><tr><td>log g (dex)</td><td>4.76 ± 0.04<br>4.81 ± 0.07</td><td>[8]<br>[9]</td></tr><tr><td>[Fe/H] (dex)</td><td>-0.14 ± 0.09<br>-0.15 ± 0.16</td><td>[8]<br>[9]</td></tr><tr><td>M<sub>★</sub> (M<sub>☉</sub>)</td><td>0.510 ± 0.051</td><td>[8]</td></tr><tr><td>R<sub>★</sub> (R<sub>☉</sub>)</td><td>0.500 ± 0.047</td><td>[8]</td></tr><tr><td>L<sub>★</sub> (L<sub>☉</sub>)</td><td>0.043 ± 0.009</td><td>[8]</td></tr><tr><td>P<sub>★,rot</sub> (days)</td><td>28.0 ± 2.9<br>30.8 ± 0.3 (Hα)<br>30.3 ± 0.2 (CaII infrared triplet)<br>30 (CaII H&K)<br>30.6 ± 0.3 (GP fit of RVs)</td><td>[10]<br>[11]<br>[11]<br>[12]<br>[7]</td></tr><tr><td>⟨log R'<sub>HK</sub>⟩</td><td>-5.10 ± 0.06</td><td>[10]</td></tr><tr><td>(U,V,W) (km s<sup>-1</sup>)</td><td>(57.42 ± 0.03, -8.11 ± 0.02,<br>-3.54 ± 0.09)</td><td>[7]</td></tr><tr><td>Kinematical age (Gyr)</td><td>> 0.8</td><td>[7]</td></tr><tr><td>Conservative HZ<sup>(a)</sup> (au)</td><td>[0.207,0.411]</td><td>[13]</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>For planet mass $m_p = 5 M_⊕$

**References.** [1] Gaia Collaboration (2016); [2] Gaia Collaboration (2021); [3] Lindegren et al. (2021); [4] Bailer-Jones et al. (2018b); [5] Koen et al. (2010); [6] Cutri et al. (2003); [7] This work; [8] Maldonado et al. (2015); [9] Schweitzer et al. (2019); [10] Suárez-Mascareño et al. (2015); [11] Fuhrmeister et al. (2019); [12] Astudillo-Defru et al. (2017a); [13] Kopparapu et al. (2013, 2014).

Sects. 4 and 6.1, our derived values for $P_{\star,\text{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

<footnote n="2">² http://phl.upr.edu/projects/habitable-exoplanets-catalog (updated to 6 Dec 2021).</footnote>
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 cover the time span between 7 April 1997 and 14 December 2014. We used the publicly available data derived by Tal-Or et al. (2019)<fnref n="3" />.

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<fnref n="4" />. 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 $\lambda \ge 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 ($\gamma$) and an uncorrelated jitter term ($\sigma_{jit}$) as free parameters for each dataset separately. We also exploited two other RV datasets extracted with alternative template-matching pipelines as a safety 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 $\sigma_{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)<fnref n="5" />.

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 $\sigma_{RV}$ uncertainty of 8.7 m s⁻¹, and are therefore not particularly useful for detecting and characterising RV signals with a semi-amplitude of few m s⁻¹, which are those we are interested in. Originally, CARMENES collected 321 spectra of Gl 514, 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 SERVAL (Zechmeister et al. 2018). During a few nights, CARMENES observed the star more than once. We binned the corresponding RVs on a nightly basis<fnref n="6" />, and then discarded 14 measurements with $\sigma_{RV} > 3$ m s⁻¹ (nearly twice the mean value of $\sigma_{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 $\sigma_{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

<footnote n="3">³ https://cdsarc.unistra.fr/viz-bin/cat/J/MNRAS/484/L8</footnote>

<footnote n="4">⁴ http://archive.eso.org/wdb/wdb/adp/phase3_main/form</footnote>

<footnote n="5">⁵ Also available at https://www2.mpio-hd.mpg.de/homes/trifonov/HARPS_RVBank.html</footnote>

<footnote n="6">⁶ As shown below, we discuss signals with periods of several days, including activity and planets, therefore the binning does not affect the analysis.</footnote>
<table><caption>**Table 2.** General properties of the RV time-series analysed in this work.</caption><thead><tr><th>Instrument</th><th>Time span (days)</th><th>No. RVs</th><th>Pipeline</th><th>RMS (m s<sup>-1</sup>)</th><th>&sigma;<sub>RV</sub> (m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>HIRES</td><td>6095</td><td>104</td><td>Tal-Or et al. (2019)</td><td>3.8</td><td>1.6</td></tr><tr><td rowspan="3">HARPS<sub>pre-2015</sub></td><td rowspan="3">3963</td><td rowspan="3">142</td><td>TERRA</td><td>2.8</td><td>0.7</td></tr><tr><td>NAIRA</td><td>2.9</td><td>0.6</td></tr><tr><td>Trifonov et al. (2020), SERVAL</td><td>3.0</td><td>0.9</td></tr><tr><td rowspan="3">HARPS<sub>post-2015</sub></td><td rowspan="3">266</td><td rowspan="3">20</td><td>TERRA</td><td>2.5</td><td>0.7</td></tr><tr><td>NAIRA</td><td>2.8</td><td>0.6</td></tr><tr><td>Trifonov et al. (2020), SERVAL</td><td>2.6</td><td>0.9</td></tr><tr><td>CARMENES-VIS</td><td>1909</td><td>274</td><td>SERVAL</td><td>2.6</td><td>1.6</td></tr></tbody></table>

**Notes.** Total time spans – HARPS: 4352 days; HARPS+CARMENES: 6154 days; HIRES+HARPS+CARMENES: 8761 days.

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Time-series of the HIRES (orange dots), HARPS<sub>TERRA</sub> (red and black dots), and CARMENES-VIS (green dots) RVs analysed in this work. The average has been subtracted from each dataset.</figcaption></figure>

stray-light contamination, which in turn badly affects the overall quality of the data.

## 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_{\star,rot}$ and its first harmonic (see Table 1).

### 4.2. HARPS

We derived activity diagnostics based on the CaII H&K, $H\alpha$, 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_{\star,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_{\star,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
<table><caption>**Table 3.** Peaks with the highest power observed in the GLS periodograms of the spectroscopic chromospheric activity indicators extracted from HIRES, HARPS, and CARMENES-VIS spectra (see Figs. 2, 3, and 4).</caption><thead><tr><th>Activity diagnostic</th><th>Frequency (days<sup>-1</sup>)</th><th>Period (days)</th><th>Note</th></tr></thead><tbody><tr><td><strong>HIRES</strong></td><td></td><td></td><td></td></tr><tr><td>S-index</td><td>0.00136<br>0.00159<br>0.0325<br>0.0673</td><td>734<br>629<br>30.7<br>14.9</td><td>FAP = 0.01%<br>FAP = 0.01%; alias of the main peak frequency (f<sub>alias</sub> = 0.000229 day<sup>-1</sup>, or P<sub>alias</sub> = 4361 days)<br>FAP ~ 0.1%<br>FAP ~ 1%</td></tr><tr><td><strong>HARPS</strong></td><td></td><td></td><td></td></tr><tr><td>CCF FWHM<sup>(a)</sup></td><td>0.00126<br>0.0304<br>0.00402<br>0.000705<br>0.0332</td><td>796<br>32.9<br>249<br>1417<br>30.1</td><td>FAP = 0.01%<br>FAP = 0.01%<br>FAP = 0.01%<br>FAP = 0.01%<br>FAP = 0.01%</td></tr><tr><td>CaII H&amp;K index</td><td>0.001022<br>0.03049<br>0.0333<br>0.00139</td><td>978<br>32.8<br>30.0<br>717</td><td>FAP = 0.01%<br>FAP = 0.01%<br>FAP = 0.01%; 1-yr alias of the previous frequency<br>FAP ~ 0.1%</td></tr><tr><td>Hα index</td><td>0.002038<br>0.00256<br>0.0007<br>0.0303</td><td>491.8<br>390.3<br>1427<br>33</td><td>FAP = 0.01%<br>FAP = 0.01%; alias of the main peak frequency (f<sub>alias</sub> = 0.00303 day<sup>-1</sup>, or P<sub>alias</sub> = 329.7 days)<br>FAP = 0.1%; 1-yr alias of the main peak frequency<br>FAP = 0.01%<br>FAP = 0.01%<br>FAP &lt; 0.1%</td></tr><tr><td><strong>CARMENES-VIS<sup>(b)</sup></strong></td><td></td><td></td><td></td></tr><tr><td>Ca-IRT index</td><td>0.00157<br>0.00165<br>0.0625</td><td>636<br>606<br>16</td><td>FAP = 0.01%<br>FAP = 0.1%<br>1% &lt; FAP &lt; 10%<br>1% &lt; FAP &lt; 10%</td></tr><tr><td>dLW</td><td>0.00175<br>0.0294</td><td>570<br>34</td><td>FAP = 0.01%<br>FAP = 1%</td></tr><tr><td>CRX</td><td>0.0024<br>0.00565<br>0.0164</td><td>410<br>177<br>61</td><td>FAP = 0.01%<br>likely related to the data sampling<br>1% &lt; FAP &lt; 10%</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>The periodogram was calculated on the residuals after removing a long-term trend. <sup>(b)</sup>The peak frequencies are calculated from pre-whitened data, except for the CRX index.

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<fnref n="7" />. Differently from the commonly used generalised Lomb-Scargle periodogram (GLS; Zechmeister et al. 2009), for each tested frequency, the MLP shows the difference $d(ln L)$

<footnote n="7">⁷ see https://ascl.net/1906.004</footnote>
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

**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.

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.

## 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_{*,\text{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_{*,\text{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.

<figure><img src="image_7.png" /></figure>

**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<fnref n="8" />. Signals

<footnote n="8">⁸ Since no significant signals have been detected at high frequencies, we show the periodograms up to 0.1 day⁻¹ for more clarity.</footnote>
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

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.

## 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

<figure><img src="image_9.png" /><figcaption>**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-whitened RVs after removing the 354-day signal.</figcaption></figure>

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.

## 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 HARPS_TERRA 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
<figure><img src="image_10.png" /></figure>

**Fig. 8.** Maximum-likelihood periodograms of the joint HARPS (pre- and post-2015) and CARMENES-VIS RVs (panel 1: HARPS<sub>TERRA</sub>; panel 2: HARPS<sub>NAIRA</sub>; panel 3: HARPS<sub>Trifonov et al.</sub>). 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.

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

<figure><img src="image_11.png" /></figure>

**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.

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.

## 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_{\star,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 $\Delta \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 $\omega_{\star,b}$, using the parametrization $\sqrt{e_b} \cos \omega_{\star,b}$ and $\sqrt{e_b} \sin \omega_{\star,b}$. 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

<table><caption>**Table 4.** Priors used for modelling the HIRES, HARPS, and CARMENES VIS radial velocities.</caption><thead><tr><th>Parameter</th><th>Prior</th></tr></thead><tbody><tr><td colspan="2"><strong>Stellar activity – QP kernel</strong></td></tr><tr><td><em>h</em> (m s<sup>-1</sup>)</td><td><em>U</em>(0, 10)<sup>(a)</sup></td></tr><tr><td><em>θ</em> (d)</td><td><em>U</em>(20, 50)</td></tr><tr><td><em>λ</em><sub>QP</sub> (d)</td><td><em>U</em>(0, 1000)</td></tr><tr><td><em>w</em></td><td><em>U</em>(0, 1)</td></tr><tr><td colspan="2"><strong>Stellar activity – QPC kernel</strong></td></tr><tr><td><em>h</em><sub>1</sub> (m s<sup>-1</sup>)</td><td><em>U</em>(0, 10)<sup>(a)</sup></td></tr><tr><td><em>h</em><sub>2</sub> (m s<sup>-1</sup>)</td><td><em>U</em>(0, 10)<sup>(a)</sup></td></tr><tr><td><em>θ</em> (d)</td><td><em>U</em>(20, 50)</td></tr><tr><td><em>λ</em><sub>QPC</sub> (d)</td><td><em>U</em>(0, 1000)</td></tr><tr><td><em>w</em></td><td><em>U</em>(0, 1)</td></tr><tr><td colspan="2"><strong>Stellar activity – dSHO kernel</strong></td></tr><tr><td>log <em>A</em></td><td><em>U</em>(0.05,10)<sup>(a)</sup></td></tr><tr><td><em>θ</em> (d)</td><td><em>U</em>(20, 50)</td></tr><tr><td>log <em>Q</em><sub>0</sub></td><td><em>U</em>(-10, 10)</td></tr><tr><td>log Δ<em>Q</em></td><td><em>U</em>(-10, 10)</td></tr><tr><td><em>f</em></td><td><em>U</em>(0, 10)</td></tr><tr><td colspan="2"><strong>First Keplerian</strong></td></tr><tr><td><em>K</em><sub>b</sub> (m s<sup>-1</sup>)</td><td><em>U</em>(0, 5)</td></tr><tr><td><em>P</em><sub>b</sub> (d)</td><td><em>U</em>(0,200)<sup>(b)</sup><br><em>U</em>(100,200)<sup>(c)</sup></td></tr><tr><td><em>T</em><sub>conj,b</sub> (BJD-2 450 000)</td><td><em>U</em>(8500,8750)</td></tr><tr><td>√<em>e</em><sub>b</sub> cos ω<sub>★,b</sub></td><td><em>U</em>(-1, 1)</td></tr><tr><td>√<em>e</em><sub>b</sub> sin ω<sub>★,b</sub></td><td><em>U</em>(-1, 1)</td></tr><tr><td colspan="2"><strong>Second Keplerian (innermost)</strong></td></tr><tr><td><em>K</em><sub>c</sub> (m s<sup>-1</sup>)</td><td><em>U</em>(0, 5)</td></tr><tr><td><em>P</em><sub>c</sub> (d)</td><td><em>U</em>(0, 100)</td></tr><tr><td><em>T</em><sub>conj,c</sub> (BJD-2 450 000)</td><td><em>U</em>(8500,8650)</td></tr><tr><td>√<em>e</em><sub>c</sub> cos ω<sub>★,c</sub></td><td><em>U</em>(-1, 1)</td></tr><tr><td>√<em>e</em><sub>c</sub> sin ω<sub>★,c</sub></td><td><em>U</em>(-1, 1)</td></tr><tr><td colspan="2"><strong>Third Keplerian (outermost)<sup>(d)</sup></strong></td></tr><tr><td><em>K</em><sub>d</sub> (m s<sup>-1</sup>)</td><td><em>U</em>(0, 5)</td></tr><tr><td><em>P</em><sub>d</sub> (d)</td><td><em>U</em>(200,4350)</td></tr><tr><td><em>T</em><sub>conj,d</sub> (BJD-2 450 000)</td><td><em>U</em>(4500,9000)</td></tr><tr><td>√<em>e</em><sub>d</sub> cos ω<sub>★,d</sub></td><td><em>U</em>(-1, 1)</td></tr><tr><td>√<em>e</em><sub>d</sub> sin ω<sub>★,d</sub></td><td><em>U</em>(-1, 1)</td></tr><tr><td colspan="2"><strong>Additional sinusoid<sup>(e)</sup></strong></td></tr><tr><td><em>K</em><sub>365-d</sub> (m s<sup>-1</sup>)</td><td><em>U</em>(0, 10)</td></tr><tr><td><em>P</em><sub>365-d</sub> (d)</td><td><em>U</em>(340,370)</td></tr><tr><td><em>T</em><sub>0,365-d</sub> (BJD-2 450 000)</td><td><em>U</em>(8500,8900)</td></tr><tr><td colspan="2"><strong>Instrument-related</strong></td></tr><tr><td>γ (m s<sup>-1</sup>)</td><td><em>U</em>(-20,20)<sup>(a)</sup></td></tr><tr><td>σ<sub>jit</sub> (m s<sup>-1</sup>)</td><td><em>U</em>(0, 10)<sup>(a)</sup></td></tr></tbody></table>

**Notes.** (a) The same prior was used for the corresponding parameter related to HIRES, HARPS (pre- and post-2015), and CARMENES-VIS data. (b) Prior used when modelling only one Keplerian. (c) Prior used when including a second Keplerian to model planetary signals with $P < 200$ days. (d) Priors for a planetary signal with period longer than 200 days, which we modelled combining HIRES+HARPS+CARMENES RVs. (e) Priors used to model the yearly signal only seen in the CARMENES-VIS RVs.
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 514 b. 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<sub>1p</sub>-ln Z<sub>0p</sub> > 5); the planetary-like signal is detected in all the three cases with K<sub>b</sub> ~ 1.2 m s<sup>-1</sup>(5.5-6σ significance), an orbital period of P<sub>b</sub> ~ 140 days, determined with an error bar lower than 1 day, and a 3σ significant eccentricity e<sub>b</sub> ~ 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<sup>9</sup>. 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. 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 corre-
lated 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 simi-
lar 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 (θ = 37.0+5.4-14.4 and λQPC ~ 47 days;
Affer et al. 2019), K2-3 (θ = 40.4+1.1-1.9 and λQPC ~ 80 days;
Damasso et al. 2018), GJ 3998 (θ = 31.8+0.6-0.5 and λQPC ~ 69 days;
Affer et al. 2016), or Gl 15A (θ = 46.7+4.8-4.3 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 HARPS<sub>TERRA</sub>+CARMENES RVs
(keeping all the priors unchanged) after removing the last
100, 150, and 200 CARMENES measurements. This cor-
responds to a reduction of 23%, 34.4%, and 45.9% of the
total HARPS<sub>TERRA</sub>+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

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Spectroscopic orbit of Gl 514 b (blue solid line) based on the best-fit solution shown in Table 5 (case *e*<sub>b</sub> ≠ 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.</figcaption></figure>

semi-amplitude K<sub>b</sub>, orbital period P<sub>b</sub>, and eccentricity e<sub>b</sub> 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 HARPS<sub>TERRA</sub>+CARMENES RVs. The
posteriors of K<sub>b</sub> are all in agreement within the uncertainties,
and those of P<sub>b</sub> become narrower and move closer to that of
the whole dataset with the increasing number of RVs. The
posteriors for e<sub>b</sub> 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.

6.1.2. Testing the two-planet model

Given the large number of HARPS<sub>TERRA</sub>+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.4 m s<sup>-1</sup>, 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<sub>*,rot</sub>. As signals with periods
longer than P<sub>*,rot</sub> 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 orig-
nal dataset, are not particularly informative when searching for
signals with P > P<sub>*,rot</sub>. From this analysis, we expect to get
at best hints for a possible sub m s<sup>-1</sup> 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

<footnote n="9">⁹ This result is generally confirmed when using the two alternative HARPS RV datasets (see Appendix D).</footnote>
<table><caption>**Table 5.** Best-fit parameters obtained for the model with the QPC kernel applied to HARPS<sub>TERRA</sub> and CARMENES-VIS RVs.</caption><thead><tr><th>Fitted parameter</th><th colspan="2">Best-fit value <sup>(a)</sup></th></tr><tr><th></th><th>e<sub>b</sub> = 0</th><th>e<sub>b</sub> ≠ 0</th></tr></thead><tbody><tr><td>h<sub>1,HARPS</sub> (m s<sup>-1</sup>)</td><td>2.1<sup>+0.4</sup><sub>-0.3</sub></td><td><b>2.2</b><sup>+0.4</sup><sub>-0.3</sub></td></tr><tr><td>h<sub>2,HARPS</sub> (m s<sup>-1</sup>)</td><td>1.4 ± 0.4</td><td><b>1.4</b>±0.4</td></tr><tr><td>h<sub>1,CARMENES</sub> (m s<sup>-1</sup>)</td><td>1.1 ± 0.3</td><td><b>1.0</b><sup>+0.3</sup><sub>-0.2</sub></td></tr><tr><td>h<sub>2,CARMENES</sub> (m s<sup>-1</sup>)</td><td>1.3 ± 0.3</td><td><b>1.3</b>±0.3</td></tr><tr><td>θ (d)</td><td>30.7 ± 0.3</td><td><b>30.6</b>±0.3</td></tr><tr><td>λ<sub>QPC</sub> (d)</td><td>121<sup>+30</sup><sub>-28</sub></td><td><b>122</b><sup>+27</sup><sub>-24</sub></td></tr><tr><td>w</td><td>0.34<sup>+0.05</sup><sub>-0.04</sub></td><td><b>0.34</b><sup>+0.05</sup><sub>-0.04</sub></td></tr><tr><td>γ<sub>HARPS pre-2015</sub> (m s<sup>-1</sup>)</td><td>-0.3 ± 0.5</td><td><b>-0.4</b>±0.5</td></tr><tr><td>γ<sub>HARPS post-2015</sub> (m s<sup>-1</sup>)</td><td>0.0 ± 1.1</td><td><b>0.0</b>±1.1</td></tr><tr><td>γ<sub>CARMENES</sub> (m s<sup>-1</sup>)</td><td>0.1 ± 0.3</td><td><b>0.1</b>±0.3</td></tr><tr><td>σ<sub>jit, HARPS pre-2015</sub> (m s<sup>-1</sup>)</td><td>0.8 ± 0.1</td><td><b>0.7</b>±0.1</td></tr><tr><td>σ<sub>jit, HARPS post-2015</sub> (m s<sup>-1</sup>)</td><td>1.5<sup>+0.5</sup><sub>-0.4</sub></td><td><b>1.6</b><sup>+0.5</sup><sub>-0.4</sub></td></tr><tr><td>σ<sub>jit, CARMENES</sub> (m s<sup>-1</sup>)</td><td>0.9 ± 0.2</td><td><b>0.9</b>±0.2</td></tr><tr><td>K<sub>b</sub> (m s<sup>-1</sup>)</td><td>0.91 ± 0.18</td><td><b>1.15</b><sup>+0.21</sup><sub>-0.19</sub></td></tr><tr><td>P<sub>b</sub> (d)</td><td>139.90<sup>+0.68</sup><sub>-0.65</sub></td><td><b>140.43</b>±0.41</td></tr><tr><td>T<sub>conj,b</sub> (BJD-2 450 000)</td><td>8670.29<sup>+6.43</sup><sub>-6.69</sub></td><td><b>8696.21</b><sup>+7.66</sup><sub>-13.71</sub></td></tr><tr><td>√e<sub>b</sub> cos ω<sub>★, b</sub></td><td>—</td><td><b>-0.604</b><sup>+0.267</sup><sub>-0.143</sub></td></tr><tr><td>√e<sub>b</sub> sin ω<sub>★, b</sub></td><td>—</td><td><b>-0.209</b><sup>+0.258</sup><sub>-0.261</sub></td></tr><tr><td>K<sub>365-d, CARMENES</sub> (m s<sup>-1</sup>)</td><td>1.2 ± 0.3</td><td><b>1.2</b>±0.3</td></tr><tr><td>P<sub>365-d, CARMENES</sub> (d)</td><td>356.98<sup>+7.92</sup><sub>-8.15</sub></td><td><b>357.76</b><sup>+7.29</sup><sub>-8.15</sub></td></tr><tr><td>T<sub>0,365-day CARMENES</sub> (BJD-2 450 000)</td><td>8784.4<sup>+18.7</sup><sub>-19.4</sub></td><td><b>8785.93</b><sup>+17.56</sup><sub>-19.13</sub></td></tr><tr><td colspan="3"><strong>Derived parameter</strong></td></tr><tr><td>eccentricity, e<sub>b</sub></td><td>—</td><td><b>0.45</b><sup>+0.15</sup><sub>-0.14</sub></td></tr><tr><td>arg. of periapsis, ω<sub>★, b</sub></td><td>—</td><td><b>-2.49</b><sup>+5.35</sup><sub>-0.47</sub></td></tr><tr><td>min. mass, m<sub>b</sub> sin i<sub>b</sub> (M<sub>⊕</sub>)</td><td>4.7±1.0</td><td><b>5.2</b>±0.9</td></tr><tr><td>semi-major axis, a<sub>b</sub> (au)</td><td>0.421<sup>+0.014</sup><sub>-0.015</sub></td><td><b>0.422</b><sup>+0.014</sup><sub>-0.015</sub></td></tr><tr><td>periapsis (au)</td><td>—</td><td><b>0.231</b><sup>+0.058</sup><sub>-0.060</sub></td></tr><tr><td>apoapsis (au)</td><td>—</td><td><b>0.612</b><sup>+0.064</sup><sub>-0.061</sub></td></tr><tr><td>equilibrium temperature, T<sub>eq,b</sub> (K)</td><td>196 ± 10<sup>(b)</sup></td><td>Orbit-averaged<sup>(c)</sup>: <b>202</b>±<b>11</b><br/>Apoapsis: <b>162</b><sup>+12</sup><sub>-11</sub><br/>Periapsis: <b>264</b><sup>+45</sup><sub>-31</sub></td></tr><tr><td>Insolation flux<sup>(d)</sup>, S<sub>b</sub> (S<sub>⊕</sub>)</td><td>0.24<sup>+0.06</sup><sub>-0.05</sub></td><td>Orbit-averaged: <b>0.28</b><sup>+0.07</sup><sub>-0.06</sub><br/>Apoapsis: <b>0.114</b><sup>+0.037</sup><sub>-0.030</sub><br/>Periapsis: <b>0.79</b><sup>+0.72</sup><sub>-0.31</sub></td></tr><tr><td>ln Z</td><td>-958.0</td><td><b>-955.2</b></td></tr><tr><td>ln Z<sub>1p</sub>-ln Z<sub>0p</sub></td><td>+3.3</td><td><b>+6.1</b></td></tr></tbody></table>

**Notes.** The model with $e_b \neq 0$ (values in bold) is our adopted solution. <sup>(a)</sup>The uncertainties are calculated as the 16th and 84th percentiles of the posterior distributions. <sup>(b)</sup>Derived from the relation $T_{eq} = T_{eff} \cdot \sqrt{\frac{R_*}{2a_b}} \cdot (1 - A_B)^{0.25}$, assuming Bond albedo $A_B = 0$. <sup>(c)</sup>This is the average equilibrium temperature based on the stellar flux received by the planet averaged over the eccentric orbit. This flux-averaged temperature scales with the eccentricity as $(1 - e_b^2)^{-1/8}$ with respect to the value for a circular orbit. <sup>(d)</sup>For the circular orbit, it is derived from the equation $S_b = \frac{L_*}{L_\odot} \cdot (\frac{\mu w}{a_b})^2$. For the eccentric orbit, the temporal average insolation flux scales with the eccentricity as $\frac{1}{\sqrt{1 - e_b^2}}$ with respect to the value for a circular orbit with the same semi-major axis (e.g. Williams & Pollard 2002).
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Posterior distributions of the semi-amplitude $K_b$ (top panel), orbital period $P_b$ (middle panel), and eccentricity $e_b$ (bottom panel) of the candidate planet Gl514b derived from a GP QPC analysis of the HARPS<sub>TERRA</sub>+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 HARPS<sub>TERRA</sub>+CARMENES dataset. Here, $\Delta \ln Z$ denotes the Bayesian evidence difference between a model with and without a Keplerian included.</figcaption></figure>

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.

<table><caption>**Table 6.** Logarithmic Bayesian evidence $\ln Z_{2\text{pl}}$ for the GP models with two Keplerians that we tested on HARPS<sub>TERRA</sub>+CARMENES VIS RVs.</caption><thead><tr><th>GP kernel</th><th>ln Z<sub>2pl</sub></th><th>Δ ln Z<sub>2pl-1pl</sub><sup>(a)</sup></th></tr></thead><tbody><tr><td>QP</td><td>-960.4</td><td>0</td></tr><tr><td>QPC</td><td>-953.9</td><td>+1.3</td></tr><tr><td>dSHO</td><td>-956.2</td><td>+0.5</td></tr></tbody></table>

**Notes.** <sup>(a)</sup> $\ln Z_{1\text{pl}}$ is the Bayesian evidence for the model with the same GP kernel and a single Keplerian.

**Test 1: QP kernel.** For planet Gl 514 b we find $K_b = 1.28^{+0.21}_{-0.18}$ m s<sup>-1</sup>, $P_b = 140.24^{+0.44}_{-0.56}$ days, and $e_b = 0.46^{+0.11}_{-0.12}$, and those of the possible innermost companion are $K_c = 0.74 \pm 0.31$ m s<sup>-1</sup>, $P_c = 63.67^{+0.29}_{-0.11}$ days, and $e_c = 0.53^{+0.18}_{-0.31}$. 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 \pm 0.17$ m s<sup>-1</sup>, $P_b = 140.40^{+0.27}_{-0.39}$ days, and $e_b = 0.47^{+0.12}_{-0.08}$, while those of the possible innermost companion are $K_c = 0.92^{+0.19}_{-0.30}$ m s<sup>-1</sup>, $P_c = 63.64^{+0.20}_{-1.93}$ days, and $e_c = 0.58^{+0.13}_{-0.20}$. The Bayesian evidence is $\ln Z = -953.9$, which is slightly higher than the Bayesian evidence for the model that includes only one Keplerian ($\Delta \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 $K_b = 1.22^{+0.22}_{-0.20}$ m s<sup>-1</sup>, $P_b = 140.55^{+0.44}_{-0.54}$ days, and $e_b = 0.50^{+0.12}_{-0.13}$, while for the possible innermost companion we find $K_c = 0.88 \pm 0.30$ m s<sup>-1</sup>, $P_c = 63.76^{+0.15}_{-0.25}$ days, and $e_c = 0.51^{+0.18}_{-0.26}$. 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 ($\Delta \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<sup>-1</sup> 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 sam-
pling from the posterior (Foreman-Mackey et al. 2013), we got
100 000 posterior samples (each with a different orbital config-
uration) for each of the four scenarios. To these, we applied
a ‘filter’ in order to select only the configurations that guar-
antee the dynamical stability of the system over 10⁵ orbits of
Gl 514 b. We adopted the following stability criteria: (i) avoid-
ing 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, differ-
ent 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 $K = 0.4_{-0.2}^{+0.3}$ m s⁻¹, $P = 1579_{-61}^{+1872}$ 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 ẋ = -0.00026 ± 0.00017 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 514 b.

# 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

<figure><img src="image_14.png" /></figure>

**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.

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σ<sub>sky</sub> 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 \sim 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 mean-
ingless, 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 simula-
tions, moving the time of central transit $T_{\text{conj}}$ recursively two
days ahead of the $T_{\text{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 recov-
ered the transits for which $T_{\text{conj, injected}} - T_{\text{conj, recovered}} < 0.1$ d.
We recovered all the injected transits with $R \ge 1.0 R_⊕$. All of
<figure><img src="image_15.png" /><figcaption>**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 - 2457000). Panel c shows the BLS periodogram. The red triangle indicates the period associated with the highest S/N.</figcaption></figure>

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 \pm 0.9 M_⊕$ would correspond to a radius $R_b = 2.1^{+0.9}_{-0.6} R_⊕$, following the probabilistic mass-radius relation of (Chen & Kipping 2017). We show that we are able to detect a single transit

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Posterior distributions for the conjunction times $T_{conj,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.</figcaption></figure>

due to a planet with a radius as large as that predicted for Gl 514 b within $1\sigma$, 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.

## 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$^{10}$. 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

<footnote>$^{10}$ For reference, a threshold of RUWE $\ge 1.4$ (Lindegren et al. 2018) is typically used to identify astrometrically variable stars.</footnote>
<figure><img src="image_17.png" /><figcaption>**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).</figcaption></figure>

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*<sub>J</sub> can be ruled out at the 1-σ level. A companion with the same value of *m*<sub>p</sub> sin *i*<sub>p</sub> 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<sup>-1</sup>, 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*<sub>jup</sub>, with a constraint on the true inclination angle ≥1°).

**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 \pm 0.9 M_\oplus$ (with an 87% probability that the true mass is within a factor of two larger) moving on an eccentric orbit ($e_b = 0.45^{+0.15}_{-0.14}$) with orbital period $P_b = 140.43 \pm 0.41$ days and semi-major axis $a_b = 0.422^{+0.014}_{-0.015}$ 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 514 b 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*<sub>jup</sub> planet at a distance of ~3–10 au, 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 514 b 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 HARPS<sub>TERRA</sub> 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<sup>11</sup>, 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

<footnote n="11">¹¹ We note that Wittenmyer et al. (2019) concluded that planet candidates with eccentricity $e \ge 0.5$ are unlikely to be impostors.</footnote>
<figure><img src="image_18.png" /></figure>

**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.

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.

## 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<fnref n="12" />, 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 \pm 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<fnref n="13" /> (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⁻⁹, a threshold that the Gl 514 system could realistically exceed (see Fig. 1 in Kasper et al. 2021).

<footnote n="12">¹² We used a list of planets with mass or minimum mass <10 $M_⊕$ that we compiled querying the NASA Exoplanet Archive, and taking the missing data from the catalogue maintained by the Planetary Habitability Laboratory or from the most recent references. For the mass and radius of some planets, we used updated values: Trappist-1 d, Trappist-1 e, Trappist-1 f, and Trappist-1 g; Agol et al. (2021); Proxima b: Suárez Mascareño et al. (2020); K2-18 b: Benneke et al. (2019). For the effective temperature and luminosity of Gl 229 A c we used values from Schweitzer et al. (2019). For Gl 163 c we used the effective temperature determined by Tuomi & Anglada-Escudé (2013). Gl 832 c is likely an artefact of the stellar activity (Suárez Mascareño et al. 2017; Gorrini et al. 2022).</footnote>

<footnote n="13">¹³ For Gl 514 b we plotted the orbit-averaged insolation.</footnote>
<figure><img src="image_19.png" /><figcaption>**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.</figcaption></figure>

## 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 ($\Delta v \sim 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, $\epsilon$, 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
<figure><img src="image_20.png" /><figcaption>**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. $M_A = v_{rel}/v_{Alfv} \le 1$; top panel) at distances equal to the semi-major axis and periapsis (vertical dashed line).</figcaption></figure>

advisable in order to better constrain the radio emission in this system.

**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 Centro 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 Ministerio 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, Landesfernwarte Königstuhl, Institut de Ciències de l'Espai, Institut für Astrophysik Göttingen, Universität Complutense de Madrid, Thüringer Landesfernwarte 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 Ministerio 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-100, PID2019-109522GB-C5[1:4], and PGC2018-098153-B-C33, and the Centre of Excellence “Severo Ochoa” and “María 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 Astrobiología (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-NAF 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 María de Maeztu CEX2020-001058-M, and by the Generalitat de Catalunya/CERCA programme. M. Perger also acknowledges support from Spanish grant PID2020-120375GB-100 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 Mascareño 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-100 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.

## References

* Affer, L., Micela, G., Damasso, M., et al. 2016, A&A, 593, A117

* Affer, L., Damasso, M., Micela, G., et al. 2019, A&A, 622, A193

* Agol, E., Dorn, C., Grimm, S. L., et al. 2021, Planet. Sci. J., 2, 1

* Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O'Neil, M. 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38, 252

* Anglada-Escudé, G., & Butler, R. P. 2012, ApJS, 200, 15

* Anglada-Escudé, G., López-Morales, M., & Chambers, J. E. 2010, ApJ, 709, 168

* Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437

* Artigau, É., Kouach, D., Donati, J.-F., et al. 2014, SPIE Conf. Ser., 9147, 914715

* Astudillo-Defru, N., Bonfils, X., Delfosse, X., et al. 2015, A&A, 575, A119

* Astudillo-Defru, N., Delfosse, X., Bonfils, X., et al. 2017a, A&A, 600, A13

* Astudillo-Defru, N., Forveille, T., Bonfils, X., et al. 2017b, A&A, 602, A88

* Bailer-Jones, C. A. L., Rybizki, J., Andrae, R., & Fouesneau, M. 2018a, A&A, 616, A37

* Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018b, AJ, 156, 58

* Barnes, R., Luger, R., Deitrick, R., et al. 2020, PASP, 132, 024502

* Benatti, S., Damasso, M., Borsa, F., et al. 2021, A&A, 650, A66

* Benneke, B., Wong, I., Piaulet, C., et al. 2019, ApJ, 887, L14

* Boisvert, J. H., Nelson, B. E., & Steffen, J. H. 2018, MNRAS, 480, 2846

* Bolmont, E., Libert, A.-S., Leconte, J., & Selsis, F. 2016, A&A, 591, A106

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

* Bonfils, X., Astudillo-Defru, N., Díaz, R., et al. 2018, A&A, 613, A25

* Bouchy, F., Pepe, F., & Queloz, D. 2001, A&A, 374, 733

* Brandt, T. D. 2021, ApJS, 254, 42

* Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

* Caballero, J. A., Guàrdia, J., López del Fresno, M., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Observatory Operations: Strategies, Processes, and Systems VI, eds. A. B. Peck, R. L. Seaman, & C. R. Benn, 9910, 99100E

* Callingham, J. R., Vedantham, H. K., Shimwell, T. W., et al. 2021, Nat. Astron., 5, 1233

* Cegla, H. 2019, Geosciences, 9, 114

* Chen, J., & Kipping, D. 2017, ApJ, 834, 17

* Chen, S. S., Donoho, D. L., & Saunders, M. A. 1998, SIAM J. Scientific Comput., 20, 33

* Cincotta, P. M., & Simó, C. 2000, A&As, 147, 205

* Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica D, 182, 151

* Cloutier, R., Astudillo-Defru, N., Doyon, R., et al. 2017, A&A, 608, A35

* Crossfield, I. J. M., Petigura, E., Schlieder, J. E., et al. 2015, ApJ, 804, 10

* Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog; II/246

* Damasso, M., Bonomo, A. S., Astudillo-Defru, N., et al. 2018, A&A, 615, A69

* Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45

* Endl, M., Cochran, W. D., Tull, R. G., & MacQueen, P. J. 2003, AJ, 126, 3099

* Faria, J. P., Suárez Mascareño, A., Figueira, P., et al. 2022, A&A, 658, A115

* Feng, F., Butler, R. P., Shectman, S. A., et al. 2020, ApJS, 246, 11

* Feroz, F., Balan, S. T., & Hobson, M. P. 2011, MNRAS, 415, 3462

* Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, Open J. Astrophys., 2, 10
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220

Fuhrmeister, B., Czesla, S., Schmitt, J. H. M. M., et al. 2019, A&A, 623, A24

Gagné, J., & Faherty, J. K. 2018, ApJ, 862, 138

Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456

Gomes da Silva, J., Figueira, P., Santos, N., & Faria, J. 2018, J. Open Source Softw., 3, 667

Gorrini, P., Astudillo-Defru, N., Dreizler, S., et al. 2022, A&A, 664, A64

Hara, N. C., Boué, G., Laskar, J., & Correia, A. C. M. 2017, MNRAS, 464, 1220

Hara, N. C., Boué, G., Laskar, J., Delisle, J. B., & Unger, N. 2019, MNRAS, 489, 738

Hara, N. C., Bouchy, F., Stalport, M., et al. 2020, A&A, 636, A6

Hara, N. C., Delisle, J.-B., Unger, N., & Dumusque, X. 2022, A&A, 658, A177

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

Hsu, D. C., Ford, E. B., & Terrien, R. 2020, MNRAS, 498, 2249

Johnson, D. R. H., & Soderblom, D. R. 1987, AJ, 93, 864

Kaltenegger, L., Madden, J., Lin, Z., et al. 2019, ApJ, 883, L40

Kane, S. R., & Gelino, D. M. 2012, Astrobiology, 12, 940

Kasper, M., Cerpa Urra, N., Pathak, P., et al. 2021, The Messenger, 182, 38

Kervella, P., Arenou, F., Mignard, F., & Thévenin, F. 2019, A&A, 623, A72

Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7

Koen, C., Kilkenny, D., van Wyk, F., & Marang, F. 2010, MNRAS, 403, 1949

Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131

Kopparapu, R. K., Ramirez, R. M., SchottelKotte, J., et al. 2014, ApJ, 787, L29

Kosiarek, M. R., Crossfield, I. J. M., Hardegree-Ullman, K. K., et al. 2019, AJ, 157, 97

Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369

Kürster, M., Trifonov, T., Reffert, S., Kostogryz, N. M., & Rodler, F. 2015, A&A, 577, A103

Lanza, A. F. 2009, A&A, 505, 339

Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2

Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2

Mahadevan, S., Ramsey, L., Bender, C., et al. 2012, SPIE Conf. Ser., 8446, 84461S

Maldonado, J., Affer, L., Micela, G., et al. 2015, A&A, 577, A132

Melrose, D. B., & Dulk, G. A. 1982, ApJ, 259, 844

Méndez, A., & Rivera-Valentín, E. G. 2017, ApJ, 837, L1

Morton, T. D., Bryson, S. T., Coughlin, J. L., et al. 2016, ApJ, 822, 86

Nardiello, D. 2020, MNRAS, 498, 5972

Nardiello, D., Bedin, L. R., Nascimbeni, V., et al. 2015, MNRAS, 447, 3536

Nardiello, D., Libralato, M., Bedin, L. R., et al. 2016, MNRAS, 455, 2337

Nardiello, D., Borsato, L., Piotto, G., et al. 2019, MNRAS, 490, 3806

Nardiello, D., Piotto, G., Deleuil, M., et al. 2020, MNRAS, 495, 4924

Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96

Pérez-Torres, M., Gómez, J. F., Ortiz, J. L., et al. 2021, A&A, 645, A77

Perger, M., García-Piquer, A., Ribas, I., et al. 2017, A&A, 598, A26

Perger, M., Anglada-Escudé, G., Ribas, I., et al. 2021, A&A, 645, A58

Pinamonti, M., Damasso, M., Marzari, F., et al. 2018, A&A, 617, A104

Pinamonti, M., Sozzetti, A., Maldonado, J., et al. 2022, A&A, 664, A65

Quintana, E. V., Barclay, T., Raymond, S. N., et al. 2014, Science, 344, 277

Quirrenbach, A. 2022, Res. Notes AAS, 6, 56

Quirrenbach, A., Amado, P. J., Ribas, I., et al. 2018, in Ground-based and Air-borne Instrumentation for Astronomy VII, eds. C. J. Evans, L. Simard, & H. Takami, 10702, International Society for Optics and Photonics (SPIE), 246

Rein H., & Liu, S. F. 2012, A&A, 537, A128

Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376

Ricker, G. R., Vanderspek, R., Winn, J., et al. 2016, SPIE Conf. Ser., 9904, 99042B

Rodigas, T. J., & Hinz, P. M. 2009, ApJ, 702, 716

Sabotta, S., Schlecker, M., Chaturvedi, P., et al. 2021, A&A, 653, A114

Saur, J., Grambusch, T., Duling, S., Neubauer, F. M., & Simon, S. 2013, A&A, 552, A119

Schweitzer, A., Passegger, V. M., Cifuentes, C., et al. 2019, A&A, 625, A68

Shulyak, D., Reiners, A., Nagel, E., et al. 2019, A&A, 626, A86

Suárez-Mascareño, A., Rebolo, R., González Hernández, J. I., & Esposito, M. 2015, MNRAS, 452, 2745

Suárez Mascareño, A., Rebolo, R., González Hernández, J. I., & Esposito, M. 2017, MNRAS, 468, 4772

Suárez Mascareño, A., Faria, J. P., Figueira, P., et al. 2020, A&A, 639, A77

Tal-Or, L., Trifonov, T., Zucker, S., Mazeh, T., & Zechmeister, M. 2019, MNRAS, 484, L8

Torres, G., Kipping, D. M., Fressin, F., et al. 2015, ApJ, 800, 99

Trifonov, T., Kürster, M., Zechmeister, M., et al. 2017, A&A, 602, L8

Trifonov, T., Tal-Or, L., Zechmeister, M., et al. 2020, A&A, 636, A74

Tuomi, M., & Anglada-Escudé, G. 2013, A&A, 556, A111

Turnpenney, S., Nichols, J. D., Wynn, G. A., & Burleigh, M. R. 2018, ApJ, 854, 72

Vedantham, H. K., Callingham, J. R., Shimwell, T. W., et al. 2020, Nat. Astron., 4, 577

Williams, D. M., & Pollard, D. 2002, Int. J. Astrobiol., 1, 61

Wittenmyer, R. A., Wang, S., Horner, J., et al. 2013, ApJS, 208, 2

Wittenmyer, R. A., Bergmann, C., Horner, J., Clark, J., & Kane, S. R. 2019, MNRAS, 484, 4230

Zarka, P. 2007, Planet. Space Sci., 55, 598

Zechmeister, M., Kürster, M., & Endl, M. 2009, A&A, 505, 859

Zechmeister, M., Anglada-Escudé, G., & Reiners, A. 2014, A&A, 561, A59

Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12

Zechmeister, M., Dreizler, S., Ribas, I., et al. 2019, A&A, 627, A49

<footnote n="1">¹ INAF – Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese, Italy
e-mail: mario.damasso@inaf.it</footnote>

<footnote n="2">² Institut de Ciències de l'Espai (ICE, CSIC), Campus UAB, Carrer de Can Magrans s/n, 08193 Bellaterra, Spain</footnote>

<footnote n="3">³ Institut d'Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain</footnote>

<footnote n="4">⁴ Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="5">5 Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="6">⁶ INAF – Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122, Padova, Italy</footnote>

<footnote n="7">⁷ Instituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain</footnote>

<footnote n="8">⁸ Observatoire Astronomique de l'Université de Genève, 51 Chemin des Pegasi, 1290 Versoix, Switzerland</footnote>

<footnote n="9">⁹ Landessternwarte, Zentrum für Astronomie der Universität Heidelberg, Königstuhl 12, 69117 Heidelberg, Germany</footnote>

<footnote n="10">¹⁰ Centro de Astrobiología (CSIC-INTA), Carretera de Ajalvir km 4, 28850 Torrejón de Ardoz, Madrid, Spain</footnote>

<footnote n="11">¹¹ Departamento de Matemática y Física Aplicadas, Universidad Católica de la Santísima Concepción, Concepción, Chile</footnote>

<footnote n="12">¹² Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="13">¹³ Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="14">¹⁴ Centro de Astrobiología (CAB, CSIC-INTA), Depto. de Astrofísica, ESAC campus, 28692 Villanueva de la Cañada (Madrid), Spain</footnote>

<footnote n="15">¹⁵ Physikalisches Institut, University of Bern, Gesellsschaftstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="16">¹⁶ European Southern Observatory, Alonso de Cordova 3107, Vitacura, Santiago, Chile</footnote>

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

<footnote n="18">¹⁸ Observatorio de Calar Alto, Sierra de los Filabres, 04550-Gérgal, Almería, Spain</footnote>

<footnote n="19">¹⁹ Thüringer Landessternwarte Tautenburg, Sternwarte 5, 07778 Tautenburg, Germany</footnote>

<footnote n="20">²⁰ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="21">²¹ Departamento de Física de la Tierra y Astrofísica & IPARCOS-UCM (Instituto de Física de Partículas y del Cosmos de la UCM), Facultad de Ciencias Físicas, Universidad Complutense de Madrid, 28040 Madrid, Spain</footnote>

<footnote n="22">²² Institut für Astrophysik, Georg-August-Universität, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany</footnote>

<footnote n="23">²³ Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany</footnote>
# 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.mEDIA-hd.mpg.de/homes/trifonov/HARPS_RVBank.html. Radial velocities derived from HIRES spectra are available at https://cdsarc.unistra.fr/viz-bin/cat/J/MNRAS/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:

$$k_{QP}(t, t') = h^2 \exp \left[ - \frac{(t-t')^2}{2\lambda_{QP}^2} - \frac{\sin^2(\pi(t-t')/\theta)}{2w^2} \right] + (\sigma_{RV}^2(t) + \sigma_{jit}^2) \delta_{t,t'}. \quad (B.1)$$

Here, $t$ and $t'$ represent two different epochs of observations, $\sigma_{RV}$ is the RV uncertainty, and $\delta_{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 $\sigma_{jit}$ which is added in quadrature to the formal uncertainties $\sigma_{RV}$. The GP hyper-parameters are $h$, which denotes the scale amplitude of the correlated signal; $\theta$, 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 $\lambda_{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

$$k_{QPC}(t, t') = \exp\left(-2\frac{(t-t')^2}{\lambda_{QPC}^2}\right) \cdot \left[ h_1^2 \exp\left(-\frac{1}{2w^2}\sin^2\left(\frac{\pi(t-t')}{\theta}\right)\right) + h_2^2 \cos\left(\frac{4\pi(t-t')}{\theta}\right) \right] + (\sigma_{RV}^2(t) + \sigma_{jit}^2) \delta_{t,t'}. \quad (B.2)$$

Again, $t$ and $t'$ represent two different epochs of observations; $h_1$ and $h_2$ are scale amplitudes; $\theta$ 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; $\lambda_{QPC}$ is defined as $2 \cdot \lambda_{QP}$, better representing the average lifetime of the activity-related features responsible for the stellar correlated signal in the RVs (Perger et al. 2021); $\sigma_{RV}$ and $\sigma_{jit}$ are the RV uncertainty and instrument-dependent jitter, respectively, and $\delta_{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)<fnref n="14" />. The power spectral density corresponding to this kernel is

$$S(\omega) = \sqrt{\frac{2}{\pi} \frac{S_1 \omega_1^4}{(\omega^2 - \omega_1^2)^2 + 2\omega_1^2 \omega^2}} + \sqrt{\frac{2}{\pi} \frac{S_2 \omega_2^4}{(\omega^2 - \omega_2^2)^2 + 2\omega_2^2 \omega^2 / Q^2}}, \quad (B.3)$$

<footnote n="14">¹⁴ https://github.com/dfm/celerite/blob/main/celerite/terms.py</footnote>
where

$$S_1 = \frac{A^2}{\omega_1 Q_1 (1 + f)}, \qquad (B.4)$$

$$S_2 = \frac{A^2}{\omega_2 Q_2 (1 + f)} \cdot f, \qquad (B.5)$$

$$\omega_1 = \frac{4\pi Q_1}{P_{\text{rot}} \sqrt{4Q_1^2 - 1}}, \qquad (B.6)$$

$$\omega_2 = \frac{8\pi Q_2}{P_{\text{rot}} \sqrt{4Q_2^2 - 1}}, \qquad (B.7)$$

$$Q_1 = \frac{1}{2} + Q_0 + \Delta Q, \qquad (B.8)$$

$$Q_2 = \frac{1}{2} + Q_0. \qquad (B.9)$$

The parameters in [B.4-B.9], where the subscripts 1 and 2 refer to the primary ($P_{\star,rot}$) and secondary ($P_{\star,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 $\omega_i$ are derived through Eq. [B.4-B.9]: the variability amplitude A, the stellar rotation period $P_{\star,rot}$, the quality factor $Q_0$, the difference $\Delta 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 sinusoid 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.
<figure><img src="image_21.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. C.2.** Phase-folded RVs of HARPS<sub>TERRA</sub> 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.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. C.3.** Sections of the residuals of the whole set of HARPS<sub>TERRA</sub>+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.**</figcaption></figure>
<figure><img src="image_24.png" /></figure>

**Fig. C.4.** GLS periodograms of the HARPS<sub>TERRA</sub>+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.
<table><caption>**Table C.1.** Best-fit parameters obtained for the model with a QP kernel applied to HARPS<sub>TERRA</sub> and CARMENES-VIS RVs.</caption><thead><tr><th rowspan="2">Fitted parameter</th><th colspan="2">Best-fit value<sup>a</sup></th></tr><tr><th>e<sub>b</sub>=0</th><th>e<sub>b</sub>≠0</th></tr></thead><tbody><tr><td>h<sub>HARPS</sub> [m s<sup>-1</sup>]</td><td>2.5 ± 0.3</td><td>2.6 ± 0.3</td></tr><tr><td>h<sub>CARMENES</sub> [m s<sup>-1</sup>]</td><td>1.7 ± 0.2</td><td>1.7 ± 0.2</td></tr><tr><td>θ [d]</td><td>30.7 ± 0.3</td><td>30.7 ± 0.3</td></tr><tr><td>λ<sub>QP</sub> [d]</td><td>63 ± 15</td><td>63<sup>+15</sup><sub>-12</sub></td></tr><tr><td>w</td><td>0.31 ± 0.04</td><td>0.31 ± 0.03</td></tr><tr><td>γ<sub>HARPS pre-2015</sub> [m s<sup>-1</sup>]</td><td>-0.5 ± 0.5</td><td>-0.5 ± 0.5</td></tr><tr><td>γ<sub>HARPS post-2015</sub> [m s<sup>-1</sup>]</td><td>0.1 ± 1.2</td><td>0.1 ± 1.2</td></tr><tr><td>γ<sub>CARMENES</sub> [m s<sup>-1</sup>]</td><td>0.1 ± 0.4</td><td>0.2 ± 0.4</td></tr><tr><td>σ<sub>jit, HARPS pre-2015</sub> [m s<sup>-1</sup>]</td><td>0.8 ± 0.1</td><td>0.7 ± 0.1</td></tr><tr><td>σ<sub>jit, HARPS post-2015</sub> [m s<sup>-1</sup>]</td><td>1.6<sup>+0.5</sup><sub>-0.4</sub></td><td>1.7<sup>+0.5</sup><sub>-0.4</sub></td></tr><tr><td>σ<sub>jit, CARMENES</sub> [m s<sup>-1</sup>]</td><td>1.0 ± 0.2</td><td>0.9 ± 0.2</td></tr><tr><td>K<sub>b</sub> [m s<sup>-1</sup>]</td><td>1.01 ± 0.18</td><td>1.23 ± 0.20</td></tr><tr><td>P<sub>b</sub> [d]</td><td>139.75<sup>+0.65</sup><sub>-0.61</sub></td><td>140.37<sup>+0.44</sup><sub>-0.46</sub></td></tr><tr><td>T<sub>conj, b</sub> [BJD-2450000]</td><td>8668.68<sup>+6.67</sup><sub>-6.29</sub></td><td>8692.00<sup>+10.38</sup><sub>-15.44</sub></td></tr><tr><td>√e<sub>b</sub> cos ω<sub>★, b</sub></td><td>-</td><td>-0.527<sup>+0.306</sup><sub>-0.180</sub></td></tr><tr><td>√e<sub>b</sub> sin ω<sub>★, b</sub></td><td>-</td><td>-0.306<sup>+0.268</sup><sub>-0.221</sub></td></tr><tr><td>K<sub>365-d, CARMENES</sub> [m s<sup>-1</sup>]</td><td>1.2 ± 0.4</td><td>1.3 ± 0.4</td></tr><tr><td>P<sub>365-d, CARMENES</sub> [d]</td><td>359.77<sup>+6.94</sup><sub>-8.59</sub></td><td>359.81<sup>+6.65</sup><sub>-8.69</sub></td></tr><tr><td>T<sub>0,365-d CARMENES</sub> [BJD-2450000]</td><td>8784.5 ± 22.6</td><td>8786.22<sup>+20.34</sup><sub>-21.93</sub></td></tr><tr><td colspan="3"><strong>Derived parameter</strong></td></tr><tr><td>eccentricity, e<sub>b</sub></td><td>-</td><td>0.42<sup>+0.14</sup><sub>-0.13</sub></td></tr><tr><td>arg. of periapsis, ω<sub>★, b</sub></td><td>-</td><td>-2.30<sup>+5.27</sup><sub>-0.66</sub></td></tr><tr><td>ln Z</td><td>-962.9</td><td>-960.4</td></tr><tr><td>ln Z<sub>1p</sub>-ln Z<sub>0p</sub></td><td>+2.7</td><td>+5.2</td></tr></tbody></table>

**Notes.** (a) The uncertainties of the parameters are calculated as the 16th and 84th percentiles of their posterior distributions.
<table><caption>**Table C.2.** Best-fit parameters obtained for the model with a rotational (dSHO) kernel applied to HARPS<sub>TERRA</sub> and CARMENES-VIS RVs.</caption><thead><tr><th rowspan="2">Fitted parameter</th><th colspan="2">Best-fit value<sup>a</sup></th></tr><tr><th>e<sub>b</sub>=0</th><th>e<sub>b</sub>≠0</th></tr></thead><tbody><tr><td>log A<sub>HARPS</sub></td><td>0.9 ± 0.1</td><td>0.9 ± 0.1</td></tr><tr><td>log A<sub>CARMENES</sub></td><td>0.6 ± 0.1</td><td>0.6 ± 0.1</td></tr><tr><td>θ [d]</td><td>31.0 ± 0.5</td><td>30.9<sup>+0.4</sup><sub>-0.5</sub></td></tr><tr><td>log Q<sub>0</sub></td><td>2.11<sup>+0.39</sup><sub>-0.42</sub></td><td>2.25<sup>+0.39</sup><sub>-0.42</sub></td></tr><tr><td>log ΔQ</td><td>-3.80<sup>+4.41</sup><sub>-4.16</sub></td><td>-3.77<sup>+4.47</sup><sub>-4.23</sub></td></tr><tr><td>f</td><td>3.70<sup>+2.57</sup><sub>-1.42</sub></td><td>3.79<sup>+2.54</sup><sub>-1.46</sub></td></tr><tr><td>γ<sub>HARPS pre-2015</sub> [m s<sup>-1</sup>]</td><td>-0.2 ± 0.2</td><td>-0.1 ± 0.2</td></tr><tr><td>γ<sub>HARPS post-2015</sub> [m s<sup>-1</sup>]</td><td>-0.2 ± 0.6</td><td>-0.3 ± 0.6</td></tr><tr><td>γ<sub>CARMENES</sub> [m s<sup>-1</sup>]</td><td>0.2 ± 0.2</td><td>0.2 ± 0.2</td></tr><tr><td>σ<sub>jit, HARPS pre-2015</sub> [m s<sup>-1</sup>]</td><td>0.8 ± 0.1</td><td>0.8 ± 0.1</td></tr><tr><td>σ<sub>jit, HARPS post-2015</sub> [m s<sup>-1</sup>]</td><td>1.0<sup>+0.4</sup><sub>-0.3</sub></td><td>1.1<sup>+0.5</sup><sub>-0.4</sub></td></tr><tr><td>σ<sub>jit, CARMENES</sub> [m s<sup>-1</sup>]</td><td>0.6 ± 0.3</td><td>0.6<sup>+0.2</sup><sub>-0.3</sub></td></tr><tr><td>K<sub>b</sub> [m s<sup>-1</sup>]</td><td>0.88 ± 0.16</td><td>1.16<sup>+0.25</sup><sub>-0.21</sub></td></tr><tr><td>P<sub>b</sub> [d]</td><td>140.23<sup>+0.67</sup><sub>-0.64</sub></td><td>140.72<sup>+0.44</sup><sub>-0.64</sub></td></tr><tr><td>T<sub>conj,b</sub> [BJD-2450000]</td><td>8674.28<sup>+5.72</sup><sub>-6.24</sub></td><td>8698.24<sup>+6.49</sup><sub>-11.29</sub></td></tr><tr><td>√e<sub>b</sub> cos ω<sub>*,b</sub></td><td>-</td><td>-0.638<sup>+0.207</sup><sub>-0.114</sub></td></tr><tr><td>√e<sub>b</sub> sin ω<sub>*,b</sub></td><td>-</td><td>-0.078<sup>+0.267</sup><sub>-0.300</sub></td></tr><tr><td>K<sub>365-d, CARMENES</sub> [m s<sup>-1</sup>]</td><td>1.2 ± 0.2</td><td>1.2 ± 0.2</td></tr><tr><td>P<sub>365-d, CARMENES</sub> [d]</td><td>354.58<sup>+7.68</sup><sub>-7.10</sub></td><td>355.21<sup>+7.03</sup><sub>-7.12</sub></td></tr><tr><td>T<sub>0,365-d CARMENES</sub> [BJD-2450000]</td><td>8782.83<sup>+14.95</sup><sub>-14.81</sub></td><td>8782.51<sup>+14.59</sup><sub>-14.82</sub></td></tr><tr><td colspan="3"><strong>Derived parameter</strong></td></tr><tr><td>eccentricity, e<sub>b</sub></td><td>-</td><td>0.48<sup>+0.14</sup><sub>-0.16</sub></td></tr><tr><td>arg. of periapsis, ω<sub>*,b</sub></td><td>-</td><td>-2.30<sup>+5.27</sup><sub>-0.66</sub></td></tr><tr><td>ln Z</td><td>-958.2</td><td>-956.7</td></tr><tr><td>ln Z<sub>1p</sub>-ln Z<sub>0p</sub></td><td>+11.9</td><td>+13.4</td></tr></tbody></table>

**Notes.** (a) The uncertainties of the parameters are calculated as the 16th and 84th percentiles of their posterior distributions.
<table><caption>**Table D.1.** Results for all the different GP-based fits performed on the combined HARPS+CARMENES dataset using the alternative NAIRA and Trifonov et al. (2020) NZP radial velocity time-series.</caption><thead><tr><th>GP kernel</th><th>NAIRA RVs</th><th>Trifonov et al. (2020) RVs</th></tr></thead><tbody><tr><td rowspan="5">Quasi-periodic (QP)</td><td>Circular orbit:</td><td>Circular orbit:</td></tr><tr><td>K<sub>b</sub> [m s<sup>-1</sup>] = 0.94<sup>+0.19</sup><sub>-0.22</sub><br>P<sub>b</sub> [d] = 139.69<sup>+0.72</sup><sub>-0.80</sub><br>T<sub>conj,b</sub> [BJD-2450000] = 8669.18<sup>+7.02</sup><sub>-6.81</sub><br>ln Z = -966.2; Δ ln Z = +1.3</td><td>K<sub>b</sub> [m s<sup>-1</sup>] = 0.91<sup>+0.17</sup><sub>-0.18</sub><br>P<sub>b</sub> [d] = 139.66<sup>+0.70</sup><sub>-0.68</sub><br>T<sub>conj,b</sub> [BJD-2450000] = 8668.56<sup>+6.56</sup><sub>-6.31</sub><br>ln Z = -970.0; Δ ln Z = +2.5</td></tr><tr><td>Eccentric orbit:</td><td>Eccentric orbit:</td></tr><tr><td>K<sub>b</sub> [m s<sup>-1</sup>] = 1.24<sup>+0.23</sup><sub>-0.21</sub><br>P<sub>b</sub> [d] = 140.38<sup>+0.46</sup><sub>-0.43</sub><br>T<sub>conj,b</sub> [BJD-2450000] = 8697.15<sup>+7.54</sup><sub>-13.15</sub><br>e<sub>b</sub> = 0.48<sup>+0.14</sup><sub>-0.13</sub><br>ln Z = -962.9; Δ ln Z = +4.6</td><td>K<sub>b</sub> [m s<sup>-1</sup>] = 1.23<sup>+0.26</sup><sub>-0.22</sub><br>P<sub>b</sub> [d] = 140.35<sup>+0.42</sup><sub>-0.34</sub><br>T<sub>conj,b</sub> [BJD-2450000] = 8699.09<sup>+6.47</sup><sub>-11.74</sub><br>e<sub>b</sub> = 0.52 ± 0.14<br>ln Z = -966.8; Δ ln Z = +5.7</td></tr><tr><td>Circular orbit:</td><td>Circular orbit:</td></tr><tr><td>Quasi-periodic cosine (QPC)</td><td>K<sub>b</sub> [m s<sup>-1</sup>] = 0.84<sup>+0.17</sup><sub>-0.19</sub><br>P<sub>b</sub> [d] = 139.91<sup>+0.72</sup><sub>-0.77</sub><br>T<sub>conj,b</sub> [BJD-2450000] = 8671.24<sup>+6.85</sup><sub>-7.35</sub><br>ln Z = -958.9; Δ ln Z = +2.0</td><td>K<sub>b</sub> [m s<sup>-1</sup>] = 0.81<sup>+0.17</sup><sub>-0.18</sub><br>P<sub>b</sub> [d] = 139.87<sup>+0.74</sup><sub>-0.76</sub><br>T<sub>conj,b</sub> [BJD-2450000] = 8670.79<sup>+6.98</sup><sub>-7.16</sub><br>ln Z = -965.2; Δ ln Z = +3.6</td></tr><tr><td>Eccentric orbit:</td><td>Eccentric orbit:</td></tr><tr><td></td><td>K<sub>b</sub> [m s<sup>-1</sup>] = 1.15<sup>+0.23</sup><sub>-0.20</sub><br>P<sub>b</sub> [d] = 140.40<sup>+0.46</sup><sub>-0.40</sub><br>T<sub>conj,b</sub> [BJD-2450000] = 8699.02<sup>+6.04</sup><sub>-10.78</sub><br>e<sub>b</sub> = 0.50<sup>+0.13</sup><sub>-0.14</sub><br>ln Z = -955.5; Δ ln Z = +5.4</td><td>K<sub>b</sub> [m s<sup>-1</sup>] = 1.14<sup>+0.27</sup><sub>-0.21</sub><br>P<sub>b</sub> [d] = 140.37<sup>+0.38</sup><sub>-0.33</sub><br>T<sub>conj,b</sub> [BJD-2450000] = 8700.17<sup>+5.73</sup><sub>-11.12</sub><br>e<sub>b</sub> = 0.53<sup>+0.14</sup><sub>-0.15</sub><br>ln Z = -962.8; Δ ln Z = +6.0</td></tr><tr><td>dSHO (rotational)</td><td>Circular orbit:</td><td>Circular orbit:</td></tr><tr><td></td><td>K<sub>b</sub> [m s<sup>-1</sup>] = 0.77±0.17<br>P<sub>b</sub> [d] = 140.43±0.88<br>T<sub>conj,b</sub> [BJD-2450000] = 8675.34<sup>+6.13</sup><sub>-6.69</sub><br>ln Z = -960.7; Δ ln Z = +9.7</td><td>K<sub>b</sub> [m s<sup>-1</sup>] = 0.88<sup>+0.17</sup><sub>-0.16</sub><br>P<sub>b</sub> [d] = 140.42±0.68<br>T<sub>conj,b</sub> [BJD-2450000] = 8674.48 ± 5.89<br>ln Z = -967.5; Δ ln Z = +10.7</td></tr><tr><td>Eccentric orbit:</td><td>K<sub>b</sub> [m s<sup>-1</sup>] = 1.19<sup>+0.33</sup><sub>-0.26</sub><br>P<sub>b</sub> [d] = 141.11<sup>+0.41</sup><sub>-0.71</sub><br>T<sub>conj,b</sub> [BJD-2450000] = 8702.57<sup>+5.08</sup><sub>-8.99</sub><br>e<sub>b</sub> = 0.56<sup>+0.14</sup><sub>-0.16</sub><br>ln Z = -958.4; Δ ln Z = +13.0</td><td>Eccentric orbit:</td></tr></tbody></table>

Note: $\Delta \ln Z = \ln Z_{1\text{pl}} - \ln Z_{0\text{pl}}$ is the difference between the Bayesian evidences $\ln Z$ of the models with and without a planetary signal included.

## 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 $\Delta \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 ($\Delta \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 ($\ln Z_{e_b \neq 0} - \ln Z_{e_b=0} \ge 2.4$ in five cases over six). The eccentricity is retrieved with significance within ~3-4$\sigma$. 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.
<figure><img src="image_25.png" /><figcaption>Fig. E.1. $l_1$ periodogram with free offsets and covariance models with the highest cross-validation score.</figcaption></figure>

## 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 HARPS_TERRA and CARMENES-VIS RVs with other methods. We first use the $l_1$ periodogram (Hara et al. 2017), a tool based on a sparse recovery technique called the basis pursuit algorithm (Chen et al. 1998). The $l_1$ 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

$$V_{k,l} = \delta_{k,l} \cdot (\sigma_{RV}^2 + \sigma_W^2 + \sigma_C^2) + \sigma_R^2 \exp\left[-\frac{(t_k - t_l)^2}{2\tau_R^2}\right] + \\ + \sigma_{act}^2 \exp\left[-\frac{(t_k - t_l)^2}{2\tau_{act}^2} - \frac{1}{2} \sin^2\left(\frac{\pi(t_k - t_l)}{P_{act}}\right)\right] \quad (E.1)$$

where $\sigma_{RV}$ is the nominal measurement uncertainty; $\sigma_W$ is an additional white noise jitter term; $\sigma_C$ is a calibration noise term; $\delta_{k,l}$ is equal to one if measurements $k$ and $l$ are taken within the same night, and zero otherwise; $\sigma_R$ and $\tau_R$ define a correlated term to model contributions due to granulation (Cegla 2019) or instrumental effects, as defined in Hara et al. (2020); $\sigma_{act}$, $\tau_{act}$, and $P_{act}$ are the hyper-parameters of a quasi-periodic covariance term to model stellar activity.

For $\sigma_W$, $\sigma_R$, $\sigma_{QP}$, we use the grid of values 0.5, 1., 1.5 m/s. $\sigma_C$ is fixed to 0.5 m/s. $\tau = 3$ or 6 days, $P_{act}$ is fixed to 31.4 days based on the analysis of activity indicators and $\tau_{act} = 60$ days. We try every combination of these values and find that the model with highest cross-validation has $\sigma_W = 1.5$ m/s, $\sigma_R = 0.5$ m/s, $\sigma_{QP} = 1$ m/1, $\tau_R = 6$ days. We use a free offset for each of the three datasets: pre- and post-2015 HARPS and CARMENES. The corresponding $l_1$ 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 \cdot 10^{-5}$, $2 \cdot 10^{-3}$ 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 $l_1$ 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 $l_1$.

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 $\chi^2$ of a linear base model $H$, and the $\chi^2$ of a model $K(\omega, t_0, \tau)$ defined as the linear model of $H$ plus an apodized sinusoid $e^{-(\omega t_0)^2 / 2\tau^2} \cdot (A \cos \omega t + B \sin \omega t)$. The ASP is defined as the $\chi^2$ difference between the two models

$$z(\omega, t_0, \tau) = \chi_H^2 - \chi_{K(\omega, t_0, \tau)}^2 \quad (E.2)$$

Like in the $l_1$ periodogram, our base model includes the offsets of the three datasets. For the apodized sinusoid, we consider a grid of timescales $\tau = 10 \cdot 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(\omega, t_0, \tau)$ as defined in Eq. (E.2), maximised over $t_0$. 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 $\tau = 10 \cdot T_{obs}$). Using $t_{(\tau, \omega)}$ to denote the value of $t_0$ maximising the value of the periodogram (E.2) for a given $\omega$ and $\tau$, we compute the distribution

$$D_z = z(\omega, t_{(\tau, \omega)}, \tau) - z(\omega, t_{(\tau', \omega)}, \tau'), \quad (E.3)$$

with the hypothesis that model $K(\omega, t_{(\tau, \omega)}, \tau, 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 $\chi^2$ distribution, with mean and variance given by an
<figure><img src="image_26.png" /><figcaption>**Fig. E.2.** Left panel: ASP corresponding to different apodization timescales: $τ = 10 · T_{obs}$, $T_{obs}/3$, $T_{obs}/9$ and $T_{obs}/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).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**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).</figcaption></figure>

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.

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.
<figure><img src="image_28.png" /></figure>

## Appendix F: Results relative to RV models including two Keplerians

In this section, we show plots related to the tests performed on the HARPS<sub>TERRA</sub>+CARMENES RVs discussed in Sect. 6.1.2, where we examined the hypothesis of a two-planet system.
<figure><img src="image_29.png" /></figure>
<figure><img src="image_30.png" /><figcaption>**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.</figcaption></figure>