© R. Allart et al. 2022

1 Introduction

The Earth’s atmosphere is rich in spectral features. The molecules H₂O, O₂, O₃, CH₄, and CO₂ are the main molecules from 0.3 to 5 µm in the Earth’s transmission spectrum (also called the telluric spectrum). When astrophysical objects are observed with ground-based instruments, their light is absorbed by the Earth’s atmosphere, and thus their spectrum is contaminated by the telluric spectrum. At high spectral resolution, the telluric lines are well resolved and thousands of them can be distinguished. Moreover, at specific wavelengths, they completely absorb the light of the observed object, thus losing any astro-physical content. The absorption of the telluric spectrum varies as a function of the volume mixing ratio along the line of sight, which is controlled at first order by airmass. The exception for this is the water vapor, for which water column density and total content in the line of sight changes inhomogeneously across the sky and over a timescale of hours due to the weather conditions (Li et al. 2018).

With increasing instrumental performances and a larger telescope collecting area, developing methodologies to measure high-fidelity unbiased spectra is becoming a real challenge for several astrophysical topics such as an accurate stellar library, stellar activity indicators, precise stellar radial velocity (RV), exoplanet atmospheric characterization, or the measurement of the fine structure constant. Therefore, accurate correction of the Earth’s transmission spectrum is crucial in the process of obtaining high-fidelity results. Several methods are available for removing the telluric spectrum from complex model-based algorithms (e.g., Molecfit; Smette et al. 2015; Kausch et al. 2015), TAPAS (Bertaux et al. 2014), and TELFIT (Gullikson et al. 2014) to the observation of a spectrophotometric standard star close in time and in the sky to the scientific observation (Vidal-Madjar et al. 1986). Ulmer-Moll et al. (2019) compared these different approaches for near-infrared archival data obtained with CRIRES. The authors concluded that the model-based approach better corrects water telluric lines while the standard star is more appropriate for correcting the O₂ lines. Such methods can be applied to different astrophysical objects. Other studies are more focused on the detection of exoplanets through RV and on the characterization of their atmospheres. For example, Artigau et al. (2014) proposed an approach based on principal component analysis (PCA) built over a library of spectrophotometric standard stars for archival HARPS data. They obtained a better RV accuracy by including a wider telluric-corrected spectral domain when compared with the restricted telluric-free domain used. Empirical approaches, of which PCA is one (e.g., SYSREM; Mazeh et al. 2007), are also often used for exoplanet atmospheric studies and take advantage of the different RV components of the planet, the star, and the Earth (Snellen et al. 2010). In addition to PCA-based algorithms, forward modeling algorithms are used to study exoplanet atmospheres (e.g., Allart et al. 2017). However, forward modeling algorithms are not yet optimized to automatically produce precise RVs as their minimization algorithm can be biased by stellar contamination depending on the combination of the spectral type, systemic velocity, and barycentric Earth radial velocity (BERV). Allart et al. (2017) demonstrated this with Molecfit, which is very sensitive to the selection of telluric regions that are to be fitted, in order to perform the most optimal correction.

Precise RVs are obtained, with stable high-resolution spectrographs, by cross-correlating stellar spectra with a binary mask (Baranne et al. 1996; Pepe et al. 2002) or through template-matching algorithms (e.g., Anglada-Escudé & Butler 2012; Astudillo-Defru et al. 2015; Zechmeister et al. 2018). Telluric lines spoil RV measurements when they are not removed or corrected for, independently of the method used to measure RVs. Exclusion of the telluric lines can be quickly done, either by masking the complete bands or by masking each line manually on each spectrum (e.g., Astudillo-Defru et al. 2015). Only the telluric lines of O₂, H₂O, and OH affect measurements in the visible wavelength. Since these cover a comparatively small fraction of the wavelength range, discarding them from the analysis has a comparatively small effect. However, with the development of high-resolution spectrographs in the near-infrared (CARMENES Quirrenbach et al. 2014), SPIRou (Donati et al. 2020), and NIRPS (Wildi et al. 2017) or installed on larger collecting areas (ESPRESSO (Pepe et al. 2021)), the need for an automatic and homogeneous telluric removal tool became apparent. Moreover, Cunha et al. (2013), Lisogorskyi et al. (2019), and Xuesong Wang et al. (2019) noticed that micro-telluric lines (unresolved or very shallow lines, with a depth up to 2%), which are generally not excluded when building the cross-correlation function (CCF), can introduce systematics in the RV at a level of 10–20 cm s⁻¹, which is comparable to the expected RV precision of the state-of-the-art optical spectrographs (ESPRESSO (Pepe et al. 2021), EXPRES (Jurgenson et al. 2016)) and to the precision required in order to detect an Earth-Sun system analog.

In recent years, several techniques have been developed to overcome the limitation of (micro-) telluric lines to reach extremely precise RVs. Similar to Artigau et al. (2014), Leet et al. (2019) propose the self-calibrating, empirical, light-weight linear regression telluric (SELENITE) model. They use the variation in water telluric line depths as a function of the precipitable water vapor and the variation in non-water lines with the airmass to model the telluric spectrum. They do not require any line list database, but instead they use the observations of a few dozen B stars to build their SELENITE telluric database. The authors applied their correction to EXPRES data and reached average residuals similar to model-based corrections discussed in Ulmer-Moll et al. (2019). Artigau et al. (2021) propose a two-step approach where they apply the first correction with a TAPAS model scaled to the water content and airmass. Then, they use the PCA-based approach discussed in Artigau et al. (2014) on the residuals of the TAPAS correction. The authors applied their technique to SPIRou data (in the near-infrared, which is more polluted by telluric lines in comparison to the visible) of TOI-1278. Finally, Bedell et al. (2019) and Cretignier et al. (2021) also propose data-driven approaches with the Wobble and YARARA algorithms that use the time series of a given star to disentangle the contribution of the telluric spectrum from the stellar spectrum. Therefore, their methods can only be applied as a postprocessing step as they require a large number of spectra of the same object with a large enough variation in the BERV. With additional analysis treatment, the authors are able to reduce the RV root mean square (rms) by 10–20% in comparison to the classical reduction for HARPS observations.

Figure 1 shows where the telluric lines fall on top of G- and M5-type stars observed with the ESPRESSO spectrograph. Both H₂O and O₂ bands are present across the full spectral range of ESPRESSO but their intensities increase toward redder wavelengths. The ESPRESSO data reduction software (DRS), which builds upon the HARPS DRS, excludes the wavelength ranges (gray bands) of the telluric lines and the active stellar lines (Ca ii H & K, Na doublet, and Balmer lines) when the CCFs are built. This drastic measure masks up to ~27% of the spectral range of ESPRESSO, which is essential for analyzing late-type stars as their blackbody maxima are toward redder wavelengths. Moreover, lines at redder wavelengths are known to be less impacted by stellar activity (Suárez Mascareño et al. 2020). Micro-telluric water bands are not excluded when the CCFs are built (50005900 Å). The impact of these lines on the RV content of HARPS was considered negligible in comparison to the number of stellar lines present at the same positions. However, the collecting area of the Very Large Telescope (VLT), the wider spectral range to redder wavelengths, and the improved instrumental stability of ESPRESSO increase their impact, arguably to a detectable level, as demonstrated by Cunha et al. (2013).

