© The Authors 2024

1 Introduction

The last two decades have seen a stupendous rise in the development of ground-based spectrographs, fuelled by their application in both the planetary and stellar communities (Pepe et al. 2014). The efficiency, spectral resolution, and wavelength coverage of available instrumentation have increased to the point where several spectrographs across the world now allow for observations in the optical and near-infrared (NIR) domains with resolving powers up to more than 100000 (e.g., ESPRESSO, Pepe et al. 2021; NIRPS, Bouchy et al. 2017). Another key point that qualifies these spectrographs to characterize exoplanetary systems (besides their high spectral resolution) is that they are fiber-fed and extremely stabilised in temperature and pressure, which ensures a high stability for their point-spread function. This allows for the lines from species in stellar and planetary atmospheres to be resolved with a high precision, thus determining their chemical composition (by identifying which transition the line corresponds to), physical properties (by analyzing the line profile shape), and dynamical motions (by analyzing the Doppler shifts of the lines with respect to their expected rest position). While spectrographs in the exoplanet field were originally designed to search for planets through the Keplerian motion of their star (Mayor & Queloz 1995; Baranne et al. 1996), transiting planets further unlocked their potential. Spectra measured during transit do indeed mix information from the unocculted photosphere, from the local photosphere fully occulted by the opaque planetary layers, and from the planet atmospheric limb filtering the local photospheric light.

The first exploitation of an exoplanet spectroscopic transit (Queloz et al. 2000) repurposed a velocimetric technique originally applied to eclipsing binaries. As a transiting body occults a local region of the photosphere, it removes its light from the observed stellar spectrum. This distortion of the disk-integrated stellar lines, known as the Rossiter-McLaughlin (RM) effect (Rossiter 1924; McLaughlin 1924), traces the trajectory of the body across the stellar disk, thus constraining the velocity field of the photosphere and the orbital architecture of the system. The latter is an important tracer for exoplanets, as their spin-orbit angle is inherited from the protoplanetary disk or shaped by late dynamical migration (see Albrecht et al. 2022 and Bourrier et al. 2023 for reviews of these processes). Measurements of the RM effect evolved over the years, from analyzing the anomalous radial velocity (hearafter, rv) deviation from the Keplerian motion of the star (e.g., Ohta et al. 2005; Giménez 2006; Hirano et al. 2011; Boué et al. 2013) to modeling the distorted disk-integrated stellar lines using Doppler tomography (e.g., Collier Cameron et al. 2010; Bourrier et al. 2015; Temple et al. 2019), and, finally, extracting the planet-occulted stellar lines to analyze their rv centroids via the reloaded RM effect (Cegla et al. 2016) or model their profiles via the RM effect Revolutions (Bourrier et al. 2021). The possibility to analyze the local stellar lines occulted by a transiting planet translates, effectively, into a spatial mapping of the photosphere (e.g., Dravins et al. 2017, 2021).

At wavelengths where the atmospheric components of an exoplanet absorb the stellar light, the transit depth appears larger than in white light (Seager & Sasselov 2000). This is the conceptual idea that Charbonneau et al. (2002) used to derive the transmission spectrum of an exoplanet atmosphere through its spectroscopic transit. While the medium spectral resolution (R ~ 5500) of the Hubble Space Telescope Imaging Spectrograph (HST/STIS) did not allow them to resolve the narrow planetary absorption lines, they measured a larger transit depth in the region of the sodium doublet that probably arises from the line wings (e.g., Carteret et al. 2024). Transit spectroscopy then expanded from space-borne, medium-resolution (e.g., Vidal-Madjar et al. 2003; Sing et al. 2008, 2015; Deming et al. 2013; Stevenson et al. 2014) to ground-based, high-resolution (Redfield et al. 2008; Snellen et al. 2008). While absolute transit depth is lost from the ground, the higher spectral resolution (R ≳ 60 000–100 000) allows for the planetary atmospheric lines to be resolved, while avoiding their blurring by the planet orbital motion, so that their core can be detected through comparisons with the surrounding continuum. High-resolution spectroscopy with stabilized spectrographs gives access to atmospheric signatures from both non-transiting and transiting planets, since the ability to perform precise rv measurements allows us to track the Doppler shift of individual spectral lines along the orbital motion. Molecular bands imprinted in the bright dayside spectrum of the planet, for example from H₂O and CO, could thus be detected in the NIR (Brogi et al. 2012; Birkby et al. 2013; Lockwood et al. 2014; Snellen et al. 2014), while iron emission lines were recently detected in the optical spectrum of ultra-hot gas giants (Pino et al. 2020; Nugroho et al. 2020; Yan et al. 2020). During transit, measurements of absorption lines blueshifted with respect to the planet orbital motion – again, from molecules in the NIR first (Snellen et al. 2010) and then from atoms in the optical (Wyttenbach et al. 2015; Louden & Wheatley 2015) were attributed to winds in the planetary atmosphere. Spectrally resolving the planetary lines provides further constraints on the thermal and dynamical broadening mechanisms in the atmosphere (Seidel et al. 2019, 2020; Wyttenbach et al. 2020). Combined with the possibility to resolve absorption lines temporally during transit, high-resolution spectroscopy thus allows mapping spatially and dynamically the planetary atmosphere (Louden & Wheatley 2015; Allart et al. 2018; Ehrenreich et al. 2020; Gandhi et al. 2022; Seidel et al. 2023).

Characterizing exoplanetary atmospheres from the ground, which was long thought to be unfeasible due to time-variable telluric contamination, became possible thanks to high spectral resolution. The fast orbital motion of close-in planets allows us to disentangle in rv space the fast-moving planetary lines from telluric lines and, outside of the transit, from the unocculted stellar lines. Analyses remain complicated during transit, when disk-integrated stellar lines may be distorted by both planet-occulted stellar lines and planetary atmospheric lines. On the one hand, atoms in the atmosphere of the hottest gas giants may absorb the same transitions as the stellar atmosphere and contaminate stellar lines traditionally used for RM analysis (e.g., Bourrier et al. 2020; Ehrenreich et al. 2020). More critically, transmission spectra of planetary atmospheres are usually corrected for the stellar lines they filter using the out-of-transit stellar spectrum. Because the disk-integrated spectrum is often not a good estimate of local stellar lines (e.g., Czesla et al. 2015; Yan et al. 2017; Casasayas-Barris et al. 2020), planet-occulted line distortions (POLDs) are introduced at every transition absorbed in the stellar spectrum (Dethier & Bourrier 2023). The rv tracks of the planetary orbital motion and of the occulted stellar surface generally overlap over a substantial fraction of the transit (e.g., Casasayas-Barris et al. 2022), so that planetary lines are contaminated by POLDs in all but the hotter stars, which have few absorption lines. Various approaches have been taken to correct transmission spectra for POLDs a posteriori (Yan et al. 2017; Casasayas-Barris et al. 2020, 2021; Borsa & Zannoni 2018; Borsa et al. 2021; Mounzer et al. 2022), or to model them together with planetary signatures (Dethier & Bourrier 2023). As a complementary approach, we propose building data-driven estimates of the planet-occulted stellar lines, by cleanly separating spectral contributions from the planet-occulted stellar regions and from the planetary atmosphere.

This requires spectra to be corrected for instrumental and environmental effects, and made as close as possible to how they would be measured directly out of the planetary system. Instrumental pipelines, such as the Data Reduction Software (DRS) of the ESPRESSO and HARPS spectrographs, provide a first layer of calibrations (in wavelength and flux) and corrections (for the grating response, cosmic rays, etc). This reduction, however, is generally performed independently on individual exposures, making it difficult to separate systematic noise from stellar and planetary signals. Processing multiple exposures from a given star together offers the advantage of using the stellar spectrum as a common reference, to better identify and clean spurious features (e.g., YARARA, Cretignier et al. 2021; APERO, Cook et al. 2022). Transit observations offer the additional advantage of providing time series of consecutive exposures, typically within a given night. In that case, the stability of the stellar spectrum and observing conditions over the course of a few hours can be exploited to perform more advanced corrections (e.g., Snellen et al. 2010; SLOPPY, Sicilia et al. 2022). However, a variety of methods are employed in the community to reduce and analyze high-resolution transit spectroscopy datasets (e.g., Wyttenbach et al. 2017; Casasayas-Barris et al. 2019; Hoeijmakers et al. 2020), making comparisons between the results of these analysis difficult.

To make the most out of transit observations and ensure the reproducibility of results, we developed a novel workflow, the Advanced and Neat Techniques for the Accurate Retrieval of Exoplanetary and Stellar Spectra (ANTARESS), which is described in two main papers focusing on planet-occulted and planetary spectra, respectively. In Sect. 2, we introduce the general scope and concept of the workflow. Section 3 presents our description of the planetary system, in particular, the stellar surface. In Sect. 4, we explain how the workflow can be used to generate realistic mock datasets. The methods devised to correct and process high-resolution spectroscopy datasets are introduced in Sects. 5 and 6, while Sect. 7 describes generic methods used across the workflow. In Sect. 8, we present the methods dedicated to the analysis of the processed datasets and their application to archival data of two planets representative of high-resolution transit studies. Following a discussion of the performance of the workflow in Sect. 9, we present our conclusions in Sect. 10.

2 ANTARESS workflow

2.1 Concept

The ANTARESS workflow consists in a set of methods to process spectral time series from an exoplanetary system, obtained consecutively in one or several epochs. The main goals are to retrieve spectra from local regions of the stellar surface or the transmission and emission spectra from the planet atmosphere; however, the workflow can also be used more generally to clean and format stellar spectra. Depending on the objectives, the time series can thus overlap with the transit or eclipse of exoplanets, cover various parts of their orbit, or even monitor a star individually. The workflow was designed to process the spectral time series from various instruments homogeneously, to remain as close as possible to the original data, and to pay particular attention to the propagation of uncertainties and spectro-temporal resampling. We are presenting ANTARESS in several papers, describing the theoretical concepts behind the methods and their application to archival datasets to assess their performance and showcase new results. The present study describes the corrections that must be applied to the time series and the steps needed to extract the spectra of stellar regions occulted by a transiting planet. A companion paper will describe the subsequent methods developed to extract and analyze planetary spectra. Once the full ANTARESS workflow is described, we will release the implementation we are developing in an open-source, user-friendly Python pipeline.

Multiple epochs and instruments can be processed with the ANTARESS workflow to compare their outputs, analyze them jointly, or combine them and increase the detectability of extracted signals. Methods must be applied according to the sequences illustrated in Fig. 1, which depend on the input type of the data (e.g., spectral corrections are not relevant for cross-correlated spectra) and on the goal of the study (e.g., stellar spectra need to be cross-correlated for the analysis of the RM effect). Methods are grouped in three main categories, represented by colored boxes in Fig. 1. Input data are first corrected for instrumental and environmental effects (red boxes), then processed in a way that allows for the extraction of the stellar surface and planetary spectra (blue boxes). Finally, the data are analyzed to derive the orbital architecture of the system and properties of the stellar and planetary atmospheres (green boxes).

CORRECTIONS methods are applied to all epochs from a given instrument together, as some corrections are better constrained using multiple spectral time series. They are applied first in the sequence so that datasets from various instruments and epochs are made comparable, and can be processed and analyzed in the same way. The existing set of corrections was developped for ground-based, stabilized high-resolution spectrographs, and some corrections are specific to a given instrument (e.g., to deal with ESPRESSO interference patterns). The current workflow is thus applicable to the spectrographs listed in Table 1. The concepts behind the PROCESSING and ANALYSIS methods, however, are generic and applicable to time series from any medium- or high-resolution spectrograph. The ANTARESS workflow could thus be applied to additional instruments (Table 1), provided that the necessary correction methods are developed. PROCESSING methods aim at extracting specific types of spectral profiles and converting them in the format required for the various ANALYSIS methods available in the workflow. Profiles from a given epoch or from several epochs of a common instrument can be binned temporally and spatially at different steps of the processing, to create new series of profiles that can be analyzed at higher signal-to-noise ratio (S/N), or a single master profile representative of the star or planetary atmosphere. Most ANALYSIS methods are either aimed at fitting the extracted spectral profiles or their derived properties, and can be applied to various outputs of the PROCESSING workflow. Profiles and properties can be analyzed jointly over multiple epochs of different instruments to get the highest possible precision from their fit.

While the ANTARESS workflow eventually aims at retrieving the orbital architecture of planetary systems as well as properties of the stellar and planetary atmospheres, it requires a number of bulk and orbital properties for planets and their host stars to be fixed. The main properties that need to be defined to apply the workflow can be found in Table 2.

2.2 Data format and notations

The ANTARESS workflow can be applied to extracted echelle order (S2D) or order-merged (S1D) spectra, as well as cross-correlated spectra (CCFs). However, our approach is most relevant when applied to S2D spectra, which are closest to the measurements made with instrument detectors. We recommend that the S1D spectra and CCFs only be used as input for preliminary analyses, as we will show that the properties they yield for the stellar surface (Sect. 9) and planetary atmosphere can be biased compared to processing S2D spectra. Hereafter we use the generic term profile to refer both to spectra, defined as a function of wavelength, and CCFs, defined as a function of rv. Our methods are applicable to input spectral data defined in the air (typically from optical ground-based spectrographs like ESPRESSO) or vacuum (typically from space-borne and IR ground-based spectrographs but also, e.g., for CARMENES and EXPRES).

Throughout the paper, the symbol E refers to an epoch and the symbol t to the time at which a given exposure was taken; F refers to the flux emitted by the star or a region of the star in the direction of the observer, measured as a function of wavelength λ and at a distance D (Table 3). For the sake of clarity we write the spectral bin width dλ(λ, E, t) as dλ. We assume that spectra input in the workflow are defined as a function of wavelength ${\lambda }_{{\mathrm{B}}_{\odot }}$ in the solar system barycentric rest frame. The symbol $\overline{\lambda }$ refers to wavelengths bins with a much larger spectral width than the instrumental pixel size. Some methods Doppler-shift profiles between two rest frames using the relevant rv, defined as negative when an object at rest in one of the frames is moving toward the other. Finally, the true flux that would be measured without any perturbations can be written as: $$\mathrm{F}(\lambda ,\mathrm{E},\mathrm{t})={\mathrm{F}}_{\star }(\lambda ,\mathrm{E}){\delta }_{\mathrm{p}}(\lambda ,\mathrm{E},\mathrm{t}),$$(1)

where F_* is the stellar spectrum, and δ_p the sum of all contributions from the planet, namely, the ratio of the in-transit flux to the stellar flux, ${\delta }_{\mathrm{p}}^{\mathrm{tr}}$, the fraction of starlight reflected by the planet, ${\delta }_{\mathrm{p}}^{\text{refl}}$, and the ratio of the thermal planetary flux to the stellar flux, ${\delta }_{\mathrm{p}}^{\text{th}}$. The flux coming from the planet can thus also be written as $F_{\mathrm{p}}(\lambda ,E,t)=F_{\star }(\lambda ,E)({\delta }_{\mathrm{p}}^{\text{refl}}(\lambda ,E,t)+{\delta }_{\mathrm{p}}^{\text{th}}(\lambda ,E,t))$.

2.3 Resampling, error propagation, and fitting

The standard approach in the literature of transmission spectroscopy is to resample transit spectra on the same wavelength table during the reduction. However, several steps require to shift spectra between different rest frames (typically the star and planet frames), with shifts specific to each exposure. This combination of shifts and resampling has two main consequences: loss of spectral resolution and the introduction of correlated uncertainties. We thus propose two improvements to the processing of transit spectra. First, each spectrum is processed on its original spectral table as long as possible, namely, its spectral table is shifted without resampling. Second, when resampling has to be applied (typically to combine spectra together) we calculate the resulting correlations and propagate them afterwards.

We represent a spectrum as an array of flux densities over pixels, to which we associate a covariance matrix. Diagonal values of the matrix, noted as σ²(λ, E, t), correspond to the variance of the signal measured in the pixel at wavelength λ, for the exposure at time t of epoch e. While the ANTARESS workflow can be applied to input data affected by correlated noise, spectra provided by the DRS of standard spectrographs are usually uncorrelated and associated with diagonal covariance matrices. If uncertainties are not provided with measured spectra, a reasonable assumption for high-S/N data (≳50) consists in defining a photon-noise variance proportional to the number of photoelectrons received during an exposure.

The resampling of a spectrum on a different spectral table is performed through the application of an interpolation algorithm, which accounts for error propagation (Appendix A.1). To limit the spreading of correlation between the resampled pixels we propose a linear and a cubic interpolator, both of which preserve the integral of the flux on each original pixel (Fig. A.1). Fast processing and analysis can be carried out using linear interpolation, but the preferred approach to limit blurring is cubic interpolation.

Models can be fitted to the data products of the ANTARESS workflow using various ANALYSIS methods (Fig. 1), following the method described in Appendix A.2. Preliminary fits can be performed using Levenberg-Marquardt minimization and the variance of the data alone, but final fits should be performed with Markov chain Monte Carlo (MCMC) sampling and account for covariances to fully benefit from our approach and avoid biases.

2.4 Application

To illustrate the capabilities of our method we apply it throughout this paper to archival ESPRESSO transit observations of HD 209458b and WASP-76b. The different properties of these two systems make the retrieval of the spatially resolved stellar spectra of the sodium doublet (for HD 209458b), and of the stellar CCFs of neutral iron (for WASP-76b), particularly interesting. HD 209458b orbits a G-type star on an aligned orbit and has no atmospheric signature in the core of the sodium doublet (Casasayas-Barris et al. 2021; Dethier & Bourrier 2023). WASP-76b orbits an F-type star on a highly misaligned orbit and the stellar light is absorbed by neutral iron in its atmosphere (Ehrenreich et al. 2020). General properties of the systems required for our analysis are given in Table 2.

For each planet we exploit two transit datasets obtained with ESPRESSO on 20 July 2019 (epoch 1) and 11 November 2019 (epoch 2) for HD 209458b and on 03 September 2018 (epoch 1) and 31 October 2018 (epoch 2) for WASP-76b. Information about the scheduling, instrumental settings, and night conditions for these epochs can be found in Casasayas-Barris et al. (2021) for HD 209458b and in Ehrenreich et al. (2020) for WASP-76b. For an optimal application of the ANTARESS workflow to transit datasets we emphasize the importance of securing observations both before and after the transit. This allows building a reference spectrum for the unocculted star, and characterizing possible short-term variations due to the star, the Earth atmosphere, and the instrument that can then be corrected for over the entire dataset. One of the interests of using ESPRESSO data as a case in point is that it is representative of a series of stabilized, high-resolution echelle spectrographs used to observe stars and transiting planets to which our methods can be applied (Table 1). ESPRESSO is installed at the Very Large Telescope (VLT) at ESO’s Paranal site. The light is dispersed on 85 spectral orders¹ from 380 to 788 nm. Our present analysis starts with S2D spectra extracted from the detector images, corrected for, and calibrated by version 3.0.0 of the DRS pipeline (Pepe et al. 2021).

3 Description of planetary systems

3.1 Coordinates

We use specific coordinates to describe the position of the planets and stellar surface regions in the ANTARESS methods. The two main reference frames (Fig. B.1), with origin at the star center, are the sky-projected stellar rest frame (x is the node line of the stellar equator, y the sky-projected stellar spin, z the line-of-sight (LOS); Fig. 2) and the star rest frame (x is the stellar equator, y the stellar spin, z the LOS). We refer to Appendix B.1 for details about our other coordinate systems.

Particular care must be taken about the definition of transit windows. We consider that an exposure is in-transit if the planetary disk overlaps with the stellar disk during at least a fraction of the exposure. Most exoplanet host stars rotate slowly enough that they can be considered spherical, in which case the condition for overlap is r_sky < R_*+R_p, where r_sky is the farthest distance between the planet and star centers during the exposure, projected onto the plane of sky. The workflow can account for host stars that rotate fast enough to be oblate. In that case the transit condition is best solved numerically. We assess whether the planet limb intersects with the sky-projected photosphere, namely, whether the x_*sky and y_*sky coordinates of the discretized planet limb provide a solution to the quadratic describing the oblate photosphere (Eqs. (13) and (14) in Barnes 2009²).

3.2 Planetary atmospheric masking

If the transiting planet possesses an atmosphere the absorption lines of its constituent species may be present in the spectra processed from the full stellar disk. While extracting those lines is one of the end goals of the workflow, in earlier steps they act as contaminants of the stellar lines. We thus propose a procedure to exclude spectral ranges contaminated by narrow absorption lines from the planetary atmosphere. This requires a list of wavelengths for the electronic transitions to be masked and a rv window corresponding to the absorbed range in the planet rest frame (assumed common to all lines). First, we calculate the $z_{\text{sky}}^{\prime }$ component of the planet velocity vector in the sky-projected orbital reference frame, by applying the changes from Eqs. (B.3) and (B.4) to the orbital velocity coordinates (Eq. (B.2)). Then, we calculate an effective planetary orbital rv for each exposure, rv_p/_*, by oversampling $z_{\text{sky}}^{\prime }$ to account for its possibly asymmetrical variations around the exposure center. The excluded rv window, shifted to the star rest frame using the rv_p/_* values, is used directly with CCF profiles or converted into spectral ranges for each electronic transition to be masked. Pixels can be masked at various steps of the workflow, and are propagated afterward through our resampling algorithm.

We illustrate the interest of atmospheric masking throughout the processing of the WASP-76b dataset, excluding conservatively the range [−25; 10] km s⁻¹ contaminated by planetary iron lines (Ehrenreich et al. 2020). These lines are present in the F9 stellar mask used by Ehrenreich et al. (2020), through the ESPRESSO DRS, to compute the CCFs of WASP-76. To be comparable with their analysis, we use the F9 mask linelist to set the excluded transitions.

3.3 Description of the stellar surface

The PROCESSING and ANALYSIS methods require the calculation of various properties associated with the planet-occulted and disk-integrated stellar surfaces, in particular their broadband emission (e.g., Sect. 6.3) and spectral profiles (e.g., Sects. 4 and 8). Since photospheric properties are not spatially uniform and the regions occulted at the stellar limbs or by multiple planets may have complex shapes, we use a numerical approach to calculate accurately these properties.

3.3.1 Local profiles

Spectral lines emitted by an “elementary” region of the photosphere are defined by their intrinsic profile, the rv of the region (to Doppler-shift the profile) and its broadband intensity (to scale the profile in flux).

