© The Authors 2025

1 Introduction

The search for exoplanets around nearby stars is now mainstream science, with roughly 6000 planets known to date¹. Within this sphere, the interest for the lowest-mass exoplanets was recently spurred by the capabilities of the ESPRESSO spectrograph (Pepe et al. 2021). Designed with the objective of detecting an Earth-mass planet within the habitable zone (HZ) of a Sun-like star, the instrument was built to reach a radial velocity (RV) precision of down to 10 cm/s over timescales of up to several years. The first studies in the field, such as that of Faria et al. (2020) and Lillo-Box et al. (2020), led to very detailed planetary system characterisations. The data obtained in the context of the consortium’s Guaranteed Time of Observations (GTO) enabled results such as those reported in Suárez Mascareño et al. (2020) and Faria et al. (2022), whereby a system of two planets was detected around Proxima Cen and the M star’s rotational-modulated activity was modelled and corrected down to a level of 30 cm/s, close to the photon noise of the measurements.

As discussed in Pepe et al. (2021), the GTO team, the associated observations, and their analysis are organised into working groups (WG). WG1 is dedicated to the detection and characterisation of Earth-like planets via blind RV searches. Proxima was one of ~25 stars observed intensively within ESPRESSO GTO WG1 to cover periodicities up to the stars’ HZs. For the preparatory work and original list of stars, we refer to Hojjatpanah et al. (2019).

In this work, we selected four ESPRESSO GTO WG1 stars for a detailed analysis: HD 10700 (τ Ceti), HD 20794 (e Eridani), HD 102365, and HD 304636. These specific stars were selected for having enough data points for an RV characterisation, while also being representative of our blind RV campaign. We first used ESPRESSO and ancillary data to characterise each star and then followed up with RV time series calculations and their analysis. We study the presence and characteristics of signals from the shortest to the longest timescales: pulsations, granulation, stellar activity, and, finally, the planetary signals up to the HZ of the host stars. We discuss the currently achievable level of radial velocity precision, our analysis choices, and the current limitations, along with potential avenues to explore.

In Section 2, we present a summary of the literature on the targets. In Section 3, we perform our own stellar characterisation using a range of methods. In Section 4, we describe the ESPRESSO data acquisition and RV derivation. We discuss the intra-night RV variations in Section 5, along with the night-to-night and longer-term RVs and their associated quantities variations in Section 6. In Section 7, we characterise the longterm variations with a planetary and stellar activity model. We discuss our results in Section 8 and wrap up our conclusions in Section 9.

2 Targets and previous works

2.1 HD 10700 (τ Ceti)

The star HD 10700, commonly know as τ Ceti, is a very bright (V = 3.50) star located only at 3.65 pc from the Sun. It is stable photometrically and has been considered to be representative of chromosphericaly inactive stars, as underlined by Gray & Baliunas (1994). The same work reports no evidence for a rotational modulation on temperature, granulation or chromospheric activity, along with only a weak magnetic cycle at a period of roughly 11 yr.

The star was target of an RV asteroseismology campaign with HARPS, providing accurate values for stellar mass, radius, and luminosity (Teixeira et al. 2009). Recently, Korolik et al. (2023) used long-baseline optical interferometric to provide additional constraints on stellar parameters, deriving a rotation period, P_rot, of 46±4 d and a stellar axis inclination i of 7±7^°. A brief summary of the parameters available in the literature is given in Table 1.

The star’s brightness and perceived stability made it a prime target for RV studies. The work of Pepe et al. (2011), using HARPS (Mayor et al. 2003) data accumulated over more than 7 years, reported an RV dispersion of only 1.5 m/s for intra-night measurements and a long-term scatter of 92 cm/s when nightly averaged RVs are considered. This study reported that several signals seemed to be present in the data but none was confidently detected.

The work of Tuomi et al. (2013) put together HARPS together with AAPS and HIRES data to find evidence for five planets at periods of 13.9, 35.4, 94, 168, and 640 d with a minimum mass of 2.0, 3.1, 3.6, 4.3, and 6.6 M_⊕, respectively, and amplitudes between 58 and 75 cm/s. The same team published an updated analysis in Feng et al. (2017b) reporting the detection of only four planets, with periods of 20.0, 49.4, 162, and 636 d and a minimum mass of 1.75, 1.83, 3.93, and 3.93 M_⊕, having amplitudes between 35 and 55 cm/s².

2.2 HD 20794 (e Eridani)

The bright (V = 4.26) star HD 20794, also known as e Eridani, has been intensively studied in exoplanet searches. It was subject to an intense RV monitoring with HARPS and had three super-Earths announced in Pepe et al. (2011). The data were later revisited by Feng et al. (2017a), who confirmed two of the previously presented planets and reported a different characterisation of the third one. The stellar characterisation available is summarised in Table 1. A recent analysis of ESPRESSO data on HD 20794 is presented in Nari et al. (2025). One particular study that is especially pertinent to our analysis is the commissioning work on HARPS, presented in Mayor et al. (2003), which reported a root-mean-square (rms) RV dispersion of HD 20794 of 80 cm/s, ascribed to p-modes.

2.3 HD 102365

HD 102365 is a bright (V = 4.89) star that stood out among ESPRESSO GTO WG1 targets for its brightness, low rotational velocity, and low activity level (as measured by its $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ index). The main stellar properties are compiled in Table 1.

Tinney et al. (2011) detected a planet around the star using 149 UCLES (Diego et al. 1990) RVs. A period inspection was done with Lomb-Scargle (LS) tools and a significant peak with a false alarm probability (FAP) of 0.1% was found at 122 d. The best-fit parameters led to a planetary orbital period of 122.1±0.3 d, an amplitude of 2.40±0.35 m/s, and an eccentricity of 0.34±0.14; these values corresponded to a lower mass limit m. sin i of 16.0±2.6 M_⊕.

Laliotis et al. (2023) performed an archival RV analysis on several stars, covering HD 102365. The authors updated the orbital parameters using additional observations; the previously published planetary period and eccentricity are recovered, but the semi-amplitude is revised down to 1.38±0.23 m/s, a value 58% lower than announced by Tinney et al. (2011). Finally, Hill et al. (2023) consider HD 102365 b as the most favourable HZ non-transiting exoplanets for a spectroscopic characterisation, underlining the proposed planet’s importance.

2.4 HD 304636

HD 304636 is listed as an M0 from Hipparcos colours (Koen et al. 2010) that is comparatively bright for its assigned spectral type (V = 9.465). Its parameters were derived in Maldonado et al. (2019b) using ratios of pseudo-equivalent widths of spectral features sensitive to the effective temperature and the stellar metallicity. The literature values on stellar characterisation are presented in Table 1.

The TESS transit survey (Ricker et al. 2014) identified a planet candidate around HD 304636 with a periodicity of 7.57 d and a transit depth of 0.29 parts per thousand (ppt), corresponding to a radius of ~0.87 R_⊕. The star was then listed as a TESS Object of Interest, TOI-741. In a recent work DiTomasso et al. (2025) used HARPS RV observations to attempt a detection of the candidate. While a detection was not possible, the authors put an upper limits of 4.0 M_⊕ (at 3σ) on the mass of the planet (if indeed present).

3 Stellar host characterisation

The modelling of stellar phenomena capable of creating RV signals requires an accurate determination of the stellar parameters. We use different methodologies and ancillary data to do so. After critically assessing and comparing the derived values, we conclude by summarising the best estimation for each star’s parameters.

3.1 Spectroscopic analysis

We used ESPRESSO spectra to derive the stellar atmospheric parameters T_eff, log g, microturbulence ξ_t, and [Fe/H]. The individual spectra were co-added to create a master spectrum with S/N up to 2500 for each star. In this regime, the uncertainty on the derived parameters is not limited by the spectra S/N, but by the intrinsic precision and accuracy of the method instead.

For the first three stars, we used the methodology described in Sousa et al. (2021) and references therein. In brief, we used ARES³ (Sousa et al. 2007, 2015) to measure the equivalent widths (EW) of iron lines on the combined spectra followed by a minimisation process to find ionisation and excitation equilibrium and converge to the best set of spectroscopic parameters. This process makes use of a grid of Kurucz model atmospheres (Kurucz 1993) and the radiative transfer code MOOG (Sneden 1973).

For M-stars like HD 304636 the stellar parameters cannot be derived with the same technique. We used the ODUSSEAS code (Antoniadis-Karnavas et al. 2020)⁴ that employs a machine learning approach on pseudo-EW of more than 4000 stellar absorption lines to predict the stellar effective temperature and metallicity. We applied the calibration presented in Antoniadis-Karnavas et al. (2024), which was demonstrated to provide more accurate and consistent results than its previous iteration. In parallel, we used the Bayesian tool SteParSyn described in Tabernero et al. (2022).

We used the derived T_eff and GAIA DR3 (Gaia Collaboration 2023) data and followed the procedure described in Sousa et al. (2021) to calculate the trigonometric surface gravity log g_trig. The work of Pecaut & Mamajek (2013) was used to assign a spectral type to the T_eff of each star. The spectroscopic parameters T_eff, log g, microturbulence ξ_t, and [Fe/H], plus the associated log g_trig and spectral type are presented in Table 2.

3.2 Spectral energy distribution analysis

We used the ARIADNE⁵ package described in Vines & Jenkins (2022) to perform a spectral energy distribution (SED) fitting of stellar models and with it determine stellar parameters from photometric data. The program cross-matches source coordinates over a wide range of catalogues to retrieve photometry and associated uncertainties; it then fits the SED across the available photometric bands using several stellar atmosphere models and employing Bayesian model averaging (BMA). With this approach, the final posterior is created by drawing from the posterior of the different models, weighted by their evidence. The parameters posteriors are then wider due to the systematic differences between the models, and error bars drawn from confidence intervals are a better representation of the method accuracy. Model averaging is done over PHOENIX (Husser et al. 2013), BT-Settl (Allard et al. 2012), ATLAS9 (Castelli & Kurucz 2003), and SYNTHE Kurucz (Kurucz 1993) models. ARIADNE then queries Gaia DR3 for parameters, and provides the best estimate for several other stellar parameters. Since our aim is to obtain an independent stellar characterisation, we used the default (weakly informative) priors on the parameters to fit.

The input photometry and model priors are detailed in Appendix A; the best-fit stellar parameters T_eff; logg [Fe/H]; distance, d; luminosity, L; radius, R; and mass, M are presented in Table 3.

3.3 Mass and radius computations from evolutionary models