The aforementioned telluric correction techniques are not practical for integration into an automatic pipeline to reduce individual frames and to be used by the whole community. In addition, having an additional database of telluric standards and complex forward models that require manual intervention are clear limitations. Therefore, we explore in this paper how a simple telluric model could rapidly correct the spectra of astro-physical objects with a minimal number of assumptions. The telluric model is based on a line-by-line radiative transfer code that assumes a single atmospheric layer. It only requires information on the sky observations, knowledge of the instrumental resolution, and a line list database such as HITRAN (Rothman et al. 2013). This algorithm will be implemented in the official DRS of ESPRESSO. Therefore, the model is built for the S2D_BLAZE_A spectra of ESPRESSO from which the CCFs are derived. The S2D_BLAZE_A product is the echelle-extracted spectra in the Earth’s rest frame produced by the ESPRESSO DRS before removing the blaze to preserve the relative lines’ strengths. We present a special focus on the telluric correction of water lines, which are the most numerous and variable, and therefore the most problematic telluric lines in the spectral range of ESPRESSO.

In Sect. 2, we describe our one-layer telluric model. In Sect. 3, we describe how the telluric model parameters are adjusted to the data. In Sect. 4, we apply the telluric correction to different spectral types. Section 5 presents the implications of the telluric correction of water for the RV time series of Tau Ceti and Proxima obtained with ESPRESSO. Afterward, we open the discussion to telluric correction of the O₂ molecular bands with its implication for other molecules and spectral domains in Sect. 6. We conclude in Sect. 7.

Fig. 1 Telluric and stellar spectra over the ESPRESSO spectral range. The gray bands are the regions excluded from the mask used in the ESPRESSO DRS. Top panel: telluric spectrum produced with Molecfit (Smette et al. 2015; Kausch et al. 2015). The H2O lines are in blue, while the O2 lines are in red. Middle panel: ESPRESSO spectrum of the G-type star Tau Ceti. Bottom panel: ESPRESSO spectrum of the M5-type star Proxima.

2 Model description

To build the Earth’s transmission spectrum, one needs to consider that its atmosphere is a gas that attenuates the stellar light as a function of the optical path length. Following the description on the HITRAN website¹, the telluric spectrum T(ν) follows the Beer-Lambert law and can be expressed as a function of the sum of the opacity function t(v) of each telluric line. The telluric spectrum is built in wavenumber units, ν, in cm⁻¹: $$T\left(v\right)=e^{-\sum \tau \left(v\right)}.$$(1)

The opacity function of a single line depends on the absorption cross-section σ(ν) in cm² molecule⁻¹, the number density n(z) in molecule cm⁻³, and the path length l in centimeters: $$T\left(v\right)={\displaystyle {\int }_0^l\sigma }\left(v\right)\cdot n\left(z\right)\cdot \text{d}z.$$(2)

If a single atmospheric layer at constant pressure and a homogeneous distribution of water vapor is considered, we can see that $${\displaystyle {\int }_0^{l}n}\left(z\right)\cdot \text{d}z=n\cdot l,$$(3)

and $${\text{IWV}}_{\text{LOSS}}=\frac{n\cdot l}{n_{{\text{H}}_2\text{O}}},$$(4)

where $n_{{\text{H}}_2\text{O}=3.34\times {10}^{22}}$ molecules cm⁻³, is a constant expressing the number of water molecules per unit volume of liquid water, and IWV_LOS is the integrated water vapor in the line of sight, which represents the condensed water column in millimeters.

Therefore, Eq. (1) can be rewritten as $$T\left(v\right)={\text{e}}^{-{\text{IWV}}_{\text{LOSS}}\cdot n_{{\text{H}}_2\text{O}}\cdot \sum \sigma \left(v\right)}.$$(5)

The absorption cross-section of each line can be expressed as the line intensity $s_{ij}^{*}$ times the line profile, f(v): $$\sigma \left(v\right)=S_{ij}^{*}\cdot f\left(v\right),$$(6)

where $s_{ij}^{*}$ is the line intensity at a given temperature, T, and is linked to the reference line intensity S_ij at ambient temperature (296 K): $$S_{ij}^{*}=S_{ij}\cdot \frac{Q_{296\text{K}}}{Q_T}\cdot \frac{exp\left(-c_2\cdot E”/T\right)}{\exp \left(-c_2\cdot E”/296\right)}\cdot \frac{1-\exp \left(-c_2\cdot v_{ij}/T\right)}{1-\exp \left(-c_2\cdot v_{ij}/296\right)}.$$(7)

Similarly, the line positions, $v_{ij}^{*}$, are subject to pressure variations, where is the pressure shift and P is the gas pressure: $$v_{ij}^{*}=v_{ij}+{\delta }_{\text{air}}\cdot P.$$(8)

To build a telluric spectrum that we will be able to use for the echelle-extracted spectra produced by the ESPRESSO DRS (S2D hereinafter), we build the model by spectral order. The wavenumber grid has a resolution of 3 000 000 and the grid is larger than the S2D grid by 40 km s⁻¹ (1Å at 7500 Å) on each side of the orders to account for the BERV that shifts the spectrum on the detectors by, at most, 30 km s⁻¹. The wavenumber grid is then converted into angstrom and convolved with a Gaussian instrumental profile to produce the spectrum as would be seen by ESPRESSO. The width of the convolving Gaussian corresponds to the width of the Gaussian on the previously produced spectral resolution map (discussed in Sect. 2.4). Finally, the last step consists in resampling the convolved telluric model to the S2D grid.

2.1 Line profiles

Different line profiles to model the telluric lines have been considered: Gaussian, Lorentzian, and Voigt profiles.

2.1.1 Gaussian profile

The Gaussian profile is the theoretical profile describing the thermal broadening of a line and is expressed as $$f_{\text{Gaussian}}\left(v\right)=\sqrt{\frac{\ln \left(2\right)}{\pi \cdot {\alpha }^2}}\cdot \exp \left(-\frac{{\left(v-v_{ij}^{*}\right)}^2\cdot \ln \left(2\right)}{{\alpha }^2}\right),$$(9)

where α is the Gaussian half width at half maximum (HWHM) that depends on the Avogadro number N_A = 6.02 × 10²³ mol⁻¹, the Boltzmann constant k_b = 1.38 × 10⁻²³ kg m² s⁻², the gas temperature T, the speed of light c, and the molar mass of water, M_mol = 18.01 g mol⁻¹: $$\alpha =\frac{v_{ij}}{c}\cdot \sqrt{\frac{2\cdot N_A\cdot k_{\text{b}}\cdot T\cdot \ln \left(2\right)}{M_{\text{mol}}}}.$$(10)

The Gaussian profile depends mainly on the temperature and not on the pressure. The Gaussian profile is expected to dominate in the upper part of the Earth’s atmosphere, where pressure broadening is negligible.

2.1.2 Lorentzian profile

The Lorentzian profile is the theoretical profile that describes the pressure broadening of a line: $$f_{\text{Lorentzian}}\left(v\right)=\frac{1}{\pi }\cdot \frac{\gamma }{{\gamma }^2+{\left(v-v_{ij}^{*}\right)}^2},$$(11)