The rv of a region at position (x_*sky,y_*) (see Sect. 3.1) is defined as: $$r{v}_{\text{surf\hspace{0.17em}}/\star }=x_{\star \mathrm{s}\mathrm{k}\mathrm{y}}{v}_{\mathrm{e}\mathrm{q}}\sin i_{\star }\frac{\mathrm{\Omega }\left(y_{\star }\right)}{{\mathrm{\Omega }}_{\mathrm{e}\mathrm{q}}}+{\displaystyle \sum _j}{c}_{\mathrm{c}\mathrm{b},\mathrm{j}}{\mu }^j,$$(2)

with $\mathrm{\Omega }\left(y_{\star }\right)={\mathrm{\Omega }}_{\mathrm{eq}}(1-{\alpha }_{\mathrm{dr}}{y}_{\star }^2-{\beta }_{\mathrm{dr}}{y}_{\star }^4)$ the stellar rotation rate as a function of latitude (equal to Ω_eq at the stellar equator), described by a Sun-like differential rotation law, and c_cb,j the coefficients of a polynomial law describing convective blueshift variations (see Cegla et al. 2016 for more details).

We define white-light and chromatic intensities for each region. The former, based on average stellar properties over the full spectral band covered by the processed data, are used to scale cross-correlated profiles. The latter, which account for variations in the broadband stellar emission, are used to scale spectral profiles in a given band. Intensities are modulated by limb-darkening, described in our formalism through uniform, linear, quadratic, nonlinear, power-2, or solar laws with coefficients specific to the spectral band considered. They can further be modulated by gravity darkening, which has the spectral balance of a black body at photospheric temperatures calculated with the corrected formalism of Barnes (2009). We caution the importance of calculating μ = cos θ accurately in the case of an oblate photosphere, as θ is the angle between the LOS and the local normal to the photosphere, which differs in this case from the stellar radius. We thus calculate μ numerically without approximations, using the projection of the local normal³ on the LOS.

We propose three possibilities to define intrinsic line profiles:

• Analytical: calculated for properties that vary across the photosphere according to parametric laws of a chosen spatial coordinate (e.g., a linear variation of contrast with μ). These variations and the model profile can be informed by measured intrinsic profiles (Sects. 8.2.3 and 8.3). Analytical profiles can be broadened by macroturbulence, using a radial-tangential (Gray 1975, 2021) or anisotropic Gaussian (Takeda & UeNo 2017) kernel.

• Data-driven: intrinsic profiles extracted from transit time series are aligned into a common rest frame (Sect. 6.5) and spatially resampled at higher S/N (Sect. 7.3) into a single master profile, or into a series function of the resampling coordinate that can then be interpolated in the chosen region.

• Theoretical: using PYSME (Wehrhahn et al. 2023), informed by a 1D stellar atmospheric model matching the type of the host star and its microturbulence velocity, and by a VALD linelist (Piskunov et al. 1995; Kupka et al. 2000; Ryabchikova et al. 2015) generated with the host properties. Solar abundances are used by default (Asplund et al. 2009) but can be varied globally or for specific species to better match the lines of interest. The spectral lines of chosen species can further be calculated using pre-computed grids of non local thermal equilibrium (NLTE) departure coefficients (see the PYSME website and Amarsi et al. 2020, 2022). Theoretical profiles are calculated over a grid in μ to sample center-to-limb variations and be interpolated in the chosen region. This approach could be extended to 3D stellar atmospheric models in future versions of the workflow.

Analytical models are typically good proxies for simple profiles, as in CCFs, but only allow defining a single line. Data-driven and theoretical profiles are better options to define complex spectral lines or a full spectrum, and by construction the former directly traces the line profiles of the target star. Using observed data to provide direct estimates of the spatially resolved stellar profiles, with no assumption on their shape, is one of the main interests of the ANTARESS workflow.

3.3.2 Disk-integrated and planet-occulted profiles

The model disk-integrated stellar profile outside of the transit is defined by integrating the elementary profiles over a cartesian grid discretizing the entire spherical or oblate sky-projected photosphere (Fig. 2). Planet-occulted profiles are integrated over the region transited during the full duration of an exposure. Disk-integrated profiles within the transit are then calculated directly as the difference between the out-of-transit disk-integrated profile and the planet-occulted profiles in each exposure. If relevant for comparison with observations (Sects. 4 and 8), those various profiles are further convolved with the response of the corresponding instrument.

The extraction and interpretation of planet-occulted stellar profiles is central to the ANTARESS workflow. The corresponding models are thus calculated over a finer grid than the one used for the star, each occulting planetary disk being discretized with a resolution adjusted to its size in the processed spectral band. This spatial oversampling allows better resolving the partial occultation of the stellar limb and accounting for spatial variations of photospheric properties across occulted areas (Sect. 2). A temporal oversampling can further be applied when the chord transited during an exposure is longer than the spatial scale of stellar surface variations, by positioning the planetary grid at regular intervals along the chord to account for the blurring introduced by the planet motion.

Stellar line profiles from planet-occulted regions can be calculated at either low or high precision (Fig. 3). In the first case, brightness-averaged properties are calculated over each planet-occulted grid along the chord, and their mean over all grids define the exposure planet-occulted line profile. In the second case, elementary line profiles are summed over the cells of each planet-occulted grid, and then averaged over all grids. Our tests suggest that models at low precision without oversampling approximate well enough the planet-occulted line profile for slow rotators, planets with small transit depth, and short exposure times whereas models at high precision with oversampling are required for giant planets transiting fast-rotating or bright stars (when blurring due the planet size and motion is strong and detectable), and for multiple transiting planets. Details can be found in Appendix B.2.

4 Mock datasets

The description of the stellar surface (Sect. 3.3) and noise properties (Sect. 2.3) in the ANTARESS methods can be used to generate realistic datasets for the planetary system of interest, which can then be processed like real observational datasets for tests and predictions. Mock datasets can be used in particular to evaluate the detectability of RM signatures, as illustrated with the simultaneous transits of two exoplanets in Fig. 4.

The method is generic and allows creating a series of disk-integrated profiles, defined by their start and end BJD times, as they would be observed with one of the spectrographs implemented in Table 1. Local profiles are defined in spectral or CCF format (Sect. 3.3) and integrated over the stellar disk numerically (Sect. 3.3.2), accounting for planetary occultation during in-transit exposures.

These disk-integrated profiles are generated in the stellar rest frame over a spectral grid at the chosen sampling. They are then Doppler-shifted to the heliocentric rest frame, using the stellar Keplerian motion induced by all planets in the system and a chosen systemic velocity (Sect. 6.2), before being convolved with the spectrograph response. If the spectral sampling of the mock dataset is low, the disk-integrated profiles are first calculated at high resolution and resampled after instrumental convolution to avoid resolution loss.

Profiles are defined in spectral flux density F_mock, namely, the number of photoelectrons per spectral bin measured during an exposure, multiplied by the exposure time, by a mean flux density over the epoch, and by a constant instrumental flux calibration factor, g_mock. The last two fields allow introducing variations in the mock stellar emission measured with different instruments or epochs. The final profiles comparable to observations $F_{\text{mock}}^{\text{meas}}$ are drawn from a Poisson distribution with number of events set by F_mock. Uncertainties on $F_{\text{mock}}^{\text{meas}}$ are then attributed as ${\sigma }_{\text{mock}}^{\text{meas}}=\sqrt{g_{\text{mock}}{F}_{\text{mock}}^{\text{meas}}}$ (Sect. 5.1).

5 Corrections methods

Input spectra are assumed to be corrected for standard instrumental effects (flat field, blaze, background) and are put into a common format at the start of the workflow (Sect. 2.2). Our methods exploit housekeeping and complementary data⁴ to monitor the quality of the observations and inform the corrections. If available, spectra can be replaced by their sky-corrected values in the relevant spectral orders. Exposures and orders considered too noisy to be exploited are removed, and spurious spectral features (e.g., bad detector pixels or poorly-corrected telluric lines) are permanently masked, with thresholds depending on the dataset and objectives.

The current ANTARESS workflow (Fig. 1) can correct for telluric contamination, spectro-temporal flux balance variations, cosmic rays and spurious peaks, as well as ESPRESSO interference patterns. Additional methods could be added to deal with other instrument-specific corrections (e.g., fringing in near-IR spectrographs, breathing in HST/STIS) and stellar noise (e.g., spots, pulsations in early-type stars). In the following we use F^meas and F^corr to refer to spectra prior to, and after, a given correction. We note that most of the proposed correction methods rely on a fixed model to correct F^meas, and therefore do not propagate possible uncertainties associated with the model. Corrected spectra are made comparable between epochs of a given instrument, but still differ between instruments because of instrumental responses. This is taken into account when jointly interpreting datasets from different spectrographs together (Sect. 8).

5.1 Spectral flux calibration

Input spectra F^meas and their associated uncertainties σ^meas are generally provided in flux units, related to the number of photoelectrons N^meas cumulated during an exposure per unit of wavelength as: $$\begin{array}{rl}& F^{\text{meas\hspace{0.17em}}}(\lambda ,E,t)=g_{\mathrm{c}\mathrm{a}\mathrm{l}}(\lambda ,E,t)N^{\text{meas\hspace{0.17em}}}(\lambda ,E,t),\\ & {\sigma }^{\text{meas\hspace{0.17em}}}(\lambda ,E,t)=g_{\mathrm{c}\mathrm{a}\mathrm{l}}(\lambda ,E,t)\sqrt{N^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}(\lambda ,E,t)}.\end{array}$$(3)

At this stage, no resampling has been applied to the spectra and there is no covariance between their bins. The g_cal profiles represent calibrations applied by the instrument’s DRS to the raw photoelectrons count, typically to correct for the blaze function (i.e., the transmission of the spectrograph grating). While the application of g_cal in DRS sets spectra closer to their true spectral flux balance, it artificially changes the relative noise level between different spectral regions.

Let us consider a flat spectrum, measured with lower count levels at the edges of a CCD (e.g., N/10, with photon noise uncertainty $\sqrt{N/10}$) compared to its middle $(N,\sqrt{N})$ due to the blaze function. Correcting for the blaze yields similar flux levels (∝ N) across the detector, but uncertainties at the edge and middle become proportional to $\sqrt{10N}$ and $\sqrt{N}$, respectively, so that the error ratio between edge and middle is 10 times larger after flux calibration (see Fig. 5). This illustrates how combining different spectral regions in calibrated flux units, in particular when binning spectra over wide bands (Sect. 5.3) or cross-correlating them (Sect. 7.5), artificially inflates uncertainties compared to operating in raw flux units, decreasing the precision on derived properties.

To circumvent this issue, we estimate a flux calibration profile <g>_cal using the spectra themselves (Appendix C.1). This profile (Fig. 5) only approximates the original calibration from the DRS (in particular it is overestimated in regions of low S/N because we neglect readout noise), but allows temporarily scaling back spectra into raw flux units in the relevant methods (Sect. 7.2), while maintaining the relative flux balance between exposures to avoid introducing biases.

5.2 Telluric correction

In spectra obtained from the ground, the stellar or planetary lines of interest can overlap with absorption (or emission) lines from Earth’s atmosphere, especially broad molecular absorption bands. The ANTARESS workflow does not yet correct for telluric emission lines, whose amplitude and profile vary in complex ways during a night. Fiber-fed spectrographs often monitor the sky simultaneously with the target, so that emission lines can potentially be removed in sky-corrected spectra.

Different methods exist to correct for telluric absorption in the literature: contemporaneous spectral measurements of spectrophotometric standard stars close to the target (e.g., Vidal-Madjar et al. 1986), model-based algorithms (e.g., MOLECFIT, Smette et al. 2015; Kausch et al. 2015; tapas, Bertaux et al. 2014; TELFIT, Gullikson et al. 2014), or empirical methods (e.g., SYSREM, Mazeh et al. 2007; the PCA approach from Artigau et al. 2014; YARARA, Cretignier et al. 2021). We adapted in the ANTARESS workflow the Automatic Telluric Correction (ATC, Allart et al. 2022), used for the DRS of state-of-the-art spectrographs (ESPRESSO, NIRPS, HARPS). This model-based approach is a simple line-by-line radiative transfer code, which calculates the telluric spectrum of a single atmospheric layer based on the sky conditions and the physical properties of molecular lines (HITRAN database, Rothman 2021). The model accounts for the most relevant telluric molecules (H₂O, O₂, CH₄, CO₂), and could process additional species if needed. The ATC provides similar performances as other established methods for exoplanet transmission spectroscopy (SYSREM, Snellen et al. 2010; MOLECFIT, e.g. Allart et al. 2017; Seidel et al. 2019; Hoeijmakers et al. 2020; Sedaghati et al. 2021) but can be applied automatically to individual spectra, which better suits the ANTARESS approach.

Properties of the theoretical atmospheric layer are derived for each molecule by fitting a CCF of the model telluric spectrum to the equivalent CCF of an exposure spectrum, using housekeeping data from the spectrograph facility to inform the fit. Once atmospheric properties are derived for all molecules, we generate a global telluric spectrum δ_⊕ at the resolution and sampling of the observed spectrum, corrected as: $$F^{\text{corr\hspace{0.17em}}}({\lambda }_{{\mathrm{B}}_{\odot }},E,t)=\frac{{\mathrm{F}}^{\text{meas\hspace{0.17em}}}({\lambda }_{{\mathrm{B}}_{\odot }},\mathrm{E},\mathrm{t})}{{\delta }_{\oplus }({\lambda }_{{\mathrm{B}}_{\odot }},\mathrm{E},\mathrm{t})}.$$(4)

The method is detailed in Appendix C.2 and illustrated in Fig. 6. The WASP-76b and HD 209458b datasets were corrected for H₂O and O₂ telluric absorption. We note that the derived δ_⊕ are also used to exclude telluric lines from CCF masks built with the workflow (Sect. 8.1).

5.3 Flux balance corrections

This subset of methods deal with changes in the relative flux balance of spectra over time, which can be described in the most general case through a coefficient c: $${\mathrm{F}}^{\text{meas\hspace{0.17em}}}({\lambda }_{{\mathrm{B}}_{\odot }},\mathrm{E},\mathrm{t})=\mathrm{c}({\overline{\lambda }}_{{\mathrm{B}}_{\odot }},\mathrm{E},\mathrm{t}){\mathrm{F}}_{\star }({\lambda }_{{\mathrm{B}}_{\odot }},\mathrm{E}){\delta }_{\mathrm{p}}({\lambda }_{{\mathrm{B}}_{\odot }},\mathrm{E},\mathrm{t}).$$(5)

The planetary contribution δ_p could be isolated if c had a repeatable temporal behaviour (e.g., the breathing effect of the HST, whose achromatic flux variations are phased with the telescope orbit, Sing et al. 2019; Bourrier et al. 2017) or if they could be measured independently (e.g., using a comparison star). Since the current ANTARESS workflow focuses on ground-based observations, c represents at first order the “color effect” caused by Earth atmospheric diffusion (e.g., Bourrier & Hébrard 2014). The spectral balance of c varies over scales of tens to hundreds of Å, with variations more regular as a function of light frequency ν. Furthermore, c evolves during an epoch (as the target moves through the sky) as well as between epochs (as atmospheric conditions and the sky position of the star evolves). As a result, low-frequency variations of planetary origin cannot be separated from atmospheric diffusion and are removed by our correction. An estimate of δ_p must be re-injected later-on in spectra corrected for flux balance (Sect. 6.3). Importantly, the correction has no impact at short spectral scales. Because highresolution spectra are spectro-photometrically accurate at high spectral frequency, the true profiles of narrow stellar and planetary lines can thus be retrieved once spectra have been scaled at low frequency (Sect. 6.4). We note that our description for c may also account for chromatic fiber losses possibly caused by a wavelength-dependent seeing profile and its projection on the spectrograph fiber.

We reset the flux balance of the F^meas spectra to that of a reference spectrum expressed as F_ref(λ, E) = C_ref(λ, E)F_*(λ, E), considering that this reference does not necessarily have the same balance as the true stellar spectrum. An analytical model R^theo describing the flux balance variations between F^meas and F_ref (Fig. 7) is derived in ν space as described in Sect. C.3, and applied as: $$\begin{array}{rl}{F}^{\text{corr\hspace{0.17em}}}({\lambda }_{{\mathrm{B}}_{\odot }},E,t)& =F^{\text{meas\hspace{0.17em}}}({\lambda }_{{\mathrm{B}}_{\odot }},E,t)\frac{R^{\text{norm\hspace{0.17em}}}(E,t)}{R^{\text{theo\hspace{0.17em}}}({\lambda }_{{\mathrm{B}}_{\odot }},E,t)},\\ \text{\hspace{0.17em}where\hspace{0.17em}}{\mathrm{R}}^{\text{norm\hspace{0.17em}}}(\mathrm{E},\mathrm{t})& =\frac{\sum _{{\lambda }_{{\mathrm{B}}_{\odot }}}{F}^{\text{meas\hspace{0.17em}}}({\lambda }_{{\mathrm{B}}_{\odot }},E,t)d{\lambda }_{{\mathbf{B}}_{\odot }}}{\sum _{{\lambda }_{{\mathrm{B}}_{\odot }}}{F}_{\text{ref\hspace{0.17em}}}({\lambda }_{{\mathrm{B}}_{\odot }},E)d{\lambda }_{{\mathrm{B}}_{\odot }}}.\end{array}$$(6)

Corrected spectra have the same broadband profile as F_ref but keep the same global flux level as the original spectra to remain as close as possible to the measured counts (absolute flux scaling is performed in Sect. 6.3).

We first set all spectra in a given epoch to the same flux balance, taking their median spectrum as reference (see Appendix C.3 and Fig. 8). The correction for the dominant color effect is performed over all spectral orders together. A second-order correction can be applied in individual orders, with c tracing local variations in Earth absorption or in the instrumental response (blaze function, flat field). This correction must be used with caution, as it may also capture small-scale features of instrumental, stellar, or planetary origin, such as the broad wings of planetary absorption lines or the ESPRESSO wiggle pattern (Sect. 5.6). Because the WASP-76 and HD 209458 datasets are dominated by this latter effect, we do not run them through the spectral order correction.

When several epochs are processed, spectra must be set to the correct relative balance of the stellar spectrum between epochs so that line profiles in different spectral bands can be properly compared and possibly combined. We thus apply a second correction derived from the ratio between the median reference spectrum in each epoch and a final reference (Appendix C.3), which can be common to all epochs (if the stellar flux balance remains stable, Fig. 8) or account for variations between epochs. The final reference should further be set to a spectrum representative of the absolute stellar emission in each epoch if one wants to extract low-frequency emission spectra from the planet.

5.4 Cosmics correction

This method exploits the fact that an epoch is made of a series of consecutive exposures, assumed to have the same broadband spectral profile thanks to the previous steps, except for slight variations due to stellar activity and planetary absorption/emission. This makes it possible to compare a given exposure spectrum with the mean spectrum over adjacent exposures, F_comp, following the procedure described in Appendix C.4, and to identify and replace spectral bins contaminated by cosmic rays. A bin is considered to be cosmics-affected if: $$\begin{array}{rl}{F}^{\text{meas}}({\lambda }_{\star },& E,t)-F_{\text{comp}}({\lambda }_{\star },E)\\ & >{\alpha }_{\text{cosm}}max({\sigma }_{{\mathrm{F}}_{\text{comp}}}({\lambda }_{\star },E),{\sigma }^{\text{meas}}({\lambda }_{\star },E,t)).\end{array}$$(7)

Where α_cosm is a chosen threshold, adjusted depending on the quality of each dataset so as not to exclude noisy pixels as cosmics. The condition on σ^meas accounts for the processed spectrum possibly being much noisier than adjacent exposures in a given bin – in which case F^meas(λ_*, E, t) can deviate significantly from F_comp(λ_*, E) due to photon noise alone. This can occur over pixels corrected for deep telluric lines that shift during an epoch. Because cosmic rays can generate charges across several pixels before being absorbed, due to their trajectory and local charge bleeding of the detector, we allow pixels adjacent to a cosmics-flagged bin to be replaced as well. Corrected flux values are set to F_comp(λ_*, E) and their variance⁵ to σ_Fcomp(λ_*, E)².

5.5 Persistent peaks masking

This method deals with flux peaks induced, for example, by hot detector pixels that were not masked by the instrument DRS, or telluric emission lines that were not removed using sky-corrected data. In contrast to transient cosmics, the affected bins show spurious flux values persisting over time. We thus identify and mask pixels that deviate from the stellar continuum by more than a threshold over all exposures (typically hot pixels) or over consecutive exposures (typically telluric emission lines, strengthening at twilight as airmass increases) in an epoch.

Since no spurious peaks were found in the ESPRESSO data, we show in Fig. 10 an application to CARMENES transit spectra of WASP-156b (Bourrier et al. 2023) where strong telluric emission lines varying in strength but persisting over consecutive exposures could be identified.

5.6 Wiggle correction

ESPRESSO spectra are modulated by “wiggles” (Fig. 11), sinusoidal patterns arising from interferences within the Coudé Train optics that affect the four VLT units. Most observations are affected by low-frequency wiggles (amplitudes up to several %, periods between ~10–100 Å; Allart et al. 2020; Tabernero et al. 2021), with the possible occasional contribution of high-frequency wiggles (amplitude ~0.1%, period ~1 Å; Casasayas-Barris et al. 2021). The wiggles’ amplitude and period vary with wavelength and over timescales of minutes, making it difficult to determine their behaviour. Nonetheless, since wiggles modify the local flux balance between exposures and thus bias transmission spectra, several empirical corrections have been attempted:

• Cubic splines (Tabernero et al. 2021; Azevedo Silva et al. 2022) fitted over a specific spectral range.

• Sines (Borsa et al. 2021; Casasayas-Barris et al. 2021; Benatti et al. 2021; Seidel et al. 2021, 2022) combined over different spectral ranges with specific period, phase, and amplitude.

• A broad Gaussian smoothing filter (Merritt et al. 2020; Kesseli et al. 2022; Maguire et al. 2023; Zhang et al. 2022) or Savitzky-Golay filter (Allart et al. 2020) to capture low-frequency wiggles, and a narrow median filter to capture high-frequency wiggles (Casasayas-Barris et al. 2022).

