© The Authors 2024

1 Introduction

The field of exoplanet science has been evolving very quickly in recent years towards the detection and characterisation of Earth-like exoplanets thanks to the combined effort of space missions such as Kepler (Borucki et al. 2009), TESS (Ricker et al. 2015), and CHEOPS (Benz et al. 2021), and ground-based highresolution ultrastable spectrographs, such as HARPS (Mayor et al. 2003), HARPS-N (Cosentino et al. 2012), CARMENES (Quirrenbach et al. 2016), and ESPRESSO (Pepe et al. 2021). In particular, the exoplanet community is already finding potentially habitable Earth-like planets (e.g. Gillon et al. 2017; Dittmann et al. 2017; Zechmeister et al. 2019; Lillo-Box et al. 2020; Suárez Mascareño et al. 2023; Cadieux et al. 2024), paving the path towards the detection of an Earth twin, the ultimate goal of the ESPRESSO project in the long term, and other projects such as the Terra Hunting Experiment with the upcoming HARPS3 spectrograph (Thompson et al. 2016). The discovery of the temperate Earth-mass planet Proxima Centauri b (Anglada-Escudé et al. 2016) orbiting the closest star to our Sun has propelled exoplanet studies to focus the search on Earth-like planets in the habitable zone around stars of the solar neighbourhood. These temperate Earth-like planets will be the main targets of future facilities in the next decade, such as ANDES (Marconi et al. 2022) at the Extremely Large Telescope (ELT) in Cerro Armazones within the European Southern Observatory (ESO), with the goal of studying their atmospheres to search for biomarkers using both transmission and reflected light spectroscopy (Palle et al. 2023).

During the last decade, the blind radial velocity (RV) search for these Earth-mass exoplanets quickly shifted to the continuous monitoring of M dwarfs, with the development of new instruments in the near-infrared, such as CARMENES (Ribas et al. 2023), SPIRou (Donati et al. 2020), and NIRPS (Bouchy et al. 2017). M dwarfs are the most common stellar type in the Galaxy, representing about 80% of the stars in the solar neighbourhood (Reylé et al. 2021). M dwarfs are cooler, intrinsically less luminous, and less massive than Sun-like stars, and have habitable zones closer to their host star, making them ideal targets for blind RV searches of Earth-like planets. The ESPRESSO spectrograph at the Very Large Telescope (VLT, ESO) has had a significant impact in exoplanet science since it began regular operations at Paranal Observatory in October 2018 (Pepe et al. 2021; González Hernández et al. 2018). ESPRESSO has demonstrated unprecedented capabilities, aiming at 10 cm s⁻¹ RV precision. ESPRESSO has confirmed, for instance, the Earth-mass planet Proxima b and discovered the sub-Earth Proxima d, a 0.26 M_⊕ planet (approximately twice the mass of Mars) orbiting Proxima Centauri, from the measurement of a small RV semi-amplitude of 39 ± 7 cm s⁻¹ (Suárez Mascareño et al. 2020; Faria et al. 2022). ESPRESSO is opening a new frontier at sub-m/s precision, making it possible to discover and characterise Earth- and sub-Earth-mass and sub-Earth size exoplanets in the solar neighbourhood. ESPRESSO has detected, for instance, the 0.4 M_⊕ planet L98-59 b (half of the mass of Venus) orbiting an M3V star (Demangeon et al. 2021), and one super-Earth and two super-Mercuries HD 23472 d,e,f with masses of 0.54–0.76 M_⊕ orbiting a K4V star (Barros et al. 2022).

Barnard’s star (GJ 699) is the second closest stellar system to our Sun, after the α Centauri stellar system, and has been investigated in great detail since its discovery (Barnard 1916). It is the nearest single star to our Sun, the closest M dwarf after Proxima Cen, at a distance of about 1.8 parsecs, and is the star with the highest proper motion (Gaia Collaboration 2021), causing significant Doppler shifts due to secular acceleration (Kürster et al. 2003). France et al. (2020) measured the X-ray flux of GJ699 with the Chandra satellite in the energy range 0.3-10keV at F_x ~ 4.8 × 10⁻¹⁴ erg cm⁻² s⁻¹ (log₁₀(L_X[erg s⁻¹]) = 25.3; L_X/L_bol = 1.6 × 10⁻⁶). The X-ray luminosity is within a factor of two of previous ROSAT data (log₁₀(L_X[erg s⁻¹]) = 25.6). This low X-ray luminosity with log₁₀(L_X/L_bol) ~ −5.8 (France et al. 2020) indicates a low level of current magnetic activity (Stelzer et al. 2013). Previous spectroscopic works revealed a low level of chromospheric activity with (Suárez Mascareño et al. 2015; Astudillo-Defru et al. 2017a; Toledo-Padrón et al. 2019), suggesting a slow rotation period of P_ROT ~ 140 d, which was estimated using Equation (1) in Suárez Mascareño et al. (2018). This value is in agreement with the photometric value of P_ROT ~ 130 d derived from HST photometry (Benedict et al. 1998). Toledo-Padrón et al. (2019) reported a rotation period of P_ROT = 145 ± 15 d from the time-series analysis of spectroscopic activity indexes, and also found evidence of a long-term activity cycle of Barnard’s star from a time series of the CaHK index and ASAS-SN m_V photometry with a periodicity of P_CYC ~ 3225–3850 d. Reiners et al. (2022) estimated a relatively low surface average magnetic field strength at ⟨B⟩ ~ 0.43 ± 0.08 kG from spectral line fitting of the Zeeman broadening covering a wide range of different Landé-g values, consistent with the star’s low magnetic activity level, whereas Cristofari et al. (2023) found a lower value of ⟨B⟩ ~ 0.21 ± 0.08 kG. Donati et al. (2023) measured longitudinal magnetic fields using SPIRou data and investigated its temporal variations to infer a rotation period of P_ROT = 136 ± 16 d in GJ 699, in agreement with previous estimates.

Ribas et al. (2018) reported the discovery of a 3.3 M_⊕ super-Earth-like planet candidate orbiting Barnard’s star with an orbital period of 233 d. This result was challenged by Lubin et al. (2021), who argues that the signal is transitory in nature and is connected to stellar activity, as the planet candidate period is close to 1 yr alias of the rotation period (Toledo-Padrón et al. 2019). More recently, Artigau et al. (2022) used SPIRou data to show that a model including a 233d planetary signal with a RV semi-amplitude of K_p = 1.2 m s⁻¹ is disfavoured when compared to a flat model.

Here we present the ESPRESSO observations of Barnard’s star (GJ 699), showing a sub-m/s precision that reveals the presence of a short-period sub-Earth planet and three additional sub-Earth planet candidates. ESPRESSO data allow us to also evaluate the presence of the super-Earth planet candidate reported in Ribas et al. (2018).

2 Observations

The ESPRESSO consortium is the collaborative effort of Switzerland, Italy, Portugal, and Spain, with ESO as an associate partner, to develop, build, and scientifically exploit the ESPRESSO¹ instrument (Pepe et al. 2021). The ESPRESSO project has mainly been dedicated to the search for and characterisation of exoplanets (e.g. Lillo-Box et al. 2021; Faria et al. 2022; Suárez Mascareño et al. 2023; Lavie et al. 2023; Castro-González et al. 2023; Suárez Mascareño et al. 2024) and exoplanet atmospheres (e.g. Ehrenreich et al. 2020; Borsa et al. 2021; Azevedo Silva et al. 2022), and the measurement of fundamental constants of the Universe (e.g. Martins et al. 2022; Murphy et al. 2022).

Barnard’s star (GJ 699) is a (main) target of the guaranteed time observations (GTOs) of the ESPRESSO instrument. It was monitored over four years from May 2019 to July 2023. The main goal of the ESPRESSO GTO has been to search for rocky planets in the habitable zone (HZ) of nearby stars. Barnard’s star is considered a primary target due to its proximity to our Sun, its relatively low magnetic activity, and the possibility to search for Earth-like planets within its HZ.

This work has also made use of public HARPS, HARPS-N, and CARMENES data, with some of the HARPS and HARPS-N spectra taken by consortium members as part of the follow-up of Barnard’s star, as we describe below.

2.1 ESPRESSO

The Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO; Pepe et al. 2021) is a fibre-fed, cross-dispersed, high-resolution échelle spectrograph located in the Combined Coudé Laboratory (CCL) at the incoherent focus, where the Front-End unit can combine the light from up to four Unit Telescopes (UT) of the Very Large Telescope (VLT) at Paranal Observatory (ESO, Chile). The so-called Coudé train optical system feeds the light of each UT to the spectrograph. The Front-End corrects the light beam for atmospheric dispersion with the atmospheric dispersion corrector (ADC) and ensures the light is centered in the fibre with two independent pupil and field stabilisation units. The light from the target and the sky enter the instrument simultaneously through two separate fibres. ESPRESSO, unlike any other ESO instrument, is able to operate simultaneously with either one UT or several of the four 8.2-m UTs. The light of one or several UTs is fed through the Front-End unit into optical fibres that scramble the light within the Fiber-Link unit and provide excellent illumination stability to the spectrograph, using octogonal (1-UT) or square (4-UT) fibres. The instrument, aiming at a long-term 10 cm s⁻¹ RV stability, is temperature-controlled and pressure-stabilised within a vacuum vessel (VV). The reference fibre fed simultaneously with stabilised Fabry-Pérot unit allows the tracking of instrument drifts down to the cm s⁻¹ level. In the most used singleHR (1-UT) mode, a fibre of 140 μm core, equivalent to 1″ on the sky, provides a resolving power of R ~ 138 000 in the wavelength range 378–789 nm, sampling properly the resolution element with 4.5 pixels in two different detector binning setups HR11 and HR21.

We obtained 157 ESPRESSO observations of Barnard’s star from May 2019 to July 2023. Nine of them were taken before the intervention done at the end of June 2019 to upgrade the fibre link, which increased the photon-detection efficiency reaching more than 10% at seeing better than 0.75″ (Pepe et al. 2021). This intervention introduced an RV offset, leading us to consider two separate E18 and E19 datasets at about BJD[d] = 2458 660. In March 2020, operations at Paranal were interrupted due to the COVID-19 pandemic and ESPRESSO was taken out of operations from March 2020 to December 2020, leading to a large gap after the first year of the ESPRESSO observing campaign. However, following our analysis of the ESPRESSO data of different RV standard stars, we conclude that the change in one of the calibration lamps after the ramp-up of the instrument at the end of 2020 does not justify any additional RV offset (see e.g. Figueira et al., in prep.). With the ESPRESSO pipeline² version 3.0.0, the wavelength calibration and chromatic drift account for the change of lamp that originally created the need to separate the E19 and E21 (see e.g. Faria et al. 2022). Thus, now there is no offset or appreciable difference, and we do not need to separate these datasets. There has been another significant intervention in the instrument in May 2022 at about BJD[d] = 2459720 to repair the blue cryostat. In fact, both the blue and the red cryostat were changed, but the analysis of RV standards does not justify the need to split the E19 dataset into two (see Section 5).

The wavelength calibration done by the DRS uses Th-Ar lamp combined with Fabry-Pérot (FP) etalon exposures. Due to the observed brightness of Barnard’s star (m_V = 9.5, see Table 1), the ESPRESSO observations were carried out with the FP as simultaneous calibration in the reference fibre B using the HR11 binning (R ~ 138 000) and a typical exposure time of 900 s, with four spectra taken with 1200 s (two in E18 and two in E19) and two spectra taken with 550 and 426 s. The ESPRESSO data covers a time baseline of 1532.7 d (4.2 yr) from BJD[d] = 2458 606.79918 (May 2019) to 2460 139.51282 (July 2023).

Fig. 1 shows three different RV computations of the ESPRESSO data reduced using the ESPRESSO pipeline version 3.0.0. In the top panel, we see the RVs provided by the ESPRESSO Data Reduction Software (DRS; Di Marcantonio et al. 2018), which uses the cross-correlation technique; in the middle panel, the RVs computed using the S-BART (Silva et al. 2022) code³, a semi-Bayesian radial velocity computation through template matching (TM); and in the bottom panel, the RVs extracted using the line-by-line (LBL) technique applied to these ESPRESSO observations (Artigau et al. 2022). DRS, TM and LBL RVs are not significantly different from each other. We also depict the generalized Lomb-Scargle (GLS; Zechmeister & Kürster 2009) periodograms of the three sets of RVs where we only see slightly different power of some peaks, mostly related to stellar activity. The root-mean-square (RMS) of the RVs are very similar, 1.81, 1.86 and 1.78 m s⁻¹ for DRS, TM and LBL, respectively, whereas the mean/median uncertainty of the RVs are 16.5/14.4, 11.0/9.6 and 10.3/9.2 cm s⁻¹ for DRS, TM and LBL, respectively. Only eight E18 RVs (instead of nine) are shown since the LBL code crashed for one spectrum with BJD[d] = 2458 642.75641. We also remove one spectrum with BJD[d] = 2 459 867.56166 that has bad quality RV, FWHM measurements, clearly off by about 10 σ from the median values.

We note that the ESPRESSO pipeline version 3.0.0 does not include yet the telluric correction from Allart et al. (2022). Thus DRS RVs have not been computed after telluric correction, which may explain the larger DRS uncertainties compared to TM and LBL techniques. S-BART masks the tellurics at a 1% threshold, which is a quite conservative mask, thus discarding a considerable amount of RV content. S-BART first constructs a synthetic spectra with the TelFit code (Gullikson et al. 2014), using the weather conditions of the observing block with the highest relative humidity. It determines the continuum level through a median filter and finds where the spectra is more than 1% away from the continuum, thus masking only deeper features. Then it flags those regions as tellurics accounting for barycentric Earth radial velocity (BERV) changes before continuing with the RV computation (Silva et al. 2022). The LBL code also discards features affected by telluric contamination. Given the small differences between the TM and LBL uncertainties, we decided to adopt the S-BART TM RVs as our preferred RV measurements for the rest of the paper. This amounts to nine E18 and 147 E19 data points. Since only a few points are taken the same night, we decided to bin RV and FWHM time-series are subsequently with a 1 d step, and after, two E19 points are discarded with RV uncertainties larger than 50 cm s⁻¹, the final sample of ESPRESSO data includes nine E18 and 140 E19 points (see Fig. 2). Full width at half maximum (FWHM) measurements from cross-correlation functions (CCFs) were automatically provided by the ESPRESSO DRS (see e.g. Fig. 3). The RMS of E18 and E19 RVs computed with S-BART are 0.96 and 1.87 m s⁻¹. The RMS of the FWHM of E18 and E19 are 2.46 and 3.34 m s⁻¹. The statistics of the uncertainties of RV and FWHM measurements are summarised in Table 2.

We also applied the telluric correction using the code from Allart et al. (2022) only to the 156 useful ESPRESSO spectra. We recomputed a new set of RVs (labelled as TMtc RVs) using the template matching S-BART code, but this time masking out those regions at a 60% threshold, to avoid including deeper features affected by tellurics that may have not been properly modelled. We use this set of TMtc RVs later in this work (see Section 6.6) but we keep TM RVs as our main ESPRESSO RV dataset, which we compare with ESPRESSO DRS and LBL RVs, and that we also use together with the HARPS, HARPS-N and CARMENES datasets described in Sections 2.2 and 2.3. The TMtc RVs are very similar to the TM RVs, with a minor improvement, showing an RMS of 0.97 and 1.83 m s⁻¹ in TMtc RVs compared with an RMS of 0.96 and 1.87 m s⁻¹ in TM RVs. The mean/median uncertainties of TMtc RVs are 0.10/0.08 and 0.10/0.09 m s⁻¹, thus very similar to those of TM RVs (see Table 2). This final ESPRESSO dataset of CCF FWHM and TM RV measurements is available at the CDS portal, together with the FWHM and RV HARPS, HARPS-N and CARMENES datasets described in Sections 2.2 and 2.3.