where γ is the Lorentzian (pressure-broadened) HWHM that depends on the pressure P, the temperature T, the partial pressure P_self, and the reference temperature T_ref = 296 K: $$\gamma ={\left(\frac{T_{\text{ref}}}{T}\right)}^{n_{\text{air}}}\cdot \left({\gamma }_{\text{air}}\cdot (P_{\text{mean}}-P_{\text{self}}+{\gamma }_{\text{self}}\cdot P_{\text{self}}\right).$$(12)

In this case, P_self can be neglected and the HWHM can be simplified to $$\gamma ={\left(\frac{296}{T}\right)}^{n_{\text{air}}}\cdot \left({\gamma }_{\text{air}}\cdot P_{\text{mean}}\right).$$(13)

The Lorentzian profile depends mainly on the pressure and not on the temperature. The Lorentzian profile is expected to dominate in the lower part of the Earth’s atmosphere, where pressure is higher and, coincidently, where most of the water vapor is located.

Fig. 2 Reversed cumulative distribution of the water vapor line intensities in the water band ranging from 3800 to 8000 A.

2.1.3 Voigt profile

The Voigt profile is simply the convolution of the Gaussian and Lorentzian profiles. Formally, it is the profile that should be used to reproduce telluric lines because all lines have a contribution of thermal and pressure broadening: $$f_{\text{Voigt}}\left(v\right)-f_{\text{Lorentzian}}\left(v\right)*f_{\text{Gaussian}}\left(v\right).$$(14)

2.2 HITRAN line list

We use the HITRAN (Rothman et al. 2013; Gordon et al. 2017) database to build the telluric spectrum. Lines from 3800 to 8000 Å for H₂¹⁶O that have a line intensity stronger than 10⁻²⁶ cm⁻¹/(molecule cm⁻²) are selected. If too few lines are included, the model will not be representative of the observation. When too many lines are included, the weakest lines that do not contribute to the telluric spectrum increase the computation time, as they are the most numerous (see Fig. 2).

The cut-off at 10⁻²⁶ cm⁻¹/(molecule cm⁻²) was obtained after including more lines or less lines in the model (Fig. 3), and by looking at the dispersion and peak-to-valley (P2V) amplitude of the residuals. It ensures that we model all telluric lines that would be seen in spectra with S/N > 10⁵ per pixel, much higher than typical observations, while keeping computing times reasonable.

The following parameters of each line have been extracted: v_ij (vacuum wavenumber in cm⁻¹), S_ij (line intensity in cm⁻¹/(molecule cm⁻²)), γ_air (air-broadened HWHM in cm⁻¹ atm⁻¹), γ_seif (self-broadened HWHM in cm⁻¹ arm⁻¹), n_air (temperature dependence exponent for γ_airin cm⁻¹), (air-pressure-induced line shift in cm⁻¹ arm⁻¹) and E” (lower state energy in cm⁻¹). Another important parameter from the HITRAN database is the table of the total partition sum value, Q_T, for H₂¹⁶O, which assigns a Q_T value to each temperature.

Fig. 3 Impact of the number of lines selected in the telluric model. Top panel: average water telluric line built for the model, with line intensity thresholds at 10−24, 10−25, 10−26, 10−27, 10−28, and 10−29 cm−1/(molecule cm−2). Middle panel: residuals of the average lines with respect to including all the lines. See Sect. 3 for more details on how the average line is built. Bottom panel: zoom on the residuals obtained with the model that includes the 10−26, 10−27, 10−28 cm−1/(molecule cm−2) strongest water lines. The dispersion and P2V amplitude of the residuals are reported in the legend.

2.3 Sky conditions

The ESPRESSO FITS header contains the following parameters linked to the target, the local weather station, and the ESO radiometer (Kerber et al. 2012): the local temperature T in K, the ground pressure P in hPa (it has to be divided by 1013.2501 to have it in atm), and the integrated water vapor² IWV measured at the zenith in millimeters, which represents the condensed water column. The IWV in the line of sight can be approximated by $${\text{IWV}}_{\text{LOS}}=\text{IWV}\cdot \text{Airmass}.$$(15)

2.4 Resolution map

The next step consists in convolving the telluric model with the appropriate instrument line spread function (LSF) to reproduce, in the best possible way, the observed spectrum. The instrumental LSF can be approximated by a Gaussian with its full-width half maximum (FWHM) equal to the instrumental resolution. However, the instrumental resolution is not constant across the spectrum. It depends on the optics used to project the light on the detectors, and on potential mechanical tensions on the detectors. To characterize the instrumental resolution of ESPRESSO across the detectors (Pepe et al. 2021), the light of a thorium-argon lamp, used for calibration, is injected into the science fiber. Thorium lines are spread across both detectors and are present at several pixel positions in each order. As they are ultra-narrow lines, they are not resolved, and thus their FWHM is a direct measure of the instrumental resolution. Consequently, it is possible to map the detector and derive a resolution map of ESPRESSO for the different observing modes. To obtain a resolution value at each pixel of each order, we applied a two-dimensional polynomial fit, where x corresponds to the pixel of each order y: $$R=a\cdot x+b\cdot y+c\cdot x\cdot y+d\cdot x^2+e\cdot y^2+f\cdot x^2\cdot y^2+g.$$(16)

This parametrization was found to best reproduce the instrumental resolution measured on the unresolved Thorium lines. The coefficients can be found in Appendix A.

Figure 4 shows the resolution map (orders versus pixels) for the high-resolution mode using the 1 x 1 binning. As described in Pepe et al. (2021) and due to the high resolution of ESPRESSO combined with the large collecting area of the VLT, a pupil sheer is used to reduce the spectrograph size. Therefore, each physical order is duplicated on the detectors into two slices. The dots are the measured resolution on the unresolved thorium lines, while the background map is the fit applied to those lines. The blue detector corresponds to slices from 0 to 89, while the red detector is for slices 90 to 169. The median resolution is about 138000, but variations are visible between both detectors. The red detector has a higher resolution than the blue one globally but shows some variations toward the bottom left edge. The resolution of the blue detector decreases toward the top left and bottom right edges. Those variations reflect the instrumental LSF but also opto-mechanic adjustments applied to the detectors. It is also important to note that there is some spectral overlap between consecutive orders, and thus the right edge of one order probes the same wavelength as the left part of the order above it. In other terms, one wavelength can be present in two consecutive orders (four slices) with different flux values, but also with different resolutions.

Fig. 4 Resolution map (orders versus pixels) for blue and red detectors in the HR1 × 1 mode. The dots are the measured resolution on individual thorium lines at known positions in pixels and in orders. The background map is the polynomial fit applied to the data.

3 Telluric model fit to ESPRESSO data

The algorithm developed to automatically correct the telluric lines consists of adjusting a minimal number of parameters able to describe the spectra in a short computation time³. To do so, the parameters are adjusted on an average telluric line through a CCF, which works as a de facto average line. The use of the CCF has the advantage of reducing stellar contamination, which leads to a more robust determination of the telluric parameters. The produced average line is specific to each molecule and is built over what we refer to hereafter as the subset of telluric lines. The subset is built such that it consists of the strongest unsaturated lines in a defined wavelength range and such that no lines are closer than 80 km s⁻¹ to each other. In the case of H₂O in the spectral range of ESPRESSO, we consider the subset to be composed of the 20 strongest lines (see Appendix B) after excluding the lines at 7186.51,7208.41, and 7236.72 A, which are removed as they can be saturated depending on the atmospheric water content.

The algorithm starts by identifying the orders or slices (with a spectral coverage reduced by 160 km s⁻¹ on each edge⁴) where the subset of telluric lines falls. In the case of water vapor for ESPRESSO, the orders are 77 and 78 (slices 154, 155, 156, and 157). For each order or slice, the HITRAN parameters are extracted for the lines of the subset present in the wavelength range and for all lines in a passband of 80 km s⁻¹ around them. Then, the telluric spectrum is computed as in Sect. 2. The next step consists in computing the CCF with a telluric binary mask covering the wavelength range. The mask is composed of the subset of telluric lines, for which their weight is fixed to one. The CCF is computed over both the observed stellar spectrum and the telluric model, afterward called observed and modeled telluric CCFs. These observed and modeled telluric CCFs are co-added over the selected orders or slices to obtain a master observed telluric CCF and a master modeled telluric CCF. Both CCFs are then normalized by their continuum. A LevenbergMarquardt algorithm is then applied to minimize the residuals between the normalized master observed telluric CCF and the normalized master modeled telluric CCF (the free parameters of the model are discussed below). Finally, the best-fit parameters are used to compute the telluric model over the full spectral range of the observation, as in Sect. 2.

3.1 Impact of line profile

It is essential to assess the impact of the line profile and to select the one that best models the data with the fewest parameters. The fitted parameters are T and IWV_LOS for the Gaussian profile (Eq. (9)), and T, P_mean, and IWV_LOS for both the Lorentzian (Eq. (11)) and Voigt (Sect. 2.1.3) profiles.

Figure 5 compares the normalized master observed telluric CCF obtained with the subset of telluric lines (in blue) to the master modeled telluric CCFs based on the three different profiles studied. The Gaussian profile is not able to reproduce the observed average water telluric line because it does not account for the pressure broadening. The two other models reproduce the data well, with residuals below the 2% P2V amplitude (the residual shape is discussed in Sect. 4.1). This threshold of 2% P2V will be used to assess the data correction quality in the following tests. This choice represents the threshold for micro-telluric lines that are not historically masked by the DRS of HARPS or ESPRESSO. In the case of water vapor, the contribution of the Gaussian profile in the Voigt profile is negligible, and to follow Occam’s razor principle, the Lorentzian profile was selected as it is the simplest one.

Fig. 5 Impact of the line profile to best model the data. Top panel: CCFs computed over the 20 strongest water lines. The blue CCF is the master observed telluric CCF based on ESPRESSO data, while the orange, green, and red CCFs are the master modeled telluric CCFs assuming Gaussian, Lorentzian, or Voigt profiles, respectively. Bottom panel: residuals only for the Lorentzian and Voigt profiles. The Gaussian profile is not shown due to its larger residuals.

3.2 Impact of the temperature

Figure 6 (top panel) shows the impact of the temperature on the Lorentzian line shape, while P_mean and IWV are respectively fixed to half the ground pressure and IWV_LOS. The temperature was set at the lowest (273 K), the median (282 K), and the highest (293 K) temperatures registered since the start of ESPRESSO operations. The temperature does not have any impact on the line shape and it can be fixed to the ground temperature measured at the time of the observations.

3.3 Impact of the pressure

Figure 6 (middle panel) shows the impact of the pressure on the Lorentzian line shape, while T and IWV are respectively fixed to the ground temperature and IWV_LOS. The pressure was set to the ground pressure, half of the ground pressure, and twice the ground pressure. When the pressure decreases, the line shape is narrower and deeper. Contrary to the temperature, the pressure plays an important role in the line shape and the algorithm fits it.

3.4 Impact of the integrated water vapor

Figure 6 (bottom panel) shows the impact of the IWV on the Lorentzian line shape, while T and P_mean are respectively fixed to the ground temperature and half the ground pressure. The IWV was scaled according to airmass, IWV_LOS, half of IWV_LOS, and twice IWV_LOS. When the water content increases, the line is deeper. The IWV affects the line as a scaling factor, and it is thus necessary to fit it.

Fig. 6 Impact of T (top), P (middle), and IWVLOS (bottom) on the telluric line shape. Top: master modeled telluric CCF at an atmospheric temperature of 282 K (blue), 293 K (orange) and 273 K (green). Pmean and IWV are fixed to the same values for the three models. Middle: master modeled telluric CCF at ground pressure (blue), half of the ground pressure (orange), and twice the ground pressure (green). T and IWVLOS are fixed to the same values for the three models. Bottom: master modeled telluric CCF at IWVLOS (blue), half of IWVLOS (orange), and twice IWVLOS (green). T and Pmean are fixed to the same values for the three models.

4 Application of the telluric correction for different spectral types

We first performed a test on the spectrophotometric standard star, HR 1544. As it is an A-type star, contamination by stellar lines is negligible as only a few broad lines are present, while telluric lines are narrower. We then tested the telluric correction for other spectral types up to late-type stars to see if the presence of stellar lines biases the telluric correction.

4.1 HR 1544

HR 1544 is an A-type star and is part of the spectrophotometric standard catalog of ESO for ESPRESSO. We selected a spectrum at a high S/N (~350) with an integrated water vapor of 3.56 mm measured at the zenith by one of the radiometers installed on Paranal (Kerber et al. 2012), which is above the median value of 2 mm⁵. Figure 7 shows this spectrum around the subset of lines selected to perform the optimization of the parameters before and after applying the telluric correction, as well as the computed telluric spectrum. The selected lines acquired under this weather condition have an absorption depth ranging from 50 to 80%. Some systematic residuals are present after the telluric correction around the deepest lines but with a decreasing residual amplitude for the weakest telluric lines. The lines shown are the strongest in the ESPRESSO spectral band, except for the lines at 7186.51, 7208.41, and 7236.72 Å. Figure 8 displays them alongside their strong residuals. We note that their residuals have different structures from one line to another and that it was not possible to properly model them, leading to their exclusion. The lack of proper modeling due to (quasi-)saturation can be explained by several factors. First, the one-layer model described here is too simple for strong telluric lines. For example, asymmetric broadening could introduce P-Cygni-like profiles. Secondly, some physics of the lines is missing to reproduce their core. It could come from the inclusion of wind that shifts the line core, or some HITRAN parameters that are inaccurate. Finally, it could be a systematic error caused by the noncommutativity of multiplication and convolution, which can be significant for deep telluric lines.

We used the subset of selected lines shown in Fig. 7 to create the master observed telluric CCF, the master modeled telluric CCF, and the master observed telluric-corrected CCF shown in Fig 9. They represent the average line profile in the 80 km s⁻¹ surrounding each of the selected lines. One can see that the Levenberg–Marquardt minimization converges to residuals below the 2% P2V amplitude on the average line. The latest value is in agreement with typical telluric correction from the literature (e.g., Sameshima et al. 2018; Ulmer-Moll et al. 2019; Baker et al. 2020).

4.2 From F-type stars to K-type stars

Figure 9 also shows the observed telluric CCF, the master modeled telluric CCF, and the master observed telluric-corrected CCF for WASP-76, WASP-127, and GJ 9827, which are F-, G-, and K-type stars, respectively. The three stars are part of the observations of the ESPRESSO consortium with an individual S/N of ~ 130 and with an IWV of 0.7, 7, and 5 mm respectively. It is noticeable that the IWV not only plays an important role in the absorption line, as expected, but it also leads to stronger residuals in the very core of telluric lines. Nonetheless, the residuals are within or close to the 2% P2V amplitude and in agreement with the continuum dispersion. The later type stars have more and stronger stellar lines that slightly affect the CCF continuum, such as at ~30 km s⁻¹ for WASP-127 or ~25 km s⁻¹ for GJ 9827, but do not prevent the recovery of the water line profile and the correction of telluric contamination.

4.3 The case of M-type stars

M-type stars are rich in molecular bands in their visible spectrum, causing absorption of the stellar continuum to an extent that what we measure is de facto a pseudo-continuum. It is thus challenging to disentangle stellar and telluric lines. Figure 9 also includes the observed telluric CCF, the master modeled telluric CCF, and the master observed telluric-corrected CCF for GJ 436 and Proxima, which are M2 and M5-type stars, respectively. It is not possible to discriminate between a telluric residual and the dispersion of the pseudo-continuum in the master observed telluric-corrected CCF. The spectrum of GJ 436 was obtained in extreme conditions with an IWV of 12 mm, leading to an absorption depth of about 90%, and thus a low S/N in the telluric line cores. On the contrary, the spectrum of Proxima was obtained under dry conditions with an IWV of 0.3 mm, leading to a telluric correction at the level of the pseudo-continuum dispersion.

From these examples, it is possible to conclude that the telluric model presented here can be applied automatically to all observations of ESPRESSO⁶, and it is not biased by the stellar spectrum and contamination between stellar and telluric lines. Nonetheless, some telluric residuals can be present in the core of the strongest lines but over a very narrow spectral range in comparison to the whole spectral domain impacted by water vapor telluric lines.

Fig. 7 Stellar spectrum of HR 1544 before (black) and after (gray) telluric correction around the subset of the 20 selected lines. Each panel is for a different line, with its position indicated by the vertical dashed blue line. The telluric spectrum is shown as the dashed red curve.

Fig. 8 Same legend as for Fig. 7 but for the three lines excluded from the subset.

5 Application of the telluric correction for radial velocity

To assess the impact of the telluric correction on the RVs obtained with ESPRESSO, we corrected the spectra of Tau Ceti and Proxima from telluric contamination by dividing each spectrum by its best-fit telluric template. Then, the spectra were cross-correlated with the DRS masks as the data without telluric correction. The next step was to see if a broader spectral coverage could be used.

It is essential to understand how ESPRESSO stellar masks are built. Each mask corresponds to a given spectral type (or sub spectral type for M dwarfs) but they are not specific to a given star. The masks are built from high S/N spectra of inactive stars observed with ESPRESSO. Firstly, the line locations are identified by searching for local minima. Then, spectral regions affected by magnetic activity stellar lines (such as Ca ii H & K, Na, and H-α) and by strong telluric bands (5871–6005 Å, 6271–6341 Å, 6449–6605 Å, 6861–7416 Å, 7588–7751 Å, and 7860–8000 Ä) are excluded. As the masks need to be the same, independently of the star or epoch of the year, margins of 64km s⁻¹ (1.5 Å at 7010 Å) for active stellar lines and 214km s⁻¹ (5 Å at 7010 Å) for telluric bands, to include systemic velocity and barycentric shifts, are excluded.

Therefore, to understand how the telluric correction could help to reduce these excluded regions, we produced new masks based on telluric-free spectra for Tau Ceti and Proxima. Then, we modified the exclusion regions as defined by the mask in different steps to understand how the new spectral regions modify the RV time series. The different RV time series produced for each star are named DRS, Microtell, Blue, Orders77_78, Orders77, 15_Slines, 1_Sline, and Ideal. The associated masks aer shown in Figs. C.1 and C.2.

DRS is the standard RV time serie with no telluric correction and no modification of the mask. It is important to note that this is not necessarily the best mask as, for example, micro-telluric lines can have an impact on the shape of the CCF.

Microtell is the RV time serie obtained with the data corrected for the telluric lines but using the same mask as DRS. These data show the impact of micro-telluric lines on the shape of the CCFs, and therefore on the extracted RVs.

Blue is the RV time serie obtained with the data corrected for the telluric lines but the mask includes the wavelength regions with the telluric water bands around the Na and H-α lines, and the reddest water band.

Orders77_78 is the RV time serie obtained with the data corrected for the telluric lines but the mask includes the wavelength regions with all the telluric water bands. It is important to note, however, that both orders 77 and 78 (slices 154, 155, 156, and 157), where water telluric lines are the strongest, and regions containing O₂ lines are still excluded, as in the standard mask.

Orders77 is the RV time serie obtained with the data corrected for the telluric lines and the mask includes the wavelength regions with all the telluric water bands, but the physical order 77 (slices 154, 155) is still excluded in addition to the oxygen bands.

15_Slines is the RV time serie obtained with the data corrected for the telluric lines and the mask includes the wavelength regions with all the telluric water bands except the 15 strongest water lines with an exclusion window of 214 km s⁻¹.

1_Sline is the RV time serie obtained with the data corrected for the telluric lines and the mask includes the wavelength regions with all the telluric water bands except the strongest water line at 7186.51 A with an exclusion window of 214 km s⁻¹ around each line.

Ideal is the RV time serie obtained with the data corrected for the telluric lines and the mask includes the wavelength regions with all the telluric water bands. It corresponds to the mask that should be used if there were no telluric water lines or if their correction was assumed to be perfect (i.e., affecting the photon noise without introducing systematic in the RVs).

We also note that ESPRESSO had several technical interventions introducing possible offsets at 2458662 and 2459202 BJD (Faria et al. 2022).

Fig. 9 Master observed telluric CCF (in black), the master modeled telluric CCF (in red, dashed) and the master observed telluric-corrected CCF (in gray) for HR 1544, WASP-76, WASP-127, GJ 9827, GJ 436, and Proxima. The CCFs are performed over the subset of selected water lines. The two horizontal blue dotted lines encompass the 2% residual P2V.

5.1 Tau Ceti

Tau Ceti (V-mag = 3.5) is one of the most studied stars and has been heavily observed by high-precision RV campaigns with HARPS and now ESPRESSO (e.g., Dumusque et al. 2011; Feng et al. 2017). Between July 2018 and December 2022, 1982 data points across 202 nights were acquired with ESPRESSO. The observational strategy for Tau Ceti consists in taking several very short exposures per night. Therefore, the results presented below are nightly binned. Tau Ceti RVs are known to be highly stable with no clear sign of activity nor any massive planets (more massive than super-Earths). It is thus a perfect target for studying the impact of the telluric correction.

Figure 10 shows the best-fit parameters of IWV and pressure with the corresponding chi-square for Tau Ceti. Our computed IWV values are in agreement with the radiometer values, and the slight differences (rms of the difference is ~0.4 mm) can be explained by different water column distributions between the line of sight and the zenith. The derived pressure varies from half of the ground pressure up to the ground pressure⁷, with seasonal variations. Finally, we note that the quality of the fits is stable over time, confirming the capability of our technique to be applied automatically. Figure 11 shows the difference between the RV time series obtained with the different masks after the telluric correction is applied, and the RV time series obtained with the DRS. One can see that the Microtell RV time series has some systematic differences in comparison to the DRS RV time series. These variations are due to the presence of microtelluric lines that induce RV variations of about 58 cm s⁻¹ from P2V with an rms of 10 cm s⁻¹. When phased as a function of the BERV, the residuals also show a periodic contribution (Fig. 12) as expected by the Earth’s orbital motion around the Sun. The same structures are also present with the other masks. We confirmed the micro-telluric contribution by computing the RVs for the DRS and Microtell RV time series in the spectral range from 4950–5855 Å, where the micro-telluric lines are located (see Fig. 1). We observed the exact same systematics but with an amplified signal up to 1.4 ms⁻¹. It is thus proof that Tau Ceti spectra are impacted by micro-telluric contamination. Moreover, the Order77, 15_Slines, l_Sline, and Ideal RV time series have increasing systematics when a broader spectral range is included. They are probably due to telluric residuals around the 7250Å water band. The large difference between the l_Sline and Ideal RV time series is caused by a single telluric line residual (at 7186.51 Å). This shows how important it is to exclude telluric lines or to be very confident in the level of correction achievable. We also note an RV offset for the last epoch of observations after 2459202 BJD, which coincides with a technical intervention as mentioned earlier. To conclude, the mask Orders77_78 is selected for G-type stars, such as Tau Ceti, as it is a good trade-off between reducing bias from telluric residuals and the gain in photon noise (~2%, Table 1). The slight gain in photon noise is caused by the blackbody maximum position at blue wavelengths for G-type stars, where few spectral regions are excluded.

Fig. 10 Properties of the telluric correction best-fit. First panel: barycentric correction as a function of time for Tau Ceti observations. Second panel: IWV as a function of time. In blue is the airmass scaled IWVlos as measured by the radiometer and in orange is the best-fit IWV value. Third-panel: pressure as a function of time. In blue is the measured ground pressure (it is the upper limit to the fitted pressure) and in orange is the best-fit pressure value. Last panel: best chi-square values of the fit of the automatic telluric correction.

Fig. 11 Tau Ceti RV time series differences between DRS and, from top to bottom: Microtell, Blue, Orders77_78, Order77, 15_Slines, 1_Sline, and Ideal.

Fig. 12 Tau Ceti RV time series differences between DRS and Microtell as a function of BERV and colored by the IWV content.

5.2 Proxima

Proxima is also one of the most studied stars and has been heavily observed by high-precision RV campaigns with HARPS and ESPRESSO. Proxima is the closest star to the Sun and is orbited by an Earth-like planet at a period of 11 days (Anglada-Escudé et al. 2016), and by a planet candidate at a large orbital distance (Damasso et al. 2020). Between February 2019 and August 2021, 117 data points were acquired and led to the independent confirmation of Proxima b by ESPRESSO (Suárez Mascareño et al. 2020), and the detection of Proxima d at an orbital period of 5.12 days (Faria et al. 2022). As the star is an M5 dwarf, its spectral energy distribution peaks toward red wavelengths, where the strongest telluric lines are located. Therefore, it is a perfect target for verifying if the telluric correction can lead to improved precision by including more stellar lines in the CCFs.

Figure 13 shows the best-fit parameters of IWV and pressure with the corresponding chi-square for Proxima. We confirm the same behaviors for the IWV (rms of the difference is ~0.7 mm), the pressure, and the fit goodness as for Tau Ceti. Figure 14 shows the difference between the RV time series obtained with the different masks after the telluric correction is applied, and the RV time series obtained with the DRS. One can see that the Microtell and Blue RV time series are similar to the DRS RV time series. This is expected as the micro-telluric lines and the spectral range gained with the Blue mask are in the blue part of the ESPRESSO range, where there is very little flux for a star such as Proxima. The Orders77_78 RV time series differs from the DRS RV time series by an offset for the RVs obtained after 2458650 BJD but no systematics are introduced. This RV offset coincides with a technical intervention on ESPRESSO. On the contrary, Order77, 15_Slines, l_Sline, and Ideal RV time series have clear systematics compared to the DRS RV time series, which increase when a greater spectral range is included. Even if the DRS time series is not the reference, introducing systematics up to 4 m s⁻¹ is not acceptable. Therefore, the telluric correction is not good enough in the core of the 7200 Å water band for it to be used for precise RV studies. Nonetheless, from Table 2, the gain in photon noise between the time series DRS and Orders77_78 is about 25%, while the gain between the time series Orders77_78 and Ideal is negligible. This factor of 1.33 in photon noise between 28 and 21 cm s⁻¹ corresponds to a factor of 1.78 in exposure time. In other words, using the Orders77_78 RVs is equivalent to increasing the exposure time of the DRS RVs by 1.78 to reach a photon noise of 21 cm s⁻¹. On the contrary, reaching the same photon noise as the DRS RVs (28 cm s⁻¹) is equivalent to a gain of 44% in exposure time on the Orders77_78 RVs. To conclude on the mask selection, the best mask to use for Proxima is Orders77_78.

Fig. 13 Properties of the telluric correction best-fit. First panel: barycentric correction as a function of time for Proxima observations. Second panel: IWV as a function of time. In blue is the airmass scaled IWVlos as measured by the radiometer and in orange is the best-fit IWV value. Third panel: pressure as a function of time. In blue is the ground pressure (it is the upper limit to the fitted pressure) and in orange is the best-fit pressure value. Last panel: best chi-square values of the fit of the automatic telluric correction.

Fig. 14 Proxima RV time series differences between DRS and, from top to bottom: Microteil, Blue, Orders77_78, Order77, 15_Slines, l_Sline, and Ideal.

5.3 Impact on the parameters of Proxima b and d

To better understand the impact of the telluric correction on the RV delivered by the ESPRESSO DRS, we measured the impact of the telluric correction on the CCFs for the recovered orbital parameters of Proxima b and d, considering the same model used in Faria et al. (2022). We compared the standard RVs produced by the DRS (hereafter referred to as DRS RVs) to the RVs computed with the Orders77_78 mask on the telluric-corrected spectra (hereafter referred to as telluric-corrected RVs). This analysis simultaneously models the RV and FWHM of the CCFs, to constrain the component of the RV signal that is caused by stellar activity. We consider a Gaussian process (GP) with shared hyperparameters between the RVs and FWHM, with the exception of the amplitudes. In the RVs, the model further assumes two Keplerian signals. Following Faria et al. (2022), we subdivide the ESPRESSO data set into three subsets, ESP18, ESP19, and ESP21, and we consider possible RV offsets between these subsets. All priors are the same as those considered in Faria et al. (2022).

The two signals with periods of 11.19 and 5.12 days are recovered with very high significance in the telluric-corrected RVs. When compared to the DRS RVs (Fig. 15), the constraints on the orbital parameters of Proxima b are similar, but the semi-amplitude and eccentricity of Proxima d are found to be smaller, although within 1 sigma (Table 3). This decrease goes in the same direction as the result obtained with the template-matching RVs analyzed by Faria et al. (2022), which uses the real RV content of each line but masks the telluric lines. We do not perform a direct comparison between the telluric-corrected CCF RVs and the template-matching RVs as the gain in photon noise and increase in parameters constraints arise from two different aspects⁸.

Considering the maximum-likelihood solution, the rms of the RV residuals reaches 34 cm s⁻¹ in the telluric-corrected RVs, compared to about 40 cm s⁻¹ in the DRS RVs. This improvement is also visible in each individual subset of the ESPRESSO data (Table 4).

Another comparison, which tries to include the full information in the posterior distributions, uses the so-called information, a measure of how much the prior is compressed to reach the posterior. Since the priors are the same in the analysis of the telluric-corrected RVs and DRS RVs, a larger compression means that the data provide more constraints on the model parameters. We obtain a value for the information of I = 61 nats⁹ for the telluric-corrected RVs and I = 59 nats for the DRS RVs. Although the values are close, they again reflect the improvement coming from using the telluric-corrected RVs.

6 Application to the O₂ molecular bands

The optical part of the telluric spectrum, in addition to the water bands discussed above, is dominated by four oxygen bands: δ (~579 nm), γ (~629 nm), B (~688 nm) and A (~765 nm). From Fig. 1, one can see that the B- and A-bands are dominated by saturated lines. There is no gain in correcting saturated or nearly saturated lines, which have an absorption depth above 90%. Therefore, we consider continuing to mask these two bands for RV purposes but applying our telluric model to the γ-band, where the gain in spectral coverage, and thus RV content, is nonzero.

The physical parameters for ¹⁶O₂ are from HITRAN (Rothman et al. 2013; Gordon et al. 2017). We found that a Lorentzian profile was not sufficient to accurately model the oxygen lines, whereas a Voigt profile with a free thermal broadening for the line core was more efficient. As opposed to water vapor, oxygen is present higher up in the Earth’s atmosphere, where the pressure broadening becomes negligible. We therefore added the temperature to the integrated oxygen column density and pressure, as a free parameter in the minimization fit. Equation (5) can thus be modified as follows: $$T\left(v\right)={\text{e}}^{-{\text{IOC}}_{\text{LOC}}\cdot n{o}_2\cdot \sum \sigma \left(v\right)},$$(17)

where $n_{{\text{O}}_2}=5.65\times {10}^{18}$ molecules cm⁻³ and IOC is the integrated oxygen column density at the zenith, and is about 7 km. On the contrary to the IWV, the IOC_LOS varies mainly as a function of the airmass. This correlation explains the strong dependence of the RVs that are measured on O₂ lines on the airmass measured with HARPS (Figueira et al. 2010).

Figure 16 shows the master observed telluric CCF, the master modeled telluric CCF, and the master observed telluric-corrected CCF, all obtained with the subset of O₂ lines in the γ-band for HR 1544 and Proxima. In both cases, telluric residuals are below the 2% P2V amplitude or within the continuum dispersion as in the case of water lines. This proves that the telluric model described here can be applied to different molecules.

Fig. 15 Posterior distributions for the orbital period, semi-amplitude, and eccentricity of the 11-day (top) and 5-day (bottom) signals. The dashed blue histograms correspond to the analysis of the DRS RVs, while the red histograms correspond to the analysis of the telluric-corrected RVs (TEL).

7 Conclusion

The telluric correction presented here is based on a simple line-by-line telluric model, which requires the physical parameters of the lines (e.g., from the HITRAN database Rothman et al. 2013), the sky conditions, and a detailed map of the instrument’s spectral resolution as a function of the position on the detector. Then, the model is adjusted to the observations through a Levenberg–Marquardt minimization by building an average line profile over a subset of selected lines with strong absorptions.

While this model can be applied to any molecules, our focus in this paper is the study of H₂O lines in the spectral range of ESPRESSO (3800–7800 Å). The water vapor lines are best represented by Lorentzian profiles, which are mainly dominated by the integrated water vapor that controls the line strength, but also by the pressure that controls the line width. We demonstrated that the telluric correction works for all stellar spectral types with residuals below a 2% P2V amplitude. Similar performances have been obtained with other methods (e.g., Artigau et al. 2014; Smette et al. 2015; Ulmer-Moll et al. 2019; Leet et al. 2019; Bedell et al. 2019; Cretignier et al. 2021), but our model does not require a telluric database, can be directly applied to a single observation, and is simple enough to be implemented in the DRS of state-of-the-art spectrographs to deliver high-fidelity spectra and RVs.

Telluric lines with absorption depths greater than 70–80% exhibit residual structures with P-Cygni-like profiles. Several reasons can explain these structures, such as a lack of sophistication in our one-layer model (even if more complex models show similar results (Smette et al. 2015; Allart et al. 2020)), inexact parameters in the HITRAN database, the presence of additional phenomena that deformed the line shape (e.g. winds), or systematic errors due to the noncommutative arithmetic between convolution and division. We note that Earth scientists use more complex line shapes to model the telluric lines but the level of complexity and number of parameters to take into account are incompatible with the simplicity that we aimed for here.

The main goal of this telluric correction is to provide more precise RVs, unbiased by micro-telluric lines (Cunha et al. 2013). The example of Tau Ceti showcased that micro-telluric lines can introduce systematics up to 58 cm s⁻¹ P2V with an rms of 10 cm s^−l. It is therefore critical to properly correct them when searching for Earth-mass planets. We also compared the RVs produced by the standard masks of the ESPRESSO DRS with different custom masks, testing the impact of the inclusion of a larger spectral region. To stick to the original idea of the ESPRESSO DRS, we kept these masks independent of the BERV and systemic velocities of the stars. Our best custom mask, Orders77_78, while excluding orders 77 and 78 (where a couple of telluric water lines still have strong residuals), enables a gain of ~500 Å (only ~ 14% of the spectral range is then excluded, compared to ~27% with the DRS mask) without introducing systematics. In the case of Proxima, this additional wavelength coverage leads to a gain of 25% in photon noise or a factor of 1.78 in exposure time. This can lead to the detection of fainter RV signals but also a more efficient use of telescope time.

We applied the same analysis as Faria et al. (2022) to the telluric-corrected CCF RVs of Proxima and we confirmed the presence of the planet Proxima b and planet candidate Proxima d. While the parameters of Proxima b are similar in both analyses, the semi-amplitude and eccentricity of Proxima d are smaller and their distributions go in the direction of the template matching analysis. Although the improvements brought by those two techniques converge to similar results for the parameters of planet candidate d, we want to note that those two methods do optimize completely different things. On one side, template matching uses the real RV content of each line and the most RV information possible, as in M dwarfs there is RV information everywhere, even in the pseudo-continuum. On the other side, the method proposed in this paper corrects for telluric absorption, and therefore a gain in RV information that leads to a better estimation of the planetary signal. We also want to note that improvements in the definition of the weights of each stellar line for the CCF technique are underway (Bourrier et al. 2021).

We extended the application of our correction to the O₂ γ-band to display that our model can correct the telluric lines of molecules other than H₂O. In the case of the γ-band, the telluric residuals are within the requirements to perform extremely precise RV measurements, even if the spectral range freed remains modest (~80 Å).

Adapting the model to other instruments, such as HARPS and HARPS-N instruments, will lead to a substantial increase in RV content. A gain in photon noise of 10% could be reached for M5-type stars, and a ~60 cm s⁻¹ P2V micro-telluric amplitude could be removed. Furthermore, near-infrared high-resolution spectrographs, such as NIRPS, are the most important spectographs for applying a well-functioning telluric correction. At these wavelengths, H₂O is not the only main absorber, and CO₂, O₂, and CH₄ absorb on large spectral ranges, and with deep lines. It is, therefore, critical to optimize the telluric correction to facilitate the maximum extraction of RV information.

The impact of the telluric correction for all the stars observed by the ESPRESSO consortium still needs to be investigated in detail but we expect that it will lead to more precise and less biased discoveries. The impact on other science cases, such as stellar activity indicators, atmospheric characterization, or the measurement of the fine structure constant, still has to be studied. Removing micro-telluric lines is potentially very important for measurements of the fine structure constant with ESPRESSO: micro-telluric lines which fall on top of the relevant transitions in the quasar spectrum could spoil the measurement if not removed or recognized. Therefore, we can be optimistic regarding the benefits of applying our correction methods to different science cases.

Fig. 16 Master observed telluric CCF (in black), master modeled telluric CCF (in red, dashed) and master observed telluric-corrected CCF (in gray) for HR 1544 (top) and Proxima (bottom). The CCFs are performed over the subset of selected 02 lines in the γ-band. The two horizontal blue dotted lines encompass the 2% residual P2V.

Acknowledgements

We thank the referee for their useful comments that improve the overall paper quality. R.A. thanks Etienne Artigau for useful discussions on the behavior of the telluric lines, and the comparison with his telluric correction for the SPIRou instrument. R.A. thanks Alain Smette and Martin Turbet for useful discussions on the telluric line behavior and line shape. R.A. is a Trottier Postdoctoral Fellow and acknowledges support from the Trottier Family Foundation. This work was supported in part through a grant from FRQNT. 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. The authors acknowledge the financial support of the SNSF. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement SCORE No. 851555). A.M.S. acknowledges support from the Fundação para a Ciência e a Tecnologia (FCT) through the Fellowship 2020.05387.BD. and POCH/FSE (EC). 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 Inter-nacionalização by these grants: UID/FIS/04434/2019; UIDB/04434/2020; UIDP/04434/2020; PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113; PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953; PTDC/FIS-AST/28987/2017 & POCI-01-0145-FEDER-028987. J.I.G.H., R.R., and A.S.M. acknowledge financial support from the Spanish Ministry of Science and Innovation (MICIN) project PID2020-117493GB-I00. J.I.G.H. and R.R. also acknowledge financial support from the Government of the Canary Islands project ProID2020010129. A.S.M. acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) under 2018 Juan de la Cierva program IJC2018-035229-I. N.J.N. acknowledges support from the following projects: UIDB/04434/2020 & UIDP/04434/2020, CERN/FIS-PAR/0037/2019, PTDC/FIS-OUT/29048/2017, COMPETE2020: POCI-01-0145-FEDER-028987 & FCT: PTDC/FIS-AST/28987/2017, PTDC/FIS-AST/0054/2022. D.M. is also supported by the INFN PD51 INDARK grant.