While these methods can be efficient at removing wiggles over local spectral bands, they consider each exposure independently and apply blind corrections without understanding of the underlying wiggle pattern. Such separate, piecewise corrections could introduce variations between transmission spectra, and inadvertently remove stellar and planetary signals. Given the importance of ESPRESSO for the astrophysical community, we thus devised a method to describe the spectro-temporal behaviour of wiggles and correct them homogeneously. This method is generic enough that it can be adapted to similar patterns, such as HARPS-N ripples (Thompson et al. 2020), HARPS interference pattern (Cretignier et al. 2021), or fringing on NIR spectrographs (e.g., Guilluy et al. 2023).

Through trial and error we determined that wiggles are best described as a beat pattern between a dominant sinusoidal component with frequency ~3.8 × 10¹⁰ s (periods⁶ between ~13–54 Å) and a weaker one at ~3.1 × 10¹⁰ s (periods between ~16–65 Å). The components’ amplitude and frequency can be defined as polynomials of light frequency. The chromatic coefficients of these polynomials, as well as the phase reference of the components, can be expressed as a linear combination of the telescope cartesian pointing coordinates, defined by its azimuth and altitude. The full model is thus controlled by a set of pointing parameters for each component (see Appendix C.6.1). These parameters change when the LOS crosses the meridian, but the model remains continuous in value and derivative. We emphasize that the telescope guide star should not be changed during an epoch, as we found that it resets all pointing parameters (see WASP-76, Fig. C.5).

Because our wiggle model still depends on many free parameters, its characterization for a given dataset has to be performed semi-automatically. Our procedure (see Appendix C.6.2) follows iterative steps that aims at initializing the model close to its best fit. First, we screen transmission spectra manually to identify the spectral ranges that constrain the wiggles and to assess the strength of the two components. Then we sample the chromatic variations of the components’ frequency and amplitude in each exposure, and determine their optimal polynomial models. Those polynomials are used to initialize the spectral wiggle model fitted to each exposure. These fits yield time series of polynomial coefficients for the amplitude and frequency of the two components, as well as phase values, which are fitted in turn with the pointing coordinate function. Finally, the resulting set of pointing parameters is used to initialize the spectro-temporal wiggle model fitted to the full dataset. Original flux spectra are corrected for the wiggles through division with this global best-fit model (Fig. 11). The efficiency of the wiggle correction is illustrated in Fig. 12.

Our exploration of the wiggles hints at trends between targets, offering hope that our model can be further simplified with generic laws linked to the telescope/instrument properties or the observed sky region. Our findings that wiggles are dominated by two components with reproducible frequencies already provide information to identify the exact optical elements that are responsible for the interference pattern. A systematic exploration of the wiggle properties over all ESPRESSO transit datasets, which can be done in a homogeneous way with ANTARESS, would help to further characterize their behavior (Allart et al. 2020).

6 Processing methods

At this stage of the workflow (Fig. 1) spectral data are corrected for instrumental and environmental variations. We now describe the methods developed to extract disk-integrated and planet-occulted profiles in each epoch, and to process them into various formats. Keeping with the spirit of ANTARESS, data is processed as S2D spectra until otherwise required. They can then be converted into 1D spectra (Sect. 7.4) or cross-correlated with the appropriate mask (Sect. 7.5), and resampled into new time series or single master profiles (Sect. 7.3). The processing pathway depends on the objective and corresponding ANALYSIS methods (Sect. 8), such as measuring global photospheric properties from a master disk-integrated stellar spectrum, or performing a RM analysis using time series of planet-occulted CCFs. We will describe methods to further process S2D spectra until the extraction of planetary atmospheric profiles in a companion paper.

6.1 Stellar line detrending

This method corrects time series of disk-integrated stellar profiles for systematic trends identified between the line properties (rv residuals, FWHM, or constrast) and environmental/housekeeping parameters (e.g., time or S/N). This is done by fitting polynomials, sinusoids, or their combination, to the out-of-transit line properties (Sect. 8.2.3) and using the derived model to correct the entire time series, in particular during transit (see details in Appendix D.1). This approach works best when stellar baseline can be measured both before and after the transit.

We highlight the importance of CCFs, whose properties can be measured at high precision, to devise detrending models that can then be applied to individual spectral lines or full spectra (for a direct detrending of CCFs, see Bourrier et al. 2022, 2023). Using this approach (not detailed here) to correct the S2D spectra of WASP-76 and HD 209458’s first epochs for correlations of their CCFs with phase and S/N reduced the out-of-transit dispersion of their rv residuals from 2.8 to 1.9 m s⁻¹ (WASP-76) and 1.25 to 1.13 m s⁻¹ (HD 209458), and contrast (normalized to its mean) from 383 to 258 ppm (WASP-76) and 275 to 174 ppm (HD 209458).

6.2 Alignment in star rest frame

Disk-integrated profiles must be aligned in the rest frame of the star, taking care of shifting their individual spectral grid rather than resampling them on a common grid (Sect. 2.3). Complementary spectral profiles, such as telluric spectra (Sect. 5.2) or flux calibration profiles (Sect. 5.1), follow the same procedure so that they can be used for weighing (Sect. 7.2). Wavelengths in the star rest frame λ_* are calculated using the formula of the classical Doppler shift, for a stationary receiver located at the solar solar System barycenter and the star as the source moving in an arbitrary direction: $${\lambda }_{\star }(t)=\frac{{\lambda }_{{\mathrm{B}}_{\odot }}}{(1+\frac{r{v}_{\star /B_{\star }}(t)}{c})(1+\frac{r{v}_{B_{\star }/B_{\odot }}}{c})}$$(8)

where rv_*/B_* is the star Keplerian rv relative to the planetary system barycenter and $r{v}_{B_{\star }/B_{\odot }}$ is the systemic rv between the stellar and solar System barycenters. If disk-integrated profiles are in CCF format, the sum of these two rv is instead subtracted from their rv grids.

The Keplerian rv (Eq. (B.1)) should account for all bodies that induce a substantial reflex motion over the course of the epoch. This motion is small during the transit of planets on low-eccentricity orbits and remains below 500 m s⁻¹ for most known host stars, which is lower than the best resolutions available (1.6 km s⁻¹ for ESPRESSO in ultra high-resolution mode, Pepe et al. 2021). There are nonetheless a few known planets inducing larger motions, and accounting for this correction may become important with next-generation spectrographs.

The systemic rv can be known a priori from velocimetry, but it must be measured from each dataset. First, because it can vary between instruments and epochs. Second, to work in the rest frame of a specific stellar line, which may deviate from the overall stellar rest frame depending on the line shape and formation layer. The measurement of ${\mathrm{rv}}_{{\mathrm{B}}_{\star }/{\mathrm{B}}_{\odot }}$ is done by isolating a given stellar line or converting disk-integrated spectra into CCFs (Sect. 7.5), aligning the profiles by correcting for the Keplerian motion, binning them over out-of-transit exposures (Sect. 7.3), and deriving the centroid of the master stellar line (Sect. 8.2.1).

We neglected the Lorentz factor $1/\sqrt{1-(v/c{)}^2}$ that accounts for relativistic effects in the Doppler shift formula, as the contribution from the systemic motion is constant, that from the Keplerian motion amounts to ~0.04 cm s⁻¹ for a velocity of 500 m s⁻¹.

6.3 Broadband flux scaling

This method resets profile time series to their true relative flux level. This step if not necessary if absolute fluxes can be measured or if the planet contribution is not removed by the flux balance correction (Sect. 5.3), but is generally critical for ground-based observations.

At this stage the F^meas spectra have the low-frequency profile of the chosen reference F_ref but retain their original flux level. We first normalize them to a common integrated flux T F_ref, set to the median flux over the epoch or to a physically-motivated value specific to each epoch. Then, we scale the normalized profiles using light curves that match lc, the flux modulation observed from the system that we calculate using the BATMAN package (for simple transit light curves, Fig. 13), the ANTARESS stellar surface model (to account for simultaneous transits and complex photospheres, Fig. 14), or a measured/theoretical light curve imported into the workflow (to account for complex features, specific to a given epoch). The real light curves trace temporal variations of the broadband stellar intensity and planetary absorption/emission, and our scaling light curves can thus be calculated chromatically over ${\overline{\lambda }}_{\star }$. Final rescaled profiles (see the full procedure in Appendix D.2) can be expressed as: $$\begin{array}{rl}{F}^{\mathrm{s}\mathrm{c}}({\lambda }_{\star },E,t)& =F({\lambda }_{\star },E,t)C_{\text{ref\hspace{0.17em}}}^{\text{norm\hspace{0.17em}}}({\overline{\lambda }}_{\star },E),\\ \text{\hspace{0.17em}with\hspace{0.17em}}{C}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}({\overline{\lambda }}_{\star },E)& =\frac{{\mathrm{T}\mathrm{F}}_{\mathrm{r}\mathrm{e}\mathrm{f}}}{\sum _{{\lambda }_{\star }}{F}_{\mathrm{r}\mathrm{e}\mathrm{f}}({\lambda }_{\star },E)d{\lambda }_{\star }}.\end{array}$$(9)

Where the F spectral time series have the low-frequency profile and global flux level of F_ref, and the relative flux level of the true spectra.

Both the HD 209458b and WASP-76b datasets are scaled using BATMAN light curves. While the medium-resolution sodium signature in the transmission spectrum of HD 209458b (Charbonneau et al. 2002) may trace broad wings from the planetary atmospheric lines (Vidal-Madjar et al. 2011), the signatures in the line cores arise from an incorrect estimate of the planet-occulted stellar line (Casasayas-Barris et al. 2021; Dethier & Bourrier 2023). In-transit spectra around the sodium doublet thus mainly decrease in flux due to the planet atmospheric continuum in this range, and we apply an achromatic scaling using the constant transit depth and limb-darkening from Table 2. In the case of WASP-76b, we study spectra over the entire ESPRESSO optical range. In-transit scaling must thus account for variations in limb-darkening and planetary absorption at low spectral resolution over this range. We use the non-linear limb-darkening coefficients and transit depths derived from HST/STIS spectra (Fu et al. 2021), which both show strong variations between 300 and 600 Å (Fig. 13).

6.4 Differential profiles

This method allows retrieving profiles that contain the starlight occulted by a planet⁷ and filtered by its atmosphere. At this stage, disk-integrated profiles are aligned in the star rest frame, have a common spectral flux balance and correct relative flux levels. They can be decomposed between regions occulted or not by the planet as: $$\begin{array}{r}{F}^{\mathrm{s}\mathrm{c}}({\lambda }_{\star },E,t)=C_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}({\overline{\lambda }}_{\star },E)(F_p({\lambda }_{\star },E,t)\\ +{\displaystyle \sum _{i\notin i_{\mathrm{o}\mathrm{c}\mathrm{c}}}}{f}_i({\lambda }_{\star },E)S_i+{\displaystyle \sum _{i\in i_{\mathrm{o}\mathrm{c}\mathrm{c}}}}{f}_i({\lambda }_{\star },E)(S_i-S_i^p({\lambda }_{\star },E)),\end{array}$$(10)

where F_p is the light emitted/reflected by the planet (Eq. (1)) and f_i the surface flux density from region i of the stellar disk, which has a sky-projected surface S_i along the LOS. $S_i^p$ is the portion of S_i hidden by the planet, which can be written as the sum of the equivalent surfaces occulted by the opaque planetary disk $S_i^{p[thick]}({\overline{\lambda }}_{\star },E)$ and by the optically thin atmospheric annulus $S_i^{p[thin]}({\lambda }_{\star },E)$. The spectral dependence of these surfaces traces broadband variations of the atmospheric continuum and high-frequency absorption by the annulus. The surfaces are also time-dependent because of the partial occultation of the star during ingress and egress, of the Doppler shift of the atmospheric signal due to the planetary orbital motion, and of the different layers of planetary atmosphere that occult the star at different times of the transit (e.g., Ehrenreich et al. 2020).

Differential profiles are extracted as the difference between the unocculted star and exposure profiles (Fig. 15): $$\begin{array}{rl}& F^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}({\lambda }_{\star },E,t)=F_{\star }^t({\lambda }_{\star },E)-F^{\mathrm{s}\mathrm{c}}({\lambda }_{\star },E,t)\\ & =C_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}({\overline{\lambda }}_{\star },E)(-F_p({\lambda }_{\star },E,t)+{\displaystyle \sum _{i\in i_{\mathrm{o}\mathrm{c}\mathrm{c}}}}{f}_i({\lambda }_{\star },E)S_i^p({\lambda }_{\star },E)).\end{array}$$(11)

In contrast to standard approaches, we calculate a reference for the star, $F_{\star }^t$, directly on the spectral grid of each exposure to avoid unnecessary resampling and subsequent resolution loss or introduction of correlated noise (see the detailed procedure in Appendix D.3). Differential profiles are also extracted outside of the transit to retrieve the light potentially coming from the planet, F_p, and to analyze the noise distribution of the data (in the absence of planetary signal and stellar variability, out-of-transit F^diff values should distribute as Gaussian noise).

Differential profiles can then be processed in two different ways. Either all planetary contributions are removed to convert differential profiles into intrinsic stellar profiles, from which are derived stellar surface properties and the orbital system architecture (as described in the following sections), or differential profiles can be corrected for planet-occulted stellar lines to isolate $S_i^{p[\mathrm{thin}]}$ and F_p.

6.5 Intrinsic profiles

This method converts in-transit differential profiles F^diff into intrinsic stellar profiles F^intr. First, differential profiles are reset to a common flux level by correcting for broadband flux variations (Sect. 6.3). Then spectral ranges contaminated by the planetary atmosphere, if present, are identified and masked (Sect. 3.2). Finally, intrinsic spectra are corrected for chromatic deviations around the average rv of the occulted regions induced by broadband variations in the planet continuum. The resulting profiles can be expressed as: $$\begin{array}{rl}& F^{\text{intr}}({\lambda }_{\star },E,t)=C_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\text{norm}}({\overline{\lambda }}_{\star },E)\\ & \times \frac{I_{\mathrm{o}\mathrm{c}\mathrm{c}}({\lambda }_{\star },E,t){\alpha }_{\mathrm{o}\mathrm{c}\mathrm{c}}({\overline{\lambda }}_{\star },E,t)S_{\mathrm{o}\mathrm{c}\mathrm{c}}^{\mathrm{p}[\text{\hspace{0.17em}thick]}}({\overline{\lambda }}_{\star },E,t)}{1-lc({\lambda }_{\star },E,t)}.\end{array}$$(12)

Here, I_occ is the specific intensity emitted by the occulted region, and α_occ its broadband intensity variations. Because the spectral flux scaling lc can be expressed as a function of α_occ and the planetary continuum, $S_{\mathrm{occ}}^{\mathrm{p}[\text{thick]}}$, the above expression can be simplified to (for full details see Appendix D.4): $$F^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}}({\lambda }_{\star },E,t)=F_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}({\overline{\lambda }}_{\star },E)I_{\mathrm{o}\mathrm{c}\mathrm{c}}({\lambda }_{\star },E,t).$$(13)

Here, $F_{\text{ref}}^{\text{norm}}$ only contains components constant in time over the transit and with low-frequency spectral variations. Intrinsic profiles are thus independent of planetary occultation (at both low-and high- spectral resolution) and broadband stellar flux variations, and only trace variations in the spectral profile I_occ of the occulted stellar lines. In this way, the ANTARESS workflow allows mapping spatially the spectral emission from exoplanet hosts, providing direct observational constraints for theory and models of stellar atmospheres (Sect. 8.3.2).

6.5.1 Cross-correlated differential and intrinsic spectra

Intrinsic spectra can be cross-correlated with a weighted mask to create high-S/N CCFs (Sect. 7.5), from which series of stellar surface properties along the transit chord can be derived. The method can be applied with masks used by standard spectrograph DRSs, or preferably with a mask generated self-consistently from the processed data (see Sect. 8.1). “White” intrinsic CCFs built from a mask encompassing the full spectrum (Fig. 17) can first be used to study the stellar surface velocity field and system orbital architecture (Sect. 8.3.1). Intrinsic CCFs built from masks specific to a given species can then be used to study how the properties of the stellar atmospheric layers populated by this species deviate from the average photospheric properties probed by the white CCFs (Sect. 8.3.2). We emphasize that, in contrast to disk-integrated CCFs distorted by planet-occulted lines and planetary atmospheric lines, CCFs computed from intrinsic spectra only trace the stellar emission and can be directly interpreted with stellar line models with fewer risks of biases. Particular care is taken to deal with planetary contamination, as described in Appendix D.4.

When computing intrinsic CCFs we cross-correlate out-of-transit differential spectra with the same mask. In-transit differential CCFs, computed as the difference between intrinsic CCFs and their best-fit model (Sect. 8.2), are then rescaled by (1 – lc_CCF(rv, E, t))/lc_CCF(rv, E, t) (where lc_CCF(rv) is the CCF of lc(λ_*), Sect. 7.5) to be made comparable to out-of-transit differential CCFs. The resulting series of differential CCFs allow studying the quality of the fit to the intrinsic stellar line, and the noise distribution of the data over a full epoch (Sect. 8.3.1).

6.5.2 Alignment of intrinsic profiles

Intrinsic profiles are extracted in the star rest frame, with stellar lines Doppler-shifted due to the motion of photosphere. This method aligns them in a common frame, in which they can be resampled temporally/spatially (Sect. 7.3) to be analyzed at higher S/N (Sect. 8.2.2) or used to extract atmospheric profiles. We take care of shifting the individual spectral grid of each profile rather than resampling them on a common grid (Sect. 2.3). Wavelengths in the rest frame of the stellar surface λ_surf are calculated for the receiver located at the star rest frame and the local photosphere as the source: $${\lambda }_{\text{surf\hspace{0.17em}}}(t)=\frac{{\lambda }_{\star }}{(1+\frac{{\mathrm{r}\mathrm{v}}_{\text{surf\hspace{0.17em}}}({\lambda }_{\star },\mathrm{t})}{c})},$$(14)

where rv_surf(λ_*,t) is the rv of the planet-occulted region, subtracted from the rv grid of intrinsic profiles if they are in CCF format. We neglect relativistic effects due to the motion of the stellar photosphere, as they are negligible.

Values of rv_surf can be set to the shifts derived for each intrinsic line profile (Sect. 8.2.2). This approach is accurate but may be imprecise, or even impossible if the line cannot be fitted due to a low S/N or inconvenient planet masking (Sect. 3.2). Alternatively, rv_surf can be calculated for all planet-occulted regions with ANTARESS stellar surface model (Sect. 3.3), accounting for rv contributions relevant for the studied line (e.g., convective blueshift in the line spectral range).

7 Generic methods

We describe here generic methods that are used in several steps of the ANTARESS workflow.

7.1 Stellar continuum

The stellar continuum is required in particular by the PEAK MASKING (Sect. 5.5) and CCF MASKS (Sect. 8.1) methods. Appendix E.1 summarizes our approach to derive the continuum from the processed data (Fig. 18), using a simplified SNAKE sequence (smoothing, finding local maxima, measuring α-shape, interpolating envelop) from the RASSINE algorithm (see Cretignier et al. 2020b for details).

7.2 Spectral weight profiles

ANTARESS methods perform weighted averages of profile time series (Sect. 7.3), or overlapping orders in a given profile (Sect. 7.4), using the squared inverse of the flux error in each pixel as weight. This accounts for different pixel flux precision between profiles due to low-frequency variations in their overall flux level, and to high-frequency variations in their spectral flux distribution. However, because the uncertainty associated with a measured pixel flux is a biased estimate of its true error, we devised an original approach to estimate its value (Appendix E.2.1). Weight profiles on disk-integrated and intrinsic spectra measured over Δt are then defined as (Appendix E.2.2): $$\begin{array}{rl}{W}^{\text{sc}}\left({\lambda }_{\star }\right)& =\frac{1}{{\sigma }^{\text{sc,true}}{\left({\lambda }_{\star }\right)}^2},\\ W^{\text{intr}}(\lambda )& =\frac{(1-lc(\lambda ){)}^2}{{\sigma }_{\star }^{\text{meas}}(\lambda {)}^2+{\sigma }^{\text{sc,true}}(\lambda {)}^2},\\ \text{\hspace{0.17em}with}{\sigma }^{\text{sc,true}}\left({\lambda }_{\star }\right)& =lc\left({\lambda }_{\star }\right)\sqrt{{\tilde{g}}_{\mathrm{c}\mathrm{a}\mathrm{l}}\left({\lambda }_{\star }\right)F_{\star }^{\text{meas}}\left({\lambda }_{\star }\right)C_{\text{corr}}\left({\lambda }_{\star }\right)/\mathrm{\Delta }t}.\end{array}$$(15)

The disk-integrated spectra are averaged in the star rest frame, while intrinsic spectra are averaged in the star (λ = λ_*) or photosphere rest frame (λ = λ_surf). The flux calibration profile ${\tilde{g}}_{\text{cal}}$ (Sect. 5.1), color effect (through C_corr, Sect. 5.3), and flux scaling lc(λ_*) (Sect. 6.3) contribute to low-frequency noise variations, weighing disk-integrated and intrinsic spectra in the same direction through σ^sc,^true. The ${\tilde{g}}_{\text{cal}}$ profiles have a strong impact when averaging overlapping orders (Fig. 19).

High-frequency noise variations arise from tellurics (through C_corr, Sect. 5.2) and disk-integrated stellar lines (through $F_{\star }^{\text{meas}}$). Tellurics weigh in the same direction disk-integrated and intrinsic spectra, with lower weights for deeper lines. Because tellurics shift with the Earth’s motion in the star or surface rest frames, the noise and associated weight of a telluric-corrected pixel evolves over time. Stellar lines remain fixed in the star rest frame and only contribute to the weighing of intrinsic spectra in the surface rest frame. In that case, the disk-integrated lines in which intrinsic lines are imprinted shift between exposures and can induce strong noise variations in a given pixel (Fig. 20). Counter-intuitively, W^intr is inversely proportional to $F_{\star }^{\text{meas}}$ because the unocculted stellar spectrum acts as a (noisy) background light source when computing its difference with in-transit spectra. The parts of resulting intrinsic spectra that fall at the bottom of disk-integrated stellar lines are thus measured more precisely than those falling in the wings or continuum, since the flux and associated noise are lower there. When planetary atmospheric lines contaminate intrinsic spectra, the masked ranges also have a strong impact on the weighing due to their large shifts during transit (Sect. 7.3).

