# ATMOSPHERIX: I- an open source high-resolution transmission spectroscopy pipeline for exoplanets atmospheres with SPIRou

Baptiste Klein<sup>1,2</sup> Florian Debras<sup>1,★</sup> Jean-François Donati<sup>1</sup> Thea Hood<sup>1</sup> Claire Moutou,<sup>1</sup>

Andres Carmona,<sup>3</sup> Merwan Ould-elkhim,<sup>1</sup> Bruno Bézard,<sup>4</sup> Benjamin Charnay,<sup>4</sup> Pascal Fouqué,<sup>1</sup>

Adrien Masson,<sup>4</sup> Sandrine Vinatier,<sup>4</sup> Clément Baruteau<sup>1</sup>,<sup>1</sup> Isabelle Boisse,<sup>5</sup> Xavier Bonfils,<sup>3</sup>

Andrea Chiavassa,<sup>6</sup> Xavier Delfosse,<sup>3</sup> William Dethier,<sup>3</sup> Guillaume Hebrard,<sup>7,8</sup> Flavien Kiefer,<sup>7</sup>

Jérémy Leconte,<sup>9</sup> Eder Martioli,<sup>7,10</sup> Vivien Parmentier,<sup>2,6</sup> Pascal Petit,<sup>1</sup> William Pluriel,<sup>11</sup> Franck Selsis,<sup>9</sup>

Lucas Teinturier,<sup>4,12</sup> Pascal Tremblin<sup>1</sup>,<sup>13</sup> Martin Turbet,<sup>9,12</sup> Olivia Venot<sup>14</sup> and Aurélien Wyttenbach<sup>3</sup>

Accepted 2023 August 24. Received 2023 August 22; in original form 2023 March 17

## ABSTRACT

Atmospheric characterization of exoplanets from the ground is an actively growing field of research. In this context, we have created the ATMOSPHERIX consortium: a research project aimed at characterizing exoplanets atmospheres using ground-based high-resolution spectroscopy. This paper presents the publicly available data analysis pipeline and demonstrates the robustness of the recovered planetary parameters from synthetic data. Simulating planetary transits using synthetic transmission spectra of a hot Jupiter that were injected into real SPIRou observations of the non-transiting system Gl 15 A, we show that our pipeline is successful at recovering the planetary signal and input atmospheric parameters. We also introduce a deep learning algorithm to optimize data reduction which proves to be a reliable, alternative tool to the commonly used principal component analysis. We estimate the level of uncertainties and possible biases when retrieving parameters such as temperature and composition and hence the level of confidence in the case of retrieval from real data. Finally, we apply our pipeline onto two real transits of HD 189733 b observed with SPIRou and obtain similar results than in the literature. In summary, we have developed a publicly available and robust pipeline for the forthcoming studies of the targets to be observed in the framework of the ATMOSPHERIX consortium, which can easily be adapted to other high resolution instruments than SPIRou (e.g. VLT-CRIRES, MAROON-X, ELT-ANDES).

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

## 1 INTRODUCTION

More than 5000 exoplanets were discovered in the last decade, paving the way for statistical exploration of the orbital and physical properties of planetary systems (e.g. Udry & Santos 2007; Fulton

et al. 2017; Debras, Baruteau & Donati 2021). One of the main objective for the next decade is now to understand thoroughly the physical nature of individual planets. This necessarily requires an in-depth study of their atmosphere in order to lift degeneracies between seemingly identical planets in terms of mass and radius (e.g. Valencia et al. 2013). JWST and Ariel (Tinetti et al. 2021) space missions will play a key role in that venture by providing high quality observations for a large number of exoplanets. However, space-based observations

© 2023 The Author(s).

<footnote n="1"><sup>1</sup>IRAP, Université de Toulouse, CNRS, UPS, F-31400 Toulouse, France</footnote>

<footnote n="2"><sup>2</sup>Department of Physics, University of Oxford, Oxford OX13RH, UK</footnote>

<footnote n="3"><sup>3</sup>Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="4"><sup>4</sup>LESIA, Observatoire de Paris, Université PSL, Sorbonne Université, Université Paris Cité, CNRS, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="5"><sup>5</sup>Aix Marseille Universite, CNRS, Laboratoire d'Astrophysique de Marseille UMR 7326, F-13388 Marseille, France</footnote>

<footnote n="6"><sup>6</sup>Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Lagrange, CS F-34229 Nice, France</footnote>

<footnote n="7"><sup>7</sup>Institut d'astrophysique de Paris, UMR7095 CNRS, Université Pierre & Marie Curie, 98bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="8"><sup>8</sup>Observatoire de Haute-Provence, CNRS, Université d'Aix-Marseille, Saint-Michel-l'Observatoire, F-04870, France</footnote>

<footnote n="9"><sup>9</sup>Laboratoire d'astrophysique de Bordeaux, Université de Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, F-33615 Pessac, France</footnote>

<footnote n="10"><sup>10</sup>Laboratório Nacional de Astrofísica, Rua Estados Unidos 154, 37504-364, Itajubá - MG, Brazil</footnote>

<footnote n="11"><sup>11</sup>Observatoire astronomique de l'Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="12"><sup>12</sup>Laboratoire de Météorologie Dynamique/IPSL, CNRS, Sorbonne Université, Ecole Normale Supérieure, Université PSL, Ecole Polytechnique, Institut Polytechnique de Paris, F-75005 Paris, France</footnote>

<footnote n="13"><sup>13</sup>Université Paris-Saclay, UVSQ, CNRS, CEA, Maison de la Simulation, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="14"><sup>14</sup>Université de Paris Cité and Univ Paris Est Creteil, CNRS, LISA, F-75013 Paris, France</footnote>

<footnote>*E-mail: florian.debras@irap.omp.eu</footnote>
of planet atmospheres have limits which are best overcome from the ground using high-resolution spectroscopy (HRS) with numerous large telescopes. Through cross-correlation HRS, relying on the statistical comparison between an absorption or emission spectrum of the planet atmosphere and theoretical models, one can extract the planetary signal which is typically 10–100 times weaker than the noise. Since the first successful characterization of an exoplanet atmosphere with HRS, a decade ago by Snellen et al. (2010), this technique has been substantially refined (e.g. Brogi et al. 2012; Birkby et al. 2013; de Kok et al. 2013; Brogi et al. 2016; Alonso-Floriano et al. 2019; Brogi & Line 2019; Giacobbe et al. 2021; Guilluy et al. 2022), and has acquired the necessary maturity to become a reliable complement to forthcoming space-based missions (B Brogi et al. 2017; Brogi & Line 2019; Kasper et al. 2023).

SPIRou (Donati et al. 2020), a high-resolution near-infrared (nIR) spectropolarimeter mounted at the Canada–France–Hawaii Telescope, is perfectly suited for this task. Thanks to its broad continuous wavelength coverage of the nIR, from 0.9 to 2.5 µm, SPIRou has the ability to resolve a high number of molecular lines in the emission or transmission spectra of planetary atmospheres. More specifically, the observations are divided into 50 overlapping diffraction orders spanning the *Y*, *J*, *H*, and *K* bands at a resolving power of ~70 000 (~2.28 km s⁻<fnref n="1" /> velocity bin). SPIRou has already been successfully used to detect water and carbon monoxide by performing emission spectroscopy of τ Boo b (Pelletier et al. 2021) and transmission spectroscopy during two transits of HD 189 733 b (Boucher et al. 2021) as well as to detect Helium on several targets (Allart et al. 2023, Masson et al. in preparation). In most cases, the absorption lines of the planet atmosphere were found to be Doppler shifted compared to theoretical predictions which remains to be understood in the light of atmospheric circulation (see e.g. Flowers et al. 2019).

With the aim of optimizing the capabilities of SPIRou for the characterization of the atmosphere of exoplanets, we have gathered a large, French-led community of observers and theoreticians, specializing in exoplanet atmospheres and stellar observations and simulations, under the ATMOSPHERIX programme. This programme aims at observing a wide range of exoplanets over several years in order to (i) constrain the composition of their atmosphere, (ii) probe the pressure–temperature (PT) profile and the amplitude of atmospheric winds through Doppler spectroscopy, and (iii) characterize the extended atmosphere through the He I metastable triplet at 1083 nm. Additionally, long-term repeated observations of a sample of planets will allow us to better understand variability in exoplanet atmospheres (e.g. Armstrong et al. 2016; Komacek & Showman 2020; Cho, Skinner & Thrastarson 2021).

This paper introduces a series of studies of the atmosphere of transiting planets observed with nIR high-resolution spectrographs as part of the ATMOSPHERIX programme. In this study, we present our publicly available code.<fnref n="1" /> to extract a planet transmission spectrum from a time-series of nIR high-resolution spectra. The extraction pipeline is applied to (i) sequences of synthetic transmission spectra of a hot Jupiter that mimics the properties of HD 189733 b and that was injected into SPIRou observations of the bright quiet M dwarf Gl 15 A and (ii) on the two same transits of HD 189 733 b as Boucher et al. (2021). The Gl 15 A input data sets are described in Section 2. Section 3 details the algorithm to extract the planet atmosphere signal from the observed sequence of spectra and infer the planet atmosphere parameters in a statistically robust way. We then present

our retrieval methods on synthetic data in Section 4 and their application on real data in Section 5. We discuss our results and their implications for real targets in Section 6 and conclude in Section 7. A companion paper (Debras et al., submitted to MNRAS) studies the biases and degeneracies in the planet atmosphere parameters retrieved with the pipeline presented here.

## 2 OBSERVATIONS AND PLANET INJECTION

We simulate nIR spectroscopic observations of a planetary transit by injecting a synthetic planet atmosphere spectrum into a sequence of spectra of the bright M dwarf Gl 15 A, collected with SPIRou in October 2020 (see Table 1). We first describe the stellar data before detailing how we injected the planet.

### 2.1 Input stellar spectra

