# Revisiting the multi-planetary system of the nearby star HD 20794

N. Nari, X. Dumusque, N. C. Hara, A. Suárez Mascareño, M. Cretignier, J. I. González Hernández, A. K. Stefanov, V. M. Passegger, R. Rebolo, F. Pepe, et al.

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

► To cite this version:

N. Nari, X. Dumusque, N. C. Hara, A. Suárez Mascareño, M. Cretignier, et al. Revisiting the multi-planetary system of the nearby star HD 20794. Astronomy & Astrophysics - A&A, 2025, 693, pp.A297. 10.1051/0004-6361/202451769. hal-05486794

HAL Id: hal-05486794

https://hal.science/hal-05486794v1

Submitted on 23 Feb 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Revisiting the multi-planetary system of the nearby star HD 20794

## Confirmation of a low-mass planet in the habitable zone of a nearby

## G-dwarf ⭐, ⭐

N. Nari<sup>1,2,3,★★</sup>, X. Dumusque<sup>5</sup>, N. C. Hara<sup>5,6</sup>, A. Suárez Mascareño<sup>2,3</sup>, M. Cretignier<sup>7</sup>, J. I. González Hernández<sup>2,3</sup>, A. K. Stefanov<sup>2,3</sup>, V. M. Passegger<sup>22,2,3,23</sup>, R. Rebolo<sup>2,3,4</sup>, F. Pepe<sup>5</sup>, N. C. Santos<sup>8,9</sup>, S. Cristiani<sup>12,13</sup>, J. P. Faria<sup>5,8,9</sup>, P. Figueira<sup>5,8</sup>, A. Sozzetti<sup>14</sup>, M. R. Zapatero Osorio<sup>11</sup>, V. Adibekyan<sup>8,9</sup>, Y. Alibert<sup>18,19</sup>, C. Allende Prieto<sup>2,3</sup>, F. Bouchy<sup>5</sup>, S. Benatti<sup>16</sup>, A. Castro-González<sup>11</sup>, V. D'Odorico<sup>12,13</sup>, M. Damasso<sup>14</sup>, J. B. Delisle<sup>5</sup>, P. Di Marcantonio<sup>12</sup>, D. Ehrenreich<sup>5</sup>, R. Génova-Santos<sup>2,3</sup>, M. J. Hobson<sup>5</sup>, B. Lavie<sup>5</sup>, J. Lillo-Box<sup>11</sup>, G. Lo Curto<sup>15</sup>, C. Lovis<sup>5</sup>, C. J. A. P. Martins<sup>8,10</sup>, A. Mehner<sup>15</sup>, G. Micela<sup>16</sup>, P. Molaro<sup>12</sup>, C. Mordasini<sup>19</sup>, N. Nunes<sup>17</sup>, E. Palle<sup>2,3</sup>, S.P. Quanz<sup>20,21</sup>, D. Ségransan<sup>5</sup>, A. M. Silva<sup>8,9</sup>, S. G. Sousa<sup>8,9</sup>, S. Udry<sup>5</sup>, N. Unger<sup>5</sup>, and J. Venturini<sup>5</sup>

*(Affiliations can be found after the references)*

Received 2 August 2024 / Accepted 15 October 2024

## ABSTRACT

**Context.** Close-by Earth analogs and super-Earths are of primary importance because they will be preferential targets for the next generation of direct imaging instruments. Bright and close-by G-to-M type stars are preferential targets in radial velocity surveys to find Earth analogs. Their brightness allows us to achieve the best precision on RV measurements and search for signals with amplitudes of less than 1 m s⁻<fnref n="1" />.

**Aims.** We present an analysis of the RV data of the star HD 20794, a target whose planetary system has been extensively debated in the literature. The broad time span of the observations makes it possible to find planets with signal semi-amplitudes below 1 m s⁻<fnref n="1" /> in the habitable zone.

**Methods.** We analyzed RV datasets spanning more than 20 years. We monitored the system with ESPRESSO. We joined ESPRESSO data with the HARPS data, including archival data and new measurements from a recent program. We applied the post-processing pipeline YARARA to HARPS data to correct systematics, improve the quality of RV measurements, and mitigate the impact of stellar activity.

**Results.** We confirm the presence of three planets, with periods of 18.3142 ± 0.0022 d, 89.68 ± 0.10 d, and 647.6<sup>+2.5</sup><sub>-2.7</sub> d, along with masses of 2.15 ± 0.17 M⊕, 2.98 ± 0.29 M⊕, and 5.82 ± 0.57 M⊕ respectively. For the outer planet, we find an eccentricity of 0.45<sup>+0.10</sup><sub>-0.11</sub>, whereas the inner planets are compatible with circular orbits. The latter is likely to be a rocky planet in the habitable zone of HD 20794. From the analysis of activity indicators, we find evidence of a magnetic cycle with a period of ~3000 d, along with evidence pointing to a rotation period of ~39 d.

**Conclusions.** We have determined the presence of a system of three planets orbiting the solar-type star HD 20794. This star is bright (V=4.34 mag) and close (d = 6.04 pc), and HD 20794 d resides in the stellar habitable zone, making this system a high-priority target for future atmospheric characterization with direct imaging facilities.

**Key words.** techniques: radial velocities – techniques: spectroscopic – planets and satellites: detection – planets and satellites: terrestrial planets – stars: activity – stars: individual: HD 20794

## 1. Introduction

Exoplanetary science has experienced a rapid and successful development since the discovery of the first planet orbiting

around a main sequence star other than the Sun (Mayor et al. 1995). This detection was obtained through the radial velocity (RV) technique, which is still one of the most successful techniques for discovering and characterizing exoplanets. This seminal discovery has opened the field to more than 5600 confirmed planets up to July 2024, according to the NASA Exoplanet Archive<fnref n="1" /> (Akeson et al. 2013).

At the onset of the exoplanet field, the main issue for the detection of new planets was related to the precision of available instruments. The High Accuracy Radial velocity Planet Searcher (HARPS) instrument (Mayor et al. 2003) represents a groundbreaking advancement in the field. It is installed at the 3.6 m Telescope in La Silla, Chile and it has been the first instrument

<footnote>* Program IDs for the ESPRESSO observations we used in the analysis are: 106.21M2.003, 106.21M2.004, 108.2254.001, 108.2254.006, and 110.24CD.001, 110.24CD.003, 1102.C-0744, 1102.C-0744, 1102.C-0744, 1104.C-0350, 1104.C-0350, 1104.C-0350, 1104.C-0350, 1104.C-0350, and 1104.C-0350.</footnote>

<footnote>** Programs IDs for the HARPS observations we used in the analysis are: 060.A-9709, 072.C-0488, 083.C-1001, 084.C-0229, 086.C-0230, 087.C-0990, 088.C-0011, 089.C-0050, 090.C-0849, 091.C-0936, 092.C-0579, 093.C-0062, 094.C-0797, 095.C-0040, 096.C-0053, 60.A-9036, 60.A-9709, and 110.24BB.001.</footnote>

<footnote>*** Corresponding author; nicola.nari@lightbridges.es</footnote>

<footnote n="1">¹ https://exoplanetarchive.ipac.caltech.edu/</footnote>
able to reach an RV precision of the level of 1 m s⁻¹. To reach this goal, the instrument is temperature-controlled and pressure-stabilized. Over the years, this instrument has achieved many results of high scientific interest, such as the discovery of the first super-Earth (Santos et al. 2004). This instrument remained the state-of-the-art RV instrument until 2018 when the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO) (Pepe et al. 2014) began operations. The latter is installed at the 8.2 m Very Large Telescope (VLT) in Paranal, Chile. Again, ESPRESSO is temperature-controlled and pressure-stabilized. All these factors contribute to reaching a precision on a single measurement better than 10 cm s⁻¹ (Pepe et al. 2021). ESPRESSO was the basis for some of the most exciting discoveries in the field of exoplanets, for example, the confirmation of Proxima d, a small, sub-Earth planet orbiting around Proxima, (Suárez Mascareño et al. 2020; Faria et al. 2022). It has been involved in other ground-breaking discoveries in the field of exoplanetary atmospheres (Ehrenreich et al. 2020) and the mass determination of sub-Earth radius planets detected by transits (Demangeon et al. 2021). Combining the high-accuracy measurements of HARPS and ESPRESSO, while exploiting the long temporal baseline of the former, we can search for signals with amplitudes of below 1 m s⁻¹ at long orbital periods. According to the NASA Exoplanet Archive only 25 planets out of more than 5600 have a measured mass of less than 10 M⊕ and an orbital period of >50 d.

HD 20794 has been part of long RV surveys dedicated to the search for low-amplitude long-period signals around solar-type stars, with hundreds of nights of observations with HARPS and ESPRESSO spanning more than 20 years available. The recent publication of a candidate at ~640 d period (Cretignier et al. 2023) gave rise to a new campaign for a dense monitoring of the star. The low stellar activity level and the brightness of HD 20794 has made this target one of the most well-suited candidates for this purpose. Furthermore, HD 20794 is only 6.04 pc away, and long-period planets orbiting this star can become targets for future facilities dedicated to the characterization of exoplanetary atmosphere, both from the ground and space as ANDES at ELT (Palle et al. 2023), LIFE (Quanz et al. 2022), and HWO (Mamajek & Stapelfeldt 2024). The work is structured as follows. In Sect. 2, we report the observations obtained for HD 20794. In Sect. 3, we report stellar characteristics and state-of-the-art knowledge on the system. In Sect. 4, we report the modeling of the RV time-series. In Sect. 5, we discuss the results and implications of our work. In Sect. 6, we summarize our analysis and findings.

# 2. Observations

## 2.1. ESPRESSO

HD 20794 has been observed as part of the ESPRESSO Guaranteed Time Observations (GTO), within working group 1 (WG1), from October 2018 to March 2023. The main objective of WG1 is the research of Earth-like planets around G-to-M type stars within the habitable zone (HZ).

ESPRESSO (Pepe et al. 2021) is a fiber-fed spectrograph installed at VLT in Paranal, Chile, and it can be connected to each of the 8.2 m telescope units. It spans a wavelength range between 378.2 and 788.7 nm and has a resolving power of ~140 000 in high-resolution mode. Designed to achieve 10 cm s⁻¹ precision in RV for solar-type stars, ESPRESSO operates within a temperature- and pressure-controlled vacuum

vessel for best instrumental stability. The simultaneous reference technique (Baranne et al. 1996) using a Fabry-Pérot etalon (Wildi et al. 2010; Pepe et al. 2021) allows us to correct for any residual instrumental drift up to 10 cm s⁻¹.

The ESPRESSO data reduction software (DRS) provides extracted and wavelength-calibrated spectra, with different byproducts, such as the cross-correlation function (CCF), CCF-derived RV, and telemetry data. The CCF technique measures RV by convolving stellar spectra with a mask from a similar spectral type template star (Fellgett 1953; Baranne et al. 1979, 1996). For this study, we derived RVs with sBART (Silva et al. 2022), a template matching tool within a semi-Bayesian framework. sBART derives RVs by comparing each spectrum to a template spectrum built from observations.

In June 2019, an intervention on the ESPRESSO fiber link improved the efficiency by 50%, introducing an offset in the data. We treated observations before BJD 2458 721 as a separate dataset. An interruption of operations due to COVID-19 and a lamp change introduced an additional offset, corrected in DRS version 3.0.0, eliminating the need for additional dataset divisions. We obtained 34 observations over 6 nights before the intervention and 661 observations over 65 nights after. We discarded pre-intervention data due to its limited usefulness.

After nightly binning and excluding outliers, that is, nights with an error on the measurement exceeding three times the median error or nights with a difference from the median RV larger than three times the standard deviation of the measurements, we retained 63 nights, covering 1307.66 days from BJD 2458 721.89 to BJD 2460 029.49. The root mean square (RMS) of the post-intervention dataset, referred to as E19, is 0.84 m s⁻¹, with a mean photonic error of 0.10 m s⁻¹ per measurement. Nightly binning reduced the mean photonic error to 0.03 m s⁻¹. This is attributed to binning multiple measurements per night, usually ten or more. However, nightly stacking does not reduce the instrumental noise on a timescale beyond the intra-night timescale. Due to the brightness of HD 20794, we used short exposure times of 60 seconds.

## 2.2. HARPS

HARPS (Mayor et al. 2003) is a fiber-fed spectrograph installed at the Observatory of La Silla, in Chile. It started operations in 2003 and is now one of the most successful ground-based facilities for the discovery of exoplanets. Designed to obtain long-term stability, the instrument has shown a precision on a single measurement better than 50 cm s⁻¹ for bright stars and a long-term stability of 1 m s⁻¹ for quiet targets such as HD 20794 (Cretignier et al. 2021a, 2023). The spectrograph is equipped with two fibers and is mechanically stabilized. The instrument operates within a temperature- and pressure-controlled vacuum vessel to improve instrumental stability. One of the two fibers provides starlight to the instrument, while the other one is used for simultaneous calibrations with a Thorium-Argon lamp or Fabry-Pérot etalon as reference. It covers a spectral range of 378–691 nm and has a spectral resolving power of R ~ 115 000.

The instrument underwent a fiber update in 2015, which introduced an offset in the RV (Lo Curto et al. 2015). For this reason, we consider the data taken before and after the update as independent datasets (henceforth, H03 and H15) and consider them as if the data were acquired by different instruments. We have applied the ESPRESSO pipeline version 3.2.0, which has been adapted to HARPS to reduce HARPS spectra. The pipeline provides extracted and wavelength-corrected spectra. The new pipeline automatically corrects for drifts due to lamp changes
and aging and for a linear trend visible in the H15 dataset after BJD 59 915 in the previous pipeline. HD 20 794 has been part of an intensive observational campaign to confirm the planetary nature of the long-period candidate announced in Cretignier et al. (2023). In our analysis, we have considered HARPS spectra available on DACE<fnref n="2" /> (Buchschacher et al. 2015). We collected 765 visits of HD 20 794 made by HARPS, 531 nights before the fiber link update, and 234 nights after the fiber link update. For our analysis we extracted the RV with the line-by-line (LBL) (Dumusque 2018) post-processing pipeline YARARA (Cretignier et al. 2021a, 2023). The YARARA pipeline automatically rejects some nights due to anomalous CCF and S/N anomalies. We removed those nights from our analysis and remained with 514 nights binned from 6224 single spectra for H03 and 232 nights binned from 1409 spectra for H15.