Low-frequency noise variations are generally accounted for in the literature by using the S/N, which approximates the true error on the disk-integrated flux over wide spectral bands. The impact of tellurics and stellar lines, however, is often ignored despite the biases it can introduce in averaged intrinsic stellar (Appendix G.3) and planetary spectra.

7.3 Temporal and spatial resampling

This method is used to calculate the weighted mean of a series of profiles, in new bins along a chosen dimension x: $$F(\lambda ,x,E)=\frac{{\displaystyle {\sum }_t}F(\lambda ,x(t),E)W(\lambda ,x(t),E)\mathrm{\Delta }x(t)}{{\displaystyle {\sum }_t}W(\lambda ,x(t),E)\mathrm{\Delta }x(t)}.$$(16)

Here, F and W are the flux and weight profiles associated to the exposure at time t (Sect. 7.2). Details on the resampling procedure are given in Appendix E.3.

Disk-integrated profiles are binned over orbital phase (temporal resampling), typically to define a reference spectrum for the unocculted star. Intrinsic profiles can be binned over phase to build a similar series at higher S/N, which may allow analyzing stellar lines in individual new exposures (Sect. 8.2.2) at the expense of blurring. Alternatively, intrinsic profiles can be binned over r_sky (spatial resampling) to build a series that maps the emission of the photosphere more directly. Indeed, with the assumption that stellar lines only vary with distance to the star center, intrinsic profiles at opposite phases around the time of inferior conjunction are symmetrical and can be combined into a higher-S/N profile in which planet-contaminated ranges are filled in. This approach is illustrated in Fig. 21. Intrinsic profiles resampled as a function of r_sky, or into a single master, can be used to tile ANTARESS stellar grid (Sect. 3.3.1) or to correct differential profiles for stellar contamination.

7.4 Conversion from 2D to 1D spectra

This method is applied to disk-integrated or intrinsic stellar spectra once they have gone through their latest processing stage (Fig. 1), to delay as much as possible resolution loss and the introduction of correlated errors. A converted disk-integrated spectrum can be seen in Fig. 18, and the procedure is detailed in Appendix E.4. The ANTARESS workflow resumes with the 1D spectral time series after conversion. One of the interests of this method is also to generate order-merged 1D spectra, with associated uncertainties, which can be used beyond ANTARESS (comparison with stellar models, template-matching to derive rv, analysis of interstellar medium absorption, etc).

7.5 Cross-correlation

This method generates CCF profiles by cross-correlating spectra with a binary mask that contains the central wavelengths and weights of a set of lines. Only mask lines that contribute to the cross-correlation in all exposures are kept, so that CCFs remain comparable. The method can be applied to disk-integrated or intrinsic stellar spectra, after which the ANTARESS workflow resumes with the CCF time series. In that case, we recommend generating masks specific to the processed datasets (Sect. 8.1). The method is also used to compute telluric CCFs from measured and theoretical spectra (Sect. 5.2). The procedure is detailed in Appendix E.5, with examples of CCFs shown in Fig. 21.

It is advised to compute CCFs at the rv resolution corresponding to the pixel width of the original spectra, to limit correlations introduced between CCF pixels. Nonetheless, the ANTARESS workflow propagates these correlations into the covariance matrices associated with CCFs, and accounts for them in their analysis, so that smaller rv steps can be used to gain in spectral resolution (Appendix G.1).

8 Analysis methods

We describe here the use and scope of the ANALYSIS methods, which exploit the data resulting from a given sequence of PROCESSING methods in spectral or CCF format, over a broad spectral range or a specific spectral line, and in individual, resampled, or joined exposures.

8.1 Stellar mask generator

Custom masks tailored to a given star have been shown to yield CCFs of improved quality (e.g., Bourrier et al. 2023). Thus, we propose a method, detailed in Appendix F.1, to generate cross-correlation masks directly from the processed data. It was adapted to the ANTARESS workflow from a version of the KITCAT method, and we refer to Cretignier et al. (2020a) for the underlying concepts. The linelist is built by identifying and selecting absorption lines in a 1D, normalized master spectrum, calculated over all epochs of a given instrument to keep their CCFs comparable and avoid biases. Masks are however specific to each instrument because their different spectral coverage and line spread functions change the shape and distribution of stellar lines. Each mask line is associated with a weight representative of the photonic error on its rv shift.

While masks are better-defined from high-S/N disk-integrated spectra, it may be interesting to derive them from intrinsic spectra in the case of fast-rotating stars, when disk-integrated lines are blended and difficult to isolate. Masks built from an entire spectrum allow generating white CCFs to analyze global properties of the photosphere (e.g., its rotational velocity, Sect. 8.2.1). By cross-matching the white linelist with a VALD linelist specific to the target star (Piskunov et al. 1995; Kupka et al. 2000; Ryabchikova et al. 2015), it is further possible to define masks specific to a given species, which trace different layers in the stellar atmosphere⁸. The interest of this approach is enhanced by ANTARESS’s method to define intrinsic CCFs as a function of stellar surface coordinates (Sect. 7.3). Series of such intrinsic CCFs, built for a suitable set of species would, in effect, map the stellar atmosphere vertically and horizontally.

8.2 Analysis of individual stellar profiles

We propose methods to derive properties from disk-integrated (Sect. 8.2.1) or intrinsic (Sect. 8.2.2) stellar profiles in individual exposures. Analyzing original time series allows assessing their quality and identifying variations in the derived stellar properties (Sect. 8.2.3). Analyzing new series, resampled within an epoch or over several epochs (Sect. 7.3), allows performing these analysis with higher precision at the expense of temporal/spatial resolution.

Properties are directly measured on a stellar line (e.g., bissector, equivalent width) or derived by fitting the stellar profile with one of the models described in Appendix F.2.1. Models can be defined analytically, to estimate first-order properties controlling the shape and position (e.g., FWHM, contrast, rv) of a single stellar line, or numerically to reproduce complex line shapes or to fit a portion of the stellar spectrum. When the analysis aims at deriving correct physical properties or comparing line properties between different instruments, the model stellar profile is convolved with the relevant instrumental responses.

8.2.1 Disk-integrated stellar lines

Most disk-integrated lines of G and K-type stars are well fitted with Gaussian profiles, especially when cross-correlated over many transitions (Fig. 22), while lines with damping wings or asymmetrical shapes can be reproduced with Voigt or skewed Gaussian profiles. Most M dwarf lines are blended and difficult to fit individually, but result in sidelobed CCFs that are well reproduced with Double-Gaussian profiles (Bourrier et al. 2018). The disk-integrated lines of early-type stars may be shaped by rotational broadening or even gravity darkening, and cannot be reproduced with these simple analytical profiles. In that case, the lines can still be fitted numerically to constrain local line properties, the photosphere geometry, and its intensity variations. Fitting the disk-integrated CCF of fast-rotating stars may constrain v sin i_* independently of the value derived from the RM analysis (e.g., Bourrier et al. 2020).

Numerical profiles can further be used to constrain the abundance of a line species. To illustrate this we fitted the sodium doublet of HD 209458 in the out-of-transit 1D spectrum (Sect. 7.4) averaged over both epochs (Sect. 7.3). Profiles were tiled over ANTARESS stellar grid (Sect. 3.3.1) using the rotational velocity and limb-darkening properties from Table 2, and series of theoretical intrinsic profiles calculated with a plane-parallel MARCS2012 stellar atmosphere model (ξ = 1 km s⁻¹), a VALD linelist for HD 209458, an overall metallicity [Fe/H] = 0.02, and the NLTE departure coefficients for sodium from Amarsi (2022). The best fit (Fig. 23) is obtained for a model sodium abundance of 6.05. While the NLTE grid reproduces better the overall shape of the sodium doublet (see also Casasayas-Barris et al. 2021; Dethier & Bourrier 2023), the core and wings of the lines cannot be well fitted at the same time, leaving residuals on the order of 2%. This discrepancy could be due to the contribution of the chromosphere in the core of the sodium lines (Bruls & Rutten 1992), highlighting the need for refinements in stellar atmosphere models that can be driven by ANTARESS measurements of spatially resolved line profiles (Sect. 8.3.2).

Stellar properties are typically derived from out-of-transit exposures to prevent biases induced by the transiting planet. Series of line shape properties (Sect. 8.2.3) allow detrending the disk-integrated profiles from systematics over the full epoch (Sect. 6.1), and rv time series allow us to check the quality of the Keplerian model and evaluating the presence of an RM anomaly. In-transit disk-integrated profiles are then used to extract planet-occulted profiles, rather than to analyze their anomalous properties. However, useful information can still be derived from those profiles if they can be masked for the spectral ranges contaminated by the planetary atmosphere (Sect. 3.2) and planet-occulted lines (Appendix F.2.2). The masked disk-integrated stellar lines are then not distorted by the planet anymore, and can be fitted with the same model as for the unocculted star. This allows deriving stellar properties homogeneously over the entire epoch, which can be of particular interest to study stellar variability at short timescales (e.g., pulsations, Wyttenbach et al. 2020). This approach is only applicable with disk-integrated lines broad enough to remain defined despite the masking. We illustrate this with the rv derived from WASP-76 CCFs, masked for planetary absorption lines (Fig. 22). The orbital motion of the planet combined with the width of its absorption signature result in many in-transit profiles being too contaminated to be accurately fitted (Fig. 17).

8.2.2 Intrinsic stellar lines

Intrinsic stellar lines are measured at lower S/N than disk-integrated lines and have simpler shapes, especially when averaged over several electronic transitions, such that a Gaussian profile typically reproduces them well (Fig. 24). For planets that are large or occult broad areas during an exposure, compared to the spatial scale of line profile and rv field variations over the photosphere, intrinsic profiles should however be fitted with numerical models (Sect. 3.3.1).

By construction (Sect. 6.5), ranges contaminated by the planetary atmosphere are excluded from intrinsic profiles (Fig. 17). Our fitting method, based on a Bayesian approach (Appendix F.2.2), allows retrieving information on the local stellar line even when it is partially masked (Fig. 25) or hidden in the noise (typically at the stellar limbs), and assessing which exposures can be exploited in time series analysis (Sect. 8.3). For both WASP-76 and HD 209458 only the very first and last exposures at the limbs had to be excluded, plus those exposures of WASP-76 in which the intrinsic line is fully contamined by the planet (Fig. 17).

Series of properties derived from individual intrinsic profiles can be analyzed with a dedicated method (Sect. 8.2.3) to map the local stellar photosphere. The line centroids trace the surface rv field, controlled by the stellar rotational velocity, convective blueshift, and differential rotation. Morphological properties, such as width and contrast, trace center-to-limb variations in their spectral profile. Average properties of local stellar lines, or a species abundance, can further be derived by fitting a master of all intrinsic profiles. Comparing the width of master intrinsic and disk-integrated CCFs, or fitting the latter using the stellar grid tiled with an intrinsic profile series, provides complementary constraints on stellar rotation. Since intrinsic lines may vary in shape along the transit chord, it is nonetheless preferred to fit them jointly with a global model (Sect. 8.3).

8.2.3 Line profile variations

Once properties have been derived from stellar lines in individual exposures, they can be studied as a series to characterize their variations. The method proposed here allows quantifying anomalous trends in disk-integrated properties to inform line detrending (Sect. 6.1), defining the models best describing intrinsic properties along the transit chord (Sect. 8.3), and identifying deviations to the series from noise-related outliers or activity-related features (spots, pulsations, etc). We detail in Appendix F.2.3 the various models that can be fitted to line property series, as a function of ambient, data-related, or stellar variables for disk-integrated properties, or as a function of stellar coordinates for intrinsic properties.

White-light analysis. Fig. 25 shows the intrinsic line contrast, FWHM, and rv in white CCFs from the regions occulted by the two planets. The surface rv of HD 209458 are consistent with those from Casasayas-Barris et al. (2021), but more similar between the two epochs at the limbs. The surface rv of WASP-76 are overall consistent with those from Ehrenreich et al. (2020), but differ sufficiently to impact the global fits to intrinsic profiles (see the next section). For both systems we found no evidence for convective blueshift. While there is a hint of differential rotation for HD 209458, it is driven by measurements at the limb that are uncertain due to imprecise transit timing and impact parameter (Casasayas-Barris et al. 2021). Hereafter, we thus assume that HD 209458 and WASP-76 rotate as solid bodies. Their line shape properties show clear variations along the transit chords, well described by polynomials of the distance to star center (linear for HD 209458, quadratic for WASP-76). The contrast of the line decreases toward the limbs as we probe cooler stellar layers while the presence of convective cells, which scan different velocities across the different depths of the layers, is likely the reason for the increase in line FWHM toward the limbs. We note that WASP-76 surface properties can be retrieved in more exposures than disk-integrated properties (Fig. 22), because intrinsic lines are not broadened by stellar rotation and overlap less with planet-contaminated ranges. Nonetheless, some intrinsic profiles that are partially masked cannot be accurately fitted in individual exposures, highlighting the interest for a joint fit to the profile series. One of the aims of analyzing line profile variations is to initialize this joint fit, from which are derived final properties for the stellar photosphere and planet architecture.

Spectral analysis. Fig. 26 shows the contrast, FWHM, and rv of the intrinsic sodium doublet from HD 209458. Gaussian profiles were fitted to CCFs of the doublet, using a mask with equal weights to match the similar line depths (Fig. 27). The sodium line width and contrast are larger than for iron lines in white CCFs (Fig. 25), and vary oppositely with distance to the star center (they get narrower and deeper toward the stellar limb). Furthermore, while the sodium line rv follow the same trend as the average photosphere they display a significant redshift common to both epochs. The same analysis performed on deep and narrow stellar lines in the vicinity of the sodium lines shows rv series consistent with the photospheric model, demonstrating that this redshift traces processes behind the formation of the doublet.

8.3 Analysis of joint intrinsic stellar profiles

This method builds on the RMR technique (Bourrier et al. 2021), whose philosophy is to fit jointly series of intrinsic profiles (Appendix F.3.1) rather than fitting individual profiles and their derived properties. The ANTARESS workflow offers two major updates: model intrinsic profiles are calculated with higher precision (Sect. 3.3.2) and can be fitted to both CCF and spectral profiles. The latter point allows constraining first the overall properties of the stellar surface (Sect. 8.3.1), and then investigating how a given spectral line compares with them (Sect. 8.3.2).

8.3.1 White-light analysis (RM revolutions)

Rather than applying the RMR to profiles derived from input disk-integrated CCFs (e.g., Bourrier et al. 2021, 2022), the ANTARESS workflow exploits CCFs calculated from extracted intrinsic spectra (Sect. 6.5.1), which trace directly the properties of planet-occulted regions. White intrinsic CCFs, calculated at high S/N with a mask of quiet and well-defined stellar lines (Appendix F.1), can constrain the architecture of the system (sky-projected λ or 3D ψ spin-orbit angle, stellar inclination i_*, scaled semi-major axis a_p/R_*, orbital inclination i_p), the global velocity field of the photosphere (stellar rotational velocity, v_eq sin i_*, differential rotation, convective blueshift), and the average shape of its absorption lines (contrast, FWHM, etc.). We detail in Appendix F.3.2 how and when these various properties can be constrained.

Here, we apply a RMR fit to the joined epochs of HD 209458 and WASP-76, using intrinsic CCFs calculated with ESPRESSO F9 mask. Based on the analysis in Sect. 8.2.3 we assume solid-body rotation for the stars, so that λ and v_eq sin i_* are the only free parameters of the surface rv model. We further define the contrast and FWHM of the theoretical line profiles as polynomials of distance to star center, modulated by coefficients specific to each epoch. Spatial oversampling is applied to the model exposures (Sect. 3.3.2), although the blurring induced by the planets’ motion and the photospheric areas they occult is negligible. Results from the best-fits, which yield reduced χ² = 1.1, are highlighted in Fig. 28. The stability of the stars, the quality of the processed data, and the good agreement between CCFs and their best-fit is further visible through the homogeneity and low dispersion of their residuals.

HD 209458. The best fit yields $\lambda =1.055{}_{-0.077}^{+{0.073}^{\circ }}$ and v_eq sin i_* = 4.272±0.007 km s⁻¹ with a similar precision as Casasayas-Barris et al. (2021) (λ = 1.58 ± 0.08° and v_eq sin i_* = 4.228±0.007km s⁻¹), who used the same datasets and system properties but applied a reloaded RM analysis (Cegla et al. 2016) to input disk-integrated CCFs. Our results are likely consistent, considering that more realistic error bars on λ and v_eq sin i_* (accounting for the uncertainty on the transit timing and impact parameter) are on the order of 1° and 100 m s⁻¹, respectively (Casasayas-Barris et al. 2021).

WASP-76. The configuration of the system, with an impact parameter close and consistent with 0 (Ehrenreich et al. 2020), is such that the planet-occulted surface rv and their derivative approximate during transit (when orbital phase ϕ ~ 0, see Eq. (2)): $$\begin{array}{rl}& r{v}_{\text{surf\hspace{0.17em}}/\star }(\varphi )=v_{\mathrm{e}\mathrm{q}}\sin i_{\star }{a}_p/R_{\star }\cos (\lambda )2\pi \varphi \\ & \dot{r}{v}_{\text{surf\hspace{0.17em}}/\star }(\varphi )=v_{\mathrm{e}\mathrm{q}}\sin i_{\star }{a}_p/R_{\star }\cos (\lambda )2\pi .\end{array}$$(17)

Thus, v_eq sin i_* cos(λ) is constrained in the same way by the slope and level of the rv series and there is a strong degeneracy between λ and v_eq sin i_* (Fig. F.2). This is in part why we derive $(\lambda =-{37.1}_{-20.6}^{+{12.6}^{\circ }};v_{\mathrm{eq}}\sin i_{\star }={0.86}_{-0.19}^{+0.13}\mathrm{km}{\mathrm{s}}^{-1})$ while Ehrenreich et al. (2020) derive (λ = 61±7°; v_eq sin i_* = 1.5±0.3 km s⁻¹), since both yield v_eq sin i_* cos(λ) ~ 0.71 km s⁻¹. Furthermore our derived rv series is slightly offset compared to Ehrenreich et al. (2020) (Fig. 25), crossing 0 km s⁻¹ after mid-transit rather than before, which requires the planet to transit the blueshifted stellar hemisphere longer (λ in between −180 and 0°, Fig. 29) than the redshifted one (λ in between 0 and 180°). These differences in rv series likely arise from our use of CCFs calculated from intrinsic spectra, rather Ehrenreich et al. (2020)’s use of intrinsic CCFs derived from disk-integrated CCFs (Sect. G.6). Finally, in an attempt to break the degeneracy, Ehrenreich et al. (2020) set a strong prior on v_eq sin i_* (1.61±0.28 km s⁻¹) derived as the quadratic difference between the FWHM of the stellar disk-integrated and planet-occulted master CCFs. Combining the average FWHM of our planet-occulted lines with the width of ESPRESSO line spread function (LSF), we derive a similar v_eq sin i_* = 1.75 km s⁻¹. However, our in-depth analysis show that the FWHM of WASP-76 intrinsic lines vary by at least 2 km s⁻¹ between the center and limbs of the star (Fig. 25). As a result, the rotational velocity estimated by Ehrenreich et al. (2020) from the average planet-occulted line, and in turn their spin-orbit angle, was biased.

To try and mitigate the degeneracy between λ and v_eq sin i_* we applied a RMR fit with ap/R_* and i_p free to vary, using photometry priors from Ehrenreich et al. (2020) $(a_p/R_{\star }={4.08}_{-0.06}^{+0.02}$, $i_p={89.623}_{-0.034}^{+{0.005}^{\circ }})$. We derived values consistent with the priors $(a_p/R_{\star }={4.059}_{-0.033}^{+0.044};i_p={89.604}_{-0.015}^{+{0.024}^{\circ }})$, showing that the intrinsic CCFs alone cannot constrain these orbital properties. Refining the orbital architecture of the WASP-76 system thus requires both additional transit spectroscopy and photometry. We note that our results change neither the conclusion from Ehrenreich et al. (2020) that WASP-76b is on a misaligned orbit, nor their analysis of the planet atmospheric signal since their surface rv model, used to mask the planet-occulted stellar lines, remain within 100 m s⁻¹ from our best-fit.

8.3.2 Spectral analysis

Analyzing intrinsic line series in specific transitions allows measuring properties of the corresponding species and its formation layer. We illustrate this approach with 1D intrinsic spectra of HD 209458 in the sodium doublet (Sect. 7.4). Numerical profiles were calculated like the white RMR fit (Sect. 8.3.1) but using theoretical intrinsic spectra, calculated with the disk-integrated properties and a free sodium abundance (Sect. 8.2.1). These theoretical spectra were Doppler-shifted according to the derived photospheric rv model (Sect. 8.3.1), allowing for a free offset as hinted in Sect. 8.2.3.

The best fit (Fig. 30) is obtained for an abundance of 6.056±0.007 and a Doppler offset of 849±27 m s⁻¹, consistent with the values derived from disk-integrated lines. This suggests that the sodium line profiles retrieved along the transit chord of HD 209458b are representative of the photosphere in the observed epochs, with the star displaying few active/spotted regions in this spectral band. While the NLTE pySME/MARCS stellar model matches well the sodium line shape it does not account for a Doppler offset like the one we measured. Furthermore, the theoretical model reproduces the decrease in line width observed from center to limb (Sect. 8.2.3) but simulates a decrease in contrast opposite to the measured variation (Fig. 27). These discrepancies emphasize the interest of confronting stellar models (e.g., Plez 2012; Wehrhahn et al. 2023) to measurements of local stellar lines, as provided by the ANTARESS workflow. We also note that our sodium abundance for HD 209458 is lower than the solar value (~6.17) used in theoretical transmission spectra accounting for stellar contamination (Casasayas-Barris et al. 2021; Dethier & Bourrier 2023). The extraction and characterization of planet-occulted stellar lines is thus also of strong interest to the interpretation of transits.

9 Performance analysis

We investigated the impact of various methods in the ANTARESS workflow on the quality of the outputs, in particular intrinsic CCFs and their RMR fits. We summarize here the results of this investigation, detailed in Appendix G:

• Resampling: we propose the first workflow that resamples spectra on a common grid as late as possible, reducing blurring and correlated noise to a minimum. Resampling has negligible impact on methods that modify the low-frequency stellar spectral profile, such as the color effect or wiggles, but it smears sharp spectral features at various steps of the workflow. While variations in stellar line properties and RMR results remain marginal for the studied datasets, they would be stronger and more significant for stars with larger Keplerian motion, narrower or asymmetrical stellar lines, and higher S/N spectra. This makes our approach particularly relevant for next-generation spectrographs on the ELT.

• Covariance: the optimal processing starts introducing correlations between pixels when extracting differential profiles, in which case the maximum covariance in cross-correlated intrinsic stellar lines remains small enough that it does not impact their analysis. Even if profiles are resampled earlier in the workflow, covariances are propagated efficiently as banded matrices and accounted for in fits, which is a major novelty of our approach. This allows oversampling CCFs and better resolving the stellar line profiles, which may be of interest with lower instrumental resolution than ESPRESSO and when lines are asymmetrical.

• Weights: neglecting the contributions of instrumental flux calibration and telluric lines (for disk-integrated spectra), or disk-integrated stellar lines (for intrinsic spectra), to the true error estimated in a given pixel can bias significantly its weighted spatial/temporal mean.

• Flux balance: a spectrum-wide correction of the color effect is critical for broadband spectral analysis and the computation of unbiased stellar CCFs. Capturing smaller-scale flux balance variations, while it may further improve the shape of CCFs, is however not critical for rv derivation and the interpretation of the RM effect.

• Wiggles: a global, homogeneous correction is necessary to remove biases when analyzing ESPRESSO spectra. However, wiggles appear to smooth out when computing stellar CCFs. While their correction may still improve the shape of the CCFs, it is not critical for rv derivation and the interpretation of the RM effect.

• Intrinsic CCFs: our workflow offers the possibility to keep data in spectral format as long as desired, so that CCFs can be directly calculated from intrinsic spectra rather than derived from disk-integrated CCFs. We recommend using the former, as the latter distorts the shape of intrinsic stellar lines and may bias the interpretation of the RM effect.

• Chromatic scaling: accounting for broadband variations in planetary absorption and stellar emission when rescaling stellar spectra may remove biases from their CCF profiles and the interpretation of the RM effect. We recommend using chromatic scaling with high-S/N data and planets with strong broadband atmospheric features.

• CCF masks: generating custom masks from the processed data allows building CCFs of comparable or higher quality than standard masks, and performing self-consistent RM analysis.

10 Summary and conclusions

The ANTARESS workflow consists of a set of rigorous methods designed to process high-resolution transit datasets homogeneously, ensuring a reproducible extraction and analysis of stellar and planetary spectra. Particular care has been taken to remain as close as possible to the original data, by reverting to raw flux units when relevant, propagating covariances, and avoiding spectral resampling unless necessary.

We introduce a set of methods meant to correct echelle spectra obtained from the ground for environmental and instrumental effects, which are mainly related to telluric absorption and flux balance variations induced by Earth’s atmosphere, cosmic rays, and optical interferences. A first-order correction for these effects is sufficient to analyze stellar line rv and perform RM analysis, but a finer correction reduces biases in spectral profiles and increases CCF quality. We then describe our methods for correcting trends in disk-integrated stellar lines, so that they are comparable and can be processed similarly by aligning disk-integrated spectra in the star rest frame, rescaling them to their relative chromatic flux level, and extracting planet-occulted spectra.

The workflow can be used to analyze stellar lines, in which case the planet is considered as a contaminant. Planet-occulted spectra are corrected for flux and rv modulations induced by the planetary continuum and then they are masked for planetary atmospheric lines. The resulting series of intrinsic spectra directly trace the properties of the photosphere along the transit chord, without any residual contamination by the planet. We highlight that the flux contrast between planet-occulted regions and the full stellar disk is better at the bottom of disk-integrated stellar lines, making slow-rotating stars with deep lines optimal cases to measure intrinsic stellar lines at high precision. We propose several methods to analyze disk-integrated and intrinsic stellar lines, using simple analytical profiles or a numerical model of the photosphere. This model, which is used as part of the data processing and to generate mock datasets, is generic enough that the workflow can be applied to any type of stars, from cool and slow-rotating dwarfs to oblate and gravity-darkened giants.

The analysis of intrinsic stellar profiles is first performed on individual exposures, to assess their variations across the stellar surface. Then, the profile series is fitted with a joint model informed by the first step, to derive properties at higher precision. This two-step RMR analysis can be applied to “white” CCFs computed over the full spectra, to constrain the rv field of the photosphere and the system orbital architecture. Once these global properties have been derived, they can be fixed and the analysis can then be applied to individual spectral lines, or even to CCFs over a subset of lines, to derive properties of a specific species or stellar atmospheric layer. We emphasize that the workflow allows us to cross-correlate intrinsic spectra directly, avoiding biases associated with the processing of disk-integrated CCFs.

To illustrate the use of the workflow, we analyzed archival ESPRESSO datasets from the iconic stars HD 209458 and WASP-76. A white-line analysis of their intrinsic spectra revealed average stellar lines that get shallower and broader toward the limbs. We confirmed the published values for the sky-projected spin-orbit angle and stellar rotational velocity of HD 209458 (Casasayas-Barris et al. 2021) and unveiled biases in the derivation of these properties for WASP-76 by Ehrenreich et al. (2020). While our revision does not change the planet atmospheric signature derived by the latter authors, it highlights the need for a robust processing of transit spectroscopy datasets. We then carried out a spectral analysis of HD 209458’s sodium doublet, detecting center-to-limb variations of the intrinsic lines opposite to those of the average photospheric lines. The intrinsic sodium abundance is consistent with the disk-integrated one, but lower than the solar abundance used in previous studies of the planetary atmosphere. Finally, we found divergences in shape and position compared to theoretical spectra, highlighting the interest of the workflow to confront stellar atmospheric models with spatially resolved stellar spectra.

The ANTARESS methods are generic and can be used outside of the workflow. For example, in this work we show how a proper calculation of weighted means between different spectra is critical to avoid biases. Combined with our method to convert 2D spectra into 1D spectra and CCFs, it allows for the calculation of accurate master profiles of the disk-integrated and local photosphere. Such profiles can be used for comparison with stellar models or to generate binary CCF masks specific to a given star. The workflow is modular, so that new methods can be added in a straightforward way. Because the concepts behind the processing and analysis methods do not depend on a given instrument, their use can be extended to other ground-based or space-borne, high-resolution spectrographs as long as the relevant correction and formatting methods are added. In this work, we present the application of the ANTARESS workflow to planet-occulted stellar spectra. In a companion paper, we will describe the subsequent methods aimed at extracting and analyzing series of atmospheric spectra uncontaminated by the star to sample the planetary limbs spatially and dynamically.

Acknowledgements

We highlight our appreciation of the referee’s work, considering the depth in which they read this comprehensive manuscript and the attention they gave to understanding its contents. We thank G. Fu and N. Nikolov for providing the low-resolution transmission spectrum of WASP-76b, K. Lind and A. Amarsi for their help with stellar NLTE grids for sodium, J. Hoeijmakers for discussions about transmission spectroscopy, A. Deline for inspiring discussions about wiggles, D. Luc for comments about coding, Y. Zhao for his help with GSL. 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. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project Spice Dune, grant agreement No. 947634; project SCORE, grant agreement No 851555). X.D acknowledges the support from the Swiss National Science Foundation under the grant SPECTRE (No. 200021_215200). 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 the Fonds de Recherche du Québec – Nature et Technologies (FRQNT). This work was funded by the Institut Trottier de Recherche sur les Exoplanètes (iREx). H.M.C. acknowledges funding from a UKRI Future Leader Fellowship, grant number MR/S035214/1. This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. ANTARESS makes use of the following packages: ASTROPY, Astropy Collaboration (2022); BATMAN, Kreidberg (2015); EMCEE, Foreman-Mackey et al. (2013); LMFIT, Newville et al. (2016); MATPLOTLIB, Hunter (2007); NUMPY, Harris et al. (2020); PANDAS, pandas development team (2020); PATHOS, McKerns et al. (2012); PYASTRONOMY, Czesla et al. (2019); PYSME, Wehrhahn et al. (2023); SCIPY, Virtanen et al. (2020).

Appendix A Resampling, error propagation, and fitting

A.1 Resampling and error propagation

We developed for ANTARESS an open-source C library with python bindings, called BINDENSITY⁹. It implements two resampling algorithms as well as basic operations on spectra, with error propagation using banded covariance matrices. We defined a linear and a cubic interpolator of the cumulative flux, which corresponds to assuming a uniform and a quadratic density over each pixel, respectively. In both cases, the cumulative flux is assumed to be continuous, and the integral of the flux on each original pixel is preserved. The linear interpolator can be used for fast, preliminary analyses, but the cubic interpolator should be preferred to limit blurring sharp spectral features. We do not propose higher order interpolators to prevent the risk of overfitting and creating spurious signals in the resampled spectra.

In the linear case, this rule is sufficient to determine the two coefficients of the interpolator, since it defines the value of the cumulative at each edge of the original pixels. With this linear interpolator, only new pixels that overlap the same original pixel will be correlated (see Fig. A.1). If the size of the new pixels is similar to the size of the original pixels, it typically generates correlation only between two consecutive pixels, and the covariance matrix will be tridiagonal (the main diagonal and the first sub/super-diagonal are non-zero). If several successive linear interpolations are performed, the covariance matrix’s bandwidth (number of non-zero sub/super diagonals) keeps increasing at each step. This bandwidth growth should be avoided because of the increased computational cost, but also because it is associated to a loss of resolution. This is why we avoid as much as possible to perform successive changes of rest frames, but instead apply spectral shifts directly to the spectra.

In the cubic case, one needs to specify two additional constraints to determine the four coefficients of the interpolating polynomial on each bin. A widely spread method is to use cubic splines, which require the interpolator to be C² (twice continuously differentiable). The polynomial coefficients can then be obtained by solving a tridiagonal system of equations, which can be implemented with dedicated efficient algorithms, but which introduces correlations between all pixels in the spectrum. Moreover, assuming the cumulative flux to be C² means that the flux density is continuously differentiable. Such an hypothesis might not be physically justified if spectra present breaks, in particular in saturated absorption lines. We thus decided to relax the C² condition and to require the derivative at each edge between two pixels to be equal to the average density of the two pixels. Our cubic interpolator is thus only C¹ in cumulative flux, which means that the flux density is continuous but not differentiable. Using these modified constraints, solving for the polynomial coefficients is straightforward, and the interpolation introduces correlations between new pixels that overlap the same original pixel or two consecutive original pixels. If the sizes of the new and original pixels are similar, the resulting covariance matrix would typically have two sub/super-diagonals (see Fig. A.1).

A.2 Fitting

Model fitting in the ANTARESS methods is performed using the log-probability function of a set of parameters θ: $$\ln p(\theta \mid d,M)=\ln \varpi (\theta )+\ln \mathcal{L}(\theta ),$$(A.1)

where m represents the family of models fitted to the data d, ϖ(θ) = p(θ|m) the prior probability of the model parameters, and $\mathcal{L}(\theta )=p(d\mid m,\theta )$ the likelihood. Under the assumption of Gaussian noise we have: $$\ln \mathcal{L}(\theta )=-\ln (\sqrt{\det (2\pi \mathrm{\Sigma })})-\frac{1}{2}{r}^t{\mathrm{\Sigma }}^{-1}r,$$(A.2)

where r = d – m(θ) is the vector of residuals and ∑ the noise covariance matrix.

We use the Cholesky decomposition of the banded covariance matrix C implemented in the workflow, so that: $$\ln \mathcal{L}(\theta )=-{\displaystyle \sum _n}\ln (\sqrt{2\pi }L[n,n])-\frac{1}{2}{u}^{t}u.$$(A.3)

Here, L is the lower triangular matrix satisfying C = LL^t, and u = L⁻¹ r. It is banded with the same width as C, so that the cost of all calculations scales with the number of spectral pixels.

If we assume that the data is dominated by white noise and correlations can be neglected, then ∑ reduces to the variance σ² along its diagonal and $$\ln \mathcal{L}(\theta )=-{\displaystyle \sum _n}\ln (\sqrt{2\pi }\sigma (n))-\frac{1}{2}{\displaystyle \sum _n}{\left(\frac{r(n)}{\sigma (n)}\right)}^2.$$(A.4)

