© The Authors 2023

1 Introduction

Since the early 2000s, the HARPS spectrograph (Pepe et al. 2000; Mayor et al. 2003), installed on the ESO-3.6m telescope at La Silla observatory, Chile, has seen numerous breakthroughs in the detection and characterisation of new worlds. Initially with a Guaranteed Time Observation (GTO) programme, and then an ESO Large Programme, the instrument dedicated its time to observe a large sample of stars. This blind planet search was performed on a well-chosen sample of main sequence and quiet stars. Among the principal results of this HARPS survey was the discovery of a new category of exoplanets, the subNeptune and Super-Earth planets on close-in orbits (Mayor & Udry 2008; Lovis et al. 2009). A significant number of detections later suggested that these planets are among the most abundant in the exoplanet population (e.g. Mayor et al. 2011; Lovis et al. 2011; Lo Curto et al. 2013; Bonfils et al. 2013; Astudillo-Defru et al. 2017b; Delisle et al. 2018; Unger et al. 2021). Their over-abundance was strengthened by the Kepler space telescope, which scrutinised a small area of the North sky for 3.5 yr in search for transiting planets (Latham et al. 2011; Borucki et al. 2011; Fabrycky et al. 2014).

The successful HARPS story motivated the development of a ‘HARPS twin’ in the northern hemisphere, HARPS-N (Cosentino et al. 2012, 2014). This high resolution spectrograph is mounted on the 3.6m Telescopio Nazionale Galileo (TNG), located at the Roque de Los Muchachos observatory on the La Palma Island, Spain. The instrument is embedded in a vacuum-controlled tank, ensuring a high level of stability in temperature and air pressure. It covers the wavelength range 383–693 nm, and reaches a spectral resolution of R = 115 000. The instrument was built with two primary science goals. (1) To ensure a synergy with the Kepler (Borucki et al. 2010), then K2 (Howell et al. 2014), TESS (Ricker et al. 2015) and CHEOPS (Benz et al. 2021) missions. This is achieved via a precise RV follow-up of promising transiting candidates. Ultimately, such a strategy provides the community with bulk density measurements of planets, which are essential to precisely characterise the planet population and constrain formation and evolution processes (e.g. Rajpaul et al. 2017; Malavolta et al. 2018; Bonomo et al. 2019; Cloutier et al. 2020; Mortier et al. 2020; Lacedelli et al. 2021). (2) To undertake a blind RV search for low-mass exoplanets, similar to the historical program of HARPS but in the northern sky. To achieve this, the Rocky Planet Search (RPS) program was started with HARPS-N. It consists in a sample of 58 bright stars widely spread in right ascension, many of which presenting low levels of activity and hence amenable for precise planet characterisation (see Motalebi et al. 2015, for a presentation of the RPS program). Unlike the HARPS sample however, we highlight that some targets with moderate activity are also part of this sample, to serve as test cases for the development of tools aimed at mitigating stellar activity (e.g. Cretignier et al. 2020a; Collier Cameron et al. 2021; de Beurs et al. 2022) in parallel with the constant RV monitoring of the Sun (Collier Cameron et al. 2019; Dumusque et al. 2021).

In order to fulfil these scientific objectives and discover low mass exoplanets, stellar activity has to be treated carefully. The RV amplitude induced by small planets on host main-sequence quiet stars is comparable to the contaminating effect of stellar activity. In order to mitigate the latter, we first adapted the observational strategy of the RPS sample. Following the prescriptions of Dumusque et al. (2011), we systematically took two exposures of 15 min each per night, with at least two hours between each exposure. This allows us to mitigate noise emanating from stellar oscillations and granulation on different scales (Dumusque et al. 2011; Chaplin et al. 2019). Coupling this observational strategy with intense monitoring led to the detection of four planets around HD 219134 (Motalebi et al. 2015). It is one of the few systems where a planet was first discovered via the RV technique, and subsequently found to transit its host star (HD 219134 b, with the use of the Spitzer space telescope). Two years later, a second planet was also observed to transit (Gillon et al. 2017).

In addition to the above, significant efforts have been dedicated to develop tools for stellar activity mitigation and RV time-series analysis (e.g. Delisle et al. 2020a,b; Collier Cameron et al. 2021; Cretignier et al. 2021; Hara et al. 2022; de Beurs et al. 2022). Simultaneously, the RPS program recently benefited from new Data Reduction Software (DRS) inspired from the ESPRESSO data reduction pipeline, and further augmented by the removal of some specific HARPS-N systematics via a study of the solar data (Dumusque et al. 2021). The convergence of all these efforts motivated us to revisit the RPS sample, in search for new exoplanet candidates. A first outcome of this work is the validation of HD 79211 b (DiTomasso et al. 2023), an exoplanet in a 24.4 days orbit and with a minimum mass of 10.6 Μ_⊗. More recently, the RV analyses of HD 166620 and HD 144579 demonstrate the potential for HARPS-N to reach sub-metre per second detectability in blind searches (Anna John et al. 2023).

In this paper, we present the results of the analysis of three systems, namely HD 99492, HD 147379, and HD 190007, each of which is already known to host a planet candidate. We aim to present a homogeneous procedure to analyse the data, using the new techniques introduced above. The present paper is organised as follows. Section 2 reports the observations available for each star, and focuses on the estimation of the stellar parameters. Section 3 describes the various tools as part of our strategy to analyse the data. Then, Sects. 4–6 present our analyses of the data and their modelling with planetary Keplerian orbits. Finally, we conclude and discuss these results in Sect. 7.

2 Observations and stellar parameters

2.1 Observations

2.1.1 HD 99492

HD 99492, also denoted as 83 LeoB, was observed with HARPS-N between January 15, 2014 and June 17, 2022 using the observational strategy mentioned above (per night, two 15 min exposures separated by at least 2 h). A total of 202 nightly binned spectra were recorded. For each binned spectrum, we extracted the RV information from cross-correlation of the stellar spectra with a K2 mask. We also derived different stellar activity indicators such as the bisector span, the full width at half maximum (FWHM) of the cross-correlation function (CCF), the S-index and H_α. 178 Keck/HIRES (Vogt et al. 1994) publicly available RV measurements were also taken over a time-span of 6365 days, between January 13, 1997 and June 19, 2014. In the top panel of Fig. 1, we present the RV time-series of the combined HIRES and HARPS-N datasets, including an adjusted offset.

HD 99492 has 71 publicly available HIPPARCOS photometric measurements, and 368 observations from the Automatic Photoelectric Telescope (APT; Henry 1999; Kane et al. 2009) at Fairborn observatory, Arizona (USA). This ensemble of photometric data did not allow Kane et al. (2016) to constrain the stellar rotation period, and is not used in our analyses. Additionally, TESS observed the star in sectors 45 and 46 (i.e. in November and December 2021) with a 2-min cadence.

Data from the third Gaia Early Data Release (EDR3; Gaia Collaboration 2020) constrain the parallax of HD 99492 to 55.062 ± 0.030 masec, leading to an updated distance estimate of 18.16 ± 0.01 pc from the Sun. The star forms a gravitational binary with HD 99491, also named 83 Leo A. With a reported angular separation between HD 99492 and HD 99491 of 28.3 asec, their mutual distance projected on the sky is estimated to 513.9 AU. We note furthermore the existence of 6 CORAVEL RV measurements of HD 99492 taken between the years 1983 and 1993, and 13 RV of HD 99491 gathered between 1980 and 1993, with a mean RV error on the individual measurements of ~300 m s⁻¹. The precision on those data is not significant enough to detect a slope in the RV due to the gravitational influence between the two stars (Halbwachs et al. 2018).

Fig. 1 RV time-series of the three stars analysed in this study. Top: HD 99492 (HARPS-N and HIRES data). Middle: HD 147379 (HARPS-N, SOPHIE, CARMENES and HIRES data). Bottom: HD 190007 (HARPS-N, HIRES and APF data).

2.1.2 HD 147379

HD 147379, also named Gl 617 A, was observed with HARPSN between February 18, 2013 and August 13, 2020. A total of 165 nightly binned spectra were recorded. The same information as for HD 99492 was extracted from the spectra, this time employing a M0 mask. The star benefits also from publicly available observations performed with three other spectrographs. HIRES gathered 30 RV measurements over about 13 yr (from May 15, 2000 to May 21, 2013); the CARMENES instrument (Quirrenbach et al. 2018) obtained 114 measurements over ~600 nights (from January 10, 2016 to September 2, 2017); the SOPHIE+ spectrograph¹ (Perruchot et al. 2011; Bouchy et al. 2013), finally, brings an additional 163 high precision RV measurements obtained over nearly 6 yr (from September 22, 2011 to June 17, 2017). The RV time-series of the combined instruments is presented in the middle panel of Fig. 1. For the sake of clarity, we present a time-series centred on the CARMENES, SOPHIE, and HARPS-N data. The full RV time-series is presented in Fig. A.1.

HD 147379 was also extensively observed in 28 TESS sectors with a 2-min cadence (sectors 14–21, 23–26, 40–41, and 47–60). The total observing window covers 1279 days, from July 18, 2019 to January 18, 2023, and contains 494 593 flux measurements. Also, the star benefits from 103 HIPPARCOS photometric measurements, obtained over a time-span of 2.5 yr. Finally, it was observed by KELT (Pepper et al. 2007) over slightly less than 2 yr, for a total of 3304 measurements.

The Gaia EDR3 estimates the parallax of HD 147379 to 92.877 ± 0.015 masec, which converts to a distance of 10.767 ± 0.002 pc from the Sun. The star has a nearby companion with a similar proper motion, namely Gl 617 B. This gravitationally bound binary presents an angular separation of 64.4 asec, which corresponds to an on-sky projected distance of 693.4 AU.

2.1.3 HD 190007

HD 190007 was observed with HARPS-N between July 21, 2015 and May 19, 2019, for a total of 37 nightly binned spectra. We derived the RV from cross-correlation of the stellar spectra with a K2 mask. These data complement the set of publicly available RV measurements from the HIRES and APF/Levy (Vogt et al. 2014) spectrographs. The former gathered 33 nightly binned RV over more than 16 yr (from June 19, 1998 to September 10, 2014), while the latter obtained 89 nightly binned high resolution spectra over more than 6yr (from July 9, 2013 to October 2, 2019). The RV time-series of the combined instruments is presented in Fig. 1, bottom panel.

HD 190007 was observed with the APT at Fairborn observatory for over 20 yr, gathering a total of 1092 flux measurements. Furthermore, the star was followed up by TESS in sector 54 (in July 2022), with a 2-min cadence.

2.2 Stellar parameters and activity

Table 1 reports on the main stellar parameters gathered from the literature for the three stars. We stress that the uncertainties on the effective temperature T_eff and the stellar mass M_* are intriguingly small. As Tayar et al. (2022) showed, important factors such as the systematic uncertainties on the fundamental stellar properties, and the scatter in the results from different stellar evolution models are often ignored. Due to these factors, they estimate a minimum uncertainty of 2.0% ± 0.5% on T_eff, to be added in quadrature to the reported uncertainties. Regarding the stellar mass estimation, the authors provide maximal fractional offsets between different model grids in the space of luminosity and effective temperature (cf. their Fig. 5). For the three stars of this work, the models variance adds uncertainties of between 5% and 10% on M_*. This additional error budget transposes directly to the mass estimates of exoplanets (a similar uncertainty increase is also applicable on the radius, but this is not relevant for this work). Finally, we also note that the distances reported in Table 1 are obtained from the parallax measurements reported in the Gaia EDR3 catalogue.

2.2.1 HD 99492

HD 99492 is a K-type main-sequence bright star (V = 7.6). Its stellar companion 83 Leo A is a late G-type star (V = 6.5). HD 99492 shows a low activity level, notably expressed by its $\text{log\hspace{0.17em}}\left(R{\prime }_{HK}\right)=-4.93$ (Marcy et al. 2005).

Kane et al. (2016) detect a magnetic cycle on HD 99492, with a periodicity between 3000 and 5000 days². In the timespan of the HARPS-N observations, which is 2270 days, we see a quadratic-like trend in the data. A similar trend is observed in the time-series of our various spectroscopic activity indicators, and we therefore attribute this trend to stellar activity (cf. Fig. A.2).