Fig. 1 ESPRESSO RV measurements (right) and GLS periodograms (left) of GJ 699 after subtracting the median of each dataset before (E18) and after (E19) the intervention in June 2019. Also shown are RV measurements from the ESPRESSO Data Reduction Software (DRS; top), from the S-BART template matching (TM) code (middle), and from the line-by-line LBL code (bottom).

Fig. 2 Radial-velocity measurements (right) and GLS periodograms (left) of GJ 699 of ESPRESSO (top), HARPS and HARPS-N (middle), and CARMENES (bottom).

Fig. 3 ESPRESSO FWHM measurements (top), RV measurements (middle), and RV residuals (bottom) from the SHO (PROT and PROT/2) GP model, and GLS periodograms (left) of GJ 699. The uncertainties include the jitter term coming from the global model A in Table 3.

2.2 HARPS and HARPS-N

The High Accuracy Radial Velocity Planet Searcher (HARPS; Mayor et al. 2003) is a fibre-fed échelle high-resolution (R ~ 115 000) spectrograph installed in 2003 at the 3.6m telescope in La Silla Observatory (ESO, Chile). It covers the wavelength range 378–691 nm, and it is contained in a vacuum vessel to minimise the temperature and pressure variations that may cause spectral drifts. HARPS spectra used in these work can be downloaded from the ESO archive⁴ from different ESO programs⁵, and cover a time baseline of 2095.2 d (5.7 yr) from BJD[d] = 2457934.65911 (June 2017) to 2460029.83365 (March 2023). All the HARPS data used in this work were taken after the fibre-link upgrade in 2015, and thus we label these data as H15. The exposure time varies from 600 to 1800 s with a typical exposure of 900 s. Wavelengths are calibrated using a Th-Ar lamp combined with FP etalon exposures (Wildi et al. 2010). Spectra taken in 2017 were taken without any reference calibration in fibre B, and from 2018 with FP simultaneous reference in fibre B. We use a total of 114 useful HARPS spectra, which after 1D binning are turned into 105 HARPS data points.

The High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPS-N; Cosentino et al. 2012) is a fibre-fed échelle high-resolution (R ~ 115000) spectrograph installed in 2012 at the 3.6m Telescopio Nazionale Galileo (TNG) in the Observatorio del Roque de los Muchachos (ORM, La Palma, Spain). Having very similar instrument specifications as HARPS, it can also reach an RV stability better than 1 m s⁻¹ (Pepe et al. 2014), and covers the wavelength range 383–693 nm. HARPS-N spectra used in this work cover a time baseline of 1496.9 d (4.1 yr) from BJD[d] = 2457 626.41888 (August 2017) to 2459 123.34337 (September 2020). HARPS-N spectra can be accessed at the TNG archive⁶ from different Spanish CAT programs⁷. As for HARPS data, the wavelength calibration is done using a Th-Ar lamp combined with FP etalon exposures, with the science spectra taken with FP simultaneous reference in fibre B. We label these data as HAN. The total number of HARPS-N spectra is 133 and after binning using 1 d step we finally have 58 HARPS-N data points. The exposure time was 900 s before 2020 and the data taking during the COVID-19 pandemic was taken with three spectra per night of 1200 s each, except for one night with three spectra of 1800 s. The HARPS-N data taken during the COVID-19 pandemic intended to cover the gap of ESPRESSO observations during the ESO Paranal Observatory shutdown in 2020.

Both H15 and HAN spectra were reduced with the standard DRS pipelines at both instruments and the RVs were extracted using the S-BART code (see Fig. 2). The CCF FWHM measurements of HARPS and HARPS-N spectra were computed by adding a colour-correction factor order by order, following Suárez Mascareño et al. (2023). After discarding RV and FWHM measurements with uncertainties larger than 0.85 and 2.5 ms⁻¹, respectively, we are left with 92 data points in H15 and 56 data points in HAN. The RMS of H15 and HAN RVs computed with S-BART are 1.97 and 1.90 m s⁻¹. The RMS of H15 and HAN FWHM measurements are 5.8 and 3.4 m s⁻¹ (see Table 2).

2.3 CARMENES

The Calar Alto high-Resolution search for M-dwarfs with Exoearths with Near-infrared and optical Échelle Spectrographs (CARMENES; Quirrenbach et al. 2016) are visual (VIS) and near-infrared (NIR) vacuum-stabilised spectrographs covering 520–960 nm and 960–1710 nm with a spectral resolution of 94600 and 80400, respectively. They are located at the 3.5 m telescope of the Centro Astronómico Hispano en Andalucía (CAHA) at Observatorio de Calar Alto (Almería, Spain). The wavelength calibration is performed by combining hollow cathode (U-Ar, U-Ne, and Th-Ne) and Fabry-Pérot etalon exposures. The instrument drift during the nights is tracked with the FP in the simultaneous calibration fibre. We downloaded the CARMENES public data (Ribas et al. 2023) of Barnard’s star, which is the CARMENES VIS data of the RV survey within the GTO programme (CARMENES Data Release 1)⁸. CARMENES RV measurements were obtained using the template-matching SERVAL algorithm (Zechmeister et al. 2018). We use the RVs corrected for nightly zero point (NZP) offsets. The CARMENES spectrograph is usually wavelength calibrated each afternoon and nightly instrumental drifts are measured with the FP etalon, but stellar RVs from the same night often share common systematic effects, producing NZP offsets generally of a few ms⁻¹ with a median error bar of 0.9 ms⁻¹ (Ribas et al. 2023). The useful 520 CARMENES spectra corrected for NZP used in this work cover a time baseline of 1751.5 d (4.8 yr) from BJD[d] = 2457422.74662 (February 2016) to 2459 174.25596 (November 2020). After binning using 1 d step we end up with 501 CARMENES data points. The FWHM measurements of CARMENES CCF profiles were obtained by Lafarga et al. (2020) using the RACCOON code, and provided in Ribas et al. (2023). We discarded those RV and FWHM measurements with uncertainties larger than 2.5 m s⁻¹ and 26 m s⁻¹, respectively, to remove only a few points and to slightly clean the dataset. The final number of CARMENES data points, labelled as CAR, is 479. The RMS of the RV measurements is 2.18 m s⁻¹ (see Fig. 2). The RMS of FWHM measurements is 10.9 m s⁻¹ (see Table 2).

3 Stellar properties

Barnard’s star (GJ 699) is a bright (m_V = 9.5) very nearby M 3.5-M 4 dwarf star located at 1.8 parsecs from the Sun (Gaia Collaboration 2021). The main stellar properties are provided in Table 1. We adopted the weighted mean mass estimated (M_* = 0.162 ± 0.007 M_⊙) from the three mass determinations (based on the mass-radius relation, the spectroscopic log g and 2MASS K_s photometry) in Schweitzer et al. (2019). We checked that the updated parallax in Gaia Collaboration (2021) does not change the values given in Schweitzer et al. (2019). We used the 156 ESPRESSO spectra of Barnard’s star to create a master mean spectrum (see Fig. A.1). We used this master ESPRESSO spectrum to derive the stellar parameters (T_eff, log g and [Fe/H]) and the total line broadening velocity, υ_br, using the STEPARSYN code⁹ described in Tabernero et al. (2022). The derived stellar parameters, given in Table 1, are compatible with those used in Schweitzer et al. (2019) which are the spectroscopic parameters (T_eff = 3273 ± 51 K, log g = 5.11 ± 0.07, [Fe/H] = −0.15 ± 0.16) determined using CARMENES VIS data (Passegger et al. 2018). These are also compatible with those derived in Marfil et al. (2021) using both CARMENES VIS and NIR data (T_eff = 3254 ± 32 K, log g = 5.13 ± 0.12, [Fe/H] = −0.57 ± 0.10) and Jahandar et al. (2023) using SPIRou data (T_eff = 3231 ± 21 K, log g = 5.08 ± 0.15, [Fe/H] = −0.39 ± 0.03).

The luminosity (L_* [10⁻³ L_⊙] = 3.558 ± 0.072) of GJ 699 from Schweitzer et al. (2019) and the spectroscopic effective temperature (T_eff [K] = 3195 ± 28) derived from the ESPRESSO master spectrum of GJ 699 was used to estimate the habitable zone (HZ)¹⁰. We find an inner boundary 0.049 AU (recent Venus) and an outer boundary 0.129 AU (early Mars) (Kopparapu et al. 2014) of the HZ, corresponding to orbital periods of 9.84 and 41.88 d, respectively. The inner edge of the HZ for worlds with very little water content (with 1% relative humidity and albedo A = 0.2) could extend inwards to 0.036 AU (Zsom et al. 2013), or 0.026 AU in the case of high albedo (A = 0.8), which corresponds to orbital periods of 6.25 and 3.85 days, respectively.

4 Stellar activity

Stellar activity is possibly the main source of false positive planetary detections from RV time series. Activity signals and their aliases, although not necessarily persistent in time, caused by the presence of long-lived large spots (or spot groups) on the stellar surface can often create periodic signals that can easily mimic planetary signals (Queloz et al. 2001; Robertson et al. 2014; Suárez Mascareño et al. 2015, 2017). However, in many cases, it is possible to track and study the behaviour of stellar activity with time series of activity indexes and the changing shapes of computed cross-correlation functions, simultaneous to the RV time series.

4.1 GP model

To evaluate the behaviour of stellar activity, we model the time series of each activity indicator using the Gaussian process framework (GP; e.g. Rasmussen & Williams 2006). The GP framework is commonly used in the analysis of stellar activity in RV times series (e.g. Haywood et al. 2014; Faria et al. 2016). The stellar noise is described with a covariance function dependent on a set of parameters, some of them typically associated with a physical quantity. GP models can be used for instance to model the activity signal without requiring a detailed knowledge of the distribution, temperature contrast and lifetime of active regions on the stellar surface. GP models are flexible to model quasi-periodic signals, accounting for changes in the amplitude, phase, or even small period changes. This flexibility can however easily over-fit the data, sometimes suppressing possible planetary signals. Recently, there have been efforts to better constrain the GP models using the variability of stellar activity indicators, such as training with photometric data or activity indicators (e.g. Haywood et al. 2014) or simultaneous modeling of activity proxies and radial velocity measurements with shared hyper parameters (e.g. Suárez Mascareño et al. 2020; Faria et al. 2022). A more sophisticated approach is the use of multi-dimensional GPs, which join the fit of all time series under a single covariance matrix (Rajpaul et al. 2015; Barragán et al. 2022; Delisle et al. 2022). This implementation assumes that there is an underlying function governing the behaviour of stellar activity, G(t), which manifests itself in each time series as a linear combination of itself and its gradient, G′(t), and their amplitudes for each time series j (Δ TS_j), with j = 0, ..., N for N times series, following the FF′ formalism (Aigrain et al. 2012), as described in Equation (1)). (1)

Suárez Mascareño et al. (2020) showed a very good correlation between FWHM of the CCF and the activity-induced RV in the analysis of Proxima using ESPRESSO data. This offers a compelling new approach to study stellar activity signals and separate them from planetary signals in M-dwarfs.

Following Suárez Mascareño et al. (2023), we use the S+LEAF code (Delisle et al. 2022), which extends the formalism of semi-separable matrices introduced with Celerite (Foreman-Mackey et al. 2017) to allow fast evaluation of GP models even in the case of large datasets. The S+LEAF code allows the simultaneous fit of a GP to several time series, based on a linear combination of the GP and its derivative, with different amplitudes for each time series (see Equation (1)). The S+LEAF code supports a wide variety of GP kernels with fairly different properties. After testing several kernel functions, based on the shape of posterior sample distributions, we chose a combination of two simple harmonic oscillators (SHO) at the first and second harmonics of the rotation period, P_ROT and P_ROT/2. The selected kernel is defined as: (2)

with τ = t_n − t_n−1, representing the time-lag between measurements.

Following Equation (1), the activity induced signal in every specific time series j is: (3)

where G_SHO,i and is the realisation of a GP with kernel k_SHO,i and its first derivative. Following Foreman-Mackey et al. (2017), the k_SHO,i kernel is defined as: (4)

with η = (1 − (2L/P_i)⁻²)^1/2, controlling the damping of the oscillator. This kernel has a power spectrum density: (5)

where ω is the angular frequency, ω_i is the undamped angular frequency for each component (ω_i = 2 π / P_i), S_i is the power at ω = ω_i, and Q_i is the quality factor. The parameters S_i, P_i and Q_i are sampled in the covariance matrix, related to the amplitude (C_i), rotation period (P = P_ROT) and timescale of evolution (L = T_ROT) as shown in Eq. (6). (6)

The covariance matrix also includes a term of uncorrelated noise (σ), independent for every instrument. This term is added in quadrature to the diagonal of the covariance matrix to account for all unmodelled noise components, such as uncorrected activity or instrumental instabilities.

The amplitudes C_i in Equations (4) and (6) are related to the amplitude of the underlying function, not to any of the specific time series. We chose to adopt S_i = 1, thus fixing their power at ω = 0. Thus, the amplitudes of every component will be governed by the parameters A_ih shown in Euqtion (3).

We model the data using Bayesian inference via nested sampling (Skilling 2004), which in turn allows efficient exploration of large parameter spaces as well as obtaining Bayesian evidence from the model (i.e. marginal likelihood, ln ). We used the Dynesty code (Speagle 2020), which employs multi-ellipsoidal decomposition (Feroz et al. 2009) to more efficiently sample large prior volumes. We use the default configuration, which uses a random walk or random cut sampling strategy (Handley et al. 2015a,b) depending on the number of free parameters. We set the number of live points equal to 100 · N_par, and the number of slices equal to 2 · N_par, with N_par, the number of free parameters of the global model, including GP parameters.

Fig. 4 Analysis of the FWHM of the ESPRESSO CCF. Panels a and b: FWHM time-series with the best-model fit. The data is split into two panels because of a large period with no observations between the two campaigns. The shaded area shows the variance of the GP model. Panel c: GLS periodogram of the CCF FWHM data. The red vertical dashed line shows the most significant period. Panel d: relationship between the CCF RV and CCF FWHM data. The best fit is shown when the slope is ≥3σ different from zero. Panels e and f: residuals of the CCF FWHM after subtracting the best model fit. Panel g: GLS periodogram of the residuals. Panel h: comparison of the CCF RV and gradient of the CCF FWHM model.

4.2 Activity indicators

Following Lovis et al. (2011); Suárez Mascareño et al. (2015, 2018); Toledo-Padrón et al. (2019), we measure from the ESPRESSO spectra the following activity indexes: S -index or S _MW, defined similar to the original Mount Wilson index, measured from the line core fluxes of Ca II H&K stellar lines relative to the continuum fluxes, and the H_α and Na I-indexes, from the stellar lines H_α and Na I doublet, all sensitive to chromospheric activity. As stellar line shape varies with magnetic activity, we also built time series of quantities extracted from the shape of the cross-correlation function: the full width at half maximum, the bisector span (BIS) and the CCF contrast. All these measurements are simultaneous to the RV measurements as they are extracted from the same spectra.