The MCMC approach implemented in the ANTARESS pipeline samples the posterior probability distributions (PDF) of the free parameters describing m using EMCEE (Foreman-Mackey et al. 2013). We propose standard priors to set on the parameters (Uniform and Gaussian) as well as original priors: “Gaussian Halves” (defined as two half-Gaussian profiles with the same center and different widths, or “Complex” (defined through a custom function that uses combinations of the model parameters). The number of walkers and the burn-in phase is adjusted based on the degrees of freedom of the considered problem and the convergence of the chains. Best-fit values for the parameters are set to the median of their PDFs, and their 1σ uncertainty ranges are defined using highest density intervals to better account for multimodal and asymmetric PDFs.

Appendix B Description of planetary systems

B.1 Coordinates

We describe here the orthonormal reference frames centered on the star that are used to define the planet coordinates and introduce various orbital properties of interest. The first frame has its x axis as the major axis of the orbit oriented toward the periapsis and its y axis obtained by rotating x within the plane along the orbital motion: $$\begin{array}{rl}& x_{\text{orb}}(t)=\frac{a}{R_{\star }}(\cos (E(t))-e),\\ & y_{\text{orb}}(t)=\frac{a}{R_{\star }}\sqrt{1-e^2}\sin (E(t)).\end{array}$$(B.1)

Here, E is the eccentric anomaly, derived numerically at a given orbital phase (or time t). The corresponding orbital velocity coordinates are: $$\begin{array}{rl}{x}_{\text{orb}}^{\mathrm{\prime }}(t)& =-\frac{2\pi }{P}{\left(\frac{a}{R_{\star }}\right)}^2\frac{1}{D_{\mathrm{p}}(t)}\sin (E(t)),\\ y_{\text{orb}}^{\mathrm{\prime }}(t)& =\frac{2\pi }{P}{\left(\frac{a}{R_{\star }}\right)}^2\frac{1}{D_{\mathrm{p}}(t)}\sqrt{1-e^2}\cos (E(t)),\\ \text{\hspace{0.17em}with\hspace{0.17em}}{D}_{\mathrm{p}}(t)& =\frac{a}{R_{\star }}(1-e\cos (E(t))).\end{array}$$(B.2)

Rotating the orbital coordinate system around the orbital plane normal by the argument of periastron ω, so that x becomes the node line, yields: $$\begin{array}{rl}& x_{\mathrm{L}\mathrm{O}\mathrm{S}}(t)=-x_{\mathrm{o}\mathrm{r}\mathrm{b}}(t)\cos (\omega )+y_{\mathrm{o}\mathrm{r}\mathrm{b}}(t)\sin (\omega ),\\ & y_{\mathrm{L}\mathrm{O}\mathrm{S}}(t)=x_{\mathrm{o}\mathrm{r}\mathrm{b}}(t)\sin (\omega )+y_{\text{orb\hspace{0.17em}}}(t)\cos (\omega ).\end{array}$$(B.3)

Rotating the coordinate system around the node line by the orbital inclination i_p then yields the sky-projected “orbital” reference frame in which y is along the sky-projected orbital normal and z is along the LOS: $$\begin{array}{rl}& x_{\text{sky}}(t)=x_{\mathrm{L}\mathrm{O}\mathrm{S}}(t),\\ & y_{\text{sky}}(t)=-y_{\mathrm{L}\mathrm{O}\mathrm{S}}(t)\cos \left(i_p\right),\\ & z_{\text{sky}}(t)=y_{\mathrm{L}\mathrm{O}\mathrm{S}}(t)\sin \left(i_p\right).\end{array}$$(B.4)

Rotating the latter reference frame around the LOS by the sky-projected obliquity λ of the planet yields the sky-projected “stellar” frame, where x is the node line of the stellar equator and y the sky-projected stellar spin: $$\begin{array}{rl}& x_{\star \mathrm{s}\mathrm{k}\mathrm{y}}(t)=x_{\text{sky}}(t)\cos (\lambda )-y_{\text{sky}}(t)\sin (\lambda ),\\ & y_{\star \mathrm{s}\mathrm{k}\mathrm{y}}(t)=x_{\text{sky}}(t)\sin (\lambda )+y_{\text{sky}}(t)\cos (\lambda ),\\ & z_{\star \mathrm{s}\mathrm{k}\mathrm{y}}(t)=z_{\text{sky}}(t).\end{array}$$(B.5)

Here, we caution that some equations in the workflow (e.g, that of the rv of planet-occulted regions) require the coordinates of the sky-projection of the planet onto the stellar surface. In that case, z_*sky is calculated as $\sqrt{1-x_{\star \mathrm{s}\mathrm{k}\mathrm{y}}^2-x_{\star \mathrm{s}\mathrm{k}\mathrm{y}}^2}$ for a spherical star and by solving the quadratic describing the photosphere for an oblate star. A final rotation around the star node line by i_*, inclination of the stellar spin with respect to the LOS, yields the position of the planet or occulted region in the star rest frame, in which x now follows the stellar equator, and y is the actual stellar spin: $$\begin{array}{rl}& x_{\star }(t)=x_{\star \mathrm{s}\mathrm{k}\mathrm{y}}(t),\\ & y_{\star }(t)=y_{\star \mathrm{s}\mathrm{k}\mathrm{y}}(t)\sin \left(i_{\star }\right)+z_{\star \mathrm{s}\mathrm{k}\mathrm{y}}(t)\cos \left(i_{\star }\right),\\ & z_{\star }(t)=-y_{\star \mathrm{s}\mathrm{k}\mathrm{y}}(t)\cos \left(i_{\star }\right)+z_{\star \mathrm{s}\mathrm{k}\mathrm{y}}(t)\sin \left(i_{\star }\right).\end{array}$$(B.6)

Another useful coordinate is the Keplerian rv of the star, relative to the barycenter of the planetary system: $${\mathrm{r}\mathrm{v}}_{\star /{\mathrm{B}}_{\star }}(\mathrm{t})={\displaystyle \sum _{\text{body\hspace{0.17em}}\mathrm{p}}}{\mathrm{K}}_{\mathrm{p}}[\cos ({\theta }_{\mathrm{p}}(\mathrm{t})+{\omega }_{\mathrm{p}})+\mathrm{e}\cos \left({\omega }_{\mathrm{p}}\right)],$$(B.7)

where K_p is the Keplerian rv semi-amplitude induced by a given body orbiting the star and θ_p the true anomaly of its orbit.

B.2 Planet-occulted profiles

We describe here the detailed calculation of line profiles. Without oversampling of the transit chord, properties of the planet-occulted region are calculated for a single position of the planetary grid corresponding to the center of the chord. With oversampling, the planetary grid is positioned at regular intervals k along the exposure chord. The disk of planet p is discretized at each position k, calculating the relevant stellar surface properties x_i in each cell $i_{k,\overline{\lambda }}$ of the grid. Only cells that fall within the sky-projected photosphere, and outside of the planetary disks already processed, are considered. We note the chromatic dependance of the planet-occulted cells and properties, linked to broadband variations of the stellar intensity and planet radius.

Planet-occulted line profiles can then be calculated at two precision levels (Fig. 3):

• Low: properties of the region occulted during the time window $\overline{t}$ are approximated as the brigthness-weighted average of all x_i: $$x(\overline{\lambda },\overline{t})=\frac{{\displaystyle {\sum }_{k,p}{\sum }_{i_{k,\overline{\lambda }}}}{x}_i{I}_i(\overline{\lambda })S_i}{{\displaystyle {\sum }_{k,p}{\sum }_{i_{k,\overline{\lambda }}}}{I}_i(\overline{\lambda })S_i},$$(B.8)

  where $I_i(\overline{\lambda })$ and S_i designate a cell intensity level and area. The intrinsic profile associated with the occulted region is then defined as $I(x(\overline{\lambda },\overline{t}),\lambda )$, with I one of the models described in Sect. 3.3.1. Because the photospheric rv field may vary substantially over the occulted region, I is broadened with a rotational kernel defined by the range of velocities covered by the exposure chord. The intrinsic profile is then continuum-scaled using the chromatic planet-occulted flux calculated with Eq. B.8 for x_i = 1, and Doppler-shifted using the chromatic planet-occulted $rv(\overline{\lambda },\overline{t})$.

• High: the line profile of each oversampled region is calculated as $F_k(\lambda )={\displaystyle \sum _{i_k}}{F}_i(\lambda )$, where F_i(λ) = I(x_i, λ)S_i is the local profile from a given cell. The exposure profile is then calculated as the average of the F_k over all regions.

Chromatic planet-occulted profiles can only be defined in low precision, as the local profiles summed in high precision already depend on wavelength.

Appendix C Spectral correction models

Here, we describe the models used by the CORRECTIONS methods of the ANTARESS workflow.

C.1 Spectral flux calibration model

Because DRS conversion and correction tables are not necessarily available, we estimate the detector flux calibration profile directly from the input spectra. Assuming that g_cal varies slowly with wavelength, Eq. 3 yields: $$g_{\mathrm{c}\mathrm{a}\mathrm{l}}(\overline{\lambda },E,t)=\frac{{\displaystyle {\sum }_{\overline{\lambda }}}{\sigma }^{\text{meas\hspace{0.17em}}}(\lambda ,E,t{)}^2}{{\displaystyle {\sum }_{\overline{\lambda }}}{F}^{\text{meas\hspace{0.17em}}}(\lambda ,E,t)}.$$(C.1)

Values of $g_{\mathrm{cal}}(\overline{\lambda },E,t)$ measured in each spectral order are fitted with a piece-wise model made of three polynomials continuous in value and derivative. They cover the central part of the order, as well as its blue and red edges where the count level decreases sharply, ensuring that the shape of g_cal is properly captured. Outliers, typically in regions of deep telluric or stellar lines, are automatically excluded from the fit. The derived profiles show slight variations in time, especially in regions of low S/N. Our final estimate of the flux calibration profile is thus <g>_cal, the median of the calibration models fitted to all exposures of a processed instrument.

C.2 Telluric spectra model

Properties of the Earth atmospheric model are derived for each molecule by fitting its output telluric spectrum to the spectrum measured in each exposure, using Levenberg-Marquard minimization. Direct measurements of the ambient temperature, airmass, and integrated water vapor at zenith are available at the facilities of ESPRESSO, HARPS, HARPS-N, CARMENES, and NIRPS. They are used to fix the model temperature, and to initialize the fitted pressure and column density of the species along the LOS, for spectra obtained with these instruments. Uniform priors can further be set on these fitted properties, which can be useful in low-S/N datasets where some exposures are poorly adjusted.

First, the spectral grid of the observed spectrum is temporarily shifted back to the Earth rest frame, and if necessary converted into vacuum wavelengths. Then, a high-resolution telluric model is calculated in wavenumber space at each step of the minimization using a subset of deep telluric absorption lines. The wavelength, width, and intensity of the lines are adjusted depending on the fitted atmospheric properties, and used to compute the line profiles (Lorenzian for H₂O, Voigt for other species; see Allart et al. 2022). The optical depth spectrum of the species is obtained through multiplication of its column density by the co-added absorption line profiles. The resulting telluric absorption spectrum is convolved with the instrumental response and sampled to the instrumental wavelength grid. 2D maps providing the width of a Gaussian resolution kernel as a function of detector pixel and spectral order were already available for the different modes and epochs of ESPRESSO and NIRPS. We derived maps for HARPS and HARPS-N following the same approach as in Allart et al. (2022), applying 2D polynomial fits over unresolved thorium lines measured across the instrument detector. For other instruments, a constant resolution is used. The comparison between model and data is performed over a CCF calculated with a binary mask composed of the selected telluric lines, with weights unity.

An interest of this approach is the possibility to determine whether the correction of a specific molecule is warranted, by evaluating the amplitude of its telluric CCF in the raw and corrected data (Fig. 6). We improve upon the original ATC by using the ANTARESS method to compute CCFs (Sect. 7.5), applying accurate spectral resampling and error propagation (Sect. 2.3), and deriving more realistic uncertainties through the use of the covariance matrix (Sect. 2.3, Fig. 6).

Once best-fit atmospheric properties are found for all molecules, we co-add their high-resolution optical depth spectra using the full HITRAN linelist (which include micro-tellurics), calculate the resulting telluric absorption spectrum, convolve it with the instrumental response, and resample it over the spectral grid of the processed exposure. Corrected spectra are obtained through division with this final telluric absorption spectrum. We set a threshold on telluric line depth above which the telluric-absorbed flux is considered too low and noisy to be retrieved. The threshold is typically set to 90% (i.e. 10% of the flux being transmitted) but is adjustable depending on the dataset. Flagged pixels are masked.

C.3 Flux balance model

We describe here the approach used to model flux balance variations (Sect. 5.3).

We first define the spectra F_ref that are used as reference for the flux balance of the star. A selection of spectra in each epoch are aligned in the star rest frame (shifting them by the stellar Keplerian and systemic rv motion), resampled on a common spectral grid, and normalized in each order to their average flux level over the selection. Reference spectra in each epoch are then calculated as the median of the processed spectra, to avoid distortions induced by spurious features (e.g., cosmics, Sect. 5.4). Reference spectra for each instrument are calculated as the mean of their epoch-specific median spectra, naturally giving more weight to epochs with higher flux levels. If a single instrument is processed, using this instrument reference to set the final flux balance for all its spectra limits the amplitude of the correction and keeps corrected spectra closer to their original balance. If several instruments are processed, their common flux balance can be set to one of the instrument-specific references, or to measured/theoretical stellar spectra independent of the processed datasets.

We then determine the color balance variations for each exposure. First, the median reference spectrum is shifted to the rest frame of the exposure spectrum and resampled over its spectral grid. Then, the reference and exposure spectra are converted back from extracted F to raw N photoelectron counts (Sect. 5.1) and binned at a lower resolution in each order. The ratio R^meas between the binned reference and exposure spectra is fitted with a smoothing spline or polynomial R^theo. A preliminary fit is performed to automatically identify and exclude outliers, and the final fitted model (Fig. 7) is an estimate of: $$\begin{array}{rl}{R}^{\text{meas\hspace{0.17em}}}({\overline{\nu }}_{{\mathrm{B}}_{\odot }},E,t)& =\frac{N^{\text{meas\hspace{0.17em}}}({\overline{\nu }}_{{\mathrm{B}}_{\odot }},E,t)}{N_{\mathrm{r}\mathrm{e}\mathrm{f}}({\overline{\nu }}_{{\mathrm{B}}_{\odot }},E)}\\ & =\frac{{\delta }_{\mathrm{p}}({\overline{\nu }}_{{\mathrm{B}}_{\odot }},E,t)c({\overline{\nu }}_{{\mathrm{B}}_{\odot }},E,t)}{C_{\mathrm{r}\mathrm{e}\mathrm{f}}({\overline{\nu }}_{{\mathrm{B}}_{\odot }},E)}.\end{array}$$(C.2)

Here, we assume that the reference spectrum and the true stellar spectrum deviate at low resolution, so that C_ref and N_* remain separated when binning N_ref = C_refN_*. Light frequency ν is used as fit coordinate to stretch the data and its variations over a more regular grid. Finally, we extrapolate $R^{\text{theo}}({\overline{v}}_{{\mathrm{B}}_{\odot }},E,t)$ over the native spectral grid of the exposure to define its correction. The same approach is used when determining the color balance variations between the median reference and the final reference for each epoch (Sect. 5.3). We note that this final reference is not used directly to correct flux balance variations within each epoch because we found it introduces high-frequency variations in the correction of individual spectra.

When correcting for flux balance variations we recommend that orders and ranges that are too contaminated by telluric absorption (typically toward longer wavelengths in optical spectra and between the J and H bands in NIR spectra), or where noise is too large to fit properly the color variations (typically the bluest orders in low-S/R optical data), are excluded. The bin size, as well as the smoothing factor or polynomial degree, must be adapted for each dataset to minimize the broadband dispersion of F^corr/F_ref while preventing the overcorrection of spectra and the introduction of artificial variations. Bin size must be smaller than the typical range of variations of the color balance c and planetary contribution δ_p but large enough to smooth out high-frequency features, especially those caused by variations in narrow stellar and planetary lines. Variations at the level of each order linked to the physical size of the detector in wavelength space (Sect. 5.3) must be sampled with a smaller binwidth but using a low-order polynomial to only capture order-wide variations.

C.4 Cosmics

We detail here the processing of the datasets for the cosmics correction method. Each epoch is treated independently. The number of adjacent exposures used for comparison can be adjusted depending on the dataset, with one half selected before the processed exposure and the other half after. In cases where the planet induces a strong RM effect, we advise choosing a small number of adjacent exposures so that their spectra remain similar (ie, with similar RM distortions). Spectra used for comparison are temporarily aligned in the stellar rest frame using the Keplerian rv model for the system (Sect. 6.2), rv provided by the instrument DRS, or delta-rv measured via cross-correlation of the adjacent exposure spectra with the processed one. In general, we recommend using the Keplerian model since it is not sensitive to the photonic rv dispersion, but the other approaches can be useful in high-S/R datasets where the rv dispersion around the Keplerian motion is caused by stellar activity (e.g., granulation). The aligned spectra of adjacent exposures are resampled on the spectral grid of the processed exposure at time t, and set to the same global flux level. We can then calculate the mean F_comp(λ, E) and standard deviation σ_Fcomp(λ_*, E) for each bin, over all adjacent exposures, and search for cosmics as described in Sect. 5.4.

C.5 Persistent peaks

To identify persistent peaks, spectra are first temporarily shifted back to the Earth rest frame so that detector features and telluric lines are aligned. We calculate a master stellar spectrum as the median of all spectra in an epoch to estimate the stellar continuum (Sect. 7.1) and calculate the residuals between the spectrum of each exposure and this continuum. Then we identify bins where residuals deviate by more than a threshold times their flux error, and flag all surrounding bins falling within a chosen window (with typical width that of the peaks to be excluded). Finally, a bin is masked if it is flagged in all exposures or in a chosen number of consecutive exposures.

C.6 Wiggle model

C.6.1 Description

From the analysis of several datasets besides those in the present study, we determined that the wiggles are best described as the sum of multiple sinusoidal components k (typically two): $$W(\nu ,t)=1+{\displaystyle \sum _k}{A}_k(\nu ,t)\sin (2\pi \int (F_k(\nu ,t)d\nu )-{\mathrm{\Phi }}_k(t)).$$(C.3)

The amplitude A_k and frequency¹⁰ F_k of each component show fewer variations as a function of light frequency ν = c/λ than wavelength, and are described as: $$\begin{array}{rl}& A_k(\nu ,t)={\displaystyle \sum _{i=0}^{da,k}}{a}_{\text{chrom\hspace{0.17em}},\mathrm{k},\mathrm{i}}(t){(\nu -{\nu }_{\mathrm{r}\mathrm{e}\mathrm{f}})}^i,\\ & F_k(\nu ,t)={\displaystyle \sum _{i=0}^{df,k}}{f}_{\text{chrom\hspace{0.17em}},\mathrm{k},\mathrm{i}}(t){(\nu -{\nu }_{\mathrm{r}\mathrm{e}\mathrm{f}})}^i,\end{array}$$(C.4)

where the polynomials are defined relative to v_ref, corresponding to 6000 Å, to decrease leverage in the fits. The chromatic amplitude and frequency coefficients a_chrom,k,i(t) and f_chrom,k,i(t), as well as the phase reference Φ_k(t), are expressed as a linear combination of the cartesian pointing coordinates g({α_point,k}, t) = $\begin{array}{cc}{\alpha }_{\mathbf{u}/\text{west}}+{\alpha }_{\mathrm{x}}x(t)+{\alpha }_{\mathrm{y}/\text{west}}y(t)+{\alpha }_{\mathrm{z}/\text{west}}z(t)& \text{if\hspace{0.17em}}\theta (\mathrm{t})>{\theta }_{\text{mer}},\\ {\alpha }_{\mathrm{u}/\text{east}}+{\alpha }_{\mathrm{x}}x(t)+{\alpha }_{\mathrm{y}/\text{east}}y(t)+{\alpha }_{\mathrm{z}/\text{east}}z(t)& \text{if\hspace{0.17em}}\theta (\mathrm{t})\le {\theta }_{\text{mer}},\\ \text{with\hspace{0.17em}}x(t)=\sin (\theta (t)),& \\ y(t)=\cos (\theta (t)),& \\ z(t)=\sin (\varphi (t)),& \end{array}$(C.5)

where θ and ϕ are the azimuth and altitude of the telescope. We found that the model coefficients change when the LOS crosses the meridian (θ_mer = 180°), but that the model remains continuous in value (α_y/west = α_y/east + sin(ϕ_mer)(α_z/west – α_z/east) – (α_u/west – α_u/east)) and derivative (α_x/west = α_x/east = α_x). The full model is thus controlled by a set of pointing parameters {a_point,k,i}, {f_point,k,i}, and {Φ_point,k,i} for each component k.

C.6.2 Fitting procedure

Each epoch is processed independently. First, 2D spectra are aligned in the star rest frame using a Keplerian RV model (Sect. 6.2), so that stellar lines are removed and the wiggles isolated when calculating transmission spectra. A master spectrum is computed for the epoch, using all exposures to smooth out as much as possible the wiggles. The master and spectra are converted back from extracted to raw photoelectron counts (Sect. 7.1), and downsampled over a common spectral grid. The downsampling binwidth is large enough to remove possible correlations between bins but small enough to sample all wiggle components. In-transit planetary lines are ignored in our analysis, considering that they are narrower than the period of the wiggles and typically smaller in amplitude. A 1D transmission spectrum is computed for each exposure by dividing their spectrum with the master in each spectral order, and grouping them together. Transmission spectra for all exposures are shifted into the Earth rest frame, which we expect traces more directly the wiggles, and run through the following steps:

1. Screening: required to exclude spectral ranges of poor quality (e.g., from telluric contamination) or unaffected by the wiggles (bluest orders) from the analysis. A periodogram cumulated over all exposures allows identifying which wiggle components are detectable (Fig. C.1). The dominant and weak components always contain most of the wiggle power, but additional components may be present and included in the model if relevant (although its power does not justify a correction we identified a third component with frequency ~ 156 × 10¹⁰ s in the HD 209458 datasets, which matches the high-frequency signal noted by Casasayas-Barris et al. 2021).

2. Chromatic sampling: we run a sliding window over each transmission spectrum, identify the strongest peak in a periodogram calculated at each window position, and fit a sine to the window spectrum using the peak frequency (Fig. C.2). Because the wiggles can be approximated by sines with constant frequency over narrow bands, this step allows sampling the frequency F_k(v) and amplitude A_k(v) of a given component in each exposure. The sampling windows must be broad to contain several cycles of a wiggle component, so that we make the windows overlap between successive positions to oversample the measurements. Once a component is processed, its piecewise model built from all windows is used to correct temporarily the transmission spectrum so that the next component can be sampled.

3. Chromatic analysis: the frequency and amplitude of each component sampled in the previous step are fitted with polynomial models of ν (Fig. C.3). This allows determining their degree and deriving guess values for the chromatic coefficients a_chrom,k,i and f_chrom,k,i in each sampled exposure. In most datasets that we analyzed, the amplitude and frequency could be described as a quadratic and linear function of ν, respectively decreasing and increasing toward shorter wavelengths. For high-S/R wiggles, which can be characterized over the full ESPRESSO range, higher degrees may be necessary to describe the wiggle variations in the bluest orders.

4. Exposure fit: the spectral wiggle model W(v) is initialized using the results from the previous step and fitted to all exposures individually (Fig. 11). This step provides accurate estimates of the a_chrom,k,i(t), f_chrom,k,i(t), and Φ_k(t) coefficients.

5. Pointing analysis: the time series of chromatic amplitude and frequency coefficients, as well as the phase, derived independently from the fit to each exposure are fitted with the pointing coordinate function (Eq. C.6.1), to evaluate the pointing parameters {a_point,k,i}, {f_point,k,i}, and {Φ_point,k,i} (Fig. C.4).

6. Global fit: the full spectro-temporal wiggle model W(v, t) is initialized using the results from the previous step and fitted to all exposures together. Given the large number of free parameters and the complexity of a model made of cumulated sines, guess values close to the best-fit are essential to the convergence of the fit.

Appendix D Detailed processing methods

D.1 Stellar line detrending

We describe here our approach to correct individual stellar lines (i.e., an isolated spectral line or a cross-correlated line) or a full spectrum for systematic trends characterized using the out-of-transit ross-correlated line properties:

• rv correction: profiles are shifted using the rv detrending model, following the same approach as in Sect. 6.2.

• Contrast correction: profiles are temporarily normalized and offset to a null continuum. A “vertical stretching” is then applied by dividing the flux with the contrast detrending model, before returning profiles to their original continuum. Normalization is performed using the average continuum flux for single line profiles, or the stellar spectral continuum (Sect. 7.1) for full spectra.

• FWHM correction: profiles are temporarily aligned on a null velocity (spectral profiles are first passed into RV space relative to the line wavelength). A “horizontal stretching” is then applied by dividing the aligned rv grid with the FWHM detrending model, before profiles are shifted back to their original frame. This correction can only be performed on single line profiles. FWHM and contrast models are normalized over each epoch to maintain the properties of the corrected profiles to their original mean value.

D.2 Broadband flux scaling

The F^meas spectra posterior to the FLUX BALANCE correction have the same low-frequency profile as the chosen reference F_ref: $$\begin{array}{rl}{F}^{\text{meas\hspace{0.17em}}}({\overline{\lambda }}_{\star },E,t)& =F_{\text{ref\hspace{0.17em}}}({\overline{\lambda }}_{\star },E)R^{\text{norm\hspace{0.17em}}}(E,t)\\ & =C_{\text{ref\hspace{0.17em}}}({\overline{\lambda }}_{\star },E)F_{\star }({\overline{\lambda }}_{\star },E)R^{\text{norm\hspace{0.17em}}}(E,t).\end{array}$$(D.1)

First, we divide all spectra by their corresponding exposure duration, converting them in the flux density units to which F will correspond to hereafter. Then we correct spectra for R^norm (ratio between the total flux of the spectrum and that of the reference, Eq. 6) through multiplication by $$C_{\mathrm{F}}(E,t)=\frac{{\mathrm{T}\mathrm{F}}_{\mathrm{r}\mathrm{e}\mathrm{f}}}{{\displaystyle {\sum }_{{\lambda }_{\star }}}{F}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}({\lambda }_{\star },E,t)d{\lambda }_{\star }}.$$(D.2)

The integrated flux of the resulting profiles is thus normalized to TF_ref, which can either be set to the median flux over the epoch < t_k,∑_λ* F^meas (λ_*, E, t)dλ_* >, or to F_ref ∑_λ* dλ_* with F_ref an arbitrary flux density. For consistency, the integration is performed over the same spectral range used to calculate R^norm. We note that regions affected by the planet should not be excluded from this spectral range, since measured light curves (see below) are integrated over bands that include those regions.

Then we re-inject the broadband contributions of the star and planet via a correction c(λ_*, E, t), so that: $$c({\lambda }_{\star },E,t)F^{\text{meas\hspace{0.17em}}}({\lambda }_{\star },E,t)=F({\lambda }_{\star },E,t).$$(D.3)

Here, the time series of F spectra have the correct relative flux level. To define the final flux scaling we can write the above equation at low frequency: $$c({\overline{\lambda }}_{\star },E,t)C_{\mathrm{r}\mathrm{e}\mathrm{f}}({\overline{\lambda }}_{\star },E){\displaystyle \sum _{{\lambda }_{\star }}}{F}_{\star }({\lambda }_{\star },E)d{\lambda }_{\star }={\displaystyle \sum _{{\lambda }_{\star }}}F({\lambda }_{\star },E,t)d{\lambda }_{\star }.$$(D.4)

We remind that the factor C_ref represents the unknown deviation from the true stellar spectrum. Because this deviation profile has low-frequency variations and is affecting all exposures in an epoch in the same way, it only matters when extracting planetary emission using absolute fluxes or if the stellar and reference spectra deviate differently between epochs (see Sect. 5.3). In general, we can thus ignore C_ref in the final correction lc, which is defined as: $$lc({\overline{\lambda }}_{\star },E,t)=\frac{{\displaystyle {\sum }_{{\lambda }_{\star }}}F({\lambda }_{\star },E,t)d{\lambda }_{\star }}{{\displaystyle {\sum }_{{\lambda }_{\star }}}{F}_{\star }({\lambda }_{\star },E)d{\lambda }_{\star }}.$$(D.5)

The scaling light curve $lc({\overline{\lambda }}_{\star },E,t)$ can be defined chromatically, namely, over a set of bands at ${\overline{\lambda }}_{\star }$, to account for broadband variations in the planetary absorption/emission and stellar emission (Fig. 13). We propose three methods to calculate lc:

• import: we import a measured or theoretical light curve into the ANTARESS workflow. This approach allows accounting for time-variable features, possibly specific to a given epoch, such as spot crossings, planetary thermal emission, or reflection of starlight.

• BATMAN: we use this package from Kreidberg (2015) to compute a transit light curve with simple shape, controlled by the transit depth of a single planet, one of the available limb-darkening laws (uniform, linear, quadratic, square-root, logarithmic, exponential, power2, nonlinear) and its associated coefficients, and time resolution.

• numerical: we use ANTARESS stellar surface model (Sect. 3.3). The light curve is calculated by subtracting from the disk-integrated emission the flux from regions occulted by the planet(s) in each exposure, calculated following the approach described in Sect. 3.3.2. This method allows accounting for planets transiting simultaneously and for complex stellar photospheres (e.g., oblate gravity-darkened stars, or spotted stars).

Those three methods allow defining light curves specific to a given epoch, to account for changes in the planetary or photo-spheric properties (in particular variations in the stellar emission due to unocculted spots). In all cases, it is important to define the light curve at high temporal resolution so that it can be averaged within the time window of each exposure and avoid blurring.

We note that algorithms like BATMAN or ANTARESS’s numerical approach actually return ${\delta }_{\mathrm{p}}({\overline{\lambda }}_{\star },E,t)$ (Sect. 2.1), which is calculated assuming a constant stellar emission and planetary absorption within a band. This quantity is not strictly equivalent to the integrated flux ratio between each exposure and the star in Eq. D.5, since F differs from F* at high resolution (due to the planet-occulted and atmospheric-filtered stellar lines, Eq. 1) and F* thus does not simplify in this ratio. This does not matter, however, as long as the theoretical light curves match the ratio.

Finally, we interpolate the chromatic $lc({\overline{\lambda }}_{\star },E,t)$ over the high-resolution spectral grid of each exposure to calculate the rescaled spectra as: $$\begin{array}{rl}{F}^{\mathrm{s}\mathrm{c}}({\lambda }_{\star },E,t)& =F^{\text{meas\hspace{0.17em}}}({\lambda }_{\star },E,t)C_{\mathrm{F}}(E,t)lc({\lambda }_{\star },E,t)\\ & =F({\lambda }_{\star },E,t)C_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\text{norm\hspace{0.17em}}}({\overline{\lambda }}_{\star },E)\\ \text{\hspace{0.17em}where\hspace{0.17em}}{\mathrm{C}}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}({\overline{\lambda }}_{\star },\mathrm{E})& =\frac{T{F}_{\mathrm{r}\mathrm{e}\mathrm{f}}}{{\displaystyle {\sum }_{{\lambda }_{\star }}}{F}_{\mathrm{r}\mathrm{e}\mathrm{f}}({\lambda }_{\star },E)d{\lambda }_{\star }},\end{array}$$(D.6)

D.3 Differential profile extraction

We detail here the procedure to retrieve differential profiles. For each exposure at time t we build a master disk-integrated profile of the unocculted star, $F_{\star }^t$, as the weighted average (Sect. 7.3) of a selection of out-of-transit F^sc profiles representative of the star during transit. Instead of resampling the out-of-transit profiles on a common grid, averaging them into a single master, and resampling this master on the individual spectral grid of each exposure, we directly resample the out-of-transit profiles on the spectral grid of the processed exposure before they are averaged. The master is thus built from profiles that are only resampled once in our approach, preventing blurring and the introduction of correlated noise (Sect. 2.3).