From the estimation $\text{log\hspace{0.17em}}\left(R{\prime }_{HK}\right)=-4.93$, Marcy et al. (2005) estimate the stellar rotation period to be ~45 days. They used the empirical relations from Wright et al. (2004) for this estimate, but stress that HD 99492 has a B – V index that lies outside of the verified calibration domain. Kane et al. (2016) analysed both the APT and HIPPARCOS photometry, with no conclusive result. From our 202 nightly binned HARPS-N spectra, we detect significant periodic signals at around 40–45 days in the periodograms of different stellar activity indicators, such as the FWHM and S-index (cf. Fig. A.2). They strengthen the previous analyses, and indicate a stellar rotation period between 40 and 45 days.

Finally, we considered the TESS photometric measurements in search for periodic signals. However, we found that the 2-min cadence Simple Aperture Photometry (SAP) fluxes of the two available sectors (45 and 46), which are derived by the TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016), are strongly contaminated by the Moon. Therefore to treat the systematics, we extracted our custom light-curve from the calibrated 2-min cadence target pixel files using lightkurve³. For each sector, the extraction of the light-curve was performed via aperture photometry employing the TESS SPOC pipeline mask. The latter is shifted from the photometric centre of the source, so to limit the contamination from the nearby companion HD 99491, which lies 28″ away from HD 99492 and is hence blended with the latter on the TESS detector. During the light-curve extraction, we corrected the fluxes from the background systematics, and notably we partly accounted for the contamination of HD 99491. This was done by building a matrix of the out-of-aperture pixels (also called a design matrix), and performing a principal component analysis to retrieve the five main trends that we then removed from the pixels inside the aperture. We also retrieved the cotrending basis vectors (CBV) computed by the Pre-search Data Conditioning (PDC; Smith et al. 2012; Stumpe et al. 2014) unit of the TESS SPOC pipeline. These are vectors that contain the most common systematic trends found for each TESS CCD, and which are produced for every sector. They are divided into three categories: the spike CBV contain the short impulsive systematics, the single-scale CBV contain all systematics in a single set of CBV, and multi-scale CBV are spread into three different band-passes. We found that the combination of single-scale and spike CBV is best suited to correct our light-curve from systematics while preserving stellar signals with time-scales of a few tens of days. We detrended the light-curve with these CBV, jointly with the background subtraction mentioned above. We then removed outlier fluxes standing beyond five standard deviations (5σ) of the smoothed light-curve, we normalised the flux data in each sector, and merged the two sectors. Finally, we undertook a periodic signal search in the resulting light-curve, via a generalised Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009). We found a significant broad periodic signal peaking at ~21 days, which is consistent with half of the expected stellar rotation period. We present those results in Fig. 2.

Fig. 2 HD 99492: periodic signal search in TESS sectors 45 and 46. Top: TESS light-curve extracted using lightkurve, and corrected from systematics (2 sectors). Bottom: GLS periodogram computed between 2 and 40 days of the 1-day binned light-curve.

2.2.2 HD 147379

HD 147379, or equivalently Gl 617 A, is a M-type main-sequence star with high proper motion (V = 8.9 mag). It is also part of a gravitationally bound binary with Gl 617 B, which is another M-type main-sequence star fainter than the former (V = 10.6 mag). HD 147379 presents a reasonably low activity index too, with $\text{log\hspace{0.17em}}\left(R{\prime }_{HK}\right)=-4.75$ (Hobson et al. 2018).

A correlation is visible in the HARPS-N spectroscopic data of HD 147379 between the RV and FWHM or S-index. We observe a long-period variation in those time-series, which is attributed to the stellar magnetic activity. The HARPS-N data spans 2734 days, and like HD 99492, does not cover a full magnetic cycle. Only a quadratic trend is observed (cf. Fig. A.3).

Based on an empirical relation between X-ray activity and stellar rotation, Reiners et al. (2018) estimate a stellar rotation period of P_rot = 31 ± 20 days. Furthermore, with the use of the empirical relation of Astudillo-Defru et al. (2017a) between $\log \left(R{\prime }_{HK}\right)$ and P_rot for M dwarfs, Hobson et al. (2018) evaluate the rotation period of HD 147379 to be 28.8 ± 6.1 days. These authors also analysed the HIPPARCOS photometry – 103 measurements over 2.5 yr – but without conclusive results except a hint of signal at ~21 days. Pepper (2018) performed a photometric analysis of HD 147379 with KELT data, and found P_rot = 22 days. All the previous studies are therefore consistent with each other. We analysed the stellar activity indicators time-series of HARPS-N to search for any hint of stellar rotation. The most significant periodic signals were seen in the S-index, which reveal a cycle of 21.5 days. This strongly supports the results of previous analyses on this star.

We analysed the extensive TESS 2-min cadence photometric observations (28 sectors) in search for any photometric periodic variation. To proceed, we extracted the light-curve of each of the 28 sectors following the same procedure as described in Sect. 2.2.1. Notably, we note again that the SPOC pipeline mask conveniently avoids a nearby source on the detector, and we opted for this mask to perform aperture photometry. We performed a periodic signal search in our custom light-curve corrected from systematics, by computing a GLS periodogram. Both the TESS time-series and the periodogram are presented in Fig. 3. We found two significant signals at 21 and 10.4 days, which we associate with the stellar rotation period and its half, respectively. This analysis is consistent with the period found in the spectroscopic observations and in Pepper (2018). We therefore conclude that the rotation period is highly likely to be around 21 days.

2.2.3 HD 190007

HD 190007 is a bright K-type main-sequence star (V = 7. 46). The star belongs to the family of BY Draconis variables (Kazarovets et al. 2003). It is indeed moderately active ($\text{log\hspace{0.17em}}\left(R{\prime }_{HK}\right)=-4.65$) and shows a clear photometric quasi-periodic variation. It is slightly metal-rich, with an age that remains uncertain. However, Burt et al. (2021) suggest that the star is at least 1 Gyr old, based on data from the Kepler mission.

The radial velocities combined from the different instruments, which span more than 20 yr, present a linear trend of ~1.4 m s⁻¹ yr⁻¹. This is in line with the measurement of a statistically significant (S /N > 3) proper motion anomaly between HIPPARCOS and Gaia DR3 mean epoch (Kervella et al. 2022). Both the RV and the HIPPARCOS-Gaia absolute astrometry hint at the existence of a long-period companion, and future observations will help shed light onto this hypothesis.

Both from the spectroscopic stellar activity indicators and the analysis of the available photometry, Burt et al. (2021) report a rotation period of P_rot ~ 29 days. This is in agreement with the previous estimation of Olspert et al. (2018) using Gaussian Process (GP) modelling, which found P_rot = 27.68 days. We observe a similar periodicity in the spectroscopic indicators when we add our HARPS-N measurements to the existing HIRES and APF data (cf. Fig. A.4).

We also analysed the TESS photometric data of this star (sector 54) in search for potential periodicity. We were limited by the small time-span of the observations (~26 days). A simple inspection of the normalised SAP light curve, consisting of 13 238 flux measurements, reveals a long-term variation covering nearly a full cycle over the sector. We present the normalised SAP light curve in Fig. 4. This is in agreement with previous analyses, suggesting a stellar rotation period around 30 days. Yet, the time span of the TESS observations is too short to firmly confirm it.

Fig. 3 HD 147379 TESS photometry. Top: TESS light-curve extracted using lightkurve, and corrected from systematics (28 sectors). Bottom: GLS periodogram computed between 2 and 40 days of the 1-day binned light curve.

Fig. 4 HD 190007 TESS SAP normalised time-series (sector 54).

3 Data analysis tools and general strategy

In this work, we employed a systematic approach to analyse the observations, making use of a number of versatile tools. We aim to demonstrate the ability of these new tools to correct for stellar activity effects and to provide unbiased planetary parameters. Initially, we explored the spectroscopic data – RV and stellar activity indicator time-series – with the Data & Analysis Centre for Exoplanets (DACE). This web-platform, hosted at the University of Geneva, is dedicated to extra-solar planet data visualisation, exchange and analysis⁴. Among many other functionalities, it shows the various spectroscopic time-series, but also the correlation plots between the RV and indicators. Additionally, it offers the possibility to perform interactive fits of the time-series, optionally with Keplerians. Therefore, we systematically used DACE to get a first valuable view of the data, and identify potential stellar activity signatures. As a second step, we initiated in-depth data analysis. We first corrected the HARPS-N data from instrumental systematics and some stellar activity features thanks to the YARARA software (Cretignier et al. 2021). Then, the datasets were investigated in light of the correlated noise model SPLEAF (Delisle et al. 2020b). Below, we briefly summarise these two tools. Finally, we explored the model parameters via a Markov chain Monte Carlo (MCMC) algorithm. For the latter, we employed samsam, a scaled adaptive metropolis algorithm (Delisle 2022).

3.1 YARARA

YARARA (Cretignier et al. 2021) is a post-processing method that aims to improve the spectra extracted by the classical DRS (Dumusque et al. 2021) by removing the extra signatures introduced by the instrument and the Earth atmosphere, leading ultimately to an improved RV precision. The method involves correcting a residual spectral time-series with respect to a master spectrum (built to be free from the systematics), where the corrections are made of multi-linear regressions in the wavelength domain. The spectra are corrected for: cosmics, tellurics, ghosts, fringing, instrumental defocus, activity and residuals outliers. The cleaned RV are extracted with a CCF using a tailored line selection based on the master spectrum itself (Cretignier et al. 2020a) in order to increase the extraction of the Doppler information content, whereas the merged-order 1d spectra are continuum normalised using RASSINE (Cretignier et al. 2020b). A few adjustments were implemented in the code in order to increase the orthogonality between the stellar activity and instrumental signatures. The main reason for this update was to be able to introduce back the stellar activity component at later stage (see Sect. 3.2 for a justification). In the old version of the code, both effects were simultaneously fitted using the S-index, the CCF FWHM and the CCF contrast. Those two moments contain redundant stellar activity information. The new version of the pipeline uses an improved version of the S-index using CaIIH&K lines better corrected for ghosts, whereas the fitted CCF moments are now orthogonal to that component thanks to a Gram-Schmidt orthogonalisation algorithm (see Appendix C).

YARARA itself is not able to simultaneously fit in the wavelength domain the planetary signals and the systematics. However, since most systematics are fixed in the terrestrial restframe, the code assumes that the planetary signals are destroyed (or strongly mixed) in the wavelength domain once such a change of rest-frame is performed. In Cretignier et al. (2021), the authors show that injected planetary signals are not absorbed more than 20% (and only for a planet close to a 1-yr harmonic). In order to further avoid such absorption, large RV amplitude signals can be pre-fitted by shifting the spectra according to a pre-determined Keplerian solution. In this work, that preliminary Keplerian solution was obtained via a RV fit with DACE which matches the orbital elements of the known planets. This Keplerian solution was then re-injected in the RV obtained on the residual spectra, and an improved Keplerian solution was obtained. This process was called in an iterative way until convergence was reached, which usually occurs after a single iteration.

Once spectra are corrected for systematics, new proxies orthogonal to a pure Doppler shift can be extracted in the timedomain (Cretignier et al. 2022). Under the assumption that the deformations of the line profiles are mainly driven by the line profile itself, the residual spectra δf (λ) = f(λ) - f₀(λ) can be expressed rather as δ f (∂f₀/∂λ, f₀), with f₀ the master spectrum of the star – an object that the authors called a ‘shell’. An advantage of that space is that a pure Doppler shift possesses a unique signature that can be fitted out analytically with a first order approximation. Other distortions, orthogonal to the pure Doppler shift, can then be extracted with a principal component analysis and their associated time-coefficients are used in a multi-linear model to correct for the RV time-series. A cross-validation algorithm is used to determine the number of components to select by randomly rejecting 10% of the observations. Such a strategy was shown successful to remove the rotational period of stars moderately active and/or potential remaining instrumental systematics uncorrected by YARARA.

This additional correction process (Cretignier et al. 2022) is subsequently designated as YARARA-Shell, as opposed to YARARA-V1 which consists of YARARA (Cretignier et al. 2021) further augmented with the few adjustments described above to increase the orthogonality between the stellar activity and instrumental signatures.

3.2 SPLEAF