YARARA improves the precision of RV by applying a correction to the 1-D spectra derived from the pipeline of the instrument and normalized by RASSINE (Cretignier et al. 2021b). RASSINE is a Python code for the normalization of merged 1D spectra. YARARA corrects systematics such as ghosts, telluric lines, stitching effects, cosmic rays, and more. The outcome from this first correction is referred to as YV1, following the definition of Cretignier et al. (2023). A second correction, named YV2, is then performed in the time domain. This correction consists of linearly decorrelating the RVs from the principal component coefficients obtained from a PCA performed on the SHELL spectral representation (Cretignier et al. 2022), in addition to the principal component coefficients obtained from a PCA performed on the LBL RVs (Cretignier et al. 2023). The SHELL spectral representation is a representation for the spectra where we consider the flux and its gradient ($\left(\frac{df_0}{dt}, f\right)$, instead of the flux and the wavelength ($\lambda, f$).

For our analysis, we used YV2 RVs for HARPS before fiber update (H03), and YV1 RVs for HARPS after fiber update (H15). This decision comes from the fact that YV2 is less effective when applied to H15 due to fewer observations and inhomogeneous sampling. After nightly binning observations and applying a $3\sigma$ cut to remove outliers, with the same procedure we followed for ESPRESSO, we retain 512 and 231 nights available for H03 and H15, respectively. Thanks to the YARARA correction we can derive HARPS RVs with a root mean square of the binned data of 1.14 and 1.21 m s⁻¹ for H03 and H15, respectively. The mean error in the measurement of binned points is 0.11 m s⁻¹ for H03 and 0.12 m s⁻¹ for H15. Here, we are considering only the photonic error on the binned measurements. We show the RVs time series in Fig. 1.

## 2.3. TESS

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) observed HD 20794 with a 2 min cadence in sectors 3, 4, 30, and 31. The data were processed with the SPOC (Science Processing Operations Center) pipeline (Jenkins et al. 2016) and the preliminary research for transiting planets was done with an adaptive, wavelet-dependent algorithm for transit (Jenkins et al. 2010). We analyzed the four time series independently and all together. We excluded flux measurements from sector 4 taken between BJD 2 458 420 and 2 458 424 for a systematic present in the data. We applied a BLS periodogram (Kovács et al. 2002; Hartman & Bakos 2016) to the time series, considering the four sectors independently and together, but we did not find any evidence of a transit. The RMS of TESS measurements is 0.08 ppt

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** RV dataset of HD 20794 for the different datasets. The HARPS RVs have been extracted using YARARA (Cretignier et al. 2021a, 2023). For ESPRESSO, RVs have been extracted with sBART (Silva et al. 2022).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Tess photometry and TESS target pixel file for HD 20794. Panel (a): TESS photometry for HD 20794. We do not find any evidence of transits in the dataset. Panel (b): TESS target pixel file for HD 20794. No nearby sources are detected by the Gaia DR3 catalog up to a contrast magnitude of +6 within the TESS aperture.</figcaption></figure>

and it goes down to 0.044 ppt when we bin the data in 30 min bins. TESS photometry is visible in Fig. 2a. We have used the publicly available code **tpfplotter**<fnref n="3" /> (Aller et al. 2020) to plot the target pixel file. The target pixel file is visible in Fig. 2b. Any other sources in the field are fainter by at least six magnitudes.

# 3. HD 20794

## 3.1. HD 20794: stellar parameters

HD 20794 ($\alpha$ = 03:19:55.669, $\delta$ = -43:04:11.29) is a bright star (V-mag ~4.34) of the spectral type G6V, located at a distance

<footnote n="2">² https://dace.unige.ch</footnote>

<footnote n="3">³ https://github.com/jlillo/tpfplotter</footnote>
<table><caption>**Table 1.** Stellar parameters of interest for HD 20794.</caption><thead><tr><th>Parameter</th><th>HD 20794</th><th>Ref</th></tr></thead><tbody><tr><td>α</td><td>3:19:55.67</td><td>1</td></tr><tr><td>δ</td><td>-43:04:11.29</td><td>1</td></tr><tr><td>Parallax (mas)</td><td>165.524 ± 0.078</td><td>1</td></tr><tr><td>d (pc)</td><td>6.0414 ± 0.0028</td><td>1</td></tr><tr><td>μ<sub>α</sub> cos δ (mas yr<sup>-1</sup>)</td><td>3035.017 ± 0.081</td><td>1</td></tr><tr><td>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>726.964 ± 0.101</td><td>1</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>5368 ± 64</td><td>0</td></tr><tr><td>log<sub>10</sub> g (cgs)</td><td>4.39 ± 0.03</td><td>0</td></tr><tr><td>Spectral type</td><td>G6V</td><td>4</td></tr><tr><td>RV (km s<sup>-1</sup>)</td><td>87.8820 ± 0.0001</td><td>1</td></tr><tr><td>[Fe/H] (dex)</td><td>-0.42 ± 0.02</td><td>0</td></tr><tr><td>M<sub>★</sub> (M<sub>☉</sub>)</td><td>0.79 ± 0.01</td><td>0</td></tr><tr><td>R<sub>★</sub> (R<sub>☉</sub>)</td><td>0.93 ± 0.03</td><td>0</td></tr><tr><td>log<sub>10</sub> R'<sub>HK</sub></td><td>-4.98 ± 0.02</td><td>3</td></tr><tr><td>B (mag)</td><td>5.130 ± 0.014</td><td>2</td></tr><tr><td>V (mag)</td><td>4.336 ± 0.009</td><td>2</td></tr><tr><td>L (L<sub>☉</sub>)</td><td>0.6869 ± 0.0026</td><td>5</td></tr><tr><td>ξ (km s<sup>-1</sup>)</td><td>0.555 ± 0.050</td><td>0</td></tr><tr><td>Inner optimistic HZ (au)</td><td>0.6322 ± 0.0020</td><td>0</td></tr><tr><td>Inner conservative HZ (au)</td><td>0.8329 ± 0.0033</td><td>0</td></tr><tr><td>Outer conservative HZ (au)</td><td>1.4502 ± 0.0086</td><td>0</td></tr><tr><td>Outer optimistic HZ (au)</td><td>1.5236 ± 0.0091</td><td>0</td></tr></tbody></table>

**Notes.** References: 0 - This work, 1 - Gaia Collaboration (2020), 2 - Høg et al. (2000), 3 - Buchschacher & Alesina (2019), 4 - Keenan & McNeil (1989), 5 - Gaia Collaboration (2018).

of 6.04 pc from the Sun (Gaia Collaboration 2020). We give a list of stellar parameters in Table 1. Surface gravity $log_{10} g$, the effective temperature $T_{eff}$, and the metallicity $[Fe/H]$ were derived directly from the analysis of the ESPRESSO spectra with the ARES + MOOG method described in Sousa et al. (2015). This method combines the ARES code (Sousa et al. 2007) for the derivation of the equivalent width of lines, with the MOOG model for the atmospheric derivation of abundances (Sneden 1973). ARES + MOOG is a specific realization of the equivalent width method for deriving stellar parameters. The equivalent width of a selected number of stellar lines is derived directly from the spectrum and it is used together with an atmospheric model to calculate the individual line abundances. The stellar parameters are found when the ionization and excitation balances of all lines are achieved; otherwise, the atmospheric model is changed to refine the fit. The ARES + MOOG method is designed to be automatically applied to a large sample of stars, so the line list was set to be wide and stable at the same time (Sousa et al. 2008). This method is more effective when the star analyzed is similar to the Sun. The ARES code derives the equivalent width and MOOG is used to calculate line abundances in the context of local thermodynamic equilibrium (LTE) through an interpolation with a grid of Kurucz Atlas-9 parallel model atmosphere (Kurucz 1993). A total of 14 spectra were added to reach the S/N of ~2500 necessary to extract the parameters. The stellar mass, $M_{★}$, and the stellar radius, $R_{★}$, are estimations coming from a calibration based on the derived spectroscopic parameters: $log_{10} g$, $T_{eff}$, and $[Fe/H]$ (Torres et al. 2010).

The star has color B-V ~ 0.7 (Høg et al. 2000), and a mass of 0.79 ± 0.01 M☉. We measure a $log_{10} R'_{HK} = -4.98 ± 0.02$. Another source for deriving the stellar parameters is the

eDR3 Gaia catalog (Gaia Collaboration 2021). The luminosity of the star is 0.6869 ± 0.0026 L☉ (Gaia Collaboration 2018). We can use the stellar mass and luminosity to derive the HZ of the system following Kopparapu et al. (2014). We refer to the HZ limits defined in Kopparapu et al. (2014) as “recent-Venus” and “early-Mars” as the optimistic inner and outer edge of the HZ, respectively. We refer to the HZ limits defined as “runaway greenhouse” and “maximum greenhouse” as the conservative inner and outer edge of the HZ. The conservative HZ for HD20794 spans between 0.833 ± 0.003 au and 1.450 ± 0.008 au, while the optimistic HZ is comprised between 0.632 ± 0.002 au and 1.524 ± 0.009 au. The corresponding orbital period for the inner optimistic and conservative HZ is 207 ± 2 d and 313 ± 3 d respectively, while for the outer HZ, we obtain values of 718 ± 8 d and 773 ± 9 d for conservative and optimistic calculations.

## 3.2. Planetary system

HD 20794 is known for hosting a multi-planetary system (Pepe et al. 2011; Feng et al. 2017; Basant et al. 2022; Cretignier et al. 2023). Through the years, multiple analyses have been conducted on this star. Different analyses use non-identical datasets and different methods for red-noise correction. As a result, they bring on outcomes that are not fully compatible. Considering the periods of the planets detected in the different works, we refer to the planets as the 18-d, 40-d, 89-d, 147-d, and 640-d planet.

Pepe et al. (2011) reports the discovery of three planets: HD 20794 b with a period of 18.315 ± 0.008 d and amplitude of 0.83 ± 0.09 m s⁻¹; HD 20794 c with a period of 40.114 ± 0.053 d and an amplitude of 0.56 ± 0.10 m s⁻¹; and HD 20794 d with a period of 90.31 ± 0.18 d and an amplitude of 0.85 ± 0.10 m s⁻¹. Feng et al. (2017) considered a larger dataset of HARPS observations, as compared to Pepe et al. (2011). They used the template matching tool TERRA (Anglada-Escudé & Butler 2012) to extract the RVs. A wavelength-dependent noise model was applied to the time series to correct spurious signals in the RVs. For the correction details, we refer to the original paper. Feng et al. (2017) confirmed the detection of planets at 18-d and 89-d, but did not find evidence of the planet at 40-d. Furthermore, they found an additional candidate at orbital period of 147.0 ± 1.1 d, with an amplitude of 0.69 ± 0.14 m s⁻¹. Feng et al. (2017) reported a significant signal at P = 331.4<sup>+5.1</sup><sub>-3.0</sub> d, but since this is so very close to a year and there is the presence of signals with the same periodicity in activity indicators, it does not allow for the RV signal to be claimed as a safe detection. Recently, Cretignier et al. (2023) revisited the system and considered the historical time series of HARPS before the fiber intervention the instrument underwent in 2015. The main change in this analysis was the usage of the YARARA tool to extract velocities (Cretignier et al. 2021a, 2023). The signals at 18 d and 89 d have been recovered with a high level of confidence. A new candidate is also recovered with a period of 644.6 ± 8.8 d and an amplitude of 0.61 ± 0.06 m s⁻¹. Cretignier et al. (2023) did not find evidence of any additional signals. The possible presence of an additional planet with an orbital period between 549 and 733 d is also reported in Basant et al. (2022). We report a summary of the results from the literature in Table B.1.

## 4. Analysis

Here, we report the analysis of the HARPS and ESPRESSO datasets. We refer to Appendix A for details about the methods
and tools used in our analysis to conduct parameter estimation and model comparison.

## 4.1. Stellar activity

We analyzed different proxies for activity, such as the full width at half maximum (FWHM), bisector timespan (BIS), S-index, $H_α$, Na I D, and the contrast of the CCF to retrieve the stellar characteristics. The periodicities we found in the activity proxies (if they were also found in RVs) would require specific attention to be paid to determine their nature, so that we could avoid misunderstanding stellar-related signals from signals that are Keplerian in origin.

When we modeled the activity, we considered both activity indicators derived from the spectra before and after the YARARA correction. YARARA can correct signatures of activity and we are interested in the inference of stellar parameters, such as the period of an eventual magnetic cycle and the period of stellar rotation. We corrected the activity indicators for changes in temperature of the echelle gratings of the instruments with an approach similar to Suárez Mascareño et al. (2023). We found evidence of a magnetic cycle in FWHM with a period of $3020_{-50}^{+111}$ d. We were also able to find evidence of a long-term magnetic cycle in the analysis of other activity indicators, such as the S-index, BIS, and contrast.

We recovered the rotation period from FWHM and BIS, for the latter in the YARARA-corrected dataset. We found a rotation period of $35.0_{-2.5}^{+3.2}$ d derived from FWHM and a rotation period of $38.8_{-2.6}^{+2.4}$ d derived from BIS. The latter is particularly important because this period is compatible with the period of the 40-d planet detected in Pepe et al. (2011). We conclude this signal is related to stellar activity. Considering the main focus of the analysis of HD 20794 is the analysis of RVs, we refer to Appendix E for a detailed description of the impact of stellar activity on the detection of exoplanets and the full analysis of the different activity indicators. Figure E.11 shows time series and generalized Lomb-Scargle (GLS) periodograms for different activity indicators extracted from spectra that were not corrected by YARARA. Figure E.12 shows the time series and GLS periodograms for different activity indicators from spectra corrected by YARARA.

## 4.2. Radial velocities analysis

Figure 1 shows the RVs of HD 20794 used in this analysis. The full RV dataset spans 7496 d, roughly corresponding to 20 years. After removing outliers and binning data as explained in Appendix A we have a total of 806 nights of observations, divided into 512 nights for HARPS before the fiber update (H03), 231 for HARPS after the fiber update (H15), and 63 nights for ESPRESSO after the intervention on the fiber link (E19). First, we conduct searches for planets in the ESPRESSO dataset alone. Then, we included the HARPS dataset directly corrected for YARARA to make a full analysis of the system.

## 4.3. ESPRESSO

We first consider the analysis of the system with the ESPRESSO dataset alone. We refer to Sect. 2.1 for a description of this dataset. We show the ESPRESSO dataset in Fig. 3a and the ESPRESSO dataset GLS periodogram in Fig. 3b.

The GLS periodogram of Fig. 3b shows the most prominent peak at 19.23 d. This peak does not show a false alarm probability (FAP) of <10%, but it is still remarkable that it stands as

<figure><img src="image_5.png" /><figcaption>Fig. 3. ESPRESSO observations of HD 20794. Panel a: E19 RVs for HD20794. Panel b: GLS periodogram of the ESPRESSO RVs for HD 20794. The main peak is around 19.23 d, a year alias of 18.3 d signal, and a second peak at 115 d, a year alias of the peak at ~87 d.</figcaption></figure>