We model each of the ESPRESSO time series individually with the adopted GP formalism and the results are shown in Fig. 4 for the FWHM of the ESPRESSO CCF, and Figs. C.1 and C.2 for the other activity indicators. The timescale of evolution of activity signals typically spans between one and two rotations, sometimes longer for M-dwarfs (Giles et al. 2017). We leave the rotation period, P_ROT, and the timescale, T_ROT, free in a wide range, with the priors (2,1000) d and (4, 4000) d, respectively, using a log-scale to allow for a long tail towards long timescales in the persistence of signals. We use log-normal priors for the amplitudes and jitter terms, centered on ln(RMS) of the data and with a sigma of ln(RMS) of the data. When using a GP with a completely free amplitude and jitter parameters, on data that includes multiple signals, there is a large risk that the GP absorbs all variations present in the data. Constraining the parameters in this way ensures a smooth GP model, preventing it from over-fitting variations at short timescales without fully excluding any region of the parameter space.

The GP analysis on the time series of the ESPRESSO FWHM measurements provides a and a timescale . A similar result was found for the bisector span ( and ) and the H_α-index ( and ). These values are consistent with previous P_ROT = 145 ± 15 d (Toledo-Padrón et al. 2019), mostly based on the time series of H_α-index measurements with a 15-yr baseline derived from seven different high-resolution spectrographs. In all cases the timescale is shorter but consistent with the rotation period. The CCF contrast ( and ), S _MW ( and ) and the NaI-index ( and ), show longer periods and larger uncertainties, thus marginally consistent. The period measured in the GP analysis of the different activity indicators is close to the period shown in the GLS of the different time series in Figs. C.1 and C.2. The structure of peaks in the GLS shows some complexity possibly related to the differential rotation, and with half a rotation, and sometimes the 1-yr alias at 240 d of the rotation period at 145 d. The highest peak of the GLS of BIS and CCF contrast falls at about the rotation period, whereas the Hα-index peaks at about half the rotation period. S _MW, Na I-index and FWHM show the highest peak at about 240 d.

We see clear correlations of the RV measurements with negative slope for H_α index, and positive slopes for FHWM and BIS measurements, in Figs. 4, C.1 and C.2. All fits performed to measure the slopes have both horizontal and vertical uncertainties taken into account. Clearly, the cleanest model with minimum residuals is provided by the FWHM, which shows a positive correlation with the CCF RVs. We therefore choose the FWHM in the joint analysis together with RVs to try to search for planetary signals in the RV time series while modeling simultaneously the stellar activity using both the FWHM and RV times series (see e.g. Fig. 3). The difference between the correlation RV vs. FWHM seen in this work and that of Suárez Mascareño et al. (2020) could be related to the specific nature of the active regions. The signature in RV of spot-induced variations causes a correlation between the δ/δt FWHM (or δ/δt Flux) and the RV. Variations caused by plages, however, cause a correlation between the FWHM (or Flux) and RV (formulas 11 and 12 in Aigrain et al. 2012, and Figure 3 in Dumusque et al. 2014). In this second case, the dominant bulk of the change in RV is due to inhibition of convective blueshift. In the case of Proxima, the rotation signal is very clear in photometric time-series, and not so much in chromospheric indicators, indicating spot-dominated variations. In the case of Barnard’s star, the rotation signal is easy to detect in chromospheric indicators (such as H alpha), but not in photometry (Toledo-Padrón et al. 2019). This, combined with the correlation between FWHM and RV, hints at inhibition of convective blueshift being the cause of the RV variations.

4.3 Photometry

We used the ASAS-SN Sky Patrol online tool¹¹ to inspect what photometry was available taken with the same time baseline that the ESPRESSO and other spectroscopic measurements shown in Fig. 2. We were able to acquire 3163 data points that cover the period January 2015 to January 2024 (approx. 2457000-2 460 300 HJD) and that contain measurements of two different pass-bands: V (N = 722) and g (N = 2441).

The proper motion of GJ 699 exceeds 10 arcsec/yr, while the ASAS-SN detector resolution is 8 arcsec/pixel. This means that the centroid of the star will drift by at least one pixel each year. We account for this in the following manner. Firstly, we took the interval 2457000-2460300 HJD which approximately reflects the ASAS-SN coverage of GJ 699. Then, we did split it in times-tamps at every 100 d, including the interval endpoints. For each timestamp, we: (i) compute the expected astrometric position of GJ 699; (ii) use this position to produce the light curve through the ASAS-SN Aperture Photometry Pipeline; (iii) clip this light curve to the region defined by a range of ± 50 d within the times-tamp. This query resulted in 34 individual light curves. Then, we concatenated all individual light curves into a common one, and computed the BJD from HJD values using the online tool¹² provided by Eastman et al. (2010). The g magnitude measurements were collected using several different cameras (see Fig. 5). We discarded the group of those g magnitude measurements higher than 10.5 which were far from the median, leaving a final set of N = 2215. We binned the data using a 1 d step. We also discarded the g magnitudes from bl (N = 161) and bH (N = 97) cameras which show a slope versus BJD not following the rest of the g measurements, leaving a final set of 1 d binned V (N = 268) from the bd camera (we discarded the only two V_bh points) and g_bt (N = 201), g_bp (N = 104) and g_bD (N = 205) magnitudes (see Fig. 5).

We use the ASAS-SN photometry to verify the stellar activity behaviour we see in the RV and FWHM measurements (see Section 6.1). We model the photometry with two double sinusoidal models as in Equation (7): (7)

one to account for the long-term cycle and another to model the rotation modulation, where ω₂ = 2ω₁ = 2π f / P, with f = 1/t and T_0,i=1,2 = t_mid + P · ϕ_i, with t_mid, the mid-time of the observation baseline and ϕ_i, the phase of the sinusoidal function. We left A₁, A₂, ϕ₁, ϕ₂ and the period P (equal to either the cycle period, P_CYC, or the rotation period, P_ROT) as free parameters, together with offsets and jitter terms to each of the magnitudes in the likelihood function. In Fig. 5, we display the ASAS-SN photometry versus BJD with fitted model. We also show GLS periodograms before and after subtracting the model, with the long-term signal and the rotation signal and its 1-yr alias detected. The posterior distributions point to a long-term cycle of P_CYC ~ 3440 d and a rotation signal of P_ROT ~ 144 d, when assuming Gaussian priors centered on the expected cycle and rotation periods based on the results obtained in the activity analysis presented in Section 6.1 and supported by previous independent results (Toledo-Padrón et al. 2019). The adopted time baseline of the photometry analysed here does not allow to get a better result with wide priors on the long-term cycle. We refer to Toledo-Padrón et al. (2019) for a deeper analysis of photometric data of Barnard’s star with a longer baseline of 15 years.

Fig. 5 ASAS-SN photometry in the same temporal baseline as the RV data together with the double sinusoidal models of the long-term cycle and rotation signals (top left) and the GLS periodogram before and after subtracting the long-term cycle (top right). The posterior and prior distributions of these models including the long-term cycle (bottom left) and rotation periods (bottom right).

5 Telemetry data

We evaluate the temperature and pressure within the vacuum vessel (VV) to track the stability of the instrument during the ESPRESSO observations of Barnard’s star. In Fig. 6, we display the échelle grating temperature, T_ech, and the vacuum vessel temperature as a function of BJD and barycentric radial velocity correction. The échelle temperature at a median temperature of 18,695 K, shows jumps of −51, −16 and +48 mK for the E18, E19 and E22 datasets. Here, this E22 refers to the jump in échelle temperature happening at about BJD[d] = 2459720 in May 2022. We check that this behaviour is shared by other temperature sensors within the VV, where, after the intervention in May 2022, we see jumps in temperature in the échelle, in the red cross disperser (RCD) and in the blue cross disperser (BCD) of approximately +50 mK, +80 mK and +120 mK, respectively. This corresponds to another significant intervention in the instrument where both the blue and the red cryostat were changed. However, the analysis of RV standards does not justify the need to split the E19 dataset into two datasets, but we separate with E22 in Fig. 6 for the sake of clarity. We show in the upper-right panel the VV pressure versus BJD keeping reasonably constant at a level of about 0.6 × 10⁻⁶ mbar and specially sensitive to the interventions coming down to stability after these interventions, in particular just after May 2022. We also display the median values of T_ech in the upper-left panel as horizontal dashed lines and subtract them to display the ΔT_ech vs. BJD in the next bottom panel. We perform a sinusoidal fit using the highest peak at 180 d (half a year) of the red GLS periodogram in the right panel to illustrate the seasonal or year dependence. The black GLS periodogram shows how this yearly dependence would disappear from ΔT_ech vs. BJD. In the next panel down, we show again more clearly this yearly dependence versus BERV, and again we perform a second-order polynomial function. We show the corresponding GLS periodogram of ΔT_ech vs. BJD in the right panel again with the result before (red line) and after (black line) subtracting this fit. We do not see any significant peak at periods shorter than 10 days.

However, this dependence of the ΔT_ech, which varies smoothly with BERV within the range of [−10,+15] mK, with a 4.9 mK RMS around this fit, does not apparently affect the RV and the FHWM as seen in the bottom four panels of Fig. 6. We do not see any trend of the RV and FWHM versus BERV or ΔT_ech. As discussed later in Section 6.1, the variations found in FHWM and RV measurements are closely related to stellar activity, and the seasonal effect that we see in ΔT_ech versus BERV, is either too small or removed with the drift corrections provided by the simultaneous FP calibration in fibre B during each observing exposure with the science target in fibre A. The drift corrections are very stable over the whole set of 4.2 yr of observations from May 2019 to July 2023, with a mean blue detector drift of −0.30 m s⁻¹ with a RMS of 0.84 m s⁻¹, and a mean red detector drift of −0.08 m s⁻¹ with a RMS of 0.33 m s⁻¹. The behaviour of both detectors seems better after the intervention in May 2022, with a mean drift and RMS of −0.26 m s⁻¹ (0.22 m s⁻¹) and −0.08 m s⁻¹ (0.22 m s⁻¹), for the blue and red detector, respectively.

Fig. 6 ESPRESSO RV and FWMH measurements, and échelle temperature sensor and vacuum vessel pressure sensor measurements versus BJD and barycentric RV correction. The GLS periodograms (left panels) show the result before (red line) and after (black line) subtracting the fitted function (dashed red line) in the left panels.

6 Analysis

Using the ESPRESSO RV and FWHM data, we build a procedure to search and confirm candidate planets orbiting Barnard’s star.

6.1 Stellar activity

The magnetic activity of the star is typically the dominant source of RV variations seen in most M dwarf stars, even the most apparently quiet stars such as Barnard’s star. The RMS of the ESPRESSO TM RV is 1.86 m s⁻¹, and the peak-to-peak RV variation goes from −4 to +4 m s⁻¹ (see Fig. 1). Following Suárez Mascareño et al. (2020); Faria et al. (2022); Suárez Mascareño et al. (2023) we run a simultaneous model using the ESPRESSO FWHM and RV data only including the GP described in Section 4.1, together with offsets and jitter parameters for the two ESPRESSO datasets E18 and E19. We note that the priors adopted for the RV jitter parameters were chosen to be higher than the expected solution to avoid over-fitting in the different model runs, thus forcing the runs if required to converge to lower jitters in the posteriors. In any case, the posterior distributions were significantly narrower than the prior distributions and typically within the prior range. Fig. 3 depicts the result of this first activity-only model (model A in Table 3). The total number of parameters includes four offsets, four jitters, and ten GP parameters including rotation, timescale and eight GP amplitudes. This includes both FWHM and RV parameters. The GP uses the two simple harmonic oscillators (SHO) at the first and second harmonics of the rotation period, P_ROT and P_ROT/2, and the realisation of the kernel G and its first derivative G′ as described in Section 4.1. We choose this implementation using P and P/2 as we see significant power in the GLS of the RVs at half of the rotation period, not only in ESPRESSO data but also in HARPS, HARPS-N and CARMENES data (see Fig. 2). We use wide priors for both the rotation period, P_ROT, and the timescale, T_ROT, with values (50, 300) and (3, 2) days, respectively. These two hyper-parameters are shared for the RV and FWHM simultaneous analysis. The posterior distributions have a narrow Gaussian shape centered on and . The timescale is quite short in comparison with the well defined rotation period. It is also shorter than those values found when analyzing individually each activity indicator separately where we found rotation periods in the range T_ROT = [138,275] d, and timescales in the range T_ROT = [101,249] d. This may indicate that the activity in RV and in the FWHM have different timescales but share the same rotation period. In the right panels of Fig. 3 one can see that the GLS periodograms of the FWHM and RV measurements in this GP-only model, after subtracting the offsets and adding the corresponding jitter terms to both datasets E18 and E19, concentrate the peaks in longer and shorter periods than 100 days for the FWHM and RV respectively. RV data exhibit half a rotation periodicity whereas FWHM data demonstrate the rotation, as shown before in the analysis of the FWHM only with a solution of P_ROT ~ 159 d and T_ROT ~ 101 d.

We also tested other GP implementations to model the activity using both the FWHM and RV measurements. We run the one-dimensional GP with a quasi-periodic (QP) kernel implemented within the George package (Foreman-Mackey et al. 2014). We also run the one-dimensional GP with a quasi-periodic and cosine (QPC) kernel, which integrates the period P in the quasi-period function and the period P/2 in the cosine function (Perger et al. 2021). However, for these two implementations, although the result of the median models was similar to the multi-dimensional GP with the SHO P and P/2 kernel, the posterior distributions were wider and required a normal Gaussian prior to provide a similar solution. We also tested other kernels within the multi-dimensional implementation, such as ESP and MEP kernels, and similarly to QP and QPC kernels, they all provide worse Bayesian evidence, with Δ ln < −4 in all cases.

The ESPRESSO RV measurements reveal a long-term variation clearly seen in the E19 dataset in Figs. 1 and 2. These variations, difficult to see by eye in other datasets such as HARPS, HARPS-N, and CARMENES due to instrument limitations, may indicate the presence of a long-term activity signal that we also see in the GLS periodograms in Figs. 1, 2 and 3. We use a longer dataset composed of ESPRESSO, HARPS, and HARPS-N data to verify this possibility. Thus, we run a simultaneous model of FWHM and RV measurements including the GP model previously described and adding a double sinusoidal model as used in the analysis of the ASAS-SN photometric data given in Equation (7). This global model is displayed in Fig. D.3, corresponding to model J1 in Table 3, including the double sinusoidal long-term cycle model, the GP model and a Keplerian model discussed in Section 6.3. The wide prior distributions, with values P_CYC (800, 5000), P_ROT (50, 300), and T_ROT ;(3, 2), and relatively narrow posterior distributions of the activity parameters are displayed in Fig. 7. This run gives the median values of the long-term cycle, , the rotation, , and the timescale, . These spectroscopically derived values perfectly match previous determination of the long-term cycle period and the stellar rotation period (Toledo-Padrón et al. 2019). We note a slightly longer timescale in this model although still at about one third of the rotation period. The global analysis including the CARMENES data without including any Keplerian model (model K1 in Table 3) provides almost the same result, with median values of the long-term cycle, , the rotation, , and the timescale, .

Fig. 7 Prior and posterior distributions of the long-term cycle, the GP rotation period, and the timescale from the global analysis of FWHM and RV measurements of ESPRESSO, HARPS, and HARPS-N (model J1 in Table 3).