Correlated noise, or Gaussian Process (GP) models have been the object of increasing research among the exoplanet community over the past 10 yr. They account for physical processes that might be poorly understood or poorly constrained, such as stellar activity (Haywood et al. 2014). Indeed, instead of a deterministic model, the GP regression technique aims to parametrise the covariance between the measurements. There exist different types of GP, which differ from each other by their kernels, that is by the functional form of the covariance matrices. SPLEAF (Delisle et al. 2020b) refers to a class of semi-separable covariance matrices, which builds on the benefits of the celerite kernel (Ambikasaran et al. 2015; Foreman-Mackey et al. 2017). SPLEAF can take into account the calibration noise, and due to the semi-separable form of the covariance matrix, it provides low computation costs. Once the GP regression is properly performed, it should account for anything that is not due to the deterministic model nor the measurements noise. Yet in practice, there might be risks of over-fitting, and the planetary signals may be (partly) absorbed by the GP. To limit the risk, Rajpaul et al. (2015) proposed to train the GP simultaneously on the RV and on activity indicators time-series, by modelling the activity-induced RV signals as linear combinations of the GP (G(t)) and its derivative G(t) (Aigrain et al. 2012) $$\text{Δ}RV=V_{c}G\left(t\right)+V_r\dot{G}\left(t\right)$$(1) $$\text{Δ}\alpha =L_{c}G\left(t\right)$$(2) $$\text{Δ}\beta =B_{c}G\left(t\right)+B_r\dot{G}\left(t\right)$$(3)

for V_c, V_r, L_c, B_c, B_r free parameters, and α and β referring to some activity indicators. Eq. (1) was first proposed by Aigrain et al. (2012) to account for the effect of stellar spots on the RV variations, where G(t) captures the suppression of the convective blueshift in the spots, and $\dot{G}\left(t\right)$ provides the effect of stellar rotation. Rajpaul et al. (2015) extended this framework to activity indicators (Eqs. (2) and (3)), where the use of one equation or another depends on the nature of the indicator. For instance, the S-index is a proxy for the proportion of the visible stellar disk covered by active regions, and is therefore accounted for by Eq. (2). On the other hand, the bisector span indicator, which informs about the asymmetry of the spectral lines, also depends on the radial velocity of the stellar surface at the location of the spots, and is hence described by Eq. (3). The planetary signal being only in the RV time-series, it is less likely to be absorbed by the GP if the latter is simultaneously trained on one or several activity indicators. As a drawback though, the computation cost is significantly larger, with a covariance matrix of size 2n × 2n where n is the number of measurements and in the case where only one activity indicator time-series is used in parallel of the RV.

Delisle et al. (2022) generalised the SPLEAF model to account for different time-series (potentially taken at different times) as modelled by Eqs. (1)–(3), while insuring a computation cost that scales linearly with the total number of measurements. In this work, we used this generalised SPLEAF correlated noise model, which is publicly available⁵. The kernel we employed for the quasi-periodic part is an approximation of the squared-exponential periodic (SEP) kernel, where the covariance has the following expression: $$k\left(\text{Δ}t\right)={\sigma }^2\text{\hspace{0.17em}exp\hspace{0.17em}}\left(-\frac{\text{Δ}{t}^2}{2{\rho }^2}-\frac{{\text{sin}}^2\left(\frac{\pi \text{Δ}t}{P}\right)}{2{\eta }^2}\right).$$(4)

It consists of a sinusoid on top of a decreasing exponential, with Δt being the time interval between two measurements. Therefore, this kernel correlates one measurement with the others according to this functional form. The hyperparameters σ, ρ, P and η describe respectively the amplitude of the correlated noise, the rate of exponential decay of the correlation, the period of variation of the sinusoid and the length-scale of the periodic component, that is its ability to capture rapidly varying features (smaller values pointing towards sharp variations). These hyperparameters are constrained from the various time-series. SPLEAF develops the periodic component of the SEP kernel in series, and keeps the two strongest harmonics. In practice, we also added a white noise term to each time-series, which consists in an additional term on the diagonal of the correlation matrix.

In this work, we systematically applied the generalised SPLEAF model - denoted as the SPLEAF model hereafter – on the RV and S-index time-series. As was described above, we find in the S-index time-series of the three stars traces of the stellar magnetic cycles and stellar rotation. Therefore, the hyper-parameters of the correlated noise model, which are estimated from the analysis of these different datasets, are expected to converge to values informing us about the stellar activity. We mostly set uninformative priors on the model parameters, that is priors following a uniform law with wide boundaries. In Table 2, we synthesise the priors we used in most of our analyses. We explicitly note in the text when certain priors are further constrained.

We express some caution about using the YARARA dataset with a correlated noise model. Indeed, we noticed a convergence issue of the fit with any dataset that would already be (partially) corrected for the stellar activity features. The reason is that while the activity indicator time-series are still left intact, the RV time-series are corrected for the trends observed in the former. Therefore, the correlated noise model trains successfully on the activity indicators, but struggles to transpose the results to the RV time-series. In conclusion, the correlated noise model does not provide satisfactory results with YARARA-V1 and YARARA-Shell. Nevertheless, YARARA also corrects the various time-series for the instrumental systematics. The use of those datasets is preferred, so to avoid observing spurious periodic signals. Therefore, in the following analyses making use of SPLEAF, we used a modified version of YARARA-V1 where the stellar activity was re-injected into the RV time-series, but not the instrumental systematics (see Appendix C for an explanation of the process to separate the various components). In other words, we used an alternative version of YARARA where only the instrumental systematics were corrected for. That allows for the convergence of the fit with the correlated noise model.

Fig. 5 HD 99492: YARARA-V1 S-index dataset. Top: S-index time-series. Middle: GLS periodogram of that time-series (k0 stands for no Keplerian in the model). Bottom: GLS periodogram of the residual S-index, after subtracting a Keplerian model with P ~ 3000 days (k1). The red bands locate the periodicities seen in this activity indicator.

4 HD 99492 data analysis

HD 99492 is known to be orbited by a planet. Marcy et al. (2005) first announced the detection with 35 HIRES RV measurements, and reported an orbital period of 17 days and a minimum mass of ~25 M_⊕. The solution was refined by Meschiari et al. (2011) with an additional 58 HIRES velocities, and the authors proposed another planet candidate at an orbital period of ~5000 days. A few years later and with 130 HIRES RV, the latter signal was shown to be related to stellar magnetic activity (Kane et al. 2016).

We added our 202 nightly binned HARPS-N spectra processed with YARARA-V1, and performed an exploratory analysis on DACE. First, we focused on stellar activity signatures. In this dataset, the S-index time-series displays a significant long-term trend that peaks at a periodicity of ~3000 days in a GLS periodogram. The latter signal is associated with the magnetic cycle described in Kane et al. (2016). Figure 5 presents this time-series (top plot) and the corresponding periodogram (middle plot). At short periods, we observe a strong signal at 42.9 days, which appears once we fit the long-period variation with a Keplerian. This is illustrated in the bottom plot of Fig. 5. Furthermore, we also detect its 1 yr aliases at 47.3 and 38.4 days. The S-index indicator reveals a notable correlation with the RV, with a correlation coefficient of r = 0.50. We detect a similar 42.9 days periodicity in other activity indicator time-series such as the CCF contrast and CCF FWHM. We also systematically identify the 1 yr aliases of this signal. Therefore, all these observations constitute strong evidence for the 42.9 days signal to be associated with the stellar rotation period. This estimate is consistent with previous works (cf. Sect. 2.2).

The HIRES publicly available data also contain the S-index activity indicator. In this dataset, we observe a clear quasi-periodic variation covering several magnetic cycles of the star, whose correlation with the RV is weaker than in the HARPS-N dataset. After fitting the most significant periodic signal in the S-index with a Keplerian (which converges to ~3000 days), we did not find any significant residual periodicity at ~43 days. In fact, there is no remaining significant signal below 800 days. We searched for the presence of the stellar rotation in the other activity indicator available (i.e. the H_α index). Again, we did not find any significant signal. As a result, we note that the rotation of the star does not present a measurable signature in the HIRES data. A potential cause of this non detection is the larger scatter observed in the S-index time-series of HIRES combined with a more sparse sampling. Finally, combining HARPS-N and HIRES datasets, we set tighter constraints on the period of the stellar magnetic cycle. With now an extended baseline of 9286 days, we constrain the main period of the S-index time-series to around 3300 days.

Taking the above information into account, we then undertook a search for planets. We performed the latter using different models and datasets, so as to compare several approaches. In the first, we used the HARPS-N data as derived from the DRS together with the HIRES dataset. In a second approach, we analysed the YARARA-V1 dataset of HARPS-N together with HIRES, and included a correlated noise model. As a third approach, we repeated this procedure on the YARARA-V1 data only. Finally, in a fourth approach we analysed the YARARA-Shell data without a correlated noise model. We describe our planet search analyses below.

4.1 Approach 1: DACE - HARPS-N DRS + HIRES

As a first preliminary analysis, we looked for periodic signals in the nightly binned RV time-series of the combined HARPS-N + HIRES datasets, using the DRS version of the HARPS-N data. While we noticed inconsistencies in the H_α time-series of the HIRES data, the S-index is consistent throughout the entire timespan. We selected this activity indicator to detrend the RV via the inclusion of a scaling parameter. Additionally, we included a quadratic drift to remove what remains of long-term variations. After adding these components in the model, we searched for periodic signals in the periodogram. The signal of the planet at 17 days is very strong, with a FAP of 7.1 × 10⁻⁷⁷. We included a Keplerian at that period into our model, and computed the periodogram of the residuals from the fit. This periodogram shows a very significant second signal. It has a period of 95.5 days and FAP of 2.1 × 10⁻²¹, and is not associated with any of the activity indicators. This favours the planet hypothesis, and we fitted that signal with another Keplerian. A third signal is revealed in the new periodogram of the residuals, with a period of 13.9 days and FAP of 0.23%. Not only its FAP is above our defined detection threshold of 0.1 %, but also, this period is about a third of the expected stellar rotation period. Hence, it has to be interpreted carefully.

We explored the two Keplerians model via a MCMC algorithm, which allows us to constrain the planet parameters. The parameters of HD 99492 b are in accordance with previous publications, with results converging to P_b = 17.0492 ± 0.0007 days and K_b =7.30 ± 0.23 m s⁻¹. Concerning the new planet candidate, we found a moderate eccentricity e_c = 0.237 ± 0.080, while the orbital period and the RV semi-amplitude are estimated to P_c = 95.373 ± 0.050 days and K_c = 3.08 ± 0.27 m s⁻¹.

Fig. 6 HD 99492: periodograms of the combined YARARA-V1 and HIRES RV time-series. The periodograms are computed from a fit of the time-series with a detrend + drift model. Top: no Keplerian. Middle: periodogram of the residual time-series, after removing a Keplerian at ~17 days. Bottom: periodogram of the residuals after removing Keplerians at ~17 and 95 days.

4.2 Approach 2: SPLEAF - HARPS-N YARARA-V1 + HIRES

We repeated the process highlighted in Sect. 4.1 with HIRES and the YARARA-V1 data reduction of HARPS-N. Again, the signals at 17 and 95.5 days are unambiguous. Figure 6 illustrates our successive fitting process. The top plot presents the periodogram of the RV time-series (the highest peak has a FAP of 1.1 × 10⁻⁸⁰), while the middle plot shows the periodogram of the residuals after fitting the RV with a Keplerian at 17 days. The highest peak in this middle plot, located at P = 95.4 days, has a FAP of 3.2 × 10⁻²². The bottom plot displays the periodogram of the residuals after performing a two Keplerian fit, including one at 17 days and one at 95.4 days. Now, the remaining signal at 13.9 days is no longer detected in the latter plot, but we observe another signal at 14.5 days with a FAP of 0.13%. Its significance is still below our threshold of 0.1% in FAP. Both the low significance and the different period detected for this third signal, in addition to its periodicity of about a third of the expected stellar rotation period, cast doubt on its planetary origin.

