Issue 
A&A
Volume 690, October 2024



Article Number  A79  
Number of page(s)  36  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202451311  
Published online  01 October 2024 
A subEarthmass planet orbiting Barnard’s star^{★}
^{1}
Instituto de Astrofísica de Canarias,
38205
La Laguna, Tenerife,
Spain
^{2}
Departamento de Astrofísica, Universidad de La Laguna,
38206
La Laguna, Tenerife,
Spain
^{3}
Consejo Superior de Investigaciones Científicas,
Spain
^{4}
Observatoire de Genève, Département d’Astronomie, Université de Genève,
Chemin Pegasi 51b,
1290
Versoix,
Switzerland
^{5}
Instituto de Astrofísica e Ciências do Espaco, CAUP, Universidade do Porto, Rua das Estrelas,
4150762
Porto,
Portugal
^{6}
Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre,
4169007
Porto,
Portugal
^{7}
Centro de Astrofísica da Universidade do Porto, Rua das Estrelas,
4150762
Porto,
Portugal
^{8}
Centro de Astrobiología (CAB), CSICINTA, ESAC campus,
Camino Bajo del Castillo s/n,
28692
Villanueva de la Cañada (Madrid),
Spain
^{9}
INAF  Osservatorio Astronomico di Trieste,
via G. B. Tiepolo 11,
34143
Trieste,
Italy
^{10}
IFPU–Institute for Fundamental Physics of the Universe,
via Beirut 2,
34151
Trieste,
Italy
^{11}
INAF  Osservatorio Astrofisico di Torino,
Strada Osservatorio 20,
10025,
Pino Torinese (TO),
Italy
^{12}
ESO  European Southern Observatory,
Av. Alonso de Cordova 3107, Vitacura,
Santiago,
Chile
^{13}
Département de Physique, Institut Trottier de Recherche sur les Exoplanètes, Université de Montréal, Montréal,
Québec,
H3T 1J4,
Canada
^{14}
INAF  Osservatorio Astronomico di Palermo,
Piazza del Parlamento 1,
90134
Palermo,
Italy
^{15}
Instituto de Astrofísica e Ciências do EspaÇo, Faculdade de Ciências da Universidade de Lisboa, Campo Grande,
1749016
Lisboa,
Portugal
^{16}
Light Bridges S. L.,
Avda. Alcalde Ramírez Bethencourt, 17,
35004
Las Palmas de Gran Canaria, Canarias,
Spain
^{17}
Departamento de Física de la Tierra y Astrofísica & IPARCOSUCM (Instituto de Física de Partículas y del Cosmos de la UCM), Facultad de Ciencias Físicas, Universidad Complutense de Madrid,
28040
Madrid,
Spain
^{18}
Center for Space and Habitability, University of Bern,
Gesellschaftsstrasse 6,
3012
Bern,
Switzerland
^{19}
Weltraumforschung und Planetologie, Physikalisches Institut, University of Bern,
Gesellschaftsstrasse 6,
3012
Bern,
Switzerland
^{20}
Subaru Telescope, National Astronomical Observatory of Japan,
650 N Aohoku Place,
Hilo,
HI
96720,
USA
^{21}
Hamburger Sternwarte,
Gojenbergsweg 112,
21029
Hamburg,
Germany
^{★★} Corresponding author; email: jonay.gonzalez@iac.es
Received:
29
June
2024
Accepted:
26
August
2024
Context. ESPRESSO guaranteed time observations (GTOs) at the 8.2m VLT telescope were performed to look for Earthlike exoplanets in the habitable zone of nearby stars. Barnard’s star is a primary target within the ESPRESSO GTO as it is the second closest neighbour to our Sun after the α Centauri stellar system.
Aims. We present here a large set of 156 ESPRESSO observations of Barnard’s star carried out over four years with the goal of exploring periods of shorter than 50 days, thus including the habitable zone (HZ).
Methods. Our analysis of ESPRESSO data using Gaussian process (GP) to model stellar activity suggests a longterm activity cycle at 3200 d and confirms stellar activity due to rotation at 140 d as the dominant source of radial velocity (RV) variations. These results are in agreement with findings based on publicly available HARPS, HARPSN, and CARMENES data. ESPRESSO RVs do not support the existence of the previously reported candidate planet at 233 d.
Results. After subtracting the GP model, ESPRESSO RVs reveal several shortperiod candidate planet signals at periods of 3.15 d, 4.12 d, 2.34 d, and 6.74 d. We confirm the 3.15 d signal as a subEarth mass planet, with a semiamplitude of 55 ± 7 cm s^{−1}, leading to a planet minimum mass m_{p} sin i of 0.37 ± 0.05 M_{⊕}, which is about three times the mass of Mars. ESPRESSO RVs suggest the possible existence of a candidate system with four subEarth mass planets in circular orbits with semiamplitudes from 20 to 47 cm s^{−1}, thus corresponding to minimum masses in the range of 0.17–0.32 M_{⊕}.
Conclusions. The subEarth mass planet at 3.1533 ± 0.0006 d is in a closeto circular orbit with a semimajor axis of 0.0229 ± 0.0003 AU, thus located inwards from the HZ of Barnard’s star, with an equilibrium temperature of 400 K. Additional ESPRESSO observations would be required to confirm that the other three candidate signals originate from a compact shortperiod planet system orbiting Barnard’s star inwards from its HZ.
Key words: techniques: radial velocities / techniques: spectroscopic / planets and satellites: terrestrial planets / stars: activity / stars: lowmass / stars: individual: GJ 699
Based [in part] on Guaranteed Time Observations collected at the European Southern Observatory under ESO programmes 1102.C0744, 1104.C0350, 106.21M2.001, 106.21M2.004, 106.21M2.006, 108.22GM.001, 108.2254.001, 108.2254.003, 108.2254.004, 108.2254.006, 110.24CD.001, 110.24CD.003 by the ESPRESSO Consortium.
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
The field of exoplanet science has been evolving very quickly in recent years towards the detection and characterisation of Earthlike 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 groundbased highresolution ultrastable spectrographs, such as HARPS (Mayor et al. 2003), HARPSN (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 Earthlike planets (e.g. Gillon et al. 2017; Dittmann et al. 2017; Zechmeister et al. 2019; LilloBox 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 Earthmass planet Proxima Centauri b (AngladaEscudé et al. 2016) orbiting the closest star to our Sun has propelled exoplanet studies to focus the search on Earthlike planets in the habitable zone around stars of the solar neighbourhood. These temperate Earthlike 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 Earthmass exoplanets quickly shifted to the continuous monitoring of M dwarfs, with the development of new instruments in the nearinfrared, 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 Sunlike stars, and have habitable zones closer to their host star, making them ideal targets for blind RV searches of Earthlike 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^{−1} RV precision. ESPRESSO has confirmed, for instance, the Earthmass planet Proxima b and discovered the subEarth Proxima d, a 0.26 M_{⊕} planet (approximately twice the mass of Mars) orbiting Proxima Centauri, from the measurement of a small RV semiamplitude of 39 ± 7 cm s^{−1} (Suárez Mascareño et al. 2020; Faria et al. 2022). ESPRESSO is opening a new frontier at subm/s precision, making it possible to discover and characterise Earth and subEarthmass and subEarth size exoplanets in the solar neighbourhood. ESPRESSO has detected, for instance, the 0.4 M_{⊕} planet L9859 b (half of the mass of Venus) orbiting an M3V star (Demangeon et al. 2021), and one superEarth and two superMercuries 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 Xray flux of GJ699 with the Chandra satellite in the energy range 0.310keV at F_{x} ~ 4.8 × 10^{−14} erg cm^{−2} s^{−1} (log_{10}(L_{X}[erg s^{−1}]) = 25.3; L_{X}/L_{bol} = 1.6 × 10^{−6}). The Xray luminosity is within a factor of two of previous ROSAT data (log_{10}(L_{X}[erg s^{−1}]) = 25.6). This low Xray luminosity with log_{10}(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; AstudilloDefru et al. 2017a; ToledoPadró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). ToledoPadrón et al. (2019) reported a rotation period of P_{ROT} = 145 ± 15 d from the timeseries analysis of spectroscopic activity indexes, and also found evidence of a longterm activity cycle of Barnard’s star from a time series of the CaHK index and ASASSN 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_{⊕} superEarthlike 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 (ToledoPadró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 semiamplitude of K_{p} = 1.2 m s^{−1} is disfavoured when compared to a flat model.
Here we present the ESPRESSO observations of Barnard’s star (GJ 699), showing a subm/s precision that reveals the presence of a shortperiod subEarth planet and three additional subEarth planet candidates. ESPRESSO data allow us to also evaluate the presence of the superEarth 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^{1} instrument (Pepe et al. 2021). The ESPRESSO project has mainly been dedicated to the search for and characterisation of exoplanets (e.g. LilloBox et al. 2021; Faria et al. 2022; Suárez Mascareño et al. 2023; Lavie et al. 2023; CastroGonzá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 Earthlike planets within its HZ.
This work has also made use of public HARPS, HARPSN, and CARMENES data, with some of the HARPS and HARPSN spectra taken by consortium members as part of the followup 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 fibrefed, crossdispersed, highresolution échelle spectrograph located in the Combined Coudé Laboratory (CCL) at the incoherent focus, where the FrontEnd unit can combine the light from up to four Unit Telescopes (UT) of the Very Large Telescope (VLT) at Paranal Observatory (ESO, Chile). The socalled Coudé train optical system feeds the light of each UT to the spectrograph. The FrontEnd 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.2m UTs. The light of one or several UTs is fed through the FrontEnd unit into optical fibres that scramble the light within the FiberLink unit and provide excellent illumination stability to the spectrograph, using octogonal (1UT) or square (4UT) fibres. The instrument, aiming at a longterm 10 cm s^{−1} RV stability, is temperaturecontrolled and pressurestabilised within a vacuum vessel (VV). The reference fibre fed simultaneously with stabilised FabryPérot unit allows the tracking of instrument drifts down to the cm s^{−1} level. In the most used singleHR (1UT) 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 photondetection 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 COVID19 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 rampup 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^{2} 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 ThAr lamp combined with FabryPé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 crosscorrelation technique; in the middle panel, the RVs computed using the SBART (Silva et al. 2022) code^{3}, a semiBayesian radial velocity computation through template matching (TM); and in the bottom panel, the RVs extracted using the linebyline (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 LombScargle (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 rootmeansquare (RMS) of the RVs are very similar, 1.81, 1.86 and 1.78 m s^{−1} 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^{−1} 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. SBART masks the tellurics at a 1% threshold, which is a quite conservative mask, thus discarding a considerable amount of RV content. SBART 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 SBART 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 timeseries are subsequently with a 1 d step, and after, two E19 points are discarded with RV uncertainties larger than 50 cm s^{−1}, the final sample of ESPRESSO data includes nine E18 and 140 E19 points (see Fig. 2). Full width at half maximum (FWHM) measurements from crosscorrelation functions (CCFs) were automatically provided by the ESPRESSO DRS (see e.g. Fig. 3). The RMS of E18 and E19 RVs computed with SBART are 0.96 and 1.87 m s^{−1}. The RMS of the FWHM of E18 and E19 are 2.46 and 3.34 m s^{−1}. 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 SBART 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, HARPSN 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^{−1} in TMtc RVs compared with an RMS of 0.96 and 1.87 m s^{−1} in TM RVs. The mean/median uncertainties of TMtc RVs are 0.10/0.08 and 0.10/0.09 m s^{−1}, 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, HARPSN and CARMENES datasets described in Sections 2.2 and 2.3.
Stellar properties of GJ 699.
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 SBART template matching (TM) code (middle), and from the linebyline LBL code (bottom). 
Statistics of difference datasets.
Fig. 2 Radialvelocity measurements (right) and GLS periodograms (left) of GJ 699 of ESPRESSO (top), HARPS and HARPSN (middle), and CARMENES (bottom). 
Fig. 3 ESPRESSO FWHM measurements (top), RV measurements (middle), and RV residuals (bottom) from the SHO (P_{ROT} and P_{ROT}/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 HARPSN
The High Accuracy Radial Velocity Planet Searcher (HARPS; Mayor et al. 2003) is a fibrefed échelle highresolution (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^{4} from different ESO programs^{5}, 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 fibrelink 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 ThAr 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 (HARPSN; Cosentino et al. 2012) is a fibrefed échelle highresolution (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^{−1} (Pepe et al. 2014), and covers the wavelength range 383–693 nm. HARPSN 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). HARPSN spectra can be accessed at the TNG archive^{6} from different Spanish CAT programs^{7}. As for HARPS data, the wavelength calibration is done using a ThAr 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 HARPSN spectra is 133 and after binning using 1 d step we finally have 58 HARPSN data points. The exposure time was 900 s before 2020 and the data taking during the COVID19 pandemic was taken with three spectra per night of 1200 s each, except for one night with three spectra of 1800 s. The HARPSN data taken during the COVID19 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 SBART code (see Fig. 2). The CCF FWHM measurements of HARPS and HARPSN spectra were computed by adding a colourcorrection 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^{−1}, 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 SBART are 1.97 and 1.90 m s^{−1}. The RMS of H15 and HAN FWHM measurements are 5.8 and 3.4 m s^{−1} (see Table 2).
2.3 CARMENES
The Calar Alto highResolution search for Mdwarfs with Exoearths with Nearinfrared and optical Échelle Spectrographs (CARMENES; Quirrenbach et al. 2016) are visual (VIS) and nearinfrared (NIR) vacuumstabilised 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 (UAr, UNe, and ThNe) and FabryPé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)^{8}. CARMENES RV measurements were obtained using the templatematching 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^{−1} with a median error bar of 0.9 ms^{−1} (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^{−1} and 26 m s^{−1}, 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^{−1} (see Fig. 2). The RMS of FWHM measurements is 10.9 m s^{−1} (see Table 2).
3 Stellar properties
Barnard’s star (GJ 699) is a bright (m_{V} = 9.5) very nearby M 3.5M 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 massradius 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^{9} 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^{−3} 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)^{10}. 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 longlived 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 crosscorrelation 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 quasiperiodic signals, accounting for changes in the amplitude, phase, or even small period changes. This flexibility can however easily overfit 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 multidimensional 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 activityinduced 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 Mdwarfs.
Following Suárez Mascareño et al. (2023), we use the S+LEAF code (Delisle et al. 2022), which extends the formalism of semiseparable matrices introduced with Celerite (ForemanMackey 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 timelag 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 ForemanMackey et al. (2017), the k_{SHO,i} kernel is defined as: (4)
with η = (1 − (2L/P_{i})^{−2})^{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 multiellipsoidal 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 timeseries with the bestmodel 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); ToledoPadró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 Iindexes, 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 crosscorrelation 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 Mdwarfs (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 logscale to allow for a long tail towards long timescales in the persistence of signals. We use lognormal 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 overfitting 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 (ToledoPadrón et al. 2019), mostly based on the time series of H_{α}index measurements with a 15yr baseline derived from seven different highresolution spectrographs. In all cases the timescale is shorter but consistent with the rotation period. The CCF contrast ( and ), S _{MW} ( and ) and the NaIindex ( 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 1yr 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 Iindex 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 spotinduced 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 timeseries, and not so much in chromospheric indicators, indicating spotdominated 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 (ToledoPadró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 ASASSN Sky Patrol online tool^{11} 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. 24570002 460 300 HJD) and that contain measurements of two different passbands: V (N = 722) and g (N = 2441).
The proper motion of GJ 699 exceeds 10 arcsec/yr, while the ASASSN 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 24570002460300 HJD which approximately reflects the ASASSN coverage of GJ 699. Then, we did split it in timestamps 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 ASASSN Aperture Photometry Pipeline; (iii) clip this light curve to the region defined by a range of ± 50 d within the timestamp. 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^{12} 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 ASASSN 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 longterm cycle and another to model the rotation modulation, where ω_{2} = 2ω_{1} = 2π f / P, with f = 1/t and T_{0,i=1,2} = t_{mid} + P · ϕ_{i}, with t_{mid}, the midtime of the observation baseline and ϕ_{i}, the phase of the sinusoidal function. We left A_{1}, A_{2}, ϕ_{1}, ϕ_{2} 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 ASASSN photometry versus BJD with fitted model. We also show GLS periodograms before and after subtracting the model, with the longterm signal and the rotation signal and its 1yr alias detected. The posterior distributions point to a longterm 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 (ToledoPadró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 longterm cycle. We refer to ToledoPadrón et al. (2019) for a deeper analysis of photometric data of Barnard’s star with a longer baseline of 15 years.
Fig. 5 ASASSN photometry in the same temporal baseline as the RV data together with the double sinusoidal models of the longterm cycle and rotation signals (top left) and the GLS periodogram before and after subtracting the longterm cycle (top right). The posterior and prior distributions of these models including the longterm 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 upperright panel the VV pressure versus BJD keeping reasonably constant at a level of about 0.6 × 10^{−6} 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 upperleft 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 secondorder 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^{−1} with a RMS of 0.84 m s^{−1}, and a mean red detector drift of −0.08 m s^{−1} with a RMS of 0.33 m s^{−1}. The behaviour of both detectors seems better after the intervention in May 2022, with a mean drift and RMS of −0.26 m s^{−1} (0.22 m s^{−1}) and −0.08 m s^{−1} (0.22 m s^{−1}), 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^{−1}, and the peaktopeak RV variation goes from −4 to +4 m s^{−1} (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 overfitting 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 activityonly 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, HARPSN 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 hyperparameters 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 GPonly 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 onedimensional GP with a quasiperiodic (QP) kernel implemented within the George package (ForemanMackey et al. 2014). We also run the onedimensional GP with a quasiperiodic and cosine (QPC) kernel, which integrates the period P in the quasiperiod 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 multidimensional 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 multidimensional 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 longterm 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, HARPSN, and CARMENES due to instrument limitations, may indicate the presence of a longterm 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 HARPSN 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 ASASSN 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 longterm 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 longterm cycle, , the rotation, , and the timescale, . These spectroscopically derived values perfectly match previous determination of the longterm cycle period and the stellar rotation period (ToledoPadró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 longterm cycle, , the rotation, , and the timescale, .
Bayesian evidence of different models.
Fig. 7 Prior and posterior distributions of the longterm cycle, the GP rotation period, and the timescale from the global analysis of FWHM and RV measurements of ESPRESSO, HARPS, and HARPSN (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 activityonly 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 falsealarm 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 activityonly 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 activityonly 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^{−6}, 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 semiamplitude velocity, k_{p}, with a uniform prior (0,5) m s^{−1}. 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 activityonly 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 activityonly model D, although in some cases, as in this particular case, it is another model (model A) the reference activityonly 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 activityonly model (model A in Table 3).
Fig. 9 Evolution of the falsealarm 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^{−1}, 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 semiamplitude of k_{p} = 54 cm s^{−1}. 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 multisinusoidal 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 activityonly model resembles the GPonly model displayed in Fig. 3 (model A in Table 3). The GLS periodogram of the RV residuals after subtracting the multisinusoidal 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 overfitting. 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 activityonly 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 activityonly 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 activityonly model with injected white noise at 30 cm s^{−1} 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 semiamplitude from the global analysis of FWHM and RV measurements of the blind search of the 233d candidate planet using ESPRESSO, HARPS, HARPSN and CARMENES data (left, model O1 in Table 4), a blind search (centreleft, model N2 in Table 4) and a guided search (centreright, 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 superEarth
Ribas et al. (2018) reported the detection of a candidate superEarthlike planet orbiting in an slightly eccentric orbit with a period of 232.8 ± 0.4 d, with a semiamplitude velocity of k_{p} = 1.2 ± 0.1 m s^{−1}. 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 (ToledoPadró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, HARPSN, 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 semiamplitude velocity are depicted in Fig. 11. This model is not satisfactory since the semiamplitude, , 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 activityonly 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 semiamplitude 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 semiamplitude 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^{−1}) 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 semiamplitude velocity of k_{p} = 1.2 m s^{−1}. We run a blind search with priors of P_{orb} (200, 300) d and k_{p} (0, 5) m s^{−1}. The posteriors are P_{orb} = 235 ± 6 d and k_{p} = 1.4 ± 0.4 m s^{−1} (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 activityinduced 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^{−1}, 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^{−1} (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.
Bayesian evidence of models evaluating the 233 d candidate.
Fig. 12 ESPRESSO FWHM measurements (top), RV measurements (middle), RV residuals (next bottom) from longterm cycle (red thick dashed line) and SHO (P_{ROT} and P_{ROT}/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 subEarth mass planet
ESPRESSO RVs show a longterm variation, in particular, in the E19 dataset, which may be interpreted as a longterm 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) (lognormal 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 HARPSN (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 activityonly 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 HARPSN 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 activityonly 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 semiamplitude 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 longterm 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^{−1} after subtracting the activity model, ending at 0.45 ms^{−1} 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^{−1}, compared with the median SBART RV uncertainties of 0.10 ms^{−1}. It is remarkable that the RMS of the RV residuals is 0.45 m s^{−1}, greater than the RMS of the FWHM residuals of 0.35 m s^{−1}, thus requiring an additional 0.56 m s^{−1} RV jitter while a 0.30 m s^{−1} 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^{−1} (~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 subEarth 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 semiamplitude 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 activityonly 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. Activityonly 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 activityonly model DT, corresponding to 0.25 and 0.05% FAP. In these cases, the semiamplitude 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 HARPSN, 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 SBART code (Silva et al. 2022), to explore any possible chromatic effects on the 3.15 d signal. We compute RVs using either only the bluedetector or the reddetector 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 reddetector 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 bluedetector 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 semiamplitude velocity these models give k_{p} = [0.40, 0.54] m s^{−1}. We consider these results reasonable given the low semiamplitude 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 semiamplitude. 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 HARPSN data together reinforces the confirmation of the 3.15 d signal as a subEarth mass planet. Figure D.3 shows the result of a model consisting of longterm 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 activityonly 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 HARPSN spectrographs. Figure D.1 displays the RV data versus orbital phase, where the RMS of the original ESPRESSO, HARPS and HARPSN RVs of 1.93 m s^{−1} goes down to 0.97 m s^{−1}, after subtracting the activityonly model, and down to 0.91 ms^{−1} 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 ESPRESSOonly 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^{−1} (6σ detection), with m_{p} sin i = 0.33 ± 0.06 M_{⊕}, consistent with ESPRESSOonly 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 activityonly 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, HARPSN and CARMENES data goes from the original values at 2.08 m s^{−1} down to 1.49 m s^{−1}, after subtracting the activityonly model, and down to 1.47 m s^{−1}, 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.
TM RV and CCF FHWM measurements of ESPRESSO datasets.
Fig. 13 Corner plot with the posterior distributions of the orbital parameters of the subEarthmass planet of GJ 699 with a 3.15 d orbital period, with prior and posterior parameters given in Table 6. 
Prior and posteriors of the activity+planet ESPRESSO model.
Parameters of planet Barnard b.
Fig. 14 RV curve of the subEarthmass 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 multiplanetary 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^{−1} 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 semiamplitude 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 twoplanet model runs by adding the datasets of HARPS and HARPSN (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 singleplanet model in each case, but not significant. It is remarkable that for the M2 model the periods and the semiamplitudes, 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 ESPRESSOonly 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 ESPRESSOonly 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 loguniform priors for periods P_{orb} (2, 40) d days, and uniform prior on semiamplitudes k_{p} (0, 5) m s^{−1}. Fig. 15 shows the results of the − log_{10} 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 injectionrecovery test. We construct activityonly RVs by generating random samples from the GP model, and adding white noise with the same standard deviation of the residuals of the 4planet model (see Appendix E). We then inject 100000 random sinusoidal signals with periods between 1 and 1100 d, semiamplitudes between 1 cm s^{−1} and 11 m s^{−1}, 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 hyperparameters as measured for our fourplanet 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^{−1} for periods from 1.2 d up to ~20 d. Our sensitivity drops to ~1 m s^{−1} at 35 d and then quickly plummets to 2.5 m s^{−1} 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 GPonly models, which try to account for all variations at low frequencies up to the characteristic value specified by the kernel. To extract possible lowamplitude 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 multiplanetary system.
The subEarth 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 subEarth 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 8m class telescope. The endeavor of filling the region in Fig. 17 below the 1 m s^{−1} precision limit at periods larger than 10 d towards finding longperiod Earthlike 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 Earthlike planets in the future with groundbased 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 HipparcosGaia 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 HipparcosGaia 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 propermotion difference technique, where the catalogue propermotion values for Hipparcos and Gaia are averaged over the timespan of the observations, thus with smearing orbital effects for shorter periods. The last peak in sensitivity loss appears, in fact, at the exact timespan 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 timespan from DR2 to DR3 over which the Gaia propermotion 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 superEarths with true masses around 10 M_{⊕} can be ruled out by the lack of statistically significant PMA. The combination of ESPRESSO RVs and HipparcosGaia absolute astrometry thus allows to infer that superEarths or larger planets are not present in the GJ 699 system out of ~10 AU.
Fig. 15 FIP periodogram of the 4planet model of the ESPRESSO data of GJ 699. The periods of the detected signals are indicated in red solid circles. The − log_{10} 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: semiamplitude k_{p} versus period P_{orb} 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 semiamplitudes 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 solartype stars, together with those discovered and confirmed using ESPRESSO (green circles). Confirmed planet Barnard b (red), and planet candidates (purple) from the 4planet 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 semiamplitude of 1 cm s^{−1} (green line) and 10 cm s^{−1} (blue line), respectively, and null eccentricity. 
Fig. 18 HipparcosGaia PMA sensitivity to companions of a given mass (in M_{⊕}) as a function of the orbital semimajor axis (in AU) orbiting Barnard’s star. The solid black line identifies the combinations of mass and separation explaining the observed propermotion anomaly (PMA) at the mean epoch of Gaia DR3 (Eq. (15) in Kervella et al. 2019). The shaded lightblue 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^{−1} RMS with about 8 m s^{−1} peaktopeak RV measurements, which mainly reflect the magnetic activity of this old middleMtype dwarf star. We modelled these RV measurements together with CCF FWHM measurements, with a 15 m s^{−1} peak to peak and 3.3 m s^{−1} RMS, using a longterm cycle described with a double sinusoidal model, and a rotationinduced activity model, which is described using a Gaussian process approach. The activity cycle is better constrained using additional HARPS, HARPSN, and CARMENES data, which help to extend the baseline of the observations up to eight years. We obtained a welldescribed activity model with a cycle period of and a rotation period of , which are consistent with previous results in the literature (e.g. ToledoPadrón et al. 2019).
The high quality of the ESPRESSO data means that we can evaluate the presence of the candidate superEarthlike 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^{−1}, with all the models giving solutions with semiamplitudes 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 semiamplitude of 55 ± 7 cm s^{−1}, thus uncovering a subEarthmass planet of 0.37 ± 0.05 M_{⊕}, which is about half of the mass of Venus or three times that of Mars. This subEarthmass 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 fourplanet system with semiamplitudes of k_{p} = 47, 41, 35, and 20 cm s^{−1}, which would correspond to a system of four subEarthmass 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 semimajor 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 fourplanet 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 Earthand subEarthmass 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/vizbin/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 PID2020117493GBI00. 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. Cofunded 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 PID2022137241NBC41 y PID2022137241NBC44. 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 HARPSN 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/Cofunded 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 PID2019107061GBC61. 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 SBART (Silva et al. 2022). Second RV computation with LBL (Artigau et al. 2022). Extensive use of the DACE platform^{13} 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 (ForemanMackey et al. 2017) and George (ForemanMackey 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 signaltonoise 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^{−1}). 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^{−1}), the lower semiamplitude 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 longterm 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 multisinusoidal function that uses up to six sinusoids with periods in the range [59,255] d, and RV residuals after subtracting this activityonly 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 activityonly 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 activityonly model with injected white noise and injecting a planet signal at 3.15 d, and RV residuals after subtracting this activityonly 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). 
Bayesian evidence of different models
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 crosscorrelation 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, pvalue (the lower the value, the better the correlation), and the slope we obtain when doing a leastsquares fit.
Spearman’s correlation coefficient, and measured slopes, between activity indicators, their gradient, and the RV measurements.
Fig. C.1 Analysis of the ESPRESSO S_{MW}, Hα and Na I spectroscopic indexes, a, b: S_{MW}, Hα and Na Iindex timeseries with the bestmodel 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 S_{MW}, Hα and Na Iindex data. The red vertical dashed line shows the most significant period, d: Relationship between the S_{MW}, Hα and Na Iindex data. The best fit is shown when the slope is ≥3σ different from zero, e, f: Residuals of the S_{MW}, Hα and Na Iindex after subtracting the best model fit. g: GLS periodogram of the residuals, h: Comparison of the CCF RV and gradient of the S_{MW}, Hα and Na Iindex model. 
Fig. C.2 Analysis of the BIS and contrast of the ESPRESSO CCF. a, b: BIS and contrast timeseries with the bestmodel 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 HARPSN 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 subEarth 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 subEarth 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, HARPSN 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 subEarthmass planet of GJ 699 with a 3.15 d orbital period together with ESPRESSO, HARPS and HARPSN 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 subEarthmass planet of GJ 699 with a 3.15 d from the model using ESPRESSO, HARPS and HARPSN 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 HARPSN datasets using SHO (P_{ROT} and P_{ROT}/2) GP model (grey solid line) to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the longterm 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 HARPSN data together, including offsets, jitter, the parameters of GP and the longterm cycle, and the parameters of the Keplerian model revealing the subEarthmass 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, HARPSN and CARMENES datasets using SHO (P_{ROT} and P_{ROT}/2) GP model (grey solid line) to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the longterm 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. 
TM RV and CCF FHWM measurements of all datasets
Appendix E Candidate multiplanet 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^{−1}. The next peak in the RV residuals is 4.12 d, so we then run a twoplanet 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 singleplanet model F1. In the RV residuals appears the peak at 2.34 d significantly above the 0.1% FAP line.
Bayesian evidence of different models
We run a threeplanet 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 twoplanet model F2. We also check the possibility of a wider prior on eccentricity (e with and (0, 0.3)) of the threeplanet 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 shortperiod 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 fourplanet 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 fourplanet 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^{−1} in the RVs after subtracting the activity cycle+rotation model to 0.49, 0.39, 0.29, and 0.24 m s^{−1} 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 semiamplitude velocities of run F4 are k_{p} = 47, 41, 35 and 20 cm s^{−1}, 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^{−1} on semiamplitude 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 semiamplitude 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 1yr alias (6.61 d and 6.86 d) of the main period 6.74 d. The semiamplitudes 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 multiplanet system
We evaluate the stability of the 4planet candidate system using the SPOCK^{14} tool (Tamayo et al. 2020), together with the Nbody Regressor (Hussain & Tamayo 2020) and the REBOUND^{15} code (e.g. Rein & Tamayo 2017). We found that assuming the zero eccentricity^{16} for all planets, the system remains stable more than 10^{9} 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 fourplanet system looks particularly compact, similar to other systems such as e.g. the TRAPPIST sevenplanet 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 Earthmass 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 semimajor 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 lowmass, 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 multiplanet systems suggest a loglinear 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^{9} 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 sevenplanet system TRAPPIST1 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 threeplanet 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 (AstudilloDefru 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 threeplanet 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 fiveplanet system at Δ(R_{H}) < 13 (Dreizler et al. 2024). On the other hand, the twoplanet 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 subEarth 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 fourplanet system would require further investigation.
Fig. E.1 ESPRESSO RV measurements versus BJD after removing the longterm cycle and the GP model (Left), GLS periodograms (middle), and RV curves of the subEarthmass 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, HARPSN and CARMENES datasets using SHO (P_{ROT} and P_{ROT}/2) GP model to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the longterm 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 semiamplitude 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 (centerleft), planet d (centerright), and planet e (right). For all planets we use an uniform prior (0, 5) m s^{−1} on semiamplitude velocity. 
Parameters of the candidate planetary system
The candidate planetary system orbiting our second closest neighbour Barnard’s star consists of four subEarth 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 fourplanet 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
 Agol, E., Dorn, C., Grimm, S. L., et al. 2021, Planet. Sci. J., 2, 1 [Google Scholar]
 Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147 [Google Scholar]
 Allart, R., Lovis, C., Faria, J., et al. 2022, A&A, 666, A196 [Google Scholar]
 AlonsoFloriano, F. J., Morales, J. C., Caballero, J. A., et al. 2015, A&A, 577, A128 [Google Scholar]
 AngladaEscudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437 [Google Scholar]
 Artigau, É., Cadieux, C., Cook, N. J., et al. 2022, AJ, 164, 84 [Google Scholar]
 AstudilloDefru, N., Delfosse, X., Bonfils, X., et al. 2017a, A&A, 600, A13 [Google Scholar]
 AstudilloDefru, N., Forveille, T., Bonfils, X., et al. 2017b, A&A, 602, A88 [Google Scholar]
 Azevedo Silva, T., Demangeon, O. D. S., Santos, N. C., et al. 2022, A&A, 666, L10 [Google Scholar]
 Barnard, E. E. 1916, AJ, 29, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Barragán, O., Aigrain, S., Rajpaul, V. M., & Zicher, N. 2022, MNRAS, 509, 866 [Google Scholar]
 Barros, S. C. C., Demangeon, O. D. S., Alibert, Y., et al. 2022, A&A, 665, A154 [Google Scholar]
 Benedict, G. F., McArthur, B., Nelan, E., et al. 1998, AJ, 116, 429 [Google Scholar]
 Benz, W., Broeg, C., Fortier, A., et al. 2021, Exp. Astron., 51, 109 [Google Scholar]
 Borsa, F., Allart, R., CasasayasBarris, N., et al. 2021, A&A, 645, A24 [EDP Sciences] [Google Scholar]
 Borucki, W., Koch, D., Batalha, N., et al. 2009, IAU Symp., 253, 289 [Google Scholar]
 Bouchy, F., Doyon, R., Artigau, É., et al. 2017, The Messenger, 169, 21 [NASA ADS] [Google Scholar]
 Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012, ApJ, 757, 112 [Google Scholar]
 Brandt, T. D. 2021, ApJS, 254, 42 [Google Scholar]
 Cadieux, C., Plotnykov, M., Doyon, R., et al. 2024, ApJ, 960, L3 [NASA ADS] [CrossRef] [Google Scholar]
 CastroGonzález, A., Demangeon, O. D. S., LilloBox, J., et al. 2023, A&A, 675, A52 [Google Scholar]
 Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996, Icarus, 119, 261 [Google Scholar]
 Cosentino, R., Lovis, C., Pepe, F., et al. 2012, Proc. SPIE, 8446, 84461V [Google Scholar]
 Cristofari, P. I., Donati, J. F., Moutou, C., et al. 2023, MNRAS, 526, 5648 [Google Scholar]
 Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
 Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, Astrophysics Source Code Library [record ascl:1906.010] [Google Scholar]
 Delisle, J. B., Unger, N., Hara, N. C., & Ségransan, D. 2022, A&A, 659, A182 [Google Scholar]
 Demangeon, O. D. S., Zapatero Osorio, M. R., Alibert, Y., et al. 2021, A&A, 653, A41 [Google Scholar]
 Di Marcantonio, P., Sosnowska, D., Cupani, G., et al. 2018, SPIE Conf. Ser., 10704, 107040F [Google Scholar]
 Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Nature, 544, 333 [Google Scholar]
 Donati, J. F., Kouach, D., Moutou, C., et al. 2020, MNRAS, 498, 5684 [Google Scholar]
 Donati, J. F., Lehmann, L. T., Cristofari, P. I., et al. 2023, MNRAS, 525, 2015 [CrossRef] [Google Scholar]
 Dreizler, S., Luque, R., Ribas, I., et al. 2024, A&A, 684, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132 [Google Scholar]
 Eastman, J., Siverd, R., & Gaudi, B. S. 2010, PASP, 122, 935 [Google Scholar]
 Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83 [Google Scholar]
 Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, Nature, 580, 597 [Google Scholar]
 Faria, J. P., Haywood, R. D., Brewer, B. J., et al. 2016, A&A, 588, A31 [Google Scholar]
 Faria, J. P., Suárez Mascareño, A., Figueira, P., et al. 2022, A&A, 658, A115 [Google Scholar]
 Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601 [Google Scholar]
 ForemanMackey, D., Hoyer, S., Bernhard, J., & Angus, R. 2014, https://doi.org/10.5281/zenodo.11989 [Google Scholar]
 ForemanMackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220 [Google Scholar]
 France, K., Duvvuri, G., Egan, H., et al. 2020, AJ, 160, 237 [Google Scholar]
 Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504 [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [Google Scholar]
 Giles, H. A. C., Collier Cameron, A., & Haywood, R. D. 2017, MNRAS, 472, 1618 [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Demory, B.O., et al. 2017, Nature, 542, 456 [Google Scholar]
 Gladman, B. 1993, Icarus, 106, 247 [Google Scholar]
 González Hernández, J. I., Pepe, F., Molaro, P., & Santos, N. C. 2018, in Handbook of Exoplanets, eds. H. J. Deeg, & J. A. Belmonte (Berlin: Springer), 157 [Google Scholar]
 Gratia, P., & Lissauer, J. J. 2021, Icarus, 358, 114038 [NASA ADS] [CrossRef] [Google Scholar]
 Gullikson, K., DodsonRobinson, S., & Kraus, A. 2014, AJ, 148, 53 [Google Scholar]
 Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015a, MNRAS, 450, L61 [Google Scholar]
 Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015b, MNRAS, 453, 4384 [Google Scholar]
 Hara, N. C., Boué, G., Laskar, J., Delisle, J. B., & Unger, N. 2019, MNRAS, 489, 738 [Google Scholar]
 Hara, N. C., Unger, N., Delisle, J.B., Díaz, R. F., & Ségransan, D. 2022, A&A, 663, A14 [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Hussain, N., & Tamayo, D. 2020, MNRAS, 491, 5258 [Google Scholar]
 Jahandar, F., Doyon, R., Artigau, É., et al. 2023, arXiv eprints [arXiv:2310.12125] [Google Scholar]
 Kervella, P., Arenou, F., Mignard, F., & Thévenin, F. 2019, A&A, 623, A72 [Google Scholar]
 Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7 [Google Scholar]
 Kirkpatrick, J. D., Gelino, C. R., Cushing, M. C., et al. 2012, ApJ, 753, 156 [Google Scholar]
 Koen, C., Kilkenny, D., van Wyk, F., & Marang, F. 2010, MNRAS, 403, 1949 [Google Scholar]
 Kopparapu, R. K., Ramirez, R. M., SchottelKotte, J., et al. 2014, ApJ, 787, L29 [Google Scholar]
 Kürster, M., Endl, M., Rouesnel, F., et al. 2003, A&A, 403, 1077 [Google Scholar]
 Lafarga, M., Ribas, I., Lovis, C., et al. 2020, A&A, 636, A36 [Google Scholar]
 Laskar, J., & Petit, A. C. 2017, A&A, 605, A72 [Google Scholar]
 Lavie, B., Bouchy, F., Lovis, C., et al. 2023, A&A, 673, A69 [Google Scholar]
 LilloBox, J., Figueira, P., Leleu, A., et al. 2020, A&A, 642, A121 [Google Scholar]
 LilloBox, J., Faria, J. P., Suárez Mascareño, A., et al. 2021, A&A, 654, A60 [Google Scholar]
 Lovis, C., Ségransan, D., Mayor, M., et al. 2011, A&A, 528, A112 [Google Scholar]
 Lubin, J., Robertson, P., Stefansson, G., et al. 2021, AJ, 162, 61 [Google Scholar]
 Marconi, A., Abreu, M., Adibekyan, V., et al. 2022, SPIE Conf. Ser., 12184, 1218424 [NASA ADS] [Google Scholar]
 Marfil, E., Tabernero, H. M., Montes, D., et al. 2021, A&A, 656, A162 [Google Scholar]
 Martins, C. J. A. P., Cristiani, S., Cupani, G., et al. 2022, Phys. Rev. D, 105, 123507 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 Murphy, M. T., Molaro, P., Leite, A. C. O., et al. 2022, A&A, 658, A123 [Google Scholar]
 Palle, E., Biazzo, K., Bolmont, E., et al. 2023, arXiv eprints [arXiv:2311.17075] [Google Scholar]
 Passegger, V. M., Reiners, A., Jeffers, S. V., et al. 2018, A&A, 615, A6 [Google Scholar]
 Pepe, F., Ehrenreich, D., & Meyer, M. R. 2014, Nature, 513, 358 [Google Scholar]
 Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96 [Google Scholar]
 Perger, M., AngladaEscudé, G., Ribas, I., et al. 2021, A&A, 645, A58 [EDP Sciences] [Google Scholar]
 Quanz, S. P., Ottiger, M., Fontanet, E., et al. 2022, A&A, 664, A21 [Google Scholar]
 Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279 [Google Scholar]
 Quirrenbach, A., Amado, P. J., Caballero, J. A., et al. 2016, Proc. SPIE, 9908, 990812 [Google Scholar]
 Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, MNRAS, 452, 2269 [Google Scholar]
 Rajpaul, V., Aigrain, S., & Roberts, S. 2016, MNRAS, 456, L6 [Google Scholar]
 Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning (Cambridge, Massachusetts: The MIT Press) [Google Scholar]
 Rein, H., & Tamayo, D. 2017, MNRAS, 467, 2377 [NASA ADS] [Google Scholar]
 Reiners, A., Shulyak, D., Käpylä, P. J., et al. 2022, A&A, 662, A41 [Google Scholar]
 Reylé, C., Jardine, K., Fouqué, P., et al. 2021, A&A, 650, A201 [Google Scholar]
 Ribas, I., Tuomi, M., Reiners, A., et al. 2018, Nature, 563, 365 [Google Scholar]
 Ribas, I., Reiners, A., Zechmeister, M., et al. 2023, A&A, 670, A139 [Google Scholar]
 Rice, D. R., & Steffen, J. H. 2023, MNRAS, 520, 4057 [NASA ADS] [CrossRef] [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telesc. Instrum. Syst., 1, 014003 [Google Scholar]
 Robertson, P., Mahadevan, S., Endl, M., & Roy, A. 2014, Science, 345, 440 [Google Scholar]
 Schweitzer, A., Passegger, V. M., Cifuentes, C., et al. 2019, A&A, 625, A68 [Google Scholar]
 Silva, A. M., Faria, J. P., Santos, N. C., et al. 2022, A&A, 663, A143 [Google Scholar]
 Skilling, J. 2004, AIP Conf. Ser., 735, 395 [Google Scholar]
 Sobell, M. G. 2015, A Practical Guide to Ubuntu Linux (Noida: Pearson Education) [Google Scholar]
 Speagle, J. S. 2020, MNRAS, 493, 3132 [Google Scholar]
 Stelzer, B., Marino, A., Micela, G., LópezSantiago, J., & Liefke, C. 2013, MNRAS, 431, 2063 [Google Scholar]
 Stock, S., Kemmer, J., Reffert, S., et al. 2020, A&A, 636, A119 [Google Scholar]
 Suárez Mascareño, A., Rebolo, R., González Hernández, J. I., & Esposito, M. 2015, MNRAS, 452, 2745 [Google Scholar]
 Suárez Mascareño, A., Rebolo, R., González Hernández, J. I., & Esposito, M. 2017, MNRAS, 468, 4772 [Google Scholar]
 Suárez Mascareño, A., Rebolo, R., González Hernández, J. I., et al. 2018, A&A, 612, A89 [Google Scholar]
 Suárez Mascareño, A., Faria, J. P., Figueira, P., et al. 2020, A&A, 639, A77 [Google Scholar]
 Suárez Mascareño, A., GonzálezÁlvarez, E., Zapatero Osorio, M. R., et al. 2023, A&A, 670, A5 [Google Scholar]
 Suárez Mascareño, A., Passegger, V. M., González Hernández, J. I., et al. 2024, A&A, 685, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tabernero, H. M., Marfil, E., Montes, D., & González Hernández, J. I. 2022, A&A, 657, A66 [Google Scholar]
 Tamayo, D., Cranmer, M., Hadden, S., et al. 2020, Proc. Natl. Acad. Sci., 117, 18194 [Google Scholar]
 Thompson, S. J., Queloz, D., Baraffe, I., et al. 2016, SPIE Conf. Ser., 9908, 99086F [Google Scholar]
 ToledoPadrón, B., González Hernández, J. I., RodríguezLópez, C., et al. 2019, MNRAS, 488, 5145 [Google Scholar]
 Van Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: CreateSpace) [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48 [Google Scholar]
 Wildi, F., Pepe, F., Chazelas, B., Lo Curto, G., & Lovis, C. 2010, SPIE Conf. Ser., 7735, 77354X [NASA ADS] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [Google Scholar]
 Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12 [Google Scholar]
 Zechmeister, M., Dreizler, S., Ribas, I., et al. 2019, A&A, 627, A49 [NASA ADS] [EDP Sciences] [Google Scholar]
 Zsom, A., Seager, S., de Wit, J., & Stamenković, V. 2013, ApJ, 778, 109 [Google Scholar]
All Tables
Spearman’s correlation coefficient, and measured slopes, between activity indicators, their gradient, and the RV measurements.
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 SBART template matching (TM) code (middle), and from the linebyline LBL code (bottom). 

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

In the text 
Fig. 3 ESPRESSO FWHM measurements (top), RV measurements (middle), and RV residuals (bottom) from the SHO (P_{ROT} and P_{ROT}/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 timeseries with the bestmodel 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 ASASSN photometry in the same temporal baseline as the RV data together with the double sinusoidal models of the longterm cycle and rotation signals (top left) and the GLS periodogram before and after subtracting the longterm cycle (top right). The posterior and prior distributions of these models including the longterm 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 longterm cycle, the GP rotation period, and the timescale from the global analysis of FWHM and RV measurements of ESPRESSO, HARPS, and HARPSN (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 falsealarm 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 semiamplitude from the global analysis of FWHM and RV measurements of the blind search of the 233d candidate planet using ESPRESSO, HARPS, HARPSN and CARMENES data (left, model O1 in Table 4), a blind search (centreleft, model N2 in Table 4) and a guided search (centreright, 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 longterm cycle (red thick dashed line) and SHO (P_{ROT} and P_{ROT}/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 subEarthmass 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 subEarthmass 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 4planet model of the ESPRESSO data of GJ 699. The periods of the detected signals are indicated in red solid circles. The − log_{10} 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: semiamplitude k_{p} versus period P_{orb} 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 semiamplitudes 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 solartype stars, together with those discovered and confirmed using ESPRESSO (green circles). Confirmed planet Barnard b (red), and planet candidates (purple) from the 4planet 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 semiamplitude of 1 cm s^{−1} (green line) and 10 cm s^{−1} (blue line), respectively, and null eccentricity. 

In the text 
Fig. 18 HipparcosGaia PMA sensitivity to companions of a given mass (in M_{⊕}) as a function of the orbital semimajor axis (in AU) orbiting Barnard’s star. The solid black line identifies the combinations of mass and separation explaining the observed propermotion anomaly (PMA) at the mean epoch of Gaia DR3 (Eq. (15) in Kervella et al. 2019). The shaded lightblue 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 multisinusoidal function that uses up to six sinusoids with periods in the range [59,255] d, and RV residuals after subtracting this activityonly 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 activityonly 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 activityonly model with injected white noise and injecting a planet signal at 3.15 d, and RV residuals after subtracting this activityonly 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 S_{MW}, Hα and Na I spectroscopic indexes, a, b: S_{MW}, Hα and Na Iindex timeseries with the bestmodel 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 S_{MW}, Hα and Na Iindex data. The red vertical dashed line shows the most significant period, d: Relationship between the S_{MW}, Hα and Na Iindex data. The best fit is shown when the slope is ≥3σ different from zero, e, f: Residuals of the S_{MW}, Hα and Na Iindex after subtracting the best model fit. g: GLS periodogram of the residuals, h: Comparison of the CCF RV and gradient of the S_{MW}, Hα and Na Iindex model. 

In the text 
Fig. C.2 Analysis of the BIS and contrast of the ESPRESSO CCF. a, b: BIS and contrast timeseries with the bestmodel 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 subEarthmass planet of GJ 699 with a 3.15 d orbital period together with ESPRESSO, HARPS and HARPSN 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 subEarthmass planet of GJ 699 with a 3.15 d from the model using ESPRESSO, HARPS and HARPSN 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 HARPSN datasets using SHO (P_{ROT} and P_{ROT}/2) GP model (grey solid line) to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the longterm 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 HARPSN data together, including offsets, jitter, the parameters of GP and the longterm cycle, and the parameters of the Keplerian model revealing the subEarthmass 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, HARPSN and CARMENES datasets using SHO (P_{ROT} and P_{ROT}/2) GP model (grey solid line) to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the longterm 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 longterm cycle and the GP model (Left), GLS periodograms (middle), and RV curves of the subEarthmass 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, HARPSN and CARMENES datasets using SHO (P_{ROT} and P_{ROT}/2) GP model to describe the activity caused by stellar rotation and including a double sinusoidal model to describe the longterm 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 semiamplitude 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 (centerleft), planet d (centerright), and planet e (right). For all planets we use an uniform prior (0, 5) m s^{−1} on semiamplitude velocity. 

In the text 
Fig. E.4 Same as Fig. 17 with the confirmed planet Barnard b (red), and planet candidates (purple) of the fourplanet 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.