precisely the 1 year alias of 18-d planet seen in all the literature works. We can see a double peak at 87.5 and 115 d, with one being the 1 year alias of the other. Again, even if the peak has a FAP of >10%, the signal at 87.5 d has a period comparable with the 89-d signal found in previous literature works.

The first model we considered is a model only concerning a zero-point of the velocity. Then, we tried a blind search for the planets but we do not find a convergent result. The following step is conducted as an informed search for a sinusoidal signal, with a uniform prior on the amplitude between 0 and $3 \text{ m s}^{-1}$. We consider a normal prior on the period centered at 18.3 d, the period of planet b, and a width of 0.2 d. This is necessary to avoid a conflict between the signal and its aliases in the determination of the period. We find an amplitude for the signal of $0.51 \pm 0.15 \text{ m s}^{-1}$, with a period of $18.33_{-0.04}^{+0.03}$ d.

To compare the different models, we used the criterion of the natural logarithm of the evidence (lnZ) associated with a model (Appendix A). The difference in evidence with the no-planet model is $\Delta \ln Z = +1.7$. This is not sufficient to claim detection. Once we subtract the sinusoidal we have fitted in the one-planet model, we obtain the GLS periodogram of the residuals shown in Fig. 4a. In the GLS periodogram, we can see a weak peak at the period of ~87 d and its 1-year alias at 113 d. These two peaks have a FAP that is slightly lower than 10%.

When we modeled the second signal, we considered an informed search on the period. We followed the same strategy of the one-planet model, performing an informed search with an additional sinusoidal. For this second sinusoidal we consider a normal prior on the period centered at 89.5 d, with a width of 5 d. We kept the same prior as before for the signal at 18.3 d. The model with two sinusoidals has a $\Delta \ln Z = +4.1$ with respect to the model with no planets and a $\Delta \ln Z = +2.4$ to the model
<figure><img src="image_6.png" /></figure>

**Fig. 4.** Analysis of the ESPRESSO residual time series. Panel a: GLS periodogram of the residuals of ESPRESSO RVs for HD 20794 once we subtract the signal at 18.3 d. We can see a double peak at 87 and 113 d, one peak being the 1 year alias of the other. Panel b: GLS periodogram of the residuals of ESPRESSO RVs for HD 20794 after subtracting the signals at 18.3 d and 87 d. We can see some peaks at 4 and 40 d. The peak at 4 d probably comes from the sampling, while the peak at 40 d could be a signature of the stellar rotation period.

with one sinusoid. The model is slightly favored both on the one-sinusoidal model and the no-planet model. For the additional signal we find an amplitude of $K2 = 0.53 \pm 0.12$ m s⁻¹ and a period of $P2 = 87.0_{-0.9}^{+1.0}$ d. We plot the GLS periodogram of the residuals once we also subtract the signal of the second planet in Fig. 4b.

In the GLS periodogram, we can see two peaks with a FAP <10%, at 4.35 d and ~40 d. The period at 40 d is of particular interest because it corresponds to the orbital period of the 40-d planet found in Pepe et al. (2011), and, following our analysis of the activity, it is likely caused by stellar rotation. We followed the same approach as before: fitting a sinusoidal with a normal prior on the period of each of the two signals in addition to the previous ones. In the case of the 4.35-day signal, we find worse evidence in comparison to the two-sinusoidal model, with a $\Delta \ln Z = -0.3$ in disfavor of the three-sinusoidal model. We discarded this model. In the case of the 40-d signal, we find evidence of a $\Delta \ln Z = +0.8$ compared to the two-planet model, which means a slight improvement. The amplitude related to the third signal is $K = 0.46_{-0.15}^{+0.13}$ m s⁻¹, with a period of $40.31_{-0.18}^{+0.24}$ d. Considering that the stellar origin is the most likely origin of the 40-d signal, we tried to model the signal at 40 d with a GP, but the fit did not converge, so we discarded that model. We tried different kernels as SHO and MEP (see Appendix E for details on the definition of the two kernels). Even if we were not able to model this signal with a GP, we considered the signal as likely to be coming from stellar rotation, because we do not see any signature of it in the full analysis. Rotation-related signals are not stable over long timespan and this is the case for the 40-d signal.

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

**Fig. 5.** Phase-folded plots and model of planets recovered in the ESPRESSO dataset. Panel a: phase-folded plot of 18 planet (on the left) and 87-d planet (on the right). Panel b: two-planet model for the ESPRESSO dataset alone. We subtracted a zero-point velocity from the dataset.

The paucity of measurements in the ESPRESSO dataset alone can also make the GP fit a non-trivial task, considering the low amplitude of activity signals in HD 20794. Taking into account what we know from the analysis of activity indicators and the slight improvement in terms of $\ln Z$ once we added a third signal, we considered the best model for the ESPRESSO dataset alone to be the two-planet model with planetary periods of $18.33_{-0.04}^{+0.03}$ d and $87.0_{-0.9}^{+1.0}$ d. For this model, we found a minimum mass for the planet at $18.33$ d of $1.91_{-0.47}^{+0.45}$ M⊕ and a minimum mass for the planet at $87.0$ d of $3.11 \pm 0.72$ M⊕. We show a phase-folded plot of the planets in Fig. 5a.

A plot of the two-planet model imposed on the ESPRESSO dataset, where a zero-point velocity was subtracted, is shown in Fig. 5b. A full analysis of the system based on the ESPRESSO dataset alone would require a larger dataset and a longer timespan of observations considering we are not able to properly recover the planets, even if the photon noise we reach is at the level of $10$ cm s⁻¹ for a single exposure.

## 4.4. HARPS + ESPRESSO

Once we had explored the possibilities of the ESPRESSO dataset alone, we added the HARPS dataset corrected by YARARA to our analysis. The full dataset is visible in Fig. 1. The offset between H03 and H15 was calculated for this target at $17.0 \pm 1.7$ m s⁻¹ (Lo Curto et al. 2015). In our analysis, we fit for the offset as a free parameter. The instrument efficiency has improved by 33–40% after the change of fibers, while the
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** GLS periodogram of the full dataset of RVs for HD 20794. The strongest peak is at 18.3 d. Also, peaks at 89.6 and 650 d are visible. The periods found in Cretignier et al. (2023) are highlighted by red vertical lines.</figcaption></figure>

resolution remained constant (Lo Curto et al. 2015). The GLS periodogram of the full dataset is visible in Fig. 6. A peak at 18.3 d is visible. Also, peaks at 89.6 and 650 d are visible in the GLS periodogram, all with FAP < 0.1%.

### 4.4.1. One-planet model

For the analysis of the presence of planetary signals in the data, we followed the same approach we followed for the ESPRESSO dataset alone, namely: adding the signals one by one. Considering the larger significance of the peaks and the large number of observations, we go first for a blind search, considering a model with a single sinusoidal with a uniform prior for the period in log space between 2 and 3000 d. The blind search in log space for the period of the signal avoids sampling with the same density regions of the parameter space at short and long orbital periods. We have a uniform prior on the amplitude between 0 m s⁻¹ and 3 m s⁻¹. We ran our analysis multiple times to investigate the stability and robustness of the result.

We find that the result on lnZ in the blind search is not constant within the error associated with the measurement by our script. This issue with the calculation of the evidence in nested sampling was reported in Nelson et al. (2020). To mitigate this effect, we ran each model five times. We discuss the mean of the best three runs, the best value and the standard deviation of the measurements of lnZ to infer the best model. The decision to consider only the best three runs comes from the model being intrinsically degenerate trying to fit a multiple planetary system with a single planet model that for construction does not have any preference among the signals. Through this procedure, we mitigate the impact of lnZ outliers in the calculation of the mean lnZ. In the multiple runs, we found two different solutions, the fit converging to both 18.3 and 89.6 d. We consider as a reference in evidence the evidence associated with a model with no planets. When we consider the mean of the evidence of the one-planet model compared to the no-planet model we find a Δ lnZ = +38.0 in favor of the model with one planet. If we compare the best result of the one-planet model with the best result of the no-planet model, the first is favored by a Δ lnZ = +40.9. The standard deviation of the repeated measurement of lnZ in the one-planet model is 3.0, while for the no-planet model, the standard deviation is 0.2. The difference in lnZ is more than 30. We summarize in Table 2 the different models in the analysis with their significance. We consider the best solution for our one-planet model

<table><caption>**Table 2.** Evidence of different models for HD 20794.</caption><thead><tr><th>Model</th><th>Δ lnZ-mean</th><th>Δ lnZ-best</th></tr></thead><tbody><tr><td>No-planet</td><td>-108.3</td><td>-108.6</td></tr><tr><td>One-sinusoidal</td><td>-69.3</td><td>-60.2</td></tr><tr><td>Two-sinusoidal</td><td>-38.1</td><td>-34.5</td></tr><tr><td>Three-sinusoidal</td><td>-8.7</td><td>-7.1</td></tr><tr><td>Three-Keplerian</td><td>0.0</td><td>0.0</td></tr><tr><td>Four-Keplerian</td><td>+2.0</td><td>+6.9</td></tr></tbody></table>

**Notes.** We consider as the zero-point of each column separately the lnZ associated with a three-Keplerian model. The values in the table reflect the comparison between different models for the planetary system around HD 20794.

the solution with the best lnZ. We find a planet with an amplitude of K = 0.61 ± 0.05 m s⁻¹ and a period of P = 18.314 ± 0.002 d. The GLS periodogram of the residuals after subtracting the signal at 18.313 d is visible in Fig. 7a.

In the GLS periodogram of the residuals, it is possible to see a strong peak at 89.66 d, which is the period of a planet found in all previous works on the system. There is also a second prominent peak at ~650 d and a peak at ~111 d.

### 4.4.2. Two-planet model

We continue in the analysis and test a model comprising a second sinusoidal signal. In this model, we carried out a blind search for both the periods of the two planets, with overlapping priors in log space between 2 and 3000 d. We keep the same setup for the other parameters. The two-planet model is favored in terms of the mean evidence of repeated runs by a Δ lnZ = +32.2 compared to the one-planet model, and by a Δ lnZ = +70.2 compared to the no-planet model. If we compare the best models, we have a Δ lnZ = +33.3 with respect to the one-planet model and a Δ lnZ = +74.2 with respect to the no-planet model in favor of the two-planet model. The standard deviation on the evidence of the repeated measurement is 3.7. In the multiple runs, we find results converging to three different planets, at 18.3, 89.6, and 650 d. Considering the best run in terms of lnZ we find a short-term signal with an amplitude K= 0.60 ± 0.05 m s⁻¹, and a period of P = 18.313 ± 0.002 d, and a second signal recovered with an amplitude K2 = 0.46 ± 0.05 m s⁻¹ and a period P2 = 89.65±0.09d. We show the GLS periodogram of the residuals in Fig. 7b, where we can see a strong peak at 651.79 d and two additional peaks at ~85 d and ~111 d, one the 1-year alias of the other. The long-period peak is compatible in period with the long-term planet found in Cretignier et al. (2023).

### 4.4.3. Three-planet model

The GLS periodogram of the residuals in the two-planet configuration indicates the presence of an additional planet at an orbital period of 651 d. This peak corresponds to the candidate planet detected in Cretignier et al. (2023). To model this third signal, we followed the same procedure as we did for the others. We considered three sinusoids with the same blind priors on the period, spanning in log space all periods between 2 d and 3000 d, and a uniform prior on the amplitudes between 0 and 3 m s⁻¹. The average evidence of the three-planet model is favored, compared to the two-planet model in terms of lnZ by Δ lnZ = +29.4, while the best models differ by Δ lnZ = +27.4, in favor of the three-planet model. The repeated runs have a standard deviation
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Analysis of the residuals time series after planetary signals subtraction. Panel a: GLS periodogram of the residuals after subtracting the 18.3 d signal. A peak at ~89 d and another one at ~650 d are visible. Panel b: GLS periodogram of the residuals after subtracting the 18.3 d and 89.6 d signals. A peak at ~650 d is visible. Panel c: GLS periodogram of the residuals after subtracting the 18.314, 650.9, and 89.65 d signals. We see some additional peaks in the periodogram of the residuals with FAP < 0.1% at periods of ~85 and 111 d.</figcaption></figure>

of 2.0. The large difference in the evidence compared to the two-planet model points toward the direction of this model as the best to describe the dataset. The periods of the planets are always recovered with the same three periods, even if the order those are recovered can differ from run to run. We refer to the planets at 18.3 d, 89.6 d, and 650 d as HD 20794 b, HD 20794 c, and HD 20794 d, respectively. For HD 20794 b, we find an amplitude of $K_b = 0.60 \pm 0.05$ m s⁻¹ and a period of $P_b = 18.314 \pm 0.002$ d. For HD 20794 c, we find an amplitude of $K_c = 0.52 \pm 0.05$ m s⁻¹ and a period of $P_c = 89.65 \pm 0.08$ d. For HD 20794 d, we find an amplitude of $K_d = 0.46 \pm 0.05$ m s⁻¹ and a period of $P_d = 650.9_{-4.9}^{+5.0}$ d. Once we subtract this model, we have the GLS periodogram of the residuals shown in Fig. 7c.

We see two peaks with FAP < 0.1% at 85.6 and ~111.7 d. These two signals are one the one-year aliases of the other. We tested a model with three Keplerians instead of three sinusoids

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Phase-folded plots: HD 20794 b (left); HD 20794 c (middle); and HD 20794 d (right).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Three-planet model for the full dataset (left) and zoom on the model on the ESPRESSO dataset (right).</figcaption></figure>

to compare with Cretignier et al. (2023). To sample for eccentricity and argument of the pericenter, we did not directly use the parameters, but combinations of them: $\sqrt{e} \cos(\omega)$ and $\sqrt{e} \sin(\omega)$. We used a normal prior centered at 0 with a width of 0.3. We imposed the condition $e \le 1$. The mean evidence of the Keplerian model was favored, compared to the model with circular orbits by a $\Delta \ln Z = +8.7$, with a $\Delta \ln Z = +7.1$ for the best model evidence. HD 20794 b has an amplitude $K_b = 0.62 \pm 0.05$ m s⁻¹ and a period $P_b = 18.314 \pm 0.002$ d; HD 20794 c has an amplitude $K_c = 0.50$ m s⁻¹ and a period $P_c = 89.67_{-0.10}^{+0.11}$ d; HD 20794 d has an amplitude $K_d = 0.57 \pm 0.07$ m s⁻¹ and a period $P_d = 647.5_{-2.8}^{+2.6}$ d. We can only put upper limits on $e_b < 0.13$ (84th percentile) for HD 20794 b, and $e_c < 0.16$ (84th percentile) for HD 20794 c; whereas for HD 20794 d, we find a value for the eccentricity from the best run $e_d = 0.45_{-0.11}^{+0.10}$. The eccentricity of HD 20794 is significantly different from zero, according to Cretignier et al. (2023). We see peaks at ~85 d and 111 d in the periodogram of the residuals again. In Fig. 8 we show the phase-folded plots of the three planets from the Keplerian model. In Fig. 9 we show the three-Keplerian model. A zoom on the model on the ESPRESSO dataset is also shown.