In order to shed light on the origin of this third signal, we scrutinised the HIRES and HARPS-N datasets in a correlated noise model via SPLEAF (cf. Sect. 3.2), using both the RV and S-index time-series. At this stage, we used the alternative YARARA-V1 dataset where the stellar activity was injected back into the RV. We did not set any constraining prior on the period of the correlated noise, nor on any other parameter (cf. Table 2). However, we initialised the correlation period to 43 days, and the decay time-scale to 500 days. The large value for the latter is aimed at mitigating the magnetic cycle, which takes place on a long time-scale. In order to facilitate the fit convergence, we fitted the correlated noise model parameters on a step-by-step strategy, successively adding an additional parameter. Then, we progressively incorporated Keplerians until no more significant signal was observed. In addition to these parameters, we also included an offset for each instrument and a white noise term, both included for each time-series – namely the RV and S-index time-series. Finally, we added a RV linear drift in the model, which provides the best results in terms of signal significance and fit convergence compared to a quadratic drift or no drift. As a result of our fit, we find a first significant periodic signal at 17.05 days, with a negligible FAP of 4.5×10⁻⁵². After adding a Keplerian in the model and proceeding with a new fit, we detect another significant signal with a period of 95.3 days and FAP = 2.5 × 10⁻¹⁸. Therefore, we incorporated a second Keplerian in our model, and performed a new fit. In the periodogram of the residuals, there remains a significant power around 3500 days which we attribute to a remnant of the stellar magnetic cycle, while the decay time-scale of the correlated noise converged to 194 days. Concerning the correlation period, it peaks at 44.8 days, which is compatible with the stellar rotation period. Then, we explored the parameter space of the two-Keplerian model with a MCMC algorithm. We performed 1M iterations, and obtained an effective sample of 4288 independent solutions. The period of the correlated noise is estimated to be 44.9 ± 0.5 days. The outermost planet, at an orbital period of 95.4 days, has an estimated RV semiamplitude of 2.7 ± 0.2 m s⁻¹, which corresponds to a planet minimum mass of 17.4 Μ_⊕. Its eccentricity is moderate, with e_c = 0.112 ± 0.086.

4.3 Approach 3: SPLEAF – HARPS-N YARARA-V1

In a third approach, we employed again SPLEAF for a correlated noise model but analysed only the HARPS-N data – YARARA-V1 with activity injected back into the time-series. We undertook the same process as described above for the second approach. Once again, we find two significant signals at 17 and 95.5 days. Their FAP are 3.0 × 10⁻³⁵ and 2.5 × 10⁻²⁰ respectively. A third residual signal with a period of 14.7 days just reaches our detection threshold, and has a FAP of 0.05%. It has a small semi-amplitude of 1.0 m s⁻¹. This period is different from the previously found signals at 13.9 and 14.5 days, but is also compatible with a third of the stellar rotation period. We note that the period of the correlated noise is estimated to 43.1 days. We remind the reader that SPLEAF approximates the periodic component of the SEP kernel by a development to the second order. To investigate further the stellar rotation origin for the 14.7 days signal, we developed the periodic component to a higher order, so as to account for the third harmonic of the stellar rotation period. We applied this refined kernel to the same dataset. After the removal of two Keplerians at 17 and 95.5 days, the residual signal at 14.7 days is still significant with a FAP of 0.07%. Developing the kernel up to one additional order so as to entail the fourth harmonic yielded similar results. The residual 14.7 days signal has a FAP of 0.08% with this refined development. However, in all cases, the period of the correlated noise unambiguously converges towards the stellar rotation period. Therefore, developing the periodic component of the kernel to a higher order should absorb the harmonics of the stellar rotation. Consequently, the survival of the 14.7 days signal irrespective of the kernel development suggests this signal is not a direct outcome of the third harmonic of stellar rotation. Instead, its origin is likely different.

We undertook an MCMC exploration of the model with three Keplerians. As the RV semi-amplitude of the third signal is small and in order to facilitate the MCMC convergence, we fixed the eccentricity and argument of periastron of the third Keplerian to 0. We performed 1M iterations, and obtained a total of 4856 independent solutions. Out of this exploration, we find RV semi-amplitudes of K_b = 6.99 ± 0.17, K_c = 2.83 ± 0.20, and K_d = 0.94 ± 0.17 m s⁻¹. The orbital eccentricity of planet c is now smaller: e_c = 0.052 with a 68.27% confidence interval of [0.019, 0.102].

We also explored the parameter space of a model with only two Keplerians, given the doubt that we cast on the origin of the 14.7 days signal. We performed a new MCMC exploration of the model with two Keplerians, again with 1M iterations. This time, we obtained a sample of 5440 independent solutions. The correlated noise parameters display similar values at the end of the MCMC exploration, compared to the model with three Keplerians. The same is observed for planets b and c RV semi-amplitudes, and the orbital eccentricity of planet c is once again constrained to be small. For this two Keplerian model, we estimate a Bayesian information criterion (BIC) of −327.1, while the BIC of the three Keplerian model amounts to −349.5. Therefore, the BIC indicates a slight preference for the model with three Keplerians. As such, we further tested the planetary nature of the 14.7 days signal, by investigating its resistance to the removal of RV measurements. If the signal was due to a planet, its power would monotonously decrease with the number of removed observations. Instead, a signal of stellar origin would react irregularly to the suppression of RV measurements. We performed different tests, following various patterns of data removal (random, lowest S/N, highest S/N, data in quadrature of the 14.7 days signal) and removing up to 10% of the time-series. We found that while the signals at 17 and 95 days successfully passed our tests, the significance of the 14.7 days signal did not decrease monotonously with the number of removed observations. Instead, its power in the periodogram evolved irregularly with the number of removed observations. From this analysis, we conclude that the current observations do not support a third planetary signal. On the other hand, the resilience of the 14.7 days signal to the degree of development of SPLEAF indicates that it does not emanate from the third harmonic of stellar rotation. Instead, it could be the result of a combination between stellar rotation (or its second harmonic) and the spectral window of the observations.

Fig. 7 HD 99492 c posterior distributions for our different models and datasets. The grey, blue, orange, and green distributions refer to our first, second, third, and fourth approaches to derive the posteriors, respectively. In the case of the third approach, we further show the distributions resulting from two different models. The plain line distribution was obtained from a model with two Keplerians, while the dashed line corresponds to a model with three Keplerians.

4.4 Approach 4: HARPS-N YARARA-Shell

Finally, we investigated the YARARA-Shell dataset which, as compared to YARARA-V1, further corrects from stellar activity features. We designed a model that includes a white noise term, but does not include correlated noise. After an iterative search, we found two significant signals at the periods already mentioned above. The first signal, at 17.05 days, has a FAP of 5.1 × 10⁻⁶², while the second signal at 95.3 days has a FAP of 3.5 × 10⁻²⁸. In the residuals of this two planet model, we did not find any significant period. The 14.7 days signal has now a FAP of 9.7%. This observation further comforts us in the nondetection of a third Keplerian signal with the current dataset. The analysis of YARARA-Shell, following the analysis of YARARA-V1 with SPLEAF, hereby demonstrates the power of using complementary analysis techniques to shed light onto the nature of periodic signals.

We explored the two-planet model via the MCMC algorithm introduced above. We performed 1M iterations, and obtained a set of 4095 independent solutions that constitute our posterior distribution. While the RV semi-amplitude of the outer planet Kc is very similar to its estimate obtained using approach 3, the orbital eccentricity of that planet converges towards a larger value of about 0.1.

In Fig. 7, we synthesise the results from the four approaches described above. It shows the posterior distributions obtained with the different approaches, and projected on the orbital parameters of planet c. Additionally, we also distinguish between the two and three planets fits performed in the third approach. We conclude that the inclusion of the HIRES RV measurements increases the estimation of the eccentricity of planet c. This is probably due to a less efficient mitigation of stellar activity effects on the HIRES measurements. Indeed, as was already mentioned, we did not find a clear periodicity of ~43 days in the S-index of HIRES. We note that the analysis without the HIRES data provides the best estimations on the planet parameters. Furthermore, using the YARARA-V1 data alone, the estimates we obtain for the parameters of planet c are very similar in the models with two or three Keplerians. The third periodic signal at 14.7 days is absent from the YARARA-Shell data. This latter dataset, however, leads to a larger estimate of the eccentricity ec. As pointed out in Hara et al. (2019), complex noise patterns - if not modelled - can boost the eccentricity estimations. This additional eccentricity induced by stellar activity likely also drives the shift observed in the distribution of ω_c from YARARA-Shell. Indeed, small orbital eccentricities do not constrain well the argument of periastron, as illustrated with the orange distributions. Any additional eccentricity would hence carry a significant influence on the distribution of ω. In conclusion, we favour the approach leading to the smallest eccentricities. We opt for the approach using the YARARA-V1 dataset only together with a correlated noise analysis, namely approach 3 with a two-planet model. These results are synthesised in Table 3. They suggest that planet c has a minimal mass of 17.9 ± 1.3 Μ_⊕, which is equivalent to the mass of Neptune. So far, few planets in that mass range have been found with large (>50 days) orbital periods. We further discuss this result in Sect. 7. We present in Fig. 8 the plot of the orbits in the Cartesian space, and the RV phase-folded on the period of the two planets.

Fig. 8 HD 99492 orbital solutions for planets b and c derived from the HARPS-N YARARA-V1 dataset. Left: Cartesian representation of the planetary system for a sample of 205 MCMC solutions. The star is in red, the orbit of HD 99492 b in blue and the orbit of HD 99492 c in orange. Additionally, the green area marks the beginning of the conservative habitable zone, which continues outwards. This was calculated from the prescription of Kopparapu et al. (2013). Right: RV measurements phase-folded on the best fit Keplerian solutions, together with the model curves.

4.5 Independent confirmation using TWEAKS

The HARPS-N CCF were independently analysed for planetary signals using TWEAKS. This pipeline was particularly designed for attaining a sub-m s⁻¹ detection threshold at long orbital periods, by combining the wavelength-domain and time-domain stellar activity mitigation (Anna John et al. 2022, 2023). We first conducted a blind search of the radial velocities, using a model with up to five unidentified Keplerian signals. For this, we used the KIMA nested-sampling package (Faria et al. 2018). Using SCALPELS (Collier Cameron et al. 2021), which involves doing a principal-component analysis on the autocorrelation function of the CCF, time-domain activity-decorrelation vectors were produced. These basis vectors were then used for the spectral line-shape decorrelation (Collier Cameron et al. 2021) in KIMA, as Anna John et al. (2022) reported that de-trending the RV for line shape variations using the SCALPELS basis vectors yields a model that is significantly better than a model that does not account for these stellar activity signatures.

The joint posteriors show clear detection of two Keplerian signals at orbital periods 17.05 days and 95.2 days. We conducted a False Inclusion Probability (FIP) analysis (Hara et al. 2022) in frequency space, with the bin size set to the Nyquist frequency resolution over the entire data duration to search for multiple planet signals simultaneously. While a strong FIP of 0.13 is found for the 17.05 days signal, the 95.2 days signal is detected even more strongly with a FIP of 0.04; in other words, 96% of the mutually independent KIMA posterior solutions favour a planet detection at an orbital period of 95.2 days. We present these results in Fig. A.5. The RV semi-amplitudes obtained from this analysis were combined with the stellar mass of 0.85 Μ_⊙, leading to a minimum mass determination of 27.13 ± 1.18 Μ_⊕ and 19.85 ± 2.19 Μ_⊕ respectively. We also validated these detections by confirming the coherency of the signals across different data subsets. Unlike the signals occurring from sampling patterns, stellar activity or aliases, both the 17.05 and 95 days signals are strongly detected in all the subsets with sigma ≥9σ. Finally, as opposed to these two planetary signals, the signal at 14d is not consistent in time. It is only detected in the second half of the data. This further motivates us to discard this signal as a potential additional planet.

4.6 Transit search

We analysed the TESS photometry in search for transit signals. To treat the instrumental systematics, we extracted our custom light-curve from the 2-min cadence pixel files with lightkurve, following the process explained in Sect. 2.2.1. To aim at the transit search, we detrended our light curve from a combination of multi-scale and spike CBV. Indeed, we found that this model performs better at minimising the Combined Differential Photometric Precision (CDPP) metric, which is a measure of the remaining scatter in the light curve expressed in parts-per-million (ppm). Our corrected light curve presents a CDPP of 123 ppm.

As a second step, we undertook the removal of the remaining stellar systematics on our corrected light-curve. We applied a detrending via splines fitting using keplersplinev2⁶, and letting free the parameter describing the time-scale of variation of the spline, only imposing a lower boundary of 0.5 day. The fit converged to a time-scale of 0.74 day, and the modelled trend was removed from our light-curve. The resulting detrended light-curve is presented in Fig. B.1 (top panel).