We derived the stellar mass and radius using the Bayesian code PARAM⁶ (da Silva et al. 2006; Rodrigues et al. 2014, 2017). A grid of stellar evolutionary tracks was matched to the derived quantities of T_eff, [Fe/H] from spectroscopic analysis and luminosity from ARIADNE as presented in the previous subsection. The grid of stellar evolutionary isochrones was taken from PARSEC⁷ v2.1 (the PAdova and TRieste Stellar Evolution Code; Bressan et al. (2012). The median and the 68% confidence intervals (C.I. 68%) of the posterior distribution are used to derive the most probable values for mass and radius, presented in Table 4.

3.4 Rotational periods from photometry and chromospheric activity indexes

On top of the fundamental stellar parameters listed, an estimation of the rotational period is key for the characterisation of RV activity signals. We analysed the public ASAS-SN (Shappee et al. 2014) and TESS photometry (Ricker et al. 2014) to search for periodic signals that could represent the stellar rotation modulation of the four targets. We used the tpfplotter⁸ (Aller et al. 2020) and TESS-cont⁹ (Castro-González et al. 2024) to ensure that most of the photometric flux comes from the target stars. The TESS Simple Aperture Photometry (SAP) and the systematics-corrected Presearch Data Conditioned Simple Aperture Photometry (PDCSAP) were downloaded from the Mikulski Archive for Space Telescopes (MAST)¹⁰.

The ASAS-SN data were downloaded through the Sky Patrol web interface¹¹. The stars in our sample have relatively large proper motions (pm) of 1–3”yr, and the instrumental FWHMs are comparable to the ASAS-SN photometric apertures (i.e. ~2 pixels, with a pixel scale or 8”). To avoid flux losses, we downloaded the image subtraction (IS) photometry in short chunks based on the pm-corrected coordinates, similarly to Castro-González et al. (2023). Given the brightness of the stars, we also downloaded the ASAS-SN photometry extracted through the novel Machine Learning (ML) approach, which is aimed at improving the photometric quality of saturated sources (Winecki & Kochanek 2024). We ran the GLS periodogram on the four data sets (TESS/SAP, TESS/PDCSAP, ASAS-SN/IS, ASAS-SN/ML) for the four stars, and found no significant activity signals in any of them. This result is consistent with the low chromospheric activity previously reported in the literature (Sect. 2). The null results on HD 304636 are in line with the recent work of DiTomasso et al. (2025).

The rotational period of a star can be calibrated by the Ca II H+K chromospheric emission as measured by the $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ index. This relation was explored in the original study of Noyes et al. (1984b) and an updated version was presented in Mamajek & Hillenbrand (2008) that gained significant traction in the community. This study calibrates the Rossby number, R₀, as function of the activity index, and associates it to the rotational period, P_rot, via the relation R₀ = P_rot/τ_c, where τ_c is the convective turnover of the star. Since τ_c is a function of effective temperature, we then have ${\mathrm{P}}_{\text{rot\hspace{0.17em}}}(\log (R_{HK}^{\mathrm{\prime }}),{\mathrm{T}}_{\text{eff\hspace{0.17em}}})$.

Suárez Mascareño et al. (2016) calibrated photometrically determined rotational periods against $\log \left(R_{HK}^{\mathrm{\prime }}\right)$, considering 125 stars either as a monolithic set or as subsets of G, K, and M dwarfs. The authors estimated the accuracy of their calibration from the achieved dispersion, reaching values of the order of 20% for low-activity stars $(\log \left(R_{HK}^{\mathrm{\prime }}\right)<4.5)$. The spectral type of our stars are well defined from the study in Section 3.1 and Section 3.2. For HD 10700 and HD 102365, we used the calibration for K dwarfs, for HD 20794 the calibration for G dwarfs, and for HD 304636 the calibration for M dwarfs. HD 102365 was paired with the K dwarfs due to its late spectral type within the G and due to the low-number statistics of the Suárez Mascareño et al. (2016) G dataset. For $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ we use the value derived by Hojjatpanah et al. (2019) and available for the four stars. We considered an uncertainty 20% for the calculated rotational period, as used by these authors. Thus, we achieved a value of 37.7±7.5 d for HD 10700, 27.7±5.5 d for HD 20794, 36.4±7.3 d for HD 102365, and 30.3±6.1 d for HD 304636. We caution, however, that this determination depends strongly on the assumed spectral type, as explained in Appendix B.

3.5 Final stellar parameters and their accuracy

The stellar parameters derived from ESPRESSO spectra are compatible with those available in the literature (Table 1) and in particular with the ESPRESSO GTO preparatory work done in Hojjatpanah et al. (2019). The higher resolution of ESPRESSO and higher S/N of the spectra of this work allowed for slightly more accurate determinations. We note the match between ODUSSEAS and SteParSyn results for HD 304636; the difference in formal uncertainty between the two methods is due to ODUSSEAS representing accuracy variance while SteparSyn representing only the formal error. Since we are interested in accuracy we use the values and uncertainties delivered by ODUSSEAS.

A comparison between the spectral parameters derivation and the SED fitting results shows agreement (within formal uncertainties) for the effective temperature and metallicity. A comparison between SED fitting and isochrones show agreement in radii but not in mass. The luminosity values derived in our work for HD 10700 and HD 102365 also match the values derived by Teixeira et al. (2009) and Pepe et al. (2011) very well, and are well within the formal error bars.

For the analysis in the remainder of the paper, if a parameter can be determined via multiple methods, we selected the one that provides the highest accuracy. For the effective temperature, we chose values derived by spectroscopic analysis, for the luminosity the values from ARIADNE SED, and for the mass and radius the values derived using PARAM and PARSEC. We note, however, that the different quantities are determined at a different accuracy level, which, in turn, depends on the method chosen. In Tayar et al. (2022) the authors compared fundamental parameters obtained through different methods and concluded on the existence of a systematic uncertainty floor of approximately 2% in effective temperature, 2% in luminosity, 4% in radius, and 5% in mass. The systematic floor level on luminosity is smaller than the ones derived in our work and, thus, it is not expected to be a limitation. However, for effective temperature it is comparable or larger than the ones we derived; for mass and radii they are appreciably larger than derived by several of our methods and, thus, we might be limited by systematic errors. We only find evidence of such an issue with respect to mass.

For the stellar rotation period we consider the values calculated in Sect. 3.4. To our knowledge, there is no prior study on the accuracy of activity-indexes calibration in delivering rotational periods. The relative uncertainty of 20% on the rotational period used by Suárez Mascareño et al. (2016) is conservative when compared with other derivations and is obtained directly from the scatter of the calibration. It should, nonetheless, be used with care as a wrong spectral type or variable activity index can lead to different rotational periods in a way that cannot be included in the uncertainty value.

Using the stellar effective temperature, luminosity and stellar mass above, we followed Kopparapu et al. (2013) and Kopparapu et al. (2014)¹², to compute the conservative HZ limits for a planet with 1 M_⊕. The values calculated using these methods are summarised in Table 5 and used in the remainder of the paper.

4 ESPRESSO data, radial velocities and indicators

ESPRESSO is a high-resolution echelle spectrograph stabilised in temperature and pressure, and with mechanically fixed optics. The instrument operates on the ESO VLT and is physically located in the Combined Coudé lab. Light feeding is done from any of the four UT telescopes (or any combination of these) with a Coudé train and enters the spectrograph via two parallel fibres; the first points to the centre of the telescope field-of-view and the second points either to the neighbouring sky or to a simultaneous calibration source. The spectra from the two fibres are recorded simultaneously in interweaving orders. A pupil slicer splits the interference orders into two parallel slices; the spectra are imaged onto and collected by two detectors, located inside two independent cryostats. The complete wavelength coverage goes from 370 to 788 nm. This design allows for an extreme RV stability down to a target precision of 10 cm/s. We refer to Pepe et al. (2021) for further details.

4.1 Data acquisition, calibration, and pipeline reduction

Observations were acquired in high-resolution (HR) mode, with light fed by the 1″-diameter fibre, and delivering a resolution of ~140 000; the detectors were read without applying spatial or spectral direction binning (sampling 1x1). Given the brightness of the stars and thus a comparatively low contribution of readout noise to the error budget, we selected a FAST readout speed.

ESPRESSO’s wavelength calibration uses a Thorium-Argon lamp and a Fabry-Perot (FP) illuminated by a flat spectral source (Wildi et al. 2010). Combining the absolute reference from the atomic lines with the regular grid of FP lines enables the derivation of an homogeneous and precise wavelength solution across ESPRESSO’s wavelength range.

On top of precise wavelength solution methodology, ESPRESSO uses a simultaneous reference technique, described in detail for its precursor ELODIE in Baranne et al. (1996). For precise RVs the second fibre is illuminated by the FP to measure the instrument’s drift relative to the last wavelength solution; this provides for a way to quantify and correct for any instrument’s internal stability variation. This drift value is calculated independently for the two detectors.

The observations were acquired with a cadence of up to 10% of the star’s HZ inner edge period, as a trade-off between time investment and ability to detect planets up to the HZ. All observations were reduced with the ESPRESSO pipeline version 3.3.1¹³. The pipeline provides a complete and automatic reduction chain, from the raw files to the post-processing analysis such as the radial velocities and associated quantities. The main improvements of the major version 3 are: the calculation of time-dependent chromatic variation of the etalon spacing and potential lamp flux variations (Schmidt et al. 2022), along with improved masks weights that take into account the line’s full-width at half maximum (FWHM) to calculate the Doppler information content (instead of the contrast only). In version 3.3.1, the pipeline also accounts for the influence of the variable Sun-Earth distance and observer’s velocity in the barycentric correction computation (an effect with an amplitude smaller than 20 cm/s).

As an additional step, we calculated the parallax-dependent term of the barycentric correction to apply to RVs (see. e.g. Eastman et al. 2010). We did so by subtracting from the barycorrpy¹⁴ (Kanodia & Wright 2018) computation (using the Gaia DR3 parallax) the correction value calculated by the ESPRESSO pipeline. For reference, the value of this term is below 1 cm/s for HD 10700, peak-to-peak. Nonetheless, we added this differential barycentric correction term to the RVs, a posteriori, for the sake of completeness in correction¹⁵.

Finally, for each acquired spectra the ESPRESSO pipeline also produces an extracted spectrum corrected for telluric H₂O absorption, computed as described in Allart et al. (2022). This telluric-corrected spectrum is also provided with associated RVs and ancillary quantities.

4.2 Instrument interventions and technical issues

In June 2019 a new fibre feed was installed in ESPRESSO. While the instrument focus did not need an adjustment, we cannot exclude the possibility that the intervention introduced an offset in the RVs calculation. As such, we consider throughout the paper two data sets ESPR18 and ESPR19 corresponding to before and after the intervention, respectively The two optomechanical configurations are treated then as two different instruments, which introduces an offset parameter on our RV analysis. For more on the topic we refer to the equivalent situation on HARPS, described in Lo Curto et al. (2015).

Previous works, such as that of Faria et al. (2022), considered the ESPR21 dataset, corresponding to the data acquired after the COVID 19 observations interruption. This dataset was differentiated by the use of a different FP light source with a markedly different flux distribution than its predecessor. However, since pipeline version 3, calibration lamp flux differences are taken into account and absorbed by the improved wavelength calibration; a dataset separation is thus no longer necessary.

From January to April 2021, ESPRESSO atmospheric dispersion corrector (ADCs) were affected by a critical communication issue. The ADC prisms were not positioned at the correct configuration, which led to the atmospheric dispersion not being corrected. A comparison of RVs taken with the ADC operating correctly showed that the compromised dataset was affected by an additional scatter ≳1 m/s. Programs aiming at high spectral fidelity are also compromised by systematic errors on the flux and line-spread function determinations. The severity of these effects will vary from spectrum to spectrum, and without having access to the ADC’s true position a correction is not possible. We used a zero value of the header ADC2 keywords on right ascension, pressure and temperature to identify the observations affected, that were then discarded¹⁶.

Since we are aiming at the highest precision in RV, we restrict our analysis to spectra in which the pipeline general quality control, flux correction of atmospheric absorption, and drift correction were all correctly executed¹⁷.

We identified three nights corrupted by localised issues or events. On the 30^th October 2018 at 2UT a pressure gauge stopped working and there was a sudden increase of pressure within ESPRESSO vacuum vessel. This affected a series of HD 10700 observations taken around 3UT (MJD = 58421.65 d), that was discarded from our analysis. The pressure gauge was replaced on the following day and did not affect any other measures. On the 28^th July 2022, around 06 UT a block of HD 10700 observations was interrupted halfway, with the astronomer on duty signalling an interruption due to weather. A check on seismic activity logs of that day showed an earthquake with magnitude 6.2 and epicentre on Calama, only ~280 km from the observatory, triggering multiple aftershocks. We discard also this incomplete series of observations (MJD = 59788.7 d).

ESPRESSO was designed to centre the image of the star on the entrance of the fibre after the spectrally dispersed image is spatially corrected by the ADC. The operational limit of this device is at a maximum airmass of 2.2; observations at higher airmass will contain an uncorrected component, with the extreme blue and red wavelengths being slightly off-centred. Up to an airmass value of 2.3, this uncorrected term remains very small, with the maximum difference between the bluest and reddest star image centre of 0.15″. In this case, and for the specific case of ESPRESSO design, the loss in flux is of 1.5% in the blue and 0.5% in the red, and the RV effect is expected to be much smaller than 1 cm/s (Wehbe et al. 2020). As such, we choose to keep the spectra taken up to an airmass of 2.3 in the analysis, affecting only 12 individual integrations spread over 3 nights. However, we choose to discard a series of observation taken in the airmass range of 2.5–2.9, on the 7^th December 2020 at 6UT (MJD = 59190.25 d, fully outside the operation range of the ADC and for which the effect in RV can be larger than 10 cm/s, due to a considerable reduction in the flux of the bluest orders.

4.3 RV calculation

The RV used in this paper were calculated from the ESPRESSO pipeline extracted spectra, and followed two parallel methodologies:

• Cross-correlation function method (CCF), in which a mask representative of the spectral lines for a spectral type is cross-correlated with the measured spectra. For the original reference on the method and upgrade to weighted masks as a function of the spectral type, we refer to Baranne et al. (1996) and Pepe et al. (2002), respectively.

• Correlation with a template mask, where a template is first constructed by stacking a large number of spectra form the star, and the correlation is calculated between the individual spectra and this empirical template. We use the methodology described and demonstrated for ESPRESSO in Silva et al. (2022), named S-BART.

The specificities of the application of S-BART to ESPRESSO high-precision RV are given in Appendix C. For the following analysis, it is important to retain that by default S-BART discards wavelength regions contaminated by telluric features deeper than 1%.

On top of the standard S-BART run, we calculate the RVs using the spectra taken from either the blue or the red detector only. These chromatic RVs enable additional diagnosis for activity and the study of detector-based systematic effects. In the remainder of the paper the RV time series corresponding to the cross-correlation, S-BART with two detectors, S-BART blue detector, and S-BART red are labelled CCF, TM, TMb, and TMr, respectively.

The ESPRESSO spectra, CCF RVs and associated quantities were retrieved via the DACE interface¹⁸. To execute the DACE queries we used the arvi¹⁹ package, which allows for command line manipulation and data download, and performs convenient tasks such as secular acceleration correction using Simbad stellar parameters and to calculate the photon-noise weighed average RV. The complete time series for the four stars, including reduced and ancillary data is made publicly available to all via the DACE interface.

4.4 Photospheric line profile and chromospheric activity indicators

On the CCF one can measure directly the contrast, FWHM, and bisector inverse slope (BIS). The BIS was early on identified as a useful line profile indicator for tracking photospheric stellar activity (Queloz et al. 2001), and has seen extensive use. For stabilised spectrographs in which the instrumental profile can be assumed as fixed, the FWHM can be used as an indicator of stellar photospheric activity as well. This has been demonstrated for ESPRESSO on the works of Suárez Mascareño et al. (2020) and Faria et al. (2022) for the M dwarf Proxima.

The most extensively used chromospheric indicator is the $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ calculated from the Mount Wilson S index data and originally developed with the objective of studying long-term stellar activity (Noyes et al. 1984a). We used the ACTIN package (Gomes da Silva et al. 2018) to implement both $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ and H_α. For a detailed discussion of $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ as used in ACTIN, we refer to Gomes da Silva et al. (2021); for the H_α implementation and in particular the choice of wavelength bandwidth around the line, we refer to Gomes da Silva et al. (2022).

For the stars that have more than one spectrum per night, the activity indexes value is calculated independently for each spectra and the uncertainty-weighted average value taken as representative for the night, in parallel with what is done on RV.

5 Intra-night radial velocity variations

When observing very bright stars at a sub-m/s RV precision, it is common to acquire multiple consecutive short integrations. This allows us to cover the typical timescales of stellar pulsations and average out their RV signal, while keeping the photoelectrons on individual exposures below CCD saturation level. We employed this observational strategy for the G–K stars HD 10700, HD 20794, and HD 102365 studied in this work. We obtained between 7 and 15 spectra at very high S/N (ranging between 200 and 660), leading to photon noise RV uncertainties below 30 cm/s and well adapted to the search for low-mass planets at high RV precision. Observations of the fainter M dwarf HD 304636 employ a single integration per night, and the star is not discussed in this section.

The acquisition of several consecutive RV points each night allows for an evaluation of the RV precision on the short timescales they cover, ranging from minutes to around half an hour. In Fig. 1 we plot the distribution of photon-noise uncertainties and the observed intra-night scatter (as root-mean-square or rms) for the four different RV time series CCF, TM, TMb, and TMr for the three stars. We considered only nights when the complete number of exposures was acquired, when the sequence was not interrupted by weather or operational issues. The cumulative distribution functions are fitted by a cumulative log-normal distribution.

The time series CCF, TM, TMb, and TMr RVs display a very low photon noise, between 10 and 20 cm/s. The intra-night scatter is significantly larger and dependent on the star, being of ~40 cm/s for HD 10700, ~50 cm/s for HD 20794, and ~60 cm/s for HD 102365.

5.1 P-mode pulsations modelling

The main source of RV scatter for main sequence G-K stars on timescales of minutes is arguably the p-mode waves. These stellar pulsations are characterised by the quantities maximum pulsation frequency, ν_max, and large frequency separation, Δν, known to scale with stellar fundamental parameters (e.g. García & Ballot 2019) as $$\frac{v_{max}^{\mathrm{\prime }}}{v_{max}}=\frac{M^{\mathrm{\prime }}}{M}{\left(\frac{R^{\mathrm{\prime }}}{R}\right)}^{-2}{\left(\frac{T_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{\prime }}}{T_{\mathrm{e}\mathrm{f}\mathrm{f}}}\right)}^{-1/2},$$(1) $$\frac{\mathrm{\Delta }{\nu }^{\mathrm{\prime }}}{\mathrm{\Delta }\nu }={\left(\frac{M^{\mathrm{\prime }}}{M}\right)}^{1/2}{\left(\frac{R^{\mathrm{\prime }}}{R}\right)}^{-3/2},$$(2)

as a function of stellar mass M, radius R, and effective temperature, T_eff. The maximum pulsation frequency is simply the frequency of the strongest pulsation mode and the large frequency separation is, to a first approximation, twice the frequency separation between consecutive modes.

For RV pulsation signals, the mode amplitude envelope is well approximated by a Gaussian function and its FWHM scales linearly with large frequency separation (Lund et al. 2017). For our cases we consider FWHM_env = 10 Δν, in line with theoretical works such as Chaplin et al. (2019) (see their Fig. 2). Pulsation studies using photometry showed independently that FWHM_env = ν_max/2 (Campante et al. 2016), again in line with our assumptions.

The RV modes’ envelope has an amplitude, A_env, that depends on stellar luminosity, L, and mass, M, (see Samadi et al. 2007) as $$\frac{A_{env}^{\mathrm{\prime }}}{A_{env}}={\left(\frac{L^{\mathrm{\prime }}M}{L{M}^{\mathrm{\prime }}}\right)}^{0.7}.$$(3)

Dedicated asteroseismology campaigns have allowed for the accurate measurement of these asteroseismic quantities for main sequence stars, which we have used as our benchmark; we chose α Cen A and μ Ara, studied in detail in Bouchy & Carrier (2002) and Bouchy et al. (2005), respectively. The stellar parameters and average reference values are presented in Table D.1.

We simulated the p-mode signals contribution to the RV measured on each star; we use the simulation to calculate both the intra-night scatter generated and represented in Fig. 1 and the additional variance introduced on the nightly averaged RVs as it will be of interest to long-term RV analysis. We proceeded as follows.

For the star’s maximum pulsation frequency, ν_max – the strongest p-mode – we created a sinusoidal RV signal of amplitude A_env with a random phase. Then, for the neighbouring modes of frequencies ± Δν/2, we calculated the mode amplitude from the envelope value at the corresponding frequency and simulate two sinusoidal signals, each with a random phase. We repeated this procedure for increasingly higher-order modes m of frequency ± m Δν/2 stopping when the mode amplitude as drawn from the envelope is smaller than 1 cm/s.

The RV signal of the calculated modes was co-added and its value projected on a temporal grid representing the detector integration sequence for each star in one night, with the integrations separated by the detector readout time and the individual photon noise drawn from the observed distribution median e_ph.noise. The simulations were then repeated 50 000 times to obtain the same number of individual sequences, or nights. For each sequence, we calculated the RV offset introduced by the pulsation signal on this particular observation as the RV weighted mean and the intra-night dispersion from the intra-night rms. The dispersion of the RV offset, e_puls, is representative of the uncertainty introduced by the pulsation on the average value, while the median of the intra-night dispersion distribution, σ_puls, can be directly compared to the observed one. For visual representations of key steps of the simulation, we refer to Appendix D.2; namely, we depict how the frequencies are drawn, an example of pulsations during one night and a distribution of scatter and distance to the mean when repeating the process a large number of times.

For the stellar parameters we use the values presented in Table 5; the measured and simulated quantities are presented in Table 6, namely the distribution median photon noise, e_ph.noise, the measured scatter, σ_meas, the final uncertainty introduced by pulsations on the nightly average, e_puls, and the dispersion created by pulsation on the intra-night set σ_puls. The results are found to be repeatable and thus consistent down to 1–2 cm/s.

The first takeaway is that e_puls is of 5–10 cm/s. This is expected, in line with the common assumption that averaging the p-modes over a timescale of a few pulsation periods is enough to correct for their contribution to the nightly averaged RV. We show here it is also achieved at a precision of better than 10 cm/s.

The simulations predict a scatter σ_puls that follows the observed σ_meas. The pulsation-induced scatter for HD 10700 matches the observations to 1–2 cm/s, and for HD 20794 to 4–5 cm/s. For the earlier spectral type (G5) HD 102365 there is a larger and clearly non-negligible difference of 10 cm/s between predicted and observed values. ESPRESSO is expected to have a RV precision of better than 10 cm/s over one hour, and if instrument stability was the source of scatter it would be present in all stars, having been easily detected for at least HD 10700. The discrepancy can be attributed to poor stellar parameter determination or an additional source of scatter.

Fig. 1 Distribution of the photon noise uncertainty values (top) and observed nightly scatter (bottom) for the stars HD 10700 (left), HD 20794 (centre), and HD 102365 (right). The solid-line histograms represent the data and the dashed lines the cumulative log-normal distribution fits.

5.2 Granulation signal modelling

The convective phenomena on stellar photospheres are commonly referred to as granulation. Contrary to pulsations, the physical mechanisms governing their RV signals remains poorly understood. Solar observations enabled the measurement of the RV signal of convection cells, or granules, and of the coordinated motion of groups of granules, often know as super-granules (Simon & Leighton 1964), which create an additional RV component know as super-granulation. While the granules lifetime and renovation time has long been estimated at ~10 min (e.g. Bahng & Schwarzschild 1961; Giannattasio et al. 2013), super-granulation timescales are still the object of active research, with Sowmya et al. (2023) showing the values range from ~24 h for the quiet to ~54 h for the active regions. For a review on super-granulation, including studies on timescales, we refer to Rincon & Rieutord (2018).

There has been significant focus on the photometric observation of granulation on giant stars, as their extended convective region and higher brightness allow for a more accurate signal characterisation than on main sequence stars. A star’s power spectrum density (PSD) is fitted with Harvey-type functions PSD(ν)=A/(1 + (τ ν)^α)), in which A is the amplitude, τ is the characteristic timescale, and α is the power-law slope.