Spectral ranges contaminated by the planetary atmosphere (i.e., where F_p and $S_i^{p[thin]}$ are not null, Eq. 10) are identified as described in Sect. 3.2 and masked so that the master represents the pure stellar emission (see Sect. 7.3). Future versions of the ANTARESS workflow could account for possible variations in the disk-integrated profile of the unocculted star, for example in the presence of spots moving over the duration of an epoch, so that $F_{\star }^t$ remains equivalent to the disk-integrated profile occulted at time t.

D.4 Intrinsic profile extraction

We detail here the procedure to retrieve intrinsic profiles, which are derived from differential profiles as: $$\begin{array}{rl}{F}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}}({\lambda }_{\star },E,t)& =\frac{F^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}({\lambda }_{\star },E,t)}{1-lc({\lambda }_{\star },E,t)}\\ & =C_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}({\overline{\lambda }}_{\star },E)\frac{{\displaystyle {\sum }_{i\in i_{\text{occ\hspace{0.17em}}}}}{f}_i({\lambda }_{\star },E)S_i^{\mathrm{p}[\text{\hspace{0.17em}thick\hspace{0.17em}}]}({\overline{\lambda }}_{\star },E)}{1-lc({\lambda }_{\star },E,t)}.\end{array}$$(D.7)

Here, the spectral ranges absorbed by the optically thin planetary atmosphere have been masked, unless intrinsic profiles are spectral and converted into CCFs later on in the workflow. Indeed, if the planetary atmosphere absorbs at transitions present in the CCF mask used to cross-correlate intrinsic spectra then its lines shift over time in the star rest frame because of the planet orbital motion (Sect. 3.2). In that case, masking those contaminating lines when defining intrinsic spectra (Sect. 6.5) would exclude mask lines from the cross-correlation in some exposures but not others, yielding CCFs that are not equivalent over the time series. To circumvent this issue, the masking of the planetary signal is performed in velocity space, after intrinsic spectra have been converted into CCFs (Fig. 17).

The spatial accuracy at which stellar line profiles can be sampled depends on the blurring induced by the planet size and orbital motion. If we can assume that the stellar emission does not vary substantially over the planet-occulted region during an exposure, then the above equation simplifies as $$F^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}}({\lambda }_{\star },E,t)\sim C_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}({\overline{\lambda }}_{\star },E)\frac{f_{\mathrm{o}\mathrm{c}\mathrm{c}}({\lambda }_{\star },E,t)S_{\mathrm{o}\mathrm{c}\mathrm{c}}^{\mathrm{p}[\text{\hspace{0.17em}thick]\hspace{0.17em}}}({\overline{\lambda }}_{\star },E,t)}{1-lc({\lambda }_{\star },E,t)}.$$(D.8)

The surface density flux from the region occulted during the exposure at time t writes as: $$f_{\mathrm{o}\mathrm{c}\mathrm{c}}({\lambda }_{\star },E,t)=I_{\mathrm{o}\mathrm{c}\mathrm{c}}({\lambda }_{\star },E,t){\alpha }_{\mathrm{o}\mathrm{c}\mathrm{c}}({\overline{\lambda }}_{\star },E,t),$$(D.9)

where I is the specific intensity emitted by the occulted region and α its broadband intensity variations (such as limb- and gravity-darkening). During transit, the spectral flux scaling lc over a given band at ${\overline{\lambda }}_{\star }$ corresponds to (Sect. 6.3): $$1-lc({\overline{\lambda }}_{\star },E,t)=\frac{\overline{I}(v){\alpha }_{\mathrm{o}\mathrm{c}\mathrm{c}}({\overline{\lambda }}_{\star },E,t){\overline{S}}_{\mathrm{o}\mathrm{c}\mathrm{c}}(E,t)}{{\displaystyle {\sum }_{{\lambda }_{\star }}}{F}_{\star }({\lambda }_{\star },E)d{\lambda }_{\star }},$$(D.10)

where ${\overline{S}}_{\text{occ}}$ and $\overline{I}$ are the achromatic planetary surface and stellar intensity values that allow reproducing the integrated flux ratio ∑_λ* F(λ_*, E, t)dλ_*/ ∑_λ* F_*(λ_*, E)dλ_* measured in this band, as noted in Sect. D.2. While ${\overline{S}}_{\text{occ}}$ does not necessarily match the true $S_{\mathrm{o}\mathrm{c}\mathrm{c}}^{\mathrm{p}[\text{\hspace{0.17em}thick]\hspace{0.17em}}}$, because of the spectral variations of I_occ(λ_*, E, t) over the band, we assume that the ratio ${\overline{S}}_{\text{occ}}/S_{\text{occ}}^{\mathrm{p}[\text{thick}]}$ remains constant during transit, so that: $$\begin{array}{rl}{F}^{\text{intr}}({\lambda }_{\star },E,t)& =C_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}({\overline{\lambda }}_{\star },E){\displaystyle \sum _{{\lambda }_{\star }}}{F}_{\star }({\lambda }_{\star },E)d{\lambda }_{\star }\times \\ & \frac{I_{\mathrm{o}\mathrm{c}\mathrm{c}}({\lambda }_{\star },E,t){\alpha }_{\mathrm{o}\mathrm{c}\mathrm{c}}({\overline{\lambda }}_{\star },E,t)S_{\mathrm{o}\mathrm{c}\mathrm{c}}^{\mathrm{p}[\text{\hspace{0.17em}thick]\hspace{0.17em}}}({\overline{\lambda }}_{\star },E,t)}{\overline{I}(v){\alpha }_{\mathrm{o}\mathrm{c}\mathrm{c}}({\overline{\lambda }}_{\star },E,t){\overline{S}}_{\mathrm{o}\mathrm{c}\mathrm{c}}(E,t)},\end{array}$$(D.11)

Can be rewritten as: $$F^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}}({\lambda }_{\star },E,t)=F_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}({\overline{\lambda }}_{\star },E)I_{\mathrm{o}\mathrm{c}\mathrm{c}}({\lambda }_{\star },E,t).$$(D.12)

Here, all the components that are constant over time throughout the transit and with low-frequency spectral variations have been combined into $F_{\text{ref}}^{\text{norm}}$.

Planet-occulted stellar spectra are distorted by broadband spectral variations in the planet size. Indeed, the planet occults a larger area of the stellar photosphere – with a broader rv distribution and a different average rv – at wavelengths where it appears larger, for example due to Rayleigh scattering. We use ANTARESS stellar surface model to calculate the surface rv associated with each planet-occulted region as a function of wavelength (Sect. 3.3), using the same chromatic planetary transit depths and stellar intensity variations used for broadband flux scaling (Sect. 6.3). Intrinsic spectra are then corrected for chromatic deviations around the “white” rv of planet-occulted regions, so that differences in the rv of intrinsic stellar lines across the spectrum only trace variations due to the stellar atmosphere. WASP-76b yields maximum chromatic deviations on the order of 5 m s⁻¹ at ingress/egress, which is much smaller than the RV precision that can be reached on stellar lines in a local spectral band with ESPRESSO. This effect is thus likely negligible with present instrumentation, but may need to be taken into account with high-resolution spectrographs on the ELT.

Appendix E Generic methods

E.1 Stellar continuum

We detail here our method, adapted from Cretignier et al. (2020b), to determine the spectral continuum of the target star using the processed data. A master spectrum of the unocculted star is first built from each epoch or a combination of epochs, resampling spectrally (Sect. 2.3) and then temporally (Sect. 7.3) a selection of out-of-transit exposures. The master spectrum is smoothed with a Savitzky-Golay filter to remove noise with higher frequency than a chosen scale. Local maxima are identified as bins whose smoothed flux value is the highest in their closest neighbourhood, defined as a window of chosen width. The resulting spectrum of local maxima is stretched so that variations in the flux direction become smaller than those in the wavelength direction. An alpha-shape algorithm is then applied over the stretched spectrum to remove maxima formed by blended lines, setting the radius of the rolling pin to a larger size than the typical spectral line width. Finally, outlying maxima are removed and the remainder used to define the continuum of the master spectrum via linear interpolation. Smoothing and local maxima windows, stretching parameter, and rolling pin radius can be adjusted depending on the instrumental resolution and stellar line properties to better control the continuum determination.

E.2 Spectral weight profiles

E.2.1 True error estimate

The flux in a given pixel is estimated as F^meas(λ), which follows a Poisson distribution that can be approximated by a normal distribution with mean the true flux F^true(λ) and standard deviation σ^true(λ), if enough photoelectrons are measured. Because the uncertainty σ^meas(λ) associated with a measured pixel flux is proportional to $\sqrt{F^{\text{meas}}(\lambda )}$, it is a biased estimate of the true error and should thus not be used to weigh the flux. Indeed, pixels where F^meas(λ) is lower (resp. larger) than F^true(λ) due to statistical variations would be associated with a lower (resp. larger) uncertainty than their true error, resulting in artifically larger (resp. lower) weights and a weighted mean that underestimates (resp. overestimates) the true flux.

We thus weigh the flux in a given pixel using an estimate of its true error. Our underlying assumption is that the measured flux averaged over a large spectral band, or over a large number of exposures, is close enough to its true value that its measured error can be considered as an estimate of its true error. The definition of the weights then depends on the type of profile, disk-integrated or intrinsic, it is applied to. Since time series of disk-integrated and intrinsic spectra are aligned in the star or surface rest frame before being averaged, all medium- and high-frequency components involved in the weight profiles are shifted in the same way as their associated spectrum. For example telluric spectra originally defined in the barycentric rest frame are shifted to the same frame as the averaged exposures. Due to these shifts, the common flux calibration profile defined in the detector rest frame becomes specific to each exposure.

E.2.2 Definition of weight profiles

At their latest processing stage disk-integrated spectra write as (see Eq. 9): $$F^{\mathrm{s}\mathrm{c}}({\lambda }_{\star },E,t)=F^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}({\lambda }_{\star },E,t)C_{\text{corr\hspace{0.17em}}}({\lambda }_{\star },E,t)lc({\lambda }_{\star },E,t),$$(E.1)

where C_corr groups all spectral corrections that were applied to the measured spectra (Sect. 5), except for cosmics and wiggles. For the sake of clarity we only report the dependency on λ_* in the following expressions. Uncertainties on the measured flux density can be expressed as a function of the raw photoelectron count (Eq. 3), yielding the following uncertainty on the true flux F^sc,true measured during Δt: $${\sigma }^{\text{sc,true}}\left({\lambda }_{\star }\right)={\tilde{g}}_{\text{cal\hspace{0.17em}}}\left({\lambda }_{\star }\right)\sqrt{N^{\text{true\hspace{0.17em}}}\left({\lambda }_{\star }\right)}{C}_{\text{corr\hspace{0.17em}}}\left({\lambda }_{\star }\right)lc\left({\lambda }_{\star }\right)/\mathrm{\Delta }t.$$(E.2)

The spectra F^sc/lc have the same low-resolution profile as the spectrum of the unocculted star F_*, allowing us to estimate the raw number of photoelectrons as: $$N^{\text{true}}\left({\lambda }_{\star }\right)\sim \frac{F_{\star }^{\text{meas}}\left({\lambda }_{\star }\right)\mathrm{\Delta }t}{{\tilde{g}}_{\text{cal}}\left({\lambda }_{\star }\right)C_{\text{corr}}\left({\lambda }_{\star }\right)},$$(E.3)

where the stellar spectrum $F_{\star }^{\text{meas}}$, averaged over many out-of-transit exposures, is assumed to be a good estimate of $F_{\star }^{\text{true}}$. Replacing N^true with the above expression in Eq. E.2.2 then yields the weight profile on a disk-integrated spectrum defined in Eq. 15.

Then, based on Eq. 11 and 12, the error on the true intrinsic flux is: $$\begin{array}{rl}{\sigma }^{\text{intr,true}}(\lambda )& =\frac{{\sigma }^{\text{diff,true}}(\lambda )}{1-lc(\lambda )}\\ \text{\hspace{0.17em}where\hspace{0.17em}}{\sigma }^{\text{diff,true}}(\lambda )& =\sqrt{{\sigma }_{\star }^{\text{true}}(\lambda {)}^2+{\sigma }^{\text{sc,true}}(\lambda {)}^2}.\end{array}$$(E.4)

Since the measured unocculted stellar spectrum is assumed to approximate the true spectrum (Sect. E.2.2), we only use its variance with ${\sigma }_{\star }^{\text{true}}\sim {\sigma }_{\star }^{\text{meas}}$, yielding the weight profile on an intrinsic spectrum defined in Eq. 15.

E.3 Temporal and spatial resampling

Temporal/spatial resampling is performed over a selection of exposures, typically within the out-of-transit window for disk-integrated profiles and within the transit window for intrinsic profiles. All profiles are resampled on a common spectral grid prior to being binned. The start and end positions of input exposures along the binning dimension, x_start(t) and x_end(t), are defined using ANTARESS numerical grid (Sect. 3.3). The start and end positions of new exposures, $x_{\text{start}}^{\mathrm{in}}$ and $x_{\text{end}}^{\mathrm{in}}$, are defined for each exposure or by specifying a x range and number of exposures. The effective start and end positions of new exposures are then calculated as $min(max{}_t(x_{\text{start\hspace{0.17em}}}(t),x_{\text{start\hspace{0.17em}}}^{\mathrm{i}\mathrm{n}}))$ and $max(min{}_t(x_{\text{end\hspace{0.17em}}}(t),x_{\text{end\hspace{0.17em}}}^{\mathrm{i}\mathrm{n}}))$, considering all input exposures that overlap with the new exposure. The central positions of new exposures are calculated as the mean of their start/end positions. Input exposures are weighted along the binning dimension by Δx, the overlap between this exposure’s window and the window of the new exposure. We prefer to set x = r_sky rather than μ (corresponding to $\sqrt{1-r_{\text{sky}}^2}$ for a spherical star) as a binning coordinate for intrinsic profiles because it samples more regularly the stellar surface (Sect. 8.2.3). For a transiting planet on a circular orbit, with i_p ~90° and orbital phases ϕ close to 0, r_sky approximates as 2 π ϕ a_p/R_*. For the same exposure time, the r_sky ranges covered by the planet at two different positions are thus more similar than in μ.

The temporal/spatial resampling method does not interpolate input profiles along x. Instead, Eq. 16 assumes a uniform flux density over the window of each exposure. This allows resampling the flux in a given spectral bin even if it is undefined in several input exposures. Undefined spectral bins in input profiles, for example masked because of contamination by the planetary atmosphere (Sect. 3.2), have their flux and weights set to zero so that they do not contribute to the mean. The covariance matrix of the mean profile is calculated as described in Sect. 2.3, based on Eq. 16. We do not account for covariance over the temporal and spatial dimension, as successive exposures are assumed to be independent.

E.4 Conversion from 2D to 1D spectra

Spectra are converted over a 1D spectral grid common to all exposures, specific to a given instrument, and with uniform resolution in ln λ. Spectra from overlapping orders have to be equivalent in flux to prevent the introduction of biases, highlighting the importance of a proper flatfielding, deblazing, and flux balance correction. We first resample spectra in each spectral order over the common 1D grid (Sect. 2.3). Spectral weight profiles, mainly shaped here by the flux calibration profiles (Sect. 7.2), are then calculated in each order and normalized by the total weight in each pixel common to overlapping orders. Flux spectra and covariance matrixes from all spectral orders are finally multiplied by their corresponding weight profiles, and coadded.

E.5 Cross-correlation

We detail here ANTARESS cross-correlation method. It can be applied over all spectral orders or a selection, for example to generate a chromatic set of CCFs (Cristo et al. 2022). CCFs are calculated over a rv grid of constant step and chosen range.

In a first step we identify lines of the binary mask that can contribute to the cross-correlation in all processed exposures. Only lines for which the requested rv range (once converted into a wavelength range around this line) is fully within the spectral range of the processed spectrum (which is not the case, for example, if the line is at the edge of a spectral order) and contains no undefined pixels are kept. The spectrum of each exposure is then cross-correlated with the filtered mask as follows: $$CCF(rv,E,t)={\displaystyle \sum _l}F({\lambda }_l,E,t)\frac{w_l}{{⟨g_l⟩}_{\mathrm{c}\mathrm{a}\mathrm{l}}}.$$(E.5)

For each mask line l the flux spectrum and covariance matrix are resampled on the CCF grid, converted from rv to wavelengths ${\lambda }_l={\lambda }_l^0\sqrt{(1+\mathrm{rv}/\mathrm{c})/(1-\mathrm{rv}/\mathrm{c})}$ relative to the line transition ${\lambda }_l^0$. They are then scaled back approximately to raw count units using the mean of the flux calibration profile over the line grid <g_l>_cal (Sect. 5.1), so as to naturally give more weight to lines measured with a larger number of raw photoelectons without modifying the local flux balance of the spectrum. The local flux spectrum is then multiplied by the line weight w_l and added to the CCF flux grid, while the covariance matrix is multiplied by the squared line weight and stored. Once all lines have been processed the module co-adds their weighted covariance matrixes, using the largest one to set the minimum number of diagonals in the CCF matrix (Sect. 2.3). When calculating CCFs of stellar spectra, the associated stellar spectral weight (Sect. 7.2) and broadband flux scaling grid (Sect. 6.3) are cross-correlated in the same way.

Appendix F Analysis methods

F.1 Stellar mask generator

Prior to using this method the ANTARESS workflow must have been set to align spectra in a common rest frame, using either the systemic¹¹ plus Keplerian RV for disk-integrated spectra (Sect. 6.2), or the surface rv along the transit chord for intrinsic spectra (Sect. 6.5.2). Aligned spectra must then have been converted from 2D to 1D (Sect. 7.4) and resampled into a master spectrum (Sect. 7.3). No in-transit exposures should be used when computing a disk-integrated master, to avoid biases due to planet-occulted stellar lines.

The method starts by normalizing the master spectrum with the estimated stellar continuum (Sect. 7.1), smoothing it to remove high-frequency noise, and oversampling it on a regular spectral grid. Local maxima and minima are identified as bins whose flux value is the highest (resp. lowest) in their closest neighbourhood, defined as a fraction of the lines FWHM. Lines can then be defined as the unique successions of distinct maximum-minimum-maximum, with the line central wavelength set to the minimum position. Lines that fall within spectral ranges known to be contaminated by deep telluric lines, or ranges selected by the user, are excluded. Successive selections can then be applied to refine the linelist, with thresholds that can be adjusted depending on the quality of the data and the type of the star (the process is illustrated in Fig. 18):

1. Lines must have a depth (relative to the global continuum, and to their local continuum) within a given range, and a minimum depth and width (relative to the closest maximum) larger than a threshold.

2. An independent estimate of the line positions, from a polynomial fit to their core, must not deviate from their minimum position by more than a threshold.

3. Lines must be deeper than contaminating telluric lines by a threshold. A telluric line is a contaminant if, as it shifts between its minimum and maximum spectral positions, it crosses the range covered by a stellar line. The use of transit time series grants a finer control over this step than with long-term rv series. We retrieve the telluric spectra (Sect. 5.2) associated with the exposures used to build the master stellar spectrum, align them in the Earth rest frame, and bin them into a master telluric spectrum. Telluric lines are then identified following the same approach as for the stellar spectrum, and checked for overlap with stellar lines by considering the stellar Keplerian and Earth barycentric motion over the processed exposures.

4. Lines must not be too asymmetric, with criteria based on the width and depth differences between their maxima lower than a threshold.

5. Lines must have a low rv dispersion. A line-by-line analysis (Dumusque 2018; Bouchy et al. 2001) is applied to each exposure used in the master spectrum, cross-correlating their spectra to derive rv time series. Each line must have a rv dispersion over the series lower than the mean error by a threshold (this criteria should not be used when the number of exposures is too small), and an average rv calculated as the weighted mean over the series lower than a threshold.

6. A final selection is applied to exclude more finely stellar lines contaminated by tellurics.

Weights representative of the photonic error on the mask line positions (Bourrier et al. 2021) are defined as the squared inverse of the error on the line rvs (Bouchy et al. 2001).

Once a general list of stellar lines is defined, a cross-matching can be applied with a VALD linelist generated with the host star properties and for the instrument spectral range. We adjust the depth of the VALD lines so that they better match the stellar lines, then identify the unblended VALD lines most consistent with each mask line based on their relative difference in depth and position (Fig. 18). Once species have been attributed to the mask lines, a refined selection can be made.

F.2 Analysis of individual stellar profiles

F.2.1 Model line profiles

All stellar profiles are defined with a polynomial continuum, multiplied by spectral absorption lines calculated analytically or numerically. Analytical profiles provide a direct way to fit a single stellar line with simple shape in disk-integrated or intrinsic spectra. Numerical profiles are calculated using ANTARESS stellar surface model (Sect. 3.3), summing elementary profiles over the full star to fit disk-integrated profiles or over planet-occulted regions to fit intrinsic profiles. As described in Sect. 3.3.1 the elementary profiles can be analytical, theoretical, or measured. The last two options allow fitting several lines over a portion of the stellar spectrum (Fig. 23).

Analytical (either used directly or in the numerical model) and theoretical profiles are defined over regular pixel grids so that they can be convolved with the instrumental response function (assumed to be a Gaussian with FWHM preset for each instrument). When working with low-resolution data, profiles are calculated over an oversampled grid before convolution, and then resampled on the pixel grid of the fitted profile, to avoid blurring. Measured profiles are already broadened by instrumental convolution and kept on their grid of origin (Sect. 6.5).

F.2.2 Line profile fitting

Disk-integrated profiles from original time series are analyzed in the rest frame of the input data, where their derived centroids can be used to check the Keplerian model and measure the systemic rv of the system. After being aligned, profiles can be analyzed in the star rest frame. We highlight that masking narrow planetary atmospheric features does not remove the effect of planetary occultation on disk-integrated profiles. Indeed, the opaque planet continuum generates spectral features through its absorption of the local stellar line profile (e.g., Dethier & Bourrier 2023). These features yield the so-called “Doppler shadow” (Collier Cameron et al. 2010) that distorts in-transit line profiles and induces the RM anomaly when deriving their centroid with an undistorted model of the unocculted star. We propose an original method to circumvent this bias, which consists in excluding the range absorbed by the planet-occulted stellar line. This range can be defined in a given exposure as a rv window covering the breadth of the line, centered around the line transition Doppler-shifted by the surface rv of the planet-occulted region (Eq. 2). Masked disk-integrated stellar line profiles are not distorted anymore by the planetary continuum, and can be fitted with the unocculted stellar line model.