Appendix A Coefficient of the resolution parametrization for the HR1x1 mode

Appendix B List of the 20 water telluric lines used in the subset and the telluric CCF computation

Appendix C Masks used

Fig. C.1 Telluric spectrum in the range of ESPRESSO. In blue, the H2O lines, and in red the O2 lines. The different masks are shown in the panels and are ordered as a function of the spectral coverage included in the mask. In gray are the regions excluded from the mask.

Fig. C.2 Inset of Fig. C.1 between 6860 and 7900 Å.

References

All Tables

All Figures

Fig. 1 Telluric and stellar spectra over the ESPRESSO spectral range. The gray bands are the regions excluded from the mask used in the ESPRESSO DRS. Top panel: telluric spectrum produced with Molecfit (Smette et al. 2015; Kausch et al. 2015). The H2O lines are in blue, while the O2 lines are in red. Middle panel: ESPRESSO spectrum of the G-type star Tau Ceti. Bottom panel: ESPRESSO spectrum of the M5-type star Proxima.

In the text

	Fig. 2 Reversed cumulative distribution of the water vapor line intensities in the water band ranging from 3800 to 8000 A.
In the text

Fig. 3 Impact of the number of lines selected in the telluric model. Top panel: average water telluric line built for the model, with line intensity thresholds at 10−24, 10−25, 10−26, 10−27, 10−28, and 10−29 cm−1/(molecule cm−2). Middle panel: residuals of the average lines with respect to including all the lines. See Sect. 3 for more details on how the average line is built. Bottom panel: zoom on the residuals obtained with the model that includes the 10−26, 10−27, 10−28 cm−1/(molecule cm−2) strongest water lines. The dispersion and P2V amplitude of the residuals are reported in the legend.