In the literature one can find different implementations. The photometric granulation signal is represented as a sum of up to three Harvey functions, and different functional forms are used to represent the pulsation signal (Gaussian vs. Lorentzian), with some authors choosing not to fit the pulsation-dominated frequencies. The work of Mathur et al. (2011) reviewed the different operational choices applied to Kepler data and concluded that the results are qualitatively similar, but the quantitative differences still leave room for different interpretations on the underlying physics. The authors showed then that from basic principles one has τ ∝ 1 / ν_max, and found evidence for the relationship between the power of granulation and ${\nu }_{max},P_{gran}\propto 1/v_{max}^2$, as proposed by Kjeldsen & Bedding (2011). The Kepler mission data is used to fit the slopes of these dependencies as $\tau \propto {\nu }_{max}^{-0.89}$ and $P_{gran}\propto {\nu }_{max}^{-1.90}$, but the residuals around the fitted models exhibit non-negligible scatter.

Comparatively, there have been much fewer works attempting to characterise granulation on RV; we summarise these in Appendix E. For our simulations we follow the PSD characterisation of Al Moulla et al. (2023) on Sun RV data. The authors used two components and fix the slope α at 2. Since the importance of scaling the parameters is still a matter of debate, we performed a first set of simulations with the quantities measured for the Sun and for a second we applied the scaling relationships on Eq. (1) to derive the amplitude and granulation timescales. For the sake of simplicity and in accordance with the literature, we refer to the two components as granulation and super-granulation.