We undertook a multi-transits search in this detrended light-curve. It was carried out with the transitleastsquares software⁷. This search was done via the computation of a box least square (BLS) periodogram, with an orbital period search comprised between 2 and 20 days. We do not find any significant transit signal in this period range. From the estimated minimum mass of planet b, the mass-radius relationship for rocky planets⁸ from Otegi et al. (2020) and the stellar radius reported on exo-fop⁹, we estimated a transit depth of ~1000 ppm. We do not observe any hint of a transit signal at the period and phase of planet b estimated from our RV analysis. Concerning planet c, the expected transit time occurred at the beginning of sector 45. Under the hypothesis of a rocky composition, the planet would induce a transit depth of 850 ppm on its parent star, which we do not see in the data. For both planets, the hypothetical transit depths are large enough compared to the light-curve noise level, and they would be detected. Therefore, we rule out the transiting configuration of HD 99492 b. Concerning the outer planet, the TESS observations do not cover all the orbital phases, and given the uncertainty on the transit time it is possible that a potential transit was missed. However, it is unlikely given the low transit probability of such a long-period planet (~1%). Details of the transit searches can be found in Fig. B.1.

We also analysed the TESS sectors 45 and 46 in search for potential single transit features. To proceed, we fitted our detrended light curve based on two distinct models: no planet and one planet models. We compared the posterior probabilities of the models using the true inclusion probability (TIP) detection criterion (Hara et al. 2022), which provides us with a rigorous transit detection threshold. A TIP of 1 would favour the one planet model with a probability of 100% (TIP = 1-FIP). This framework was first designed for RV data, and was applied later on to transit analyses (Hoyer et al. 2022; Ehrenreich et al. 2023; Wilson et al., in prep.). To aid convergence, we performed the TIP analyses on one-day slices of the TESS sectors. We ran two sets of analyses employing different priors to first search for any transit signal and subsequently to probe the transiting nature of HD 99492 b: wide planetary priors in one case, priors constrained on the HD 99492 b planet parameters as reported in this study in the other case. We did not find any preference for the one-planet model as the TIP never exceeds 0.5 across the entire light curve (cf. Fig. B.2). Therefore, we conclude on the absence of transit signals in this data set.

5 HD 147379 data analysis

HD 147379 was first found to harbour a planet by the CARMENES team (Reiners et al. 2018). With the help of 114 CARMENES RV measurements and an additional 30 HIRES RV, they detected a planet with a minimum mass of ~25 Μ_⊕ in a 86.5 day orbit. This detection was independently confirmed by the SOPHIE team, with a dataset of 163 SOPHIE RV measurements (Hobson et al. 2018). A noticeable fact about HD 147379 b is that it lies in the conservative habitable zone of its host star, as defined by Kopparapu et al. (2013).

5.1 RV analysis

We analysed this system using exclusively the 165 nightly binned HARPS-N RV measurements. As with HD 99492, we first looked in the modified YARARA-V1 dataset for activity signatures - that is YARARA-V1 where the stellar activity was injected back in the RV. Besides the long-term trend due to the stellar magnetic cycle, the most significant signal is detected in the S-index indicator and peaks at 21.6 days. Furthermore, we observe a weak correlation between the RV and the S-index measurements, with a correlation coefficient of R = 0.39 (cf. Fig. A.3).

In a second step, we undertook an analysis of the RV in search for planetary signals. We employed a correlated noise model with SPLEAF, using both the RV and S-index time-series. We also included in the model an offset, a white noise term and a linear drift term¹⁰. Wide uniform priors were injected, and we initialised the correlated noise with a period of 22 days and a large correlation decay time-scale of 500 days. After performing a first fit – again through a successive addition in the model of the different correlated noise parameters – we find a significant signal at 86.5 days in the periodogram of the residuals (Fig. 9). We included the latter in our model as a Keplerian, and fitted again the RV time-series with the new model. In the periodogram of the residuals, we do not find additional significant signals and the correlated noise period is estimated to 21.8 days, which is the expected stellar rotation period. The most prominent signal that we find in the residuals has a period of 7.1 days, which is about a third of the estimated stellar rotation period, and FAP = 3.8%. This signal is pushed down to a FAP of 23% when we use the augmented version of SPLEAF, which develops further the periodic component of the kernel to account for the third harmonic. This result further suggests the stellar rotation origin of the 7.1 days signal. Furthermore, we do not find any signal around 500 days, which was reported by Hobson et al. (2018) using the SOPHIE data alone. Therefore, we stop our planet search and consider a model composed of one Keplerian only. The exploration of this model was performed with a MCMC using 500k iterations, leading to a final sample of 3298 independent solutions. The results are presented in Table 3. The stellar rotation period that we estimate (P_GP = 21.9 ± 0.4 days) is consistent with previous investigations, and our correlated noise model converges to a correlation decay time-scale of 28.3 ± 4.4 days. The minimum mass of HD 147379 b is found to be smaller than both the estimates of Reiners et al. (2018) with ${24.7}_{-2.4}^{+1.8}\text{}{M}_{\oplus }$ and Hobson et al. (2018) with 28.6 ± 1.5 M_⊕ – our result suggests 21.6 ± 1.1 Μ_⊕, or about 1.26 times the mass of Neptune, and was derived from a careful mitigation of stellar rotation features. The bottom plot of Fig. 9 presents the RV measurements folded on the orbital period and phase of HD 147379 b.

As a comparison, we searched for planetary signals with HARPS-N using the YARARA-Shell post-process. After deriving the cleaned RV, we computed a periodogram which unambiguously revealed the 86.5d planet (FAP = 3.3 × 10⁻³³). After fitting the latter with a Keplerian model and computing the periodogram of the residual RV, we did not find any signal that reaches our detection threshold of FAP = 0.1 %. The strongest peak stands at a period of 12.3 days, and has FAP = 0.5%. Because of this too large FAP and the weakness of the signal in V1-SPLEAF (FAP > 10%), we did not retain this signal as a planetary candidate. Further observations are needed to shed light on this signal.

Finally, we also analysed the entire dataset composed of the HARPS-N (DRS), CARMENES, SOPHIE and HIRES RV. This ensemble consists of 458 RV measurements¹¹. However, the CARMENES and SOPHIE public data do not contain spectroscopic activity indicators, but only the RV measurements. Therefore, we could not carry a reliable modelling of stellar activity. We only applied a linear drift fit, in order to account for an observable trend. We computed the GLS periodogram of the residuals on this combined dataset, and found the 86.5 days planet as the most significant peak. After its subtraction from the time-series and the computation of the updated periodogram, we found several periodic signals with a FAP below 0.1% at 10.6, 12.3 and 21.4 days. While the signals at 10.6 and 21.4 days are attributed to stellar rotation (they correspond to 0.5 and 1 time the stellar rotation period, respectively), the signal at 12.3 days is more unclear. Nevertheless, we note that its fit with a Kep-lerian yields a large eccentricity of 0.25. With the HARPS-N data alone, this remaining signal not only is less significant, but also leads to a large orbital eccentricity of 0.62 if we fit it with a Keplerian. As a result, we put strong caution about the planetary nature of this 12.3 days signal, and presently interpret it as an artefact of stellar activity. The Keplerian orbital parameters of HD 147379 b that we derived from the combined dataset suggest a minimum mass of m sin i = 25.2 M_⊕ (K = 5.25 m s⁻¹). In Fig. A.6, we present the RV folded on this Keplerian solution. Again, we stress that these results do not include a careful modelling of stellar activity. Therefore, our analysis of the HARPS-N YARARA-V1 data alone together with SPLEAF is favoured, and its results defines our final solution for the orbital parameters of HD 147379 b.

Fig. 9 HD 147379 planet search. Top and middle panels: GLS periodograms of the HARPS-N YARARA-V1 RV time-series, in a model with correlated noise and with no Keplerian in the model (k0) and one Keplerian (k1) at a period of 86.5 days. Bottom panel: RV measurements phase-folded on the 86.5 days signal.

5.2 Transit search

We scrutinised the full TESS light curve in search of any transiting signal. To proceed, we extracted the light curve of each of the 28 sectors following the same procedure as described in Sect. 2.2.1. After downloading the CBV of those sectors, we found that the combination of multi-scale and spike CBV provides the best corrections in terms of minimising the CDPP metric, with an average CDPP of -90 ppm over the sectors. We further detrended the resulting light curve from any remaining stellar systematics with keplersplinev2. Our detrending fit converged to a spline variation time-scale of 0.9 day.

We then undertook several successive transit searches, for orbital periods between 2 and 20 days, between 20 and 70 days, and between 70 and 100 days. In none of the BLS periodograms did we find significant peaks. Concerning HD 147379 b, using the mass-radius relationship of the rocky population from Otegi et al. (2020) - considering the unfavourable case of a high bulk density - and the stellar radius provided on exofop, we estimated the transit depth to -1200ppm. The planet would hence be detected given the large time-span of the TESS observations, and we conclude that it does not transit its parent star. We refer to Fig. B.3 for more details.

Fig. 10 HD 190007: periodograms of the full RV time-series. Top left: no Keplerian (k0), and no correlated noise (CN) in the model. Bottom left: one Keplerian (k1) in the model, at a period of 11.7 days. Top right: inclusion of a CN model, no Keplerian. Bottom right: inclusion of both CN and one Keplerian models.

6 HD 190007 data analysis

A planet candidate was recently suggested to orbit HD 190007 (Burt et al. 2021). The authors used the combined RV dataset of APF and HIRES in order to determine the planetary characteristics. They found a body that orbits its star with a period of 11.7 days and a minimum mass of 16.46 Μ_⊕. The stellar rotation is visible in the data, and they modelled the latter using a Keplerian with a period consistent with the photometric variability (~29 days).

6.1 RV analysis

We reviewed this solution by adding our set of 37 nightly binned HARPS-N measurements, and analysing the combined dataset. In total, we used 159 RV. Given the small number of HARPS-N RV measurements, we could not reliably process them with YARARA. Indeed, the algorithms behind the post-processing require several dozens of observations with different barycentric Earth RV (BERV) in order to not over-fit the data. Also, to work properly, having several observations per BERV bin element is mandatory, a condition not satisfied in this case. Instead, we used the data derived from the new DRS (Dumusque et al. 2021). After correcting for the different instrument offsets and accounting for drifts, the periodogram of the RV time-series clearly shows the presence of the 11.7 days planet. Additionally, a significant signal is visible with a period of ~28 days, consistent with the stellar rotation period. This periodogram is presented in the top left plot of Fig. 10. At the bottom left, we show the periodogram of the residual RV, after subtraction from a Keplerian model with a period of 11.7 days. Significant power is still observed at the rotation period and half of it. Burt et al. (2021) fitted the 28 days signal with another Keplerian.

Our three datasets - HARPS-N, HIRES and APF - contain the S-index indicator. We used this information to train a correlated noise model, and apply it simultaneously on the RV time-series. Again, we employed SPLEAF to model the correlated noise. Given the tight constraints on the stellar rotation period, we set a uniform prior on the correlated noise period between 25 and 30 days. We set priors on the other parameters of the correlated noise model according to Table 2. We initialised the model with a correlated noise period of 28 days and a decay time-scale of 200 days. Beside the correlated noise parameters, we also included in the model the instrumental offsets and a white noise term. We fitted the data with all the parameters of the model, and present the resulting periodogram of the RV residuals at the top right plot of Fig. 10. The planetary signal at 11.7 days is more significant. Its has a FAP of 3 × 10⁻¹⁰. Naturally, we further note that the signal at ~28 days disappeared. As a second step, we included a Keplerian with a period of 11.7 days in the model, and performed a new fit. The bottom right plot of Fig. 10 shows the periodogram of the new residuals. Without ambiguity, no significant signal remains. Therefore, the correlated noise modelling provides a solution cleaned from the stellar rotation effects. It is able to model the main stellar rotation period and its harmonic better than a single Keplerian as previously done by Burt et al. (2021). The minimum mass of the planet and the orbital eccentricity are estimated to 15.5 ± 1.3 Μ_⊕ and 0.14 ± 0.08, respectively. This is in agreement with Burt et al. (2021), which estimated a minimum mass of 16.5 ± 1.7 Μ_⊕ and an eccentricity of 0.14 ± 0.07. We present the updated planetary parameters in Table 3, and the phase-folded RV time-series in Fig. 11.

Fig. 11 HD 190007: RV measurements phase-folded on the Keplerian solution of HD 190007 b.

6.2 Transit search