Fig. 8 GLS periodograms of ESPRESSO CCF FWHM measurements (top), TM RV measurements (middle), and the window function (bottom) of GJ 699 after subtracting the median of each dataset before (E18) and after (E19) the intervention in June 2019.

6.2 Candidate planetary signals

In the RV residuals of the simplest activity-only model shown in Fig. 3 (model A in Table 3), described in Section 6.1, we find several signals in the GLS periodogram with a false-alarm probability (FAP) of less than 1%, which we cannot attribute to any activity process, and we tentatively associated them with candidate planetary signals. We identify three of them as main signals at 3.15 d (with 1.46 d and 0.76 d as 1 d aliases, and 0.59 d as 1 d alias of the 1.46 d signal), 4.12 d (with 1.32 d as 1 d alias), and 2.34 d (with 1.74 d as 1 d alias). Figure 8 shows the GLS periodogram of the ESPRESSO FWHM and RV measurements, and the window function computed as the periodogram of all values equal to 1 at the BJD values of the ESPRESSO FWHM and RV measurements. We note strong peaks at 1 d and 1 yr in the window function and we do not see any significant peak at the position of the main candidate planetary signals. The GLS of the FWHM and the RV measurements are dominated by the activity signals at about 247 d and 67 d, respectively, but there are no significant peaks at the position of the candidate planetary signals. In Fig. 3, we display the GLS after adding the jitter of the corresponding activity-only model which explain the slight differences with those in Fig. 8. The GLS after subtracting the GP model in the FWHM does not show either any significant peak in the whole period range from 1 to 10,000 d (see Fig. 3).

6.3 Evaluating the 3.15 d signal

We evaluate the strongest signal at 3.15 d in the RV residuals after subtracting the GP model (see Fig. 3). In Fig. 9 we show the evolution of both the 3.15 d signal and the 1 d alias at 1.46 d as a function of number of RV points. To do that, we run the activity-only model for both the FWHM and RV measurements from 50 RV points and adding ten points in each run until we end up with all 149 points. We clearly see a steady increase in confidence of both signals. In particular, the 3.15 d signal, reaches in the last run a FAP of about 10⁻⁶, being already very significant (FAP < 0.1%) after 110 ESPRESSO observations. In all these runs, we adopt wide priors on P_ROT (50, 300) d and on T_ROT (3, 2). We see how the estimated rotation period is converging towards the final value of of this model A in Table 3. The 3.15 d signal is certainly slightly diminished due to the GP modeling since we are not including a Keplerian model targeting the 3.15 d signal in this run.

We first include a circular planet model (model C1c in Table 3) in a new run to search for that particular signal but again with a wide prior (0.5, 50) d on orbital period, P_orb. We use the time of conjunction given by the phase, ϕ, with a prior (−0.5, 0.5), centered around the maximum time, t_max, of the observation baseline as in Eq. (8): (8)

and the semi-amplitude velocity, k_p, with a uniform prior (0,5) m s⁻¹. We detect unequivocally the 3.15 d signal as a planetary signal with a Bayesian evidence ln = −443.7 (model C1c in Table 3), and Δ ln = +5.6, with respect to the activity-only model A, which corresponds to a 1/e^+5.6 = 0.37% false alarm probability for the activity+planet model. We run this model five times to check any possible ln variance and found ln in the range [−442.7,−445.3], with a mean ln of −444.2 (σ = 1.0). This run converges to P_orb = 3.1533 ± 0.0005 d and . We note that in Table 3 all the differences Δ ln are always referred for simplicity to the activity-only model D, although in some cases, as in this particular case, it is another model (model A) the reference activity-only model.

We also run a Keplerian planet model (model C1e in Table 3) given in Eq. (9): (9)

where the true anomaly η, which is the angle between periastron and the planet, measured from the barycenter of the system (e.g. Eastman et al. 2013), is related to the solution of the Kepler’s equation, and depends on the orbital period of the planet P_orb, the orbital phase ϕ, and the eccentricity of the orbit e. We use the RadVel toolkit to implement the Keplerian model (Fulton et al. 2018). We associate the orbital phase ϕ with the time at inferior conjunction, T_c, which depends on the maximum time, t_max, of the observation baseline as in Eq. (8). Then, we convert this T_c time into time of periastron, T_p, which depends on the argument of periastron ω and the eccentricity e. Following Eastman et al. (2013), we parameterise the eccentricity as in Eq. (10): (10)

Initially we tried an eccentricity with uniform prior (0, 1), but the solution converged to a very low value consistent with zero. Thus we decided to sample and with normal priors (0, 0.3), still allowing the possibility of high eccentricity values but favoring low eccentricity values, expected for the short periods of the candidate planets. Eccentricity is typically overestimated in noisy data and datasets with unmodelled sources of variability (Hara et al. 2019). We measure a ln of −445.5 (model C1e in Table 3), which corresponds to a Δ ln of +3.8 (equivalent to 2.2% FAP), with respect to the reference activity-only model (model A in Table 3).

Fig. 9 Evolution of the false-alarm probability of the 3.15 d and 1.46 d (1 day alias) signals (top) and the GP rotation period (bottom) with the number of observations. The 0.1, 1, and 10% FAP lines are computed using the full 149 ESPRESSO dataset.

Fig. 10 GLS power spectral density of the observed signals 3.15d and 1.46d after removing the GP compared to a simulated 54 cm s−1 injected sinusoidal signal at different orbital periods.

6.4 Statistical tests to validate the 3.15 d signal

We ran several tests to validate the 3.15 d signal, using the model C1e in Table 3 as a reference model. Firstly, following Rajpaul et al. (2016), we tested the injection and recovery of different sinusoidal signals with k_p = 54 cm s⁻¹, with random phases (−0.5, 0.5), and periods uniformly distributed within 3.15 ± 2.0 d, excluding 0.1 d around 3.15 d and its 1 d alias 1.46 d. We performed 10 000 simulations and the result is displayed in Fig. 10. We randomly selected 10000 different samples among the resulting 124 524 samples of the run of model C1e. From each sample, we built the GP model using the parameters of that sample, subtracted it from the observed RVs and injected one simulated sinusoidal signal with period P_x, random phase, and a semi-amplitude of k_p = 54 cm s⁻¹. We then computed the GLS periodogram and measured the power spectral density (PSD) of the signals at the 3.15 d, 1.46 d, and P_x. The distribution of PSD of the simulated planet x appears to perfectly match that of signal 3.15 d, and both show larger values of PSD distributions than that of the 1.46 d signal. This indicates that indeed the 1.46 d signal is just the 1 d alias of the main signal at 3.15 d. We note that the PSD distribution of the 3.15 d signal has a long tail below PSD = 0.2 due to the fact that about 10% of the samples converge to a detection of 4.12 d signal or the 1.46 d signal as the main signal in the GLS periodogram, thus decreasing the PSD of the 3.15 d signal. This, and the fact that we adopted a Keplerian model in this run, may explain why the PSD of the injected planet x is larger than that of the 3.15 d signal in about 60% of the simulations.

Secondly, we ran the search of a Keplerian solution with the same priors as in model C1e in Table 3, but we randomly removed 15 points (10% of data points) from the datasets of FWHM and RV measurements. We repeated this 25 times and recovered the 3.15 d signal in 18 cases (72% with the 1 d alias 1.46 d in 3 cases) and the 4.12 d signal in 7 cases (28% with the 1 d alias 1.32 d in 1 case).

Thirdly, we ran another simulation using again the 10000 samples of the run model C1e in Table 3, where we built the GP using each sample parameter. For each iteration of these 10000 samples, we computed the GLS of the RVs, corrected for offsets, and, with the jitter added to the RV uncertainties and after subtracting the GP, measured the PSD at 3.15 and 1.46 d. The distributions of these PSD values are labelled as ‘observed’ in Fig. B.1. We see that the observed PSD distributions of the RVs measured at 3.15 d and 1.46 d have some structure at lower PSD values, which is related to a significant fraction of the samples where the solution peaks at 4.12 d (10.85%) and at 2.34 d (0.08%) instead of 3.15 d (89%). This is consistent with the aforementioned test. As expected, a blind search for the planetary signal recovers in most cases the 3.15 d signal and thus the median of the PSD distribution of the observed 3.15 d signal is larger than that of the PSD distribution of the observed 1.46 d signal. We refer to Appendix B for further discussion of Fig. B.1.

Fourthly, we run another test where we model the activity of the star using only the RVs with a multi-sinusoidal function composed of six sinusoidals with periods in the range [60,255] d (see Fig. B.2). We found that six sinusoidal functions were sufficient to be able to reproduce all possible RV variations due to stellar activity. The fitting procedure was done following a prewhitening methodology, which is based on a single sinusoidal fit every time at the period of highest peak (lowest FAP) in the GLS periodogram followed by subtraction of this fit and repeating this sequence. The fitted periods were 79.1, 254.9, 68.2, 89.8, 59.5, and 244.9 d. This activity-only model resembles the GP-only model displayed in Fig. 3 (model A in Table 3). The GLS periodogram of the RV residuals after subtracting the multi-sinusoidal fit also exhibits the detection of the 3.15 d signal above the 0.1% FAP line.

We subsequently performed an additional test, where we modelled the RVs and all activity indicators using a moving average model with an exponential decay using a 5.25 d window. The aim of this test was to run a simplified model with less freedom than GP models. To choose this 5.25 d window of the moving average, we used the FWHM, trying to avoid over-fitting. This 5.25 d exponential moving average model was then applied to the RVs and all activity indicators. The result of this test is depicted in Fig. B.3, where we are able to recover again the 3.15 d signal in the GLS periodogram of the RV residuals after subtracting this activity-only model using the exponential moving average. This 3.15 d signal does not appear in the residuals of any of the activity indicators after fitting the same activity-only model using the 5.25 d exponential moving average. All GLS periodograms of all activity indicators have all GLS peaks below 10% FAP line.

Finally, we performed a last test where we build several simulated RV time series: (i) a first flat model to evaluate the sampling, with flat result in the GLS periodogram; (ii) a second one using the activity-only model with injected white noise at 30 cm s⁻¹ level, resembling the moving average model, again resulting in a flat GLS periodogram with all peaks below the 10% FAP line; and (iii) a third one by injecting a planetary signal at 3.15 d (see Fig. B.4), thus in this case recovering the 3.15 d signal. We note the difference between the GLS periodograms of the residuals of real and simulated RVs, pointing that the real data shows additional peaks such as the 4.12 d signal although at the 10% FAP level.

All these tests allow us to confidently confirm the 3.15 d signal as a planet signal.

Fig. 11 Prior and posterior distributions of the period and semi-amplitude from the global analysis of FWHM and RV measurements of the blind search of the 233d candidate planet using ESPRESSO, HARPS, HARPS-N and CARMENES data (left, model O1 in Table 4), a blind search (centre-left, model N2 in Table 4) and a guided search (centre-right, model N1 in Table 4) of the 233d candidate planet using ESPRESSO only, and a blind search of the 233d simulated signal (right, model S1 in Table 4).

6.5 Evaluating the 233 d candidate super-Earth

Ribas et al. (2018) reported the detection of a candidate super-Earth-like planet orbiting in an slightly eccentric orbit with a period of 232.8 ± 0.4 d, with a semi-amplitude velocity of k_p = 1.2 ± 0.1 m s⁻¹. This result has been challenged by Lubin et al. (2021); Artigau et al. (2022) which argue that the signal is connected to stellar activity as the planet candidate period is close to 1 yr alias of the rotation signal (Toledo-Padrón et al. 2019). We run several models with the RV and FWMH data to search for the 233 d signal. Firstly, we use the whole dataset of ESPRESSO, HARPS, HARPS-N, and CARMENES to run a blind search for the 233 d signal (model O1 in Table 4), with a uniform prior /(200,250) on the orbital period. The posterior distribution of the orbital period and semi-amplitude velocity are depicted in Fig. 11. This model is not satisfactory since the semi-amplitude, , is consistent with zero at 1.8σ, and significantly lower than that of the reported candidate planet. The model also exhibit a Bayesian evidence Δ ln = −0.5 with respect to activity-only model, but the orbital period, , is shorter than that of the candidate planet although still consistent within the error bars.

Secondly, we run a blind search only with ESPRESSO data, with the period free (model N2 in Table 4) in the same range /(200, 250), and in this case the posterior of the orbital period is flat as the prior, and the semi-amplitude is even lower and consistent with zero (see central left panels in Fig. 11). Thirdly, we run a guided search (model N1 in Table 4) with ESPRESSO only with the period prior at (233,0.5) and the result is similar to the blind search, the planet is not detected with a worse Bayesian evidence in both cases with Δ ln ~ −2.5. One could argue that the number of data points and baseline is not enough to detect such a long orbital period, but the semi-amplitude is really significant for ESPRESSO to have missed the candidate planet.

Finally, we run a simulation by adding a planetary signal in the activity model with injected white noise using the uncertainties including the jitter terms (with a mean value of 0.7 m s⁻¹) of the previous run (model N1 in Table 4). We injected a Keplerian with the same eccentricity (very close to zero) of the Keplerian solution of model N1, a period of P_orb = 233 d, with a semi-amplitude velocity of k_p = 1.2 m s⁻¹. We run a blind search with priors of P_orb (200, 300) d and k_p (0, 5) m s⁻¹. The posteriors are P_orb = 235 ± 6 d and k_p = 1.4 ± 0.4 m s⁻¹ (with the k_p only detected at 3.5σ) which are consistent with the injected Keplerian signal but not exactly peaking at right values (see right panels in Fig. 11). This result (model S1 in Table 4) is probably related to the baseline and number of RV epochs, as well as to the search for a signal at longer period than the activity-induced signal by stellar rotation. We explore if using a conservative white noise of the input model had a significant impact, so we run one last simulation (model S2 in Table 4) similar to the previous one but with the injected white noise at 0.3 m s⁻¹, and we get a slightly improved result whose posteriors are very similar with median values P_orb = 235 ± 5 d and k_p = 1.3 ± 0.3 m s⁻¹ (at 4.3σ). Therefore, the quality of the ESPRESSO RVs would have allowed us to clearly detect such a relatively strong signal as that of the candidate planet reported in Ribas et al. (2018). Thus, we conclude that ESPRESSO data do not support the existence of the candidate planetary signal.

Fig. 12 ESPRESSO FWHM measurements (top), RV measurements (middle), RV residuals (next bottom) from long-term cycle (red thick dashed line) and SHO (PROT and PROT/2) GP model (grey solid line), and RV residuals (bottom) from the Keplerian model (grey shaded area), and GLS periodograms (left) of GJ 699. The uncertainties include the jitter term coming from the global model E1e in Table 3, with prior and posterior parameters given in Table 6.

6.6 The 3.15 d sub-Earth mass planet

