© The Authors 2025

1 Introduction

The field of exoplanet research has seen a massive growth over the past three decades, with 5849¹ detections of exoplanetary signals. Despite this impressive growth, the discovery and characterisation of ‘Earth analogs’ still remains a major challenge, driving improvements in both instrumentation and data analysis techniques. A major workhorse of the field is the radial velocity (RV) method, based on the measurement of variations of a star’s velocity along the line of sight. This is done via the measurement of displacements in the location of the spectral lines in comparison to their expected position.

Often, RV extraction is done through the cross-correlation function (CCF) method, based on the cross-correlation between stellar spectra and a (weighted) binary mask (Baranne et al. 1996; Pepe et al. 2002a). However, this is not always the optimal approach. In cooler stars (e.g. M dwarfs), the template-matching (TM) algorithm typically yields better results, leading to smaller RV scatters and improved RV precision, as seen in different implementations, such as HARPS-TERRA (Anglada-Escudé & Butler 2012), NAIRA (Astudillo-Defru 2015), and SERVAL (Zechmeister et al. 2018). Template matching is a data-driven method, based on the construction of a high-S/N stellar model from available observations of a given star, which is then aligned with the individual observations. Usually, TM methods share a common set of characteristics: (i) chromatically independent RV shifts, where each spectral order is assumed to be an independent measurement of the same underlying phenomenon; (ii) the stellar template (model) is capable of accurately describing every observation of a given star. s-BART² (Silva et al. 2022) tackled the assumption of chromaticity, as the RV signal of exoplanets is independent of the wavelength at which it is measured. This is done by imposing a wavelength-independent RV shift - common to all spectral orders - to describe the differences between the stellar model and observations. Since its inception, s-BART has used for the discovery and characterisation of several systems across different state-of-the-art spectrographs (e.g. Faria et al. 2022; Palethorpe et al. 2024; Balsalobre-Ruza et al. 2025). To refine template-matching avenues, we have seen a push towards a careful selection of the lines that are used for RV extraction. Line-by-line (LBL) methods (e.g. Dumusque 2018; Artigau et al. 2022) share a common working principle with TM methods: a high-S/N stellar model is constructed from the data and, afterwards, RVs are estimated through their analysis. This process is carried out for each individual spectral line and, finally, a selection is made in order to remove the ones that present a larger variability.

In exoplanetary RV blind-searches, or in the follow-up of transiting exoplanets, observations are typically carried out over multiple months. This allows us to properly sample the orbital period of the orbiting companions and to characterise any stellar signals. However, some science cases require a much shorter cadence of observations, such as the transmission or emission spectroscopy studies (e.g. Wyttenbach et al. 2015; Costa Silva et al. 2024), or asteroseismological analysis (e.g Bouchy et al. 2005; Campante et al. 2024). Similarly, the RV detection of planets with unusually ultra-short periods (e.g. P_orb < 12 h, Adams et al. 2021; Castro-Gonzâlez et al. 2025) can be achieved through very short observing baselines. These datasets are typically collected within a single night or over a few consecutive nights.

A priori, one would expect the RV extraction through a template-based algorithm (TM or LBL) to perform similarly for different observing strategies. However, in this paper, we show that this is not necessarily the case. In particular, we present a systematic bias in the RV extraction through template-based algorithms when the stellar model is constructed from observations collected within a short time span (e.g. within a night). This bias appears in two TM pipelines, one LBL pipeline, is present in the data of multiple state-of-the-art spectrographs, and it has different amplitudes for different stars. The bulk of our analysis is done mainly through the publicly available s-BART pipeline, using data from ESPRESSO (Pepe et al. 2021) and HARPS (Mayor et al. 2003; Pepe et al. 2002b).

We start by providing a description of our datasets in Section 2, followed by an overview of the methods for RV extraction in Section 3. Then, in Section 4, we showcase the magnitude of this spurious signal in different stars and explore its wavelengthdependence. In Section 5 we focus our analysis on ESPRESSO observations and explore the effect that the selection of observations has on the retrieved RVs. In Section 6 we hypothesise on the root cause of this systematic bias. Lastly, we wrap-up in Section 7 with a discussion on the possible avenues to mitigate this bias, as well as the impact it introduces on different science cases.

2 Overview of the observations

We use observations collected with the ESPRESSO and HARPS spectrographs. The raw data was processed with the corresponding version of the official data reduction software (DRS) of each spectrograph³. It handles the extraction of spectral orders from the raw images, performs the necessary calibrations and corrections, and provides an RV extraction through the CCF method. The CCF RVs are obtained through a Gaussian fit to the outcome of cross-correlating the stellar spectra with a weighted mask that is selected from a pre-defined list. This selection is based on the spectral type, using the one that most resembles the star that is being processed. We used HARPS DRS version 3.2.5, whilst the bulk of ESPRESSO data was reduced with pipeline version 3.0.0⁴.

Our analysis will be focused on three different datasets: aster-oseismology observations, long-period coverage of τ Ceti, and a large sample of single-night observations. We provide a detailed description of each dataset in the following subsections.

2.1 Asteroseismology datasets

Our asteroseismology datasets comprise observations of three different stars:

1. HD 160691 - This G3 star has 2183 HARPS public spectra, available through the DACE interface, under ESO program 073.D-0578(A). This dataset was the basis for the discovery of the first Neptune-like exoplanet ever announced (Santos et al. 2014). For the purposes of our analysis, we are only interested in the asteroseismology time series (i.e. those that span a few hours). We found six consecutive nights (200406-05 to 2004-06-10, Bouchy et al. 2005) and selected the last two for our analysis.

2. HD 209100 ( Indi) - This K5 star was observed with ESPRESSO for six nights on December 2022 - ESO program 109.236P.001 - allowing for the detection of solar-like oscillations. These have a peak amplitude (radial modes) of 2.6 cm s⁻¹, making this the coolest seismic dwarf ever observed (Campante et al. 2024).

3. HD 40307 - This bright K-dwarf star was observed during a five-night campaign during ESPRESSO’s commissioning, ESO program 0102.D-0346(A), with a tentative detection of p-mode oscillations (Pepe et al. 2021). However, the dataset was revisited by Campante et al. (2024), which found no evidence for their presence.

In total, we used 475 HARPS observations of HD 160691, 2084 ESPRESSO observations of Indi and 1140 ESPRESSO observations ofHD 40307. The BERV coverage of the three stars can be found in Figure A.1 and Table A.1 presents an overview of the nights in which the stars were observed, alongside their median S/N at a wavelength of 600 nm.

2.2 HD 10700 (τ Ceti)

HD 10700, commonly referred to as τ Ceti, is a bright G8 star that was intensively observed within the Guaranteed Time of Observations (GTO) of ESPRESSO⁵. The ESPRESSO spectrograph observed τ Ceti for close to five years with a median number of 15 observations per night. In total, 2442 observations were collected and a detailed analysis of this system (using CCF and s-BART RVs) can be found in Figueira et al. (2025). Unlike other ESPRESSO observations, the τ Ceti dataset was reduced with ESPRESSO DRS version 3.3.0 for the analysis of Figueira et al. (2025). Figure A.2 presents the BERV coverage of the ESPRESSO observations, highlighting in blue the nights we found a decrease in BERV (in comparison with the previous night). In this paper, we also looked at individual nights of this star, whose BERV distribution is shown in Figure A.3 and some of their properties are summarised in Table A.2.

2.3 Large sample of single-night observations