The RVs were calculated as in Al Moulla et al. (2023), decomposing the power spectra into frequencies, ν, with width, Δν, that generate signals with amplitude $\sqrt{PSD(\nu )\mathrm{\Delta }\nu }$, and randomised phase. The integration limits are set by the limits of our knowledge: since we cannot characterise granulation on timescales longer than 100 d, our lowest frequency was set at 0.1 μHz; given we cannot sample signals with higher cadence than 1 minute (our Nyquist frequency), we set the upper limit at 6 mHz. The integration is done in steps of 1 μHz. For the sake of simplicity, we perform the simulations for CCF and TM only. The results for the effect of the granulation on the long-term uncertainties, e_g, and on the nightly scatter, σ_g, are presented in Table 7.

The results with the solar and scaled values are similar, with the scaled values showing a larger variation across stars, up to 5 cm/s. The granulation component introduces an uncertainty of 19–28 cm/s on the nightly averaged RVs; on the other hand, the super-granulation component introduces an uncertainty of 58–71 cm/s. Either component introduces an intra-night scatter of 15–20 cm/s.

With this model, the granulation or super-granulation components impact on the intra-night scatter remains of around 25 cm/s across stars, and scaled / unscaled modelling choices. When one discounts the effects of pulsations, this estimation is in excess of the values measured for HD 10700 but remains compatible with the scatter measurements on the other two stars. However, the effect of super-granulation, of 65–70 cm/s, is not averaged out by our observation strategy and should be easily detectable as a parasitic effect on the night-to-night RV variation.

6 Nightly averaged time series

We perform an analysis of the nightly averaged RV time series of the targets HD 10700 and HD 102365, and on those of HD 304646, for which a single observation per night was taken. The analysis of nightly averaged time series of HD 20794 was presented in Nari et al. (2025).

6.1 Dataset properties and simple statistics

The RV time series are shown in Fig. 2 and simple statistics listed in Table 8. For each star, dataset, and RV timeseries we list the number of nights and observing time span, plus the average RV uncertainty (calculated only from the photon noise) and weighted rms of the nightly averaged RVs.

On HD 10700, CCF TM and TMb show an average uncertainty of 3 cm/s, with TMr being larger at 5 cm/s. The weighted rms is of 67–70 cm/s for the different RV calculations; the similarity between TM, TMb and TMr suggests that the dominant source of scatter is achromatic.

For HD 102365 there is a clear difference between the ESPR18 and ESPR19 datasets. While the average photon noise is similar, the rms of the nightly averaged RVs increases significantly from ESPR18 to ESPR19. In contrast, HD 304636 ESPR18 shows a significantly higher rms than ESPR19. For both stars TMr shows a significantly lower scatter than TMb. This colour-dependent rms points towards an activity-induced flux-dominated RV variation, that is known to be mitigated at redder wavelengths, where the photometric contrast (or flux) effect between active regions is smaller. A variation of the instrument stability is excluded by the contemporaneous RVs of HD 10700, that show a scatter of only 60–70 cm/s over approximately the same period.

It is worth noting that for the two G-K stars HD 10700 and HD 102365 the effect of pulsations on the nightly average was estimated to be slightly larger than when considering photon noise alone: for HD 10700 we estimate the uncertainties to increase up to 5 cm/s, and for HD 102365 up to 8 cm/s.

The temporal evolution of the stars in RV does not define a clear shape. In Fig. 2 for HD 10700 we have what seems to be a periodic signal in the relative variation Δ(TM - CFF), which we discuss in Sect. 7.2. The temporal evolution of the photospheric and chromospheric indicators is presented in Fig. 3. Their scatter is much larger than their formal uncertainties. The line asymmetry and width variations measured by the BIS and FWHM are not expected to vary for observations of very quiet stars with RV-dedicated spectrograph, that were built to have stable instrumental profiles. As such, the variability shown could actually represent an activity effect impacting the RV.

Lovis et al. (2011) had already noted that long-term RV trends caused by the stellar magnetic cycles would be detected as opposite sign trends on the FWHM and the contrast indicators; this might be at play here. However, a very slow monotonic variation in resolution over the years²⁰ can explain small-amplitude indicator variations as seen on HD 10700.

The $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ chromospheric indicator is known to show variability even for the lowest activity stars. The Sun shows a variation in $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ of 0.05 due to rotational modulation during its low-activity phase (e.g., Milbourne et al. 2019; Maldonado et al. 2019a), and low-activity solar-like stars, like the ones targeted commonly by RV searches, show similar amplitudes (Baliunas et al. 1985). Gomes da Silva et al. (2021), showed that even HD 10700, one of the lowest variability stars within their HARPS sample, has a variability of 1.6%, that corresponds to 0.017. This matches well the ESPRESSO observations presented here.

Overall, the long-term (≳ 100 d) evolution of FWHM, Contrast, BIS, and $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ on different stars and with different slopes is suggestive of a long-term activity evolution of the stars over the 4 years of our observations.

6.2 Presence of periodic signals

We computed the Lomb-Scargle periodogram (LS, Lomb 1976; Scargle 1982) in its floating-mean form (known as generalised Lomb-Scargle, GLS) as introduced by Zechmeister & Kürster (2009) and implemented in astropy. The false-alarm probability (FAP) is computed as defined on Baluev (2008). We chose three FAP threshold probabilities of 1%, 0.1% and 0.01%, and a maximum frequency of 1 d⁻¹ for the GLS and FAP calculations. The LS periodogram for the different RV sets is shown in Fig. 4, along with the different FAP. Given the observations time span is ~1500 d we consider the longest period that can be confidently detected is of the order of 1000 d; variations longer than this value are better modelled as linear or quadratic slopes and a frequency-oriented discussion is not particularly insightful.