ESPRESSO RVs show a long-term variation, in particular, in the E19 dataset, which may be interpreted as a long-term activity cycle. We then run a model for the ESPRESSO dataset (model D in Table 3) with the GP and the cycle as done in Section 6.1 but with priors on P_CYC (3250, 300) (normal prior centered on P_CYC = 3250 d with a σ = 300 d), on P_ROT (50, 300) (uniform prior on P_ROT in the range [50,300] d), and on T_ROT (3, 2) (log-normal prior centered on log T_ROT[d] = 3 with a σ(log T_ROT[d]) = 2). This is supported by models I1 and J1 in Table 3, which use a longer baseline of ESPRESSO, HARPS and HARPS-N (see Fig. 7). Model D results in a Bayesian evidence of ln = −447.2 which we adopt as our reference model for the ESPRESSO dataset. The posteriors of this model give , the rotation, , and the timescale, . We run three times this activity-only model to check any possible variance on ln values and we found for these three models equivalent ln values at −447.1 ± 0.1. We also explore other possibilities as testing the GP plus first, second and third order polynomials (models B1, B2 and B3 in Table B.1), but they provide worse Bayesian evidence with Δ ln = −6.2, −6.7 and −15.5, respectively. The same behaviour is seen with ESPRESSO, HARPS and HARPS-N dataset with the model G3 compared to I1 and I2 in Table B.1.

We run a blind search of a planet using either a circular or a Keplerian model (E1c and E1e in Table 3) with a prior P_orb (0.5, 50) d. The cycle+GP+Keplerian model, E1e, displayed in Fig. 12, uses the and with normal priors (0, 0.3). We note that even if this eccentricity prior favors low eccentricity values, it still leaves the possibility to choose high eccentricity values if needed. In Table 5 we provide RV and FWHM measurements of ESPRESSO datasets used in Fig. 12. We detect unequivocally the 3.15 d signal as a planetary signal with a Bayesian evidence ln = −438.6, and Δ ln = +8.6 with respect to the activity-only model (model D in Table 3), which corresponds to a 0.018% FAP for the activity+planet model. We repeated 10 times the run for model E1e and found some variance in the ln values in the range [−438.2,−444.7], with a mean of −442.2 (σ = 2.1). In all these runs, we clearly detect the 3.15d planet but the planet posteriors are slightly different. The runs with lower ln values produce slightly asymmetric posterior distribution for the planet semiamplitude k_p and also few per cent of the planet period samples (typically less than 5–10%) corresponds to either aliases of 3.15d or other signals (typically 4.12d and its alias). We believe this is the main reason why we find some variance in the ln values. We adopted as our final solution for model E1e one of the runs with higher ln value at −438.6 which provides a fairly symmetric posterior for the planet semi-amplitude k_p and all the posterior period samples at 3.15d (see Fig. 13).

Figure 12 depicts the FWHM and RV ESPRESSO measurements together with the activity model and the GLS periodograms before and after subtracting the long-term cycle. We note that the SHO (P and P/2) GP (model E1e, with priors and posteriors parameters given in Table 6) is efficiently absorbing all GLS peaks at periods longer than 50 d, even with the P_ROT = 136 ± 10 d, for both the FWHM and RV ESPRESSO datasets. In the lower left panels of Fig. 12 we show the RVs after subtracting the activity model and below that panel the RV residuals after subtracting the Keplerian model. In the lower right panels, the GLS periodograms show clearly the 3.15 d signal with PSD greater than 0.4 well above the 0.1% FAP line. In the GLS periodogram of the RV residuals after subtracting the Keplerian model (lower right panel of Fig. 12), the 3.15 d signal together with the 1 d alias at 1.46 d and 0.76 d and 0.59 d (1 d alias of 1.46 d) have disappeared. Other signals, already mentioned before, appear at 4.12 d (with 1 d alias at 1.32 d), 2.34 d (with 1 d alias at 1.74 d), still above the 0.1% FAP line. At periods below 1 d, a few new peaks that have grown are all 1 d alias of the signals at 4.12 d and 2.34 d and their 1 d alias, but still below the 10% FAP line. The RMS of the RVs goes from the original value at 1.86–0.60 ms⁻¹ after subtracting the activity model, ending at 0.45 ms⁻¹ in the RV residuals after subtracting the Keplerian model. The median jitter of the E18 and E19 datasets in model E1e are 0.56 m s⁻¹, compared with the median S-BART RV uncertainties of 0.10 ms⁻¹. It is remarkable that the RMS of the RV residuals is 0.45 m s⁻¹, greater than the RMS of the FWHM residuals of 0.35 m s⁻¹, thus requiring an additional 0.56 m s⁻¹ RV jitter while a 0.30 m s⁻¹ FWHM jitter in the E19 dataset (see Table 6). This suggest that this remaining RMS in RV may be possibly related to additional candidate planetary signals present in the data and not included in this model E1e.

Figure 14 displays the RV measurements versus orbital phase. The planetary signal at 3.15 d is not eccentric. The E1 e model converges to P_orb = 3.1533 ± 0.0006 d and k_p = 55 ± 7 cm s⁻¹ (~7.9σ detection) and an eccentricity less than 0.16, consistent with zero. The planet parameters, hereinafter referred as Barnard b, are provided in Table 7. Fig. 13 shows the corner plot with the minimum mass of the planet equal to m_p sin i = 0.37 ± 0.05 M_⊕, i.e. a sub-Earth mass planet with about half of the mass of Venus and about three times the mass of Mars.

We verified that we were able to recover the planet using the DRS RVs and the LBL RVs, with semi-amplitude velocity and for DRS and LBL runs (models E1eD and E1eL in Table 3), respectively, thus all compatible with TM run within the error bars (see also Appendix B and Table B.1). We run several additional models with the TMtc RVs derived using the telluric corrected ESPRESSO spectra (see Section 2.1). We first run an activity-only model with only a GP (model AT in Table 3). We note that the Bayesian evidence cannot be directly compared between models A, AD, AL and AT, but we provide Δ ln always with respect to model D in Table 3 for simplicity. Activity-only model DT shows better Bayesian evidence than model AT. We also recover the 3.15 d planet using the TMtc RVs in the runs E1cT and E1eT with Δ ln = 6.0 and 7.5 with respect to the activity-only model DT, corresponding to 0.25 and 0.05% FAP. In these cases, the semi-amplitude velocity is and for the circular and Keplerian orbits respectively. In the following we keep using the TM RVs as reference ESPRESSO RVs, since we do not see significant improvement, and we do consistently the RV computation with HARPS and HARPS-N, without using telluric corrected spectra, but masking out the regions possibly affected by telluric features taking into account BERV changes.

We also perform several runs with the S-BART code (Silva et al. 2022), to explore any possible chromatic effects on the 3.15 d signal. We compute RVs using either only the blue-detector or the red-detector spectra, or by dividing the whole ESPRESSO spectral coverage into three parts with approximately identical RV content, being the blue, green and red RVs. We run models using the red-detector RVs, with priors on orbital period P_orb (0.5, 50) d and (0.5, 20) d, providing clear detections of the 3.15 d signal, with posteriors of the planet and activity parameters with Gaussian shapes peaking at the values very similar to reference model E1e in Table 3. The blue-detector RVs, and the blue, green, and red RVs, with the prior P_orb (0.5, 20) d, were also providing a detection of the 3.15 d signal. In all these models the activity parameters give median values in the ranges for the cycle P_CYC = [3254, 3498] d, and the rotation P_ROT = [136, 144] d, whereas for the semi-amplitude velocity these models give k_p = [0.40, 0.54] m s⁻¹. We consider these results reasonable given the low semi-amplitude of the signal, which makes it difficult to be detected. The red detector seems to contain most of the RV content, according to the shape of the posterior distributions on semi-amplitude. These two runs provide a result very similar to the reference model, which seems reasonable given the late spectral type of the star that shifts the stellar flux towards the red part of the spectrum.

The analysis of the ESPRESSO, HARPS and HARPS-N data together reinforces the confirmation of the 3.15 d signal as a sub-Earth mass planet. Figure D.3 shows the result of a model consisting of long-term cycle, a GP and a Keplerian model. This run (model J1 in Table 3) uses wide priors on P_CYC (800, 5000) d, P_ROT (50, 300) d, P_orb (0.5, 50) d. This run J1 gives a Δ ln = +12.4 with respect to the activity-only model I1, equivalent to 0.0004% FAP. This model leaves in the GLS RV residuals (bottom right panel in Fig. D.3) the tentative signal at 4.12 d above 1% FAP line. We believe that this is the result of the balance between including more RV data with a longer baseline but, on the other hand, poorer quality data according to the lower RV precision of the HARPS and HARPS-N spectrographs. Figure D.1 displays the RV data versus orbital phase, where the RMS of the original ESPRESSO, HARPS and HARPS-N RVs of 1.93 m s⁻¹ goes down to 0.97 m s⁻¹, after subtracting the activity-only model, and down to 0.91 ms⁻¹ after subtracting the 3.15 d Keplerian RV curve. Figure D.2 shows the corner plot of the posteriors of the planet parameters of this run, very similar to the ESPRESSO-only run E1e. The posteriors give again a practically null eccentricity, and the values P_orb = 3.1537 ± 0.0004 d and k_p = 48 ± 8 cm s⁻¹ (6σ detection), with m_p sin i = 0.33 ± 0.06 M_⊕, consistent with ESPRESSO-only model E1e within the uncertainties. Figure D.4 depicts the prior and posteriors distributions of the 40 parameters of model J1.

Finally, we run the model of FWHM and RV measurements including the CARMENES data (see Fig. D.5) that also confirms the detection of a planetary signal at 3.15 d. This model L1 in Table 3 shows again a difference on the Bayesian evidence of Δ ln = +12.3 with respect to the activity-only model K1. We also verify that using uniform or normal priors for the cycle provides very similar results. Now, again, the 4.12 d and 2.34 d signals appear in the GLS of the RV residuals above the 0.1% FAP line. The RMS of the RV measurements in this run K1 involving ESPRESSO, HARPS, HARPS-N and CARMENES data goes from the original values at 2.08 m s⁻¹ down to 1.49 m s⁻¹, after subtracting the activity-only model, and down to 1.47 m s⁻¹, after subtracting the Keplerian 3.15 d signal. The 3.15 d signal is too weak to be independently confirmed using CARMENES data only, using H15+HAN data only, or using both datasets together. A search with a very narrow prior at 3.15 d ±0.1 d provides a tentative detection but weaker than 10% FAP. Other candidate signals at 4.12 and 2.34 d are also present, but still less significant than 10% FAP.

Fig. 13 Corner plot with the posterior distributions of the orbital parameters of the sub-Earth-mass planet of GJ 699 with a 3.15 d orbital period, with prior and posterior parameters given in Table 6.

Fig. 14 RV curve of the sub-Earth-mass planet of GJ 699 with a 3.15 d orbital period together with ESPRESSO RVs, with uncertainties with (light colour) and without the jitter term (dark colour) coming from the global model E1e, with prior and posterior parameters given in Table 6.

6.7 Candidate multi-planetary system

We investigate the possibility to detect and confirm in a blind search these additional signals seen in the GLS periodogram of the ESPRESSO RV residuals. We note that the high required jitter of 56 cm s⁻¹ invites one to believe that this is the consequence of including only one Keplerian signal in the model. Thus we first run several models including two Keplerian signals with priors on orbital period P_orb (0.5, 20) d, where we only put the restriction that planet b should have shorter period, and therefore shorter orbital distance, than planet c. In addition, we only permit those orbits where the eccentricity of planet b is such that excludes the orbital collision (e.g. Laskar & Petit 2017), given by the Eq. (11): (11)

where P_i, a_i and e_i are the orbital period, distance and eccentricity of planet b and c. We note that the prior on eccentricity e with and (0,0.3), though it favors low eccentricity orbits, it leaves some freedom to reach high eccentricity values. Models E22c and E22e in Table 3 shows a satisfactory result in terms of the shape of the posteriors of the activity and planetary parameters of planet b, but not for planet c, where the semi-amplitude gets lower median value than expected and consistent with zero in model E22e. The result is better for the circular model E22c with the posteriors with better shape and median values of P_b = 3.1533 ± 0.0006 d and P_c = 4.1249 ± 0.0014 d, and and . However, the Bayesian evidence of these models is worse than that of the models E12c and E12e. We note here that we computed, as we did for model E1e, five times the model E12c and E12e and we see there is some (slightly lower) variance in the ln in the range [−438.8,−442.5] for E12c, with a mean of −441.1 (σ = 1.5), and in the range [−439.8,−442.3] for E12e, with a mean of −440.9 (σ = 1.1). Therefore, from a blind search of two planetary signal using ESPRESSO only FWHM and RV data we cannot confirm the presence of the planet c signal at 4.12 d. We also tried additional two-planet model runs by adding the datasets of HARPS and HARPS-N (model J2) and adding the CARMENES dataset (model M2). The Bayesian evidence of these models were a bit worse (model J2) or a bit better (model M2) than the single-planet model in each case, but not significant. It is remarkable that for the M2 model the periods and the semi-amplitudes, with median values of P_b = 3.1543 ± 0.0003 d and P_c = 4.1245 ± 0.0005 d, and and , are marginally compatible with the ESPRESSO-only run E22c, within the uncertainties.

In Appendix E we describe the exercise to run a guided search of the four candidate signals at 3.15, 4.12, 2.34 and 6.74 d, corresponding to models F1, F2, F3, and F4 in Table E.1.

7 Discussion

We take advantage of the result of run F4 in Table E.1 (described in Appendix E) to test the presence of planets, and assess their significance with the False Inclusion Probability (FIP; Hara et al. 2022). This framework uses the posterior distribution of the nested sampling run, and computes the probability of having a planet within a certain orbital frequency interval based on the fraction of samples within that frequency interval. It can test the presence of an arbitrary number of planetary signals. To produce the posterior distribution, we run a ESPRESSO-only model composed of the same activity model as described before, and four circular planetary signals. For all these signals we use the same priors, with wide log-uniform priors for periods P_orb (2, 40) d days, and uniform prior on semi-amplitudes k_p (0, 5) m s⁻¹. Fig. 15 shows the results of the − log₁₀ FIP as a function of period. The results support the detection of a planet at a period of a 3.15 d, and a solid hint of the presence of a second planet at 4.12 d. This result is expected from the previous discussion in Sections 6.4, 6.6 and E.

We also estimate our detection limits for additional companions in the ESPRESSO data by performing a simple injection-recovery test. We construct activity-only RVs by generating random samples from the GP model, and adding white noise with the same standard deviation of the residuals of the 4-planet model (see Appendix E). We then inject 100000 random sinusoidal signals with periods between 1 and 1100 d, semi-amplitudes between 1 cm s⁻¹ and 11 m s⁻¹, and random phases. We reject those combinations that would create an RMS of the data significantly higher than the real RMS of our data. Then we subtract the stellar activity using the same GP hyper-parameters as measured for our four-planet model, but recomputing the activity model for each specific dataset. Then we generate the periodogram of the residuals and check the false alarm probability of the highest peak at periods within 10% of the injected period.

In Fig. 16 we display the result of this exercise. With this dataset and activity model, we are sensitive to signals of ~0.3 m s⁻¹ for periods from 1.2 d up to ~20 d. Our sensitivity drops to ~1 m s⁻¹ at 35 d and then quickly plummets to 2.5 m s⁻¹ at periods longer than 40 d. These values correspond to minimum masses of ~0.5 M_⊕ at 15 d, ~1.5 M_⊕ at 35 d and > 4 M_⊕ at periods longer than 40 d. This is a typical effect of GP-only models, which try to account for all variations at low frequencies up to the characteristic value specified by the kernel. To extract possible low-amplitude signals at longer periods we would need either a very significant amount of new data or a different strategy to mitigate stellar activity. This also shows that finding the signal of Barnard b (3.15 d) is within expectations. Finding the signals at 2.34 and 4.12 d is also within expectations, with the 4.12 d signal being easier to find. Finding the signal at 6.74 d, however, is not within expectations. Exploring the region of periods longer than 30–40 d likely requires simultaneous modeling of the GP and activity signals using a Nested Sampling, or MCMC, algorithm. We discuss in Appendix E the candidate multi-planetary system.