The last dataset in our analysis is made up of a large sample of ESPRESSO targets. The majority of them were selected from within ESPRESSO GTO’s working group 2 (atmosphere characterisation), as it collected a large number of multi-hour time series over multiple nights⁶. All in-transit observations were discarded from the analysis to avoid contamination from the spectral imprint of a transiting exoplanet. This is done through a manual identification of the observations taken outside of the planetary transit (typically known as out-of-transit observations). To those observations, we added the ESPRESSO observations of HD 40307 and Indi (described in Section 2.1).

In total, we collected 4124 ESPRESSO observations of 19 different stars that were observed over 38 nights. Figure A.4 presents an histogram of the BERV distribution within this sample, showing a good coverage of the ±30 km s⁻¹ maximum barycentric Earth radial velocity (BERV) variation over the year.

3 RV extraction methods

In this section, we present an overview of the RV-extraction process of the different algorithms that will be applied to the stellar spectra. In particular, we focus on the s-BART algorithm, using it for the bulk of our analysis.

3.1 The s-BART algorithm for template-matching RVs

The semi-Bayesian Approach for RVs with Template matching (s-BART) algorithm (Silva et al. 2022) is a TM algorithm for RV extraction, built around a common assumption: the RV shift that an orbiting companion introduces is independent of the wavelength at which it is measured. This is enforced by assuming a common RV shift to describe all differences between the stellar model (the template) and the data with which the template is compared to (the individual spectra). Furthermore, by implementing this algorithm within a semi-Bayesian context, we have access to a straightforward and consistent method to characterise the posterior distribution of the RVs. Below, we present a brief overview of the RV extraction process, referring to Silva et al. (2022) for a full description.

• We rejected wavelength regions based on the quality control (QC) flags provided by the instrument’s official pipeline. We also rejected wavelength regions that are known to be sensitive to activity and spectral orders that do not meet the internal quality control flags of s-BART.

• We constructed a transmittance profile of Earth’s atmosphere using Telfit (Gullikson et al. 2014), assuming the worst atmospheric conditions of our dataset, defined as the observation with the highest relative humidity. Then, we constructed a binary mask to reject any telluric line that has a depth greater than 1%. This is a very conservative threshold, leading to a rejection of several wavelength regions. However, it also minimises the number of micro-telluric features that are left in the stellar spectra. This binary mask was also expanded to account for the maximum yearly BERV variation of the star.

• Using all observations, we constructed an high signal-tonoise-ratio (S/N) model by stacking the individual spectra and propagating the flux uncertainty. In this process, we placed all observations in a common-rest frame using a previous RV estimate; computed through the CCF approach or an initial run of the TM algorithm, for instance. Then, we interpolated the spectra to a common wavelength grid. We also enforced the condition that any wavelength gap in any of the observations would be translated into a gap in the stellar template, ensuring that every pixel in the stellar template has a corresponding pixel in each observation.

• We applied the s-BART algorithm to align the stellar template with the spectra. The RV posterior distribution was characterised via a Laplace approximation, which assumes the posterior distribution is Gaussian in shape. The RV measurement can be obtained through the maximum of the posterior distribution, along with the associated uncertainty via its standard deviation. During the RV extraction process, we accounted for changes in the stellar continuum through the marginalisation of a first-degree polynomial to describe differences between the stellar spectra and the template.

It is also important to note that in our RV extraction, we ensured we were using the same orders for all observations. We note that this also translates into a loss of RV information, as the rejection of one order in a single observation will imply the loss of that same order in all other observations. However, this ensures that the final RV measurement of all observations will be informed by the same wavelength regions.

3.2 The SERVAL pipeline for template-matching RVs

The SpEctrum Radial Velocity AnaLyser (SERVAL) templatematching pipeline is a widely used implementation of the template matching method. It constructs the stellar template through a B-spline regression of the continuum-normalised individual observations. Then, order-wise RVs are computed through a χ² minimisation between the Doppler-shifted template and the observations. Lastly, the order-by-order radial velocities are combined with a weighted mean to produce the final RV estimate for the observation.

3.3 The ARVE algorithm for line-by-line radial velocities

Analyzing Radial Velocity Elements⁷ (ARVE, Al Moulla 2025) is a Python code that has the capability (among other functionalities) to extract line-by-line (LBL) RVs. Unlike other LBL algorithms and to speed up computations, ARVE utilises a grid of synthetic spectra, from which the properties of expected spectral lines within a provided wavelength range are retrieved. The LBL RVs are thereafter calculated through a per-lineχ²-minimisation between the observed spectrum and a template constructed from all or a subset of the observations. Specifically, as in Bouchy et al. (2001), the observed spectrum, S_obs, is described as the scaled first-order Taylor expansion of the template spectrum, S, $$S_{\mathrm{o}\mathrm{b}\mathrm{s}}=A(S+\frac{\mathrm{d}S}{\mathrm{d}\lambda }\mathrm{\Delta }\lambda ),$$(1)

where the only free variables are the scaling, A, to account for integrated-flux differences between the observation and template, and the RV shift, v, which relates to the wavelength differential, Δλ=λν/c.

Fig. 1 Mean-subtracted RV time series for one night (December 24, 2018) of ESPRESSO observations of HD 40307. The CCF RVs (blue triangles) are compared against those from s-BART (black), SERVAL (red), and ARVE (orange).

4 Systematic RV bias for short time baselines

During the development and subsequent application of s-BART to multiple targets, we found a transversal issue to templatebased techniques for RV extraction: if all observations are acquired within a short period (i.e. within a few days), then a quasi-linear systematic trend in RV appears over each individual night. This issue is present in the observations of multiple spectrographs and multiple stars, with varying amplitudes (Section 4.1). This RV bias was also reproduced using a different TM pipeline: the publicly available⁸ version of SERVAL. In Figure 1, we present a comparison of the RVs extracted for one night of HD 40307, showing that all template-based algorithms present a large peak-to-peak RV variation (~40 m s⁻¹), whilst the CCF RVs do not.

This rules out a possible pipeline error and instead reveals a systematic bias in template-matching RVs. This is, however, unexpected as s-BART (similarly to other TM pipelines) provides improvements on the cooler stars when the stellar template is constructed with data collected during a long time series (e.g. Silva et al. 2022). The fact that a LBL extraction of RVs also reveals a similar bias leads to two possible explanations: (1) the issue is connected with the spectral lines, as LBL approaches do not use the continuum; or (2) the root source of the bias is linked to the stellar template, as both TM and LBL methods utilise a data-driven, high-S/N, stellar template. As shown in Appendix B, we failed to find a clear correlation between the RV residuals and the line-shape metrics that would be derived from the CCF. In the next sections, we explore the behaviors of the observed trend.

4.1 Variation across targets

In this section, we set out to characterise the behavior of this RV trend based on a larger sample of observations (described in Section 2.3). We do not include the HARPS observations, preferring to use data from the same instrument, to avoid differences that could be caused by the instrument.

We subdivided the data of each star, so that each of the 38 nights of observations is treated fully independent from one another. Then, for each night, we applied the following prescription:

1. For each night, we constructed a telluric binary mask spanning the maximum yearly BERV variation for the target. This leads to a larger rejection of wavelength regions, but we preferred to keep a consistent telluric masking process across this paper.

2. We applied the telluric mask to remove the telluric features of each observation of the night, which were then used to construct the stellar model and extract s-BART RVs.

3. We computed the residuals between the s-BART and CCF RVs and fit a first-degree polynomial to them. We note that the linear model fitting is an simplification of the effect that we are measuring, but as a first approximation, it allows us explore the magnitude of the RV bias.

Figure 2 presents the distribution of the slopes (of the residuals) and their correlation with different metrics, with Appendix C presenting their values. The RV bias that we retrieve always shows a negative slope: independently of the star, we find a decrease of RVs over time. This result allowed us to discard the hypothesis that the TM algorithm is more sensitive to RV variations introduced by stellar activity, since there is no reason these variations would be in the same direction for all stars.