We analysed the photometry of TESS sector 54 in search for transit signals. Again, we retrieved the 2-min cadence target pixel files and downloaded the CBV of sector 54. We applied aperture photometry and corrected the light-curve for instrumental systematics via a joint fit of multi-scale and spike CBV and background subtraction at the pixels level (as explained in Sect. 2.2.1). A special care with the detrending was needed due to the rapidly varying variability of the star. When let free, the parameter describing the time-scale of variation of the spline converged to the lower limit that we imposed (i.e. 0.5 day). We undertook several BLS transit searches between orbital periods of 0.5 and 13 days, exploring various detrend time-scales in the interval [0.5, 1.5] days with the parameter bkspace of the function keplersplinev2. Irrespective of the detrend time-scale, we did not find any significant signal in the BLS periodogram.

We computed the expected transit times, depth and duration of planet b. To proceed, we used the results of our RV analysis, the stellar radius reported on exofop, and the mass-radius relationship for rocky planets taken from Otegi et al. (2020) – again, we tested the transit detectability in the unfavourable case of a high bulk density, noting however that HD 190007 b is more likely gas-rich. We estimated the transit depth to 920 ppm, which is significant compared to the scatter observed in our detrended light curve. However, we did not find hints of transits that match with the expected ephemerides of planet b. We refer to Fig. B.4 for more details.

In a second approach, we scrutinised the TESS sector 54 considering potential single transit features. To proceed, we applied the same single-transit search process as the one employed for HD 99492 (cf. Sect. 4.6) and based on the TIP framework. Again, we ran two distinct analyses with wide priors on one side, and priors constrained on the ephemerides of HD 190007 b as derived from our RV modelling on the other side. In both cases, the absence of planetary transits was unambiguously favoured. Therefore, we rule out the presence of transits in the data.

7 Discussion and conclusions

In this work, we presented the HARPS-N spectroscopic data acquired on three stars: HD 99492, HD 147379, and HD 190007, all of them already known to harbour planets. We reviewed those planetary systems with the help of advanced data analysis tools. YARARA-V1 provided us with HARPS-N datasets cleaned from any instrumental effects. Additionally, we used the generalised SPLEAF model to carefully and efficiently mitigate stellar activity. The combination of these two tools was shown to provide the best results with respect to other strategies, as discussed for HD 99492 in Sect. 4. We updated the orbit of the known planet HD 99492 b (P = 17.05 days), and unambiguously detected a second planetary companion, namely HD 99492 c, that we confirmed independently using TWEAKS. This second planet orbits its parent star in 95.2 days, and has an estimated minimum mass of 17.9 ± 1.3 Μ_⊕. The analysis of the 165 nightly binned HARPS-N measurements of HD 147379 did not lead to a new planet detection. We also did not find any transit signal in the extensive 28 TESS sectors. However, we updated the parameters of HD 147379 b, and notably found a minimum mass smaller than previous studies (Reiners et al. 2018; Hobson et al. 2018) - our new estimate, peaking at 21.6 ± 1.1 Μ_⊕, is 2.6σ away from the nearest value provided by Reiners et al. (2018). Finally, we reviewed the system HD 190007 with the addition of our 37 nightly binned HARPS-N measurements to the publicly available HIRES and APF data. We performed a correlated noise analysis to account for the strong stellar rotation signal. We updated the solution of the known planet at P = 11.7 days, and obtained a minimum mass of 15.5 ± 1.3 Μ_⊕. The results from our review of the three planetary systems are presented in Table 3. We also carried out a systematic transit search in the available TESS sectors, and did not find transit signals.

HD 99492 c and HD 147379 b have orbital periods and minimum masses that stand in an underpopulated region of the parameter space. Figure 12a presents the distribution of known exoplanets with measured minimum masses m sin i, and with large orbital periods P > 50 days. This population is composed of 665 planets, according to the Exoplanet Archive¹² as of March 22, 2023. We note that this distribution also includes long-period transiting planets whose masses are constrained with RV measurements. The positions of HD 99492 c and HD 147379 b in this plot are represented by vertical grey and green lines, respectively. For comparison, we indicate the position of Neptune, Saturn, and Jupiter in this histogram. This distribution unambiguously reveals a larger number of long-period giant exoplanets. Naturally, observational biases play an important role in this picture, since Neptune-mass exoplanets on long orbital periods are technically more challenging to detect. However, a closer inspection of this histogram may reveal a potential lack of planets between 30 and 50 Earth masses, giving rise to two distinct populations that the observational biases hardly explain.

To estimate the statistical significance of this bimodality, and hence the existence of two planet populations, we tested the similarity between this distribution and a single-mode Gaussian distribution. To proceed, we defined the large-mass sub-sample by setting a mass threshold at 50 Μ_⊕, and fitted a Gaussian law on this sub-sample. Then, we generated a random sample following this law, and with the same size as the full distribution. Finally, we computed the p-value between this Gaussian and the blue histogram. We repeated this procedure 10 000 times, and derived a distribution of p-values. In those results, the maximum p-value obtained is 3 × 10⁻⁴, which indicates that all of our 10 000 tests support the hypothesis of two distinct distributions between the planet distribution and the fitted Gaussian. As a second statistical test, we performed a Z-test between the sample of planets with m sin i < 50 M_⊕ (red shaded area), and the sample with m sin i > 50 M_⊕ (grey shaded area) to quantify their mutual difference. This test compares the means and standard deviations of two samples, and indicates if these samples significantly differ from each other. We fitted both large-mass and small-mass sub-samples with a Gaussian, and applied the Z-test on these adjusted laws. The test returned the value of 2.94, which points towards a significant difference between the two populations. We repeated the Z-test with a new separation threshold of 30 M_⊕ between the two populations, and obtained the value of 3.06, again indicating a significant difference between these two populations. Therefore, our analyses support the bimodality of the planet distribution.

We bring an additional support to the bimodality in the bottom panel of Fig. 12a. It presents, among the considered population, the percentage of planets found in multi-planet systems. The planets seem to be globally found in distinct system architectures. Despite the large uncertainties in the small-mass population caused by the small number of detections, we observe that these planets are most often found in multi-planet systems. For higher-mass planets, this ratio is significantly smaller. Those results support previous RV-detected planet statistics. Notably, Winn & Fabrycky (2015) pinpoint a bimodality in the distribution of minimum masses, without distinction of orbital period (see their Fig. 4). In that same figure, the authors illustrate, at the population level, the increase in eccentricity scatter for larger minimum masses. Besides, among the exoplanet population, lower multiplicity systems present also a larger eccentricity scatter (Wright et al. 2009; Limbach & Turner 2015). Through the link of orbital eccentricity, we hence deduce that large-mass planets tend to be found in systems with smaller multiplicity, according to the RV surveys. The results that we present in the lower panel of Fig. 12a corroborate these observations, since larger-mass planets seem to be more often found in single-planet systems. Again, while we raise a hint of a trend, we stress that a treatment of detection biases is needed to confirm this statement.

We compared these results with the short-period planet population. Figure 12b shows the population of exoplanets with measured minimum masses, but this time with a focus on short-period planets (P < 50 days). The double peak distribution is more prominent than for the long-period population, and a similar trend with the system multiplicity is observed. We undertook a new Z-test for a threshold mass of 50 M_⊕, and obtained a value of 2.78 which confirms the likely bimodality of the distribution. We note that at short orbital periods, low-mass planets appear more common than giant planets. The opposite is observed for long-period planets. This observation is in line with classical theories of planet formation, for which gas giants are formed in cooler regions of the protoplanetary disk, as opposed to smaller-mass planets. However, we stress that detection biases likely alter this picture. Furthermore, we highlight with the red histogram the short-period planet companions to the low-mass long-period planets located in the red zone of Fig. 12a. Most of these planets have smaller masses too. Conversely, we highlight with the grey histogram the short-period planet companions to the large-mass long-period planets inside the grey zone of Fig. 12a, which does not indicate a strong preference for large or small mass companions. Additionally, the red histogram constitutes most of the planet companions to low-mass long-period planets with a measured minimum mass, with 80 planets out of 94 companions in total, that is 85%. The remaining 15% are long-period planets. As a result, most of the companions to the small-mass long-period planets have small masses and short periods.

Based on all those observations, we would therefore expect HD 147379 to be orbited by at least a second planetary companion inner to HD 147379 b, with the hypothesis that this system follows the observed trend. We stress, however, that the statistics have low numbers and that the observational biases could alter this picture. Continuing observational efforts using high-resolution spectrographs will help to augment those statistics and refine the architectures of planetary systems. Finally, we note that a similar significant bimodality is observed in the true mass distribution of exoplanets, as corroborated by Bertaux & Ivanova (2022) using a numerical inversion of the m sin i distribution. This apparent desert in the exoplanet population is in line with population synthesis models (Ida & Lin 2004; Mordasini et al. 2009), and strengthen former observational hints (e.g. Arriagada et al. 2013; Winn & Fabrycky 2015).

We put these results in light of the Kepler and TESS transit surveys in Fig. 13. The histogram presents the population of confirmed Kepler and TESS planets spread over radius, in the range [3.5, 15] R_⊕ so to entail the Neptune-size regime and above. In these data, there might also be a tentative hint of two distinct populations separated around 7 R_⊕, which corresponds to a mass of 38 Μ_⊕ using the mass-radius relationship for volatile-rich planets (Otegi et al. 2020). The location of this boundary hence overlaps with the ones observed in Figs. 12a,b. Nevertheless, this separation is here too shallow for any further prospects. Moreover, there is no clear trend in the proportion of planets in multi-planet systems as a function of planet radius in this considered radius range, as is illustrated in the bottom panel. Different factors could explain this discrepancy with the RV surveys, such as a non homogeneous scaling law for the planets bulk density, or probing different types of systems from the transit and RV surveys. This is briefly discussed in Leleu et al. (2023), in the context of the mass dichotomy between RV and transit-timing variations (TTV) techniques.

This work has once again proven the importance of long-term RV surveys with stable high-precision spectrographs for the detection and characterisation of warm planets, which constitute a yet under-explored family of exoplanets. Additionally, the use of advanced tools for the removal of instrumental systematics and for modelling correlated noise sheds light on the noise pattern, and pushes further the barriers of small-mass planet detections.

Fig. 12 Exoplanet population with a minimum mass measurement m sin i. Panel a: sub-population with large orbital periods (P > 50 days). This sample contains 665 exoplanets. Top: histogram of that population, with respect to the planet minimum mass. The vertical grey and green lines identify the position of HD 99492 c and HD 147379 b, respectively. Bottom: proportion of planets that are part of multi-planet systems, among the considered sub-population. Panel b: Sub-population with small orbital periods (P < 50 days, 389 exoplanets). Top: distribution of that subpopulation according to m sin i (in orange). The red histogram indicates companions to the low-mass long-period planets (red zone in panel a). The grey histogram focuses on companions to the large-mass long-period planets (grey zone in panel a). Bottom: Proportion of planets part of multi-planet systems.

Fig. 13 Exoplanet population from the Kepler and TESS surveys, distributed in planet radius between 3.5 and 15 R⊕. This population is composed of 632 planets as of March 22, 2023. Top: histogram of the considered populations. Bottom: proportion of planets part of multiplanet systems.

Acknowledgements

This work is based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated by the Fundacion Galileo Galilei (FGG) of the Istituto Nazionale di Astrofisica (INAF) at the Observatorio del Roque de los Muchachos (La Palma, Canary Islands, Spain). The HARPS-N project was funded by the Prodex Program of the Swiss Space Office (SSO), the Harvard University Origin of Life Initiative (HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of Geneva, the Smithsonian Astro-physical Observatory (SAO), the Italian National Astrophysical Institute (INAF), University of St. Andrews, Queen’s University Belfast, and University of Edinburgh. This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. The authors acknowledge the financial support of the SNSF. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. M.C. acknowledges the SNSF support under grant P500PT_211024. R.D.H. is funded by the UK Science and Technology Facilities Council (STFC)’s Ernest Rutherford Fellowship (grant number ST/V004735/1). F.P.E. and C.L.O. would like to acknowledge the Swiss National Science Foundation (SNSF) for supporting research with HARPS-N through the SNSF grants nr. 140649, 152721, 166227 and 184618. The HARPS-N Instrument Project was partially funded through the Swiss ESA-PRODEX Programme. T.G.W. acknowledges support from STFC consolidated grant number ST/V000861/1 and UKSA grant number ST/R003203/1. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement SCORE No 851555). A.S. acknowledges financial support from the agreement ASI-INAF n.2018-16-HH.0. Finally, the authors thank the referee for their insightful comments and suggestions on the paper.