In the text

	Fig. 4 Resolution map (orders versus pixels) for blue and red detectors in the HR1 × 1 mode. The dots are the measured resolution on individual thorium lines at known positions in pixels and in orders. The background map is the polynomial fit applied to the data.
In the text

Fig. 5 Impact of the line profile to best model the data. Top panel: CCFs computed over the 20 strongest water lines. The blue CCF is the master observed telluric CCF based on ESPRESSO data, while the orange, green, and red CCFs are the master modeled telluric CCFs assuming Gaussian, Lorentzian, or Voigt profiles, respectively. Bottom panel: residuals only for the Lorentzian and Voigt profiles. The Gaussian profile is not shown due to its larger residuals.

In the text

Fig. 6 Impact of T (top), P (middle), and IWVLOS (bottom) on the telluric line shape. Top: master modeled telluric CCF at an atmospheric temperature of 282 K (blue), 293 K (orange) and 273 K (green). Pmean and IWV are fixed to the same values for the three models. Middle: master modeled telluric CCF at ground pressure (blue), half of the ground pressure (orange), and twice the ground pressure (green). T and IWVLOS are fixed to the same values for the three models. Bottom: master modeled telluric CCF at IWVLOS (blue), half of IWVLOS (orange), and twice IWVLOS (green). T and Pmean are fixed to the same values for the three models.