The sub-Earth mass planet orbiting at period 3.15 d our second closest neighbour Barnard’s star has a minimum mass of 0.37 ± 0.05 M_⊕. Fig. 17 compares the confirmed planet Barnard b with those other planets in the exoplanet database. We highlight the ESPRESSO planet discoveries and confirmations with particular emphasis on our closest neighbour star, hosting the two planets Proxima d and b (Suárez Mascareño et al. 2020; Faria et al. 2022). ESPRESSO has gone beyond the Earth mass frontier, thus starting to discover sub-Earth mass planets orbiting nearby stars to our Sun such as Proxima Cen and Barnard’s star, although still at short orbital periods. These discoveries have required significant investment with more than 100 measurements in an 8-m class telescope. The endeavor of filling the region in Fig. 17 below the 1 m s⁻¹ precision limit at periods larger than 10 d towards finding long-period Earth-like planets would require many years of continuous, well planned, sufficiently sampled observations and a deep understanding of magnetic activity of every nearby star. This would open the opportunity to sample the habitable zone of nearby stars with the aim of further investigating the atmosphere of Earth-like planets in the future with ground-based facilities such as ANDES at ELT (Marconi et al. 2022) and with space missions such as LIFE (Quanz et al. 2022).

Stringent upper limits on the presence of planetary mass companions at larger separations can be placed using the catalogues of Hipparcos-Gaia astrometric accelerations produced by Brandt (2021) and Kervella et al. (2022). For example, using the analytical formulation of Kervella et al. (2019), Fig. 18 shows the sensitivity curve for Barnard’s star given the measured PMA in the Hipparcos-Gaia DR3 catalogue of Kervella et al. (2022). The peaks, corresponding to regimes of orbital separation with reduced sensitivity, are a direct consequence of the observing window smearing effect of the orbital motion of any companion. This effect is intrinsic to the proper-motion difference technique, where the catalogue proper-motion values for Hipparcos and Gaia are averaged over the time-span of the observations, thus with smearing orbital effects for shorter periods. The last peak in sensitivity loss appears, in fact, at the exact time-span of the observations of the astrometric catalogue we consider. For Gaia DR3, this means 2.83 years, i.e. a separation of 1.08 AU for a star with the mass of Barnard’s star. Therefore, there are also changes to the precise location of the regions of sensitivity loss due to the fact that the time-span from DR2 to DR3 over which the Gaia proper-motion values are computed varies. In the region of uniform sensitivity, the gain in mass limits is a factor of several, almost one order of magnitude in comparison with Figure A.3 in Kervella et al. (2019). Thus Fig. 18 shows that in the regime of approximately uniform sensitivity 2–10 AU, the presence of super-Earths with true masses around 10 M_⊕ can be ruled out by the lack of statistically significant PMA. The combination of ESPRESSO RVs and Hipparcos-Gaia absolute astrometry thus allows to infer that super-Earths or larger planets are not present in the GJ 699 system out of ~10 AU.

Fig. 15 FIP periodogram of the 4-planet model of the ESPRESSO data of GJ 699. The periods of the detected signals are indicated in red solid circles. The − log10 FIP as a function of the center of the period is represented as a blue line. The green dashed and dotted orange lines show the 1 and 10% FIP thresholds, respectively.

Fig. 16 Detection limits. Top panel: semi-amplitude kp versus period Porb for our detection limits based on injection and recovery tests. The white region shows the signals that we would detect in a periodogram with 1 % FAP. The blue region shows signals we would detect with a FAP lower than 1 % FAP. With the current dataset and activity model, we would miss the red region signals. The top grey dashed line shows the maximum semi-amplitudes that would still be compatible with the data. The bottom grey dashed line shows the minimum amplitude that we expect to detect with ESPRESSO data and the current exposure times. The green solid circle shows the confirmed planet and the orange circles show the position of the candidate planets. Bottom panel: same as top panel but in mass versus period.

Fig. 17 Minimum mass vs. orbital period diagram for known planets from the NASA Exoplanet archive as of December 2023 orbiting solar-type stars, together with those discovered and confirmed using ESPRESSO (green circles). Confirmed planet Barnard b (red), and planet candidates (purple) from the 4-planet candidate system orbiting Barnard’s star, together with the two planets orbiting Proxima Cen, Proxima b (yellow) and d (blue) are highlighted. Planets of the solar system (grey circles) are also labelled. Inclined solid and dashed lines show the RV semi- amplitude of planets orbiting a late M dwarf star with 0.25 M⊙ and a G dwarf star with 1 M⊙ star assuming a RV semi-amplitude of 1 cm s−1 (green line) and 10 cm s−1 (blue line), respectively, and null eccentricity.

Fig. 18 Hipparcos-Gaia PMA sensitivity to companions of a given mass (in M⊕) as a function of the orbital semi-major axis (in AU) orbiting Barnard’s star. The solid black line identifies the combinations of mass and separation explaining the observed proper-motion anomaly (PMA) at the mean epoch of Gaia DR3 (Eq. (15) in Kervella et al. 2019). The shaded light-blue region corresponds to the 1σ uncertainty region.

8 Conclusions

We analysed 156 ESPRESSO spectra of Barnard’s star, the second closest stellar system to the Sun after the α Centauri stellar system, at a distance of about 1.8 parsecs. These spectra have been taken with the ESPRESSO spectrograph at the VLT as part of the GTOs over the last five years. The ESPRESSO RV data show a 1.8 ms⁻¹ RMS with about 8 m s⁻¹ peak-to-peak RV measurements, which mainly reflect the magnetic activity of this old middle-M-type dwarf star. We modelled these RV measurements together with CCF FWHM measurements, with a 15 m s⁻¹ peak to peak and 3.3 m s⁻¹ RMS, using a long-term cycle described with a double sinusoidal model, and a rotation-induced activity model, which is described using a Gaussian process approach. The activity cycle is better constrained using additional HARPS, HARPS-N, and CARMENES data, which help to extend the baseline of the observations up to eight years. We obtained a well-described activity model with a cycle period of and a rotation period of , which are consistent with previous results in the literature (e.g. Toledo-Padrón et al. 2019).

The high quality of the ESPRESSO data means that we can evaluate the presence of the candidate super-Earth-like planet at 233 d reported in Ribas et al. (2018). We are unable to recover the 233 d signal in either a blind search or a guided search with a narrow prior on the planet period (233, 0.5) d. We do not recover a clean signal above 0.50 m s⁻¹, with all the models giving solutions with semi-amplitudes consistent with zero. We also tested simulated data at the ESPRESSO precision and conclude that ESPRESSO data do not support the existence of the 233 d candidate planet.

The RV residuals of the activity model reveal several signals at periods of shorter than 10 d; in particular four signals at periods of 3.15, 4.12, 2.34, and 6.74 d, sorted here according to signal strength. We assessed the nature of the main signal at 3.15 d, and confirm it as a planetary signal. We modelled the 3.15 d signal with a Keplerian, resulting in an almost circular orbit with a semi-amplitude of 55 ± 7 cm s⁻¹, thus uncovering a sub-Earth-mass planet of 0.37 ± 0.05 M_⊕, which is about half of the mass of Venus or three times that of Mars. This sub-Earth-mass planet is located inwards from the HZ of the star, with an equilibrium temperature of 400 K – assuming an albedo of 0.3.

We are unable to confirm the other signals in a blind search. However, we ran a guided search of these signals modelled as Keplerian with narrow priors on each of the four signals at P_orb = 3.15, 4.12, 2.34, and 6.74 d, assuming very low eccentricities, recovering a candidate four-planet system with semi-amplitudes of k_p = 47, 41, 35, and 20 cm s⁻¹, which would correspond to a system of four sub-Earth-mass planets with m_p sin i = 0.32, 0.31, 0.22, and 0.17 M_⊕. All candidate planetary orbits would be located inwards from the HZ of the star, with orbital semi-major axes of between 0.019 AU and 0.038 AU. Thus, all the candidate planets would be irradiated to a greater extent than the Earth, with incident fluxes of between 2.4 S_⊕ and 10.1 S_⊕, and their equilibrium temperatures, assuming an albedo of 0.3, would be in between 440 K for the inner planet and 310 K for the outer planet.

Confirming the presence of a compact four-planet system orbiting Barnard’s star, similar to other planetary systems orbiting nearby stars, would require many more ESPRESSO observations. These observations would need to be carried out with sufficient cadence to sample the planet periods, and with a sufficiently long baseline to be able to properly model the activity of the star; in particular, those activity signals associated with stellar rotation. This result further stimulates the search for Earth-and sub-Earth-mass planets in the nearest stars of the solar neighbourhood, and encourages new detailed studies with current and future facilities, such as ANDES at ELT (Marconi et al. 2022; Palle et al. 2023).

Data availability

Full Tables 5 and D.3 are available at the CDS via anonymous cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/690/A79

Acknowledgements

J.I.G.H., A.S.M., V.M.P., N.N., A.K.S., C.A.P. and R.R. acknowledge financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) project PID2020-117493GB-I00. J.I.G.H., A.S.M., V.M.P. and R.R. acknowledge financial support from the Government of the Canary Islands project ProID2020010129. A.K.S. acknowledges financial support from La Caixa Foundation (ID 100010434) under the grant LCF/BQ/DI23/11990071. N.N. acknowledges financial support by Light Bridges S.L, Las Palmas de Gran Canaria. F.P.E. and C.L.O. would like to acknowledge the Swiss National Science Foundation (SNSF) for supporting research with ESPRESSO through the SNSF grants nr. 140649, 152721, 166227, 184618 and 215190. The ESPRESSO Instrument Project was partially funded through SNSF’s FLARE Programme for large infrastructures. 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; UIDP/04434/2020. HMT acknowledges support from the “Tecnologías avanzadas para la exploración de universo y sus componentes” (PR47/21 TAU) project funded by Comunidad de Madrid, by the Recovery, Transformation and Resilience Plan from the Spanish State, and by NextGenerationEU from the European Union through the Recovery and Resilience Facility. H.M.T. acknowledges financial support from the Agencia Estatal de Investigación (AEI/10.13039/501100011033) of the Ministerio de Ciencia e Innovación through projects PID2022-137241NB-C41 y PID2022-137241NB-C44. 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). We thank the TNG Director E. Poretti and the TNG staff for the observational effort made to guarantee the continuity of the survey by using HARPS-N during the closure of the Paranal site due to pandemic. This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. The INAF authors acknowledge financial support of the Italian Ministry of Education, University, and Research with PRIN 201278X4FL and the “Progetti Premiali” funding scheme. This work was supported by FCT -Fundação para a Ciência e a Tecnologia through national funds by these grants: UIDB/04434/2020 (DOI: 10.54499/UIDB/04434/2020); UIDP/04434/2020 (DOI: 10.54499/UIDP/04434/2020). Funded/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. 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) A.C.-G. is funded by the Spanish Ministry of Science through MCIN/AEI/10.13039/501100011033 grant PID2019-107061GB-C61. 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. Main analysis performed in Python3 (Van Rossum & Drake 2009) running on Ubuntu (Sobell 2015) systems and Mac OSX systems. Radial velocity computation with S-BART (Silva et al. 2022). Second RV computation with LBL (Artigau et al. 2022). Extensive use of the DACE platform¹³ Extensive usage of Numpy (Harris et al. 2020). Extensive usage of Scipy (Virtanen et al. 2020). All figures built with Matplotlib (Hunter 2007). Periodograms and phase folded curves built using Pyastronomy (Czesla et al. 2019). Gaussian processes modelled with S+LEAF (Delisle et al. 2022), Celerite2 (Foreman-Mackey et al. 2017) and George (Foreman-Mackey et al. 2014). Nested sampling with Dynesty (Speagle 2020).

Appendix A ESPRESSO spectrum of Barnard’s star

The ESPRESSO spectra were reduced using the DRS pipeline v3.0.0. The 156 S1D spectra were combined using the mean values at each pixel after interpolating all spectra to a common wavelength array. The resulting mean spectrum is display in Fig. A.1. The individual ESPRESSO spectra have a signal-to-noise ratio (S/N), from the flux over flux error pixel, of ~ 10, 100, 200, and more than 300 at 4000, 5000, 6000 and 7500 Å, respectively. The mean spectrum has a S/N following approximately Poisson statistics, thus, about 12.5 times the S/N of individual spectra.

Fig. A.1 Mean ESPRESSO HR11 spectrum obtained from the 156 individual S1D_A spectra reduced using the DRS pipeline v3.0.0.

Appendix B Evaluating 3.15 d signal

In order to further investigate the origin of the 3.15 d signal we perform 10 000 simulations from the 124254 samples of analysis of model including activity and planetary signal. We build a simulated RV time series by injecting white noise to the GP model using a random normal distribution with a σ equal to the RV uncertainties including the jitter term (with a median value of 0.61 m s⁻¹). In top panel of Fig. B.1 we show the PSD distributions measured at the 3.15 d and 1.46 d periods of this simulated RV time series (labelled as “simul”). We finally add to this simulated RV time series the three Keplerian signals (with the parameters from model F3 in Table E.1) associated with the three main signals (3.15, 4.12 and 2.34 d) we see above the 1% FAP line in Fig. 3, and we measure the PSD distributions (labelled as “injected”) at 3.15 d and 1.46 d, before and after subtracting the corresponding GP model of each sample (upper and lower panel in Fig. B.1). Since most of the PSD of GP simulated time series is concentrated in strong activity signals at periods larger than 60 d, the PSD at these short periods stays below the 10% FAP line (upper panel). We note that slight PSD increase of the 3.15 d injected signal as compared with the simulated PSD values. We also note that the 1 d alias 1.46 d signal is stronger in the top panel possibly due to the sampling of the time series, as this signal is also stronger in the simulated RV times series that do not contain any planet model. On the other hand, it is also remarkable that the observed signals are stronger that those of the injected signals in the simulated time series, which may be related to maybe too conservative injected white noise (RV uncertainties plus jitter with a median value of 0.61 m s⁻¹), the lower semi-amplitude of the 3.15 d signal in the three planet model and the different GP model not adapted for this three planet model.

We provide in Table B.1 the blind searches of the 3.15 d planet signal using polynomial functions to model the long-term activity signal. We show in this table different runs computed for TM, DRS and LBL datasets.

Fig. B.1 Power spectral density (PSD) distributions of GLS periodograms of 10,000 samples extracted from the run of ESPRESSO FWHM and RV measurements using a model with the GP and a Keplerian orbital model to search for the 3.15 d signal. The PSD distributions are extracted from the GLS periodograms at 3.15 d and the 1 d alias 1.46 d periods from simulated RV time series of the GP values plus white noise before (upper panel) and after (lower panel) subtracting the GP and injecting a planet signal in the simulated data. The observed PSD distributions (lower panel) are measured after subtracted the GP to the original RVs.

Fig. B.2 Activity model of the ESPRESSO RV measurements (upper panels) using a multi-sinusoidal function that uses up to six sinusoids with periods in the range [59,255] d, and RV residuals after subtracting this activity-only model (lower panels). The periods fitted are 79.1, 254.9, 68.2, 89.8, 59.5, and 244.9 d. The GLS periodograms of both the RVs and RV residuals are also displayed (left panels), with the detection of the 3.15 d signal(lower right panel).