For HD 10700 the different RV time series exhibit different significant peaks: CCF at 1.13 d, TM at 1.05 d and 1.13 d (the three at (1% FAP) and TMr a doublet at 20.07 and 21.23 d (reaching 0.1% FAP). TMb shows a significant peak at 1500 d (between 1% 0.01% FAP). No significant signals are detected on HD 102365 or HD 304636.

The peaks of the LS periodogram do not represent directly the best-fit matches to a sinusoidal signal. The amplitude of the periodogram at a given frequency depends not only on the signals present at this frequency but but also on signal spectral leakage, that depends directly on the temporal sampling. This effect is represented by the Window function, estimated by setting the measured RV signal at 1.0 (see VanderPlas 2018); the window functions for the different stars is shown on the bottom row of Fig. 4. The strongest individualised peak is at approximately one year (360 d for HD 10700, 350 d for the other two stars), only topped by the broad signals with periods longer than 2000 d. The window functions of HD 102365 and HD 304636 also show strong power at 1000 d. It is worth noting the doublet seen in HD 10700 TMr are aliases if one considers a sampling frequency of approximately year.

The FAP is known to be a reliable estimator, as long as the uncertainties and the frequency grid are correctly defined (see e.g. Delisle et al. 2020, for a discussion). Unfortunately, so far we have only considered photon-noise uncertainties for the error budget. If we repeat the GLS and FAP calculations adding quadratically 10 cm/s to the individual photon noise error bars the signals of HD 10700 become non-significant, and the signals at 310–330 d in CCF, TM, and TMb approach 1% FAP. This serves as a word of caution on the interpretation of the frequency analysis.

The LS periodogram evaluates the least-square fit of a single sinusoidal function to the data and assumes that the formal uncertainties describe the variance of the dataset. These two assumptions can be lifted by considering the ℓ₁ periodogram (Hara et al. 2017), that uses sparse recovery tools to represent a RV times series as linear combinations of a small number of sinusoidal functions. The tool also allows us to define the noise model via a covariance matrix V, and we adopt the representation with element V_kl of the form $$V_{kl}={\delta }_{kl}({\sigma }_k^2+{\sigma }_W^2)+{\sigma }_R.e^{-\frac{{(t_k-t_l)}^2}{2{\tau }^2}},$$(4)

where (on top of the individual nominal uncertainties, σ_k), we consider a white noise, σ_W, common to all data and a red noise with standard deviation, σ_R, and a timescale, τ.

To solve the basis pursuit model we used the default choice of the least-angle regression (LARS)²¹ algorithm. Significance is estimated with the FAP, defined for the sequence of peaks, namely, the first value evaluates the presence of the highest amplitude signal, the second value evaluates the presence of the two highest amplitude signals, and so on, all calculated relative to the null hypothesis of no signals present in the data. We used ℓ₁ on our RV datasets, restricting our analysis to periodicities longer than 1 d.

In a first analysis we considered a white noise component of 10 cm/s. The resulting periodograms are shown in Fig. 5, with the ten strongest peaks represented. The calculated FAP values for the three stars and different RV calculation methods are listed in Appendix F.

The detected signals have either a FAP larger than 1% or a period very similar to our temporal baseline. For HD 10700 CFF show a 3% FAP signal close to one year, while TM and TMb show a significant signal at 1423.5 d. A signal with 20 d in TM and TMb and a very similar signal in TMr show a FAP of a few percent. For HD 102365 there is no formally significant signal and for HD 304636, we would need a combination of 10 (CCF), 6 (TMb), or the 10 peaks (TMr) to reach an FAP of 1%.

The forest of peaks on the ℓ₁ periodogram is characteristic of noisy datasets that do not have a sparse representation in frequency, as discussed in Hara et al. (2020). In such a case, we should evaluate the impact of the noise model and to test models including correlated noise. We followed the cross-validation methodology described in Appendix B.1 of Hara et al. (2020) and compared models with white noise σ_W of 10, 20, 40, and 60 cm/s combined with σ_R with a value of 0 (no red noise component), 20, and 60 cm/s and timescales of 0.02 and 0.54 d. These amplitudes and timescales correspond roughly to the two timescales and amplitudes of granulation measured on the Sun. An average of the 20% of the models with the highest cross-validation ranking still led to FAP that were, at best, much larger than 1%; the only exception were the already detected signals at 1423.5 d for HD 10700.

In summary, there is only one periodic signal detected in RV with a period below 1000 d. It is present in the different template matching RV time series of HD 10700 and has an approximate period of 20 d. The signal is statistically significant when measured with an LS analysis but only if photon-noise error bars are considered, and with varying significance between 1% and 0.05%, depending on the RV reduction. On ℓ₁ the signal is present but with FAP of a few percent, and cannot be considered significant. The other signals detected are a consequence of the window function effect.

Fig. 2 RV time series using the different RV calculation methods CCF (top) and difference relative to it of TM (center), TMb, and TMr (bottom) for the stars HD 10700 (left), HD 102365 (center), and HD 304636 (right). ESPR18 data are represented by open circles while ESPR19 data are represented by filled circles.

Fig. 3 Time series for the line profile and chromospheric indicators for the stars HD 10700 (left), HD 102365 (centre), and HD 304636 (right). ESPR18 data are represented by open circles, while ESPR19 data are represented by filled circles.

Fig. 4 Lomb-Scargle periodogram for the different methods of RV calculation and associated window function for the stars HD 10700 (left), HD 102365 (center), and HD 304636 (right). The FAP levels of 1, 0.1, and 0.01% are shown as horizontal grey lines.

Fig. 5 ℓ1 periodogram for the stars HD 10700 (left), HD 102365 (centre), and HD 304636 (right). The ten strongest peaks periods are identified in the image; the associated FAP is listed in Table F.

6.3 Activity indicators study: periodicity and correlations

We perform an LS periodogram analysis for the indicators as previously done for RVs. The results are presented in Fig. 6. By construction, the window function for the indicators is almost identical to that of the RV data and was therefore omitted.

For HD 10700, the different photospheric indicators show peaks at common periods: 307 d, and 440 d, and a wide peak from 1250 d approximately (<0.01% FAP). The $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ indicator has a peak at 230 d (1% FAP) and a much wider and stronger peak from 550 to 1050 d (<0.01% FAP).

HD 102365 shows highly significant power for periods longer than 1380 d in FWHM (<0.01% FAP) and between 570 and 690 d in Ha₀₆ (between 0.1% and 1% FAP). While never significant at 1%, it is interesting to note all indicators show a peak at 125 d; $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ even shows a doublet at 119 d, which are the aliases pair for a sampling frequency of 1000 d.

For HD 304636, we have a peak in $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ and Ha₀₆ at 437 d (0.1% and 1% FAP, respectively). The FWHM and $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ show significant power above 1500 d.

The presence of significant signals on a timescale similar to the observing span motivates a linear fitting on the indicators and analysis of the residuals. For HD 102365 on $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ and Ha₀₆, we detect a significant signal at ~123 d, very similar to the announced planet period at 122 d. On HD 3046336, the period of 30.3 d is detected on $\log \left(R_{HK}^{\mathrm{\prime }}\right)$, close to the expected rotation period. We checked for correlations between nightly averaged RVs and indicators using the Pearson correlation coefficient and display these values in the radar chart in Fig. 7.

The correlation coefficients absolute values for HD 10700 and HD 304636 are below 0.4, showing weak correlations, if present at all. On the other hand, for HD 102365 the coefficients values range from 0.4 to 0.8. The only exception is Ha₁₆, that is arguably poorly adapted to this spectral type. Interestingly, for CCF, TM, and TMb the photospheric indicators FWHM and contrast show a correlation of 0.7. The explained variance is given by R², that in the simplistic case of a linear model is given by the square of the coefficient of correlation. In this approximation, the correlation coefficient value of 0.7 sets an important threshold: half of the variance is explained by the correlation. We have then clear evidence that HD 102365 RVs are contemporaneous with activity.

The study of periodicities of RVs and their correlations with activity indicators is insightful but incomplete. Activity-induced signals may not be periodic, and remain undetected in a periodicity analysis; they can also show a time (or phase) lag relative to the associated RV signal and thus lead to a null correlation coefficient while being physically connected (see e.g. Collier Cameron et al. 2019, for an example on the Sun). For a complete analysis one must model simultaneously activity and potential planetary signals.

Fig. 6 Lomb-Scargle periodogram for the activity indicators for the stars HD 10700 (left), HD 102365 (centre), and HD 304636 (right). On the upper panels we have the photospheric indicators FWHM, line contrast, and BIS, and on the bottom on we have $\log (R_{HK}^{\mathrm{\prime }})$ and the two Ha parametrisations. The FAP levels of 1, 0.1, and 0.01% are shown as horizontal grey lines.

Fig. 7 Radar chart of the correlation between nightly averaged RV and each of the indicators, for the four RV calculation methods and three stars. The different colours represent the different RV datasets, with filled or open dots representing positive or negative signed Pearson correlation coefficients.

7 Radial velocity analysis with kima

We used the kima package²² (Faria et al. 2018) to evaluate the number of planetary signals present in the data while considering a realistic modelling of activity. kima uses the diffusive nested sampling (DNS) algorithm from Brewer et al. (2011) to sample from the posterior distribution of the model parameters; this method delivers an estimate for the marginal likelihood for each model and enables direct model comparisons (see e.g. Brewer 2014; Feroz et al. 2011). In this context, the number of detected planets N_p is a parameter that can be left free or fixed for the analysis. For each model we obtain at least 50 000 effective samples from the posterior, which allows for a robust sampling of the different parameters.

7.1 Modelling with a quasi-periodic GP

kima is able to make use of a Gaussian process (GP) with a quasi-periodic kernel as a noise model to accurately represent activity-induced signals; the hyper-parameters η₁-η₄ of the GP are inferred together with the planetary orbital parameters, star background velocity, and instrument-specific noise and offset values. η₁ is the amplitude of the GP quasi-periodic kernel, η₂ is the evolution timescale, η₃ the rotational period it tries to represent and η₄ the signal’s harmonic complexity. This activity signal representation has had its origin in photometric activity modelling and has seen a large number of applications and success on RV modelling. For an early implementation see Haywood et al. (2014); for a recent revision work see Nicholson & Aigrain (2022). The choice of priors for each of the parameters is discussed in detail in Appendix G.1 and listed in Table G.1.

The best-fit parameters for the different stars and RV calculation methods is presented in Table G.2. Using the quasi-periodic kernel and the previously mentioned priors we could reproduce a rotation-modulated activity signal with a period close to the expected P_rot value, as listed in Table 5.

Letting the number of planets free between zero and one, no significant planetary signals were detected on any of the stars using any of the RV calculation methods. We repeated the analysis fixing the number of planets at zero to obtain the best characterisation of stellar activity; the posteriors for the case with variable planet number and fixing planet number at zero are almost identical. For a given star, the results across RV calculation methods are also very similar; the only exception are for η₃ and η₄ that for the CCF RV of HD 10700 are very poorly constrained. The distribution of the key parameters for the TM RV calculation methods, with and without planets, are shown in Figure 8.

The modelling on HD 10700 can be interpreted as successful for several reasons. The best-fit rotational period η₃ is of approximately 42 d, halfway between the 38 d P_rot value estimate from spectroscopic calibrations and the 46 d value from interferometric studies (see Sect. 3). Not only a rotation period but a first harmonic are detected, and an evolution timescale η₂ maximum is found at 2.5 times the rotation period, compatible with what is observed on the Sun (Collier Cameron et al. 2019). Moreover, the jitter on the long timespan dataset ESPR19 and η₁ parameter distributions are approximately normal and have values of 45–50 cm/s. However, η₄ is of 2.2 for the CCF dataset, which can lead to a degeneracy in the activity signal characterisation (as discussed on Appendix G.1). The RV offsets for the three stars are of the order of m/s, compatible with zero, and mutually compatible across RV calculation methods.

The corner plots are presented on Appendix G.5, and show correlation between different parameters. The linear slope and quadratic term of the stellar systemic trend are mutually correlated and correlated with the instrumental offset; the quadratic term is correlated with the systemic velocity. The jitter parameters for the three stars are slightly correlated with their GP parameters η₂ and η₄, with the correlations being less pronounced for HD 10700. The correlation between jitter and η₂ / η₄ occurs when the jitter is very small and the timescales very short (with η₂ similar to the rotation period), and η₄ is at its lowest limit, creating a signal with very low complexity.

The classification of a detection as significant depends on the (Keplerian) planet model being favoured, which in turn depends on a chosen metric crossing an a-priori threshold. For exoplanet detection within the Bayesian Framework, the metric is typically the Bayes factor and the threshold 150 (Kass & Raftery 1995). While no star shows a significant detection, the posterior on the number of planets of HD 10700 has a similar number of samples for N_p = 1 than for N_p = 0. The posterior of the orbital period and the joint posterior of semi-amplitude period and eccentricity period for TM are shown in Fig. G.1. There is a pile-up of samples at periods of ~1 d and ~20 d; when looking at the semi-amplitude, we see that these extend from 15 to 60 cm/s, with the eccentricity extending up to around 0.6. This cluster of posterior parameters does not correspond then to a well-characterised and well-defined signal, but are very variable parametrisations of tentative detections. The posteriors for the CCF RVs are similar, with a third clustering at a period of one year.

Fig. 8 Normalised posterior distribution of the two instrument jitter and four GP parameters as estimated by kima, for the three stars studied in the section and the TM RV calculation method. The distributions of the posteriors for GP-only (dashed line) and up to one planet runs (solid line) are both shown.

7.1.1 Changing the priors: Highly constrained planet searches

The comparison between planetary model and non-planetary model depends on the assumed priors. In particular, an orbital period constrained to regions of high likelihood (e.g., periodogram peaks) can boost the posterior evidence for the planetary model when compared to a wider parameter search. Whether the narrow priors have a physically sound justification remains a separate question.

To test the presence of planets in this skewed scenario we consider three cases of narrow searches around 20, 122, and 7.57 d for HD 10700, HD 102365, and HD 304636, respectively. We consider Gaussian priors around these periods, with a standard deviation of 20% of the chosen period value. For HD 10700 we choose the TM and for HD 102365 and HD 304636 we choose the TMr datasets as optimistic cases. The results remain largely unchanged; there is no significant planet detection even in this scenario.

7.1.2 Changing the model: Quasi-periodic joint kernel with FWHM

We tested a joint modelling of RV and FWHM, as used in Faria et al. (2022). The FWHM is known to track activity changes in a wide range of scenarios, but in our study it only shows a significant correlation with RVs for HD 102365. None of the three stars has an RV or indicator time series with a significant periodicity at the expected stellar rotational period. This model is tried thus for purely exploratory reasons, without having any a priori support from data inspection.

In our joint model the parameters η₂ and η₃ are shared between the RVs and FWHM variation, as a way of anchoring the rotational period and the evolution timescale. The parameter η_FWHM,1 is drawn from a log-Uniform between 0.1 and the maximum FWHM variation and the parameter η_FWHM,4 is drawn from the same prior as η₄. The rest of the analysis remains unchanged.

The joint model results are presented in Table G.2, and the posterior distributions are very similar to those seen in the previous section. For HD 10700 the only noticeable difference occurs for η₃, that now has a very skewed distribution towards small values. For HD 102365 η₁ and η₄ show larger values than seen in the previous section; η₄ being larger than 2.0 is associated with a higher model degeneracy. For HD 304636 the posteriors are compatible, but the parameters η₂, η₃ and η₄ show narrower posteriors, indicating this model provides increased constraints.

7.2 Telluric correction and its impact on RVs

Telluric water absorption features are present across the optical domain and as such cannot be masked out without a considerable loss in RV precision. With a small relative depth (<2%) these are often referred to as micro-tellurics, and are expected to have an impact on RVs that is larger than 10 cm/s over a wide range of conditions (target spectral types, airmass, precipitable water vapour content on the atmosphere), reaching up to 80 cm/s on extreme cases. For more on the topic, we refer to Cunha et al. (2014).

Allart et al. (2022) presented a telluric absorption correction model of ESPRESSO and tested it on a smaller dataset on HD 10700 than present here, identifying an effect on RVs that was of the order of 10 cm/s rms, could reach a maximum value of 58 cm/s, and had a significant impact on the detection of signals with a period close to one year.

Since CCF and TM are subject to different telluric masking criteria, water micro-tellurics are a prime candidate for the modulation seen in Δ (CCF - TM) (see Sect. 6.1). To test for this hypothesis and attempt a gain in RV precision we apply the CCF and S-BART RV computation to telluric-corrected spectra of HD 10700.

We tested different telluric masking conditions on S-BART; on top of the baseline choice of masking regions affected by tellurics deeper than 1%, we alternatively masked water features to deeper than 25%, 40%, and 50%. The masking of non-water features remained fixed at 1%. Since there is a larger density of water features in the redder optical wavelengths, the scatter of TMr is the most efficient metric for an improved correction. We notice that for the 25% threshold there is a reduction in 6–8 cm/s, and a negligible reduction for deeper thresholds. For 25% one also has a reduction of 2–3 cm/s rms in the TM. This is in line with the correction being at its most effective for shallow water lines, as intended. We keep then this threshold as our nominal S-BART correction threshold.

The night scatter for the telluric corrected spectra is presented in Table 9, and the rms can be compared directly with that of Table 8. The scatter on CCF_tc remains the same when compared with the uncorrected one very probably because the correlation masks were not extended in wavelength to make use of the now corrected wavelength range. An LS of the uncorrected-corrected CCF RVs and indicators show clear peaks at one year and fractions of it (Fig. G.2), illustrating nonetheless the positive impact of the correction.

The kima analysis using the same model and priors as in Sect. 7 on the telluric corrected RVs is presented in Table G.3. The telluric correction allows us to reach lower ESPR19 jitter values for the different TM runs, between 41 and 45 cm/s.

7.3 Compatibility limits

When running kima one can fix N_p to the number of detected planets plus one to calculate the parameters of undetected planets that are compatible with the data (see Standing et al. 2022). The samples with the largest planet semi-major axis amplitude, K₀, in a given period bin will yield the most massive planet compatible with the data and with the priors. The compatibility limits produced in this way have important advantages when compared with the traditional detection limits from the injection-recovery method: the different orbital parameters are marginalised by construction and the results do not depend on signal injection properties. We calculate and compare the compatibility limits for the different RV time series in Fig. 9; for HD 10700 and HD 102365 we mark the parameters of the planets announced in the literature.

For a given star the results for the different RV calculation methods are similar but not identical. For HD 10700, as we increase in period we have a first bump (increase in equivalent planetary mass) at around 20 d, the value of the posterior pileup we saw in Sect. 7.1, and the period of an announced planet. For periods longer than 300 d the CCF RV calculation method shows a lower detection capability than the others, unfortunately for periods covering the star’s HZ. However, using the telluriccorrected version of the RVs corrects for this effect, creating an almost linear limit between 100 and 1000 d. While the gain is the most pronounced for CCF, it is also present for the other reductions. For HD 102365 and HD 304636 TMr provides the lowest (i.e., most sensitive) limits. For HD 102365 this is interpreted as the red wavelengths being less sensitive to star’s activity. The different RV time series and different models yield slightly different compatibility limits, with the bumps showing a different amplitude; this points to these being created by the spurious fitting of activity signals. For HD 304636 we obtain very similar detection limits when comparing pairwise CCF − TMb and TM − TMr. We can speculate that this is probably due to the template matching being more efficient at extracting RVs from M-dwarfs, that have most of the spectral information on the red detector, with the CCF being more similar to the information on the blue detector. This could be the effect of activity at play, leading to slightly lower limits measured with TMr than TM. For this star we have a clear bump at 30–40 d, the expected rotational period of the star.

For each star we need to choose an RV calculation method to define the most representative compatibility limit. For HD 10700 we choose TM as is the one that best constraints the GP parameters, while allowing for a lower detection limit to be reached up to the HZ. For HD 102365 and HD 304636 the lowest detection limits are obtained with TMr, that we thus select. A summary of the compatibility limits mass for key periods is shown in Table 10. For HD 10700 we are below 1 M_⊕ for periods shorter than 10 d and at 1.7 M_⊕ for periods shorter than 100 d; within the HZ we are sensitive to planets with a minimum mass above 2.5 to 5 M_⊕, approximately.

Fig. 9 Compatibility limits for the three stars studied, HD 10700 (top), HD 102365 (middle), and HD 304636 (bottom). The maximum semi-major RV amplitude of the samples is plotted for each of the RV calculation methods model, and telluric correction set; HZ limits as presented in Table 5 and orbital properties of announced planed described in Sect. 2 are shown.

7.4 RV stability and precision floor

Any RV variability not parametrised by η₁ − η₄ will be captured by the jitter term of each dataset. While the dataset ESPR18 spans only 150 days, ESPR19 spans almost 4 years, with the data clustered into the months when a star is observable from Paranal observatory.

For a further inspection of ESPRESSO stability and the RV precision limit of our model we take HD 10700 telluric corrected data and subdivide ESPR19 into four yearly datasets ESPRy19 − ESPRy22, with 20 to 33 points each, timespans in the range 70–230 d, and temporal gaps between subsets from 120 to 270 days, approximately. We use the previously described quasi-periodic kernel and priors within the GP to reproduce the RVs derived with the different methods. The ad-hoc separation leads to the derivation of five jitter values and four offset ones; their posteriors are presented as box plots in Fig. 10. The dataset properties and measured values are presented in Table 11.

The result from the different RVs is, once again, compatible; we focus our analysis on the results for TM. While the jitter for ESPR18 remains at 73 cm/s, for ESPRy19 − ESPRy22 this is of 32, 38, 23, and 57 cm/s, respectively. The datasets have different number of points and span, with ESPRy20 being particularly small and short; the results on it can be affected by small-number statistics. Nonetheless, the different sub-datasets show variable jitter, up to a factor of 2 to 3.5, depending on the reduction. The reason for the difference is unclear, and can come from either a variation in instrument stability or induced stellar activity. The different reductions have a minimum yearly timescale jitter between 20 cm/s and 30 cm/s.

Fig. 10 Boxplots for the posteriors of the jitter (top) and offsets between datasets (bottom) for the different RV calculation methods with five datasets of telluric-corrected RVs of HD 10700.

8 Discussion

8.1 A forward model of pulsations

The asteroseismology of main sequence stars experienced significant developments over recent decades. We refer to García & Ballot (2019) for a comprehensive review of the topic. Our work uses simple principles and scaling equations to estimate the contribution of p-mode pulsations to averaged intra-night RVs and intra-night RV scatter. We take into consideration the time intervals covered by individual exposures and the detector readout-time gaps in the acquisition of consecutive integrations, as used in RV searches, and use individual p-mode frequencies instead of an envelope to represent the ensemble of oscillations. These two aspects are certainly behind the very good match with observations, and key differences between our others such as Chaplin et al. (2019).

Some important shortcomings apply. There have been reports of significant correlations between surface activity indicators and solar-like oscillations properties for the Sun, with the oscillation frequencies increasing during higher activity periods. The commonly accepted explanation is that magnetic fields induce perturbations in the outer layers of the star, influencing the oscillation cavity of the p-modes. The exact mechanism behind this frequency shift remains unknown; see Howe (2008) and references therein. In addition, several other p-mode characteristics such as amplitude and widths change within an activity cycle (e.g. Salabert et al. 2011). A similar mechanism may be at play and detectable at very high RV precision.

Asteroseismology enables the measurement of stellar parameters with exquisite accuracy. We could compare our own parameter estimation against those provided by asteroseismology runs for HD 10700 and reached an excellent agreement. It was precisely on this very well-characterised star that we were able to reproduce pulsations scatter to better than 10 cm/s. We cannot know if the mismatch between simulations and observed values for HD 102365 and HD 20794 comes from a higher activity impact on pulsations or from poorly determined parameters; either can be responsible for the mismatch, without one having to invoke additional sources of intra-night scatter.

In the last section we saw that for the GP model on telluric-corrected HD 10700 RVs using TM/TMb/TMr the best estimate for the jitter of ESPR19 (measured over 3.5 yr) was of approximately 40 cm/s. This value is very similar the intra-night scatter generated by pulsations. The similarity might be coincidental, but in the worst-case scenario can be the consequence of the pulsations intra-night scatter not being averaged, and thus being present in the nightly averaged dataset. To test that the pulsations were effectively averaged out, we binned the data of HD 10700 considering only 3 points per night for the time series TM/TMb/TMr (with the actual points selected being randomly drawn from the available within each night). While the photon noise of the averaged points increases, it remains below 15 cm/s. The jitter on ESPR18 and ESPR19 increases to 85 and 45 cm/s, approx, pointing to an additional source of noise of 30–40 cm/s, with the other parameters remaining very similar. These tests reinforce the idea that the stellar pulsations signal were significantly averaged out.

8.2 Granulation’s signature and impact

Our understanding of the detailed physical mechanisms creating granulation RV signatures remains poor and has historically been influenced by the phenomenological characterisation of granulation coming from photometry. In the direct comparison in Fig. 4 of García & Ballot (2019) the photometry shows a very clear shape and knees that allow for the fitting of Harvey functions for granulation, while for RV the location, amplitude and presence of the different components is still not firmly established.

A Harvey function with α = 2, as commonly used for parametrising granulation, is a Cauchy distribution. This is a well-known pathological distribution, with both expected value and variance undefined; if we draw samples from this distribution, the sample average does not converge. This provides a simple mathematical explanation as to why the many works that simulate granulation from such prescription do not measure an increase in planetary detection capability when stacking up a very large number of observations (e.g. Meunier et al. 2023).

The expected granulation signature is clearly not detected in the intra-night scatter of HD 10700, being still compatible with those of HD 102365 and HD 304636. On the other hand, the predicted impact of super-granulation on average RVs, of 65–75 cm/s rms, is above the long-term stability of 40 cm/s measured on HD 10700. Our jitter measurement imposes then stringent constraints on granulation signals modelling and parametrisations.

8.3 A different signature: low-amplitude white-light flares

Main sequence stars with spectral types A to M exhibit flares, associated with events of magnetic reconnection in the solar corona. It has long been know that flares are detected predominantly on M dwarfs (e.g. Yang & Liu 2019). However, a comparison of flare size distributions on different hosts show these are mutually compatible across spectral types and with the distribution measured on the Sun (Aschwanden & Güdel 2021). This points towards detectability as the main source for the difference. The study of Pietras et al. (2022) on TESS data of 300 000 stars and using a photometric cadence of 2 min detected that 7.7% of stars show flares. The distribution of the duration of the flares peaks at 50 min and has a long queue up to several hours. The amplitude of the vast majority of the detected flares is of lower than 10% of the continuum flux, with the number of flares as a function of energy following a power-law.

In their seminal work on the physics of flares, Lu & Hamilton (1991) argue that the observed power-law dependence stems from the solar coronal magnetic field being in a self-organised critical state. In this scenario, solar flares are avalanches of many small reconnection events, where the size of a given flare is determined by the number of elementary reconnection events, and explaining why the phenomenon extends over more than five orders of magnitude in peak flux.

White-light flares are expected to extend to much lower amplitude than currently detected, and at a much higher frequency. These low-amplitude flares are expected to have an RV signature, that unfortunately has not been characterised yet. The recent work of Pietrow et al. (2024) compares flaring activity with contemporaneous HARPS-N and concludes that flares, if present, should have an RV signature below 50 cm/s. The timescale of 50 min measured on TESS photometry matches the knee on the power spectrum of solar RVs as measured by Al Moulla et al. (2023), and the power spectrum of these phenomena have been characterised by a Harvey function with α = 2 (see Hudson 1991). One can speculate that small-amplitude white light flares contribute to the RV variability measured by ESPRESSO at 40 cm/s. However, without a detailed study or benchmark on the impact of small-amplitude white-light flares it remains impossible to quantify the contribution.

8.4 Stellar signals or planets

The RV analysis of periodic signals only showed a significant signal at around 20 d for HD 10700, and under specific noise assumptions (without other sources of variance than photon noise). Being located close to an harmonic of the expected rotational period, at 42 d, the physical origin of this signal remains hard to assess. The best-match planet parameters posteriors show an over density of samples at this period. However, the RV amplitude of the signal varies from 15 to 60 cm/s, and is at the edge of our compatibility limits. The posterior of the GP η₃ parameter shows also a local maximum at 20 d, which could be a spurious fit to a planetary signal or simply the harmonic of the star rotational period at 42 d. The detection in RV of the first harmonics of the rotation period is well documented in the literature (e.g., Dumusque et al. 2017; Suárez Mascareño et al. 2017; Georgieva et al. 2023) and is even expected under some stellar configurations (Boisse et al. 2011). The recent work of Korolik et al. (2023) shows that HD 10700 is observed practically pole-on, and there are examples on the literature of pole-on stars like Vega that show a first-harmonic RV signature of the rotational period (Hurt et al. 2021). However, without a detailed characterisation of the star’s active regions, the modelling of the activity signal and confirmation of this hypothesis remains out of reach. The GP delivers a very low jitter value of ~40 cm/s to ESPR19 over 3.5 yr, to which contribute uncharacterised RV effects and undetected signals.

It is possible to make a comparison of these results with those obtained by Feng et al. (2017b) using the compatibility limits shown in Figure 9. The planets announced at 20.0 d (1.75 M_⊕) and 49.4 d (1.83 M_⊕) are very close to the detection limit. The 3.93 M_⊕ planet at 162 d should have been detected with our analysis, while the announced planet at 636 d is below our detection limits. In a nutshell, our data and analysis remains compatible with Feng et al. (2017b) on the potential presence of a signal at 20 d. However, we do not find significant evidence to call it a detection, and cast some doubts on its planetary origin (even if confirmed). Our results do not support the presence of a 3.93 M_⊕ planet at 162 d as previously announced. When interpreting the discrepancy in detected planets and detectability estimates between HARPS and ESPRESSO one should consider that HARPS can reach an intrinsic precision of down to 60 cm/s (see Lovis et al. 2006, for one of the first results). As such the detection of Earth-mass planets up to hundred days remains extremely difficult when done with HARPS, even with extensive datasets and the most advanced analysis.

The detailed study of HD 102365 by Tinney et al. (2011) was motivated by the large measured RV scatter (~3 m/s) when compared to the expected RV scatter on the star. The work of Suárez Mascareño et al. (2017) yields an an activity-induced RV semiamplitude of only 41 cm/s for the mean value of $\log \left(R_{HK}^{\mathrm{\prime }}\right)$, giving the star a privileged status for RV studies.

However, our ESPRESSO data shows clear evidence for activity-induced RV. First of all, the RV scatter doubled from ESPR18 to ESPR19. The LS of BIS, $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ and Ha₀₆ show signals almost reaching an FAP of 1% ate 123 d. Very importantly, the correlation between nightly averaged indicators and RV shows that half of the RV variance can be attributed to the correlation between an indicator and RV.

Finally, our kima analysis shows no evidence for the presence of a planet of 16 M_⊕ on a 122 days period orbit. The compatibility limits, calculated under the assumption of rotation-period informed GP model, show no candidate signals at this period and amplitude compatible with the ESPRESSO RV data. Therefore, we consider the 122–123 d signal as having a stellar origin.

The M star HD 304636 does not show a clear rotational period modulation in photometry, or on any spectral activity indicator. The kima analysis lead to a non-Gaussian jitter distribution for ESPR18 and poor detection limits around the expected rotational period of 30–40 d. A natural question is if the ESPRESSO data can provide further information on the TESS candidate found at 7.57 d. We used Forecaster (Chen & Kipping 2017) that for the radius of ~0.87 R_⊕ delivers a planetary mass of ~0.58 M_⊕. This value is well below the 1.5 M_⊕ (for 5 days) and 1.7 M_⊕ (for 10 days) compatibility limits found in this work, showing that the detection of such a planet is beyond our ability.

The kima model of the three stars delivers a compatible RV slope; this common behavior suggests a long-term variation of instrumental RV. The best characterisation is provided by HD 10700, delivering an RV slope of around 30 cm/s per year with an uncertainty of 10 cm/s.

It should be stated that the different RV calculations CCF, TM, TMb and TMr have different sensitivities to telluric contamination. This comes first as a consequence of the different method’s selection or exclusion of contaminated wavelength regions, but also due to the template matching calculating a unique template for a dataset. We showed that using telluric correction provides appreciable gains for HD 10700, but a finer analysis on a wide range of targets can be done.

As a concluding remark, we note that while HD 10700 and HD 102365 have very similar activity values as measured by the most common indicators, the stars show a very different RV scatter and a different correlation level between indicators and RV. Moreover, they are modelled with kima with very different success, both in terms of the posteriors reached and the final jitter value reached. This exemplifies how the a priori measurement of activity indicators is still insufficient when selecting the most suitable stars for the detection of Earth-mass planets. This conclusion expands the results of the systematic study of Gomes da Silva et al. (2021), where it was shown that low-activity stars are not necessarily low-variability stars.

9 Conclusions

We present a study of main sequence stars using ESPRESSO GTO data. We started with a stellar characterisation using complementary methods to draw a complete picture of the stars. The RVs were derived from the ESPRESSO spectra with the CCF and TM methods and ancillary photospheric and chromospheric activity indicators.

We were able to study the intra-night RV behaviour on the three G–K stars with multiple consecutive exposures per night. We compared the measured RV signals with the expected scatter from pulsations and granulation. We reached a very good match between the predicted RV variability from oscillations and measured RV down to a level of better or equal to 10 cm/s, with an additional scatter that might be attributed to poor stellar parameter determination or the impact of activity on oscillations.

A kima analysis was used to evaluate the number of planets supported by the data, under the assumption that a quasi-periodic GP can model the activity signal. A rotationally modulated signal close to the expected period was modelled for the three stars, with the analysis for HD 10700 being particularly successful in reproducing the observed data. A planetary signal at a period of 20 d could be modelled along with it, but the evidence for its presence is not significant. Moreover, its physical origin remain uncertain due to the similarity between this period and the first harmonic of rotation. ESPRESSO data on their own do not provide conclusive evidence for planets around HD 10700, HD 102365, or HD 304636. This is at odds with the detection by Tinney et al. (2011) of a 123 d period planet with 16 M_⊕ around HD 102365, which would clearly have been revealed by our data. The latest analysis of Feng et al. (2017b) on HD 10700 HARPS archive data also predicts a planet with 3.93 M_⊕ planet at 162 d that while within our detection capabilities is not detected, with our analysis casting a different interpretation on the 20 d signal. Our RV analysis allowed us to detect planets around HD 10700 down to 1.7 M_⊕ for periods up to 100 d and of 2–5 M_⊕ within the star’s habitable zone. As data analysis progresses, these values are expected to decrease.

In this work, ESPRESSO demonstrates an on-sky RV precision of better 10 cm/s on short time-scales (<1 h) and of around 40 cm/s on a timescale of 3.5 yr on HD 10700 telluric-corrected spectra. By further further splitting the datasets, we estimated that the RV stability over a timescale of one year can be as low as 25–30 cm/s, variable from year to year.

The convective phenomena of granulation were modelled with literature prescriptions that unfortunately do not reproduce the level of scatter we observed in this work. In particular, the parametrisations of super-granulation overestimate the scatter measured in our RVs for HD 10700, which provides strong constraints for detailed modelling in future studies. The analysis presented here illustrates the significant differences among stellar RV signatures when analysed at sub-m/s level. It underlines the fact that no current activity indicator allows us to select RV-stable stars in advance for searches at better than 1 m/s.

Data availability

The RV data for the targets discussed in this work are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/700/A174.

The ESPRESSO data, reduced with the latest version of the pipeline, can also be accessed through the DACE platform via https://doi.org/10.82180/dace-3js7g1l9.

In parallel, the ESPRESSO Legacy data release can be found at https://www.doi.org/10.82180/dace-3g5jf9h8.

Acknowledgements

In memory of Denis Mégevand, project manager of ESPRESSO and dear friend, who passed away in 2022. We thank the anonymous referee for a review and comments that very much improved the readability of the manuscript. PF and JPF thank Patrick Eggenberger for very fruitful discussions on the impact on RV of p-mode oscillations, Ângela Santos for her expertise on solar physics and stellar flaring, Khaled Al-Moulla for discussing the temporal analysis of HARPS and HARPS-N data, and Alex Pietrow for a discussion on the state-of-the art on flaring observations. 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 grant 51NF40_205606. The authors acknowledge the financial support of the SNSF. FPE and CLO would like to acknowledge the Swiss National Science Foundation (SNSF) for supporting research with ESPRESSO through the SNSF grants nr. 140649, 152721, 166227, 184618 and 215190. The ESPRESSO Instrument Project was partially funded through SNSF’s FLARE Programme for large infrastructures. Co-funded by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This work was supported by FCT – Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UIDB/04434/2020; UIDP/04434/2020. A.C.-G., O.B.-R. and J.L.-B. are funded by the Spanish Ministry of Science through MCIN/AEI/10.13039/501100011033 grants PID2019-107061GB-C61. HMT, PID2023-150468NB-I00 and CNS2023-144309. 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. HMT also acknowledges financial support from the Agencia Estatal de Investigación (AEI/10.13039/501100011033) of the Ministerio de Ciencia e Innovación through projects PID2022-137241NBC41 y PID2022-137241NB-C44. J.I.G.H., A.S.M., and R.R. acknowledge financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) projects PID2020-117493GB-I00 and PID2023-149982NB-I00, and from the Government of the Canary Islands project ProID2020010129. X.D. acknowledges the support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement SCORE No 851555) and from the Swiss National Science Foundation under the grant SPECTRE (No 200021_215200). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work made use of astropy, a community-developed core Python package and an ecosystem of tools and resources for astronomy, described in astropy collaboration (2013), astropy collaboration (2018), and astropy collaboration (2022). Corner plots were created using the package corner.py, described in Foreman-Mackey (2016).