Gl 15 A has been intensively monitored with SPIRou over the last 3 yr and does not have any known transiting planet, making it an ideal target to benchmark our data analysis code. Our input observations consist in a series of 192 consecutive spectra collected with SPIRou on 2020 October 8, spanning a total of 5 h. We divided these 192 spectra into two series of 96 spectra (the odd and even file numbers, respectively) in order to ensure the robustness of our pipeline on two sets of data. The peak signal-to-noise ratio (SNR) per 2.28 km s⁻<fnref n="1" /> velocity bin ranges from 250 to 320 (median of 291), and the airmass from 1.1 and 1.3 (see panels 2 and 4 of Fig. 1). Note that the binary companion, Gl 15 B, is a M 3.5 dwarf located at 146 au from Gl 15 A (Reid, Hawley & Gizis 1995). The velocimetric effect of this binary system is neglected in the present analysis given the low acceleration that B induces on A's RV (about 2 m s⁻<fnref n="1" /> per yr; Howard et al. 2014). We also neglect the RV effect of the two recently detected planets around Gl 15 A (Howard et al. 2014; Pinamonti et al. 2018), inducing respective signatures of 1.68 and 2.5 m s⁻<fnref n="1" /> modulated with orbital periods of 11.44 and 7600 d, respectively.

For this paper, the SPIRou observations were reduced using the version 0.6.132 of APERO, the official data reduction software (DRS) of the instrument Cook et al. (2022). In short, the pipeline applies the optimal extraction method of Horne (1986) to extract each individual exposure from the H4RG detector (Artigau et al. 2018). The wavelength solution is obtained by combining calibration exposures of a UNe hollow-cathod lamp and a thermally stabilized Fabry-Periot Etalon (Bauer, Zechmeister & Reiners 2015; Hobson et al. 2021). APERO performs a correction of the telluric contamination using a method, summarized in Cook et al. (2022, see section 8), which will be presented in a forthcoming paper (Artigau et al. in preparation). This technique applies TAPAS (Bertaux et al. 2014) to pre-clean telluric absorption and the low level residuals are removed in using a data set of spectrum of hot stars observed in different atmospheric conditions to build a residual models in function of few parameters (optical depths of H₂O and of dry components). Note that the deepest telluric lines (relative absorption larger than 90 per cent) are masked out by the pipeline as the low amount of transmitted flux will most likely result in an inaccurate telluric correction. Our input sequence of spectra contains the blaze- and telluric-corrected spectra.

### 2.2 Planet injection

We then inject synthetic planet atmosphere transmission spectra on top of the APERO-provided telluric-corrected spectra of Gl 15 A. We consider the case of a hot Jupiter (based on HD 189 733 b, Bouchy et al. 2005). The injected planet spectra are generated using

<footnote n="1">¹https://github.com/baptklein/ATMOSPHERIX_DATA RED</footnote>
<table><caption>**Table 1.** Physical parameters for Gl 15 A, HD189733 b, and for the simulated hot Jupiter used in the study. When taken from the literature, the reference for each parameter is indicated in the right-hand column†. Note that references cited for planet parameters refer to HD189733 b.</caption><thead><tr><th>Stellar parameters</th><th>Value</th><th>Gl 15 A</th><th>Reference</th></tr></thead><tbody><tr><td>Mass (M<sub>☉</sub>)</td><td>0.398 ± 0.004</td><td></td><td>Cr22</td></tr><tr><td>Radius (R<sub>☉</sub>)</td><td>0.388 ± 0.013</td><td></td><td>Cr22</td></tr><tr><td>Effective temperature (K)</td><td>3603 ± 60</td><td></td><td>Cr22</td></tr><tr><td>H magnitude</td><td>4.476 ± 0.2</td><td></td><td>Cu03</td></tr><tr><td>Systemic velocity [km s<sup>-1</sup>]</td><td>11.73 ± 0.1</td><td></td><td>Fo18</td></tr><tr><td>Limb darkening (Quadratic)</td><td>0.0156, 0.313</td><td></td><td>Cl11</td></tr><tr><th colspan="4">Planet parameters</th></tr><tr><td></td><td>HD 189733 b</td><td>Synthetic planet</td><td>Reference</td></tr><tr><td>Transit depth (%)</td><td>2.2 ± 0.1</td><td>2.2</td><td>Ad19</td></tr><tr><td>Radius (R<sub>J</sub>)</td><td>1.142 ± 0.040</td><td>0.57</td><td>Ad19</td></tr><tr><td>Mass (M<sub>J</sub>)</td><td>1.13 ± 0.05</td><td>0.568</td><td>Ad19</td></tr><tr><td>g (m s<sup>-2</sup>)</td><td>22.45 ± 1.5</td><td>45.29</td><td>From planet mass and radius</td></tr><tr><td>Orbital period (d)</td><td>2.218579 ± 0.000001</td><td>2.218577</td><td>Ad19</td></tr><tr><td>Mid-transit time (BJD TBD)</td><td>2458334.990899 ± 0.0007</td><td>2459130.8962180</td><td>Ad19</td></tr><tr><td>Inclination (deg)</td><td>85.3 ± 0.2</td><td>90.0</td><td>Ad19</td></tr><tr><td>Eccentricity</td><td>~0.0</td><td>0.0</td><td>–</td></tr><tr><td>Equilibrium temperature (K)</td><td>1203 ± 39</td><td>900</td><td>Ad19 &amp; Ro21</td></tr><tr><td>Orbital semi-amplitude (km s<sup>-1</sup>)</td><td>151.2 ± 4.5</td><td>120.0</td><td>–</td></tr><tr><td>Transit duration (h)</td><td>1.84 ± 0.04</td><td>1.84</td><td>Ad19</td></tr></tbody></table>

Note. † To gain some space in the table, we use aliases for the references. Cr22, Cu03, Fo18, Cl11, Ad19, and Ro21 stand respectively for Cristofari et al. (2022), Cutri et al. (2003), Fouqué et al. (2018), Claret & Bloemen (2011), Addison et al. (2019), and Rosenthal et al. (2021).

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Continuum-normalized transit light curve of HD 189733 b (top panel), airmass (panel 2), topocentric-to-stellar rest frame RV correction (RV<sub>corr</sub>, panel 3), and peak SNR per velocity bin during the two simulated transits of the HD 189733 b analogue (panel 4). Note that RV<sub>corr</sub> contains the RV contributions of the barycentric Earth motion and of the systemic velocity of the star. On panels 1 and 4, the two different transits are respectively shown as blue dots and pink crosses. The vertical grey band (resp. vertical grey dotted line) indicates the mid-transit primary transit (resp. mid-transit time) of the simulated planet. The horizontal grey dotted line on the bottom panel indicates the average value of the peak SNR for the observed spectra.</figcaption></figure>

petitRADTRANS (Mollière et al. 2019), which gives the planet radius as a function of wavelength assuming an isothermal planet atmosphere solely containing chosen molecules at a constant volume mixing ratio. This radius is then transformed into an absorbed flux by multiplying the total flux by the ratio of planetary to stellar radius squared, called transit depth. This model is then convoluted with a Gaussian of half-width 2 SPIRou pixels (i.e. ~4.5 km s<sup>-1</sup> for a resolving power of ~70 000 in the nIR Donati et al. 2020) to account for the instrumental broadening.

Since Gl 15 A is significantly smaller than HD 189733, some adjustments of the injected planet parameters are needed to keep the simulations realistic. We decided to conserve 4 quantities: (i) the transit depth, (ii) the transit duration, (iii) the ratio between the stellar radius and the atmospheric scale height, and (iv) a consistent atmospheric temperature at the limbs. Our synthetic planet therefore has a lower mass, radius, and velocimetric amplitude than HD189733 b, but a larger surface gravity. Although not a physical planet (the orbital mass is not consistent with the gravitational mass), the injected planet represents a good observational analogue of a hot Jupiter. The stellar properties were left untouched in our simulated data. The parameters adopted for synthetic planet are given in Table 1.

The planet orbit is assumed circular with a mid-transit time corresponding to the mean values of our observation times (see the two transits in Fig. 1). For each spectrum, we generate a transit curve $F_C$ using the python module `batman` (Kreidberg 2015), assuming an aligned circular planet orbit and using the H-band quadratic limb- darkening coefficients computed in Claret & Bloemen (2011) for Gl 15 A's properties (see Table 1). From the resulting transiting curves, we compute transit window functions $W_C$ whose values range from 0 (for out-of-transit times) and 1 (at mid-transit time), using $W_C = (1 - F_C)/\max(1 - F_C)$.

We then build a sequence of planet transmission spectra by applying the following steps to the synthetic planet atmosphere spectrum. First, the simulated planet atmosphere spectrum is multiplied by the transit window function $W_C(t)$ at time $t$. Secondly, the window-
weighted simulated spectrum is shifted in the Earth rest frame by correcting for the Barycentric Earth Radial Velocity $V_{BE}(t)$, the stellar systemic velocity $V_{sys}$ and the RV signature $V_{RV}(t)$ induced by the injected planet on the host star. We then shift this spectrum by an additional $30 \text{ km s}^{-1}$ corresponding to the planetary shift in velocity during transit plus three times the SPIRou resolution to ensure that the stellar and planetary molecular lines are separated. Note that, as Gl 15 A is a M-dwarf star, it contains water and carbon monoxide in its atmosphere which complicates the planetary retrieval. This is discussed further in the companion paper. The resulting synthetic spectrum is then convolved at SPIRou's spectral resolution, and multiplied by the stellar spectrum observed at time $t$. Our input sequence of spectra is finally built by repeating the steps listed earlier to all the observed spectra.

The final intensity $I_f$ as a function of time $t$ and wavelength $\lambda$ can thus be expressed as follows