Fig. B.3 Activity model of the ESPRESSO RV measurements (upper panels) using moving average with an exponential decay, and RV residuals after subtracting this activity-only model (lower panels). The GLS periodograms of both the RVs and RV residuals are also displayed (left panels), with the detection of the 3.15 d signal(lower right panel).

Fig. B.4 Same as Fig. B.3, but using simulated RV data (upper panels) from the activity-only model with injected white noise and injecting a planet signal at 3.15 d, and RV residuals after subtracting this activity-only model (lower panels). The GLS periodograms of both the RVs and RV residuals are also displayed (left panels), with the detection of the 3.15 d signal(lower right panel).

Appendix C Activity indicators

We measured the activity indexes S_MW, Na I doublet, and Hα from the ESPRESSO spectra, and from the ESPRESSO DRS cross-correlation functions, we measured the full width at half maximum (FWHM), the contrast, and the bisector span (BIS). We modelled the time series of these activity indicators using the GP formalism described in Section 4.1 to evaluate their possible correlations with the time series of RV measurements. Figs. C.1 and C.2 show the model of each index for which we display the index (and its derivative) against the RV. The uncertainties include the jitter. The correlation trend in the index versus RV plots is only drawn in cases where the slope is significant at least at 3 σ (slope / error > 3).

We summarise the analysis in Table C.1, which includes Spearman’s correlation index, p-value (the lower the value, the better the correlation), and the slope we obtain when doing a least-squares fit.

Fig. C.1 Analysis of the ESPRESSO SMW, Hα and Na I spectroscopic indexes, a, b: SMW, Hα and Na I-index time-series with the best-model fit. The data is split into two panels because of a large period with no observations between the two campaigns. The shaded area shows the variance of the GP model, c: GLS periodogram of the SMW, Hα and Na I-index data. The red vertical dashed line shows the most significant period, d: Relationship between the SMW, Hα and Na I-index data. The best fit is shown when the slope is ≥3σ different from zero, e, f: Residuals of the SMW, Hα and Na I-index after subtracting the best model fit. g: GLS periodogram of the residuals, h: Comparison of the CCF RV and gradient of the SMW, Hα and Na I-index model.

Fig. C.2 Analysis of the BIS and contrast of the ESPRESSO CCF. a, b: BIS and contrast time-series with the best-model fit. The data is split into two panels because of a large period with no observations between the two campaigns. The shaded area shows the variance of the GP model, c: GLS periodogram of the CCF BIS and contrast data. The red vertical dashed line shows the most significant period, d: Relationship between the CCF RV and CCF BIS and contrast data. The best fit is shown when the slope is ≥3σ different from zero, e, f: Residuals of the CCF BIS and contrast after subtracting the best model fit. g: GLS periodogram of the residuals, h: Comparison of the CCF RV and gradient of the CCF BIS and contrast model.

Appendix D Analysis of all datasets

In Fig. D.3 we show the FWHM and RV measurements, including ESPRESSO, HARPS and HARPS-N data, with the analysis of the GP+cycle model together with one Keplerian model using wide priors (model J1 in Table 3). Figure D.1 depicts the RV curve in orbital phase of the sub-Earth mass planet at 3.15 d of model represented in Fig. D.3. Figure D.2 displays the corner plot of the orbit of the 3.15 d sub-Earth mass planet. In Fig. D.4 we show the prior and posteriors distributions of all parameters of the model in Fig. D.3. Finally, in Fig. D.5 we represent the FWHM and RV measurements, including ESPRESSO, HARPS, HARPS-N and CARMENES data, with the analysis of the GP+cycle model together with one Keplerian model using wide priors (model L1 in Table 3). We also provide the prior, posterior and parameter values of this model L1 in Table D.1 and the model statistics in Table D.2. In Table D.3 we provide RV and FWHM measurements of all datasets used in Fig. D.5.

Fig. D.1 RV curve of the sub-Earth-mass planet of GJ 699 with a 3.15 d orbital period together with ESPRESSO, HARPS and HARPS-N RVs with uncertainties with (light colour) and without including the jitter term (dark colour) coming from the global model.

Fig. D.2 Corner plot with the posterior distribution of the orbital parameters of the sub-Earth-mass planet of GJ 699 with a 3.15 d from the model using ESPRESSO, HARPS and HARPS-N data.

Fig. D.3 FWHM measurements (top), RV measurements (middle), RV residuals from cycle and GP model (next bottom) and RV residuals from Keplerian model (bottom), of the ESPRESSO, HARPS and HARPS-N datasets using SHO (PROT and PROT/2) GP model (grey solid line) to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the long-term cycle (red thick dashed line), and a Keplerian model (grey shaded area), and GLS periodograms (left) of GJ 699 (model J1 in Table 3). The uncertainties include the jitter term coming from the global model.

Fig. D.4 Prior and posterior distributions of all 40 parameters of model J1 in Table 3 of ESPRESSO, HARPS and HARPS-N data together, including offsets, jitter, the parameters of GP and the long-term cycle, and the parameters of the Keplerian model revealing the sub-Earth-mass planet of GJ 699 with a 3.15 d orbital period.

Fig. D.5 FWHM measurements (top), RV measurements (middle), RV residuals from cycle and GP model (next bottom) and RV residuals from Keplerian model (bottom), of the ESPRESSO, HARPS, HARPS-N and CARMENES datasets using SHO (PROT and PROT/2) GP model (grey solid line) to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the long-term cycle (red thick dashed line), and a Keplerian model (grey shaded area), and GLS periodograms (left) of GJ 699 (model L1 in Table 3 with prior and posterior parameters given in Table D.1). The uncertainties include the jitter term coming from the global model.

Appendix E Candidate multi-planet system

The fact that the candidate signals are so weak makes it difficult to detect them in a blind search, so we decided to run guided search with normal priors around the period of the signals with a σ = 0.3 d (prior (P_orb, 0.3) d). We first search for the highest peak in the GLS periodogram after subtraction of the activity model, which is the 3.15 d signal. We adopt a prior on orbital period P_orb (3.15, 0.3) d and eccentricity e with and (0, 0.05), basically an almost circular orbit, based on previous results. This is model F1 in Table E.1. The Bayesian evidence, though, may be not that informative as compared with the blind search runs because the prior on period is very narrow. However, it indicates that the global fit model reduces significantly the RMS of the RV residuals from 0.60 to 0.45 m s⁻¹. The next peak in the RV residuals is 4.12 d, so we then run a two-planet model F2 with priors centered at period 3.15 d and 4.12 d, with an increase in Δ ln = +3.9 with respect to the single-planet model F1. In the RV residuals appears the peak at 2.34 d significantly above the 0.1% FAP line.

We run a three-planet model for periods 3.15 d, 4.12 d and 2.34 d, again with narrow priors on eccentricity. The resulting run F3 has a higher Bayesian evidence by Δ ln = +3.9 over the two-planet model F2. We also check the possibility of a wider prior on eccentricity (e with and (0, 0.3)) of the three-planet model. This model tend to favor very low eccentricity (consistent with zero) for the signals 3.15 d and 4.12 d, but it goes to high eccentricity values for the signal at 2.34 d. In this case, the RV residuals do not leave any peak above the 0.1% FAP line. However, such a system with short-period candidate planets with the inner one being a high eccentric planet would not be stable in the long term. In the case of low priors on eccentricity, model F3 reveals a peak at 6.74 d in the GLS periodogram of the RV residuals.

We then run a four-planet model, with a narrow prior on eccentricity e with and (0, 0.05), for the fourth signals with narrow priors centered at 3.15 d, 4.12 d, 2.34 d and 6.74 d. This model F4 has the same Bayesian evidence as F3, only increasing by Δ ln = +0.5. We checked that using this narrow prior on eccentricity gives almost the same result as using sinusoids as circular orbits (models F1c and F4c in Table E.1), for both models F1 and F4. Fig. E.1 summarises the result of run F4, where we display just the RVs after subtracting the activity model displayed in Fig. E.2. Although F4 is a simultaneous fit of a four-planet Keplerian model, we represent in Fig. E.1 the RVs after subtracting first the activity model, and then in each row, after removing the signal at the highest peak in the GLS periodograms given in middle panels, one by one. The left panels show the RVs in phase of the given planet period after subtracting the activity model and the rest of the fitted planetary signals. In the middle panels of Fig. E.1, the GLS periodograms of the RVs (displayed in left panels) show that after subtracting the last candidate signal at 6.74 d, there are no more signals above the 10% FAP line. We also see, in the central panels of Fig. E.2, several peaks corresponding to the main signals and their 1 d alias, 3.15 d (with its 1 d alias 1.46 d and 0.76 d, and with 0.59 d as 1 d alias of 1.46 d), 4.12 d (with its 1 d alias 1.32 d and 0.80 d), 2.34 d (with its 1 d alias 1.74 d and 0.70 d), and 6.74 d (with its 1 d alias 1.17 d and 0.87 d). The RMS of the ESPRESSO RVs goes from 0.61 m s⁻¹ in the RVs after subtracting the activity cycle+rotation model to 0.49, 0.39, 0.29, and 0.24 m s⁻¹ after removing the candidate planet signals at 3.15, 4.12, 2.34, and 6.74 d (see Figs. E.1 and E.2). The semi-amplitude velocities of run F4 are k_p = 47, 41, 35 and 20 cm s⁻¹, thus providing the candidate planetary minimum masses of 0.32, 0.31, 0.22, 0.17 M_⊕ for candidate planets tentatively labelled as b, c, d and e (see Table E.2).

Finally, we also tested a run with the first confirmed planet as a circular orbit with a narrow prior (3.15, 0.3) d on orbital period and three additional planetary signals also as circular orbits but uniform priors (2, 7) d on orbital period. We use for all candidate planetary signals an uniform prior (0, 5) m s⁻¹ on semi-amplitude velocity. We recover the same solution as in models F4 and F4c in Table E.1, but the Bayesian evidence of this model is just the same as for model F1c in Table E.1, i.e. ln = −438.0. Figure E.3 shows the prior and posterior distributions of the orbital period and semi-amplitude velocities of the four candidate planets. The orbital periods are unequivocally retrieved for all planets, in particular for the first three planets, and for the fourth planet with just a few samples falling on both 1-yr alias (6.61 d and 6.86 d) of the main period 6.74 d. The semi-amplitudes are also unequivocally recovered for the first three planets at 9.2, 8.2 and 7.2 σ, whereas for the fourth planet it is only recovered at 2.4 σ.

E.1 Stability of the candidate multi-planet system

We evaluate the stability of the 4-planet candidate system using the SPOCK¹⁴ tool (Tamayo et al. 2020), together with the Nbody Regressor (Hussain & Tamayo 2020) and the REBOUND¹⁵ code (e.g. Rein & Tamayo 2017). We found that assuming the zero eccentricity¹⁶ for all planets, the system remains stable more than 10⁹ orbits of the inner candidate planet at 2.34 d, which is the maximum Nbody instability time explored by default by this particular Nbody regressor. We run this calculation using the planet masses assuming m_p = m_p sin i (thus i = 90 degrees), but also we run it again assuming by m_p = (m_p sin i) × 3. This latter conservative case would mean an orbital inclination of i = 19.5 degrees. In both cases the result is the same.

The candidate system architecture is displayed using the REBOUND plotting tool in Fig. E.5. We see all these candidatem planets in orbits inner to the habitable zone limits estimated in Section 3. The candidate four-planet system looks particularly compact, similar to other systems such as e.g. the TRAPPIST seven-planet system (Gillon et al. 2017). Dreizler et al. (2024) have recently revisited the planetary system around the Teegarden’s star with a detection and confirmation of one additional Earth-mass planet in a 26 d orbit, and additional candidate planets. These authors compare the Teegarden’s system with other compact systems with their planets orbiting within 0.1 AU, similar to the candidate planets of Barnard’s star. They highlight mutual Hill radius separation as a good indicator of the dynamic compactness of a planetary system, also called fractional orbital separation (Gladman 1993). The mutual Hill radius separation is calculated from the stellar mass, M_*, the planetary masses, m_p, and semi-major axes of two neighbouring planets j and j+1 as in eq. E.1: (E.1)

Dreizler et al. (2024) compare the values of Δ(R_H) of planets in Teegarden’s star with other six compact systems. Although the sample is small, they argue that planet systems with more massive planets have mutual Hill radius separations above Δ(R_H) ≥ 30 while the ones with low-mass, likely terrestrial, planets have Δ(R_H) ≤ 20. Weiss et al. (2018) indicated that systems with more planets tend to have smaller Δ(R_H) values, but only about 10% of the systems have a Δ(R_H) < 10. Simulations of the stability evolution of compact multi-planet systems suggest a log-linear dependency with the mutual Hill radius separation (e.g. Chambers et al. 1996). Chambers et al. (1996) found that systems with mutual Hill radius separation Δ(R_H) less than 10 are always unstable. These simulations suggest that at about 13 mutual Hill radius separations, the orbital crossing timescale seems to be longer than at least 10⁹ orbits for all systems (Gratia & Lissauer 2021; Rice & Steffen 2023). Below this limit the orbital crossing timescale changes by more than one order of magnitude for small changes in mutual Hill separation. The seven-planet system TRAPPIST-1 with planet masses from 0.33 to 1.37 M_⊕ at orbital periods from 1.5 d to 18.8 d (Agol et al. 2021) have all values Δ(R_H) below 13. Similarly, the three-planet system YZCeti with planets at P_orb = 2.02, 3.06 and 4.66 d with minimum masses 0.70, 1.14, and 1.09 M_⊕ respectively (Astudillo-Defru et al. 2017b; Stock et al. 2020) have also all mutual Hill radius separations in the range 10 < Δ(R_H) < 13 (see Fig. 6 in Dreizler et al. 2024). The three-planet system orbiting Teegarden’s star at periods 4.9 d, 11.4 d, and 26.1 d with minimum masses 1.16, 1.05 and 0.82 M_⊕ have Δ(R_H) = 19.3 and 19.7, for planet pairs bc and cd, respectively. They also suggest additional candidate planets at 7.7 d and 17 d, which would move the mutual Hill radius separations of all planet pairs of the candidate five-planet system at Δ(R_H) < 13 (Dreizler et al. 2024). On the other hand, the two-planet systems GJ 1002 with planets at 10.3 d and 21.2 d with minimum masses 1.1 and 1.4 M_⊕ (Suárez Mascareño et al. 2023) and Proxima Cen system with the planets Proxima d and Proxima b at periods 5.1 d and 11.2 d with minimum masses 0.26 M_⊕ and 1.1 M_⊕ (Faria et al. 2022), have mutual Hill radius separations larger than 13, specifically at Δ(R_H) = 17.2 and 23.0 for planet pairs bc in GJ 1002 and db in Proxima Centauri. The sub-Earth mass candidate planets in Barnard’s star system have Δ(R_H) of 13.3, 11.3 and 22.6 for planet pairs db, bc and ce, respectively. The two inner and outer pairs of planets db and ce have a Δ(R_H) > 13, whereas the middle pair of planets bc would have Δ(R_H) < 13 but larger than 10. Whether this could cause orbital instabilities to the candidate four-planet system would require further investigation.

Fig. E.1 ESPRESSO RV measurements versus BJD after removing the long-term cycle and the GP model (Left), GLS periodograms (middle), and RV curves of the sub-Earth-mass planet candidates (right) of GJ 699 with 3.15 d, 4.12 d, 2.34 d and 6.74 d orbital periods together with ESPRESSO RVs with uncertainties with (light colour) and without including the jitter term (dark colour) coming from the global model F4 in Table E.1.