Appendix A Complementary spectroscopic data

We present in this section complementary spectroscopic data, in the order with which they are mentioned in the main text. Fig. A.1 illustrates the entire high-precision RV time-series of HD 147379, built up from the combination of HIRES, CARMENES, SOPHIE and HARPS-N data. Fig. A.2, A.3, and A.4 display the various spectroscopic activity indicator data on HD 99492, HD 147379 and HD 190007, respectively. Trends and correlations with the RV are hinted. Fig. A.5 presents the TWEAKS results on HD 99492, which confirm the detection of two planet signals and help to rule out the planetary nature of the 14d signal. Finally, Fig. A.6 depicts the HD 147379 RV folded on our Keplerian solution for planet b obtained with all the available RV data. This solution results from an analysis of the RV only, and hence it is not favoured. We refer to Table 3 for the preferred updated Keplerian parameters of HD 147379 b, based on the analysis of the HARPS-N YARARA-V1 data together with SPLEAF.

Fig. A.1 Entire RV time-series of HD 147374. It combines our HARPS-N RV measurements with the public data from HIRES, CARMENES, and SOPHIE high-resolution spectrographs.

Fig. A.2 HARPS-N RV and spectroscopic activity indicators of HD 99492. This dataset is cleaned from instrumental systematics but still contains the stellar activity features, as obtained with the YARARA-V1 post-process described in Sect. 3.1 and Appendix C. The left column presents the RV and indicator time-series, while the middle column contains the periodograms of these time-series. No drift nor white noise was included prior to their computation. However concerning the activity indicators, we included a Keplerian in the model so to absorb most of the long-period power, which impacts the short-period signals. In those cases, we present the periodograms of the residuals. The red bands depict the estimated range for the stellar rotation period, obtained from both literature and a preliminary analysis of our data. The right column illustrates the correlation between the various indicators and the RV.

Fig. A.3 HARPS-N YARARA-V1 RV and spectroscopic activity indicators of HD 147379. Similar description to Fig. A.2, except with the periodogram computation. The long-period trends observed in the time-series of the spectroscopic indicators are better fitted by quadratic drifts instead of Keplerians. The corresponding periodograms were computed from the residuals of these fits.

Fig. A.4 HARPS-N YARARA-Vl, HIRES and APF RV and spectroscopic activity indicators of HD 190007. Similar description to Fig. A.2. We did not include any Keplerian nor drift prior to the computation of the periodograms. We note that only two and one chromospheric activity indicators are publicly available for the APF and HIRES data, respectively (namely the Hα and S-index indicators). Hα being absent from the HIRES data (Butler et al. 2017), the correlation coefficient between the RV and Hα only takes the HARPS-N and APF data into account. Furthermore, this correlation is stronger in HARPS-N, for which R = 0.47.

Fig. A.5 HD99492: Results from the TWEAKS pipeline. Left column: posterior distribution. Right column: FIP periodogram. From top to bottom rows: Analyses on the entire RV time-series, on the first half and on the second half of the time-series, respectively. The vertical plain lines (light green) depict the two planet detections, at 17.05 and 95.2 days, respectively. The red dashed lines help to guide the eye on the 14d signal, which is inconsistently present in the data.

Fig. A.6 Entire RV time-series of HD 147374 folded on the best-fit Keplerian solution for the 86.5d planet. This time-series combines our HARPS-N RV measurements with the public data from HIRES, CARMENES, and SOPHIE high-resolution spectrographs. The best-fit Keplerian model is illustrated with the grey curve.

Appendix B Transit searches

This section provides details on our transit searches in the three stars. Fig. B.1, B.3, and B.4 present the corrected and detrended light-curves together with the hypothetical transit models for the four planets around HD 99492, HD 147379, and HD 190007, respectively. To obtain a light-curve cleaned from systematic trends and amenable to a reliable transit search, we employed a two-stages process. First, we used lightkurve to treat the instrumental systematics mostly (see Sect. 2.2.1 for more details), and then, we detrended the light-curve from the remaining stellar activity trends with keplersplinev2 (see Sect. 4.6 for more details). These cleaned light-curves are presented in the top panel of each of the three figures. To illustrate the transit detectability of the planets, we generated hypothetical transit model curves with the assumptions of null impact parameter and null eccentricity. Furthermore, these curves are drawn from the unfavourable assumption of rocky composition, even though we expect those planets to stand in the volatile-rich regime. This decision aims to strengthen the transit detectability, by stating that the planets would be detected even under unexpectedly high bulk densities. The transit depths of the model curves were derived from the mass-radius relationship for rocky planets provided by Otegi et al. (2020), utilising the results from our RV modelling, together with the stellar radius indicated on exofop. In the bottom panels of Fig. B.1, B.3, and B.4, we present the cleaned light-curves phase-folded onto the ephemerides of the considered planets, together with the transit model curves. When relevant, we also provide zoomed-in views of these plots, centred around the expected transit time. The red dots are binned flux measurements, and allow for a better representation of the variations in the light curves.

Finally, Fig. B.2 illustrates our single-transit search in the TESS sectors 45 and 46 of HD 99492, using the TIP framework. The search was performed in consecutive 1-day temporal windows, and the plot presents one of the outcomes. The singletransit search process is explained in Sect. 4.6. The TIP compares two models, without and with one planet. In none of the one-day windows did we find a TIP higher than 0.5, which means that the zero-planet model is favoured. A similar procedure was undertaken with the TESS sector 54 of HD 190007, and the same conclusion was reached.

Fig. B.1 HD 99492: Transit search in the TESS light-curve, sectors 45 and 46. Top panel: cleaned light-curve following our custom extraction from the target pixel files. Second and third panels: Light-curve folded on the orbital period and phase of planet b predicted from our RV analysis, together with the transit model curve (the modelled transit depth is 1038ppm for a rocky composition). Bottom panel: Light-curve folded on the ephemerides of planet c, together with its theoretical transit curve for a rocky composition (the modelled transit depth is 850ppm).

Fig. B.2 Illustration of the single-transit search in the HD 99492 TESS sectors 45 and 46 data.

Fig. B.3 HD 147379: Transit search in the TESS light-curve composed of 28 sectors (sectors 14-21, 23-26, 40-41 and 47-60). Top panel: cleaned light-curve. Second and third panels: Light-curve folded on the orbital period and phase of planet b predicted from our RV analysis, together with the transit model curve (the modelled transit depth for a rocky composition is 1250ppm).

Fig. B.4 HD 190007: Transit search in the TESS light-curve, sector 54. Top panel: cleaned light-curve. Second and third panels: Light-curve folded on the orbital period and phase of planet b predicted from our RV analysis, together with the transit model curve (the modelled transit depth and duration for a rocky composition are 920ppm and 190 minutes, respectively).

Appendix C Orthogonalisation of stellar activity and instrumental systematic vectors in YARARA

During the cleaning process of spectra, YARARA is fitting simultaneously for stellar activity and instrumental systematics in the wavelength domain by fitting a vector basis in a multilinear model. In the previous version of the code (Cretignier et al. 2021), the basis was fitted using the S-index, the CCF’s contrast and the CCF’s FWHM. However, on HARPS-N, the CCF moments usually contain the superposition of 1) a PSF variation (induced by an instrumental defocus or technical operations on the instrument) and 2) a stellar activity component that affects the individual stellar lines (Cretignier et al. 2020a). Because of that, it was not possible to re-inject the stellar activity component fitted out by YARARA, without reintroducing as well the PSF defocus systematics. In order to remove part of the degeneracy, a better splitting of activity and instrumental systematics was mandatory, which was obtained by producing a better vector basis.

Since the S-index is a broad flux integration in the core of the strong CaIIH&K lines, the proxy is insensitive to small PSF defocus, opposite to the thin photospheric stellar lines that dominate the moments of the CCF. On the other hand, PSF variations are mainly symmetrical. This led us to develop a new methodology to derive better time-domain vectors orthogonal to the stellar activity S-index component. We show the results of this analysis on HD4628 as well as on HD99492. HD4628 is not presented in the present article and its analysis will be described in a forthcoming paper. However, the star is one of the most consistently monitored by HARPS-N from the very beginning. The star also possesses a clear magnetic cycle and has a similar spectral type than HD99492. This is why it was selected for the present demonstration.

