in the wavelength direction. An alpha-shape algorithm is then applied over the stretched spectrum to remove maxima formed by blended lines, setting the radius of the rolling pin to a larger size than the typical spectral line width. Finally, outlying maxima are removed and the remainder used to define the continuum of the master spectrum via linear interpolation. Smoothing and local maxima windows, stretching parameter, and rolling pin radius can be adjusted depending on the instrumental resolution and stellar line properties to better control the continuum determination.

## E.2. Spectral weight profiles

### E.2.1. True error estimate

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

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

### E.2.2. Definition of weight profiles

At their latest processing stage disk-integrated spectra write as (see Eq. 9):

$$F^{\text{sc}}(\lambda_{\star}, E, t) = F^{\text{meas}}(\lambda_{\star}, E, t)C_{\text{corr}}(\lambda_{\star}, E, t)lc(\lambda_{\star}, E, t), \quad (\text{E.1})$$

where $C_{\text{corr}}$ groups all spectral corrections that were applied to the measured spectra (Sect. 5), except for cosmics and wiggles. For the sake of clarity we only report the dependency on $\lambda_{\star}$ in the following expressions. Uncertainties on the measured flux density can be expressed as a function of the raw photoelectron count (Eq. 3), yielding the following uncertainty on the true flux $F^{\text{sc, true}}$ measured during $\Delta t$:

$$\sigma^{\text{sc, true}}(\lambda_{\star}) = \tilde{g}_{\text{cal}}(\lambda_{\star}) \sqrt{N^{\text{true}}(\lambda_{\star})} C_{\text{corr}}(\lambda_{\star}) lc(\lambda_{\star}) / \Delta t. \quad (\text{E.2})$$

The spectra $F^{\text{sc}}/lc$ have the same low-resolution profile as the spectrum of the unocculted star $F_{\star}$, allowing us to estimate the raw number of photoelectrons as:

$$N^{\text{true}}(\lambda_{\star}) \sim \frac{F_{\star}^{\text{meas}}(\lambda_{\star}) \Delta t}{\tilde{g}_{\text{cal}}(\lambda_{\star}) C_{\text{corr}}(\lambda_{\star})}, \quad (\text{E.3})$$

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

Then, based on Eq. 11 and 12, the error on the true intrinsic flux is:

$$\sigma^{\text{intr, true}}(\lambda) = \frac{\sigma^{\text{diff, true}}(\lambda)}{1 - lc(\lambda)} \quad (\text{E.4})$$
$$\text{where } \sigma^{\text{diff, true}}(\lambda) = \sqrt{\sigma_{\star}^{\text{true}}(\lambda)^2 + \sigma^{\text{sc, true}}(\lambda)^2}.$$

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

## E.3. Temporal and spatial resampling

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

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

## E.4. Conversion from 2D to 1D spectra

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

### E.5. Cross-correlation

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

In a first step we identify lines of the binary mask that can contribute to the cross-correlation in all processed exposures. Only lines for which the requested rv range (once converted into a wavelength range around this line) is fully within the spectral range of the processed spectrum (which is not the case, for example, if the line is at the edge of a spectral order) and contains no undefined pixels are kept. The spectrum of each exposure is then cross-correlated with the filtered mask as follows:

$$CCF(sv, E, t) = \sum_l F(\lambda_l, E, t) \frac{w_l}{<g_l>_{cal}}. \quad (E.5)$$

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

## Appendix F: Analysis methods

### F.1. Stellar mask generator

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

The method starts by normalizing the master spectrum with the estimated stellar continuum (Sect. 7.1), smoothing it to remove high-frequency noise, and oversampling it on a regular spectral grid. Local maxima and minima are identified as

bins whose flux value is the highest (resp. lowest) in their closest neighbourhood, defined as a fraction of the lines FWHM. Lines can then be defined as the unique successions of distinct maximum-minimum-maximum, with the line central wavelength set to the minimum position. Lines that fall within spectral ranges known to be contaminated by deep telluric lines, or ranges selected by the user, are excluded. Successive selections can then be applied to refine the linelist, with thresholds that can be adjusted depending on the quality of the data and the type of the star (the process is illustrated in Fig. 18):

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

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

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

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

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

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

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

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

### F.2. Analysis of individual stellar profiles

#### F.2.1. Model line profiles