For the majority of the targets, we find that the absolute value of the residuals RV slope is smaller than 10 m s⁻¹ h⁻¹ and the values range from 0.3 m s⁻¹ h⁻¹to 53.6 m s⁻¹ h⁻¹, with Table C.1 presenting the individual measurements. It is puzzling to consider why such a large variation is not present on the CCF RVs and does not translate into a change in the line-shape metrics that the CCF provides.

The bottom panels of Figure 2 show a scatter plot between the residual slope and different metrics: the intra-night BERV slope; the effective temperature as taken from the SweetCat catalogue (Sousa et al. 2024); the median S/N of the spectra at 600 nm⁹; and the median difference between BERV and CCF RV. Here, we find two noteworthy results: i) the RV bias is not uniform across multiple nights for the same star. An example is HAT-P-26, which was observed in two different nights with residual slopes of −27 ms⁻¹ h⁻¹ and −6.9 ms⁻¹ h⁻¹; ii) the magnitude of the bias increases as the S/N decreases. We find that for S/N levels under 100 we have a large dispersion of the retrieved RV slopes, whilst for S/N greater than 100, the maximum residual slope found was of −1.02 m s⁻¹ h⁻¹;

4.2 Wavelength-dependence

One of the core assumptions behind s-BART is that a common RV value is enough to describe all differences between the high-S/N stellar model and the spectra. This assumption is in line with possible orbiting companions (the signals that we are interested in detecting and characterising), but not with stellar signals and/or instrumental systematics.

As a simple test for wavelength dependence of the bias, we computed two different RV time series: one that only uses the spectral orders that exist in the “blue” detector and another that only uses information from spectral orders that fall in the “red” detector. We note, however, that the wavelength coverage of the blue and red detectors is not the same for the two spectrographs.

For this analysis, we used the two nights of HARPS observations of HD 160691 and two of the ESPRESSO nights of observations of Indi. We preferred to use the data from ϵ Indi as it presents a larger flux-S/N than HD 40307. A subset of the nights was chosen to illustrate the behaviour, as the results were compatible across all nights.

In Figure 3, we present the residuals between CCF and s-BART RVs for both detectors. We find that in all targets the systematic effect has a larger amplitude in the red detector. The RV time series from the blue detector still presents a bias, albeit with a smaller amplitude. Unfortunately, our data are not sufficient for us to draw conclusions from the difference in amplitude between the two spectrographs, as we find large variations even within a single spectrograph (see Section 4.1).

It is also important to note that the RV time series that utilises the full spectral information will reflect the average effect on the two detectors. Although it is not represented here, we also confirm that this RV drift over time is also seen at the level of the individual spectral orders. From these results, we conclude that the bias is present at different wavelengths, showing a small wavelength dependence in amplitude.

Fig. 2 RV bias for 3785 ESPRESSO observations of 20 different stars. The top row presents an histogram of the slopes of a first degree polynomial that was adjusted to the residuals between s-BART and CCF RVs. The bottom two rows present a scatter plot with the median value of different metrics: intra-night slope of BERV variation, effective temperature, median S/N at 600 nm, and median difference between CCF RV and BERV. All observation nights of the same star are shown with a common color. For representation purposes we do not show the two outliers (with a slope greater than −20 ms−1 h−1 ) in the correlation plots.

Fig. 3 Residuals between the s-BART and CCF RVs for two nights of HARPS observations of HD 160691 (top row) and two ESPRESSO nights of Indi (bottom row). The red and blue points represents the RV time series extracted when using the corresponding detectors.

5 A closer look into the BERV-dependency

In Section 4, we show the presence of a quasi-linear, systematic, bias in RVs collected within a short period of time (i.e. over afew days). The main goal of this section is to understand the impact of relative and absolute BERV variations in the RV bias. To measure this dependence, we narrowed our analysis to observations of τ Ceti, which was intensively observed over five years in the context of ESPRESSO’s GTO.

5.1 Night-by-night analysis of τ Ceti

Selecting a subset of the available nights of τ Ceti observations we compute s-BART RVs using two different stellar templates:

1. We assumed that that each night is fully independent from the others and we constructed individual stellar templates for each. The RV time series obtained from this extraction is denoted as “single-night.”

2. One single stellar template, constructed from all available observations, as used in Figueira et al. (2025).

To ensure consistency between the different nights, we handled (as before) the telluric features through a common mask. We created a binary telluric mask that rejects wavelength regions with a drop of transmittance greater than 1%, which is then expanded to account for the maximum yearly BERV variation (which for τ Ceti is ~29.05 km s⁻¹). This guarantees that the RV extraction is always informed by the same wavelength regions, independently of the position of the telluric features. Similarly to before, we evaluated the magnitude of the bias through a linear fit to the residuals.

Figure 4 shows that the systematic bias is closely linked to the observations that were selected to construct the stellar template. If the template is constructed from observations that are spread over the year, we do not find any structure in the residuals (top right panel), nor any meaningful slope within the error bars (median value of ~0 ms⁻¹ h⁻¹ ). In the case where the nights were assumed to be independent (left column), each with its own stellar template, there is a linear structure in the residuals, with the maximum slope reaching −1.5 ms⁻¹ h⁻¹. This is a clear indication that the origin of this systematic bias is in the construction of the stellar template and the observations that are selected to be part of it.

Once again, we can confirm that the linear bias has a negative slope and that it is manifested with different amplitudes in different nights. We can also see that even within the data of the same target, we find a difference of 1 m s⁻¹ on the residual slope. Furthermore, we also fail to find (Appendix D) any correlation with the meteorological conditions of the individual nights.

Fig. 4 Top: residuals between s-BART RVs and CCF RVs of τ Ceti when using the template constructed from a single night (left) and the one constructed from all available observations (right). Each night is represented by a unique color in the two time series. Bottom: histogram of the slopes of the first-degree polynomials adjusted to the residuals of the corresponding column.

5.2 Binning observations based on their BERV

To disentangle a dependence with time from a dependence with BERV (which, within one night would be equivalent), we bin all of our observations based on their BERV. We constructed bins of 600 m s⁻¹, a value close to the expected BERV variation found within a single night of observations. For this, we applied the following prescription:

1. To avoid having a mixture of red-shifted and blue-shifted telluric features, within a given bin, we removed all observations for which the derivative of the BERV (over multiple nights) is positive. The resulting dataset is highlighted in Figure A.2 and ensures that within a given bin, the telluric features are all drifting in the same “direction.”

2. We divided our observations into the two ESPRESSO datasets - ESPRESSO18 and ESPRESSO19 - as they were treated independently during the RV extraction process. This division of ESPRESSO data is commonly used and motivated by an instrumental intervention that translated into a change of the spectral profile¹⁰.

3. We found the minimum BERV of each of the datasets and constructed intervals of 600 m s⁻¹ starting from it.

4. We iterated over all bins and discarded those with less than 20 observations. This ensured that we had a reasonable amount of observations to construct the stellar model and to draw any meaningful conclusion on the presence of a systematic bias.

We were left with 6 ESPRESSO18 bins and 24 ESPRESSO19 bins, which we applied s-BART to, following the same prescription that has been used throughout the paper. The results from this application are shown in Figure 5, where we show the residuals as a function of hourly difference to the first observation of the corresponding night. Once again we found a systematic bias within the nights that are included within the bins, with peak-to-peak amplitudes comparable with our previous results. In contrast to previous results, we find that some of the bins present a positive slope, albeit without any correlation with the median BERV of the corresponding bin. Furthermore, we find no obvious difference between the bins that only contain data collected during the same year (squares from the bottom panel) or contain observations collected over multiple years (circles).

