# The ANTARESS workflow: I. Optimal extraction of spatially resolved stellar spectra with high-resolution transit spectroscopy

V. Bourrier, J.-B. Delisle, C. Lovis, H. M. Cegla, M. Cretignier, R. Allart, K.

Al Moulla, S. Tavella, O. Attia, D. Mounzer, et al.

► To cite this version:

V. Bourrier, J.-B. Delisle, C. Lovis, H. M. Cegla, M. Cretignier, et al., The ANTARESS workflow:

I. Optimal extraction of spatially resolved stellar spectra with high-resolution transit spectroscopy.

Astronomy & Astrophysics - A&A, 2024, 691, 10.1051/0004-6361/20249203. insu-04836790

HAL Id: insu-04836790

https://insu.hal.science/insu-04836790v1

Submitted on 15 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The ANTARESS workflow

## I. Optimal extraction of spatially resolved stellar spectra with high-resolution transit spectroscopy

V. Bourrier<fnref n="1" />,*, J.-B. Delisle<fnref n="1" />, C. Lovis<fnref n="1" />, H. M. Cegla<fnref n="1" />,<fnref n="2" />, M. Cretignier<fnref n="1" />,<fnref n="3" />, R. Allart<fnref n="4" />,<fnref n="1" />,**, K. Al Moulla<fnref n="1" />, S. Tavella<fnref n="1" />, O. Attia<fnref n="1" />, D. Mounzer<fnref n="1" />, V. Vaulato<fnref n="1" />, M. Steiner<fnref n="1" />, T. Vrignaud<fnref n="1" />,<fnref n="5" />, S. Mercier<fnref n="1" />, X. Dumusque<fnref n="1" />, D. Ehrenreich<fnref n="1" />, J. V. Seidel<fnref n="1" />,<fnref n="6" />, A. Wyttenbach<fnref n="1" />, W. Dethier<fnref n="7" />, and F. Pepe<fnref n="1" />

Received 10 January 2024 / Accepted 13 July 2024

## ABSTRACT

High-resolution spectrographs open a detailed window onto the atmospheres of stars and planets. As the number of systems observed with different instruments grows, it is crucial to develop a standard in analyzing spectral time series of exoplanet transits and occultations, for the benefit of reproducibility. Here, we introduce the ANTARESS workflow, a set of methods aimed at processing high-resolution spectroscopy datasets in a robust way and extracting accurate exoplanetary and stellar spectra. While a fast preliminary analysis can be run on order-merged 1D spectra and cross-correlation functions (CCFs), the workflow was optimally designed for extracted 2D echelle spectra to remain close to the original detector counts, limit the spectral resampling, and propagate the correlated noise. Input data from multiple instruments and epochs were corrected for relevant environmental and instrumental effects, processed homogeneously, and analyzed independently or jointly. In this first paper, we show how planet-occulted stellar spectra extracted along the transit chord and cleaned from planetary contamination provide a direct comparison with theoretical stellar models and enable a spectral and spatial mapping of the photosphere. We illustrate this application of the workflow to archival ESPRESSO data, using the Rossiter-McLaughlin effect Revolutions (RMR) technique to confirm the spin-orbit alignment of HD 209458b and unveil biases in WASP-76b's published orbital architecture. Because the workflow is modular and its concepts are general, it can support new methods and be extended to additional spectrographs to find a range of applications beyond the proposed scope. In a companion paper, we will present how planet-occulted spectra can be processed further to extract and analyze planetary spectra decontaminated from the star, providing clean and direct measurements of atmospheric properties.

**Key words.** methods: data analysis – techniques: spectroscopic – planets and satellites: atmospheres – stars: atmospheres

## 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 100 000 (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

<footnote n="1">¹ Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51b, 1290 Versoix, Switzerland</footnote>

<footnote n="2">² Physics Department, University of Warwick, Coventry CV4 7AL, UK</footnote>

<footnote n="3">³ Centre for Exoplanets and Habitability, University of Warwick, Coventry CV4 7AL, UK</footnote>

<footnote n="4">⁴ Department of Physics, University of Oxford, Oxford OX1 3RH, UK</footnote>

<footnote n="5">⁵ Département de Physique, Institut Trottier de Recherche sur les Exoplanètes, Université de Montréal, Montréal, Québec H3T 1J4, Canada</footnote>

<footnote n="6">⁶ Institut d'astrophysique de Paris, CNRS, UMR 7095, Sorbonne Université, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="7">⁷ European Southern Observatory, Alonso de Córdova 3107, Vitacura, Región Metropolitana, Chile</footnote>

<footnote n="8">⁸ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote>* Corresponding author; vincent.bourrier@unige.ch</footnote>