Using the same generic line lists as the DRS and for a given velocity grid v_i, we first derived the cross-correlation functions for all the epochs: CCF(t, The CCFs were centred by an RV amount corresponding to the centre of a Gaussian fit before being interpolated on the initial velocity grid. Such an operation allows to cancel any large Doppler shift. The median CCF was then subtracted to produce residual profiles and each of these profiles was itself median centred afterwards. The resulting CCF residuals are plotted in the top panel of Fig. C.1.

In order to remove the stellar activity dependency, the residual CCFs for each velocity grid element v_i were linearly decorrelated from the S_index(t) time-series: CCF_res(t, v_i) = CCF(t, v_i) – α_i. S_index(t). In Cretignier et al. (2021), the authors showed that the S-index was strongly correlated with the line profiles distortions (see their Fig.3), hence the same correlations are expected for the CCFs. However, the S-index usually also contains long-term variations due to stellar magnetic cycles (Baliunas et al. 1995). Such signals can be easily collinear with the long-term drift defocus observed in the HARPS and HARPS-N instruments, producing incorrect scaling coefficients α_i. In order to avoid this issue, the scaling coefficients between the S-index and the residual CCF time-series were obtained after performing a low-pass filtering. The fitted scaling coefficients between the S-index and CCF residuals were obtained on the high frequency part of the filtered vectors, which implicitly means that the correlation was done on the rotational modulation of the S-index signal and not on the magnetic cycle. The CCF_res profiles obtained can be seen in the bottom panel of the Fig. C.1.

The decorrelated residual profiles were split into a left (L(t,v_i) = CCF_res(t,v_i),v_i < 0) and a right (R(i,v_i) = CCF_res(t,v_i),v_i > 0) wings profile. Symmetric S (t, v_i) and anti-symmetric A(t, v_i) vectors were created, respectively computing S (t, v_i) = 0.5 · (L(t, −|v_i|) + R(t, |v_i|)) and A(t,v_i) = 0.5 · (L(t, −|v_i|) – R(t, |v_i|)). Those vectors help categorising symmetric versus anti-symmetric CCF variations from the PCA decomposition performed here after.

Fig. C.1 Visualisation of the CCF residuals profiles obtained before (top panel) and after (bottom panel) the decorrelation by the S-index for HD4628. The colour code is the Z-normalised S-index vector (in arbitrary units). The stellar activity signature, initially visible in the top panel by the blue-to-green colour gradient at the line centre (v = 0 km s−1), is removed, whereas the clear instrumental defocus is preserved. The instrumental time-domain vectors can be extracted after performing a PCA decomposition of those profiles (see Fig. C.2).

Fig. C.2 Extraction of the instrumental systematics basis for HD4628 by performing a PCA decomposition on the M(t, vi) matrix (as defined in the text). The first five principal components Pj(vi) are displayed (left panel) with an arbitrary offset as well as their time-domain coefficients Uj(t) and their corresponding periodogram (right panel). The grey periodograms were derived after applying a high-pass filter on the time-series, while the black ones do not correspond to any filtering. The S-index (top row) is also displayed for comparison. The ratio score R is displayed for each component. Only the first two components (blue and orange) are found symmetric (R1 = 91.0 and R2 = 15.1). The first one clearly shows the instrumental defocus during the first part of the instrument’s life (1st and 2nd seasons), whereas an instrumental refocus is also visible at the middle of the penultimate season, which corresponds to a manual intervention on the instrument on the 4th of October 2021 to fix the HARPS-N cryostat (priv. communication).

Indeed, because the PSF variation is unknown a priori, a PCA was performed simultaneously on the symmetric and anti-symmetric profiles of the joint matrix M(t, v_i) = [S (t, v_i); A(t, v_i)]. The principal components P_j(v_i) and time-coefficients U_j (t) were recovered. The first five components are represented in Fig. C.2 as well as the S-index coefficient on top (also converted into a symmetric and anti-symmetric component for comparison purpose). Each PCA component was then classified computing a ratio score R defined as: $$\{\begin{array}{l}R={\sigma }_{sym}/{\sigma }_{anti},\text{if}{\sigma }_{sym}>{\sigma }_{anti}\hfill \\ R=-{\sigma }_{anti}/{\sigma }_{sym},\text{if}{\sigma }_{sym}<{\sigma }_{anti}\hfill \end{array}$$

with σ_sym the standard deviation of P_j(v_i) on the symmetric part and σ_anti the standard deviation in the anti-symmetric side. We considered a component as symmetric if R > 2, anti-symmetric if R < −2 and noisy otherwise. We then selected the N first symmetric components until a first asymmetric or noisy component was detected.

Similar components were identified for both HD4628 and HD99492 (see Fig. C.3), even if for the latter the components are noisier as expected from the lower S/N of the observations. It can be noted how similar the activity profile and the instrumental defocus profile are, showing how challenging the decoupling of the signals would be without the S-index decorrelation. The third PCA component is anti-symmetric, even though its score is sometimes just below the threshold value (e.g R = −1.6 for HD99492). This component shows some magnetic cycle and rotational period (marked by red bands on the periodograms in Fig. C.3) modulations and correlates slightly with the VSPAN of the CCFs, which is a CCF’s moment measuring the asymmetry of the CCF profile.

The present algorithm was successfully tested on dozens of stars and showed to be efficient to produce time-domain vectors dominated by the instrumental defocus observed on HARPS-N and HARPS. Among those stars: the Sun (Dumusque et al. 2021), for which results will also be presented in a forthcoming paper. Although the solar data is not in the scope of the present paper, it is tempting to use them in order to interpret the PCA components, even if it is always a dangerous and difficult task since components are often a mixture of contaminations that change from star-to-star. The first component on the Sun contains almost exclusively the change of the v sin i (Collier Cameron et al. 2019) that may be interpreted as a change of the PSF in some way. The second solar component seems strongly correlated with the S/N of the solar observations which are related to the ageing of the transparent dome in which the solar telescope is located (Dumusque et al. 2015). Contrarily to the Sun observations which are obtained with the solar telescope, observations of other stars are obtained with the TNG which is not covered by such a dome. While such systematics are necessarily different in origin for other stellar observations, it may still indicate that this component is related to unusual S/N perhaps due to clouds or peculiar observing conditions (seeing).

Fig. C.3 Same as Fig. C.2 but applied on the HD99492 dataset.

Note that, by design, nothing prevents stellar activity to be slightly present in the derived instrumental vectors since residual symmetric activity variation uncorrelated with the S-index may be captured by the PCA. However, even if such a risk exists, such a component would be now strongly mitigated, as is depicted by the absence of rotational peak in the periodogram of the extracted components.

References

All Tables

All Figures

	Fig. 1 RV time-series of the three stars analysed in this study. Top: HD 99492 (HARPS-N and HIRES data). Middle: HD 147379 (HARPS-N, SOPHIE, CARMENES and HIRES data). Bottom: HD 190007 (HARPS-N, HIRES and APF data).
In the text

	Fig. 2 HD 99492: periodic signal search in TESS sectors 45 and 46. Top: TESS light-curve extracted using lightkurve, and corrected from systematics (2 sectors). Bottom: GLS periodogram computed between 2 and 40 days of the 1-day binned light-curve.
In the text

	Fig. 3 HD 147379 TESS photometry. Top: TESS light-curve extracted using lightkurve, and corrected from systematics (28 sectors). Bottom: GLS periodogram computed between 2 and 40 days of the 1-day binned light curve.
In the text

	Fig. 4 HD 190007 TESS SAP normalised time-series (sector 54).
In the text

	Fig. 5 HD 99492: YARARA-V1 S-index dataset. Top: S-index time-series. Middle: GLS periodogram of that time-series (k0 stands for no Keplerian in the model). Bottom: GLS periodogram of the residual S-index, after subtracting a Keplerian model with P ~ 3000 days (k1). The red bands locate the periodicities seen in this activity indicator.
In the text

	Fig. 6 HD 99492: periodograms of the combined YARARA-V1 and HIRES RV time-series. The periodograms are computed from a fit of the time-series with a detrend + drift model. Top: no Keplerian. Middle: periodogram of the residual time-series, after removing a Keplerian at ~17 days. Bottom: periodogram of the residuals after removing Keplerians at ~17 and 95 days.
In the text

Fig. 7 HD 99492 c posterior distributions for our different models and datasets. The grey, blue, orange, and green distributions refer to our first, second, third, and fourth approaches to derive the posteriors, respectively. In the case of the third approach, we further show the distributions resulting from two different models. The plain line distribution was obtained from a model with two Keplerians, while the dashed line corresponds to a model with three Keplerians.

In the text

Fig. 8 HD 99492 orbital solutions for planets b and c derived from the HARPS-N YARARA-V1 dataset. Left: Cartesian representation of the planetary system for a sample of 205 MCMC solutions. The star is in red, the orbit of HD 99492 b in blue and the orbit of HD 99492 c in orange. Additionally, the green area marks the beginning of the conservative habitable zone, which continues outwards. This was calculated from the prescription of Kopparapu et al. (2013). Right: RV measurements phase-folded on the best fit Keplerian solutions, together with the model curves.

In the text

	Fig. 9 HD 147379 planet search. Top and middle panels: GLS periodograms of the HARPS-N YARARA-V1 RV time-series, in a model with correlated noise and with no Keplerian in the model (k0) and one Keplerian (k1) at a period of 86.5 days. Bottom panel: RV measurements phase-folded on the 86.5 days signal.
In the text

	Fig. 10 HD 190007: periodograms of the full RV time-series. Top left: no Keplerian (k0), and no correlated noise (CN) in the model. Bottom left: one Keplerian (k1) in the model, at a period of 11.7 days. Top right: inclusion of a CN model, no Keplerian. Bottom right: inclusion of both CN and one Keplerian models.
In the text

	Fig. 11 HD 190007: RV measurements phase-folded on the Keplerian solution of HD 190007 b.
In the text

Fig. 12 Exoplanet population with a minimum mass measurement m sin i. Panel a: sub-population with large orbital periods (P > 50 days). This sample contains 665 exoplanets. Top: histogram of that population, with respect to the planet minimum mass. The vertical grey and green lines identify the position of HD 99492 c and HD 147379 b, respectively. Bottom: proportion of planets that are part of multi-planet systems, among the considered sub-population. Panel b: Sub-population with small orbital periods (P < 50 days, 389 exoplanets). Top: distribution of that subpopulation according to m sin i (in orange). The red histogram indicates companions to the low-mass long-period planets (red zone in panel a). The grey histogram focuses on companions to the large-mass long-period planets (grey zone in panel a). Bottom: Proportion of planets part of multi-planet systems.

In the text

	Fig. 13 Exoplanet population from the Kepler and TESS surveys, distributed in planet radius between 3.5 and 15 R⊕. This population is composed of 632 planets as of March 22, 2023. Top: histogram of the considered populations. Bottom: proportion of planets part of multiplanet systems.
In the text

	Fig. A.1 Entire RV time-series of HD 147374. It combines our HARPS-N RV measurements with the public data from HIRES, CARMENES, and SOPHIE high-resolution spectrographs.
In the text

Fig. A.2 HARPS-N RV and spectroscopic activity indicators of HD 99492. This dataset is cleaned from instrumental systematics but still contains the stellar activity features, as obtained with the YARARA-V1 post-process described in Sect. 3.1 and Appendix C. The left column presents the RV and indicator time-series, while the middle column contains the periodograms of these time-series. No drift nor white noise was included prior to their computation. However concerning the activity indicators, we included a Keplerian in the model so to absorb most of the long-period power, which impacts the short-period signals. In those cases, we present the periodograms of the residuals. The red bands depict the estimated range for the stellar rotation period, obtained from both literature and a preliminary analysis of our data. The right column illustrates the correlation between the various indicators and the RV.

In the text

	Fig. A.3 HARPS-N YARARA-V1 RV and spectroscopic activity indicators of HD 147379. Similar description to Fig. A.2, except with the periodogram computation. The long-period trends observed in the time-series of the spectroscopic indicators are better fitted by quadratic drifts instead of Keplerians. The corresponding periodograms were computed from the residuals of these fits.
In the text

Fig. A.4 HARPS-N YARARA-Vl, HIRES and APF RV and spectroscopic activity indicators of HD 190007. Similar description to Fig. A.2. We did not include any Keplerian nor drift prior to the computation of the periodograms. We note that only two and one chromospheric activity indicators are publicly available for the APF and HIRES data, respectively (namely the Hα and S-index indicators). Hα being absent from the HIRES data (Butler et al. 2017), the correlation coefficient between the RV and Hα only takes the HARPS-N and APF data into account. Furthermore, this correlation is stronger in HARPS-N, for which R = 0.47.

In the text

Fig. A.5 HD99492: Results from the TWEAKS pipeline. Left column: posterior distribution. Right column: FIP periodogram. From top to bottom rows: Analyses on the entire RV time-series, on the first half and on the second half of the time-series, respectively. The vertical plain lines (light green) depict the two planet detections, at 17.05 and 95.2 days, respectively. The red dashed lines help to guide the eye on the 14d signal, which is inconsistently present in the data.

In the text

	Fig. A.6 Entire RV time-series of HD 147374 folded on the best-fit Keplerian solution for the 86.5d planet. This time-series combines our HARPS-N RV measurements with the public data from HIRES, CARMENES, and SOPHIE high-resolution spectrographs. The best-fit Keplerian model is illustrated with the grey curve.
In the text

Fig. B.1 HD 99492: Transit search in the TESS light-curve, sectors 45 and 46. Top panel: cleaned light-curve following our custom extraction from the target pixel files. Second and third panels: Light-curve folded on the orbital period and phase of planet b predicted from our RV analysis, together with the transit model curve (the modelled transit depth is 1038ppm for a rocky composition). Bottom panel: Light-curve folded on the ephemerides of planet c, together with its theoretical transit curve for a rocky composition (the modelled transit depth is 850ppm).

In the text

	Fig. B.2 Illustration of the single-transit search in the HD 99492 TESS sectors 45 and 46 data.
In the text

	Fig. B.3 HD 147379: Transit search in the TESS light-curve composed of 28 sectors (sectors 14-21, 23-26, 40-41 and 47-60). Top panel: cleaned light-curve. Second and third panels: Light-curve folded on the orbital period and phase of planet b predicted from our RV analysis, together with the transit model curve (the modelled transit depth for a rocky composition is 1250ppm).
In the text

	Fig. B.4 HD 190007: Transit search in the TESS light-curve, sector 54. Top panel: cleaned light-curve. Second and third panels: Light-curve folded on the orbital period and phase of planet b predicted from our RV analysis, together with the transit model curve (the modelled transit depth and duration for a rocky composition are 920ppm and 190 minutes, respectively).
In the text

Fig. C.1 Visualisation of the CCF residuals profiles obtained before (top panel) and after (bottom panel) the decorrelation by the S-index for HD4628. The colour code is the Z-normalised S-index vector (in arbitrary units). The stellar activity signature, initially visible in the top panel by the blue-to-green colour gradient at the line centre (v = 0 km s−1), is removed, whereas the clear instrumental defocus is preserved. The instrumental time-domain vectors can be extracted after performing a PCA decomposition of those profiles (see Fig. C.2).

In the text

Fig. C.2 Extraction of the instrumental systematics basis for HD4628 by performing a PCA decomposition on the M(t, vi) matrix (as defined in the text). The first five principal components Pj(vi) are displayed (left panel) with an arbitrary offset as well as their time-domain coefficients Uj(t) and their corresponding periodogram (right panel). The grey periodograms were derived after applying a high-pass filter on the time-series, while the black ones do not correspond to any filtering. The S-index (top row) is also displayed for comparison. The ratio score R is displayed for each component. Only the first two components (blue and orange) are found symmetric (R1 = 91.0 and R2 = 15.1). The first one clearly shows the instrumental defocus during the first part of the instrument’s life (1st and 2nd seasons), whereas an instrumental refocus is also visible at the middle of the penultimate season, which corresponds to a manual intervention on the instrument on the 4th of October 2021 to fix the HARPS-N cryostat (priv. communication).

In the text

	Fig. C.3 Same as Fig. C.2 but applied on the HD99492 dataset.
In the text