Fig. 5 RV extraction on the ESPRESSO dataset of τ Ceti, after binning observations based on their BERV. Each bin is attributed a unique color. Top: residuals between the RV time series of s-BART and the CCF. To ease visualisation, we plot the residuals as a function of the time difference between the observation and the first observation of the night it was taken in. Bottom: slope of the first degree polynomial that is fitted to the RV residuals, as a function of the median BERV in the corresponding bin. The squares represent the bins that only contain data for one year, whilst the circles have two or more years of observations.

6 On the origin of the RV bias

Section 5 shows how the selection of observations for the construction of the stellar template plays a key role in the presence (or lack) of the RV bias. The templates constructed with the full τ Ceti dataset lead to no meaningful RV slopes in individual nights (median value of 0.0 m s⁻¹ h⁻¹), whilst the construction of the template from the data of the same night leads to a median slope of −0.53 m s⁻¹ h⁻¹. This leads us to think that the origin of this effect comes from the difference between datasets collected within a single night or over time.

If we focus on the fact that the BERV separation in singlenight observations is on the same order of magnitude as the mean pixel size of ESPRESSO we can postulate two sources of contamination on the stellar template: micro-telluric features and/or any detector fixed pattern noise. We know that in multinight observations microtellurics can - independently of the RV extraction method used - impact RV measurements at the 10 cm s⁻¹ level (Cunha et al. 2014) and bias our RV extraction. If we consider their impact on the stellar template, we see that they would be (assuming a large BERV span and an high number of observations) averaged out of the stellar template. They are, by definition, shallow features and when combining the multiple observations to construct the template, they will be averaged out of it. In the same vein, any instrumental systematic (or correlated noise) that is coherent over a timescale of a few hours will also end up averaged out of the stellar template in the case of multi-night observations.

In single-night observations, those assumptions might no longer hold. Since the BERV separation within our dataset is small, any contamination on the detector referential, either from tellurics or from artifacts on the detector, will end up approximately in the same position and overlap over time. The BERV variation during one night is smaller than the RV coverage of one pixel and than the full width at half maximum (FWHM) of the spectral lines. Thus, when stacking our observations (through a weighted mean or equivalent), they are not averaged out. In the specific case of microtelluric lines, they will be appear in the template slightly broader and shallower than in any individual observation. In the case of fixed pattern noise on the detector, the stellar template will contain the mean effect across the night, also blurred by the BERV variation. Afterwards, in the RV extraction process, those features and any other instrumental contamination could be driving the alignment between spectra and template since they will be shifting in the same "direction" in the spectra during the night.

This contamination from microtellurics and correlated noise would provide an explanation of the quasilinear shape of the bias and the fact that all observations present a negative slope (Section 4.1), as the light moves in the detector in the same direction due to the negative slope of the BERV variation over the night. It would also explain the fact that the bias is also present in SERVAL and ARVE RV measurements, as they also utilised a stellar template constructed from stacking multiple observations. It is also possible to justify the lack of effect on CCF RVs through a combination of two different factors: (i) the CCF method uses a fixed mask through time, leading to a model (the CCF mask) that is not contaminated; and (ii) the CCF mask only includes a small subset of the spectral lines (and regions) that are present in the stellar template. Overall, this hypothesis seems to be compatible with the results that we report in this work:

• In Section 4.2, with a chromatic analysis, we report a greater effect in the red wavelengths, where telluric features are more common.

• In Section 5.2, we binned τ Ceti observations based on their BERV and found the RV bias, even when combining observations of multiple years. This gives weight to the hypothesis that the RV bias is connected with the physical position of the spectral lines on the detector. We also saw that increasing the number of observations that are used to construct the stellar model (and consequently increase the BERV span) leads to a decrease of the RV slope, until we achieve convergence towards a slope of zero, as we would find in multi-night datasets.

• We find that a time separation between single-night spectra and the respective template leads to a decrease in the effect (Appendix E). This is still compatible with our two sources of contamination, since the BERV separation will: (i) place the micro-tellurics of the template further way from the ones in the spectra; and (ii) lead to a time evolution of the correlated noise sources, diluting its impact on the RV extraction as it wold no longer match the effect imprinted on the first night.

In Sections 6.1 and 6.2, we further explore the hypothesis of telluric contamination in the data, based on changes in the telluric masking and the BERV of the observations.

6.1 Impact of the telluric threshold

Here, we study the effects of telluric contamination by changing the threshold to reject telluric-contaminated regions. A decrease in this threshold will result in a higher number of (deeper) features bypassing our rejection criteria and, consequently, increase this contamination of our observations. Thus, if the RV trend is introduced by leftover tellurics in the stellar template, as described in Section 6, we would expect the amplitude to increase as our rejection criterion gets more relaxed.

To evaluate this hypothesis, we restricted our nights even further: from the subset of nights presented in Table A.2, we selected those for which the last observation was taken more than three hours after the first one. This choice leaves us with three nights of data and its selection was motivated by visualisation purposes: the systematic bias is clearer over longer time-scales. Then, we extracted RVs in the following fashion:

1. We constructed a new binary mask using a given threshold (Tr_i) for the telluric contamination. It is important to note in here that s-BART uses this threshold as a percentage of the continuum level. If Tr_i = 1%, it rejects every wavelength in which the continuum varies more than 1%. This means that higher values of Tr_i reject less telluric features than lower ones. For further details on the construction of the telluric mask, we refer to Silva et al. (2022).

2. For each night in our restricted set (i.e, the three nights selected based on having a time span greater than three hours), we used the telluric template from point 1 and we constructed a new stellar model. Then, we applied the usual s-BART RV extraction to get a new RV time series.

We repeated this process for a list of seven different thresholds - with Tr_i taking the values of {0.5, 1, 5, 10, 15, 20, 50}% - and compared the RV time series with the one that is extracted from the CCF, as shown in Figure 6.

From Figure 6, we find that a relaxation in the masking of the telluric features (increasing the threshold) leads to an increased effect in our RVs, which subsequently leads to a greater difference to the ones from the CCF. This is true for all cases, except for a telluric rejection of 0.5%, where the slope is slightly larger than for a threshold of 1%. With such a small threshold (half of s-BART default value of 1%), we end up with a massive rejection of wavelength regions, leading to a noisier time series that still leads to a slope that is compatible at the 1-σ level.

If we go on to focus on the difference between the results obtained with a telluric mask with Tr_i of 50% and the one with Tr_i of 1% we find an interesting result. On the two nights from November, we found a decrease in the RV slope of ~40%, whilst on the night from the December the telluric threshold variation only translates to a decrease of ~27%. Furthermore, the night from December also presents the smaller peak-to-peak variation, the larger time-span between start and end of observations, and the larger gap in the middle of the observations. This results seem to be compatible with what one would expect from the presence of telluric features. The wide gap in the data from December eases the “averaging-out” of the noise sources from the template. Thus, the extraction is not as sensitive (in comparison with the other nights) to their presence.

Independently of the threshold we used, we still find the quasi-linear trend between the two RV time series. In Appendix F, we present an overview of the effect that is introduced on data that was corrected from telluric features, showing results that are compatible (at the 1-σ level) to those that are present in here.

Fig. 6 Impact on the RV trend from template matching after creating a telluric mask with different thresholds for the rejection of telluric-contaminated regions. The thresholds were applied to the same transmittance spectrum, generated through Telfit (Gullikson et al. 2014). Each row represents different nights of τ Ceti, chosen from within our subset specified in Table A.2. Left: residuals to the CCF RVs, color-coded based on the transmittance threshold; Right: slope of a first-degree polynomial adjusted to the residuals as a function of the transmittance threshold.

6.2 Impact of BERV span on recovered RVs