All stellar profiles are defined with a polynomial continuum, multiplied by spectral absorption lines calculated analytically or numerically. Analytical profiles provide a direct way to fit a single stellar line with simple shape in disk-integrated or intrinsic

<footnote n="11">¹¹ Since a precise systemic rv is later derived from disk-integrated CCFs, we only use an approximate value at this stage and keep in mind that masks are defined with respect to this reference.</footnote>
<figure><img src="image_39.png" /><figcaption>**Fig. F.1.** Flowchart of CCF mask generation, illustrated with HD 209458. Thresholds applied in each step (magenta lines) can be adjusted based on the resulting line selection displayed as 2D distributions (green and red disks show lines kept and excluded, respectively), occurrence histograms, and cumulative functions of the line weights. We illustrate how telluric contamination is flagged by showing the master telluric spectrum at its minimum (blue) and maximum (red) doppler shifts during the two epochs, relative to the master stellar spectrum (black profile). The stellar line at ~7223.7 Å is excluded due to contamination by the telluric line highlighted at its mean position by the dashed grey line.</figcaption></figure>

spectra. Numerical profiles are calculated using ANTARESS stellar surface model (Sect. 3.3), summing elementary profiles over the full star to fit disk-integrated profiles or over planet-occulted regions to fit intrinsic profiles. As described in Sect. 3.3.1 the elementary profiles can be analytical, theoretical, or measured. The last two options allow fitting several lines over a portion of the stellar spectrum (Fig. 23).

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

### F.2.2. Line profile fitting

Disk-integrated profiles from original time series are analyzed in the rest frame of the input data, where their derived centroids can be used to check the Keplerian model and measure the systemic rv of the system. After being aligned, profiles can be analyzed in the star rest frame. We highlight that masking narrow planetary atmospheric features does not remove the effect of planetary occultation on disk-integrated profiles. Indeed, the opaque planet continuum generates spectral features through its absorption of the local stellar line profile (e.g., Dethier & Bourrier 2023). These features yield the so-called “Doppler shadow” (Collier Cameron et al. 2010) that distorts in-transit line profiles and induces the RM anomaly when deriving their centroid with an undistorted model of the unocculted star. We propose an original method to circumvent this bias, which consists in excluding the range absorbed by the planet-occulted stellar line. This range can be defined in a given exposure as a rv window covering the
<figure><img src="image_40.png" /><figcaption>**Fig. F.2.** Correlation diagrams for the PDFs of the RMR model parameters of the WASP-76b transits. Ci and FWHMi indicate polynomial coefficients for the line shape properties (see text). Green and blue lines show the 1 and 2 simultaneous 2D confidence regions that contain, respectively, 39.3% and 86.5% of the accepted steps. 1D histograms correspond to the distributions projected on the space of each line parameter, with the green dashed lines limiting the 68.3% HDIs. The blue lines and squares show the median values.</figcaption></figure>

breadth of the line, centered around the line transition Doppler-shifted by the surface *rv* of the planet-occulted region (Eq. 2). Masked disk-integrated stellar line profiles are not distorted anymore by the planetary continuum, and can be fitted with the unocculted stellar line model.

Intrinsic profiles from original time series are analyzed in the star rest frame, where their derived centroids can be used to determine the best model for stellar surface *rv* along the transit chord. After being aligned, profiles can be analyzed in the photosphere rest frame. We highlight the importance of using a Bayesian approach when fitting intrinsic profiles (Bourrier et al.

2021), to better explore the line model properties and their correlations, and to retrieve information from noisy or masked local stellar lines that could not be obtained through least-square minimization. Priors set on the model properties can be informed by their range of variations along the transit chord (Sect. 8.2.3) or physical arguments (e.g., local FWHM smaller than the disk-integrated value; contrast between 0 and 1; RV centroid smaller than the projected stellar rotational velocity). PDFs from the fits to intrinsic stellar lines inform on the exposures that can be exploited in time series analysis (Sect. 8.3). Exposures with flat, non-constraining PDFs for the width, contrast, and *rv* of the intrinsic line, typically at the limb of the star where the planet
only partially occults dim regions of the stellar surface, should be excluded (Fig. 17).

### F.2.3. Line profile variations

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

## F.3. Analysis of joint intrinsic stellar profiles

### F.3.1. Global model line profiles

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

### F.3.2. RMR fitting

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