Fig. E.2 FWHM measurements (top), RV measurements (middle), RV residuals from cycle and GP model (next bottom) and RV residuals from four Keplerian model (bottom), of the ESPRESSO, HARPS, HARPS-N and CARMENES datasets using SHO (PROT and PROT/2) GP model to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the long-term cycle, and four Keplerian models, and GLS periodograms (left) of GJ 699 (model F4 in Table E.1). The uncertainties include the jitter term coming from the global model.

Fig. E.3 Prior and posterior distributions of the period and semi-amplitude from the global analysis of FWHM and RV measurements of the blind search of the candidate planet system using ESPRESSO data. We use a normal prior on orbital period (3.15, 0.3) d for the confirmed planet b (left), and an uniform prior (2, 7) d for planet c (center-left), planet d (center-right), and planet e (right). For all planets we use an uniform prior (0, 5) m s−1 on semi-amplitude velocity.

The candidate planetary system orbiting our second closest neighbour Barnard’s star consists of four sub-Earth mass planets at periods 2.34 d, 3.15 d, 4.12 d, and 6.74 d with minimum masses 0.22, 0.32, 0.31, 0.17 M_⊕. Fig. E.4 compares the confirmed planet Barnard b and the candidate planets with other planets in the exoplanet database.

Fig. E.4 Same as Fig. 17 with the confirmed planet Barnard b (red), and planet candidates (purple) of the four-planet candidate system orbiting Barnard’s star, together with the two planets orbiting Proxima Cen, Proxima b (yellow) and d (blue) are highlighted.

Fig. E.5 Schematic view of the candidate planetary system of GJ 699 together with the habitable zone (light green region). The planets are depicted as green (planet b) circle, blue (planet c), red (planet d), and yellow (planet e).

References

All Tables

All Figures

	Fig. 1 ESPRESSO RV measurements (right) and GLS periodograms (left) of GJ 699 after subtracting the median of each dataset before (E18) and after (E19) the intervention in June 2019. Also shown are RV measurements from the ESPRESSO Data Reduction Software (DRS; top), from the S-BART template matching (TM) code (middle), and from the line-by-line LBL code (bottom).
In the text

	Fig. 2 Radial-velocity measurements (right) and GLS periodograms (left) of GJ 699 of ESPRESSO (top), HARPS and HARPS-N (middle), and CARMENES (bottom).
In the text

	Fig. 3 ESPRESSO FWHM measurements (top), RV measurements (middle), and RV residuals (bottom) from the SHO (PROT and PROT/2) GP model, and GLS periodograms (left) of GJ 699. The uncertainties include the jitter term coming from the global model A in Table 3.
In the text

Fig. 4 Analysis of the FWHM of the ESPRESSO CCF. Panels a and b: FWHM time-series with the best-model fit. The data is split into two panels because of a large period with no observations between the two campaigns. The shaded area shows the variance of the GP model. Panel c: GLS periodogram of the CCF FWHM data. The red vertical dashed line shows the most significant period. Panel d: relationship between the CCF RV and CCF FWHM data. The best fit is shown when the slope is ≥3σ different from zero. Panels e and f: residuals of the CCF FWHM after subtracting the best model fit. Panel g: GLS periodogram of the residuals. Panel h: comparison of the CCF RV and gradient of the CCF FWHM model.

In the text

	Fig. 5 ASAS-SN photometry in the same temporal baseline as the RV data together with the double sinusoidal models of the long-term cycle and rotation signals (top left) and the GLS periodogram before and after subtracting the long-term cycle (top right). The posterior and prior distributions of these models including the long-term cycle (bottom left) and rotation periods (bottom right).
In the text

	Fig. 6 ESPRESSO RV and FWMH measurements, and échelle temperature sensor and vacuum vessel pressure sensor measurements versus BJD and barycentric RV correction. The GLS periodograms (left panels) show the result before (red line) and after (black line) subtracting the fitted function (dashed red line) in the left panels.
In the text

	Fig. 7 Prior and posterior distributions of the long-term cycle, the GP rotation period, and the timescale from the global analysis of FWHM and RV measurements of ESPRESSO, HARPS, and HARPS-N (model J1 in Table 3).
In the text

	Fig. 8 GLS periodograms of ESPRESSO CCF FWHM measurements (top), TM RV measurements (middle), and the window function (bottom) of GJ 699 after subtracting the median of each dataset before (E18) and after (E19) the intervention in June 2019.
In the text

	Fig. 9 Evolution of the false-alarm probability of the 3.15 d and 1.46 d (1 day alias) signals (top) and the GP rotation period (bottom) with the number of observations. The 0.1, 1, and 10% FAP lines are computed using the full 149 ESPRESSO dataset.
In the text

	Fig. 10 GLS power spectral density of the observed signals 3.15d and 1.46d after removing the GP compared to a simulated 54 cm s−1 injected sinusoidal signal at different orbital periods.
In the text

Fig. 11 Prior and posterior distributions of the period and semi-amplitude from the global analysis of FWHM and RV measurements of the blind search of the 233d candidate planet using ESPRESSO, HARPS, HARPS-N and CARMENES data (left, model O1 in Table 4), a blind search (centre-left, model N2 in Table 4) and a guided search (centre-right, model N1 in Table 4) of the 233d candidate planet using ESPRESSO only, and a blind search of the 233d simulated signal (right, model S1 in Table 4).

In the text

Fig. 12 ESPRESSO FWHM measurements (top), RV measurements (middle), RV residuals (next bottom) from long-term cycle (red thick dashed line) and SHO (PROT and PROT/2) GP model (grey solid line), and RV residuals (bottom) from the Keplerian model (grey shaded area), and GLS periodograms (left) of GJ 699. The uncertainties include the jitter term coming from the global model E1e in Table 3, with prior and posterior parameters given in Table 6.

In the text

	Fig. 13 Corner plot with the posterior distributions of the orbital parameters of the sub-Earth-mass planet of GJ 699 with a 3.15 d orbital period, with prior and posterior parameters given in Table 6.
In the text

	Fig. 14 RV curve of the sub-Earth-mass planet of GJ 699 with a 3.15 d orbital period together with ESPRESSO RVs, with uncertainties with (light colour) and without the jitter term (dark colour) coming from the global model E1e, with prior and posterior parameters given in Table 6.
In the text

	Fig. 15 FIP periodogram of the 4-planet model of the ESPRESSO data of GJ 699. The periods of the detected signals are indicated in red solid circles. The − log10 FIP as a function of the center of the period is represented as a blue line. The green dashed and dotted orange lines show the 1 and 10% FIP thresholds, respectively.
In the text

Fig. 16 Detection limits. Top panel: semi-amplitude kp versus period Porb for our detection limits based on injection and recovery tests. The white region shows the signals that we would detect in a periodogram with 1 % FAP. The blue region shows signals we would detect with a FAP lower than 1 % FAP. With the current dataset and activity model, we would miss the red region signals. The top grey dashed line shows the maximum semi-amplitudes that would still be compatible with the data. The bottom grey dashed line shows the minimum amplitude that we expect to detect with ESPRESSO data and the current exposure times. The green solid circle shows the confirmed planet and the orange circles show the position of the candidate planets. Bottom panel: same as top panel but in mass versus period.

In the text

Fig. 17 Minimum mass vs. orbital period diagram for known planets from the NASA Exoplanet archive as of December 2023 orbiting solar-type stars, together with those discovered and confirmed using ESPRESSO (green circles). Confirmed planet Barnard b (red), and planet candidates (purple) from the 4-planet candidate system orbiting Barnard’s star, together with the two planets orbiting Proxima Cen, Proxima b (yellow) and d (blue) are highlighted. Planets of the solar system (grey circles) are also labelled. Inclined solid and dashed lines show the RV semi- amplitude of planets orbiting a late M dwarf star with 0.25 M⊙ and a G dwarf star with 1 M⊙ star assuming a RV semi-amplitude of 1 cm s−1 (green line) and 10 cm s−1 (blue line), respectively, and null eccentricity.

In the text

Fig. 18 Hipparcos-Gaia PMA sensitivity to companions of a given mass (in M⊕) as a function of the orbital semi-major axis (in AU) orbiting Barnard’s star. The solid black line identifies the combinations of mass and separation explaining the observed proper-motion anomaly (PMA) at the mean epoch of Gaia DR3 (Eq. (15) in Kervella et al. 2019). The shaded light-blue region corresponds to the 1σ uncertainty region.

In the text

	Fig. A.1 Mean ESPRESSO HR11 spectrum obtained from the 156 individual S1D_A spectra reduced using the DRS pipeline v3.0.0.
In the text

Fig. B.1 Power spectral density (PSD) distributions of GLS periodograms of 10,000 samples extracted from the run of ESPRESSO FWHM and RV measurements using a model with the GP and a Keplerian orbital model to search for the 3.15 d signal. The PSD distributions are extracted from the GLS periodograms at 3.15 d and the 1 d alias 1.46 d periods from simulated RV time series of the GP values plus white noise before (upper panel) and after (lower panel) subtracting the GP and injecting a planet signal in the simulated data. The observed PSD distributions (lower panel) are measured after subtracted the GP to the original RVs.

In the text

Fig. B.2 Activity model of the ESPRESSO RV measurements (upper panels) using a multi-sinusoidal function that uses up to six sinusoids with periods in the range [59,255] d, and RV residuals after subtracting this activity-only model (lower panels). The periods fitted are 79.1, 254.9, 68.2, 89.8, 59.5, and 244.9 d. The GLS periodograms of both the RVs and RV residuals are also displayed (left panels), with the detection of the 3.15 d signal(lower right panel).

In the text

	Fig. B.3 Activity model of the ESPRESSO RV measurements (upper panels) using moving average with an exponential decay, and RV residuals after subtracting this activity-only model (lower panels). The GLS periodograms of both the RVs and RV residuals are also displayed (left panels), with the detection of the 3.15 d signal(lower right panel).
In the text

	Fig. B.4 Same as Fig. B.3, but using simulated RV data (upper panels) from the activity-only model with injected white noise and injecting a planet signal at 3.15 d, and RV residuals after subtracting this activity-only model (lower panels). The GLS periodograms of both the RVs and RV residuals are also displayed (left panels), with the detection of the 3.15 d signal(lower right panel).
In the text

Fig. C.1 Analysis of the ESPRESSO SMW, Hα and Na I spectroscopic indexes, a, b: SMW, Hα and Na I-index time-series with the best-model fit. The data is split into two panels because of a large period with no observations between the two campaigns. The shaded area shows the variance of the GP model, c: GLS periodogram of the SMW, Hα and Na I-index data. The red vertical dashed line shows the most significant period, d: Relationship between the SMW, Hα and Na I-index data. The best fit is shown when the slope is ≥3σ different from zero, e, f: Residuals of the SMW, Hα and Na I-index after subtracting the best model fit. g: GLS periodogram of the residuals, h: Comparison of the CCF RV and gradient of the SMW, Hα and Na I-index model.

In the text

Fig. C.2 Analysis of the BIS and contrast of the ESPRESSO CCF. a, b: BIS and contrast time-series with the best-model fit. The data is split into two panels because of a large period with no observations between the two campaigns. The shaded area shows the variance of the GP model, c: GLS periodogram of the CCF BIS and contrast data. The red vertical dashed line shows the most significant period, d: Relationship between the CCF RV and CCF BIS and contrast data. The best fit is shown when the slope is ≥3σ different from zero, e, f: Residuals of the CCF BIS and contrast after subtracting the best model fit. g: GLS periodogram of the residuals, h: Comparison of the CCF RV and gradient of the CCF BIS and contrast model.

In the text

	Fig. D.1 RV curve of the sub-Earth-mass planet of GJ 699 with a 3.15 d orbital period together with ESPRESSO, HARPS and HARPS-N RVs with uncertainties with (light colour) and without including the jitter term (dark colour) coming from the global model.
In the text

	Fig. D.2 Corner plot with the posterior distribution of the orbital parameters of the sub-Earth-mass planet of GJ 699 with a 3.15 d from the model using ESPRESSO, HARPS and HARPS-N data.
In the text

Fig. D.3 FWHM measurements (top), RV measurements (middle), RV residuals from cycle and GP model (next bottom) and RV residuals from Keplerian model (bottom), of the ESPRESSO, HARPS and HARPS-N datasets using SHO (PROT and PROT/2) GP model (grey solid line) to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the long-term cycle (red thick dashed line), and a Keplerian model (grey shaded area), and GLS periodograms (left) of GJ 699 (model J1 in Table 3). The uncertainties include the jitter term coming from the global model.

In the text

	Fig. D.4 Prior and posterior distributions of all 40 parameters of model J1 in Table 3 of ESPRESSO, HARPS and HARPS-N data together, including offsets, jitter, the parameters of GP and the long-term cycle, and the parameters of the Keplerian model revealing the sub-Earth-mass planet of GJ 699 with a 3.15 d orbital period.
In the text

Fig. D.5 FWHM measurements (top), RV measurements (middle), RV residuals from cycle and GP model (next bottom) and RV residuals from Keplerian model (bottom), of the ESPRESSO, HARPS, HARPS-N and CARMENES datasets using SHO (PROT and PROT/2) GP model (grey solid line) to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the long-term cycle (red thick dashed line), and a Keplerian model (grey shaded area), and GLS periodograms (left) of GJ 699 (model L1 in Table 3 with prior and posterior parameters given in Table D.1). The uncertainties include the jitter term coming from the global model.

In the text

Fig. E.1 ESPRESSO RV measurements versus BJD after removing the long-term cycle and the GP model (Left), GLS periodograms (middle), and RV curves of the sub-Earth-mass planet candidates (right) of GJ 699 with 3.15 d, 4.12 d, 2.34 d and 6.74 d orbital periods together with ESPRESSO RVs with uncertainties with (light colour) and without including the jitter term (dark colour) coming from the global model F4 in Table E.1.

In the text

Fig. E.2 FWHM measurements (top), RV measurements (middle), RV residuals from cycle and GP model (next bottom) and RV residuals from four Keplerian model (bottom), of the ESPRESSO, HARPS, HARPS-N and CARMENES datasets using SHO (PROT and PROT/2) GP model to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the long-term cycle, and four Keplerian models, and GLS periodograms (left) of GJ 699 (model F4 in Table E.1). The uncertainties include the jitter term coming from the global model.

In the text

Fig. E.3 Prior and posterior distributions of the period and semi-amplitude from the global analysis of FWHM and RV measurements of the blind search of the candidate planet system using ESPRESSO data. We use a normal prior on orbital period (3.15, 0.3) d for the confirmed planet b (left), and an uniform prior (2, 7) d for planet c (center-left), planet d (center-right), and planet e (right). For all planets we use an uniform prior (0, 5) m s−1 on semi-amplitude velocity.

In the text

	Fig. E.4 Same as Fig. 17 with the confirmed planet Barnard b (red), and planet candidates (purple) of the four-planet candidate system orbiting Barnard’s star, together with the two planets orbiting Proxima Cen, Proxima b (yellow) and d (blue) are highlighted.
In the text

	Fig. E.5 Schematic view of the candidate planetary system of GJ 699 together with the habitable zone (light green region). The planets are depicted as green (planet b) circle, blue (planet c), red (planet d), and yellow (planet e).
In the text