Appendix A ARIADNE priors and input photometry

For the effective temperature, log(g), and [Fe/H] ARIADNE uses as priors the observed parameter distributions from the fifth RAVE survey data release (Kunder et al. 2017). A Normal distribution centred around the Bailer-Jones et al. (2018) value and with a standard deviation of 5 times the associated uncertainty is used as prior for the distance. As we wanted to retain uninformative priors we kept these default choices; for a full description we refer to Vines & Jenkins (2022). Since the distance to these stars is at maximum 10 pc, we fixed the extinction coefficient A_v at 0, in what is our only informed choice.

The photometry and associated uncertainties retrieved for the four stars is presented on Table A.1. Since 2MASS J, H, and K_s photometry is saturated for magnitudes below 4.5, 4, and 3.5, respectively, we exclude from the analysis of HD 10700, HD 20794, and HD 102365 the 2MASS data ²³.

Gaia DR3 photometric measurements are also affected for saturation for bright targets. Nonetheless, the consortium provides a correction formula readily applicable to our case, described in Riello et al. (2021), Appendix C1. Following this correction all Gaia magnitudes could be used.

Appendix B Rotational period derived assuming different spectral types

For the sake of completeness and comparison, we computed the rotational period using the calibrations of Suárez Mascareño et al. (2016) for each star’s estimated spectral type, closest class, plus GKM’s as a whole. For $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ we use the value derived by Hojjatpanah et al. (2019). The results are shown in Table B.1.