Throughout this section, we evaluate the plausible contamination of telluric features and correlated noise sources on the stellar template. If this were indeed the case, we would find that as we increase the temporal span of our dataset (i.e., the observations that are used to construct the stellar template), the bias would decrease. For this analysis, we went back to our subset of preselected nights (Table A.2) and we sequentially added more observations to the stellar template in the following fashion:

1. We selected the ith night and construct a stellar template using all observations taken until the end of that night (i.e. all observations since the start of the first night).

2. We extracted the RVs from all observations collected in the night of the 28th of December, 2020, whilst using the newly created stellar template. During this process, we used the telluric mask with a transmittance cutoff of 1%.

In Figure 7, we present the results from this comparison, where it is clear that the addition of new observations into the template leads to a significant decrease in the RV systematic bias. There, we see that the slope of the residuals quickly converges to a near-zero value as the stellar template includes observations from more than five nights. It is also important to note that the τ Ceti observations present some of the lowest trends in our datasets, which might mean that the results from Figure 7 represent the best-case scenario.

Fig. 7 Top: comparison of different RV time series of τ Ceti, extracted with stellar templates that sequentially utilise one more night of observations than the previous. Bottom: slope of a first degree polynomial adjusted to the residuals as a function of the number of nights in the stellar template.

6.3 A summary of the RV bias

We described in Section 6.1, a less conservative telluric mask leads to an increase in the effect. However, even with the most conservative threshold that we use (masking tellurics deeper than 1%) we were not able to fully reduce the RV slope to 0m s⁻¹ h⁻¹. However, it is also important to keep in mind that this test only probes the effect on the lines that are present on the transmittance model. Our results can be interpreted in different ways: (i) the models that we use to construct the telluric binary mask are missing lines; (ii) telluric contamination could introduce a similar effect to what we see, but it is not the root cause; or (iii) the effect that we are measuring is a combination of different sources of contamination of the stellar template. The first hypothesis is difficult to support based on our data, as both the telluric masking and telluric correction (described in Appendix F) display compatible results. The increase in the absolute RV slope as we relax the telluric threshold does not allow us to confidently claim the telluric features to be the root source of this issue. Ifwe consider a given fixed-pattern noise artifact that exists on the referential of the detector, it will introduce a flux contamination similar to the one introduced by telluric features. Over the course of one night, the physical drift of stellar lines on the detector leads to them being imprinted on the data, similarly to how telluric features are overlapped with the stellar spectra. Under this hypothesis, a given percentual masking of telluric features minimises the RV residuals introduced by them, still keeping the contamination from the fixed-pattern artifact. A relaxation of the telluric masking translates to a magnified effect since the unmasked tellurics would introduce a new systematic RV slope to the pre-existing one from a fixed-pattern artifact. Based on our current results and understanding, we postulate that this RV systematic trend has been introduced by a combination of micro-telluric features and other (unidentified) artifacts on the detector that are coherent over timescales of a few hours.

The fact that the bias decreases as we increase the time between the data and template (Figure E.2) is fully expected, since we are left with a larger separation between the systematics of the template and the data, breaking the coherence of the flux artifacts. Furthermore, a combination of observations from different nights (as shown in Figure 7) leads to a decrease of the RV trend, which is compatible with both sources of contamination. The larger BERV span results in a larger separation between the telluric lines across all observations, diluting their presence in the stellar template. Furthermore, the temporal separation can reduce the coherence of the fixed pattern systematics present on the spectra. Both of those effects contribute to a dilution of the contamination of the stellar template and, consequently, a decrease in the RV trend that we measure. In Appendix G, we discuss other possible sources of noise on single-night observations.

7 Discussion

In this paper, we present and explore a previously unidentified systematic RV bias on template-based RVs that is not found for CCF radial velocities. This effect is seen for several stars observed with different spectrographs (see Section 4.1) with residual RV (of TM compared with CCF) slopes in the interval [−0.3, −52] m s⁻¹ h⁻¹. Furthermore, we find that the amplitude of the RV residuals is not consistent across different nights of observations of the same target. The fact that this trend is found across different spectrographs, as well in different pipelines, hints towards a systematic effect at play that is introducing a bias in our RV extraction. More worrisome is that not only do we have to consider the impact on intra-night observations, but we also need to understand how it translates to RV extraction in multi-night observations. Fortunately, as shown in Section 5.1, this bias is either not present under such conditions, or exists at an amplitude below the noise floor of our instruments.

7.1 Impact on different scientific cases

In this paper, we explore a systematic RV signal that affects our intra-night measurements with state-of-the-art high-resolution optical spectrographs, such as ESPRESSO and HARPS. This contamination will have an effect on the RV measurements obtained for different scientific purposes. We discuss the impact on three of those areas below.

Detection and characterisation of exoplanets: The determination of precise RVs plays a pivotal role in the detection and characterisation of exoplanets (e.g. Castro-Gonzâlez et al. 2023; Mascareno et al. 2023; Gonzâlez Hernândez et al. 2024; Passegger et al. 2024). With these scientific goals, stellar spectra is usually acquired within a large window of BERVs and, consequently, template-matching RVs are not affected by this bias. The effect that we report has been averaged out from the stellar template (given a certain number of observations spread over different nights) to a level that cannot be measured on RV time series with current spectrographs. However, at shorter time spans of observations, this time-correlated effect could be interpreted as a planetary signal. Fully understanding this spurious effect and developing a robust enough model will be critical to exploit the stellar spectra from the next generations of spectrographs and to possibly reduce the current noise floor of our instruments.

Impact on transit and atmospheric characterization: Transit observations have become central to modern astrophysics, providing key insights into the architecture of planetary systems (Winn et al. 2010; Albrecht et al. 2022; Cristo et al. 2024) and enabling the study of exoplanet atmospheres via transmission spectroscopy (Casasayas-Barris et al. 2019; Ehrenreich et al. 2020; Azevedo Silva et al. 2022). A crucial method for determining the relative orientation between a star’s rotation axis and a planet’s orbital plane is the Rossiter-McLaughlin (RM) effect (Holt 1893; Rossiter 1924; McLaughlin 1924). Boué et al. (2013) demonstrated that the RV signal from the RM effect is sensitive to the method used for RV extraction. Therefore, understanding the application of template-based RVs in this context is essential, particularly in terms of determining how any spurious slope could affect the interpretation of the RM curve. In transmission spectroscopy, the standard approach for constructing absorption spectra involves creating a template (or master-out) from out-of-transit observations, which is then compared to in-transit spectra. However, it is possible that a spurious contamination of the master-out (and associated RV measurements) may contaminate in-transit absorption spectra, potentially introducing biases in the determination of planetary properties or leading to false species detections, especially when both the behind-the-planet stellar regions and the planet have similar velocities.

Asteroseismology impact: Prior to the mid-2000s, RV observations constituted the primary method for measuring convection-driven, solar-like oscillations in stars other than the Sun, leading to detections for a few select stars (see e.g. Kjeldsen et al. 2008). The advent of space photometry, with missions such as CoRoT (Auvergne et al. 2009), Kepler/K2 (Koch et al. 2010; Howell et al. 2014), and TESS (Ricker et al. 2014), has since revolutionised the field, enabling the asteroseismic analysis of hundreds of thousands of stars (Campante et al. 2016; Hon et al. 2021). Despite these advancements, there remained a notable detection void in cool main-sequence stars (i.e. cooler than ~5000 K). Owing to their low luminosities, K dwarfs display extremely small oscillation amplitudes (below the 10 cm/s level or, equivalently, a few parts-per-million), which are thus beyond the reach of current space-based photometric missions. A viable alternative is offered by the much higher S/N in RV observations over the typical oscillation timescales, which has been driving a resurgence of observing campaigns on K dwarfs using state-of-the-art spectrographs, such as ESPRESSO (Pepe et al. 2021) and KPF (Gibson et al. 2020), mounted on large-aperture telescopes (Campante et al. 2024; Hon et al. 2024; Li et al. 2025). However, the short-cadence (<~60 s) and continuous observations needed to properly sample these oscillations make it impractical to secure more than a handful of consecutive nights at the telescope, which necessarily translates into short observational baselines. It is clear throughout this paper that this systematic bias in the RV extraction presents a critical challenge in the utilisation of template-based RV time series in the context of asteroseismology studies.