#### 4.4.4. Four-planet model

We applied a blind search for a model with four Keplerians to search for additional signals in the dataset. We considered a Keplerian model because we found a significative eccentricity for HD 20794 d. This model has a $\Delta \ln Z = +2.0$ compared to the
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** FIP periodogram for HD 20794. We recover three signals at the same periodicities we derived in our blind search. The FIP levels for the planets at 18.3, 89.7, and 657.6 d are cut because the numerical calculation gave us back an infinite value.</figcaption></figure>

three-planet model, with a standard deviation of 4.6. The best solution for the four-planet model has a $\Delta \ln Z = +6.9$ compared to the three-planet model. We do not see the fourth signal converging clearly to any of the peaks we see in Fig. 7c. The posterior distribution of the additional signal is shared between the signal at ~85 d, the signal at ~111 d, and additional long-term signals at ~1000 d and ~1400 d. The amplitude of this hypothetical signal would be at ~30 cm s⁻¹. Furthermore, the improvement in lnZ is not significant enough to claim a new detection. Investigating the presence of additional planetary candidates, at the level of 30 cm s⁻¹, requires additional observations to be thoroughly investigated.

### 4.4.5. False inclusion probability

The analysis of the ESPRESSO dataset and the HARPS dataset corrected by YARARA revealed the presence of at least three planets in the system. Statistical tools are available to assess the robustness of the result. We analyzed the false inclusion probability (FIP) (Hara et al. 2022) to assess the number and significance of planetary signals in the data. FIP has been demonstrated to be optimal to maximize true detections for a certain tolerance to false positives (Hara et al. 2024). We take a grid of frequency intervals with a fixed length. The k element of the grid is defined as $[\omega_k - \Delta\omega/2, \omega_k + \Delta\omega/2]$, where $\omega_k$ is defined as $\omega_k = k\Delta\omega/N_{\text{oversampling}}$ and $\Delta\omega$ is defined as $\Delta\omega=2\pi/T_{\text{obs}}$. Here, $T_{\text{obs}}$ is the total time of the observations, and we take a $N_{\text{oversampling}} = 5$. We considered a Keplerian model for the planets. In Fig. 10, we show the result of our FIP analysis. In the plot, we cut the levels of FIP for the 18.3 d, the 89.6 d, and the 657.6 d signals because the numerical result of the FIP is infinite (this points toward a strong confirmation for these two signals). We see some additional peaks in the FIP at ~111 d and ~1000 d, but not exceeding the 1% FIP, that translates in a confidence in the detection lower than 1%.

Considering the issue with the stability of the evidence calculation of the nested sampling we reported in our previous analysis we repeated the FIP calculation up to five times to check the consistency through the runs. We always recover each of the three signals with confidence >99.999% in every run. The peaks related to the planets at 111 d and 1000 d exceed the 1% FIP limit just once over multiple runs. The hypothetical characteristics of the 111-d signal are discussed in Appendix D.

### 4.4.6. Final results

We conducted multiple runs to marginalize the intrinsic instability of our sampling method. The dispersion on lnZ due to the sampler is one order of magnitude less than the difference of lnZ we obtain every time we add planets b,c, and d, respectively. The FIP indicates the presence of the same three planets we have found in the blind search. An eventual additional candidate arises once we subtract the signal of the three known planets from our dataset, but both a blind search and the FIP analysis do not corroborate its presence solidly. Considering these elements we take as our best model the one with three Keplerians. Once we have assessed the significance of the three signals in our analysis framework, we can relax the constraints on searching the periods blindly. We do this to find the best parameters for the model. We want to define disjoint regions of the parameter space where to search for different signals, to avoid the nested sampling method to struggle between competing signals. We considered a normal prior on the periods based on our previous knowledge derived from the blind-search and the FIP analysis. In principle, the decision on the width of the prior can influence the result of the analysis. To study this possible source of noise, we performed different analyses with a width of the normal prior equal to 5, 10, and 20% the period of the planet. We chose, as the starting periods for the planets, 18.3, 89.6, 650 d for HD 20794 b, c, and d, respectively. We determined our results on the periods of the planets different at less than $0.1\sigma$. We concluded that the prior definition, on the limits of a prior larger than 5% of the period and narrower than 20%, does not influence the result of our analysis. We refer to the results obtained with a prior width of 10% of the period value. We obtain for HD 20794 b an amplitude of $K_b = 0.614 \pm 0.048 \text{ m s}^{-1}$, with a period $P_b = 18.3140 \pm 0.0022 \text{ d}$, and an upper limit on eccentricity, $e_b < 0.13$. This corresponds to a planet with a minimum mass of $M_{p\text{sini}} = 2.15 \pm 0.17 \text{ M} \oplus \text{orbiting at a distance of } 0.12570^{+0.00053}_{-0.00052} \text{ au}$. For HD 20794 c we obtain an amplitude of $K_c = 0.502^{+0.048}_{-0.049} \text{ m s}^{-1}$, with a period of $P_c = 89.68 \pm 0.10 \text{ d}$, and an upper limit on eccentricity $e_c < 0.16$. This corresponds to a planet with a minimum mass of $M_{p\text{sini}} = 2.98 \pm 0.29 \text{ M} \oplus \text{orbiting at a distance of } 0.3625^{+0.0015}_{-0.0016} \text{ au from the star. HD 20794 d has an amplitude of } K_d = 0.567^{+0.064}_{-0.067} \text{ m s}^{-1}$, a period of $P_d = 647.6^{+2.5}_{-2.7} \text{ d}$, an eccentricity of $e_d = 0.45^{+0.10}_{-0.11}$, a minimum mass of $M_{p\text{sini}} = 5.82 \pm 0.57 \text{ M} \oplus$, and an orbital distance of $1.3541 \pm 0.0068 \text{ au}$. Once we subtract the model with three Keplerian, the RMS of the dataset goes from 1.14 to 0.93 m s⁻¹. The RMS for H03 goes from 1.14 to 0.93 m s⁻¹, the RMS for H15 goes from 1.21 to 0.99 m s⁻¹, and the RMS for E19 goes from 0.84 to 0.72 m s⁻¹. In Table 3, we summarize the main results obtained in our analysis.

## 5. Discussion

### 5.1. Planetary system

Previous works from the literature assess the presence of a multi-planetary system orbiting around HD 20794. In our analysis, we exploited a larger number of HARPS observations that were collected to investigate the 640-d candidate reported in Cretignier et al. (2023). Furthermore, we added the ESPRESSO dataset, and exploited a new tool for extracting velocities after the correction of activity and systematics at the spectral level, YARARA. We considered a consequential blind search for new signals to avoid any possible bias from previous works on the same subject. With this approach, we detected three significant signals, corresponding to the same planets detected in Cretignier et al.
<table><caption>**Table 3.** Parameters for planets b,c,d derived in a model with norm priors.</caption><thead><tr><th>Parameter</th><th>HD 20794 b</th><th>HD 20794 c</th><th>HD 20794 d</th></tr></thead><tbody><tr><td>T0 - 2 450 000 [d]</td><td>10 383.26 ± 0.54</td><td>10 356.2 <sup>+5.2</sup><sub>-4.5</sub></td><td>9933 <sup>+22</sup><sub>-17</sub></td></tr><tr><td>P (d)</td><td>18.3140 ± 0.0022</td><td>89.68 ± 0.10</td><td>647.6 <sup>+2.5</sup><sub>-2.7</sub></td></tr><tr><td>K [ms<sup>-1</sup>]</td><td>0.614 ± 0.048</td><td>0.502 <sup>+0.048</sup><sub>-0.049</sub></td><td>0.567 <sup>+0.067</sup><sub>-0.064</sub></td></tr><tr><td>e</td><td>0.064 <sup>+0.065</sup><sub>-0.046</sub></td><td>0.077 <sup>+0.084</sup><sub>-0.055</sub></td><td>0.45 <sup>-0.11</sup><sub>+0.10</sub></td></tr><tr><td>m<sub>p</sub>sini [M<sub>⊕</sub>]</td><td>2.15 ± 0.17</td><td>2.98 ± 0.29</td><td>5.82 ± 0.57</td></tr><tr><td>a [AU]</td><td>0.12570 <sup>+0.00052</sup><sub>-0.00053</sub></td><td>0.3625 <sup>+0.0015</sup><sub>-0.0016</sub></td><td>1.3541 ± -0.0068</td></tr></tbody></table>

**Notes.** A Keplerian model is considered here.

(2023). HD 20794 b has an orbital period of 18.3140 ± 0.0022 d and its RV signal has an amplitude of 0.614 ± 0.048 m s⁻¹. These values correspond to a planet with a minimum mass of 2.15 ± 0.17 M⊕, orbiting the star at a distance of 0.12570 <sup>+0.00053</sup><sub>-0.00052</sub> au. We can only put an upper limit on the eccentricity with e<sub>b</sub> < 0.13 HD 20794 c has an orbital period of 89.67 ± 0.10 d and its RV signal has an amplitude of 0.502 <sup>+0.048</sup><sub>-0.049</sub> m s⁻¹. These values correspond to a planet with a minimum mass of 2.98 ± 0.29 M⊕, orbiting the star at a distance of 0.3625 <sup>+0.0015</sup><sub>-0.0016</sub> au. Again, we can only put an upper limit on the eccentricity, with e<sub>c</sub> < 0.16. HD 20794 d has an orbital period of 647.5 <sup>+2.5</sup><sub>-2.7</sub> d and its RV signal has an amplitude of 0.567 <sup>+0.064</sup><sub>-0.067</sub> m s⁻¹. For this planet, we find an eccentricity of e<sub>d</sub> = 0.45 <sup>+0.10</sup><sub>-0.11</sub>. These values correspond to a planet with a minimum mass of 5.82 ± 0.57 M⊕, orbiting the star at a distance of 1.3541 ± 0.0068 au. The orbit of the long-period planet is compatible with an eccentric one. This result is in agreement with the result obtained in Cretignier et al. (2023). In that work, it was derived an eccentricity for the outer planet e<sub>d</sub> = 0.40 ± 0.07. The larger dataset and the addition of ESPRESSO point toward a similar solution. All the planets of the system have a minimum mass compatible with a super-Earth. Still, we cannot infer the true mass of the planets without information on the inclination of the orbit with respect to the line of sight. Masses of HD 20794 b and HD 20794 c are also compatible with an Earth-like scenario. HD 20794 d could also be a mini-Neptune with a non-negligible H/He atmosphere. The fact that we cannot access the radii of these planets does not permit us to infer their nature. The orbital period of HD 20794 d resides both in the optimistic and conservative HZ. In Fig. 12 we show the position of the planets in the conservative and optimistic HZ. This is an interesting result because we do not have many examples of planets with M < 10 M⊕ with mass measurement from RVs in the HZ of Sun-like stars. Due to the high eccentricity of this signal we need to take into account the large variation in stellar flux received by the planet during its orbit. At the apoaster, the distance of the planet to the star is 1.96<sup>+0.13</sup><sub>-0.16</sub> AU, while at the periastron the planet-star distance is 0.75<sup>+0.15</sup><sub>-0.13</sub> AU. This means the orbit of the HD 20794 d crosses the HZ as shown in Fig. 11 for the best-fit solution. The stellar flux at the periastron is almost seven times stronger than the stellar flux at the apoaster. HD 20794 d spends ~59% of its orbit inside the optimistic HZ, and ~38% of its orbit inside the conservative HZ. Biasiotti et al. (2024) investigated the habitability of planets crossing the HZ for just a fraction of their orbit, as is the case for GJ 514 b (Damasso et al. 2022). This planet has a similar eccentricity to HD 20794 d. The level and number of tests of their analysis are beyond the scope of this work. In their work, we can see how habitability could be possible for planets

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Position of the HZ relative to the elliptical orbit of HD 20794 d. It is possible to see how the planet crosses the HZ both in the optimistic and conservative boundaries for a large proportion of time.</figcaption></figure>

on high-eccentric orbits based on stellar and planetary parameters such as stellar age, CH₄ abundance in the atmosphere, axis obliquity, ocean fraction, rotation period of the planet, and other properties. The possibility of maintaining habitable conditions even in a highly eccentric orbit raises the interest in investigating the HD 20794 system in future studies. Figure 13 offers a plot of the mass-period relation for detected planets with masses measured from RVs. Planets orbiting HD 20794 populate the lower edge of the diagram in mass at their orbital period. This brings us to highlight the importance of long-term and high-cadence campaigns for the characterization of low-mass signals, especially in the outer regions of planetary systems. Upcoming surveys as the Terra Hunting Experiment (Hall et al. 2018), will follow the same approach, paving the way to the characterization of habitable terrestrial planets orbiting Sun-like stars.

In our analysis, we did not find any evidence of the presence of the planet at 147 d reported in Feng et al. (2017). We tried an informed search for the planet with priors centered around the period of the missing planet. The addition of the signal worsened the evidence of the model compared to the three-planet model. We can see some power excess in the GLS periodogram of the ESPRESSO dataset alone at the period of ~40 d. We did not
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** HZ for HD 20794. Following Kopparapu et al. (2014) we mark the conservative HZ in green, while the optimistic HZ is marked in yellow. The red part of the plot considers all the periods remaining on the interior of the optimistic HZ, while the blue part represents the periods wider than the outer edge of the optimistic HZ.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Mass-period relationship for planets orbiting around HD 20794 compared to the planets detected in literature with the RV method. The planets orbiting HD 20794 are among the lightest planets detected for G-type stars at their respective orbital period, well below the 1 m s⁻¹ threshold.</figcaption></figure>

recover the same signal when we consider the full dataset. We found evidence of a signal of the same period related to stellar rotation in some activity indicators (see Sect. 4.1). The possibility the 40-d signal is related to rotation can explain the fact the signal is not persistent in the full dataset.

## 5.2. Detection limits