Appendix C S-BART application to ESPRESSO data for high-precision radial velocities

The specifics of our application of S-BART to ESPRESSO data can be summed up as follows:

• the method assumes a common RV value for the individual orders and slices of the instrument. In this context, the two ESPRESSO slices of each interference order are treated as two independent measurements, ignoring the correlations between the two orders flux values and their uncertainties;

• To decrease the impact of micro-telluric features on the stellar template, only observations taken at an airmass smaller than 1.6 are considered for the template construction.

• S-BART masks wavelength regions to prevent spectral contamination from Earth’s atmosphere on the spectra, that leads to parasite effects on RV. A reference synthetic transmittance profile is constructed with TELFIT (Gullikson et al. 2014) at a resolution of 140000 for the wavelength range of ESPRESSO; the airmass, stellar RV and precipitable water vapour of the observation with the highest relative humidity are assumed as worst-case scenario for each star. Wavelength bins with a transmission lower than 99% relative to the continuum and within the maximum yearly BERV variation (~30 km/s) are masked out. This mask is applied to both the construction of the stellar template and the calculation of the individual RVs. This criteria leads to discarding 36 out of the 170 slices of ESPRESSO.

• two independent templates are generated for the two independent datasets ESPR18 and ESPR19;

• When discarding wavelength domains, and to ensure that the RV measurements are informed by the same spectral regions, we consider only orders accepted for both ESPR18 and ESPR19.