7.2 Possible avenues to decrease the impact

Throughout the paper, we explore different avenues to mitigate the impact of this contamination of the stellar template. In Section 5.1, we describe how the construction of a template from a large number of observations - spread over multiple years - is one possible avenue to mitigate this systematic effect. This was also covered in Section 6.2, where the addition of new nights of observations leads to a gradual decrease of the RV slope. In the context of exoplanetary detection and characterisation of exoplanets with orbital periods with timescales greater than a few days, such observations will be collected over time, leading to a natural mitigation of this systematic bias. However, this is not the case for all observational campaigns, as some consist of a short number of nights with high cadence (e.g. asteroseismology, transit and atmospheric studies, characterisation of ultra-short period exoplanets). It is thus of paramount importance for such science cases to have observations collected over a larger window of BERVs, so that the construction of data-driven stellar models is not impacted by this effect.

The option of using observations collected after a certain time has passed (as detailed in Appendix E) also seems to lead to satisfactory results on the higher-S/N stars. However, it also seems to be accompanied by increased RV scatter. It is not clear if this excess of variability could also be a byproduct of the root source of the systematic bias. The interpretation of the lower S/N case is made more difficult by the shorter coverage ofHD 40307 (five nights, in comparison with the previous six) and the overlap in BERV seen between observations (see Appendix 2). Although we find a similar behaviour to the one from the high-S/N case, it comes with a delay of one day. Thus, we were not able to properly interpret the lower S/N case of the asteroseismology datasets.

In the case where all options above are not valid, it would also be possible to use a stellar model constructed with data from a different star. If we find a star that has similar properties to the one in which we are interested in, we can construct the stellar template with its data and then apply it to our original star (e.g. as described in Artigau et al. 2022). The core concept behind template-matching algorithms is that the stellar model should be a perfect match to the observations of our star. By using a model from a different star, we would be introducing differences between model and data, which would then introduce some extra variability among the results.

8 Conclusions

In this paper, we present a systematic RV effect that is contaminating template-matching (TM) and line-by-line (LBL) RVs in the specific case where all observations fall within a short time window (i.e. within afew days). Not only do we find this effect in different template-based pipelines (s-BART, SERVAL, and ARVE), but also on different stars and spectrographs, with slopes that vary (in our sample) in the interval [−0.3, −50] m s⁻¹ h⁻¹. This systematic effect presents a significant challenge for the usage of template-matching RVs for science cases that focus on extensive, intra-night, observations; for instance, in asteroseismology analysis, transit, and atmospheric characterisation.

From ESPRESSO and HARPS observations, we find that the contamination can be found on each detector, although it is systematically greater on the red detector. We also find that this effect either does not come into play for multi-night datasets or it exists under the current noise-floor of our instruments and analysis techniques. We managed to pin-point the source of this systematic bias to the construction of the stellar template. In particular, to the selection of observations from which it is constructed. We postulate that a combination of leftover micro-telluric features and correlated noise in the stellar template leads to the systematic bias of the RVs over the night. We found that it can be mitigated by either using a large temporal separation between the template and data (in the case where the template is also constructed from one single night) or by using a large number of observations collected over a period of several months.

Acknowledgements

We thank the helpful comments of Andrew Collier Cameron as the referee of this paper. This work was 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 also supported by FCT - Fundaçâo para a Ciência e a Tecnologia through national funds by these grants: UIDB/04434/2020 DOI: 10.54499/UIDB/04434/2020, UIDP/04434/2020 DOI: 10.54499/UIDP/04434/2020, PTDC/FIS-AST/4862/2020, UID/04434/2025. TLC is supported by Fundaçâo para a Ciência e a Tecnologia (FCT) in the form of a work contract (2023.08117.CEECIND/CP2839/CT0004). J.H.C.M. acknowledges further support from national funds through the FCT DarkMAGE project (grant ID: PTDC/FIS-AST/4862/2020) and from the project e-CHEOPS (PEA: 4000142255), funded by ESA/PRODEX. JIGH acknowledge financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) projects PID2020-117493GB-I00 and PID2023-149982NB-I00. ÉA is supported by the Trottier Family Foundation through the Trottier Institute for Research on Exoplanets (IREx). ODSD acknowledges support from e-CHEOPS: PEA No 4000142255. The INAF authors acknowledge financial support of the Italian Ministry of Education, University, and Research with PRIN 201278X4FL and the “Progetti Premiali” funding scheme. X.De acknowledges funding from the French ANR under contract number ANR-24-CE49-3397 (ORVET). This work is supported by the French National Research Agency in the framework of the Investissements d’Avenir program (ANR-15-IDEX-02), through the funding of the “Origin of Life" project of the Grenoble-Alpes University. KA acknowledges support from the Swiss National Science Foundation (SNSF) under the Postdoc Mobility grant P500PT_230225. J.L.-B. is funded by the Spanish Ministry of SCience, Innovation and Universities (MCIN/AEI/10.13039/501100011033) through grants PID2019-107061GB-C61, PID2023-150468NB-I00 and CNS2023-144309. We acknowledge financial support from the Agencia Estatal de Investigaciôn of the Ministerio de Ciencia e Innovaciôn MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through project PID2021-125627OB-C32, and from the Centre of Excellence “Severo Ochoa” award to the Instituto de Astrofisica de Canarias. X.Du acknowledges the support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement SCORE No. 851555) and from the Swiss National Science Foundation under the grant SPECTRE (No 200021_215200). This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. FPE 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. This work was financed by Portuguese funds through FCT (Fundaçâo para a Ciência e a Tecnologia) in the framework of the project 2022.04048.PTDC (Phi in the Sky, DOI 10.54499/2022.04048.PTDC). CJM also acknowledges FCT and POCH/FSE (EC) support through Investigador FCT Contract 2021.01214.CEECIND/CP1658/CT0001 (DOI 10.54499/2021.01214.CEECIND/CP1658/CT0001). ASM acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) projects PID2020-117493GB-I00 and PID2023-149982NB-I00. ARCS acknowledges support from FCT fellowship 2021.07856.BD; This work has received support in the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grant 51NF4_205606. The authors acknowledge financial contribution from the European Union - Next Generation EU RRF M4C2 1.1 PRIN MUR 2022 project 2022CERJ49 (ESPLORA).

Appendix A BERV coverage of the sample of stars

Fig. A.1 BERV values for the different nighs of the asteroseismol-ogy stars as a function of the temporal distance to the first observation of that star. Top and middle row: ESPRESSO observations of Indi and HD40307, respectively. Bottom row: HARPS observations of HD160691. The dashed red line highlights the BERV value of the last observation taken in the first night.

Fig. A.2 BERVs of the ESPRESSO observations of τ Ceti (left) and their distribution (right) in bins of 1 km/s. The black points (2122 in total) represent the observations that we keep for our analysis of Section 5.2, whilst the red points were discarded. The dashed blue line represents the BERV curve of τ Ceti.

Fig. A.3 BERV coverage of the subset of observing nights of τ Ceti we selected to explore the systematic RV bias.