We detected three planets orbiting HD 20794. The possibility of detecting a planet orbiting around a star depends on a list of factors: the quality of data, the quality of correction for stellar activity, and the quality of sampling. We want to investigate the regions of the parameter space where we are not sensitive to the detection of planets. To determine our sensitivity limits we do an injection-recovery test similar to the one in Suárez Mascareño et al. (2023). We base our knowledge of the planetary system on the result of Table 3. We subtracted planetary signals with the same ephemerides as the ones we found for HD 20794 b, c, and d. We did not consider any additional correction for activity outside of the YARARA correction for HARPS. Figure 14 gives a plot of the results of the injection-recovery test for HD 20794, where we show the sensitivity depending on the mass and period of the hypothetical planets. To perform this test we injected 90,000 sinusoidal signals with periods between 2 d and 8000 d taken from a grid of 300 different periods uniformly spaced in logarithm and a grid of 300 masses between 0.1 and 15 M⊕ uniformly spaced in logarithm. For each injection, we considered a GLS periodogram of the resulting dataset. We considered to be sensitive to injected signals when we can find a peak in the GLS periodogram with FAP < 1% at the injected period. We needed to subtract an additional Keplerian signal from the dataset fitted with a norm prior on the period centered at 111 d. This is necessary to remove this signal and its aliases from the analysis. Even if we cannot claim the planetary nature of these signals they have a FAP < 1%. Their inclusion would corrupt the test at their

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Sensitivity limits for HD 20794 in mass. HD 20794 b,c, and d are above our sensitivity limits. The dataset, thanks to the YARARA correction of HARPS and the extreme precision of the ESPRESSO dataset, shows a sensitivity limit in mass under 10 M⊕, with the possibility to detect Earth-like planets for up to ~20 d. We have a sensitivity down to 2.5 M⊕ for the inner edge of the conservative HZ and down to 4.1 M⊕ for the outer edge of the HZ. We highlight in green the conservative HZ and in yellow the optimistic HZ.</figcaption></figure>

respective periodicities. We can detect Earth or sub-Earth planets up to ~20 d orbital period. We have the sensitivity to detect potential candidates with a minimum mass superior to 10 M⊕ for all the periods. We considered the minimum mass where we did not obtain a detection for the inner and outer conservative edges of the HZ. For the inner edge of the HZ, we are sensitive to planets with $M_p$sini > 2.5 M⊕, while for the outer edge of the HZ, we are sensitive to $M_p$sini > 4.1 M⊕. We can exclude the presence of giant planets in the HD 20794 system, making its architecture different from our Solar system. HD 20794 has a metallicity $[Fe/H] = -0.42 \pm 0.02$. The frequency of gas giants in systems known for hosting a super-Earth is derived for metal-poor stars in Bryan & Lee (2024) as $P(GG|SE, [Fe/H] < 0.0) = 4.5^{+2.6}_{-1.9} \%$. Their work defines a gas giant as a planet with a mass of > 0.5 M$_{Jup}$, while a super-Earth is a planet with a mass of between 1 and 20 M⊕. Bonomo et al. (2023) search for the frequency of cold Jupiters (CJ) in the presence of short-period (SP) planets. They define as cold Jupiters planets with mass comprised between 0.3 and 13 M$_{Jup}$ and orbital distance comprised between 1 and 10 AU. They define short-period planets with masses comprised between 1 and 20 M⊕ and orbital period < 100 d. Bonomo et al. (2023) find $P(CJ|SP) = 9.3^{+7.7}_{-2.9} \%$ for a sample of 38 stars with transiting planets from Kepler (Borucki et al. 2010) and K2 (Howell et al. 2014). The absence of a cold Jupiter orbiting around HD 20794 is in line with the absence of a strong correlation between the population inner of super-Earths and the presence of an outer giant seen in statistical surveys. We repeated the test for considering the detection limits as a function of the amplitude. We considered the same period grid as before and a
grid of 300 amplitudes uniformly spaced in logarithm spanning between 0.01 and 1.5 m s⁻¹. We recovered a detection limit in the amplitude of ~30 cm s⁻¹ almost constant with the orbital period.

We have for HD 20794 a *Gaia* DR3 renormalised unit weight error (RUWE) metric = 1.97. A value of this metric larger than one could indicate the source is not single, or problematic for the astrometric solution. Furthermore, there is evidence for astrometric acceleration due to a marginally significant HIPPARCOS-*Gaia* proper motion anomaly (Brandt 2021; Kervella et al. 2022). We are sensitive to the presence of a companion with a true mass of 50 M⊕ in the range of 3–10 au. A 50 M⊕ planet at 3–10 au would correspond, for such a star as HD 20794, to a signal with an amplitude of 2.9–1.6 m s⁻¹ at an orbital period of ~2000 and ~13 000 d, respectively. Such a signal would be detectable as a Keplerian signal at periods shorter than our baseline and should generate a detectable linear trend at longer orbital periods. We do not see evidence of additional Keplerian signals with P > 2000. We considered a model with an acceleration term to model an eventual long-term trend present in the data and we found a result compatible with no acceleration. We conclude that an eventual long-term companion with the characteristics indicated should lie on a low-inclination orbit.

### 5.3. A candidate target for future atmospheric characterization

The atmospheric characterization of exoplanets is fundamental to studying the composition of planets, and to infer their habitability. The atmosphere of exoplanets is studied both from the ground and space taking advantage of the footprint the atmosphere of the planet impresses to the stellar spectra once the planet transits. The most common technique used for this investigation is transmission spectroscopy.

Future facilities such as the high-resolution spectrographs RISTRETTO (Lovis et al. 2022) at VLT and ANDES (Marconi et al. 2021, 2022; Palle et al. 2023) at ELT (Padovani & Cirasuolo 2023) suggest an alternative path for the characterization of the atmosphere of exoplanets. The large diameter of the mirror of VLT and ELT makes it possible to spatially resolve the angular distance between stars and outer planets. Disentangling the light coming from the star from the light coming from the planet will make it possible to analyze the spectroscopy of exoplanets' atmospheres directly through high-dispersion coronagraphy (HDC). RISTRETTO is expected to permit the characterization of the atmosphere of exoplanets at planet-star contrast of 10⁻⁷. The angular distance limit is given by 2λ/D, where D is the diameter of the mirror. For VLT this formula transforms in an angular resolution of ~37.5 mas for a wavelength of 760 nm (Lovis et al. 2022). We are interested in this wavelength because it is where an oxygen band resides. The angular distance of HD 20794 d from the star is ~222 mas considering as reference the semi-major axis, while it spans between ~124 mas and ~322 mas for periastron and apoaster. The angular distance for HD 20794 c and HD 20794 b are ~60 and 20 mas respectively considering a circular orbit. Following Otegi et al. (2020) we can derive a minimum radius for our planets from their minimum masses. We considered the two regimes described in this work, ρ < 3.3 and ρ > 3.3 g cm⁻³. In the rocky planet scenario we derived R<sub>b</sub> ~ 1.3 R⊕, R<sub>c</sub> ~ 1.4 R⊕, and R<sub>d</sub> ~ 1.7 R⊕. In the volatile-rich scenario, we derived R<sub>b</sub> ~ 1.1 R⊕, R<sub>c</sub> ~ 1.4 R⊕, and R<sub>d</sub> ~ 2.1 R⊕. To derive the planet-to-star contrast, we followed Lovis et al. (2017). The planet-to-star contrast ratio is ~1.8·10⁻⁸/1.3·10⁻⁸ for planet b, ~2.6·10⁻⁹/2.2·10⁻⁹ for planet c, and ~2.7·10⁻¹⁰/4.1·10⁻¹⁰ for planet d in rocky/volatile-rich case.

We are considering an Earth-like albedo. If we consider the condition for HD 20794 d at the periastron we obtain a planet-to-star contrast ratio of ~8.9·10⁻¹⁰/1.4·10⁻⁹ for the rocky/volatile-rich case. All the targets are not feasible for RISTRETTO due to limits in contrast and/or distance from the star. ANDES is supposed to span a range between 400 and 1800 nm. ELT has an angular resolution of ~13.7 mas for a wavelength of 1300 nm (Suárez Mascareño et al. 2023). This wavelength is important because it is where we can find another oxygen band. The maximum field of view of ANDES will be 100X100 mas. Due to the characteristics of the spectrograph HD 20794 b and HD 20794 c will be feasible, in terms of angular separation, for atmospheric characterization. HD 20794 d, even if it is interesting for its location in the HZ, resides at an angular distance from the star outside the field of view of ANDES, so the telescope should be pointed away from the guiding star. Furthermore, the low planet-to-star contrast ratio for HD 20794 d and HD 20794 c would make it necessary for long integration times to achieve a significant S/N ratio. HD 20794 b could be an interesting case of directly characterizing the atmosphere of a 2 M⊕ planet orbiting close to a G-type star.

HD 20794 is part of the Habitable World Observatory Exoplanet Exploration Program (ExEP) Precursor Science Stars (Mamajek & Stapelfeldt 2024), a list of amenable candidates for characterization with Habitable World Observatory (HWO). In the future, HWO will serve as a space mission with a 6-m telescope aimed at characterizing exoplanets in reflected light through direct imaging at optical/near-infrared wavelengths. The main goal of the mission will be to study the atmosphere of planets in the HZ. HWO is currently part of the Great Observatory Maturation Program (GOMAP) of NASA. The selection criteria for the target stars were based on the brightness and similarity to the Sun. HD 20794 d is a planet of particular interest for this mission because it is a light planet orbiting in the HZ of its hosting star. HWO requires targets with a separation from the star of more than 60–70 mas, while HD 20794 d is at more than 200 mas from the star, crossing the HZ.

The Large Interferometer For Exoplanets (LIFE) is a mission aimed at characterizing habitable planets through nulling interferometry at mid-infrared wavelengths (Quanz et al. 2022). LIFE will investigate the thermal emission of the planets and the goal of the mission is twofold: the discovery of new planets and an in-depth atmospheric characterization of the most promising targets. The mid-infrared wavelength regime is suitable for the detection of various biosignatures, such as ozone (O₃), methane (CH₄), nitrous oxide (N₂O), chloromethane (CH₃Cl), and phosphine (PH₃) (Schwieterman et al. 2018). These biosignatures are not necessarily produced only in biological processes but combined detection of more than one of those will improve the possibility of a biological origin. Again, HD 20794 d could be a preferential candidate considering its characteristics.

### 5.4. Stellar activity

HD 20794 is known in the literature for being a quiet star. A measurement for the stellar rotational period is missing in the previous works. The RVs do not show clear signatures of the activity. We look for additional info in the activity indicators derived from spectra not corrected with YARARA. We can find a magnetic cycle in FWHM, BIS, Contrast, and S-index with a period of ~3000 d. Further analysis of the FWHM time series with GP points toward a rotation period of 35.0±3.2 d for the star. When we analyzed the time series of activity indicators derived from spectra corrected by YARARA, we found in BIS a strong
peak at ~40 d. A fit with a GP found a rotation period of 38.8$^{+2.4}_{-2.6}$ d. The fact we can see a signal compatible with a period of ~40 d is of particular importance because this period is reported in Pepe et al. (2011) as the orbital period of a planetary signal. This work did not use a correction for the activity. Furthermore, we exploit a larger number of observations in a longer observation timespan. Activity signals are not stable over time and this can lead their signature to change in strength over time. We conclude that the possibility for the 40-d signal to be of planetary origin is unlikely, while the possibility this periodicity is related to the rotation period is the most plausible option. We refer to Appendix E for more detailed activity analysis.

## 5.5. Telemetry of the instrument

We investigated the variation in the temperature of echelle gratings of the spectrographs and how this impacts RVs and activity indicators (Suárez Mascareño et al. 2023). We looked for correlations with FWHM and RVs and we found a correlation between the temperature of the Echelle gratings and the FWHM. In HARPS, we found a periodicity of the temperature variation related to the year. In ESPRESSO, we see a linear trend in the variation of temperature. We have a Pearson linear correlation of 0.67, with a *p*-value of 2.48·10⁻⁹. We see a jump in the temperature of the echelle grating of ESPRESSO in correspondence with the intervention the instrument underwent in May 2022. If we consider measurements taken before and after the intervention as separate datasets and subtract the mean value, we can see a variation related to the year with a peak-to-peak amplitude of 25 mK. For HARPS, this variation is larger, with a peak-to-peak variation of more than 300 mK. We can see the same trends in the temperature of other detector elements, such as the dichroic, the main collimator, and the field mirror. We do not see a correlation between ESPRESSO RVs, FWHM, and BIS with the pressure of the blue and red cryostat. We considered a linear detrending between indicators and temperature for all indicators in the analysis. This detrending is fundamental to revealing the signatures of magnetic cycle and rotation in our dataset.

## 5.6. Intra-night scatter and possible origins of the dispersion of ESPRESSO dataset

During the observations of HD 20794 with ESPRESSO, a sequence of short exposures was taken within the same night. This included one night with one single exposure and an additional night with two exposures. For all the remaining nights, we had between 7 and 25 exposures per night. The sampling strategy we followed observing HD 20794 with ESPRESSO was dictated by the brightness of the target and the presence of phenomena such as granulation and oscillations, which have relatively short timescales (from minutes to several days) compared to the periodicities we are interested in the detection of exoplanets. The observing strategy, alongside the nightly binning of the observations, should permit us to mitigate the effect of *p*-mode oscillations. The median timespan between the first and the last exposure on the same night is 15 min. We observe a mean intranight dispersion of 49 cm s⁻¹, considering only the nights with more than five exposures. The maximum scatter we observe over a single night is 90 cm s⁻¹, while the minimum is 22 cm s⁻¹. This must be compared to an error on the single measurement with a mean of 10 cm s⁻¹. For a star such as HD 20794, pulsations should have a timescale of the order of minutes (Chaplin et al. 2019). The most common strategy to take out the contribution of the oscillation is to integrate the stellar spectra