Appendix D P-mode simulations

D.1 Reference values used in the scaling relationship

The reference values for asteroseismic and stellar parameters are presented in Table D.1, with the corresponding bibliography.

D.2 Simulation diagrams and examples

A representation of how the different frequencies and amplitudes generate different pulsation components is shown in Fig. D.1.

Fig. D.1 Representation on the frequencies and amplitudes drawn from the modes envelope (left), and how these plus random phases create the signals in the temporal domain (right). From top to bottom we have the main frequency νmax, plus first and second order terms. In the temporal domain, open shutter is represented in blue.

An example of a simulated night sequence for HD 102365 is represented in Fig D.2, depicting how the open shutter time covers the final pulsation signal. Once a large number of realisations of this simulation is done, one can evaluate the distribution of nightly averaged RVs and RV scatter (rms), as seen in Fig. D.3.

Fig. D.2 Example of a simulated pulsation signal for HD 102365. The open shutter time intervals are represented in blue.

Fig. D.3 Distributions of the nightly averaged RV (top) and intra-night scatter (bottom) introduced by pulsations, calculated for 50 000 runs for the star HD 102365.

Appendix E Granulation measurements and simulations

Dumusque et al. (2011) used HARPS asteroseismology data on the stars β Hyi, α Cen A, μ Ara, and τ Cet and fitted the PSD with Harvey and Lorentzian function, applying to RV the same principles studied in detail in photometry. In a more recent work Al Moulla et al. (2023) used HARPS and HARPS-N solar telescopes RV data and fitted two Harvey functions to represent granulation on the RV’s power density spectrum. While Dumusque et al. (2011) fitted three granulation signals, Al Moulla et al. (2023) fitted two, and the power-law slope α was fixed at 2.0 to avoid degeneracy and because this value reproduced well the observational data. The results are summarised in Table E.1.

From the proposed scaling relationships one gets that the amplitude and granulation timescales are proportional to Δν; yet, for α Cen A and mu Ara, two stars with very similar large frequency separations, the parametrisations of granulation obtained by Dumusque et al. (2011) is very different. This can be a consequence of poor time sampling of the phenomenon or the poor granulation modelling, but one should also consider that the scaling model could be inaccurate, as suggested already in Mathur et al. (2011).

Appendix F FAP from ℓ₁ analysis

The FAP as estimated from ℓ₁ analysis for the three stars is presented in Tables F.1, F.2, and F.3.

Appendix G kima analysis

G.1 Description of prior distributions

We use non-informative priors as much as possible, including only information available from simple summary statistics and the previously determined rotation periods listed in Table 5.

The modified log-uniform distribution is a log-uniform distribution with parameters knee a and upper limit b, introduced by Gregory (2005). Using the modified version allows us to extend the prior to zero, assigning non-zero probability to values below the knee value. For the planet semi-major amplitude we use a modified log-uniform prior with a knee at 10 cm/s and an upper limit of 10 m/s.

For the period we use a log-uniform prior from 1 to 1000 d. The eccentricity is represented with a Kumaraswamy (1980) distribution which is an approximation to the β distribution proposed in Kipping (2013). For the mean anomaly and argument of periastron we use uniform distributions between 0 and 2π.

For the stars systemic RV and the instrument-related parameters we derive the priors from the datasets: the systemic RV is taken from a uniform distribution between the RV extreme values, and the jitter from a uniform distribution between 5 cm/s to 5 m/s. The offset is drawn from a uniform distribution between minus and plus the RV span value. The linear and quadratic slopes are zero-centred normal distributions with variance given by the data dispersion.

As the GP we use kima parametrisation of the quasi-periodic kernel $$k_{QP}(t,t^{\mathrm{\prime }})={\eta }_1^2\exp [-\frac{2}{{\eta }_4^2}{\sin }^2\left(\frac{\pi (t-t^{\mathrm{\prime }})}{{\eta }_3}\right)-\frac{{(t-t^{\mathrm{\prime }})}^2}{2{\eta }_2^2}],$$(G.1)

where η₁ is the amplitude of the kernel, η₂ is the evolution timescale, η₃ is the rotational period and η₄ is the harmonic complexity. η₂ represents the timescale of coherence of the active regions. The adimensional parameter η₄ governs the short-scale variation of the GP relative to the rotation timescale η₃; a smaller value of η₄ indicates a higher number of short-scale structures within a single stellar rotation period. For further information of the quasi-periodic kernel GP and on GPs in general, we refer to Rasmussen & Williams (2006) and Roberts et al. (2012).

It is worth mentioning that this GP model has several potentially degenerate variables. Equation G.1 is basically the product of two exponential functions: the term with η₃ and η₄, designed to represent the periodic part with higher harmonic complexity, and the term with η₂, intended to represent the temporal evolution. However, for a small temporal difference relative to the putative period (t − t′)/η₃, when 2πη₂ = η₃η₄ there is a potential issue²⁴. Since η₂ ≥ η₃, for η₄ < 1 there is no confusion; however, when η₄ > 1 there is a degeneracy between the terms that intend to represent temporal evolution and periodic behaviour. Moreover, for very low or high values of harmonic complexity η₄, the final function is not quasi-periodic. For η₄ ≪ 1 the function will contain a very large number of peak-to-valley variations per timescale η₃, while for η₄ ≫ 1 the function will be constant over timespans much longer than η₃, and including the GP corresponds to simply adding a second component of white noise. The term of harmonic complexity η₄ has its origin from the temporal evolution of spot complexes as seen in photometry, but its usage in RV is far from trivial. The exact value is hard to recover even when injected on simulated RV datasets, and its interpretation at the physical level is not trivial; for a discussion on the subject, we refer to Nicholson & Aigrain (2022).

For η₁ we consider a uniform distribution between 5 cm/s and 5 m/s. For the rotation timescale η₃ we use a Normal distribution centred at P_rot and with standard deviation σ as the period’s estimated uncertainties, with the two variables as presented in Table B.1. For the evolution timescale parameter η₂ we consider a Uniform distribution between 1 and 5× the rotational period. There is little information about this parameter on stars other than the Sun, but on our star the values of η₂ are assumed to be a few times the rotation period. For η₄ we choose a uniform between 0.2 and 5.

All priors are listed in Table G.1.

G.2 Posterior samples for HD 10700 TM run

G.3 Telluric corrected HD 10700 time series

Fig. G.1 Posterior distribution for the orbital period (top) and joint posteriors semi-amplitude-period (middle) and eccentricity period (bottom) for the TM RV calculation method of HD 10700.

Fig. G.2 Lomb-Scargle periodogram for the quantities ΔRV, ΔFWHM and Δcontrast calculated by subtracting telluric corrected time series from uncorrected ones for HD 10700. The vertical dashed lines represent one year and fractions of it, from 1/2 to 1/10.

G.4 Best-fit parameters

G.5 Corner plots

Fig. G.3 Corner plot for the no-planet GP analysis done with kima for HD 10700 TM RVs.

Fig. G.4 Corner plot for the no-planet GP analysis done with kima for HD 102365 TM RVs.

Fig. G.5 Corner plot for the no-planet GP analysis done with kima for HD 304636 TM RVs.

References

All Tables

All Figures

	Fig. 1 Distribution of the photon noise uncertainty values (top) and observed nightly scatter (bottom) for the stars HD 10700 (left), HD 20794 (centre), and HD 102365 (right). The solid-line histograms represent the data and the dashed lines the cumulative log-normal distribution fits.
In the text

	Fig. 2 RV time series using the different RV calculation methods CCF (top) and difference relative to it of TM (center), TMb, and TMr (bottom) for the stars HD 10700 (left), HD 102365 (center), and HD 304636 (right). ESPR18 data are represented by open circles while ESPR19 data are represented by filled circles.
In the text

	Fig. 3 Time series for the line profile and chromospheric indicators for the stars HD 10700 (left), HD 102365 (centre), and HD 304636 (right). ESPR18 data are represented by open circles, while ESPR19 data are represented by filled circles.
In the text

	Fig. 4 Lomb-Scargle periodogram for the different methods of RV calculation and associated window function for the stars HD 10700 (left), HD 102365 (center), and HD 304636 (right). The FAP levels of 1, 0.1, and 0.01% are shown as horizontal grey lines.
In the text

	Fig. 5 ℓ1 periodogram for the stars HD 10700 (left), HD 102365 (centre), and HD 304636 (right). The ten strongest peaks periods are identified in the image; the associated FAP is listed in Table F.
In the text

	Fig. 6 Lomb-Scargle periodogram for the activity indicators for the stars HD 10700 (left), HD 102365 (centre), and HD 304636 (right). On the upper panels we have the photospheric indicators FWHM, line contrast, and BIS, and on the bottom on we have $\log (R_{HK}^{\mathrm{\prime }})$ and the two Ha parametrisations. The FAP levels of 1, 0.1, and 0.01% are shown as horizontal grey lines.
In the text

	Fig. 7 Radar chart of the correlation between nightly averaged RV and each of the indicators, for the four RV calculation methods and three stars. The different colours represent the different RV datasets, with filled or open dots representing positive or negative signed Pearson correlation coefficients.
In the text

	Fig. 8 Normalised posterior distribution of the two instrument jitter and four GP parameters as estimated by kima, for the three stars studied in the section and the TM RV calculation method. The distributions of the posteriors for GP-only (dashed line) and up to one planet runs (solid line) are both shown.
In the text

	Fig. 9 Compatibility limits for the three stars studied, HD 10700 (top), HD 102365 (middle), and HD 304636 (bottom). The maximum semi-major RV amplitude of the samples is plotted for each of the RV calculation methods model, and telluric correction set; HZ limits as presented in Table 5 and orbital properties of announced planed described in Sect. 2 are shown.
In the text

	Fig. 10 Boxplots for the posteriors of the jitter (top) and offsets between datasets (bottom) for the different RV calculation methods with five datasets of telluric-corrected RVs of HD 10700.
In the text

	Fig. D.1 Representation on the frequencies and amplitudes drawn from the modes envelope (left), and how these plus random phases create the signals in the temporal domain (right). From top to bottom we have the main frequency νmax, plus first and second order terms. In the temporal domain, open shutter is represented in blue.
In the text

	Fig. D.2 Example of a simulated pulsation signal for HD 102365. The open shutter time intervals are represented in blue.
In the text

	Fig. D.3 Distributions of the nightly averaged RV (top) and intra-night scatter (bottom) introduced by pulsations, calculated for 50 000 runs for the star HD 102365.
In the text

	Fig. G.1 Posterior distribution for the orbital period (top) and joint posteriors semi-amplitude-period (middle) and eccentricity period (bottom) for the TM RV calculation method of HD 10700.
In the text

	Fig. G.2 Lomb-Scargle periodogram for the quantities ΔRV, ΔFWHM and Δcontrast calculated by subtracting telluric corrected time series from uncorrected ones for HD 10700. The vertical dashed lines represent one year and fractions of it, from 1/2 to 1/10.
In the text

	Fig. G.3 Corner plot for the no-planet GP analysis done with kima for HD 10700 TM RVs.
In the text

	Fig. G.4 Corner plot for the no-planet GP analysis done with kima for HD 102365 TM RVs.
In the text

	Fig. G.5 Corner plot for the no-planet GP analysis done with kima for HD 304636 TM RVs.
In the text