Intrinsic profiles from original time series are analyzed in the star rest frame, where their derived centroids can be used to determine the best model for stellar surface rv along the transit chord. After being aligned, profiles can be analyzed in the photosphere rest frame. We highlight the importance of using a Bayesian approach when fitting intrinsic profiles (Bourrier et al. 2021), to better explore the line model properties and their correlations, and to retrieve information from noisy or masked local stellar lines that could not be obtained through least-square minimization. Priors set on the model properties can be informed by their range of variations along the transit chord (Sect. 8.2.3) or physical arguments (e.g., local FWHM smaller than the disk-integrated value; contrast between 0 and 1; RV centroid smaller than the projected stellar rotational velocity). PDFs from the fits to intrinsic stellar lines inform on the exposures that can be exploited in time series analysis (Sect. 8.3). Exposures with flat, non-constraining PDFs for the width, contrast, and rv of the intrinsic line, typically at the limb of the star where the planet only partially occults dim regions of the stellar surface, should be excluded (Fig. 17).

F.2.3 Line profile variations

Parametric models are used to describe the variation of a line property as a function of the chosen variable. They are fitted to the properties measured over all epochs or a selection of epochs, with parameters that can be common or independent between fitted datasets to allow for variations between epochs. Disk-integrated properties measured out-of-transit can be studied as a function of ambient (e.g., airmass, seeing measured at the instrument facility), data-related (e.g., time, mean S/N over a given selection of orders), or stellar (e.g., activity indexes, retrieved from the DACE platform) variables. They are fitted directly with the chosen parametric model, which can be a polynomial, a sine function, or their combination (e.g., additive for RVs, multiplicative for FWHM and contrast). Intrinsic properties, by definition measured in-transit, are studied as a function of stellar coordinates (e.g., center-to-limb angle, sky-projected distance from star center, latitude or longitude). They are fitted with brightness-averages over the planet-occulted regions (Sect. 3.3.2), tiled with the chosen parametric model. Intrinsic line centroids, which trace the surface RVs along the transit chord, are described with Eq. 2. Morphological line properties m are described with absolute (m(x) = ∑_i≥0 c_ixⁱ)) or modulated (m(x) = m₀(1 + ∑_i≥1 c_ixⁱ)) polynomial models. The latter possibility allows for a common dependence of the property with stellar coordinate x, with a scaling m₀ specific to each epoch.

F.3 Analysis of joint intrinsic stellar profiles

F.3.1 Global model line profiles

Series of intrinsic profiles are fitted over multiple epochs with a global model for the stellar lines, defined analytically or numerically (Sect. 8.2). Line profiles are linked between consecutive exposures through models describing their shape and position as a function of stellar coordinates, informed by the analysis of individual exposures (Sect. 8.2.3). Profiles can further be linked between epochs by using common model coefficients, although they can be varied to account for stellar variability. The line profiles are defined before convolution with each instrument’s response, so that their coefficients remain comparable between different spectrographs.

F.3.2 RMR fitting

The surface rv model controlling the theoretical line centroids (Eq. 2) depends at first order on the sky-projected spin-orbit angle λ and stellar rotational velocity v_eq sin i_*. If allowed by the data, the rv model can further depend on convective blueshift and on differential rotation, the latter breaking the degeneracy on the stellar latitudes transited by the planet. The resulting measurement of the stellar inclination i_*, combined with λ and v_eq sin i_*, respectively yields the 3D spin-orbit angle ψ and the equatorial stellar rotational velocity v_eq. Alternatively, ψ can be derived using an independent constraint on the stellar inclination (for example through asteroseismology) or the equatorial rotation period P_eq (typically using photometric modulation or activity indexes).

The coefficients describing the line profiles depend on the type of profile (e.g., FWHM and contrast for analytical profiles, abundances for theoretical profiles, etc), but their spatial dependence along the transit chord is always described using the polynomials described in Sect. F.2.3. We remind that analytical models for intrinsic lines can include macroturbulence (Sect. 3.3.1), which can be used to separate its contribution from that of thermal broadening.

By construction, measured intrinsic profiles are corrected for the broadband flux scaling associated with stellar intensity variations and the area occulted by the planet. However, an intrinsic profile can be seen as the brightness-average of all elementary stellar profiles occulted by the planet during an exposure (Sect. 3.3), so that its shape and Doppler-shift still depend on the size of the planet and on the stellar intensity properties. Model intrinsic profiles thus also depend on the planet-to-star radius ratio in the fitted band (which determines which stellar grid cells are occulted) and on the stellar gravity/limb-darkening coefficients (which determine the relative weights of the elementary profiles from each cell). In practice, intrinsic line profiles are however not measured to a sufficiently good precision, even in CCF format, to constrain these properties. On the other hand the scaled semi-major axis a_p/R_* and orbital inclination i_p can in some cases be constrained by the data, and set as free parameters in the joint fit. Indeed, those parameters can be degenerate with either λ (for aligned orbits) or v sin i_* (for polar orbits), and affect the model stellar lines by controlling the theoretical transit chord and thus the occulted stellar grid cells. Importantly, the extraction of intrinsic profiles does not require a precise knowledge of the planet-to-star radius ratio, semi-major axis, and orbital inclination, but only of the transit light curve used for broadband scaling (Sect. 6.3, 6.5).

Appendix G Performance analysis

G.1 Resampling

Spectra are measured in different exposures over slighly different wavelength grids, and the subsequent shifts required to align them in common rest frames (Sect. 6.2, 6.5.2) further differentiate their spectral grids. Resampling spectra on a common grid, as is commonly done in transmission spectroscopy analyses, blurs sharp spectral features and introduces correlations between pixels (Sect. 2.3). The ANTARESS workflow innovates by avoiding resampling as much as possible, and we evaluate here the benefits of this approach using the first epoch of HD 209458 as an example.

In our nominal approach, a model telluric spectrum is calculated on the native grid of each observed spectrum before calculating telluric CCFs of the spectra and their models (Sect. 5.2). If spectra are first resampled on a common grid, the best-fit telluric CCFs are shallower due to the blurring of the observed telluric lines, which are then undercorrected (Fig. G.1).

We fitted Gaussian profiles to deep and narrow stellar lines in the red, high S/N parts of the spectrum to assess the impact of resampling on the line shape. Although the lines were systematically shallower and broader when resampled, the difference remained marginal (1σ for the constrast, 0.5σ for the FWHM, 0.05σ for the rv). We derived the same conclusion from fits to the cross-correlated stellar line averaged over out-of-transit exposures. We thus conclude that resampling spectra after they are aligned in the star rest frame has a negligible impact on individual stellar lines. This is because the Keplerian motion of HD 209458 around the transit windows only represents a small fraction of an ESPRESSO pixel, so that the interpolation due to resampling remains limited.

The impact of resampling intrinsic spectra is mitigated by their lower S/N. Nonetheless, the cross-correlated stellar line averaged along the transit chord shows differences on the order of 0.1%, 10 m s⁻¹ and 1 m s⁻¹ for the contrast, FWHM, and rv position, with typical significance of ~1σ. This is because intrinsic spectra aligned in the rest frame of the stellar surface are corrected for its Doppler motion, which represents here a substantial fraction of a pixel. The RMR fit is similarly impacted, with differences at the ~1σ level for λ and v_eq sin i_* (1.055±0.077°, 4.272±0.007 km s⁻¹ for the nominal analysis; 1.183±0.073°, 4.262±0.007 km s⁻¹ for resampled data).

Finally, one of the interests of the ANTARESS workflow is the possibility to oversample CCFs with a rv step smaller than the instrumental pixel width, thus better resolving the line profiles. The resulting correlations between pixels are propagated thanks to the use of covariance matrixes, and accounted for when fitting the CCFs (Sect. 2.3). To illustrate this point we calculated the CCFs of HD 209458 intrinsic spectra with rv steps of 0.25 and 0.10 km s⁻¹, smaller than the nominal value of 0.50 km s⁻¹. Fits to the master CCF along the transit chord in epoch 1 show a significant refinement of the stellar line shape from 0.5 to 0.25 km s⁻¹, as it gets deeper (contrast from 0.5683±0.0009 to 0.58693±0.0010) and narrower (FWHM from 7.656±0.016 to 7.199±0.018 km s⁻¹). Using even finer rv steps does not resolve better the line shape (Fig. G.2). Oversampling the CCFs does not change the results of the RMR fit, likely because the line profiles of HD 209458 are symmetrical and their positions are already measured with high precision at the nominal rv step.

G.2 Covariance matrix

The ANTARESS workflow calculates the covariance between pixels of the processed profiles, and propagates it as a banded covariance matrix (Sect. 2.3). With the optimal processing, when 2D spectra are resampled as late as possible, no covariance exists until differential profiles are extracted (Sect. 6.4). At this stage, correlations are introduced through the difference between in-transit profiles and the master out-of-transit spectrum used as reference, since disk-integrated spectra must be aligned and resampled on a common grid before they can be averaged into a master. The master spectra of HD 209458 and WASP-76 are associated with three covariance diagonals, whose ratios to the variance (from ~3× 10⁻³ to 4× 10⁻⁵) however remain negligible. Strong correlations are then inevitably created when computing CCFs, although this is mitigated in the ANTARESS workflow by accounting for the covariance in fits. The CCFs of intrinsic spectra are associated with eight covariance diagonals, with typical ratios to the variance decreasing from 0.2 to 2× 10⁻¹⁴. However, covariances remain small enough that discarding them in the RMR fits does not change the derived values and their uncertainties for either HD 209458b or WASP-76b.

G.3 Weighted means

Particular care is taken in the ANTARESS workflow to calculate the weighted mean of spectra (Sect 7.2), using estimates of the true error on individual pixels. Here we evaluate the impact of not considering various components in the weighing of disk-integrated and intrinsic spectral profiles. Our goal is not to perform an exhaustive assessment but rather to estimate the typical biases introduced by improper weighing.

Figure G.3 illustrates the role of the instrumental flux calibration profile when averaging overlapping spectral orders in a given disk-integrated spectrum (Sect. E.2.2), typically to convert 2D echelle spectra into a 1D spectrum. We observe differences up to a few percents, increasing toward the edges of the orders, when performing a simple mean that does not account for lower photon count and associated large flux errors.

Figure G.4 illustrates the role of telluric correction when averaging disk-integrated spectra from different exposures (Sect. E.2.2). These spectra are typically averaged in the star rest frame, in which telluric lines shift over time, so that a given pixel should be weighted differently according to its level of telluric correction. The strongest biases are observed for the deepest telluric lines when they are not included in the weight profile, but even moderately strong tellurics can induce biases up to a few tens of ppm when averaging two exposures. The bias is reduced when more exposures are averaged, as a pixel absorbed by tellurics in a given exposure will weigh less in the overall mean. Nonetheless, we observe differences up to a few hundreds of ppm in regions of deep telluric lines, even when averaging all out-of-transit exposures.

Figure G.5 illustrates the role of disk-integrated stellar lines when averaging intrinsic spectra in the photosphere rest frame (Sect. E.2.2). In that case, disk-integrated stellar lines shift between different exposures, so that a given pixel should be weighted differently according to the depth of the stellar line it overlaps with. Similarly to telluric lines, the strongest biases are observed for the deepest disk-integrated stellar lines, when they are not included in the weight profile. The sodium lines are among the deepest in the disk-integrated spectrum of HD 209458, with nearly full absorption at the core of the line, and induce biases up to tens of percents even in the intrinsic profile averaged over all in-transit exposures. It is thus critical for the analysis of averaged planet-occulted stellar lines to account for disk-integrated stellar lines in weight profiles.

G.4 Flux balance corrections

Correcting for spectro-temporal variations of the flux balance is necessary to avoid biases in planetary transmission spectra. While local corrections can be applied at high spectral resolution, forcing the continuum around the line of interest to a constant value, the use of a spectrum-wide correction in the ANTARESS workflow may allow the retrieval of planetary signatures at medium spectral resolution. This type of correction also ensures the computation of unbiased stellar and planetary CCFs, for which it is not possible to correct the continuum of each independent line. Relative variations between the continua of different lines over time indeed translate into variations between their relative weights in the CCF, modifying its profile between exposures.

Here we evaluate how much a poor correction of the flux balance impacts the quality of disk-integrated and intrinsic stellar CCFs. We reprocessed the data from the FLUX BALANCE method (Sect. 5.3) onward, using low-order polynomials instead of smoothing splines to fit the color balance variations. This approach allows capturing the overall shape of the color balance, but not its variations over a few spectral orders as finely as with splines (Fig. 7). We first assessed the quality of disk-integrated CCFs by comparing the dispersion of their out-of-transit properties with the nominal analysis. The degraded flux balance correction does not have a clear impact on WASP-76 CCFs. For HD 209458, it increases the dispersion of the CCF contrast and FWHM by ~3–8% but has no impact on their RVs. This is likely because the residual, uncorrected color variations are broad enough that they are approximately constant over the breadth of a stellar line and thus do not modify their shape. Since the correction improves the disk-integrated CCFs of HD 209458 we further evaluated its impact on the intrinsic CCFs and their RMR fit. The degraded correction yields λ = 1.12±0.07° and v_eq sin i_* = 4.289±0.008 km s⁻¹, consistent with 1.06±0.07°, 4.272±0.007 km s⁻¹ from the nominal analysis.

Our analysis suggests that an accurate flux balance correction, while important for spectral analysis and to compute high-quality CCFs, is not critical for the derivation of RVs and the interpretation of the RM effect. Biases induced by the color effect (e.g., Bourrier & Hébrard 2014) mainly arise from the spectrum-wide variations reaching tens of percents, which can be captured by low-order polynomials.

G.5 Wiggle correction

The same spectral region in ESPRESSO data is affected by different wiggle patterns over time (Fig. G.6). Applying a global, homogeneous correction to all spectra acquired during an epoch (Sect. 5.6) is thus required to compare a given spectral line between different exposures, or to combine a line assumed to be stable over several exposures.

Here, we assess its impact on stellar CCFs and their interpretation. Similarly to color balance variations, but at shorter frequencies, wiggles indeed change the relative flux balance between different spectral regions over time and thus the relative weights between CCF lines. We deactivated the wiggle correction (Sect. 5.6) and reprocessed the data from the SPECTRAL LINE DETRENDING method (Sect. 6.1) onward. Correcting for the wiggles has no significant impact on WASP-76 disk-integrated CCFs. For HD 209458, it improves the dispersion of their out-of-transit contrast and FWHM by ~3–6% but has no impact on their rvs. This is likely because the wiggle pattern is too broad to affect the line shapes and thus only changes their overall flux level, as discussed for the flux balance correction (Sect. G.4). We then evaluated the impact of the wiggles on the intrinsic CCFs of HD 209458 and their RMR fit. Disabling the correction yields λ = 1.01±0.08° and v_eq sin i_* = 4.271±0.007 km s⁻¹, to be compared to 1.06±0.07°, 4.272±0.007 km s⁻¹ for the nominal analysis.

Although they have a strong impact on spectral data, wiggles thus appears to have a moderate to negligible impact on disk-integrated and intrinsic CCFs. This is likely because wiggles are smoothed out over the thousands of cross-correlated lines, so that temporal variations in the wiggle pattern have little impact on the CCF profiles. This suggests that wiggle correction is not critical to RMR fits on ESPRESSO data, although it may be safer to apply it to ensure unbiased results, especially if a small number of stellar lines are used in the CCFs.

G.6 CCF from intrinsic spectra / from disk-integrated CCF

The ANTARESS workflow offers the new possibility to apply the RMR technique (Sect. 8.3.1) to “white” CCFs calculated by cross-correlating intrinsic spectra over broad spectral ranges: $$CCF\left[F^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}}\right](rv)={\displaystyle \sum _l}\frac{F_{\star }^t\left({\lambda }_l\right)-F^{\mathrm{s}\mathrm{c}}\left({\lambda }_l\right)}{1-\mathrm{lc}\left({\lambda }_l\right)}{\alpha }_l,$$(G.1)

where α_l designates the combined weights applied to local line spectra (Sect. 7.5). It is also possible, as in past studies (e.g., Bourrier et al. 2021), to use the CCFs of disk-integrated spectra, readily provided by the DRS of high-resolution spectrographs or calculated with the ANTARESS workflow (Sect. 6.5), and extract intrinsic CCFs as: $$CC{F}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}}(rv)=\frac{CC{F}_{\star }(rv)-CC{F}^{\mathrm{s}\mathrm{c}}(rv)}{1-l{c}_{\mathrm{w}\mathrm{h}}},$$(G.2)

where lc_wh (Sect. 6.3) is either the achromatic light curve calculated over the “white” band best matching the CCF linelist (when CCFs are provided as input of ANTARESS), or the cross-correlation of the chromatic light curve used to scale disk-integrated spectra (when CCFs are calculated within our workflow). The CCF[F^intr] and CCF^intr are thus not equivalent when a chromatic light curve is used, as with WASP-76b, and more generally because the CCF_* and CCF^sc are calculated by resampling, weighing, and scaling disk-integrated CCFs in rv space.

To evaluate the bias induced by these approximations, we compared the nominal RMR fit of the CCF[F^intr] (Sect. 8.3.1) with the same fit applied to CCF^intr derived from disk-integrated CCFs calculated with the ANTARESS workflow. The quality of the fits and the precision of the results is similar for both planets $({\chi }_r^2=1.1)$, with consistent results for HD 209458b ($\lambda =1.164\pm {0.069}^{\circ },v_{\mathrm{eq}}\sin i_{\star }={4.263}_{-0.0083}^{+0.0072}\mathrm{km}{\mathrm{s}}^{-1}$ from the CCF^intr; λ = 1.06 ± 0.07°, v_eq sin i_* = 4.272±0.007 km s⁻¹ from the nominal CCF[F^intr]). While the derived rotational velocities are also consistent for WASP-76b (${0.79}_{-0.15}^{+0.12}\mathrm{km}{\mathrm{s}}^{-1}$ from CCF^intr; ${0.86}_{-0.19}^{+0.13}\mathrm{km}{\mathrm{s}}^{-1}$ from CCF[F^intr]), the fit to the CCF^intr prefers a positive solution for λ (${27.0}_{-18.6}^{+{29.4}^{\circ }}$, compared to the nominal $‐{37.1}_{-20.6}^{+{12.6}^{\circ }}$). It is noteworthy that our fit to the CCF^intr yields similar results as those derived by Ehrenreich et al. (2020) without prior on the rotational velocity $(v_{\mathrm{eq}}\sin i_{\star }={0.73}_{-0.03}^{+0.21}\mathrm{km}{\mathrm{s}}^{-1}$, $\lambda ={10.2}_{-5.7}^{+{31.0}^{\circ }}$; priv. comm.), supporting the conclusion from Sect. 8.3.1 that their different surface rv series and results arise from their use of intrinsic CCFs derived from disk-integrated CCFs.

We remind that both types of intrinsic CCFs analyzed here are derived from the same 2D spectra processed by the SPECTRAL LINE DETRENDING method (Sect. 6.1). The differences we observe are thus entirely caused by cross-correlating spectra too early in the workflow. We thus recommend analyzing CCFs calculated from intrinsic spectra rather than extracted from disk-integrated CCFs, as the latter can bias the analysis of the RM effect.

G.7 Chromatic flux scaling

The ANTARESS workflow can account for broadband variations in the planetary absorption and stellar emission when rescaling disk-integrated spectra F^sc to their correct relative flux level, and when rescaling intrinsic spectra F^intr to a common flux level (Sect. 6.3). As can be seen in Eq. G.6, using an achromatic scaling light curve instead of the true chromatic one thus biases the continuum of the planet-occulted spectrum $F_{\star }^t-F^{\text{sc}}$ and its intrinsic normalization by 1 – lc. This changes the relative weight given to stellar lines from different regions of the spectrum when cross-correlating intrinsic spectra, and may thus bias the resulting CCF properties.

To investigate this issue we performed a RMR fit of the CCF[F^intr] from WASP-76b, processed from the BROADBAND FLUX SCALING method onward with the white transit depth and limb-darkening values from Table 2. The achromatic processing yields $v_{\text{eq}}\sin i_{\star }={0.78}_{-0.10}^{+0.08}\mathrm{km}{\mathrm{s}}^{-1}$ and $\lambda =‐{23.9}_{-24.7}^{+{16.6}^{\circ }}$, to be compared with the nominal values ${0.86}_{-0.19}^{+0.13}\mathrm{km}{\mathrm{s}}^{-1}$ and $‐{37.1}_{-20.6}^{+{12.6}^{\circ }}$. While the derived properties remain consistent, they are impacted by the chromaticity of the scaling light curve. Future studies of data measured at higher precision, or for planets with strong chromatic variations of their continuum, should thus be wary of this bias.

G.8 Stellar CCF masks

To illustrate the use of the STELLAR MASK GENERATOR (Sect. 8.1) we built custom masks from the disk-integrated spectra of HD 209458 and WASP-76 (Fig. 18).

We assessed the quality of disk-integrated CCFs derived with the custom masks by comparing the dispersion of their out-of-transit properties with those of CCFs derived with the standard F9 ESPRESSO DRS mask. For HD 209458, the custom CCFs are less stable in contrast, more stable in FWHM, and comparable in RVs to the DRS CCFs. For WASP-76, the custom CCFs are less stable in contrast and FWHM, but more stable in RVs. Given that HD 209458 and WASP 76 have types F9 and F7, these results are not unexpected. Indeed, from a similar comparison over a wider stellar sample, Bourrier et al. (2023) noted that masks customized to a specific target star, rather than representative of a spectral type proxy as in ESPRESSO DRS, yield CCFs of comparable quality for F-type stars and CCFs of much better quality for G-type and especially K-type stars. As a complementary analysis, we performed a RMR fit of the CCF[F^intr] from HD 209458, cross-correlating intrinsic spectra with its custom mask (Sect. 6.5.1). The fit is of similar quality than the nominal one based on the F9 mask, and yields consistent results.

These tests show that custom masks, generated with the ANTARESS workflow from input datasets of a given star, allow building CCFs of at least comparable quality to standard masks, with the interest of performing self-consistent RM analysis.

References

---

All Tables

All Figures