for a time comparable to the pulsation timing of the star. In our case, this approach was impossible because a longer exposure time would have led to detector saturation. Instead of a longer exposure time, we exploit the presence of multiple exposures to mitigate the effect of pulsation (O'Toole et al. 2008; Dumusque et al. 2011). The mitigation of granulation effects on RVs can be obtained by observing a star multiple times per night with a span of some hours between observations. We do not have multiple observations per night for HD 20794 and we cannot even follow the granulation pattern with the multiple exposures we have on a single night, because of the short baseline (15 min) of our exposures. The granulation effect on RVs is difficult to predict, but it can induce a scatter similar to the scatter we observe in the ESPRESSO dataset once we remove the planetary signals (RMS of 72 cm s⁻¹ in the residuals). Most likely this scatter is of stellar origin, but in the time series, there is no evidence of a significant impact from the magnetic cycle of the star and rotation. Considering we are mitigating the effect of pulsations by binning, granulation remains the main effect we cannot model. Only a deep analysis of the behavior of granulation on RVs, through dedicated campaigns of consecutive exposures over a full night or longer, can help us understand the unexplained scatter we continue to see in the data.

## 6. Conclusion

The re-analysis of the system HD 20794, a bright G6V star (V = 4.34 mag) brought on the confirmation of the presence of the two planets with orbital periods of 18-d and 89-d, combining the RV measurements taken with HARPS and ESPRESSO, and thanks to the correction of the HARPS spectra with the YARARA-pipeline for extraction of velocities, which has permitted us to clean HARPS time-series from different systematic errors. These two planets were present in all the previous papers from the literature on the system and are found with more than 9σ significance.

We confirmed the presence of the 650-d planet which was reported first in Cretignier et al. (2023). This planet is of particular interest because it spends most of its orbit inside the HZ of the star. The closeness of the planetary system, summed with the distance of the star and the planet and the planet-to-star contrast ratio make this planet a good candidate for future atmospheric characterization through direct imaging facilities. We call the three planets HD 20794 b, HD 20794 c, and HD 20794 d for planets at 18-d, 89-d, and 650-d, respectively.

HD20794 b has a period of 18.3140$^{+0.0022}_{-0.0022}$ 0.003 d and a minimum mass of $M_p$ inini of 2.15 ± 0.17 $M_⊕$ with an eccentricity $e_b$ < 0.13. HD20794 c has an orbital period of 89.68 ± 0.10 d, a minimum mass of 2.98 ± 0.29 $M_⊕$, and an eccentricity $e_d$ < 0.16. HD 20794 d has an orbital period of 647.6$^{+2.5}_{-2.7}$ d, a minimum mass of 5.82 ± 0.57 $M_⊕$, and an eccentricity $e_d$ of 0.45$^{+0.10}_{-0.11}$. The amplitude related to the signals is by far under 1 m s⁻¹, 0.614 ± 0.048 m s⁻¹, 0.502$^{+0.048}_{-0.049}$ m s⁻¹, and 0.567$^{+0.067}_{-0.064}$ m s⁻¹ for planets b, c, and d respectively.

We did not find any evidence of signals related to the stellar rotation in RVs. For the activity indicators, we ran a parallel analysis with YARARA and not YARARA-corrected spectra. In the case of raw indicators, we found evidence of the presence of a magnetic cycle. We can see the cycle in the FWHM, S-index, BIS, and contrast. The timescale of the cycle is compatible with FWHM and S-index, pointing toward a cycle of ~3000 d. We can see some hints of the stellar rotation period in the analysis of FWHM with GP, where the rotation period results amount to
~40 d. In the YARARA corrected dataset, we find an excess of power at ~40 d in BIS that might be related to rotation, as we have seen from our analysis.

The low level of activity we find confirms the stability of HD 20794 and its suitability as a standard target for observations carried out by future facilities. The brightness of the target permits us to achieve a photon noise level for ESPRESSO of 10 cm s⁻¹.

*Acknowledgements.* N.N. acknowledges funding from Light Bridges for the Doctoral Thesis "Habitable Earth-like planets with ESPRESSO and NIRPS", in cooperation with the Instituto de Astrofísica de Canarias, and the use of Indefea- sible Computer Rights (ICR) being commissioned at the ASTRO POC project in the Island of Tenerife, Canary Islands (Spain). The ICR-ASTRONOMY used for his research was provided by Light Bridges in cooperation with Hewlett Packard Enterprise (HPE). J.I.G.H., A.K.S., R.R., C.A.P., N.N., V.M.P., and A.S.M. acknowledge financial support from the Spanish Ministry of Science and Innovation (MICINN) project PID2020-117493GB-I00. The project that gave rise to these results received the support of a fellowship from the "la Caixa" Foundation (ID 100010434). The fellowship code is LCF/BQ/D123/11990071. Co-funded by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This work was supported by FCT - Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 - Programa Operacional Competitividade e Internacionalização by these grants: UIDB/04434/2020 (DOI: 10.54499/UIDB/04434/2020); UIDP/04434/2020 (DOI: 10.54499/UIDP/04434/2020)". A.M.S acknowledges support from the Fundação para a Ciência e a Tecnologia (FCT) through the Fellowship 2020.05387.BD (DOI: 10.54499/2020.05387.BD). 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. This work was financed by Portuguese funds through FCT (Fundação para a Ciência e a Tecnologia) in the framework of the project 2022.04048.PTDC (Phi in the Sky, DOI 10.54499/2022.04048.PTDC). C.J.M. also acknowledges FCT and POCH/FSE (EC) support through Investigador FCT Contract 2021.01214.CEECIND/CP1658/CT0001 (DOI 10.54499/2021.01214.CEECIND/CP1658/CT0001). This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualization, exchange, and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch. X.D. acknowledges the support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement SCORE No 851555) and from the Swiss National Science Foundation under the grant SPECTRE (No 200021_215200). This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. F.P.E. and C.L.O. would like to acknowledge the Swiss National Science Foundation (SNSF) for supporting research with ESPRESSO and HARPS through the SNSF grants nr. 140649, 152721, 166227, 184618 and 215190. The ESPRESSO and HARPS Instrument Projects were partially funded through SNSF's FLARE, resp. FINES Programmes for large infrastructures. The INAF authors acknowledge financial support of the Italian Ministry of Education, University, and Research with PRIN 201278X4FL and the "Progetti Premiiali" funding scheme. This research has made extensive use of the SIMBAD database operated at CDS, Strasbourg, France, and NASA's Astrophysics Data System. This research has made use of NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. The manuscript was written using Overleaf. Extensive use of numpy (van der Walt et al. 2011) and scipy (Virtanen et al. 2020). The main analysis was performed in Python 3 (Van Rossum & Drake 2009) running on a Ubuntu system (Sobell 2015).

## References

* Aigrain, S., & Foreman-Mackey, D. 2023, ARA&A, 61, 329

* Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989

* Aller, A., Lillo-Box, J., Jones, D., Miranda, L. F., & Barceló Forteza, S. 2020, A&A, 635, A128

* Andretta, V., Doyle, J. G., & Byrne, P. B. 1997, A&A, 322, 266

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

* Baranne, A., Mayor, M., & Poncet, J. L. 1979, Vistas Astron., 23, 279

* Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373

* Basant, R., Dietrich, J., & Apai, D. 2022, AJ, 164, 12

* Biasiotti, L., Simonetti, P., Vladilo, G., et al. 2024, MNRAS, 530, 4300

* Bonomo, A. S., Dumusque, X., Massa, A., et al. 2023, A&A, 677, A33

* Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977

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

* Bryan, M. L., & Lee, E. J. 2024, ApJ, 968, L25

* Buchschacher, N., & Alesina, F. 2019, ASP Conf. Ser., 521, 757

* Buchschacher, N., Ségransan, D., Udry, S., & Díaz, R. 2015, ASP Conf. Ser., 495, 7

* Chaplin, W. J., Cegla, H. M., Watson, C. A., Davies, G. R., & Ball, W. H. 2019, AJ, 157, 163

* Cretignier, M., Dumusque, X., Hara, N. C., & Pepe, F. 2021a, A&A, 653, A43

* Cretignier, M., Francfort, J., Dumusque, X., Allart, R., & Pepe, F. 2021b, Astrophysics Source Code Library [record ascl:2102.022]

* Cretignier, M., Dumusque, X., & Pepe, F. 2022, A&A, 659, A68

* Cretignier, M., Dumusque, X., Aigrain, S., & Pepe, F. 2023, A&A, 678, A2

* Cretignier, M., Pietrow, A. G. M., & Aigrain, S. 2024, MNRAS, 527, 2940

* Damasso, M., Perger, M., Almenara, J. M., et al. 2022, A&A, 666, A187

* Delisle, J. B., Unger, N., Hara, N. C., & Ségransan, D. 2022, A&A, 659, A182

* Demangeon, O. D. S., Zapatero Osorio, M. R., Alibert, Y., et al. 2021, A&A, 653, A41

* Desort, M., Lagrange, A. M., Galland, F., Udry, S., & Mayor, M. 2007, A&A, 473, 983

* Díaz, R. F., Cincunegui, C., & Mauas, P. J. D. 2007, MNRAS, 378, 1007

* Dumusque, X. 2018, A&A, 620, A47

* Dumusque, X., Udry, S., Lovis, C., Santos, N. C., & Monteiro, M. J. P. F. G. 2011, A&A, 525, A140

* Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132

* Efron, B., Hastie, T., Johnstone, I., & Tibshirani, R. 2004, arXiv e-prints [arXiv:math/0406456]

* Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, Nature, 580, 597

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

* Fellgett, P. 1953, J. Opt. Soc. Am., 43, 271

* Feng, F., Tuomi, M., & Jones, H. R. A. 2017, A&A, 605, A103

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

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

* Gaia Collaboration 2020, VizieR Online Data Catalog: I/350

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

* Gomes da Silva, J., Santos, N. C., Bonfils, X., et al. 2011, A&A, 534, A30

* Hall, R. D., Thompson, S. J., Handley, W., & Queloz, D. 2018, MNRAS, 479, 2968

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

* Hara, N. C., Unger, N., Delisle, J.-B., Díaz, R. F., & Ségransan, D. 2022, A&A, 663, A14

* Hara, N. C., de Poyferré, T., Delisle, J.-B., & Hoffmann, M. 2024, Annal. Appl. Stat., 18, 749

* Hartman, J. D., & Bakos, G. Á. 2016, Astron. Comput., 17, 1

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

* Hég, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27

* Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398

* Jenkins, J. M., Chandrasekaran, H., McCauliff, S. D., et al. 2010, SPIE Conf. Ser., 7740, 77400D

* Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, SPIE Conf. Ser., 9913, 99133E

* Keenan, P. C., & McNeil, R. C. 1989, ApJS, 71, 245

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

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

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

* Kurucz, R. 1993, Robert Kurucz CD-ROM, 13

* Lo Curto, G., Pepe, F., Avila, G., et al. 2015, The Messenger, 162, 9

* Lovis, C., Dumusque, X., Santos, N. C., et al. 2011, arXiv e-prints [arXiv:1107.5325]

* Lovis, C., Snellen, I., Mouillet, D., et al. 2017, A&A, 599, A16

* Lovis, C., Blind, N., Chazelas, B., et al. 2022, SPIE Conf. Ser., 12184, 121841Q

* Mamajek, E., & Stapfeldt, K. 2024, AAS/Div. Extreme Sol. Syst. Abs., 56, 628,17

* Marconi, A., Abreu, M., Adibekyan, V., et al. 2021, The Messenger, 182, 27

* Marconi, A., Abreu, M., Adibekyan, V., et al. 2022, SPIE Conf. Ser., 12184, 1218424

* Mayor, M., Queloz, D., Marcy, G., et al. 1995, IAU Circ., 6251, 1

* Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20

* Meunier, N. 2021, arXiv e-prints [arXiv:2104.06072]
Nelson, B. E., Ford, E. B., Buchner, J., et al. 2020, AJ, 159, 73

Otegi, J. F., Bouchy, F., & Helled, R. 2020, A&A, 634, A43

O'Toole, S. J., Tinney, C. G., & Jones, H. R. A. 2008, MNRAS, 386, 516

Padovani, P., & Cirasuolo, M. 2023, Contemp. Phys., 64, 47

Palle, E., Biazzo, K., Bolmont, E., et al. 2023, arXiv e-prints
[arXiv:2311.17075]

Pepe, F., Lovis, C., Ségransan, D., et al. 2011, A&A, 534, A58

Pepe, F., Molaro, P., Cristiani, S., et al. 2014, Astron. Nachr., 335, 8

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

Quanz, S. P., Ottiger, M., Fontanet, E., et al. 2022, A&A, 664, A21

Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279

Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine
Learning (Cambridge: MIT Press)

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, SPIE Conf. Ser., 9143, 914320

Saar, S. H., & Donahue, R. A. 1997, ApJ, 485, 319

Santos, N. C., Bouchy, F., Mayor, M., et al. 2004, A&A, 426, L19

Schwieterman, E. W., Kiang, N. Y., Parenteau, M. N., et al. 2018, Astrobiology, 18, 663

Silva, A. M., Faria, J. P., Santos, N. C., et al. 2022, A&A, 663, A143

Sneden, C. 1973, ApJ, 184, 839

Sobell, M. G. 2015, A Practical Guide to Ubuntu Linux (London: Pearson Education)

Sousa, S. A. G., Santos, N. C., Monteiro, M. J. P. F. G., Israelian, G., & Mayor, M. 2007, in JENAM-2007, "Our Non-Stable Universe", 9-9

Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373

Sousa, S. G., Santos, N. C., Adibekyan, V., Delgado-Mena, E., & Israelian, G. 2015, A&A, 577, A67

Speagle, J. S. 2020, MNRAS, 493, 3132

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

Suárez Mascareño, A., Rebolo, R., & González Hernández, J. I. 2016, A&A, 595, A12

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

Suárez Mascareño, A., González-Álvarez, E., Zapatero Osorio, M. R., et al. 2023, A&A, 670, A5

Toner, C. G., & Gray, D. F. 1988, ApJ, 334, 1008

Torres, G., Andersen, J., & Giménez, A. 2010, A&A Rev., 18, 67

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22

Van Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: CreateSpace)

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Wildi, F., Pepe, F., Chazelas, B., Lo Curto, G., & Lovis, C. 2010, SPIE Conf. Ser., 7735, 77354X

Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577

<footnote n="1">¹ Light Bridges S.L., Observatorio del Teide, Carretera del Observatorio, s/n Guimar, 38500 Tenerife, Canarias, Spain</footnote>

<footnote n="2">² Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="3">³ Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="4">⁴ Consejo Superior de Investigaciones Científicas, Madrid, Spain</footnote>

<footnote n="5">⁵ Observatoire de Genève, Département d'Astronomie, Université de Genève, Chemin Pegasi 51b, 1290 Versoix, Switzerland</footnote>

<footnote n="6">⁶ Aix Marseille Univ, CNRS, CNES, Institut Origines, LAM, Marseille, France</footnote>

<footnote n="7">⁷ Department of Physics, University of Oxford, Oxford OX13RH, UK</footnote>

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

<footnote n="9">⁹ Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal</footnote>

<footnote n="10">¹⁰ Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="11">¹¹ Centro de Astrobiología (CAB), CSIC-INTA, ESAC campus, Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada (Madrid), Spain</footnote>

<footnote n="12">¹² INAF - Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, 34143 Trieste, Italy</footnote>

<footnote n="13">¹³ IFPU-Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy</footnote>

<footnote n="14">¹⁴ INAF - Osservatorio Astrofísico di Torino, Strada Osservatorio 20, 10025 Pino Torinese (TO), Italy</footnote>

<footnote n="15">¹⁵ ESO - European Southern Observatory, Av. Alonso de Cordova 3107, Vitacura, Santiago, Chile</footnote>

<footnote n="16">¹⁶ INAF - Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, 90134 Palermo, Italy</footnote>