<footnote>** Trottier Postdoctoral Fellow.</footnote>
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 veloc-
ity 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). Measure-
ments 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 tomogra-
phy (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, effec-
tively, 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 con-
ceptual 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 plane-
tary 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 signa-
tures 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 spec-
trum 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 emis-
sion 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 atmo-
sphere. Spectrally resolving the planetary lines provides further
constraints on the thermal and dynamical broadening mecha-
nisms in the atmosphere (Seidel et al. 2019, 2020; Wyttenbach
et al. 2020). Combined with the possibility to resolve absorp-
tion 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 tel-luric 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, transmis-sion spectra of planetary atmospheres are usually corrected for the stellar lines they filter using the out-of-transit stellar spec-trum. 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 distor-tions (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 sur-face 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 sep-arating 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. Instru-
mental 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 fea-
tures (e.g., YARARA, Cretignier et al. 2021; APERO, Cook et al.
2022). Transit observations offer the additional advantage of pro-
viding 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 developed for ground-based, stabilized high-resolution spectrographs, and some corrections are specific to a given instrument

<table><caption>**Table 1.** Spectrographs suitable to the ANTARESS workflow.</caption><thead><tr><td>Spectrograph</td><td>Reference</td><td>Status</td></tr></thead><tbody><tr><td>CARMENES</td><td>Quirrenbach et al. (2016)</td><td>Applied</td></tr><tr><td>CORALIE</td><td>Queloz et al. (2001)</td><td>Applied</td></tr><tr><td>ESPRESSO 1/4UT</td><td>Pepe et al. (2021)</td><td>Applied</td></tr><tr><td>EXPRES</td><td>Jurgenson et al. (2016)</td><td>Applied</td></tr><tr><td>GIANO</td><td>Claudi et al. (2018)</td><td>Applicable</td></tr><tr><td>HARPS</td><td>Mayor et al. (2003)</td><td>Applied</td></tr><tr><td>HARPS-N</td><td>Cosentino et al. (2012)</td><td>Applied</td></tr><tr><td>HST/STIS</td><td></td><td>Applicable</td></tr><tr><td>JWST/NIRSPEC</td><td></td><td>Applicable</td></tr><tr><td>MAROON-X</td><td>Seifahrt et al. (2018)</td><td>Applicable</td></tr><tr><td>NIRPS</td><td>Bouchy et al. (2017)</td><td>Applied</td></tr><tr><td>SOPHIE</td><td>Perruchot et al. (2008)</td><td>Applied</td></tr><tr><td>SPIRou</td><td>Donati et al. (2020)</td><td>Applicable</td></tr><tr><td>UVES</td><td>Dekker et al. (2000)</td><td>Applicable</td></tr></tbody></table>

**Notes.** The list of spectrographs that could be applicable to the workflow is not exhaustive.

(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
<figure><img src="image_2.png" /><figcaption>Fig. 1. Chart of the ANTARESS process flow.</figcaption></figure>

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 λB in the solar system barycentric rest frame. The symbol λ 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 :

$$F(\lambda, E, t) = F_{\star}(\lambda, E) \delta_{p}(\lambda, E, t), \quad (1)$$

where F$_{\star}$ 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, δ$^{th}$, the fraction of starlight reflected by the planet, δ$^{refl}$, and the ratio of the thermal planetary flux to the stellar flux, δ$^{th}$. The flux coming from the planet can thus also be written as F$_{p}$ (λ, E, t) = F$_{\star}$ (λ, E) (δ$^{refl}$ (λ, E, t) + δ$^{th}$ (λ, 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><caption>**Table 2.** Properties of the HD 209458 and WASP-76 system used in our analysis.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value<br>HD 209458</th><th>Value<br>WASP-76</th><th>Unit</th></tr></thead><tbody><tr><td colspan="5"><strong>Star</strong></td></tr><tr><td>Spectral type</td><td></td><td>F7</td><td>F9</td><td></td></tr><tr><td>Radius</td><td>R<sub>★</sub></td><td>1.160 [1]</td><td>1.77 [4]</td><td>R<sub>☉</sub></td></tr><tr><td>Microturbulence velocity</td><td>ξ</td><td>1.03 [1]</td><td>1.38 [4]</td><td>km s<sup>-1</sup></td></tr><tr><td>Effective temperature</td><td>T<sub>eff</sub></td><td>6069 [1]</td><td>6316 [4]</td><td>K</td></tr><tr><td>Surface gravity</td><td>log g</td><td>4.41 [1]</td><td>4.13 [4]</td><td>cgs</td></tr><tr><td>Metallicity</td><td>[Fe/H]</td><td>0.02 [1]</td><td>0.34 [4]</td><td></td></tr><tr><td>Projected rotational velocity</td><td>ν sin i<sub>★</sub></td><td>4.272 [0]</td><td>0.86 [0]</td><td>km s<sup>-1</sup></td></tr><tr><td>Limb-darkening coefficients</td><td>u<sub>1</sub><br>u<sub>2</sub></td><td>0.380 [1]<br>0.234 [1]</td><td>0.393 [5]/Chromatic [6]<sup>(‡)</sup><br>0.219 [5]/Chromatic [6]<sup>(‡)</sup></td><td></td></tr><tr><td>Stellar reflex motion amplitude</td><td>K<sub>★</sub></td><td>84.27 [2]</td><td>116.02 [5]</td><td>m s<sup>-1</sup></td></tr><tr><td colspan="5"><strong>Planet</strong></td></tr><tr><td>Orbital period</td><td>P</td><td>3.52474859 [2]</td><td>1.80988198 [5]</td><td>days</td></tr><tr><td>Eccentricity</td><td>e</td><td>0 [2]</td><td>0 [5]</td><td></td></tr><tr><td>Argument of periastron</td><td>ω</td><td>90 [2]</td><td>90 [5]</td><td>deg</td></tr><tr><td>Mid-transit time</td><td>T<sub>0</sub>-2 450 000</td><td>4560.80676 [1]<sup>(†)</sup></td><td>8080.626165 [5]</td><td>BJD<sub>TDB</sub></td></tr><tr><td>Orbital inclination</td><td>i</td><td>86.71 [1]<sup>(†)</sup></td><td>89.623 [5]</td><td>deg</td></tr><tr><td>Scaled semi-major axis</td><td>a<sub>p</sub>/R<sub>★</sub></td><td>8.90 [1]<sup>(†)</sup></td><td>4.08 [5]</td><td></td></tr><tr><td>Projected spin-orbit angle</td><td>λ</td><td>1.06 [0]</td><td>-37.1 [0]</td><td>deg</td></tr><tr><td>Planet-to-star radius ratio</td><td>R<sub>p</sub>/R<sub>★</sub></td><td>0.12086 [3]</td><td>0.10852 [5] / Chromatic [6]<sup>(‡)</sup></td><td></td></tr></tbody></table>

**Notes.** [0] Our work; [1] Casasayas-Barris et al. (2021); [2] Bonomo et al. (2017); [3] Torres et al. (2008); [4] Tabernero et al. (2021); [5] Ehrenreich et al. (2020); [6] Fu et al. (2021) (†) We use the values adjusted from the RM analysis in [1]. (‡) We use the values measured at low spectral resolution by [6].

<table><caption>**Table 3.** Nomenclature of the main types of profiles across the workflow.</caption><thead><tr><th>Symbol</th><th>Type</th><th>Frame</th><th>Description</th></tr></thead><tbody><tr><td>F<sup>meas</sup></td><td>Disk-integrated</td><td>Heliocentric</td><td>Profiles integrated over the stellar disk, initially provided by a DRS, then corrected by the workflow.</td></tr><tr><td>F<sup>sc</sup></td><td>Disk-integrated</td><td>Stellar</td><td>F<sup>meas</sup> profiles aligned and scaled to their correct relative flux levels.</td></tr><tr><td>F<sup>diff</sup></td><td>Differential</td><td>Stellar</td><td>Difference profiles between a reference for the unocculted star and each F<sup>sc</sup>.</td></tr><tr><td>F<sup>intr</sup></td><td>Intrinsic</td><td>Stellar, photosphere</td><td>In-transit F<sup>diff</sup> profiles, i.e., planet-occulted stellar profiles, corrected for planetary contamination.</td></tr></tbody></table>

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 $\sigma^2(\lambda, E, t)$, correspond to the variance of the signal measured in the pixel at wavelength $\lambda$, 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 through-out 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<fnref n="1" /> 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

<figure><img src="image_3.png" /></figure>

**Fig. 2.** Views of ANTARESS sky-projected stellar grid for an imaginary system in various configurations highlighting the workflow possibilities. The plots X and Y axis correspond to $x_{*\text{sky}}$ and $y_{*\text{sky}}$ (see text). The stellar equator is shown as a solid black line, and the projected stellar spin axis as a black arrow extending from the north pole. Two planets (black disks) move along their orbital trajectory, whose 3D orientation is indicated by the normal to the orbital plane (green and golden arrows). *Top*: star shown rotating differentially between equator and poles and is colored as a function of its photospheric rv field (dotted black lines show iso-rv curves). *Middle*: star is fast-rotating and oblate, colored as a function of its intensity field (darkened at the limbs and brightened at the poles by gravity-darkening). *Bottom*: zoom-in on the planet-planet occultation. The grids discretizing the star and planets are displayed with a low resolution for clarity purpose.

<footnote n="1">¹ Each order is covered by two independent slices on the detectors, which are hereafter treated as individual orders.</footnote>
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_{\star\text{sky}}$ and $y_{\star\text{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}}$ 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/\star}$, by oversampling $z'_{\text{sky}}$ 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/\star}$ 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_{\star\text{sky}}, y_{\star})$ (see Sect. 3.1) is defined as:

$$rv_{\text{surf}/\star} = x_{\star\text{sky}} v_{\text{eq}} \sin i_{\star} \frac{\Omega(y_{\star})}{\Omega_{\text{eq}}} + \sum_{j} c_{\text{cb,j}} \mu^{j}, \quad (2)$$

with $\Omega(y_{\star}) = \Omega_{\text{eq}}(1 - \alpha_{\text{dr}} y_{\star}^2 - \beta_{\text{dr}} y_{\star}^4)$ the stellar rotation rate as a function of latitude (equal to $\Omega_{\text{eq}}$ at the stellar equator), described by a Sun-like differential rotation law, and $c_{\text{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 $\mu = \cos\theta$ accurately in the case of an oblate photosphere, as $\theta$ 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 $\mu$ numerically without approximations, using the projection of the local normal<fnref n="3" /> 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 $\mu$). 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 $\mu$ 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.

<footnote n="2">² We caution that all instances of $(1-f^2)$ in Eq. (14) of Barnes (2009) must be replaced with $(1-f)^2$.</footnote>

<footnote n="3">³ Defined as $(-\partial p/\partial x_{\star\text{sky}}, -\partial p/\partial y_{\star\text{sky}}, 1)$, where the photosphere is described by $z_{\star\text{sky}} = p(x_{\star\text{sky}}, y_{\star\text{sky}})$ following Eq. (13) in Barnes (2009).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Model line profile from the region occulted during a 15 min exposure ($R_p = R_I, R_\star = R_\odot, v \sin i_\star = 25$ km s$^{-1}$). The reference black profile (high precision, 1.5 min oversampling) can be compared with the blue (low precision, rotational broadening accounted for, 1.5 min oversampling), blue dashed (low precision, rotational broadening not accounted for, 1.5 min oversampling), and dotted black (high precision, no oversampling) profiles.</figcaption></figure>

### 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_{\text{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_{\text{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_{\text{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<fnref n="4" /> 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^{\text{meas}}$ and $F^{\text{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^{\text{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).

<footnote n="4">⁴ Provided by the DRS of standard spectrographs or, for activity indexes, by the Data & Analysis Center for Exoplanets (DACE) hosted by the University of Geneva at https://dace.unige.ch/</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Mock ESPRESSO-like time series in the sodium D2 line, for HD 209458b ($R_p = 1.36 R_i$; P = 3.5 d) and a mock companion ($R_p = 2 R_i$; P = 7 d) transiting simultaneously HD 209458 (1.16 $R_\odot$) with $v \sin i_* = 15$ km s$^{-1}$. The transit configuration is illustrated in Fig. 2. Orbital phases are relative to the outer planet. Dotted and dashed lines show the contacts of the inner and outer planets, respectively. Top left: in-transit disk-integrated profile distorted by the occultation of the two planets. Bottom left: double-transit light curve. We note that photometric light curves are typically not an observable from ground-based spectroscopic datasets. Top right: RV derived from a Gaussian fit to the disk-integrated profile series (similar to the output of most spectrographs DRS), blending the well-known RM anomalies from the larger planet on a highly misaligned orbit (~flat deviation) and from the smaller planet on its aligned orbit (S-shaped deviation). The solid blue line show the mock Keplerian stellar reflex motion. Bottom right: 2D map of intrinsic profiles (see Sect. 6.5), showing the tracks of local stellar lines occulted by the planets along their respective transit cores.</figcaption></figure>

## 5.1. Spectral flux calibration

Input spectra $F^{\text{meas}}$ and their associated uncertainties $\sigma^{\text{meas}}$ are generally provided in flux units, related to the number of photoelectrons $N^{\text{meas}}$ cumulated during an exposure per unit of wavelength as:

$$ F^{\text{meas}}(\lambda, E, t) = g_{\text{cal}}(\lambda, E, t)N^{\text{meas}}(\lambda, E, t), \\ \sigma^{\text{meas}}(\lambda, E, t) = g_{\text{cal}}(\lambda, E, t)\sqrt{N^{\text{meas}}(\lambda, E, t)}. \quad (3) $$

At this stage, no resampling has been applied to the spectra and there is no covariance between their bins. The $g_{\text{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_{\text{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 ($\propto 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 $\langle g \rangle_{\text{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).
<figure><img src="image_6.png" /><figcaption>Fig. 5. Flux calibration. Top panel: disks (colored over the rainbow scale as a function of orbital phase, but mostly identical) show the $g_{cal}(\lambda, t)$ values measured at low spectral resolution in a given order, for all exposures of WASP-76 epoch 2. The black dashed line shows the median profile $\langle g \rangle_{cal}$ over both epochs. Dotted vertical lines indicate where the three piece-wise polynomials join. Bottom panel: spectrum measured during one of WASP-76 epoch 2 exposure, scaled back from the DRS calibrated counts (red, with orange uncertainties) to ANTARESS equivalent raw counts (blue, with light blue uncertainties). Plotted uncertainties are amplified by three to visualize the difference in relative errors between the order edge and center.</figcaption></figure>

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 $\delta_{\oplus}$ at the resolution and sampling of the observed spectrum, corrected as:

$$F^{\text{corr}}(\lambda_{\text{B}_{\odot}}, E, t) = \frac{F^{\text{meas}}(\lambda_{\text{B}_{\odot}}, E, t)}{\delta_{\oplus}(\lambda_{\text{B}_{\odot}}, E, t)}. \quad (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 $\delta_{\oplus}$ are also used to exclude telluric lines from CCF masks built with the workflow (Sect. 8.1).

<figure><img src="image_7.png" /><figcaption>Fig. 6. Telluric correction for WASP-76 epoch 1. Top panel: CCF of H₂O over the telluric lines used in the fit, from an observed spectrum before (red) and after (blue) correction, and from its best-fit telluric spectrum (green). Middle panel: Average atmospheric pressures derived for H₂O over the epoch. Bottom panel: Spectrum of the exposure in the top panel, before (red) and after (blue) correction, with its best-fit telluric model (green).</figcaption></figure>

### 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:

$$F^{\text{meas}}(\lambda_{\text{B}_{\odot}}, E, t) = c(\bar{\lambda}_{\text{B}_{\odot}}, E, t) F_{\star}(\lambda_{\text{B}_{\odot}}, E) \delta_{\text{p}}(\lambda_{\text{B}_{\odot}}, E, t). \quad (5)$$

The planetary contribution $\delta_{\text{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 $v$. 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 $\delta_{\text{p}}$ must be re-injected later-on in spectra corrected for flux balance (Sect. 6.3). Importantly, the
<figure><img src="image_8.png" /><figcaption>Fig. 7. Flux balance variations in WASP-76 spectra, between each exposure and their median reference in epoch 1 (top panel) and 2 (middle panel), and between the median reference and their mean over both epochs (bottom panel). Disks show the spectra-to-reference ratios measured at low spectral resolution in each exposure (top and middle panels, colored from purple to red over the rainbow scale with increasing orbital phase) or in each epoch (bottom panel, in blue for epoch 1 and green for epoch 2). Solid lines with matching colors show their best-fit models. Vertical dashed grey lines indicate the central position of one in every eight ESPRESSO slice. The medium-frequency variations captured around slice 120 are observed at high airmass and might correspond to the ozone Chappuis absorption bands at $\lambda = 5750$ and 6030 Å ($v = 52.1$ and $49.7 \times 10^{-10}$ s$^{-1}$).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Data-driven reference spectra of WASP-76 for flux balance corrections in epoch 1 (top panel) and 2 (bottom panel). The dark grey profile is the mean of the median spectra in each epoch (solid black profiles), taken as final reference for the two epochs. Spectra from individual exposures are colored over the rainbow scale with increasing orbital phase.</figcaption></figure>

correction has no impact at short spectral scales. Because high-resolution 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^{\text{meas}}$ spectra to that of a reference spectrum expressed as $F_{\text{ref}}(\lambda, E) = C_{\text{ref}}(\lambda, E)F_{\star}(\lambda, E)$, considering that this reference does not necessarily have the same balance as the true stellar spectrum. An analytical model $R^{\text{theo}}$ describing the flux balance variations between $F^{\text{meas}}$ and $F_{\text{ref}}$ (Fig. 7) is derived in $v$ space as described in Sect. C.3, and applied as:

$$F^{\text{corr}}(\lambda_{\text{B}\odot}, E, t) = F^{\text{meas}}(\lambda_{\text{B}\odot}, E, t) \frac{R^{\text{norm}}(E, t)}{R^{\text{theo}}(\lambda_{\text{B}\odot}, E, t)}, \quad (6)$$

$$\text{where } R^{\text{norm}}(E, t) = \frac{\sum_{\lambda_{\text{B}\odot}} F^{\text{meas}}(\lambda_{\text{B}\odot}, E, t) d\lambda_{\text{B}\odot}}{\sum_{\lambda_{\text{B}\odot}} F_{\text{ref}}(\lambda_{\text{B}\odot}, E) d\lambda_{\text{B}\odot}}.$$

Corrected spectra have the same broadband profile as $F_{\text{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_{\text{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:

$$F^{\text{meas}}(\lambda_{\star}, E, t) - F_{\text{comp}}(\lambda_{\star}, E) > \alpha_{\text{cosm}} \max(\sigma_{F_{\text{comp}}} (\lambda_{\star}, E), \sigma^{\text{meas}} (\lambda_{\star}, E, t)). \quad (7)$$
<figure><img src="image_10.png" /></figure>

**Fig. 9.** Cosmics processing applied to epoch 2 of WASP-76. Top sub-panels show $F(t) - F_{\text{comp}}$ divided by $\sigma_{F_{\text{comp}}}$ (black) and $\sigma^{\text{meas}}(t)$ (green), with $\alpha_{\text{cosm}}$ plotted as a red line, in a sample exposure whose flux spectrum is shown in the bottom subpanels. The red pixel in the top block (together with adjacent pixels) is corrected because it deviates with respect to both its own uncertainty and to the standard deviation over adjacent exposures. The pixel in the bottom block is not corrected because it could result from photon noise in the exposure.

Where $\alpha_{\text{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 $\sigma^{\text{meas}}$ accounts for the processed spectrum possibly being much noisier than adjacent exposures in a given bin – in which case $F^{\text{meas}}(\lambda_{\star}, E, t)$ can deviate significantly from $F_{\text{comp}}(\lambda_{\star}, 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_{\text{comp}}(\lambda_{\star}, E)$ and their variance<fnref n="5" /> to $\sigma_{F_{\text{comp}}}^2(\lambda_{\star}, E)$.

<figure><img src="image_11.png" /></figure>

**Fig. 10.** Persistent peaks masking of a CARMENES time series. Top panel: telluric emission lines (highlighted in blue in the first exposure spectrum, shown in red) were identified through comparison with the stellar continuum (black). Middle panel: zoom on the strongest telluric line from the top panel, plotted over the rainbow scale with increasing orbital phase. This specific line was flagged over the first ten exposures, after which the decreasing airmass makes it negligible. Bottom panel: zoom on one of the strongest telluric lines in the CARMENES range, masked over the entire epoch.

## 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

<footnote n="5">⁵ At this stage of the workflow spectra still display no covariance since they were processed on their own spectral grid.</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Transmission spectra of HD 209458 as a function of light frequency, highlighting different wiggle patterns at the start (exposure 14, top panel) and end (exposure 79, bottom panel) of epoch 2. Top (resp. bottom) subpanels show data before (resp. after) wiggle correction. The wiggle model, overplotted as a dashed red line, is a beat pattern between a dominant (blue) and weaker (orange) component. Right panels show periodograms of the data (horizontal dashed lines are false-alarm probability levels at 1, 5, and 10%), highlighting the removal of these components.</figcaption></figure>

%, 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<fnref n="6" /> 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

<footnote n="6">⁶ The period in wavelength space is $P(\lambda) = \lambda^2/(F(\lambda)c)$.</footnote>
<figure><img src="image_13.png" /><figcaption>Fig. 12. RMS of HD 209458 transmission spectra in epoch 2 before (red) and after (blue) wiggle correction. Green squares are the medi-ans of the propagated error tables. The bottom panel shows the ratio between RMS and median error, expected to be unity for photon-dominated noise. The correction removes the structure in this ratio, which can be attributed to systematic noise from the wiggles.</figcaption></figure>

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 param-eters, 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 iden-tify the spectral ranges that constrain the wiggles and to assess the strength of the two components. Then we sample the chro-matic variations of the components’ frequency and amplitude in each exposure, and determine their optimal polynomial mod-els. 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 tar-gets, 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 respon-sible 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 cor-rected 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 ANALY-SIS 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 spec-tra 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 pro-files for systematic trends identified between the line prop-erties (rv residuals, FWHM, or constrast) and environmen-tal/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 tran-sit (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 spec-tra of WASP-76 and HD 209458’s first epochs for correlations of their CCFs with phase and S/N reduced the out-of-transit disper-sion 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_{*}(t) = \frac{\lambda_{B_{\odot}}}{\left(1 + \frac{rv_{*}/B_{*}(t)}{c}\right)\left(1 + \frac{rv_{B_{*}}/B_{\odot}}{c}\right)} \quad (8) $$

where rv*/*B* is the star Keplerian rv relative to the planetary sys-tem barycenter and rvB*/*B* 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.
<figure><img src="image_14.png" /><figcaption>Fig. 13. Chromatic broadband scaling. Scaling values for WASP-76b exposures (epoch 1) were derived from a chromatic set of BATMAN light curves (see text).</figcaption></figure>

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 rv<sub>B*</sub>/B<sub>☉</sub> 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^{\text{meas}}$ spectra have the low-frequency profile of the chosen reference $F_{\text{ref}}$ but retain their original flux level. We first normalize them to a common integrated flux $TF_{\text{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 $\bar{\lambda}_{\star}$. Final rescaled profiles (see the full procedure in Appendix D.2) can be expressed as:

$$ F^{\text{sc}}(\lambda_{\star}, E, t) = F(\lambda_{\star}, E, t)C_{\text{ref}}^{\text{norm}}(\bar{\lambda}_{\star}, E), \\ \text{with } C_{\text{ref}}^{\text{norm}}(\bar{\lambda}_{\star}, E) = \frac{TF_{\text{ref}}}{\sum_{\lambda_{\star}} F_{\text{ref}}(\lambda_{\star}, E)d\lambda_{\star}} \quad (9) $$

Where the $F$ spectral time series have the low-frequency profile and global flux level of $F_{\text{ref}}$, and the relative flux level of the true spectra.

<figure><img src="image_15.png" /><figcaption>Fig. 14. Numerical broadband scaling. The light curve, generated with ANTARESS stellar grid, corresponds to the mock system in Fig. 2 (middle panel). Note the distorted shape of the light curve due to the misaligned orbit of the planet across the gravity-darkened star, and the asymmetric transit contacts due to its oblateness.</figcaption></figure>

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<fnref n="7" /> 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:

$$ F^{\text{sc}}(\lambda_{\star}, E, t) = C_{\text{ref}}^{\text{norm}}(\bar{\lambda}_{\star}, E)(F_{p}(\lambda_{\star}, E, t) + \sum_{i \notin \text{occ}} f_{i}(\lambda_{\star}, E) S_{i} + \sum_{i \in \text{occ}} f_{i}(\lambda_{\star}, E) (S_{i} - S_{i}^{p}(\lambda_{\star}, E))), \quad (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^{[thick]}(\bar{\lambda}_{\star}, E)$ and by the optically thin atmospheric annulus $S_i^{[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

<footnote n="7">⁷ What we call planet-occulted starlight is actually light from the unocculted stellar profile coming from the equivalent region of the star outside of the transit.</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Map of differential profiles during HD 209458b epoch 1, in the region of the sodium doublet. Values are colored as a function of flux and plotted as a function of wavelengths in the star rest frame (in abscissa) and orbital phase (in ordinate). The planet-occulted stellar lines of sodium and nickel are clearly visible during transit (contacts are shown as white dashed lines). Visual inspection of out-of-transit and in-transit continuum values do not show any strong systematic features.</figcaption></figure>

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):

$$ F^{\text{diff}}(\lambda_{\star}, E, t) = F'_{\star}(\lambda_{\star}, E) - F^{\text{sc}}(\lambda_{\star}, E, t) \\ = C_{\text{ref}}^{\text{norm}}(\bar{\lambda}_{\star}, E)(-F_{p}(\lambda_{\star}, E, t) + \sum_{i \in i_{\text{occ}}} f_{i}(\lambda_{\star}, E) S_{i}^{p}(\lambda_{\star}, E)). \quad (11) $$

In contrast to standard approaches, we calculate a reference for the star, $F'_\star$, 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^{\text{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[\text{thin}]}$ and $F_p$.

## 6.5. Intrinsic profiles

This method converts in-transit differential profiles $F^{\text{diff}}$ into intrinsic stellar profiles $F^{\text{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

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Intrinsic spectra of HD 209458 in epoch 1, showing the sodium doublet from occulted regions of the photosphere. The purple (from end of ingress) and red (from beginning of egress) spectra trace HD 209458b moving from the blueshifted to the redshifted stellar hemisphere along its aligned orbit (spectra are downsampled for clarity, with original data shown as faded profiles).</figcaption></figure>

profiles can be expressed as:

$$ F^{\text{intr}}(\lambda_{\star}, E, t) = C_{\text{ref}}^{\text{norm}}(\bar{\lambda}_{\star}, E) \\ \times \frac{I_{\text{occ}}(\lambda_{\star}, E, t) \alpha_{\text{occ}}(\bar{\lambda}_{\star}, E, t) S_{\text{occ}}^{[\text{thick}]}(\bar{\lambda}_{\star}, E, t)}{1 - lc(\lambda_{\star}, E, t)}. \quad (12) $$

Here, $I_{\text{occ}}$ is the specific intensity emitted by the occulted region, and $\alpha_{\text{occ}}$ its broadband intensity variations. Because the spectral flux scaling $lc$ can be expressed as a function of $\alpha_{\text{occ}}$ and the planetary continuum, $S_{\text{occ}}^{[\text{thick}$]}$, the above expression can be simplified to (for full details see Appendix D.4):

$$ F^{\text{intr}}(\lambda_{\star}, E, t) = F_{\text{ref}}^{\text{norm}}(\bar{\lambda}_{\star}, E) I_{\text{occ}}(\lambda_{\star}, E, t). \quad (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_{\text{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.
<figure><img src="image_18.png" /><figcaption>Fig. 17. Intrinsic CCF map of WASP-76 in epoch 2, plotted as a function of rv in the star rest frame (in abscissa) and orbital phase (in ordinate). Values are colored as a function of normalized flux, unveiling the track of local stellar lines occulted by the planet along its transit chord (green dotted lines show transit contacts). The magenta line shows the planet orbital track, along which was excluded the range contaminated by its atmosphere. CCFs were computed with ANTARESS mask built from WASP-76 data (Sect. 8.1).</figcaption></figure>

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(\lambda_*)$, 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 $\lambda_{surf}$ are calculated for the receiver located at the star rest frame and the local photosphere as the source:

$$ \lambda_{surf}(t) = \frac{\lambda_{*}}{\left(1 + \frac{rv_{surf}(\lambda_{*}, t)}{c}\right)}, \qquad (14) $$

where $rv_{surf}(\lambda_*, 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).

<figure><img src="image_19.png" /><figcaption>Fig. 18. Generation of a CCF mask for HD 209458. Top panel: Master disk-integrated spectrum over all out-of-transit 1D spectra, before (blue) and after (black) normalization with the stellar continuum (red). Shaded blue areas are excluded due to strong telluric contamination. Middle panel: final mask lines, identified in the master spectrum by their minimum (green disks) and surrounding maxima (red disks). Bottom panel: zoom highlighting the range of each mask line (shaded bands) and, when cross-matched with VALD, their origin species.</figcaption></figure>

## 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 $\alpha$-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
<figure><img src="image_20.png" /><figcaption>Fig. 19. Order-per-order weight profiles of disk-integrated stellar spectra in the first (blue) and last (red) exposures of WASP-76 epoch 1 (the lower panel shows a zoom-in of the upper panel). This shows how lower weights are given at the edges of detector orders, at the bottom of telluric lines, and in regions more strongly corrected for Earth's atmospheric diffusion (i.e., toward blue wavelengths and at the start/end of epochs). The stellar spectrum is not included, as it does not contribute to disk-integrated spectra weighing.</figcaption></figure>

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 $\Delta t$ are then defined as (Appendix E.2.2):

$$ 
\begin{aligned}
W^{\text{sc}}(\lambda_{\star}) &= \frac{1}{\sigma^{\text{sc, true}}(\lambda_{\star})^2}, \\
W^{\text{intr}}(\lambda) &= \frac{(1 - lc(\lambda))^2}{\sigma_{\star}^{\text{meas}}(\lambda)^2 + \sigma^{\text{sc, true}}(\lambda)^2},
\end{aligned} 
\quad (15) $$

$$ \text{with } \sigma^{\text{sc, true}}(\lambda_{\star}) = lc(\lambda_{\star}) \sqrt{\tilde{g}_{\text{cal}}(\lambda_{\star}) F_{\star}^{\text{meas}}(\lambda_{\star}) C_{\text{corr}}(\lambda_{\star}) / \Delta t}. $$

The disk-integrated spectra are averaged in the star rest frame, while intrinsic spectra are averaged in the star ($\lambda = \lambda_{\star}$) or photosphere rest frame ($\lambda = \lambda_{\text{surf}}$). The flux calibration profile $\tilde{g}_{\text{cal}}$ (Sect. 5.1), color effect (through $C_{\text{corr}}$, Sect. 5.3), and flux scaling $lc(\lambda_{\star})$ (Sect. 6.3) contribute to low-frequency noise variations, weighing disk-integrated and intrinsic spectra in the same direction through $\sigma^{\text{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_{\text{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). Counterintuitively, $W^{\text{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).

<figure><img src="image_21.png" /><figcaption>Fig. 20. Order-per-order weight profiles of intrinsic stellar spectra near the second (blue) and third (red) transit contacts of WASP-76 in epoch 1, as a function of wavelength in the surface rest frame. Besides the contributions highlighted in Fig. 19, we note the impact of masked planetary ranges (which redshift over the transit as the planet orbits faster than the star rotates) and disk-integrated stellar lines (which blueshift over the transit as the planet first occults the stellar hemisphere rotating toward, and then away, from the observer; see inset).</figcaption></figure>

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{\sum_t F(\lambda, x(t), E)W(\lambda, x(t), E) \Delta x(t)}{\sum_t W(\lambda, x(t), E) \Delta x(t)}. \quad (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_{\text{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_{\text{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
<figure><img src="image_22.png" /><figcaption>**Fig. 21.** Binning of WASP-76 intrinsic CCFs in epoch 1. The top panel shows the $r_{\text{sky}}$ windows covered by each exposure (offset vertically for visibility, and colored over the rainbow scale with increasing orbital phase). Middle panels show the intrinsic CCFs of exposures overlapping with the binning windows chosen as example. The bottom panel shows the resulting binned profiles, corrected for planet contamination and tracing variations in shape across the stellar surface.</figcaption></figure>

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<fnref n="8" />. 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., bisector, 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

<footnote n="8">⁸ Lines of species with differing atomic properties, or different excitation and ionization states, may form at different temperatures in the stellar atmosphere, see e.g., Al Moulla et al. (2022).</footnote>
<figure><img src="image_23.png" /></figure>

**Fig. 22.** Fit to WASP-76 disk-integrated CCFs, computed with the ESPRESSO DRS F9 mask. **Top panel:** the dashed black profile is the Gaussian best-fit to the CCF in the last exposure (blue profile, phase 0.0162) before the stellar line is too contaminated by the planetary signal (red area) to be accurately fitted. Blue areas are never contaminated by the planet, and used to measured the continuum. **Bottom panel:** rv residuals from the Keplerian model. Orange squares and green disks were derived from CCFs with no strong planetary contamination in epochs 1 and 2. Matching black symbols (shown without errors for clarity) were derived without excluding planetary ranges from the fits.

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 ($\xi$ = 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;

<figure><img src="image_24.png" /></figure>

**Fig. 23.** Fit to the disk-integrated sodium doublet of HD 209458. **Top panel:** the black profile shows the best-fit PYSME theoretical spectrum to the measured spectrum (blue profile) over the unshaded regions. Only the sodium abundance was varied in this example, so that other deep stellar lines were excluded from the fit. **Bottom panel:** residuals from the fit.

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
<figure><img src="image_25.png" /><figcaption>**Fig. 24.** Average line profile of HD 209458 in the white CCFs calculated with ESPRESSO F9 mask (epoch 2). Exposures at orbital phases -0.014 (ingress, blue profile), 0 (conjunction, green profile), and 0.014 (egress, red profile) illustrate the change in line shape between stellar limb and center. Because of the aligned orbit, profiles occulted at ingress and egress have opposite *rv* and similar shapes. Dashed black lines show the Gaussian profiles used to derive the line properties (Fig. 25).</figcaption></figure>

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 contaminated 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
<figure><img src="image_26.png" /><figcaption>**Fig. 25.** Properties of stellar regions occulted by HD 209458b and WASP-76b. Dotted vertical lines show transit contacts. Measured values (first epochs in blue, second epochs in red) are derived from best-fit profiles to each intrinsic CCF, before convolution with ESPRESSO LSF (Sect. 8.2.2). They are thus comparable to the models (dashed black curves) derived from fits to intrinsic CCF series (Sect. 8.3.1). The black curve in WASP-76 surface *rv* is the model from Ehrenreich et al. (2020). The two outliers in contrast at the limbs of HD 209458 are consistent with the FWHM and *rv* series, and thus kept in the analysis. WASP-76 local stellar lines are fully contaminated by the planet over the shaded range.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. 26.** Sodium doublet properties of stellar regions occulted by HD 209458b (same format as Fig. 25). A common model to both epochs (dashed orange curves) is preferred for all properties. The photospheric *rv* model is shown as a solid black curve.</figcaption></figure>

(sky-projected $\lambda$ or 3D $\psi$ 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 $\lambda$ 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 $\chi^2 = 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 \pm 0.007 \text{ km s}^{-1}$ with a similar precision as Casasayas-Barris et al. (2021) ($\lambda = 1.58 \pm 0.08^{\circ}$ and $v_{eq} \sin i_* = 4.228 \pm 0.007 \text{ km s}^{-1}$), 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,
<figure><img src="image_28.png" /><figcaption>**Fig. 27.** Intrinsic sodium doublet of HD 209458, for the joint fit to both epochs. Top panel: measured (blue) and model (black) profiles at mid-transit in epoch 2. Original data (semi-transparent spectrum) was downsampled (solid spectrum) for clarity. Shaded regions were excluded from the fit. Middle panel: residual spectrum. Bottom panel: zoom on the D1 line, illustrating center-to-limb variations from the NLTE pySME/MARCS model (profiles are colored in purple, green, and red at orbital phases -0.0126, -0.0002, and 0.0123).</figcaption></figure>

considering that more realistic error bars on $\lambda$ and $v_{eq} \sin i_\star$ (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 $\phi \sim 0$, see Eq. (2)):

$$ \begin{aligned} rv_{\text{surf}/\star}(\phi) &= v_{\text{eq}} \sin i_{\star} a_p / R_{\star} \cos(\lambda) 2\pi\phi \\ \dot{rv}_{\text{surf}/\star}(\phi) &= v_{\text{eq}} \sin i_{\star} a_p / R_{\star} \cos(\lambda) 2\pi. \end{aligned} \qquad (17) $$

Thus, $v_{eq} \sin i_\star \cos(\lambda)$ is constrained in the same way by the slope and level of the rv series and there is a strong degeneracy between $\lambda$ and $v_{eq} \sin i_\star$ (Fig. F.2). This is in part why we derive $(\lambda = -37.1^{+12.6}_{-20.6}^\circ$; $v_{eq} \sin i_\star = 0.86^{+0.13}_{-0.19}$ km s⁻¹) while Ehrenreich et al. (2020) derive $(\lambda = 61 \pm 7^\circ$; $v_{eq} \sin i_\star = 1.5 \pm 0.3$ km s⁻¹), since both yield $v_{eq} \sin i_\star \cos(\lambda) \sim 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 ($\lambda$ in between -180 and 0°, Fig. 29) than the redshifted one ($\lambda$ 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_\star$ (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_\star = 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 $\lambda$ and $v_{eq} \sin i_\star$ we applied a RMR fit with $a_p/R_\star$ and $i_p$ free to vary, using photometry priors from Ehrenreich et al. (2020) ($a_p/R_\star = 4.08^{+0.02}_{-0.06}$, $i_p = 89.623^{+0.005}_{-0.034}$). We derived values consistent with the priors ($a_p/R_\star = 4.059^{+0.044}_{-0.033}$; $i_p = 89.604^{+0.024}_{-0.015}$), 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.
<figure><img src="image_29.png" /><figcaption>**Fig. 28.** RMR fits of HD 209458b (shown for epoch 1) and WASP-76b (shown for epoch 2). Maps are plotted as a function of *rv* in the star rest frame (in abscissa) and orbital phase (in ordinate). Green dotted lines show transit contacts. The masked streak in WASP-76 maps is contaminated by the planetary atmosphere. **Top panels:** CCFs of intrinsic spectra, colored with flux. Green solid lines show the best stellar surface *rv* models. **Bottom panels:** residuals over the full epochs. Out-of-transit values show the CCFs of differential spectra, and in-transit values the difference between the CCFs of intrinsic spectra and their best-fit models (scaled to the out-of-transit level). The slanted in-transit features aligned with the surface *rv* model of HD 209458 may arise from stellar lines close to the CCF mask lines.</figcaption></figure>

**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.
<figure><img src="image_30.png" /><figcaption>**Fig. 29.** Sky-projected view of WASP-76 for the best-fit orbital architecture (see plot details in Fig. 2). Thin lines surrounding the orbital trajectory (thick curve) show orbits obtained for orbital inclination, semi-major axis, and sky-projected obliquity values drawn randomly within $1\sigma$ from their probability distributions.</figcaption></figure>

- 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

<figure><img src="image_31.png" /><figcaption>**Fig. 30.** Joint fit to HD 209458 intrinsic sodium doublet, shown for epoch 2 as a function of wavelength in the star rest frame (in abscissa) and orbital phase (in ordinate). Green dotted lines are transit contacts. *Top panel:* intrinsic spectra, colored with flux. Stellar lines kept to solar abundance were excluded from the fit (shaded regions). *Middle panel:* best-fit theoretical model. *Bottom panels:* residuals over the full epoch, calculated as in Fig. 28.</figcaption></figure>

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); SCI PY, Virtanen et al. (2020).

## References

* Al Moulla, K., Dumusque, X., Cretignier, M., Zhao, Y., & Valenti, J. A. 2022, A&A, 664, A34

* Albrecht, S. H., Dawson, R. I., & Winn, J. N. 2022, PASP, 134, 082001

* Allart, R., Lovis, C., Pino, L., et al. 2017, A&A, 606, A144

* Allart, R., Bourrier, V., Lovis, C., et al. 2018, Science, 362, 1384

* Allart, R., Pino, L., Lovis, C., et al. 2020, A&A, 644, A155

* Allart, R., Lovis, C., Faria, J., et al. 2022, A&A, 666, A196

* Amarsi, A. M. 2022, https://doi.org/10.5281/zenodo.7088951

* Amarsi, A. M., Lind, K., Osorio, Y., et al. 2020, A&A, 642, A62

* Amarsi, A. M., Liljegren, S., & Nissen, P. E. 2022, A&A, 668, A68

* Artigau, É., Astudillo-Defru, N., Delfosse, X., et al. 2014, SPIE Conf. Ser., 9149, 914905

* Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

* Azevedo Silva, T., Demangeon, O. D. S., Santos, N. C., et al. 2022, A&A, 666, A10

* Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373

* Barnes, J. W. 2009, ApJ, 705, 683

* Benatti, S., Damasso, M., Borsa, F., et al. 2021, A&A, 650, A66

* Bertaux, J. L., Lallement, R., Ferron, S., Boonne, C., & Bodichon, R. 2014, A&A, 564, A46

* Birkby, J. L., de Kok, R. J., Brogi, M., et al. 2013, MNRAS, 436, L35

* Bonomo, A. S., Desidera, S., Benatti, S., et al. 2017, A&A, 602, A107

* Borsa, F., & Zannoni, A. 2018, A&A, 617, A134

* Borsa, F., Allart, R., Casasayas-Barris, N., et al. 2021, A&A, 645, A24

* Bouchy, F., Pepe, F., & Queloz, D. 2001, A&A, 374, 733

* Bouchy, F., Doyon, R., Artigau, É., et al. 2017, The Messenger, 169, 21

* Boué, G., Montalto, M., Boisse, I., Oshagh, M., & Santos, N. C. 2013, A&A, 550, A53

* Bourrier, V., & Hébrard, G. 2014, A&A, 569, A65

* Bourrier, V., Lecavelier des Etangs, A., Hébrard, G., et al. 2015, A&A, 579, A55

* Bourrier, V., Ehrenreich, D., King, G., et al. 2017, A&A, 597, A26

* Bourrier, V., Lovis, C., Beust, H., et al. 2018, Nature, 553, 477

* Bourrier, V., Ehrenreich, D., Lendl, M., et al. 2020, A&A, 635, A205

* Bourrier, V., Lovis, C., Cretignier, M., et al. 2021, A&A, 654, A152

* Bourrier, V., Zapatero Osorio, M. R., Allart, R., et al. 2022, A&A, 663, A160

* Bourrier, V., Attia, O., Mallonn, M., et al. 2023, A&A, 669, A63

* Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2012, Nature, 486, 502

* Bruls, J. H. M. J., & Rutten, R. J. 1992, A&A, 265, 257

* Carteret, Y., Bourrier, V., & Dethier, W. 2024, A&A, 683, A63
Casasayas-Barris, N., Pallé, É., Yan, F., et al. 2019, A&A, 628, A9

Casasayas-Barris, N., Pallé, É., Yan, F., et al. 2020, A&A, 635, A206

Casasayas-Barris, N., Palle, E., Stangret, M., et al. 2021, A&A, 647, A26

Casasayas-Barris, N., Borsa, F., Palle, E., et al. 2022, A&A, 664, A121

Cegla, H. M., Lovis, C., Bourrier, V., et al. 2016, A&A, 588, A127

Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. 2002, ApJ, 568, 377

Claudi, R., Benatti, S., Carleo, I., et al. 2018, SPIE Conf. Ser., 10702, 107020Z

Collier Cameron, A., Bruce, V. A., Miller, G. R. M., Triaud, A. H. M. J., & Queloz, D. 2010, MNRAS, 403, 151

Cook, N. J., Artigau, É., Doyon, R., et al. 2022, PASP, 134, 114509

Cosentino, R., Lovis, C., Pepe, F., et al. 2012, SPIE Conf. Ser., 8446

Cretignier, M., Dumusque, X., Allart, R., Pepe, F., & Lovis, C. 2020a, A&A, 633, A76

Cretignier, M., Francfort, J., Dumusque, X., Allart, R., & Pepe, F. 2020b, A&A, 640, A42

Cretignier, M., Dumusque, X., Hara, N. C., & Pepe, F. 2021, A&A, 653, A43

Cristo, E., Santos, N. C., Demangeon, O., et al. 2022, A&A, 660, A52

Czesla, S., Klocová, T., Khalafinejad, S., Wolter, U., & Schmitt, J. H. M. M. 2015, A&A, 582, A51

Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, PyA: Python astronomy-related packages

Dekker, H., D'Odorico, S., Kaufer, A., Delabre, B., & Kotzlowski, H. 2000, SPIE Conf. Ser., 4008, 534

Deming, D., Wilkins, A., McCullough, P., et al. 2013, ApJ, 774, 95

Dethier, W., & Bourrier, V. 2023, A&A, 674, A86

Donati, J. F., Kouach, D., Moutou, C., et al. 2020, MNRAS, 498, 5684

Dravins, D., Ludwig, H.-G., Dahlén, E., & Pazira, H. 2017, A&A, 605, A90

Dravins, D., Ludwig, H.-G., & Freytag, B. 2021, A&A, 649, A17

Dumusque, X. 2018, A&A, 620, A47

Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, Nature, 580, 597

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fu, G., Deming, D., Lothringer, J., et al. 2021, AJ, 162, 108

Gandhi, S., Kesseli, A., Snellen, I., et al. 2022, MNRAS, 515, 749

Giménez, A. 2006, ApJ, 650, 408

Gray, D. F. 1975, ApJ, 202, 148

Gray, D. F. 2021, The Observation and Analysis of Stellar Photospheres (Cambridge University Press)

Guilluy, G., Bourrier, V., Jaziri, Y., et al. 2023, A&A, 676, A130

Gullikson, K., Dodson-Robinson, S., & Kraus, A. 2014, AJ, 148, 53

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hirano, T., Narita, N., Sato, B., et al. 2011, PASJ, 63, L57

Hoeijmakers, H. J., Seidel, J. V., Pino, L., et al. 2020, A&A, 641, A123

Hunter, J. 2007, Comput. Sci. Eng., 9, 90

Jurgenson, C., Fischer, D., McCracken, T., et al. 2016, SPIE Conf. Ser., 9908, 99086T

Kausch, W., Noll, S., Smette, A., et al. 2015, A&A, 576, A78

Kesseli, A. Y., Snellen, I. A. G., Casasayas-Barris, N., Mollière, P., & Sánchez-López, A. 2022, AJ, 163, 107

Kreidberg, L. 2015, PASP, 127, 1161

Kupka, F. G., Ryabchikova, T. A., Piskunov, N. E., Stempels, H. C., & Weiss, W. W. 2000, Baltic Astron., 9, 590

Lockwood, A. C., Johnson, J. A., Bender, C. F., et al. 2014, ApJ, 783, L29

Louden, T., & Wheatley, P. J. 2015, ApJ, 814, L24

Maguire, C., Gibson, N. P., Nugroho, S. K., et al. 2023, MNRAS, 519, 1030

Mayor, M., & Queloz, D. 1995, Nature, 378, 355

Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20

Mazeh, T., Tamuz, O., & Zucker, S. 2007, in Transiting Extrapolar Planets Workshop, eds. C. Afonso, D. Weldrake, & T. Henning, Astronomical Society of the Pacific Conference Series, 366, 119

McKerns, M. M., Strand, L., Sullivan, T., Fang, A., & Aivazis, M. A. G. 2012, Building a Framework for Predictive Science

McLaughlin, D. B. 1924, ApJ, 60, 22

Merritt, S. R., Gibson, N. P., Nugroho, S. K., et al. 2020, A&A, 636, A117

Mounzer, D., Lovis, C., Seidel, J. V., et al. 2022, A&A, 668, A1

Newville, M., Stensitzki, T., Allen, D. B., et al. 2016, Lmfit: Non-Linear Least-Square Minimization and Curve-Fitting for Python, Astrophysics Source Code Library, [record ascl:1606.014]

Nugroho, S. K., Gibson, N. P., de Mooij, E. J. W., et al. 2020, ApJ, 898, L31

Ohta, Y., Taruya, A., & Suto, Y. 2005, ApJ, 622, 1118

Pepe, F., Ehrenreich, D., & Meyer, M. R. 2014, Nature, 513, 358

Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96

Perruchot, S., Kohler, D., Bouchy, F., et al. 2008, SPIE Conf. Ser., 7014, 70140J

Pino, L., Désert, J.-M., Brogi, M., et al. 2020, ApJ, 894, L27

Piskunov, N. E., Kupka, F., Ryabchikova, T. A., Weiss, W. W., & Jeffery, C. S. 1995, A&AS, 112, 525

Plez, B. 2012, Turbospectrum: Code for spectral synthesis, Astrophysics Source Code Library, [record ascl:1205.004]

Queloz, D., Eggenberger, A., Mayor, M., et al. 2000, A&A, 359, L13

Queloz, D., Mayor, M., Udry, S., et al. 2001, The Messenger, 105, 1

Quirrenbach, A., Amado, P. J., Caballero, J. A., et al. 2016, SPIE Conf. Ser., 9908, 990812

Redfield, S., Endl, M., Cochran, W. D., & Koesterke, L. 2008, ApJ, 673, L87

Rossiter, R. A. 1924, ApJ, 60, 15

Rothman, L. S. 2021, Nat. Rev. Phys., 3, 302

Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr, 90, 054005

Seager, S., & Sasselov, D. D. 2000, ApJ, 537, 916

Sedaghati, E., MacDonald, R. J., Casasayas-Barris, N., et al. 2021, MNRAS, 505, 435

Seidel, J. V., Ehrenreich, D., Wyttenbach, A., et al. 2019, A&A, 623, A166

Seidel, J. V., Ehrenreich, D., Pino, L., et al. 2020, A&A, 633, A86

Seidel, J. V., Ehrenreich, D., Allart, R., et al. 2021, A&A, 653, A73

Seidel, J. V., Cegla, H. M., Doyle, L., et al. 2022, MNRAS, 513, L15

Seidel, J. V., Borsa, F., Pino, L., et al. 2023, A&A, 673, A125

Seifahrt, A., Stürmer, J., Bean, J. L., & Schwab, C. 2018, SPIE Conf. Ser., 10702, 107026D

Sicilia, D., Malavolta, L., Pino, L., et al. 2022, A&A, 667, A19

Sing, D. K., Vidal-Madjar, A., Désert, J.-M., Lecavelier des Etangs, A., & Ballester, G. 2008, ApJ, 686, 658

Sing, D. K., Wakeford, H. R., Showman, A. P., et al. 2015, MNRAS, 446, 2428

Sing, D. K., Lavvas, P., Ballester, G. E., et al. 2019, AJ, 158, 91

Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77

Snellen, I. A. G., Albrecht, S., de Mooij, E. J. W., & Le Poole, R. S. 2008, A&A, 487, 357

Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049

Snellen, I., Brandl, B., de Kok, R., et al. 2014, arXiv e-prints [arXiv:1404.7506]

Stevenson, K. B., Désert, J.-M., Line, M. R., et al. 2014, Science, 346, 838

Tabernero, H. M., Zapatero Osorio, M. R., Allart, R., et al. 2021, A&A, 646, A158

Takeda, Y., & UeNo, S. 2017, PASJ, 69, 46

Temple, L. Y., Hellier, C., Anderson, D. R., et al. 2019, MNRAS, 490, 2467

The pandas development team. 2020. https://doi.org/10.5281/zenodo.3509134

Thompson, A. P. G., Watson, C. A., Haywood, R. D., et al. 2020, MNRAS, 494, 4279

Torres, G., Winn, J. N., & Holman, M. J. 2008, ApJ, 677, 1324

Vidal-Madjar, A., Ferlet, R., Gry, C., & Lallement, R. 1986, A&A, 155, 407

Vidal-Madjar, A., Lecavelier des Etangs, A., Désert, J.-M., et al. 2003, Nature, 422, 143

Vidal-Madjar, A., Sing, D. K., Lecavelier Des Etangs, A., et al. 2011, A&A, 527, A110

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261

Wehrhahn, A., Piskunov, N., & Ryabchikova, T. 2023, A&A, 671, A171

Wyttenbach, A., Ehrenreich, D., Lovis, C., Udry, S., & Pepe, F. 2015, A&A, 577, A62

Wyttenbach, A., Lovis, C., Ehrenreich, D., et al. 2017, A&A, 602, A36

Wyttenbach, A., Mollière, P., Ehrenreich, D., et al. 2020, A&A, 638, A87

Yan, F., Pallé, E., Fosbury, R. A. E., Petr-Gotzens, M. G., & Henning, T. 2017, A&A, 603, A73

Yan, F., Pallé, E., Reiners, A., et al. 2020, A&A, 640, A5

Zhang, Y., Snellen, I. A. G., Wyttenbach, A., et al. 2022, A&A, 666, A47
<figure><img src="image_32.png" /><figcaption>Fig. A.1. Sketch of correlated pixels after resampling a spectrum with BINDENSITY, depending on the chosen interpolator and the new bins size. We highlight in orange all the pixels which are correlated with a reference pixel (marked in blue). The original spectrum is assumed to be uncorrelated.</figcaption></figure>

# 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<fnref n="9" />. 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^2$ (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^2$ 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^2$ 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^1$ 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 $\theta$:

$$ \ln p(\theta|d, M) = \ln \bar{\sigma}(\theta) + \ln \mathcal{L}(\theta), \qquad (A.1) $$

where $m$ represents the family of models fitted to the data $d$, $\bar{\sigma}(\theta) = p(\theta|m)$ the prior probability of the model parameters, and $\mathcal{L}(\theta) = p(d|m, \theta)$ the likelihood. Under the assumption of Gaussian noise we have:

$$ \ln \mathcal{L}(\theta) = - \ln(\sqrt{\det(2\pi\Sigma)}) - \frac{1}{2} r' \Sigma^{-1} r, \qquad (A.2) $$

where $r = d - m(\theta)$ is the vector of residuals and $\Sigma$ 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) = - \sum_{n} \ln(\sqrt{2\pi L[n, n]}) - \frac{1}{2} u' u. \qquad (A.3) $$

Here, $L$ is the lower triangular matrix satisfying $C = LL'$, and $u = L^{-1}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 $\Sigma$ reduces to the variance $\sigma^2$ along its diagonal and

$$ \ln \mathcal{L}(\theta) = - \sum_{n} \ln(\sqrt{2\pi\sigma(n)}) - \frac{1}{2} \sum_{n} \left( \frac{r(n)}{\sigma(n)} \right)^2. \qquad (A.4) $$

The MCMC approach implemented in the ANTARESS pipeline samples the posterior probability distributions (PDF)

<footnote n="9">⁹ The BINDENSITY package is installable with pip and conda and the sources are publicly available at https://gitlab.unige.ch/jean-baptiste.delisle/bindensity</footnote>
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{aligned}
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{aligned} 
\qquad \text{(B.1)}
$$

Here, *E* is the eccentricity anomaly, derived numerically at a given orbital phase (or time *t*). The corresponding orbital velocity coordinates are:

$$ 
\begin{aligned}
x'_{\text{orb}}(t) &= -\frac{2\pi}{P} \left(\frac{a}{R_{\star}}\right)^2 \frac{1}{D_{\text{p}}(t)} \sin(E(t)), \\
y'_{\text{orb}}(t) &= \frac{2\pi}{P} \left(\frac{a}{R_{\star}}\right)^2 \frac{1}{D_{\text{p}}(t)} \sqrt{1 - e^2} \cos(E(t)),
\end{aligned} 
\qquad \text{(B.2)}
$$

with $D_{\text{p}}(t) = \frac{a}{R_{\star}}(1 - e \cos(E(t)))$.

Rotating the orbital coordinate system around the orbital plane normal by the argument of periastron $\omega$, so that *x* becomes the node line, yields:

$$ 
\begin{aligned}
x_{\text{LOS}}(t) &= -x_{\text{orb}}(t) \cos(\omega) + y_{\text{orb}}(t) \sin(\omega), \\
y_{\text{LOS}}(t) &= x_{\text{orb}}(t) \sin(\omega) + y_{\text{orb}}(t) \cos(\omega).
\end{aligned} 
\qquad \text{(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{aligned}
x_{\text{sky}}(t) &= x_{\text{LOS}}(t), \\
y_{\text{sky}}(t) &= -y_{\text{LOS}}(t) \cos(i_p), \\
z_{\text{sky}}(t) &= y_{\text{LOS}}(t) \sin(i_p).
\end{aligned} 
\qquad \text{(B.4)}
$$

Rotating the latter reference frame around the LOS by the sky-projected obliquity $\lambda$ 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{aligned}
x_{\star\text{sky}}(t) &= x_{\text{sky}}(t) \cos(\lambda) - y_{\text{sky}}(t) \sin(\lambda), \\
y_{\star\text{sky}}(t) &= x_{\text{sky}}(t) \sin(\lambda) + y_{\text{sky}}(t) \cos(\lambda), \\
z_{\star\text{sky}}(t) &= z_{\text{sky}}(t).
\end{aligned} 
\qquad \text{(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_{\star\text{sky}}$ is calculated as $\sqrt{1 - x_{\star\text{sky}}^2 - x_{\star\text{sky}}^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_\star$, 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{aligned}
x_{\star}(t) &= x_{\star\text{sky}}(t), \\
y_{\star}(t) &= y_{\star\text{sky}}(t) \sin(i_{\star}) + z_{\star\text{sky}}(t) \cos(i_{\star}), \\
z_{\star}(t) &= -y_{\star\text{sky}}(t) \cos(i_{\star}) + z_{\star\text{sky}}(t) \sin(i_{\star}).
\end{aligned} 
\qquad \text{(B.6)}
$$

Another useful coordinate is the Keplerian rv of the star, relative to the barycenter of the planetary system:

$$ 
\text{rv}_{\star/B_{\star}}(t) = \sum_{\text{body p}} K_{\text{p}} [\cos(\theta_{\text{p}}(t) + \omega_{\text{p}}) + \text{e} \cos(\omega_{\text{p}})], \qquad \text{(B.7)}
$$

where $K_{\text{p}}$ is the Keplerian rv semi-amplitude induced by a given body orbiting the star and $\theta_{\text{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,\bar{\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 dependence 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 $\bar{t}$ are approximated as the brigitthness-weighted average of all $x_i$:

$$ x(\bar{\lambda}, \bar{t}) = \frac{\sum_{k,p} \sum_{i_{k,\bar{\lambda}}} x_i I_i(\bar{\lambda}) S_i}{\sum_{k,p} \sum_{i_{k,\bar{\lambda}}} I_i(\bar{\lambda}) S_i}, \qquad \text{(B.8)} $$

where $I_i(\bar{\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(\bar{\lambda}, \bar{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($\bar{\lambda}, \bar{t}$).

- High: the line profile of each oversampled region is calculated as $F_k(\lambda) = \sum_{i_k} F_i(\lambda)$, where $F_i(\lambda) = I(x_i, \lambda)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.
<figure><img src="image_33.png" /><figcaption>Fig. B.1. Main ANTARESS reference frames (green axis, orbital; blue axis, sky-projected orbital; red axis, sky-projected stellar; black, stellar). The star is at the origin. An eccentric planetary orbit is shown as a green curve, oriented in the direction of the planet motion.</figcaption></figure>

# 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_{cal}(\bar{\lambda}, E, t) = \frac{\sum_{\lambda} \sigma^{\text{meas}}(\lambda, E, t)^2}{\sum_{\lambda} F^{\text{meas}}(\lambda, E, t)}. \quad (C.1)$$

Values of $g_{cal}(\bar{\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 $\langle g \rangle_{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$_2$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_{\text{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^{\text{meas}}$ between the binned reference and exposure spectra is fitted with a smoothing spline or polynomial $R^{\text{theo}}$. A preliminary fit is performed to automatically identify and exclude outliers, and the final fitted model (Fig. 7) is an estimate of:

$$ R^{\text{meas}}(\bar{v}_{\text{B}_{\odot}}, E, t) = \frac{N^{\text{meas}}(\bar{v}_{\text{B}_{\odot}}, E, t)}{N_{\text{ref}}(\bar{v}_{\text{B}_{\odot}}, E)} = \frac{\delta_{\text{p}}(\bar{v}_{\text{B}_{\odot}}, E, t) c(\bar{v}_{\text{B}_{\odot}}, E, t)}{C_{\text{ref}}(\bar{v}_{\text{B}_{\odot}}, E)}. \quad (C.2) $$

Here, we assume that the reference spectrum and the true stellar spectrum deviate at low resolution, so that $C_{\text{ref}}$ and $N_{\star}$ remain separated when binning $N_{\text{ref}} = C_{\text{ref}}N_{\star}$. Light frequency $v$ is used as fit coordinate to stretch the data and its variations over a more regular grid. Finally, we extrapolate $R^{\text{theo}}(\bar{v}_{\text{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^{\text{corr}}/F_{\text{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 $\delta_{\text{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_{\text{comp}}(\lambda_{\star}, E)$ and standard deviation $\sigma_{F_{\text{comp}}}(\lambda_{\star}, 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(v, t) = 1 + \sum_{k} A_{k}(v, t) \sin(2\pi \int (F_{k}(v, t)dv) - \Phi_{k}(t)). \quad (C.3) $$

The amplitude $A_k$ and frequency$^{10}$ $F_k$ of each component show fewer variations as a function of light frequency $v = c/\lambda$ than

<footnote>$^{10}$ We highlight that the phase of the sine term is calculated as the integral of the frequency $F_k(v)$. This is because the rate at which a sinusoid phase changes is equal to the frequency at which this sinusoid oscillates, in other words a sinusoid frequency is the derivative of its phase.</footnote>
wavelength, and are described as:

$$A_k(v, t) = \sum_{i=0}^{da,k} a_{\text{chrom,k,i}}(t)(v - v_{\text{ref}})^i, \qquad (C.4)$$

$$F_k(v, t) = \sum_{i=0}^{df,k} f_{\text{chrom,k,i}}(t)(v - v_{\text{ref}})^i,$$

where the polynomials are defined relative to $v_{\text{ref}}$, corresponding to 6000 Å, to decrease leverage in the fits. The chromatic amplitude and frequency coefficients $a_{\text{chrom,k,i}}(t)$ and $f_{\text{chrom,k,i}}(t)$, as well as the phase reference $\Phi_k(t)$, are expressed as a linear combination of the cartesian pointing coordinates $g(\{\alpha_{\text{point,k}}\}, t) =$

$$\begin{array}{ll}
\alpha_{\text{u/west}} + \alpha_{\text{x}}x(t) + \alpha_{\text{y/west}}y(t) + \alpha_{\text{z/west}}z(t) & \text{if } \theta(t) > \theta_{\text{mer}}, \\
\alpha_{\text{u/east}} + \alpha_{\text{x}}x(t) + \alpha_{\text{y/east}}y(t) + \alpha_{\text{z/east}}z(t) & \text{if } \theta(t) \le \theta_{\text{mer}},
\end{array}
$$
with $x(t) = \sin(\theta(t))$, \\
$y(t) = \cos(\theta(t))$, \\
$z(t) = \sin(\phi(t))$, \\
(C.5)

where $\theta$ and $\phi$ are the azimuth and altitude of the telescope. We found that the model coefficients change when the LOS crosses the meridian ($\theta_{\text{mer}} = 180^{\circ}$), but that the model remains continuous in value ($\alpha_{\text{y/west}} = \alpha_{\text{y/east}} + \sin(\phi_{\text{mer}})(\alpha_{\text{z/west}} - \alpha_{\text{z/east}}) - (\alpha_{\text{u/west}} - \alpha_{\text{u/east}})$) and derivative ($\alpha_{\text{x/west}} = \alpha_{\text{x/east}} = \alpha_{\text{x}}$). The full model is thus controlled by a set of pointing parameters $\{a_{\text{point,k,i}}\}$, $\{f_{\text{point,k,i}}\}$, and $\{\Phi_{\text{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 $v$ (Fig. C.3). This allows determining their degree and deriving guess values for the chromatic coefficients $a_{\text{chrom,k,i}}$ and $f_{\text{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 $v$, 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_{\text{chrom,k,i}}(t)$, $f_{\text{chrom,k,i}}(t)$, and $\Phi_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_{\text{point,k,i}}\}$, $\{f_{\text{point,k,i}}\}$, and $\{\Phi_{\text{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 cross-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
<figure><img src="image_34.png" /><figcaption>Fig. C.1. Periodogram cumulated over all transmission spectra in HD 209458b epoch 2, before (top panel) and after (bottom panel) wiggle correction. The two components included in the model are highlighted. Although a peak to the left of the first component may be removed by including a third component to the model, we emphasize the change in power scale between the two panels.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. C.2. Chromatic sampling, illustrated with exposure at index 14 in HD 209458b epoch 2. Subpanels show the processed transmission spectrum for consecutive positions of the sampling window, with the periodogram used to determine the local frequency $A_k(v_{window})$ of the considered wiggle component (the unhatched range is set by the user to limit the search around the frequency expected from the screening step). The red model shows the sine model fitted at this frequency. The blue and orange boxes respectively show the sampling of the first and second wiggle components, whose frequencies are distinct enough that they can be processed successively. The width of the sliding window is adjusted to cover several cycles of the considered wiggle component while remaining small enough that its frequency is constant over the band.</figcaption></figure>

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{aligned}
F^{\text{meas}}(\bar{\lambda}_{\star}, E, t) &= F_{\text{ref}}(\bar{\lambda}_{\star}, E)R^{\text{norm}}(E, t) \\
&= C_{\text{ref}}(\bar{\lambda}_{\star}, E)F_{\star}(\bar{\lambda}_{\star}, E)R^{\text{norm}}(E, t).
\end{aligned} 
\quad (\text{D.1})
$$
<figure><img src="image_36.png" /><figcaption>**Fig. C.3.** Chromatic analysis. Red points in the top subpanels show measurements from the sampling step in Fig. C.2, for the amplitude (top panels) and frequency (bottom panels) of the first (blue box) and second (orange box) wiggle components. Bottom subpanels show residuals between the measurements and their best-fit model (solid black curves). Blue points are outliers automatically excluded from the fit through sigma-clipping. The wave-like variations in the measurements are an artifact of the sampling method.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>**Fig. C.4.** Pointing analysis. Diamonds (resp. squares) show the chromatic amplitude (top), frequency (middle) and phase (bottom) coefficients derived east (resp. west) of the meridian from the fit to individual exposures, plotted as a function of orbital phase. The meridian crossing is indicated as a vertical dotted line. Bottom subpanels show residuals between the measurements and their best-fit pointing coordinate model (dashed black curves). The blue and orange boxes correspond to the first and second wiggle components.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>Fig. C.5. Same as Fig. C.4, for the dominant wiggle component in WASP 76 epoch 2. These plots illustrate the reset of pointing coefficients following a change in guide star, indicated as a vertical dashed line.</figcaption></figure>

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 Rnorm (ratio between the total flux of the spectrum and that of the reference, Eq. 6) through multiplication by

$$C_F(E, t) = \frac{TF_{\text{ref}}}{\sum_{\lambda_*} F^{\text{meas}}(\lambda_*, E, t) d\lambda_*} \quad (\text{D.2})$$

The integrated flux of the resulting profiles is thus normalized to TFref, which can either be set to the median flux over the epoch < t<, Σλ<, 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 Rnorm. 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_*, E, t) F^{\text{meas}}(\lambda_*, E, t) = F(\lambda_*, E, t). \quad (\text{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(\bar{\lambda}_*, E, t) C_{\text{ref}}(\bar{\lambda}_*, E) \sum_{\lambda_*} F_*(\lambda_*, E) d\lambda_* = \sum_{\lambda_*} F(\lambda_*, E, t) d\lambda_* \quad (\text{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(\bar{\lambda}_*, E, t) = \frac{\sum_{\lambda_*} F(\lambda_*, E, t) d\lambda_*}{\sum_{\lambda_*} F_*(\lambda_*, E) d\lambda_*} \quad (\text{D.5})$$

The scaling light curve lc(λ*, E, t) can be defined chromatically, namely, over a set of bands at λ*, 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 photospheric 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 δp(λ*, 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(λ*, E, t) over the high-resolution spectral grid of each exposure to calculate the rescaled spectra as:

$$ \begin{align} F^{\text{sc}}(\lambda_*, E, t) &= F^{\text{meas}}(\lambda_*, E, t) C_F(E, t) lc(\lambda_*, E, t) \nonumber \\ &= F(\lambda_*, E, t) C_{\text{ref}}^{\text{norm}}(\bar{\lambda}_*, E) \tag{D.6} \\ \text{where } C_{\text{ref}}^{\text{norm}}(\bar{\lambda}_*, E) &= \frac{TF_{\text{ref}}}{\sum_{\lambda_*} F_{\text{ref}}(\lambda_*, E) d\lambda_*} \nonumber \end{align} $$

### 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$, as the weighted average (Sect. 7.3) of a selection of out-of-transit $F^{\text{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^{\text{[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$ 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{aligned}
F^{\text{intr}}(\lambda_{\star}, E, t) &= \frac{F^{\text{diff}}(\lambda_{\star}, E, t)}{1 - lc(\lambda_{\star}, E, t)} \\
&= C_{\text{ref}}^{\text{norm}}(\bar{\lambda}_{\star}, E) \frac{\sum_{i \in I_{\text{occ}}} f_i(\lambda_{\star}, E) S_i^{\text{[thick]}}(\bar{\lambda}_{\star}, E)}{1 - lc(\lambda_{\star}, E, t)}.
\end{aligned}
\quad (\text{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^{\text{intr}}(\lambda_{\star}, E, t) \sim C_{\text{ref}}^{\text{norm}}(\bar{\lambda}_{\star}, E) \frac{f_{\text{occ}}(\lambda_{\star}, E, t) S_{\text{occ}}^{\text{[thick]}}(\bar{\lambda}_{\star}, E, t)}{1 - lc(\lambda_{\star}, E, t)}. \quad (\text{D.8}) $$

The surface density flux from the region occulted during the exposure at time $t$ writes as:

$$ f_{\text{occ}}(\lambda_{\star}, E, t) = I_{\text{occ}}(\lambda_{\star}, E, t) \alpha_{\text{occ}}(\bar{\lambda}_{\star}, E, t), \quad (\text{D.9}) $$

where $I$ is the specific intensity emitted by the occulted region and $\alpha$ its broadband intensity variations (such as limb- and

gravity-darkening). During transit, the spectral flux scaling $lc$ over a given band at $\bar{\lambda}_\star$ corresponds to (Sect. 6.3):

$$ 1 - lc(\bar{\lambda}_{\star}, E, t) = \frac{\bar{I}(v) \alpha_{\text{occ}}(\bar{\lambda}_{\star}, E, t) \bar{S}_{\text{occ}}(E, t)}{\sum_{\lambda_{\star}} F_{\star}(\lambda_{\star}, E) d\lambda_{\star}}, \quad (\text{D.10}) $$

where $\bar{S}_{\text{occ}}$ and $\bar{I}$ are the achromatic planetary surface and stellar intensity values that allow reproducing the integrated flux ratio $\sum_{\lambda_{\star}} F(\lambda_{\star}, E, t) d\lambda_{\star} / \sum_{\lambda_{\star}} F_{\star}(\lambda_{\star}, E) d\lambda_{\star}$ measured in this band, as noted in Sect. D.2. While $\bar{S}_{\text{occ}}$ does not necessarily match the true $S_{\text{occ}}^{\text{[thick]}}$, because of the spectral variations of $I_{\text{occ}}(\lambda_{\star}, E, t)$ over the band, we assume that the ratio $\bar{S}_{\text{occ}}/S_{\text{occ}}^{\text{[thick]}}$ remains constant during transit, so that:

$$ F^{\text{intr}}(\lambda_{\star}, E, t) = C_{\text{ref}}^{\text{norm}}(\bar{\lambda}_{\star}, E) \sum_{\lambda_{\star}} F_{\star}(\lambda_{\star}, E) d\lambda_{\star} \times \frac{I_{\text{occ}}(\lambda_{\star}, E, t) \alpha_{\text{occ}}(\bar{\lambda}_{\star}, E, t) S_{\text{occ}}^{\text{[thick]}}(\bar{\lambda}_{\star}, E, t)}{\bar{I}(v) \alpha_{\text{occ}}(\bar{\lambda}_{\star}, E, t) \bar{S}_{\text{occ}}(E, t)}, \quad (\text{D.11}) $$

Can be rewritten as:

$$ F^{\text{intr}}(\lambda_{\star}, E, t) = F_{\text{ref}}^{\text{norm}}(\bar{\lambda}_{\star}, E) I_{\text{occ}}(\lambda_{\star}, E, t). \quad (\text{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 \text{ m s}^{-1}$ 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