Fig. A.4 Histogram of the BERV distribution of all stars in our ESPRESSO-WG2 sample. Multiple nights observations of the same star are shown with the same color.

Appendix B Searching for correlations with line-shape evolution

In Section 4, we focus our analysis on a asteroseismol-ogy dataset, containing ESPRESSO data of 2 stars (ϵ Indi and HD40307) and HARPS observations of another star (HD160691). Here, we present, in Figure B.1, a scatter plot between the RV residuals and the line-shape metrics. This comparison must, however, be interpreted carefully, as the variation of both the line-shape metrics and the RV residuals have different amplitudes in the different stars.

Across the three stars, we find no correlations with the FWHM and bisector span. As for the contrast, we find the same lack of correlations on Indi and HD160691 data. Strangely, that is not the case for HD40307 as we find a clear correlation between variation in contrast and RV residuals, with a correlation coefficient of 0.62. A closer look at this star reveals that the seeing gradually worsened during the night which, in turn, leads to a correlation between contrast and time. As the residuals that we measure are also correlated with time, we are left with a correlation between RV residuals and contrast.

Lastly, it is important to note that the results of this analysis does not paint a full picture of line variability. By definition, the CCF metrics are only sensitive to variations in the lines that are present on the CCF mask. It is still possible that we are being affected by variability on the lines that are not present on the mask and, consequently, we find no clear correlations.

Fig. B.1 Correlation between the s-BART and CCF RVs for HD40307 (ESPRESSO, left column), HD160691 (HARPS, middle column), and ϵ Indi (ESPRESSO, right column). In each row we plot the values from the activity indicators from the CCF: FWHM (top row), BIS SPAN (middle row), and contrast (bottom row). As we are including more than one night, we removed the median value of each indicator, so that they are centered at zero. We employ a uniform color map for all stars, representing the difference in time to the first observation of the corresponding night (the blue points always correspond to the start of the night and the yellow represents a different of 10.5 hours).

Appendix C Slopes measured in the large-sample analysis

In Table C.1, we present the slope of the first-degree polynomials that were adjusted to the large sample of stars that we selected, sorted by their absolute value. The slope was adjusted through a least-squares fit, using the inverse variance as weights¹¹. The uncertainty associated with this measurement was estimated through the diagonal of the covariance matrix, ignoring cross-correlations between parameters.

This analysis reveals the presence of large differences in the retrieved RV slope between different targets and, on top of that, between different nights of the same target (e.g. as seen in HD40307 or HAT-P-26).

Appendix D Searching for correlations with weather conditions

In this appendix, we explore a possible correlation between the meteorological conditions and the τ Ceti RVs extracted with the nightly templates, as presented in Section 5.1. From the headers of the ESPRESSO files we collect the measurements of: ambient temperature, relative humidity, and pressure; airmass; wind speed and direction. For ease of analysis we plot, in Figure D.1, the slope and the median value of each of those metrics.

Unexpectedly, we failed to find a meaningful correlation with any of the aforementioned metrics. This dampens our explanation of telluric features being the root source of this contamination. If the drift was indeed caused by micro-telluric features that are bypass our mask, we should expect a correlation with the metrics that influence their depth.

Fig. D.1 Correlation between the slopes measured in Figure 4 and the median value of meteorological metrics collected by ESO for the τ Ceti observations. In the cases where there are measurements at the start and end of the observation, we report the value at the beginning.

Appendix E Impact of the BERV separation on the extraction

Lastly, we evaluate whether this systematic contamination is coherent at longer timescales or if it exists only within a small time window. To do so, we constructed stellar templates from individual nights and then use them to extract RVs. This process was applied to the two ESPRESSO asteroseismology targets (e Indi and HD40307) and for τ Ceti observations. The former dataset provides a dense coverage of a few consecutive nights, and the later provides a larger BERV coverage, albeit with a sparser sampling within any given night.

Asteroseismology observations: To have a better understanding of the relation between the bias and BERV, we select the high-S/N observations of ϵ Indi and the comparatively lower S/N observations of HD40307. For each star, we constructed a template from the data collected on their first night and usde it to extract RVs from all others, maximising the BERV separation between the stellar template and the last observation of the dataset. Figure E.1 presents the results from this comparison, clearly showing that the systematic bias disappears as the temporal distance between data and model increases. If we start by focusing on the top row (data from ϵ Indi) we find that the decrease of the peak-to-peak amplitude of the residuals is accompanied by a slight increase in the dispersion. A comparison between the first two nights shows a decrease in the structure of the residuals. Strangely, the application to the third night leads to a slight positive trend and, after that point, it converges to values close to zero. The behaviour of HD40307 (middle row) differs from the one found on Indi, as the change to a positive trend only occurs one day later (on the 4th night). Furthermore, the convergence towards ’flat’ residuals is slower, as the trend persists in all cases. However, the peak-to-peak impact of this systematic bias is larger on this star, in comparison with Indi.

These conclusions are summarised in the bottom panels of Figure E.1, which shows the RV slope as a function of the separation between the template and the data, alongside with Pearson’s correlation coefficient (which is summarised in Table E.1). One possible explanation for the difference in behavior between the two cases is that the BERV variation for ϵ Indi is such that we there is no overlap between two consecutive nights. However, the long duration of the HD40307 nights of observations allow for a significant overlap in BERV between different days.

τ Ceti observations: The τ Ceti data allow us to expand on this analysis over a longer timescale. We selected the night of 4 December 2022 as a baseline. Then, we select all nights that were observed in a window of 50 days before and after it. In total, we collected 17 nights (including the baseline one), which were then used to construct single-night stellar templates. We used each of those templates to compute s-BART RVs for the baseline night.

Figure E.2 presents the evolution of the residual slope as a function of temporal separation between the observation and template. The slope converges to zero after 20 days of difference between the data and template and takes a positive value in some cases. The scatter in the residuals also converges to zero after about ten days of separation.

Fig. E.1 RV residuals between s-BART and CCF RVs for the ESPRESSO observations of Indi (in black, top row) and HD40307 (in blue, middle row), divided in the different consecutive nights. The data from the first night of observations was used to construct the stellar template that was then applied to the other nights. Bottom row: Slope of a first degree polynomial that fitted to the RV residuals of Indi (left) and HD40307 (middle). The right-most column presents the Pearson correlation coefficient of the residuals over time, with the error bars representing a 90% confidence interval.

Finally, it is important to note that this analysis of the three stars is not enough for us to be fully certain that a given temporal separation is enough to mitigate most of the effect of this RV contamination. We have previously seen that the amplitude of this parasitic effect has a large variation across different stars (with a more pronounced effect in lower S/N regimes). In the analysis of these three stars, we have also seen that the decrease of the RV slope occurs at different rates for different stars.

Appendix F The effect of telluric correction

The latest release of ESPRESSO pipeline has introduced telluric correction for water features, following the methodology of Allart et al. (2022). As we discussed throughout the paper, the hypothesis of telluric contamination would imply that we have micro-telluric features that are not accounted for in our synthetic transmittance profile. One avenue to test this is to apply the same s-BART recipe - that we used in this work - to the telluric-corrected spectra. It is important to note that this new RV extraction will apply the telluric-masking on the telluric-corrected spectra. This means that the only difference at the level of the stellar spectra is in the micro-telluric features that bypass our binary mask.

Fig. E.2 Top: Residuals between the s-BART and CCF RVs of τ Ceti, color-coded by the difference in days between the baseline night and the night from which the template was constructed. Bottom: Evolution of two metrics with the distance (in days) between template and data: slope of a first degree polynomial (left) and the standard deviation of the RV residuals (right).