<footnote n="17">¹⁷ Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências da Universidade de Lisboa, 1749-016 Lisboa, Portugal</footnote>

<footnote n="18">¹⁸ Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="19">¹⁹ Weltraumforschung und Planetologie, Physikalisches Institut, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="20">²⁰ ETH Zurich, Institute for Particle Physics & Astrophysics, Wolfgang-Pauli-Str. 27, 8093 Zurich, Switzerland</footnote>

<footnote n="21">²¹ ETH Zurich, Department of Earth Sciences, Sonneggstrasse 5, 8092 Zurich, Switzerland</footnote>

<footnote n="22">²² Subaru Telescope, National Astronomical Observatory of Japan, 650 N Aohoku Place, Hilo, HI 96720, USA</footnote>

<footnote n="23">²³ Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany</footnote>
## Appendix A: Methods

For the parameter estimation, we used the nested-sampling tool Dynasty (Speagle 2020). Dynasty gives an estimate of the natural logarithm of the evidence associated with a model, allowing for an easy model comparison. To reckon the different zero-points of different instruments, we always consider an offset term for each instrument in the analysis. A jitter term is added in quadrature to the nominal error of the different instruments for each time series. The jitter term takes into account all the sources of noise we are not modeling for, and the instrumental noise. To remove outliers from different datasets, we bin the observations nightly and we apply a cut on the dataset consisting of a 3 σ clipping joined with the exclusion of measurements where the error is larger than three times the median error for each dataset. To implement Gaussian processes (Sect. 4.1) in our analysis we used S+LEAF (Delisle et al. 2022). S+LEAF allows for only a certain kind of semi-separable matrix to be taken into account as a covariance matrix. In this way, the computational cost scales linearly with the dimension of the dataset, instead of with its cube, as it used to be in standard implementations. Before a fit for the jitter was available, exploratory periodograms of our datasets were generated adding in quadrature a white noise term to the error on the measurements equal to the standard deviation of the dataset.

## Appendix B: Planetary system

We show in Table B.1 the previous results obtained in the literature, compared with the results obtained in our work.

## Appendix C: L1 periodogram

A different way to assess the presence of significant signals in the dataset, instead of a GLS periodogram, is the L1 periodogram (Hara et al. 2020). GLS periodogram assesses the significance of each frequency of a grid of frequencies alone. L1 periodogram attempts to fit simultaneously for all the periodicities present in the dataset. We show in Fig. C.1 the result of the L1 periodogram applied to HD 20794. We are calculated the L1 periodogram with the LARS algorithm (Efron et al. 2004). We can see the most prominent peaks at 18.3 d, 89.5 d, and 646 d. These peaks correspond to the peaks we find in our blind-search analysis. Less significant peaks are found at ~ 85 d and 1026 d.

<figure><img src="image_17.png" /><figcaption>**Fig. C.1.** L1 periodogram of the full dataset of RVs for HD20794. We see the most significant peaks at 18.3 d, 89.6 d, and 646 d. These peaks correspond to the planets b, c, and d we are seeing in our analysis. We see other peaks at ~ 85 d and 1027 d. These peaks appear in our analysis but we cannot claim those as new detections.</figcaption></figure>

## Appendix D: 111-day signal

The periodogram of the residuals of the 3-Keplerian model shows a double peak in the residuals at 85.4 and 111.5 d. The blind search we conducted and the FIP analysis did not permit us to claim this new signal as a safe planet detection. Here, we discuss a dedicated analysis to investigate the characteristics of this signal. We tried a model with normal priors on the periods for the claimed planets and the 111-d signal. We considered both a circular and a Keplerian model. In the circular model, we obtained a solution for the 111-d signal with $K = 0.276_{-0.050}^{+0.049} \text{ m s}^{-1}$, and $P = 111.51_{-0.18}^{+0.19} \text{ d}$. This corresponds to a planet of mass $M_p \text{ sini} = 1.77 \pm 0.31 \text{ M} \oplus$ orbiting at a distance of $0.4192_{-0.18}^{+0.19} \text{ m}$. In the Keplerian model, we obtained a solution with a significant eccentricity $e = 0.40_{-0.12}^{+0.13} \text{ K}$, $K = 0.38 \pm 0.07 \text{ m s}^{-1}$, and $P = 111.52_{-0.16}^{+0.15} \text{ d}$. This corresponds to a planet mass of $M_p \text{ sini} = 2.23_{-0.34}^{+0.35} \text{ M} \oplus$ and an orbital semi-major axis of $0.4192 \pm 0.0018$. The presence of such an eccentric signal at a period close to the period of planet c would imply a crossing orbit for the two planets, making the system likely unstable. Furthermore, the amplitude of the hypothetical Keplerian signal should be detectable in our dataset considering our detection limit. The high level of eccentricity could be related to the tendency of the fitting algorithm to include all the remaining structures in the data inside the 111-d signal. We lack the sensitivity necessary for a proper characterization of this structure. Additional observations are needed to solve the planetary system of HD 20794 in the lower mass regimes.

## Appendix E: Stellar activity

The advancement in the development of ultra-stable spectrographs opened the possibility of achieving precision on the single measurement of RV well below $1 \text{ m s}^{-1}$. At this level of precision, the change in the stellar RV due to stellar physics-related phenomena becomes detectable so we need to consider it in our models because it can hide or mimic planetary signals in RV time series (Saar & Donahue 1997; Meunier 2021). Stars are affected by a plethora of non-equilibrium processes. We can find evolving oscillation patterns, granulation, rotation, flares, plages, and long-term magnetic cycles. These phenomena can affect the RV in timescales varying from minutes to years and with amplitudes ranging from cm s$^{-1}$ to tenths of m s$^{-1}$. In this section, we first outline the analysis of the stellar activity.

For oscillation and granulation, some strategies have been developed to mitigate their effect (Dumusque et al. 2011; Chaplin et al. 2019), including: (i) appropriate exposure times to average out the p-mode oscillations, (ii) multiple observations per night to mitigate granulation signal.

The stellar rotation period can differ significantly, depending on the stellar type and the stellar age (Suárez Mascareño et al. 2015). For solar-type stars, we do not have examples in the literature of rotation periods longer than 50 d. The most common tools in the field of exoplanets to deal with rotation related-stellar activity are Gaussian Processes (GP) (Rasmussen & Williams 2006; Haywood et al. 2014; Aigrain & Foreman-Mackey 2023). GPs offer a framework to represent signals that present an excess of power at certain frequencies, without being strictly periodic. Because of their flexibility, GP can adapt very well to the peculiar features of activity signals. Activity signals are not strictly periodic and are usually unstable over multiples of their characteristic timescale. GPs are flexible enough to capture
<table><caption>**Table B.1.** Planetary parameters from the literature works on HD 20794 and this work.</caption><thead><tr><th>Reference</th><th>K [ms<sup>-1</sup>]</th><th>P [d]</th><th>m<sub>p</sub>·sini [M<sub>⊕</sub>]</th><th>e</th></tr></thead><tbody><tr><td colspan="5"><strong>Pepe et al. (2011)</strong></td></tr><tr><td>18-d</td><td>0.83 ± 0.09</td><td>18.315 ± 0.008</td><td>2.7 ± 0.3</td><td>0</td></tr><tr><td>40-d</td><td>0.56 ± 0.10</td><td>40.114 ± 0.053</td><td>2.4 ± 0.4</td><td>0</td></tr><tr><td>89-d</td><td>0.85 ± 0.10</td><td>90.309 ± 0.184</td><td>4.8 ± 0.6</td><td>0</td></tr><tr><td colspan="5"><strong>Feng et al. (2017)</strong></td></tr><tr><td>18-d</td><td>0.81<sup>+0.00</sup><sub>-0.24</sub></td><td>18.33<sup>+0.01</sup><sub>-0.02</sub></td><td>2.82<sup>+0.10</sup><sub>-0.80</sub></td><td>0.27<sup>+0.04</sup><sub>-0.22</sub></td></tr><tr><td>89-d</td><td>0.60<sup>+0.00</sup><sub>-0.18</sub></td><td>88.90<sup>+0.37</sup><sub>-0.41</sub></td><td>3.52<sup>+0.48</sup><sub>-0.91</sub></td><td>0.25<sup>+0.16</sup><sub>-0.20</sub></td></tr><tr><td>147-d</td><td>0.69<sup>+0.15</sup><sub>-0.13</sub></td><td>147.02<sup>+1.43</sup><sub>-0.91</sub></td><td>4.77<sup>+0.96</sup><sub>-0.86</sub></td><td>0.29<sup>+0.14</sup><sub>-0.18</sub></td></tr><tr><td colspan="5"><strong>Cretignier et al. (2023)</strong></td></tr><tr><td>18-d</td><td>0.56<sup>+0.05</sup><sub>-0.05</sub></td><td>18.32<sup>+0.01</sup><sub>-0.01</sub></td><td>2.0<sup>+0.2</sup><sub>-0.2</sub></td><td>0.09<sup>+0.08</sup><sub>-0.06</sub></td></tr><tr><td>89-d</td><td>0.78<sup>+0.05</sup><sub>-0.05</sub></td><td>89.56<sup>+0.09</sup><sub>-0.10</sub></td><td>4.7<sup>+0.4</sup><sub>-0.4</sub></td><td>0.13<sup>+0.07</sup><sub>-0.07</sub></td></tr><tr><td>640-d</td><td>0.61<sup>+0.06</sup><sub>-0.06</sub></td><td>644.6<sup>+9.9</sup><sub>-7.7</sub></td><td>6.6<sup>+0.6</sup><sub>-0.7</sub></td><td>0.40<sup>+0.07</sup><sub>-0.07</sub></td></tr><tr><td colspan="5"><strong>This work</strong></td></tr><tr><td>18-d</td><td>0.614 ± 0.048</td><td>18.3142 ± 0.0022</td><td>2.15 ± 0.17</td><td>0.065<sup>+0.065</sup><sub>-0.046</sub></td></tr><tr><td>89-d</td><td>0.502<sup>+0.048</sup><sub>-0.049</sub></td><td>89.68 ± 0.10</td><td>2.98 ± 0.29</td><td>0.077<sup>+0.084</sup><sub>-0.055</sub></td></tr><tr><td>650-d</td><td>0.567<sup>+0.064</sup><sub>+0.067</sub></td><td>647.6<sup>+2.5</sup><sub>-2.7</sub></td><td>5.82 ± 0.57</td><td>0.45<sup>+0.10</sup><sub>-0.11</sub></td></tr></tbody></table>

**Notes.** We can see planets at 18 d and 89 d are recovered in all literature works on the system and by the present work. Planet at 40 d is recovered only in Pepe et al. (2011), the planet at 147 d is recovered only in Feng et al. (2017), and the planet at 640 d is recovered only in Cretignier et al. (2023) and this work.

this variability over time. We directly sample for the parameters of the covariance matrix between samples of the process under analysis, which correspond to our dataset. In this way, we can model the correlated noise in a Bayesian and mathematically defined framework. The largest drawback of GP is their tendency to overfit, which can be mitigated by deciding on appropriate kernels. Kernels describe the shape of the covariance matrix between observations. Different kernels have been proposed and used in the analysis of exoplanets. We tried different kernels in our analysis, which we briefly discuss here. The stochastic harmonic oscillator (SHO) is a kernel described in Foreman-Mackey et al. (2017). It requires three parameters: amplitude, rotation period, and a quality factor. The Matérn 3/2 exponential period (MEP) kernel and the exponential-sine periodic (ESP) kernel are approximations of the quasi-periodic kernel (Eq. (E.1)) as defined in Aigrain & Foreman-Mackey (2023):

$$k(\Delta t) = A \exp \left( - \frac{(\Delta t)^2}{2\rho^2} - \frac{2 \sin^2 \left( \frac{\pi \Delta t}{P} \right)}{\eta^2} \right), \quad (E.1)$$

where we can relate the parameters to physical quantities: the amplitude, *A*, the rotation period, *P*, the timescale of evolution, *ρ*, and the harmonic complexity, *η*.

Stars can also show in the RV time series signatures of magnetic cycles, which usually have periods of thousands of days for G-type stars (Lovis et al. 2011). The magnetic long-term cycle can also be seen in many activity indicators (Lovis et al. 2011; Gomes da Silva et al. 2011; Suárez Mascareño et al. 2016). Usually, this kind of noise is modeled with a sinusoid and a few harmonics at rational fractions of the cycle period. Over the years many diagnostic methods have been developed as proxies for activity, and nowadays the stellar activity analysis is a central part of the analysis performed in almost every work on extrasolar planetary systems. Activity indicators are sensitive to the change in the spectral profile stellar activity can give rise to (FWHM, bisector, contrast), or to the behavior of spectral lines,

which are good proxies for chromospheric activity ($S_{MW}$, H-$\alpha$, and Na doublet).

The FWHM is one of the most reliable indicators of stellar activity. The WHM measures the width of the CCF profile. A temporal variation of the broadening of the spectral lines is due to the change in the activity pattern of the star (Desort et al. 2007; Dumusque et al. 2014). This activity indicator has been used in many analyses as a proxy for the rotational period of the star, or alongside RV to better constrain the RV model (Suárez Mascareño et al. 2020; Faria et al. 2022; Suárez Mascareño et al. 2023).

The bisector span metric (BIS) measures the asymmetry in spectral lines. The variation in BIS is related to a change in the stellar profile due to the stellar activity happening on the star. The presence of spots on the stellar surface can create some unbalances in the flux from the star with a periodicity related to the rotational period of this. This results in a modification of the shape of the spectral lines (Queloz et al. 2001; Toner & Gray 1988). The pipelines of HARPS and ESPRESSO automatically compute this indicator.

The contrast of the CCF measures the ratio between the core of the CCF and the continuum. This indicator is sensitive to the magnetic level of the star, which can change due to phenomena as rotation and long-term cycles.

The H-$\alpha$ index is an indicator of the chromospheric activity of a star. It measures the strength of the core of the line of the H-$\alpha$ emission, or absorption for weakly active stars. This line is mostly effective for M-dwarfs, even if it can also be used for G-type stars. On G-type stars the indicator is usually sensitive to the contamination of filaments and telluric lines. The calculation of this indicator is based on the comparison between the flux in the H$\alpha$ line and the continuum, as explained in Gomes da Silva et al. (2011). The exact method to extract this indicator can slightly change from work to work. In our analysis, we average the flux in a window with a width of 1.6 Å centered at the center of the line, at 6562.808 Å (H$\alpha_{core}$) and divide it by the sum of the averaged flux in two bands, centered at 6550.87 Å (H$\alpha_V$) and 6580.31 Å
(HαR), with a width of 10.75 Å and 8.75 Å respectively. The equation for the extraction of the indicator is