$$I_f(t, \lambda) = I_i(t, \lambda) \left( 1 - W_C(t) \left( \frac{R_p(\lambda'(t))}{R_*} \right)^2 \right), \quad (1)$$

where $I_i$ is the initial intensity (i.e. the APERO reduced SPIRou observations), $R_*$ is the stellar radius<fnref n="2" />, and $R_p$ is the planetary radius (degraded at SPIRou resolution) which depends on wavelength because of the wavelength-dependent opacity of the planet atmosphere and $\lambda'$ the Doppler-shifted wavelength. As explained earlier

$$\lambda'(t) = \lambda \left( 1 - \frac{K_* \sin(2\pi\phi(t)) + V_*(t)}{c_0} \right), \quad (2)$$

$$V_*(t) = V_{BE}(t) - V_{sys} + 30\text{km.s}^{-1} \quad (3)$$

where $K_*$ is the orbital RV semi-amplitude of the star, $\phi$ the planet orbital phase centred on the mid-transit, $V_*(t)$ the non-orbital Doppler shift, and $c_0$ the speed of light.

Finally, we have also created synthetic sequences without the star, where the model is injected into a map of white noise with a variance equal to the SNR of the observations modulated by the blaze function that increases the variance at the edges of the orders. These synthetic models do not need to go through any further step of data analysis, and serve as references to identify the effects of the data analysis on the atmospheric retrieval.

# 3 DATA PROCESSING

## 3.1 Data cleaning

The extraction of the planetary signal requires to correct for the residual telluric absorption lines, the stellar spectrum and additional correlated noise induced by the instrument and the observing conditions. Following Boucher et al. (2021), we can perform an additional masking of the telluric lines with absorption deeper than 70 per cent of the continuum level.<fnref n="3" /> Diffraction orders 57–54 (i.e. ~1300 to ~1500 nm) and 42–40 (i.e. ~1800 to ~2000 nm), located within nIR water absorption bands, are discarded in what follows.<fnref n="4" />

We use a data analysis that resembles that of Boucher et al. (2021) albeit with a few differences described further. Our analysis consists of the following steps, independently applied to each of the 42 remaining diffraction orders, and illustrated on a given order in Fig. 2.

(i) We create an out-of-transit stellar reference spectrum, $I_{ref}$, by averaging the out-of-transit spectra interpolated in the stellar rest frame (the star moves by ~200 m s⁻¹ during the course of the transit). This reference spectrum, shifted back to the Earth rest frame, is linearly adjusted in flux to each observed spectrum $I_{obs}$, and $I_{obs}$ is divided by this best-fitting solution. This operation is then performed once again to the resulting spectra (hereafter reference-corrected spectra), but, this time, the out-of-transit reference spectrum is computed in the Earth rest frame in order to remove residual telluric contamination. Note that, contrary to Boucher et al. (2021), the rest of our data analysis is conducted in the Earth rest frame rather than in the stellar rest frame, so that the interpolation of the noise only affects the master (median) spectra and not individual spectra. Additionally, note that our data reduction enables the user to correct for planet-induced distortions of the stellar line profiles (e.g. centre-to-limb variations or the Rossiter–Maclaughlin effect; see Chiavassa & Brogi 2019). At each epoch, the code linearly fits a user-provided distorted stellar spectrum to the data prior to step (i) and normalizes the data. As planet-induced distortions of the stellar line profiles are not taken into account in our simulations, we will not give further details on its implementation in this paper and redirect the reader to a paper in preparation (Klein et al., in preparation). Note that, as shown in the panel [b] of Fig. 2, low-frequency variations in time and wavelength domains are still identifiable after correcting for the reference spectrum. These residual variations are most likely due to modal noise from the fibers (Oliva et al. 2019) and requires additional normalization.

(ii) We normalize each residual spectrum by a noise-free continuum estimate, computed using a rolling mean window, and remove outliers in the normalized spectra using a $5\sigma$ clipping procedure. These two steps are repeated until no outlier is flagged in the data. By measuring how the variance of the normalized spectra varies with the size of the rolling window, we find that a minimum width of ~23 km s⁻¹ (10 SPIRou pixels) is required to reliably average the spectrum noise. The exact size of the window has no more than a marginal impact on the recovered planet signature, provided that it is small-enough to correct for the low-frequency structures induced by modal noise in the data. In what follows, we fix the window size to 50 pixels (115 km s⁻¹) and discard the same amount of points at the ends of each diffraction order.

(iii) At this stage, outliers have been removed in the wavelength space, for each spectrum individually. However, some pixels (i.e. wavelength bins) might still exhibit large temporal variations (e.g. due to telluric contamination). In order to flag and remove these high-variance pixels, we compute the variance in the wavelength space (i.e. for each pixel), and perform an iterative parabolic fit with a $5\sigma$ outlier removal to the variance distribution.<fnref n="5" /> The outliers flagged in the process are masked out in the subsequent steps of the data processing.

(iv) Our data processing pipeline features an optional additional filtering step to correct for the variation of residual telluric absorption with airmass. Accordingly to Brogi et al. (2018), we fit a second

<footnote n="2">²Note that the potential wavelength-dependencies of $R_*$ are expected to be corrected in steps (i) and (ii) of our data cleaning procedure (see Section 3.1) and are therefore neglected in the simulations. We also assume that limb- darkening coefficients are wavelength-independent over SPIRou's spectral range.</footnote>

<footnote n="3">³This mask extends from the line centre until the relative absorption is lower than 5 per cent of the continuum level.</footnote>

<footnote n="4">⁴Due to their significant fraction of high-absorption/saturated telluric lines, keeping these orders in the analysis leads systematically to degraded results.</footnote>

<footnote n="5">⁵Note that, as a result of the blazed grating, the noise at the edge of each diffraction order is larger than in the centre, which justifies the choice of a parabolic fit</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Time series of spectra in Order 46 (1643 to 1695 nm) at subsequent steps of the data processing as detailed in Sections 3.1 and 3.2. From top to bottom: [a] Blaze- and telluric-corrected APERO-provided spectra (prior to step (i) in Section 3.1); [b] spectra corrected from the master-out template (step (i) in Section 3.1); [c] Normalized spectra (step (ii) in Section 3.1); and [d] PCA-corrected spectra (with 8 components removed; see Section 3.2). Note that time series of normalized spectra cleaned with the auto-encoder, visually similar to the bottom panel, is not shown here.</figcaption></figure>

order polynomial of the log of the intensity with airmass and remove it out (hence divide it out in intensity). However, the airmass is no more than an incomplete proxy of the water telluric absorption, expected to unpredictably change on short time-scales during the observations. As a consequence, several studies have chosen to bypass this step (e.g. de Kok et al. 2013; Boucher et al. 2021; Giacobbe et al. 2021) in favour of a more statistical approach (often based on PCA). In this study, we have kept the quadratic airmass detrending but it can be easily by-passed in our pipeline (and the order of the polynomial can also straightforwardly be changed).

The distribution of the variance in wavelength of the processed spectra is compared to the APERO-provided photon noise in Fig. 3. The dispersion of the processed spectra remains similar to the APERO estimates for the blue half of the spectrum, but is significantly higher in the *H* and *K* bands. This is most likely due to the fact that modal and thermal noises, stronger at redder wavelengths, are not including in the APERO estimation of the photon noise. To correct for residual correlated noise due to both imperfect corrections of the stellar and Earth absorption spectra and instrumental noise, we apply an additional data-driven procedure described in the following section.
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Variance at the centre of each diffraction order after processing the data as described in Section 3.1 before and after the principal component analysis (PCA) (or auto-encoder) reduction (resp. dark-blue stars and light-blue triangles). As a reference, the photon noise estimate provided by APERO is shown in gold filled circles. Each point and error bars give the mean and standard deviation across all spectra. Orders removed due to strong telluric contamination are indicated by the vertical grey bands. The position of the *YJHK* photometric bands is indicated by the horizontal magenta solid lines.</figcaption></figure>

## 3.2 PCA and auto-encoder

The last step of the data processing is conceptually different to the others in the sense that we aim at getting rid of the correlated variance in time in our data on which we have no physical priors. Defining a deterministic framework to do so seems impracticable as we expect this correlated noise to highly depend on the target and on the observing conditions. This step is therefore necessarily a data-driven approach. In practice, we have developed two different methods that we independently apply to the data. Having two methods to statistically reduce the correlated noise in the data provides additional robustness to any claim of planet atmosphere detection and prevents false positives. We stress that both methods are applied to the log of the reduced data where we can consider at first order that the total spectra is a linear combination of the planet's and noise.

Our first method is based on principal component analysis (PCA). PCA is a linear method, that recovers the dominant source of correlated variance from an eigenvector decomposition of the covariance matrix: the principal components are these eigenvectors sorted by decreasing eigenvalues. This technique has been extensively used and discussed in several HRS-based planet atmosphere studies (e.g. de Kok et al. 2013; Brogi & Line 2019; Damiano, Micela & Tinetti 2019; Boucher et al. 2021; Pelletier et al. 2021). Our second method uses a deep-learning approach based on an auto-encoder, and is a new method in the HRS exoplanet community. An auto-encoder is an artificial neural network which aims at reproducing the dominant features of a data set by encoding them into a much lower number of points through subsequent reduction matrices. Initially proposed by Hinton & Salakhutdinov (2006), it is now widely used in many fields of applied mathematics and in some astrophysical works as well (e.g. Yang & Li 2015; Čotar et al. 2021). In essence, both methods rely on reducing dimensionality by transforming our spectra into smaller

sets, but, unlike PCA, our auto-encoder is not linear and we lose information about how the data are coded.<fnref n="6" /> We now give details on the practical implementation of both methods in the further two paragraphs.

### 3.2.1 PCA implementation

In our data analysis pipeline, the PCA-based dimensional reduction is applied independently to each order in the time domain.<fnref n="7" /> The number of components associated with correlated noise and subsequently discarded for the analysis is tuned using the following procedure, illustrated on a given order in Fig. 4. For each order, we generate 5 to 10 sequences of spectra matching our observed wavelengths and times, but containing only uncorrelated Gaussian noise of level similar to our empirical photon noise estimate. To account for the larger noise level at the edges of the order, the sequences of noise are amplified by the normalized inverse of the square root of the blaze function. In principle, these sequences are free from correlated noise and can be used as references to tune our PCA. We apply PCA to each sequence of noise and store the largest eigenvalue, $S_{max}$. When we apply PCA to our observed sequences of processed spectra, any component associated with an eigenvalue significantly larger than $S_{max}$ (e.g. $2 \times S_{max}$, see the red-dotted line in Fig. 4) likely encloses a significant amount of correlated noise, and is discarded. For $15$ A, this procedure typically removes 4 PC in the blue part of the spectrum and 8 in the reddest part, which we attribute to the complex

<footnote n="6">⁶See this comparison between the PCA and auto-encoder: https://towardsdatascience.com/autoencoders-vs-pca-when-to-use-which-73de063f5d7.</footnote>

<footnote n="7">⁷The structure of the correlated noise is expected to vary from one order to the next and, therefore, the number of components associated with correlated noise and subsequently discarded has no reason to be the same for all orders.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Eigenvalues associated to each PCA component of SPIRou échelle orders 77 (~999 nm, light-blue symbols), 45 (~1707 nm, gold crosses) and 32 (~2400 nm, dark-red stars). The eigenvalues for each of the three orders have been vertically shifted for better clarity. For each order, larger symbols indicate components which have been flagged, in Section 3.2, as dominated by correlated noise and removed in the analysis, that is, 2, 5, and 7 components for orders 77, 45, and 32, respectively.</figcaption></figure>

structure of the stellar atmosphere and/or modal noise. For the hotter star HD 189 733 (see Section 5), we typically remove 2 components in the blue part and 5 to 7 in the red part. Note that injecting the synthetic planet signature to the noise maps has a (i) marginal impact on the eigenvalues and (ii) affects all the components by more- or-less the same factor, ensuring the planet atmosphere spectrum is not removed in the process. This effect is further discussed in Section 3.3.3. Finally, note that the weighted PCA framework of Delchambre (2014), is also implemented in our publicly available data processing code (via the wpea python module).

### 3.2.2 Auto-encoder implementation

As in Čotar et al. (2021), our implementation of the auto-encoder relies on four layers, which allows us to reduce the data dimension- ality from a few thousands pixels (the size of a corrected SPIRou order, which varies from 2000 to 4000 pixels after removal of bad pixels from telluric correction) to eight. Each SPIRou order has a different auto-encoding process, but, in the same order, the spectra have a common encoding and decoding matrix: as for PCA, the auto-encoder takes into account time correlated features. We train the network in the following way: we randomly select 70 per cent of the reduced spectra (panel 3 of Fig. 2), encode them through the four layers which reduce dimensionality (to 1024, then 256, then 64 then 8 pixels) and then reconstruct the original spectra. All these numbers are just optimization of the auto-encoding process, and can of course be changed (the number of layers as well). This creates an auto- encoder which is then applied on the 30 per cent remaining spectra in order to validate that their reconstruction is reliable as well. After 5000 iterations, we consider that the network has sufficiently learned and use the resulting encoding matrix as our final neural network. We then apply the final auto-encoder on the reduced spectra to create a reconstruction of the dominant feature and remove it from the observations. This is equivalent to the way we apply PCA, although we completely lose the information about the number of components order by order and how they are encoded.

The auto-encoder takes much longer to run than PCA (typically 1 min per order on a GPU against 0.1 s on a CPU for PCA) because of the learning curve. However, once the algorithm has learned and its

transformation matrices are defined for a given sequence of spectra, it takes only a few milliseconds to run it on a CPU.

## 3.3 Uncovering the planetary signature

### 3.3.1 Template matching

Once the reduced data have been cleaned through the PCA or auto-encoder, the planetary signal is still largely buried under the noise as can be seen on the last panel of Fig. 2. The use of a correlation function between a theoretical model and the reduced data has therefore been proposed since the first successful exoplanet atmosphere characterization by HRS of Snellen et al. (2010). As is done in the literature, we first create an atmospheric model at extremely high spectral resolution (between 300 000 and 1 million) with petitRADTRANS. We then use this model to build a sequence of synthetic spectra matching the observing epochs and wavelengths. This requires to Doppler-shift the model by the planet RV, $V_p$, computed at each observed planet orbital phase $\phi$ using

$$V_p(\phi) = K_p \sin(2\pi\phi) + V_{\text{sys}}, \quad (4)$$

for different values of the planet velocimetric semi-amplitude ($K_p$) and systemic Doppler shift ($V_{sys}$), and convoluting the shifted models with a Gaussian at SPIRou's resolving power. The synthetic sequence is then processed with some of the key steps described in the previous sections, as we expect the data analysis to affect the planetary spectra. This is described in Section 3.3.3.

Finally, we build sequences of processed synthetic spectra for a range of $K_p$ and $V_{sys}$ values, and compute the scalar product between each of these sequences and the observed spectra to create a correlation function (as in Boucher et al. (2021))

$$CCF = \sum_{i} \frac{d_i m_i}{\sigma_i^2}, \quad (5)$$

where $d_i$, $m_i$, and $\sigma_i$ are respectively the observed flux, the model value, and the flux uncertainty at pixel $i$ (corresponding to time $t$ and wavelength $\lambda$). Our correlation maps typically extend from $K_p = 0 \text{ km s}^{-1}$ to twice the theoretical value of $K_p$, computed from the masses of the star and planet and the semi-amplitude of the planet-induced stellar RV wobble. For $V_{sys}$, we typically explore a $200 \text{ km s}^{-1}$ wide window centred on 0. A detection can be claimed if the maximum of correlation between the reduced data and the model is obtained close to the injected semi-amplitude and Doppler shift. Following Boucher et al. (2021), we define $\sigma_i$ as the standard deviation of the value of the pixel $i$ weighted by the S/N of each spectrum

$$\sigma_i^2 = \sigma^2(t, \lambda) = \frac{\sum_t (d(t, \lambda) - \overline{d(\lambda)})^2}{N_{\text{spectra}}} \frac{\overline{\text{SNR}}}{\text{SNR}(t)}, \quad (6)$$

where the bar denotes a time average, $N_{spectra}$ is the number of spectra, and $d_i = d(t, \lambda)$. Finally, in order to convert the correlation values to significance of detection, we divide the former by the standard deviation of the correlation map in regions dominated by white noise (i.e. away from the planetary signal), as frequently done in the literature. Note that the cross-correlation analysis is only used for first-order searches of planet signatures, whereas a more statistically robust (but more time consuming) exploration of the parameter space is performed in the Bayesian framework described in Section 3.3.2.

In terms of speed, we tried to optimize the calculation of this correlation in the public code, and for a low resolution map ($50 \times 50$ points in $K_p$ and $V_{sys}$), it typically takes a couple of minutes per transit over the whole SPIRou domain on one processor. We have
not parallelized it as this is sufficiently efficient for the use we make of it, but it would be very straightforward to do so by splitting the calculations for different regions of the ($K_p$, $V_{sys}$) map.

### 3.3.2 MCMC and nested sampling

Finally, we have the possibility to robustly explore the parameter space in a Bayesian framework. We implemented two methods: a Markov Chain Monte Carlo algorithm (MCMC) based on the python module `emcee` (Foreman-Mackey et al. 2013) and a nested sampling algorithm based on the python module `pyMULTINEST` (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009; Buchner et al. 2014; Feroz et al. 2019). The second one is typically 50 times faster than the first one, and, in our tests, we never noticed any difference in the results of these two methods. We therefore present only results using the nested sampling algorithm in the rest of the paper, but having the possibility to use both samplers allows us to have independent avenues to validate the results. Both methods rely on the a likelihood $L$, defined in Brogi & Line (2019) and Gibson et al. (2020) by

$$L = \prod_i \frac{1}{\sqrt{2\pi}\zeta_i} \exp \left\{ -\frac{[am_i - d_i]^2}{b\zeta_i^2} \right\}, \quad (7)$$

where $\zeta_i$ accounts for the uncertainty of the $i^{th}$ pixel and a and b are scaling factors to account for incorrect modelling and incorrect estimation of the order variance, respectively. In the rest of this paper, a is set to 1. The main difference between the two approaches is that Brogi & Line (2019) uses a unique $\zeta$ value that does not depend on the pixel. They derive the likelihood relative to $\zeta$ and select the value that cancels the derivative, hence ensuring a maximum of likelihood, whereas Gibson et al. (2020) allows $\zeta$ to be defined pixel by pixel.

When applying the Brogi & Line (2019) likelihood, b is set to 1 and we calculate the optimal $\zeta$ for each spectrum (which is what the authors advise (M. Brogi, private communication)). Essentially, this is similar to say that our log-likelihood is the sum of the log-likelihoods of each spectrum, with a different $\zeta$ optimized for each spectra. Alternatively, with the Gibson likelihood, the values of $\zeta_i$ are defined from prior information: in this paper we chose $\zeta_i = \sigma_i$, as defined in equation (6) The b value is then obtained in a similar manner than the $\zeta$ in Brogi & Line (2019): we chose b that cancels the derivative of the likelihood, as explained in Gibson et al. (2020). We have the freedom to optimize this b value for (i) each spectrum, (ii) each order, (iii) each transit, or (iv) globally. We found that, for the simple tests presented in this paper, the four options provided very similar results.

The typical time to converge a nested sampling algorithm for a model with four parameters ($K_p$, $V_{sys}$, temperature and water mass mixing ratio) and 384 live points, which are the nested sampling equivalent of a walker in a usual MCMC framework, is 2–3 h on 36 processors. This gets drastically longer with more molecules as we face a memory issue, which is inherent to petitRADTRANS for now (P. Molliere, private communication). This problem could be overcome by pre-computing a grid of models and interpolating in this grid rather than calculating a model at each iteration (which is our choice here), but that becomes prohibitively complicated with too large a number of molecules (typically $\gtrsim 3$).

### 3.3.3 Degradation the model

Although the PCA or the auto-encoder mainly remove planet-unrelated noise, they do affect the planetary signature in the data. There is a lot of work in the literature to reproduce at best the

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Zoom on a planet atmosphere model degraded with Gibson et al. (2022)'s framework (see Section 3.3.3) when removing 0, 1, 3, or 5 principal components (PC).</figcaption></figure>

degradation by PCA onto the synthetic model so as to optimize the template matching function and/or likelihood calculations (Brogi & Line 2019; Gibson et al. 2020; Boucher et al. 2021; Pelletier et al. 2021). Skipping this step leads to significant errors in the retrieved atmosphere parameters (see the discussion in 6.3).

In order to be performed in the nested sampling algorithm, this degradation must be as fast as possible numerically. We have therefore implemented the fastest of these methods for PCA (i.e. that of Gibson et al. (2022)), which we detail further. We have not yet found an equivalently fast method for the auto-encoder, because of non-linearities in the process, and, therefore, we limit the statistical exploration of the parameter space to PCA-reduced data. However, using PCA we have realized that not degrading the model was not an issue for molecular detection when performing template matching: it only reduces the significance marginally. Our use of auto-encoder can therefore be applied to molecular detection through template matching, but is not yet suited for parameter retrieval.

Our implementation of the Gibson et al. (2022) method for PCA works as follows. During the data reduction, we store the removed eigenvectors (i.e. associated with correlated noise) for each order into a matrix U (calculated in log-space). We then multiply U by its pseudo-inverse $U^\dagger$ to create an orthogonal projector of the vector space defined by these eigenvectors. We then project the logarithm of our synthetic sequence model M on this vector space, and remove it from M. Our final degraded sequence $M'$ is therefore given by

$$M' = \exp (\log M - U U^\dagger \log M). \quad (8)$$

We stress again that U changes from one order to the another. As in Gibson et al. (2022), we do not need to take the weights into account as they can be naturally implemented in the weighted PCA algorithm. Fig. 5 shows the effect of such a degradation on an isothermal model of our HD189733 b-analogue, containing only water, when 1, 3, and 5 PCA components are removed of the data (for order 52 here).
## 3.4 Including rotation and winds

Brogi et al. (2016) defined a framework to include the effect of rotation and winds on a 1D transmission spectra in a phase dependent manner. This method is however quite time consuming and hence not suitable for a large parameter space exploration. In Appendix A, we show that the inclusion of rotation can be expressed as a double convolution at mid-transit. This provides a very fast first-order calculation of the effects of rotation (and eventually winds, see the appendix), albeit not as accurate as the framework of Brogi et al. (2016), since it does not take the phase dependence or limb darkening effects into account. In what follows, we rely on the equations presented in Appendix A to recover planetary rotation at first order in our nested sampling algorithm. In particular, note that by separating Appendix equation (A7) into its blueshifted and redshifted components, one could straightforwardly create a transmission spectrum where both hemisphere have different physical parameters. This hemispheric dichotomy is notably applied in a forthcoming work of the ATMOSPHERIX consortium (Hood et al., in preparation).

# 4 APPLICATION ON SIMULATED DATA

## 4.1 Simple isothermal model

Following the process described in Section 2.2, we first inject a simple, isothermal atmospheric model, containing only water with a volume mixing ratio of $10^{-2}$ and a temperature of 900 K, in the APERO-provided telluric-corrected spectra of Gl 15 A. As a first step, we run the cross-correlation analysis to the data processed using the different steps described in Section 3.1, but prior PCA (or auto-encoder) cleaning. Unsurprisingly, strong signatures at low semi-amplitudes and velocity shifts dominate the correlation map, confirming that residuals water lines from Gl 15 A and the Earth atmosphere are still prominent in the reduced spectra. Note that detrending the data with airmass is not sufficient to uncover the injected signal, which confirms that a PCA/auto-encoder treatment is needed.

In contrast, fair detections of the injected signals are obtained when the data are cleaned with PCA or auto-encoder. In terms of cross-correlation, the signal was found at a SNR of about 6, as shown in Figs 6 and 7. Both the auto-encoder- and the PCA-based treatments yield similar level of signal detection, thereby confirming that the auto-encoder is a reliable, robust approach to PCA, which would gain at being developed further. The relatively large significance of detection (of the order of what can be found on hot Jupiters; see Line et al. 2021) can be explained by the fact that the planet absorption template used in the cross-correlation is the same as the injected model.

In terms of parameter retrieval from the nested sampling algorithm, Fig. 8 shows the corner plot of our results. Our MCMC process has converged towards Gaussian-looking posterior densities, roughly matching the injected parameters. Temperature and water content are unsurprisingly degenerated, but we do recover the injected values in the $1\sigma$ ellipse of posteriors. Note that changing the white noise realization (by fitting the other synthetic transit) or the likelihood definition (see Section 3.3.2) only marginally affects the retrieved parameters, which remain consistent within ~1$\sigma$. This confirms that our analysis and parameter estimation process do not introduce strong biases in the retrieval. Additionally, we tested the effects of changes in the input water content and did not find systematic trends in the parameters recovered using different likelihood definitions: a same

<figure><img src="image_6.png" /><figcaption>**Figure 6.** Cross-correlation significance as a function of Doppler velocity and orbital semi-amplitude for the sequence of Gl 15 A spectra, in which an isothermal planet atmosphere model has been injected (see Section 4.1), and reduced using the procedure described in Section 3.1 and PCA-cleaning (Section 3.2).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Same as Fig. 6, this time for the encoder-reduced data.</figcaption></figure>

likelihood can overestimate the water content for one synthetic model and underestimate it for another model. On real data, we therefore recommend to gather results from different likelihoods to define conservative error bars in the retrieved parameters.

## 4.2 Including rotation

We have first tried to recover rotation from a model that did not include any, and we found that our nested sampling algorithm could not differentiate between a non-rotated model or models with equatorial rotation lower than $1 \text{ km s}^{-1}$. Then we have created a model with a planetary rotation of $3 \text{ km s}^{-1}$ at the equator, following appendix A (hence assuming the rotation axis is perpendicular to the orbit). When we try to recover this model with a non-rotating model, the best fit value of the water mixing ratio is decreased by a factor of ~30 and we also get lower values for the temperature. This is due to the fact that rotation decreases the strength of the absorption lines by spreading them on a larger width (see Fig. 9), whereas temperature and water content typically increase the line contrasts. As we are mostly sensitive to line amplitude and not shape, this creates a degeneracy between rotation and temperature/composition.

When we try to retrieve a model with rotation, we recover the correct parameters in the posteriors as shown in Fig. 10 but the mean recovered values for water and temperature are 3 and $2\sigma$.
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Posterior densities resulting from a *pyMULTINEST* retrieval of the data containing the injected planet atmosphere signal described in Section 4.1, with the simple isothermal model described in Section 4.1. The blue squares with error bars show the best fit values, with the ±1σ uncertainties, and the red line and crosses indicate the injected values.</figcaption></figure>

away from the maximum of posterior probability for water and temperature. No matter the rotation speed, we always obtained too low water content and too high temperature, showing that this is intrinsic to the analysis which exhibits a bias for large rotation rate. As demonstrated in Appendix B, we have performed a wide range of tests and demonstrated that it most likely comes from a lack of model degradation. Indeed, although our model is degraded consistently with the PCA applied to the data, the first phases of the data analysis (average stellar spectrum division, moving average normalization and airmass detrending which worsens this effect when included) are not applied to the models during retrievals. This mainly affect

models that are almost constant with time: very broadened (hence fast rotating) spectra or planets with low semi-amplitudes.<fnref n="8" /> We are working on finding the most efficient way to include this additional model degradation in the Nested Sampling algorithm in order not to increase too much the numerical cost. Until then, our analysis is biased to lower molecular content for models with large rotation rate.

<footnote n="8">⁸For real planet, this effect is probably smaller because of the variability of the planet spectra, due to intrinsic variability and 3D geometric projections</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Top: rotation kernel (arbitrary units) in plain orange line and instrumental profile in dashed line as a function of velocity similar to Brogi et al. (2016). Bottom: zoom on the absorption lines of a model containing only water, broadened at SPIRou resolution, and the same model when applying a rotation with an equatorial velocity of 3000 m s⁻¹.</figcaption></figure>

## 4.3 Including clouds

Our further test was to study the influence of clouds in our data analysis which are a major limitation for atmospheric analysis (e.g. Kreidberg et al. (2014)). HRS has the potentiality to see above the clouds when they are deep enough (Gandhi, Brogi & Webb 2020; Hood et al. 2020) and we first tried to recover synthetic planets with grey cloud deck at different pressure levels. When the cloud deck was below 0.1 bar, we recover the model roughly at the same amplitude than the non-cloudy model. This is consistent with the fact that we expect to probe pressure levels around 10–1000 Pa through water absorption. When we move the clouds higher up in the atmosphere, we typically lose one point of signal to noise detection per order of magnitude in pressure, until 0.1mbar where the SNR becomes lower than 2. However, some absorption lines are still theoretically observable (Gandhi et al. 2020) and combining several observations might allow to push this limit upwards.

We have then ran a retrieval including clouds on a model with no clouds and obtained that, even if clear models exhibit higher likelihood, models with clouds are not excluded by our analysis:

we obtain a degeneracy between water content and cloud coverage. This is not surprising as HRS is only sensitive to the variations of the atmospheric absorption, and not its absolute value. Additional constraints, such as LRS or fixed temperature can lift this degeneracy, as we will see in the application to real data on Section 5.

Finally, we have created a model with a grey cloud deck at 10mbar and applied our MULTINEST algorithm which results are shown on Fig. 11. Globally, the fit is poorer which is expected as the amplitude of planetary lines is reduced. We recover a very tight degeneracy between water and cloud top pressure spanning almost 5 orders of magnitude in water, showing that we lose our capability to obtain a precise water mixing ratio in that case without further constraints. This will be discussed again in Section 5 where the change in the temperature profile allows to lift part of this degeneracy. This confirms our test with the non-cloudy model: the use of HRS alone does not allow to lift the cloud-composition-temperature degeneracy, and additional information must be added. However, we note that our model does recover the injected model at the 1σ level, which validates our pipeline for (simple) cloudy models.
<figure><img src="image_10.png" /><figcaption>**Figure 10.** Corner plot showing the results of a *pyMULTINEST* retrieval with a model containing only water with a mass mixing ratio of $10^{-2}$ and a rotation with equatorial velocity of 3 km s$^{-1}$. The different figure elements are the same as in Fig. 8.</figcaption></figure>

# 5 APPLICATION ON REAL DATA

## 5.1 Short description of the data

In order to validate entirely our pipeline, we have applied it on real data with already published result for comparison. We have therefore used the two transits of HD 189 733 b observed by SPIRou as in Boucher et al. (2021) (hereafter B21). We shortly detail these data here, a larger discussion can be found in B21. The physical parameter of the planet are referenced in Table 1.

Two transits of HD 189 733 b were observed as part of the SPIRou Legacy Survey (SLS, PI: J.-F. Donati). The first transit was observed

on UT 2018 September 22 (hereafter Sep18), as part of SPIRou commissioning observations, and the second on UT 2019 June 15 (hereafter Jun19). The first data set consists of 2.5 h, divided into 36 exposures, where the first 21 are in transit and the remaining 15 are out-of-transit. The second data set consists of 50 exposures in total, where 24 are in transit, 12 before, and 14 after transit, for a total of ~3.5 h. The data were reduced using APERO version 0.6.132. In both observations, the airmass remains below 1.3 and even below 1.15 during the transit. The mean SNR per order ranges from 50 in the telluric contaminated region and in the bluest part of the instrument to 250 in the centre of the H band. Conditions were photometric for
<figure><img src="image_11.png" /><figcaption>**Figure 11.** Corner plot showing the results of a *pymultinest* retrieval with a model containing only water with a mass mixing ratio of 10⁻² and grey cloud deck at 10 mbar. The different figure elements are the same as in Fig. 8 with *P*<sub>clouds</sub> being the recovered cloud top pressure.</figcaption></figure>

both transit sequences, with an average seeing of around 0.82 as estimated from the guiding images.

## 5.2 Data analysis and retrieval

We have therefore applied our pipeline on these two transits of HD 189733 b in order to retrieve atmospheric signatures and compare with B21. In order to be as consistent as possible with their methods, we have used the additional telluric correction of B21 (i.e. masking telluric lines deeper than 70 per cent from the continuum level, as described in Section 3.1), and the detrending with airmass was only performed through PCA (see Section 3). However, as presented in

Section 3, our pipeline has some intrinsic differences with B21, notably (i) we interpolate only the reference spectrum, and not individual ones, and (ii) the number of PCA components to remove is decided automatically by the pipeline. In both data sets, the number of PCA removed ranges from 1 in the bluest orders to 4/5 in the reddest orders.

In the template matching algorithm, we have used the exact same model than the best model of B21. The resulting cross-correlation map is shown in Fig. 12 to be compared with their Fig. 5. The comparison is excellent: the maximum is recovered at the expected theoretical semi-amplitude (151 km s⁻¹) and the recovered Doppler shift is comparable within B21 with less than 500 m s⁻¹ of difference.
<figure><img src="image_12.png" /></figure>

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

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Top: cross-correlation map between the combined transits of HD 189 733 b analysed with PCA and the atmospheric model used in B21. The white dashed line show the theoretical position with no atmospheric Doppler shift. Middle: same with the auto-encoder. Bottom: cross-correlation significance from the PCA-reduced data for individual transits and both transits and an orbital semi-amplitude equal to the planetary semi-amplitude (151.2 km s⁻¹). The black dashed line is the 0 Doppler shift.</figcaption></figure>

We obtain a slightly higher maximum of correlation (SNR of 4.6 compared to 4 for B21) with a lower amplitude for a same non-planetary peak obtained in both our works at $K_p \sim 270$ km s⁻¹ and $V_{sys} \approx -75$ km s⁻¹, showing that our pipeline corrects better for spurious signatures. We also applied our auto-encoder on these data and found that the detection was slightly lower (SNR of 4.2) but exactly at the same position and the secondary spurious peak disappeared, confirming that it is not a physical signature. Additionally, the negative maximum of correlation next to the positive maximum, which is often recovered in studies using PCA disappears with the auto-encoder. This is promising towards a more global use of this technique: it shows that coupling PCA and auto-encoder can allow to disentangle between numerical and physical signatures which will be of real added values for planets with lower atmospheric detectability.

We have then used our nested sampling framework to retrieve parameters and compare with fig. 11 of B21. We created models of HD 189 733 b using petitRADTRANS with the water line list POKAZATEL from Exomol Tennyson et al. (2016); Polyansky et al. (2018) and a temperature profile from Guillot (2010) as in B21. The resulting posterior distributions is shown in Fig. 13 with the red crosses being the mean recovered values of B21. We see that we are perfectly consistent with their recovered parameters although we recover a higher temperature (which is more consistent with physical expectation of the temperature at the limbs of HD 189 733 b (e.g. Drummond et al. (2018))), a higher water content and deeper cloud top pressure. Comparing to our test on synthetic data, it is striking how well the deep cloud top pressure is recovered. This surprising good retrieval led us to consider an isothermal profile as shown in Appendix Fig. C1. As we expected, our algorithm then does not distinguish between a high cloud deck-high water content and deep cloud deck-low water content. We indeed see two Gaussian distributions in the water posterior density, as in Fig. 11. This shows how different variables are intricated, as we explore further in our companion paper, and how adding information on the temperature can lift the degeneracies in other parameters. As already mentioned, a combination with low resolution spectroscopy would also allow to solve this discrepancy: the observation of the slope of the continuum favours clear atmospheres (Sing et al. 2016; Barstow 2020).

In Appendix C, we show two other posterior distributions when including rotation: one where the rotation speed is imposed as the expected tidally locked value for HD 189 733 b (2.6 km s⁻¹ at the equator) and one where the rotation speed is left as a free parameter. In both cases, we did not include clouds as they only increase complexity and are disfavoured, as mentioned earlier. For the first one, Fig. C2, we obtain similar results than in Section 4.2: we recover a lower water MMR and a higher temperature with larger error bars. This globally confirms that the 1D, non-rotated water MMR is a good estimate as this is coherent with our analysis on simulated data. We also recover a higher $K_p$ and the systemic velocity is changed by 500 m s⁻¹, being then perfectly in line with B21. In the second one, Fig. C3, the error bars are expectedly larger as rotation adds another degeneracy but we still recover consistent parameters and show that the data are consistent with tidally locked rotation. In summary, this application to real data confirms the validity of our methods globally.

# 6 DISCUSSION

## 6.1 Errors on $K_p$

In our simulations (and in general HRS planet observations), the absorption lines of the planet atmosphere are largely drowned in the noise, which limits the achievable precision of our recovered
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Posterior density for our nested sampling algorithms applied on the two transits of HD 189733 b with a temperature profile from Guillot (2010). The red crosses are the B21 recovered values. The water quantity is in mass mixing ratio, contrary to B21 in volume mixing ratio.</figcaption></figure>

parameters. From our simulated data, for example Fig. 8, we see that the 3σ error on the velocity is of the order of a third of a SPIRou pixel, and the 3σ error on $K_p$ leads to a shift of half a pixel at the beginning and end of the transit. These are extremely simplified cases as the injected and recovered model are very similar and it is therefore expected that, for real data, the error on $K_p$ can be a factor of a few larger. It does however provide a good understanding of the precision of the method: being limited at the half pixel precision points towards the fact that we pre-dominantly recover the centre of the lines and not their shapes.

Another interesting aspect is that, although we have performed a lot of simulations, we never recovered a mean $K_p$ value that was lower

than the injected $K_p$ in synthetic data. In contrast, the broadening induced by the rotation can lead to a recovered mean $K_p$ more than 20 km s⁻¹ higher than the injected value. We attribute this effect to the data reduction process: the division by the median as well as the PCA/auto-encoder aims at suppressing signals that are almost constant in time over the whole sequence, hence that have low $K_p$ values. The PCA is also applied to the models, which reduce this effects, but not the first steps of the data analysis which have a lower but non-zero impact on the model. This artificially enhances the recovered $K_p$ and is coherent with the fact that this trend gets worse with higher rotation rates: the broadening of the lines make them more sensitive to the data analysis. The Doppler shift between
the beginning and end of transit in our fiducial sequence is about
12 km s⁻¹: a line broadened by a rotation kernel of a few km s⁻¹
will be more affected by the pipeline than non-broadened lines.
Finally, we note that Vsys is well recovered by our model, confirming
that blueshifts (or redshifts) in observed data will most likely be of
physical origin (i.e. atmosphere dynamics such as winds).

Globally, the fact that this technique performs better at high K<sub>p</sub> was known and expected: a better separation between planetary, stellar and telluric lines during the course of the transit increases the level of detection. This will be a challenge for example PLATO targets (Rauer et al. 2014) which will have no detectable Doppler variation between the star and planet during transit. Through our consortium, we have obtained part of the transit of HIP 41378f, an inflated Neptune mass planet on a 1.5 yr orbit (see e.g. Alam et al. (2022)), that shift by less than a meter per second from the stellar line over the course of the transit. It will be a good test to optimize our method for planets with low K<sub>p</sub>, or a proof of the absolute limitation of this method in the context of slowly displacing planets.

## 6.2 Water-temperature degeneracy

As we expected, there is a degeneracy between composition and temperature. Although their physical effect on the atmosphere is different, both these parameters affect strongly the amplitude of the lines: temperature through a change in pressure scale height and composition through a change in opacity (and to a lesser extent to the pressure scale height as well). It is however important to notice that our analysis is not biased: although the maximum of likelihood does not necessarily correspond to the injected model, the injected parameters do lie in the ellipse of recovered parameters. A naive way to reduce the impact of this degeneracy is to provide informative priors but one obviously has to be careful about their physical motivation. Essentially, we need other diagnostics to lift the degeneracies.

Two main diagnostics come into mind: (i) visible observations, where the strength of the absorption lines such as Ca⁺, Na or FeI is so much larger than it can provide tighter constraint on the rotation or temperature through the line shape, and (ii) combination with low resolution spectroscopy (LRS). LRS can provide a reference value for the strength of the absorption line, hence a tighter constraint on the temperature-composition degeneracy. This in turn will provide a more appropriate estimate of dynamical mechanisms in the planet.

## 6.3 Degrading the model

As we discussed in Section 1, we degrade the model because of the PCA treatment accordingly to the fast method developed by Gibson et al. (2022). However, as mentioned in Section 4.2 and discussed in Appendix B, we do not degrade the model from the data analysis steps prior to PCA which leads to errors for broadened models. This is work in progress and is not discussed further here.

Regarding the PCA degradation, we have made several tests to try to understand how this step was important. As can be seen on Fig. 5, the degradation first leads to reduce the depth of several absorption lines. When we do not degrade the model, we therefore expect to recover lower mixing ratio or lower temperature to mimic this effect. This is exactly what we obtained, and when we tried to recover only the water content the error could reach a factor of 20 with a non-degraded model. Degrading the model correctly is therefore of primordial importance.

Essentially, without this step of model degradation, our recovered values were much further away to the real data and sometimes

incompatible at the 3σ level. The correction provided by Gibson et al. (2022) allows to correct these effects, while only costing a small amount of calculation time. We are looking for such an efficient method to apply with the auto-encoder.

## 6.4 Perspectives of exploration

As the goal of this paper is mainly to present the pipeline, we have focused on a few examples but we could obviously not cover all the issues in the analysis of planetary atmospheres. Our companion paper, Debras et al. submitted, tackles the importance of biases and degeneracies in a few different situations. There are however many aspects that we did not include in these first ATMOSPHERIX papers.

Probably the most important is that we have not mentioned at all 3D effects, although they are known to impact the retrieval of atmospheric parameters (Caldas et al. 2019; Flowers et al. 2019; Pluriel et al. 2020; Falco et al. 2022; Pluriel et al. 2022). It would have been too large a task to explore these effects for a benchmark study, and we dedicate it to individual planet studies with the forthcoming works of the ATMOSPHERIX consortium.

As we mention in the introduction, the way forward in our understanding of planetary atmospheres is the combination of low- and HRS. The method to combine these observations efficiently has been presented and discussed in Brogi et al. (2017), and an application with SPIRou data has already been performed (Boucher et al. 2023). We are therefore working on a similar benchmark paper with the combination of space and ground-based data, notably in the goal of exploiting at best the JWST and Ariel capabilities.

Finally, we chose to focus on infrared transmission spectroscopy with SPIRou but this pipeline could in theory be easily adapted to emission/reflection spectroscopy or data from another instrument/wavelength range. Emission spectroscopy has already been performed with SPIRou observations (Pelletier et al. 2021) with similar methods and we have gathered and reduced visible data from MAROON-X with our pipeline, requiring marginal changes. Therefore, our pipeline is straightforwardly applicable to a much broader range of observations, which will be presented in the future.

## 7 CONCLUSION

In this paper, we have presented our publicly available pipeline to (i) generate synthetic SPIRou transmission spectroscopy data, (ii) reduce the data with state-of-the-art methods including PCA and the use of an auto-encoder, and (iii) analyse the data to recover the injected planetary signal either through cross-correlation or statistical exploration of the parameter space in a Bayesian framework. We have also included a fast way to include and retrieve (super-) rotation in planetary atmospheres.

By creating synthetic sequences, we demonstrated the validity of our pipeline and explored its limitations. We have first confirmed that the auto-encoder was a working, independent method to combine with PCA to recover planetary signals. We have shown that our method is unbiased for simple 1D models, but some issues remain in the retrieval for models with large rotation rates. We have explored the impact of clouds, showing that they can also bias the results and require additional constraints to be properly recovered.

Importantly, we have confirmed that degrading the model was needed to ensure a proper retrieval. We have implemented the Gibson et al. (2022) method for PCA and are still working on a fast method for the auto-encoder and for the non-PCA steps of the data analysis
as well. When the model is not degraded, the retrieved value of the mixing ratio in our tests could be more than 1 dex away from the input value, and the temperature about 200K wrong.

Finally, we have applied our pipeline on real SPIRou observations of HD 189 733 b and obtained slightly better results than the literature for the detection and characterization of the atmosphere. We recover water at a SNR greater than 4.5 with a volume mixing ratio in line with the literature and a temperature profile consistent with physical priors. We also show that we are consistent with a tidally locked rotation rate, an important result for hot Jupiters.

In conclusion, we have benchmarked our pipeline for atmospheric observations of exoplanets. A companion paper (Debras et al., submitted) tests its limitations for non-isothermal and multispecies models. With the ongoing observations of the ATMOSPHERIX consortium, we have proven to be ready for the challenge of atmospheric detection and characterization and will present results from SPIRou observations in a suite of papers to come.

## ACKNOWLEDGEMENTS

Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) which is operated from the summit of Maunakea by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. The observations at the Canada-France-Hawaii Telescope were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site. BK acknowledges funding from the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 865624, GPRV). FD thanks the CNRS/INSU Programme National de Planétologie (PNP) and Programme National de Physique Stellaire (PNPS) for funding support. This work was supported by the Action Spécifique Numérique of CNRS/INSU. This work was granted access to the HPC resources of CALMIP supercomputing centre under the allocation 2021-P21021. JFD, CM, IB, XB, AC, XD,GH, FK acknowledge funding from Agence Nationale pour la Recherche (ANR, project ANR-18-CE31-0019 SPlaSH). AM, BC, BB, SV acknowledge funding from Programme National de Planétologie (PNP) and the Scientific Council of the Paris Observatory. XD and AC acknowledge funding by the French National Research Agency in the framework of the Investissements d'Avenir programme (ANR-15-IDEX-02), through the funding of the 'Origin of Life' project of the Université de Grenoble Alpes.JL acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 679030/WHIPLASH), and from the french state: CNES, Programme National de Planétologie (PNP), the ANR (ANR-20-CE49-0009: SOUND). JFD and BK acknowledge funding from the European Research Council (ERC) under the H2020 research & innovation programme (grant agreement #740651 NewWorlds). Part of the work has also been carried out in the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). WP acknowledge financial support from the SNSF for project 200021_200726. PT acknowledges supports by the European Research Council under Grant Agreement ATMO 757858. EM acknowledges funding from FAPEMIG under project number APQ-02493-22 and research productivity grant number 309829/2022-4 awarded by the CNPq, Brazil. Finally, we thank the anonymous referee for valuable comments and suggestions throughout the paper.

## 8 DATA AVAILABILITY

The data are available upon request to the author and the code to analyse them is publicly available.

## REFERENCES

* Addison B., et al., 2019, PASP, 131, 115003

* Alam M. K., et al., 2022, ApJ, 927, L5

* Allart R., et al., 2023, A&A, 677, A164

* Alonso-Floriano F. J., et al., 2019, A&A, 621, 74

* Armstrong D. J., et al., 2016, Nature Astron., 1, 0004

* Artigau E., et al., 2018, in Holland A. D., Beletic J. eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10709, High Energy, Optical, and Infrared Detectors for Astronomy VIII. SPIE, Bellingham, p. 107091P

* Barstow J. K., 2020, MNRAS, 497, 4183

* Bauer F. F., Zechmeister M., Reiners A., 2015, A&A, 581, 117

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

* Birkby J. L., de Kok R. J., Brogi M., de Mooij E. J. W., Schwarz H., Albrecht S., Snellen I. A. G., 2013, MNRAS, 436, 35

* Boucher A., et al., 2021, AJ, 162, 233

* Boucher A., et al., 2023, MNRAS, 522, 5062

* Bouchy F., et al., 2005, A&A, 444, 15

* Brogi M., Line M. R., 2019, AJ, 157, 114

* Brogi M., Snellen I. A. G., de Kok R. J., Albrecht S., Birkby J., de Mooij E. J. W., 2012, Nature, 486, 502

* Brogi M., et al., 2016, ApJ, 817, L106

* Brogi M., Line M., Bean J., Désert J. M., Schwarz H., 2017, ApJ, 839, L2

* Brogi M., Giacobbe P., Guilluy G., de Kok R. J., Sozzetti A., Mancini L., Bonomo A. S., 2018, A&A, 615, 16

* Buchner J., et al., 2014, A&A, 564, 125

* Caldas A., Leconte J., Selsis F., Waldmann I. P., Bordé P., Rocchetto M., Charnay B., 2019, A&A, 623, 161

* Chiavassa A., Brogi M., 2019, A&A, 631, 100

* Cho J. Y. K., Skinner J. W., Thrastarson H. T., 2021, ApJ, 913, L32

* Claret A., Bloemen S., 2011, A&A, 529, 75

* Cook N. J., et al., 2022, PASP, 134, 114509

* Čotar K., et al., 2021, MNRAS, 500, 4849

* Cristofari P. I., et al., 2022, MNRAS, 516, 3802

* Cutri R. M., et al., 2003, VizieR Online Data Catalog, p. II/246

* Damiano M., Micela G., Tinetti G., 2019, ApJ, 878, L153

* Debras F., Baruteau C., Donati J.-F., 2021, MNRAS, 500, 1621

* Delchambre L., 2014, MNRAS, 446, 3545

* de Kok R. J., Brogi M., Snellen I. A. G., Birkby J., Albrecht S., de Mooij E. J. W., 2013, A&A, 554, 82

* Donati J. F., et al., 2020, MNRAS, 498, 5684

* Drummond B., et al., 2018, ApJ, 869, L28

* Falco A., Zingales T., Pluriel W., Leconte J., 2022, A&A, 658, 41

* Feroz F., Hobson M. P., 2008, MNRAS, 384, 449

* Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601

* Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2019, Open J. of Astrophysics, 2, 10

* Flowers E., et al., 2019, AJ, 157, 209

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

* Fouqué P., et al., 2018, MNRAS, 475, 1960

* Fulton B. J., et al., 2017, AJ, 154, 109

* Gandhi S., Brogi M., Webb R. K., 2020, MNRAS, 498, 194

* Giacobbe P., et al., 2021, Nature, 592, 205

* Gibson N. P., et al., 2020, MNRAS, 493, 2215

* Gibson N. P., et al., 2022, MNRAS, 512, 4618

* Guillot T., 2010, A&A, 520, 27

* Guillot G., et al., 2022, A&A, 665, 104

* Hinton G. E., Salakhutdinov R. R., 2006, Science, 313, 504

* Hobson M. J., et al., 2021, A&A, 648, 48

* Hood C. E., et al., 2020, AJ, 160, 198
Horne K., 1986, PASP, 98, 609

Howard A. W., et al., 2014, ApJ, 794, L51

Kasper D., et al., 2023, AJ, 165, 7

Komacek T. D., Showman A. P., 2020, ApJ, 888, L2

Kreidberg L., 2015, PASP, 127, 1161

Kreidberg L., et al., 2014, Nature, 505, 69

Line M. R., et al., 2021, Nature, 598, 580

Mollière P., et al., 2019, A&A, 627, 67

Oliva E., Rainer M., Tozzi A., Sanna N., Iuzzolino M., Brucalassi A., 2019, A&A, 632, 21

Pelletier S., et al., 2021, AJ, 162, 73

Pinamonti M., et al., 2018, A&A, 617, 104

Pluriel W., Zingales T., Leconte J., Parmentier V., 2020, A&A, 636, 66

Pluriel W., et al., 2022, A&A, 658, 42

Polyansky O. L., Kyuberis A. A., Zobov N. F., Tennyson J., Yurchenko S. N., Lodi L., 2018, MNRAS, 480, 2597

Rauer H., et al., 2014, Exp. Astron., 38, 249

Reid I. N., Hawley S. L., Gizis J. E., 1995, AJ, 110, 1838

Rosenthal L. J., et al., 2021, ApJS, 255, 8

Sing D. K., et al., 2016, Nature, 529, 59

Snellen I. A. G., et al., 2010, Nature, 465, 1049

Tennyson J., et al., 2016, J. Mol. Spectrosc., 327, 73

Tinetti G., et al., 2021, preprint (arXiv:2104.04824)

Udry S., Santos N. C., 2007, ARA&A, 45, 397

Valencia D., et al., 2013, ApJ, 775, L10

Yang T., Li X., 2015, MNRAS, 452, 158

## APPENDIX A: CONVOLUTING WINDS

In order to derive simple equations for including the effects of planet rotation and winds in our retrieval, one has to remember that the observable quantity is the transit depth, that is, the area of the projected planetary disc over the stellar surface, as a function of frequency. The elementary area of a curve $r(\theta)$ between $\theta$ and $\theta + d\theta$ is the area of the elementary triangle, hence $1/2r \times rd\theta$. Assuming that the centre of the planet is the centre of our polar coordinate system, the projected area of the planet at a frequency $v$, defined as the curve $R_p(v, \theta)$ on the stellar disc is

$$A(v) = \frac{1}{2} \int_{-\pi}^{\pi} R_p^2(v, \theta) d\theta, \qquad (A1)$$

where $-\pi$ corresponds to the evening limb, $-\pi/2$ is the South pole, 0 the morning limb, and $\pi/2$ the North pole of the planet. We can define an effective radius, $R_{\text{eff}}$, such that

$$A(v) = \pi R_{\text{eff}}^2(v) = \frac{1}{2} \int_{-\pi}^{\pi} R_p^2(v, \theta) \cdot d\theta. \qquad (A2)$$

In this appendix, we make the assumption that the planet is properly characterized as a 1D object whose radius only depends on the wavelength Doppler-shifted by rotation. However, we note that it would be very easy to decompose the planet into different regions within this framework (e.g. equatorial region and high latitudes), hence creating a pseudo-2D planet. Similarly, a latitude-dependent weight could be added on this integral to account for limb darkening effect, but we restrict to simple considerations here.

Assuming that the planet rotates as a solid body with equatorial speed $v_0$, the observed planetary radius at mid-transit in the stellar frame at a frequency $v$ is simply the integral of the planetary radius Doppler-shifted by rotation

$$R_{\text{eff}}^2(v) = \frac{1}{2\pi} \int_{-\pi}^{\pi} R_p^2 \left( v \cdot \left[ 1 + \frac{v_0 \cos \theta}{c} \right] \right) \cdot d\theta, \qquad (A3)$$

where $c$ is the speed of light. Assuming North-South symmetry in the plane, we can remove the factor of 1/2 and integrate from 0 to $\pi$,

that is, morning limb to evening limb. For computational efficiency, we wish to express this integral as a convolution and we need to express it as a function of speed instead of wavelength in order to convolve it with the instrumental profile as in Brogi et al. (2016). For a given frequency chunk of mean $v_0$ (e.g. spectrograph diffraction order), the corresponding velocity $v$ can be simply expressed as

$$v = c \left( \frac{v}{v_0} - 1 \right), \qquad (A4)$$

which leads to

$$v \left( 1 + \frac{v_0 \cos \theta}{c} \right) = v_0 \left( 1 + \frac{v + v_0 \cos \theta}{c} + \frac{vv_0 \cos \theta}{c^2} \right) \qquad (A5)$$

$$\approx v_0 \left( 1 + \frac{v + v_0 \cos \theta}{c} \right). \qquad (A6)$$

Within a SPIRou order, $v$ typically goes from $-5000 \text{ km s}^{-1}$ to $5000 \text{ km s}^{-1}$ and the equatorial rotation speed reaches no more than a few km s$^{-1}$. Therefore, the approximation made in equation (A6) is verified and the error on the frequency is lower than the spectral resolution of SPIRou at the border of the order.<fnref n="9" /> We can then write

$$R_{\text{eff}}^2(v) = \frac{1}{\pi} \int_{0}^{\pi} R_p^2(v + v_0 \cos \theta) d\theta. \qquad (A7)$$

The observed squared planet radius $R_{\text{obs}}^2$ is obtained by convolving $R_{\text{eff}}^2$ with the instrumental profile. We therefore fall back onto the numerical result of Brogi et al. (2016) by writing

$$R_{\text{obs}}^2(v) = \frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^{\infty} R_{\text{eff}}^2(t) e^{-\frac{(v-t)^2}{2\sigma^2}} dt. \qquad (A8)$$

The rotation kernel $K(v)$ at mid-transit assuming no limb darkening is therefore

$$R_{\text{observed}}^2(v) = \frac{1}{\pi} (R_p^2 * K)(v) \qquad (A9)$$

$$= \frac{1}{\pi} \left[ R_p^2(t) * \frac{1}{\sigma\sqrt{2\pi}} \int_0^{\pi} e^{-\frac{(v+t)^2}{2\sigma^2}} d\theta \right] (v). \qquad (A10)$$

Although expression (A10) is quite simple, we can still improve its computational efficiency by writing it as two convolutions, which are extremely fast numerical operations. Denoting

$$x = -v_0 \cos \theta, \qquad (A11)$$

we have

$$d\theta(x) = \frac{1}{v_0 \sqrt{1 - (x/v_0)^2}} dx, \qquad (A12)$$

which is properly defined for $\theta \in [0, \pi]$. We can therefore modify equation (A7) to obtain

$$R_{\text{eff}}^2(v) = \frac{1}{v_0 \pi} \int_{-v_0}^{v_0} R_p^2(v-x) \frac{1}{\sqrt{1-(x/v_0)^2}} dx \qquad (A13)$$

$$= \frac{1}{v_0 \pi} \left[ R_p^2(x) * \frac{1}{\sqrt{1 - (x/v_0)^2}} \right] (v), \qquad (A14)$$

where the convolution is defined between $-v_0$ and $v_0$ only.

Note that, in practice, $d\theta$ is not defined in $x = v_0$. However, since equation (A14) is a Riemann sum, only the following term of the integral cannot be computed

$$\frac{1}{v_0 \pi} \int_{v_{\text{final}}}^{v_0} R_p^2(v-x) \frac{1}{\sqrt{1-(x/v_0)^2}} dx, \qquad (A15)$$

<footnote n="9">⁹Note in passing that such approximation is also implicitly performed in Brogi et al. (2016) although not explicitly written.</footnote>
where $v_{final}$ is the edge of our array of integration. We can then assume that $R_p$ is constant in this interval [$v_{final}$, $v_0$] and simply calculate the theoretical value that is missing and correct it in our code for all frequencies. But actually, to a very good approximation, one can also assimilate this correction to a constant correction at all frequencies within an order as the variation of $R_p$ are negligible in front of the mean of $R_p$: the first order error is a constant shift of the recovered radius. As the input radius and rotated radius should have the same mean (easily verified theoretically), it is an excellent approximation to simply calculate this convolution and shift the result to its expected theoretical mean.<fnref n="10" />

Finally, the observed radius is simply the convolution by the instrumental profile as expressed earlier

$$R_{\text{obs}}^{2}(\text{v}) = \left( R_{\text{eff}}^{2}(t) * \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{t^{2}}{2\sigma^{2}}} \right) (\text{v}) \qquad (\text{A16}) \\ = \frac{1}{v_{0}\pi} \left\{ \left[ R_{p}^{2}(x) * \frac{1}{\sqrt{1 - (x/v_{0})^{2}}} \right] (t) * e^{-\frac{t^{2}}{2\sigma^{2}}} \right\} (\text{v}). \qquad (\text{A17})$$

We found that this expression was the fastest numerically speaking, and allows us to include a simple prescription for rotation in our parameter space exploration. Its comparison with Brogi et al. (2016) is excellent. Additionally, as mentioned previously, the inclusion of super-rotation is extremely easy: one just has to separate the first convolution into a sum of convolutions onto smaller latitude ranges and adding a shift which is latitude dependent. In the same vein, one can separate the planet into two hemispheres of different composition/temperature profile straightforwardly.

## APPENDIX B: BIASES WITH STRONG ROTATION AND LOW SEMI-AMPLITUDE

In Section 4.2 and Fig. 10, we have shown that the retrieved water content was $3\sigma$ smaller than the injected model. We have tried different values for the rotation rate and water content and always obtained smaller retrieved water content when rotation was included, pointing towards a bias in the analysis. We have therefore performed many tests, listed further, and came to the conclusion that this issue arises from the data analysis, prior to PCA, which is not applied to the model, and is more important for signatures that are almost constant with time. Indeed

(i) we tested whether it was a SNR issue: rotation lowers the amplitude of molecular lines and hence decreases the SNR. We have therefore used models with much lower volume mixing ratio of water but no rotation which leads to comparable SNR with our high-VMR rotated model. There was no biases in our low-VMR models: the retrieval was worsen (higher error bars) but the mean retrieved values much closer to the injected values.

(ii) We have then tested the influence of rotation and not surprisingly, the higher the rotation rate the larger the mismatch between the injected and retrieved planet properties. However, the rotation kernel is not the origin of the problem: when we simply include a Gaussian broadening, we obtain similar biases.

(iii) We have also tested whether the error on the retrieved parameters could be due to an erroneous PCA correction. We have created models with a gravity divided by 4 (hence much larger atmospheric depth), which allowed to detect the molecular signatures without using PCA. The same trend was found (although to lesser extent): broadened models were biased.

(iv) We noticed that the bias is increased for reduced data that contained a step of airmass detrending prior to PCA, compared to models where the airmass correction is performed by PCA.

(v) We have created m models with a planetary radial velocity shift divided or multiplied by 4. When we increase (decrease) the semi-amplitude, we reduce (increase) the bias. Typically, with $K_p = 400$ km s⁻¹ and a rotation speed of 3 km s⁻¹ at the equator, there is no bias anymore. Hence, this effect is more important for planetary models that don't shift much in wavelength with time.

(vi) To confirm further that this is not just a SNR issue we have included a planet with a $K_p$ multiplied by 4, increasing the SNR, but with 10 times less water in the atmosphere, decreasing the SNR. Cross-correlation detects this planet with a similar significance than our reference model. In this case, we don't retrieve biases on the recovered parameters, confirming that the limiting factor is not the level of detection.

(vii) Finally, we have included the a and b value of equation (7) in the retrieval which were set to 1 in the rest of this paper. We observed that the retrieved b was not equal to 1 but does not affect the retrieval at all (only the value of the likelihood). Alternatively, the retrieved a was very close to 1 for non-broadened models and deviates largely from 1 (converging to 0.5) for broadened models, showing that there is a discrepancy between the injected and recovered model. Interestingly, when we fit the a value, the bias is lowered but that is not a satisfactory solution to us and aim at providing the most reliable data analysis.

In summary, the bias arises for broadened models, especially when the velocity shift of the planet's atmosphere with respect to the star is lower than a few resolution elements over the course of the transit and disappears for synthetic planets with large semi-amplitudes. It is increased when an airmass detrending step is included, but does not depend on PCA. Hence everything points toward the fact that the first phases of the data analysis, which are not applied consistently to degrade the models during retrieval are responsible for this bias. We are working on efficient techniques to include them.

## APPENDIX C: ADDITIONAL POSTERIOR DENSITIES FOR HD 189733 B

<footnote n="10">¹⁰Care must be taken in that case when including super-rotation or other latitudinal dependent effect.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure C1.** Same as Fig. 13 but with an isothermal model. The temperature in red is the 1 bar temperature retrieved with the Guillot (2010) profile in B21.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure C2.** Same as Fig. 13 but with the inclusion of a rotation in our model, with equatorial speed of 2.6 km s⁻¹.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

**Figure C3.** Same as Fig. 13 but when the rotation is let as a free parameter and we display the rotation speed at the equator here. The tidally locked value is indicated by a red cross.

This paper has been typeset from a TEX/LATEX file prepared by the author.