The coefficients describing the line profiles depend on the type of profile (e.g., FWHM and contrast for analytical profiles, abundances for theoretical profiles, etc), but their spatial

dependence along the transit chord is always described using the polynomials described in Sect. F.2.3. We remind that analytical models for intrinsic lines can include macroturbulence (Sect. 3.3.1), which can be used to separate its contribution from that of thermal broadening.

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

# Appendix G: Performance analysis

## G.1. Resampling

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

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

We fitted Gaussian profiles to deep and narrow stellar lines in the red, high S/N parts of the spectrum to assess the impact of resampling on the line shape. Although the lines were systematically shallower and broader when resampled, the difference remained marginal ($1\sigma$ for the constrast, $0.5\sigma$ for the FWHM, $0.05\sigma$ for the *rv*). We derived the same conclusion from fits to the cross-correlated stellar line averaged over out-of-transit exposures. We thus conclude that resampling spectra after they are aligned in the star rest frame has a negligible impact on individual stellar lines. This is because the Keplerian motion of HD 209458 around the transit windows only represents a small fraction of an ESPRESSO pixel, so that the interpolation due to resampling remains limited.
<figure><img src="image_41.png" /><figcaption>**Fig. G.1.** Average out-of-transit spectrum of HD 209458 (epoch 1), when resampling spectra on a common wavelength grid (red) or not (blue). In the former case, undercorrections of telluric lines (highlighted by the best-fit telluric model at mid-epoch, in green) in individual exposures add up and yield biases visible by eye.</figcaption></figure>

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

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

<figure><img src="image_42.png" /><figcaption>**Fig. G.2.** Master intrinsic CCF of HD 209458, averaged along the transit chord in epoch 1. The blue, red, and green profiles were calculated with rv steps of 0.50, 0.25, and 0.10 km s⁻¹, respectively.</figcaption></figure>

## G.2. Covariance matrix

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

## G.3. Weighted means

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

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

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

Figure G.5 illustrates the role of disk-integrated stellar lines when averaging intrinsic spectra in the photosphere rest frame (Sect. E.2.2). In that case, disk-integrated stellar lines shift between different exposures, so that a given pixel should be weighted differently according to the depth of the stellar line it overlaps with. Similarly to telluric lines, the strongest biases are observed for the deepest disk-integrated stellar lines, when they are not included in the weight profile. The sodium lines are among the deepest in the disk-integrated spectrum of HD 209458, with nearly full absorption at the core of the line,
<figure><img src="image_43.png" /></figure>

**Fig. G.3.** Weighing performance (instrumental flux calibration) for disk-integrated spectra. **Top panel:** Overlapping orders in 2D flux spectra (exposure at index 10 in HD 209458 epoch 1). **Middle panel:** Zoom on the 1D spectrum averaged from the order spectra, when accounting (black) or not (magenta) for instrumental calibration in the weight profiles. **Bottom panel:** Relative difference between the magenta and black profiles, overlapped with the contribution of the instrumental calibration to the order weight profiles (arbitrary scale).

and induce biases up to tens of percents even in the intrinsic profile averaged over all in-transit exposures. It is thus critical for the analysis of averaged planet-occulted stellar lines to account for disk-integrated stellar lines in weight profiles.

### G.4. Flux balance corrections

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

Here we evaluate how much a poor correction of the flux balance impacts the quality of disk-integrated and intrinsic stellar CCFs. We reprocessed the data from the FLUX BALANCE method (Sect. 5.3) onward, using low-order polynomials instead of smoothing splines to fit the color balance variations. This approach allows capturing the overall shape of the color balance, but not its variations over a few spectral orders as finely as with

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

**Fig. G.4.** Weighing performance (tellurics) for disk-integrated spectra. **First (top) panel:** Typical telluric spectrum in HD 209458 epoch 1. **Second panel:** Relative difference between the master-out spectra accounting or not for telluric correction in weight profiles. **Third panel:** Same as the second panel, when averaging only two exposures (index 5 and 84). **Fourth (bottom) panel:** Zoom from the third panel in the region of the sodium doublet, illustrating how shifting telluric lines (reported on an arbitrary scale, in blue and red at indexes 5 and 84) bias the mean over the pixels they absorb when not included in weight profiles.

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