In the text

	Fig. 7 Stellar spectrum of HR 1544 before (black) and after (gray) telluric correction around the subset of the 20 selected lines. Each panel is for a different line, with its position indicated by the vertical dashed blue line. The telluric spectrum is shown as the dashed red curve.
In the text

	Fig. 8 Same legend as for Fig. 7 but for the three lines excluded from the subset.
In the text

	Fig. 9 Master observed telluric CCF (in black), the master modeled telluric CCF (in red, dashed) and the master observed telluric-corrected CCF (in gray) for HR 1544, WASP-76, WASP-127, GJ 9827, GJ 436, and Proxima. The CCFs are performed over the subset of selected water lines. The two horizontal blue dotted lines encompass the 2% residual P2V.
In the text

Fig. 10 Properties of the telluric correction best-fit. First panel: barycentric correction as a function of time for Tau Ceti observations. Second panel: IWV as a function of time. In blue is the airmass scaled IWVlos as measured by the radiometer and in orange is the best-fit IWV value. Third-panel: pressure as a function of time. In blue is the measured ground pressure (it is the upper limit to the fitted pressure) and in orange is the best-fit pressure value. Last panel: best chi-square values of the fit of the automatic telluric correction.

In the text

	Fig. 11 Tau Ceti RV time series differences between DRS and, from top to bottom: Microtell, Blue, Orders77_78, Order77, 15_Slines, 1_Sline, and Ideal.