$$H\alpha = \frac{H\alpha_{core}}{H\alpha_{R} + H\alpha_{V}}. \qquad (E.2)$$

The core of the emission of the H & K lines of calcium is an important indicator of chromospheric activity. We used the same approach of Lovis et al. (2011) to derive the Mount Wilson S-index (S<sub>MW</sub>). The calculation is done by defining two triangular passbands centered in the core of the emission lines, respectively at 3968.469 Å for the H line and at 3933.663 Å for the K line, with an FWHM of 1.09 Å. We normalized this for the continuum using two rectangular bandpasses of 20 Å centered at V = 3901.070 Å and R = 4001.070 Å respectively. The equation we used to derive the indicator is

$$S_{MW} = \alpha \frac{\bar{N}_{H} + \bar{N}_{K}}{\bar{N}_{R} + \bar{N}_{V}} + \beta, \qquad (E.3)$$

where Ñ<sub>H</sub>, Ñ<sub>K</sub>, Ñ<sub>R</sub>, Ñ<sub>V</sub> are the mean fluxes in each band and α and β are constants for the calibration, namely, α = 1.111 and β = 0.153. For a description of the effect of different phenomena such as spots or plages on the profile of Ca H & K lines, see Cretignier et al. (2024).

Na I D resonance lines (D1: 5895.92 Å; D2: 5889.95 Å) are visible in all kinds of stars as an important absorption feature of the spectra. In cool stars as late G, K, and M dwarfs the Sodium absorption lines are usually present as strong absorption wings and have been demonstrated as a useful tool for studying stellar atmospheres (Andretta et al. 1997; Díaz et al. 2007). A method to compute this activity indicator is described in (Díaz et al. 2007).

In Fig. E.11 we show the time series of the activity indicators together with their GLS periodograms (Zechmeister & Kürster 2009). For observations taken before BJD 2453500, we see a higher dispersion in the H03 dataset than for observations made after that date. For this reason, we discard observations made previous to this date. We can see a linear trend present in FWHM and contrast in the H03 dataset. This effect is due to a change in the focus of the instrument through the years and can be modeled with a linear trend for the H03 dataset alone. We desire to model stellar activity on itself. For this reason, we first analyze the activity indicators derived without taking into account the YARARA correction. YARARA corrects spectra also for stellar activity, hiding its signature in the indicators. The trade-off is the presence in our time series of the systematics YARARA corrects for. The analysis of activity indicators corrected for YARARA will be reported in Sec. E.8 and YARARA activity indicators and their GLS periodograms are shown in Fig. E.12.

## E.1. Temperature of the echelle gratings

Modern spectrographs need a stabilized environment in pressure and temperature to achieve the necessary precision in RV measurements. Even if the environmental conditions are strictly controlled and monitored, some noise is still present. This variability in principle can affect RVs and activity indicators time series (Suárez Mascareño et al. 2023). We controlled the variations of pressure and temperatures for different components of HARPS and ESPRESSO. We checked for correlations with activity indicators and RVs. We found a strong linear correlation (Pearson correlation = 0.67, p-value = 2.48·10<sup>-9</sup>) between FWHM and the temperature of echelle gratings for ESPRESSO. We saw a similar correlation for BIS. Also, the HARPS dataset

shows similar trends for BIS and FWHM compared to the temperature of the echelle gratings. In the following analysis, we corrected the correlation with temperature for all the activity indicators with linear trends. Figure E.1 shows the time series of echelle gratings temperature for HARPS and ESPRESSO together with their GLS periodogram. Figure E.2 shows the correlation between FWHM and variation in temperature of echelle gratings.

## E.2. FWHM

We analyzed the FWHM as it is extracted by the ESPRESSO pipeline (Pepe et al. 2021).

The first step in the analysis of FWHM is done by detrending the time series for the change of focus HARPS suffered before fiber intervention. Secondly, we modeled the correlation with the change in temperature of echelle gratings. Once we linearly detrend the FWHM for focus and temperature, we can see a peak in the GLS periodogram of the residuals at 3093.89 d as shown in Fig. E.3.

This periodicity is typical for the magnetic cycle of a G-type star. We modeled the cycle with a sinusoidal and a few of its harmonics. We found the best solution to be a model with a sinusoidal and its fourth harmonic, with a period of 1/4 of the main term period. For the cycle, we found a period of 3005<sup>+25</sup><sub>-20</sub> d. In Fig. E.4 we show a plot of the cycle model with the GLS periodogram of the residuals. The plotted time series is already detrended for temperature changes and focus, while these parameters were fitted simultaneously to the long-term cycle. We searched for the rotation period of the star with GP. We tried different GP kernels to model the suspected rotation signal. We found as the best kernel the MEP kernel. We found a rotation period of 35.0<sup>+3.2</sup><sub>-2.5</sub> d.

In the GP model, we found a period for the cycle of 3020<sup>+111</sup><sub>-50</sub> d. The period of the cycle is compatible with the one we found in the cycle-only model, even if with a larger uncertainty. The amplitude of the cycle remains different from zero at 6.6 σ. We show the posterior distribution of the amplitude of the cycle, the period of the cycle, and the rotation period in Fig. E.5

## E.3. BIS

For the analysis of BIS, we followed a procedure similar to that of the analysis of the FWHM. We considered a model with a linear detrending of the BIS against the temperature of the echelle gratings. Once we detrend the indicator for temperature of echelle gratings we can still see strong peaks at 2530 d and 850 d. These two peaks suggest the presence of a magnetic cycle for HD 20794. The ratio between the two periods is very close to 1/3, pointing toward the presence of a cycle with its harmonic at P/3. We modeled the cycle with a sinusoid with a norm prior for the period centered at 3000 d, with a width of 250 d, and a second sinusoid, with free amplitude and phase but with the period constrained to be 1/3 of the period of the main component. We found a period for the cycle of 2521<sup>+40</sup><sub>-34</sub> d. Figure E.6 shows the model for BIS over-imposed to the dataset once we have subtracted the correlation with the temperature of echelle gratings.

Once we subtract the model for BIS we only see peaks with FAP > 1 %. The most prominent peak in the residuals is at 36.79 d. This could be a signature of the rotation period of the star. We tried a fit with a GP but the result did not converge to a defined solution. In Fig. E.7 we show the posterior distribution for the amplitude and the period of the BIS cycle.
<figure><img src="image_18.png" /><figcaption>**Fig. E.1.** Variation in the temperature of echelle gratings for HARPS (upper left panel) and ESPRESSO (lower left panel) with their respective GLS periodograms. We can see a strong peak in the GLS periodogram of the HARPS dataset close to 1 year, related to the seasonal variation of the instrument.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. E.2.** Correlation between variation of the temperature of the echelle grating and the FWHM.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. E.3.** GLS periodogram of the FWHM time series after linear detrending for the temperature of echelle gratings and the change of focus in the H03 dataset. We can see the presence of a long-term magnetic cycle at more than 3000 d.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. E.4.** Model of the cycle in FWHM for HD 20794 (top). The time series is detrended for correlation with temperature and focus. GLS periodogram of the residuals once the model is subtracted (bottom).</figcaption></figure>

### E.4. S-index

From the dataset and the GLS periodogram of the residuals of Fig. E.11, we see a quadratic trend in the S-index time series. We tried to model this indicator with a quadratic trend together with a cycle, its first harmonic, a sinusoidal to model 1-year variability, and the first harmonic of that year. For the cycle period, we considered a normal prior centered at 3000 d with a width of 500 d. The period of the cycle derived for the S-index is $P = 2884^{+86}_{-90}$ d. Once we subtracted the model for the cycle, we could see some additional structure in the residuals, resulting in a peak at 527 d, but no peaks in correspondence with the orbital periods of the planets we have found in our analysis. We can see some
<figure><img src="image_22.png" /></figure>

**Fig. E.5.** Posterior distributions for parameters of the magnetic cycle and rotation period we can observe in HD 20794 for FWHM. Panel (a): Posterior distribution of the amplitude of the cycle found in FWHM. Panel (b): Posterior distribution of the period of the cycle found in FWHM. Panel (c): Posterior distribution of the rotation period found in FWHM.

excess in power at lower periods, so we tried a blind search for
the rotation period through a GP analysis. The GP analysis did
not converge to a satisfying solution, probably due to the higher
level of noise in this indicator compared to the FWHM. In Fig.
E.8 we show the model of the cycle for S-index together with the
GLS periodogram of the residuals. In Fig. E.9 we show the pos-
terior distribution for the amplitude and the period of the S-index
cycle.

E.5. Contrast

In the contrast, we can see a structure similar to the one we find in
FWHM. We see a strong linear trend in the H03 dataset, related
to the change in focus of the instrument through the years. We
can see some trends in the E19 dataset too, related to the change
in temperature of the echelle gratings. We started the analysis
by detrending the different time series with a linear trend to the
temperature of echelle gratings, with an additional term for the
H03 dataset to correct the change in focus. Once we subtracted
this model we could see peaks in the GLS periodogram of the
residuals related to the 1-year variation and probably related to
the magnetic cycle of the star. We saw peaks at ~ 1800 d and ~

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

**Fig. E.6.** Model of the cycle in BIS for HD 20794 (top). The time series is detrended for correlation with temperature. GLS periodogram of the residuals once the model is subtracted (bottom).

3200 d. We modeled this signal with a sinusoid and its first and
second harmonic and a sinusoid with normal prior on the period
centered around 1-year. We found a solution with a period for
the cycle equal to 4065<sup>+32</sup><sub>-30</sub> d. Once we subtracted this signal we
still found some structure in the residuals but we did not find any
periodicity related to the period of the detected planet or stellar
rotation.

E.6. Hα

In Hα, we can see the presence of a signal related to the year
with the first three harmonics, respectively at 1/2,1/3, and 1/4
of the year-signal. We found the best model a model with three
sinusoids, one with a period of 1 year and the other two at half
and one-third of this period. When we subtracted this model, we
did not see any peak in the GLS periodogram in correspondence
with any interesting signal in RVs, nor a region of the parameter
space suitable for a signature of stellar rotation. The modeling
with a GP did not converge to a result on the rotation period.

E.7. Na I D

In Na I D, we only find a peak in the GLS periodogram at 367
d and a second peak at the first harmonic of the year signal. We
tried to model the time series with a sinusoidal and a harmonic
at half the period. Once we subtracted this model, we can see
a quadratic trend in the full dataset. We modeled it and we did
not find any significant signal in the GLS periodogram of the
residuals. We did not go further in the analysis of this indicator.

E.8. YARARA activity indicators

In our analysis of activity, we did not only consider the activ-
ity indicators extracted from the spectra not corrected from
YARARA. We also tried to analyze the HARPS dataset alone
for activity indicators extracted from spectra corrected for
YARARA. Here we report the analysis of BIS. BIS is the only
indicator showing a prominent signal in the GLS periodogram
that could be related to the rotation of the star, while most of the
other indicators just show a signal related to the 1-year variation
or its harmonics.
<figure><img src="image_24.png" /><figcaption>**Fig. E.7.** Posterior distributions for parameters of the magnetic cycle we can observe in HD 20794 for BIS. Panel (a): Posterior distribution of the amplitude of the cycle found in BIS. Panel (b): Posterior distribution of the amplitude of the harmonic at P/3 of the period of the cycle found in BIS. Panel (c): Posterior distribution of the period of the cycle found in BIS.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. E.8.** Model of the cycle in S-index for HD 20794 (top). The time series is detrended for correlation with temperature. GLS periodogram of the residuals once the model is subtracted (bottom).</figcaption></figure>

We can see a prominent peak in BIS at 40.53 d. This period is of particular interest because it is similar to the period of the

<figure><img src="image_26.png" /><figcaption>**Fig. E.9.** Posterior distributions for parameters of the magnetic cycle we can observe in HD 20794 for S-index. Panel (a): Posterior distribution of the amplitude of the cycle. Panel (b): Posterior distribution of the amplitude of the harmonic at P/2. Panel (c): Posterior distribution of the period of the cycle.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. E.10.** Model for BIS dataset in case of YARARA correction. We can retrieve clues on the rotation period of the star around 40 d. The amplitude of the signal is very low as expected from the low level of activity we can observe for HD 20794 and the fact spectra were previously corrected with YARARA</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>**Fig. E.11.** Activity indicators time series derived before the YARARA correction. Panel (a): FWHM time series. It is possible to see trends in the different datasets. For H03 we have a trend with the change of focus of the instrument, while for H15 and E19 we have a linear correlation with the temperature of the echelle gratings. Panel (b): GLS periodogram for FWHM. Panel (c): BIS time series. It is possible to see a linear trend in E19 and H15. Panel (d): GLS periodogram for BIS. Panel (e): Time series of Contrast. Panel (f): GLS periodogram of Contrast. Panel (g): Hα time series. Panel (h): GLS periodogram of Hα. Panel (i): S-index time series. Panel (j): GLS periodogram of S-index. Panel (k): Na I D time series. Panel (l): GLS periodogram of Na I D. The red/green/blue horizontal dashed lines indicate FIP levels of 10/1/0.1 % respectively.</figcaption></figure>

40-d planet (40.114 ± 0.053) detected in Pepe et al. (2011). This planet was not confirmed in the following analyses on the target. This points toward a stellar origin for this signal. This periodicity could be related to the rotation period of the star. To test this hypothesis, we tried to model this indicator with a GP, using a MEP kernel. In our analysis, we found a non-converging result. We noticed the mean error on the measurement BIS provided by

YARARA is larger than the RMS of the measurement. Assuming the error on the BIS is overestimated we considered as error on the BIS the error we had for RV. With this assumption, we found a convergent solution for the rotation period at 38.8 $^{+2.4}_{-2.6}$ d, with a low amplitude, as expected for a quiet star such as HD 20794. We overplotted the model on the HARPS dataset in Fig. E.10.
<figure><img src="image_29.png" /><figcaption>**Fig. E.12.** Activity indicators time series derived after the YARARA correction. Panel (a): FWHM time series. Panel (b): GLS periodogram for the FWHM, a signature of the year is present. Panel (c): BIS time series. Panel (d): GLS periodogram for BIS, peak at 40.53 d. Panel (e): Time series of Contrast. Panel (f): GLS periodogram of the Contrast, 1-year signature is visible. Panel (g): $H\alpha$ time series. Panel (h): GLS periodogram of $H\alpha$, a peak corresponding to the revolution of the Earth is visible. Panel (i): S-index time series. Panel (j): GLS periodogram of S-index. Panel (k): Na I D time series. Panel (l): GLS periodogram of Na I D. The red/green/blue horizontal dashed lines indicate FIP levels of 10/1/0.1 % respectively.</figcaption></figure>