**Fig. G.5.** Weighing performance (stellar lines) for intrinsic spectra. **Top panel:** Mean of intrinsic spectra along the transit chord of HD 209458b in the region of the sodium doublet, when accounting (black) or not (magenta) for the disk-integrated spectrum in weight profiles. **Bottom panel:** Relative difference between the magenta and black profiles, overlapped with the contribution of the disk-integrated spectrum to the weight profiles (arbitrary scale, in blue and red for two exposures at ingress and egress).

splines (Fig. 7). We first assessed the quality of disk-integrated CCFs by comparing the dispersion of their out-of-transit properties with the nominal analysis. The degraded flux balance correction does not have a clear impact on WASP-76 CCFs. For HD 209458, it increases the dispersion of the CCF contrast and FWHM by ~3–8% but has no impact on their RVs.
<figure><img src="image_46.png" /><figcaption>Fig. G.6. Wiggle correction in epoch 1 of HD 209458b, for exposures at index 60 (top panel) and 88 (bottom panel). Blue and red profiles show spectra before and after correction, respectively. Flux spectra are normalized by the stellar continuum (estimated from the corrected data; Sect. 7.1) and offset for clarity. The dashed black line represents the normalized continuum unity. Transmission spectra, shown at the top of the panels and calculated with the master out-of-transit spectrum, highlight wiggle patterns strong enough to be visible in the flux spectra.</figcaption></figure>

This is likely because the residual, uncorrected color variations are broad enough that they are approximately constant over the breadth of a stellar line and thus do not modify their shape. Since the correction improves the disk-integrated CCFs of HD 209458 we further evaluated its impact on the intrinsic CCFs and their RMR fit. The degraded correction yields $\lambda = 1.12 \pm 0.07^{\circ}$ and $v_{\text{eq}} \sin i_{\star} = 4.289 \pm 0.008 \text{ km s}^{-1}$, consistent with $1.06 \pm 0.07^{\circ}$, $4.272 \pm 0.007 \text{ km s}^{-1}$ from the nominal analysis.

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

## G.5. Wiggle correction

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

Here, we assess its impact on stellar CCFs and their interpretation. Similarly to color balance variations, but at shorter frequencies, wiggles indeed change the relative flux balance between different spectral regions over time and thus the relative weights between CCF lines. We deactivated the wiggle correction (Sect. 5.6) and reprocessed the data from the SPECTRAL LINE DETRENDING method (Sect. 6.1) onward. Correcting for the wiggles has no significant impact on WASP-76 disk-integrated CCFs. For HD 209458, it improves the dispersion of their out-of-transit contrast and FWHM by ~3–6% but has

no impact on their rvs. This is likely because the wiggle pattern is too broad to affect the line shapes and thus only changes their overall flux level, as discussed for the flux balance correction (Sect. G.4). We then evaluated the impact of the wiggles on the intrinsic CCFs of HD 209458 and their RMR fit. Disabling the correction yields $\lambda = 1.01 \pm 0.08^{\circ}$ and $v_{\text{eq}} \sin i_{\star} = 4.271 \pm 0.007 \text{ km s}^{-1}$, to be compared to $1.06 \pm 0.07^{\circ}$, $4.272 \pm 0.007 \text{ km s}^{-1}$ for the nominal analysis.

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

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

The ANTARESS workflow offers the new possibility to apply the RMR technique (Sect. 8.3.1) to “white” CCFs calculated by cross-correlating intrinsic spectra over broad spectral ranges:

$$ \text{CCF}[F^{\text{intr}}](rv) = \sum_{l} \frac{F_{\star}^{l}(\lambda_{l}) - F^{\text{sc}}(\lambda_{l})}{1 - l c(\lambda_{l})} \alpha_{l}, \qquad (\text{G.1}) $$

where $\alpha_l$ designates the combined weights applied to local line spectra (Sect. 7.5). It is also possible, as in past studies (e.g., Bourrier et al. 2021), to use the CCFs of disk-integrated spectra, readily provided by the DRS of high-resolution spectrographs or calculated with the ANTARESS workflow (Sect. 6.5), and extract intrinsic CCFs as:

$$ \text{CCF}^{\text{intr}}(rv) = \frac{\text{CCF}_{\star}(rv) - \text{CCF}^{\text{sc}}(rv)}{1 - l c_{\text{wh}}}, \qquad (\text{G.2}) $$

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

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

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

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

### G.7. Chromatic flux scaling

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

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

### G.8. Stellar CCF masks

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

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