In the text

	Fig. 12 Tau Ceti RV time series differences between DRS and Microtell as a function of BERV and colored by the IWV content.
In the text

Fig. 13 Properties of the telluric correction best-fit. First panel: barycentric correction as a function of time for Proxima observations. Second panel: IWV as a function of time. In blue is the airmass scaled IWVlos as measured by the radiometer and in orange is the best-fit IWV value. Third panel: pressure as a function of time. In blue is the ground pressure (it is the upper limit to the fitted pressure) and in orange is the best-fit pressure value. Last panel: best chi-square values of the fit of the automatic telluric correction.

In the text

	Fig. 14 Proxima RV time series differences between DRS and, from top to bottom: Microteil, Blue, Orders77_78, Order77, 15_Slines, l_Sline, and Ideal.
In the text

	Fig. 15 Posterior distributions for the orbital period, semi-amplitude, and eccentricity of the 11-day (top) and 5-day (bottom) signals. The dashed blue histograms correspond to the analysis of the DRS RVs, while the red histograms correspond to the analysis of the telluric-corrected RVs (TEL).
In the text

	Fig. 16 Master observed telluric CCF (in black), master modeled telluric CCF (in red, dashed) and master observed telluric-corrected CCF (in gray) for HR 1544 (top) and Proxima (bottom). The CCFs are performed over the subset of selected 02 lines in the γ-band. The two horizontal blue dotted lines encompass the 2% residual P2V.
In the text

	Fig. C.1 Telluric spectrum in the range of ESPRESSO. In blue, the H2O lines, and in red the O2 lines. The different masks are shown in the panels and are ordered as a function of the spectral coverage included in the mask. In gray are the regions excluded from the mask.
In the text

	Fig. C.2 Inset of Fig. C.1 between 6860 and 7900 Å.
In the text