In Figure F.1 we present the residuals between the two s-BART RVs, showing an agreement at the 1-σ level. There is a slight difference between the two RV time series, but it is not significant. A follow-up analysis with telluric-correction on other molecules will be done after it is implemented in ESPRESSO’s DRS.

Fig. F.1 Comparison of the τ Ceti RV time series extracted from telluric corrected spectra for three nights explored in this paper.

Since we did not find any significant differences in the telluric-corrected spectra, it is even more unlikely that we are indeed missing telluric lines from our model.

Appendix G Other possible sources of RV biases

In this appendix, we postulate on other effects that could introduce trends within datasets that are collected within a single night:

1. A short-scale variability from the spectral continuum, that is biasing the stellar RVs. One of the main differences between CCF and TM algorithms is the regions of the stellar spectra that they use. In template-matching algorithms the stellar model will contain both spectral lines and the continuum, whilst on the CCF the model only contains the stellar lines (usually weighted by their RV information). Within s-BART we assume that the continuum of a given observation can have a linear difference to the stellar model. In other words, we marginalise over the hyper-parameters of a first-degree polynomial that describes the flux difference. If we have an higher-order variation on the flux levels over a night, they would not be accounted for. However, it is not clear to us why a change in the continuum level would translate into a multi-meter per second bias on RV time series.

2. A contamination between fiber A and fiber B of our instrument. In nominal operations¹², fiber A collects light from the science target, whilst fiber B will either be placed on the sky (to allow for the subtraction of its effect) or on a calibration light source¹³. It is not clear why such an effect would not appear on CCF RV and, on top of that, we find this bias on data collected with the two configurations of fiber B, namely, for simultaneous calibration or sky. This rules out the plausible (albeit unlikely) contamination between the two fibers.

3. The last effect that we can think of is a possible injection of a systematic effect of the stellar spectra from the calibration frames that are used by the instrument’s official pipeline. When the DRS is applied to the raw images, it selects the calibration frames closest to each image. This means that the DRS will reduce every observation (of the same night) using the same set of calibration frames, which is not the case in observations taken in different nights. A possible deviation of the instrumental properties, in comparison with the calibration frame, could be leading to an intra-night evolution of the stellar spectra. A preliminary analysis of some of the calibration frames (flat-fields and darks) revealed no obvious artifacts, but we leave a full analysis of those data products for future work. Furthermore, we find this effect not only on ESPRESSO and HARPS, but we are also aware of the same systematic bias on data from other high-resolution spectrographs¹⁴.

References

All Tables

All Figures

	Fig. 1 Mean-subtracted RV time series for one night (December 24, 2018) of ESPRESSO observations of HD 40307. The CCF RVs (blue triangles) are compared against those from s-BART (black), SERVAL (red), and ARVE (orange).
In the text

Fig. 2 RV bias for 3785 ESPRESSO observations of 20 different stars. The top row presents an histogram of the slopes of a first degree polynomial that was adjusted to the residuals between s-BART and CCF RVs. The bottom two rows present a scatter plot with the median value of different metrics: intra-night slope of BERV variation, effective temperature, median S/N at 600 nm, and median difference between CCF RV and BERV. All observation nights of the same star are shown with a common color. For representation purposes we do not show the two outliers (with a slope greater than −20 ms−1 h−1 ) in the correlation plots.

In the text

	Fig. 3 Residuals between the s-BART and CCF RVs for two nights of HARPS observations of HD 160691 (top row) and two ESPRESSO nights of Indi (bottom row). The red and blue points represents the RV time series extracted when using the corresponding detectors.
In the text

	Fig. 4 Top: residuals between s-BART RVs and CCF RVs of τ Ceti when using the template constructed from a single night (left) and the one constructed from all available observations (right). Each night is represented by a unique color in the two time series. Bottom: histogram of the slopes of the first-degree polynomials adjusted to the residuals of the corresponding column.
In the text

Fig. 5 RV extraction on the ESPRESSO dataset of τ Ceti, after binning observations based on their BERV. Each bin is attributed a unique color. Top: residuals between the RV time series of s-BART and the CCF. To ease visualisation, we plot the residuals as a function of the time difference between the observation and the first observation of the night it was taken in. Bottom: slope of the first degree polynomial that is fitted to the RV residuals, as a function of the median BERV in the corresponding bin. The squares represent the bins that only contain data for one year, whilst the circles have two or more years of observations.

In the text

Fig. 6 Impact on the RV trend from template matching after creating a telluric mask with different thresholds for the rejection of telluric-contaminated regions. The thresholds were applied to the same transmittance spectrum, generated through Telfit (Gullikson et al. 2014). Each row represents different nights of τ Ceti, chosen from within our subset specified in Table A.2. Left: residuals to the CCF RVs, color-coded based on the transmittance threshold; Right: slope of a first-degree polynomial adjusted to the residuals as a function of the transmittance threshold.

In the text

	Fig. 7 Top: comparison of different RV time series of τ Ceti, extracted with stellar templates that sequentially utilise one more night of observations than the previous. Bottom: slope of a first degree polynomial adjusted to the residuals as a function of the number of nights in the stellar template.
In the text

	Fig. A.1 BERV values for the different nighs of the asteroseismol-ogy stars as a function of the temporal distance to the first observation of that star. Top and middle row: ESPRESSO observations of Indi and HD40307, respectively. Bottom row: HARPS observations of HD160691. The dashed red line highlights the BERV value of the last observation taken in the first night.
In the text

	Fig. A.2 BERVs of the ESPRESSO observations of τ Ceti (left) and their distribution (right) in bins of 1 km/s. The black points (2122 in total) represent the observations that we keep for our analysis of Section 5.2, whilst the red points were discarded. The dashed blue line represents the BERV curve of τ Ceti.
In the text

	Fig. A.3 BERV coverage of the subset of observing nights of τ Ceti we selected to explore the systematic RV bias.
In the text

	Fig. A.4 Histogram of the BERV distribution of all stars in our ESPRESSO-WG2 sample. Multiple nights observations of the same star are shown with the same color.
In the text

Fig. B.1 Correlation between the s-BART and CCF RVs for HD40307 (ESPRESSO, left column), HD160691 (HARPS, middle column), and ϵ Indi (ESPRESSO, right column). In each row we plot the values from the activity indicators from the CCF: FWHM (top row), BIS SPAN (middle row), and contrast (bottom row). As we are including more than one night, we removed the median value of each indicator, so that they are centered at zero. We employ a uniform color map for all stars, representing the difference in time to the first observation of the corresponding night (the blue points always correspond to the start of the night and the yellow represents a different of 10.5 hours).

In the text

	Fig. D.1 Correlation between the slopes measured in Figure 4 and the median value of meteorological metrics collected by ESO for the τ Ceti observations. In the cases where there are measurements at the start and end of the observation, we report the value at the beginning.
In the text

Fig. E.1 RV residuals between s-BART and CCF RVs for the ESPRESSO observations of Indi (in black, top row) and HD40307 (in blue, middle row), divided in the different consecutive nights. The data from the first night of observations was used to construct the stellar template that was then applied to the other nights. Bottom row: Slope of a first degree polynomial that fitted to the RV residuals of Indi (left) and HD40307 (middle). The right-most column presents the Pearson correlation coefficient of the residuals over time, with the error bars representing a 90% confidence interval.

In the text

	Fig. E.2 Top: Residuals between the s-BART and CCF RVs of τ Ceti, color-coded by the difference in days between the baseline night and the night from which the template was constructed. Bottom: Evolution of two metrics with the distance (in days) between template and data: slope of a first degree polynomial (left) and the standard deviation of the RV residuals (right).
In the text

	Fig. F.1 Comparison of the τ Ceti RV time series extracted from telluric corrected spectra for three nights explored in this paper.
In the text