# Mass-loss rate and local thermodynamic state of the KELT-9 b thermosphere from the hydrogen Balmer series

A. Wyttenbach, P. Mollière, D. Ehrenreich, H. M. Cegla, V. Bourrier, C. Lovis, Laurent Pino, R. Allart, J. V. Seidel, H. J. Hoeijmakers, et al.

► To cite this version:

A. Wyttenbach, P. Mollière, D. Ehrenreich, H. M. Cegla, V. Bourrier, et al. Mass-loss rate and local thermodynamic state of the KELT-9 b thermosphere from the hydrogen Balmer series. Astronomy and Astrophysics - A&A, 2020, 638, pp.A87. 10.1051/0004-6361/201937316. hal-02903834

HAL Id: hal-02903834

https://hal.science/hal-02903834

Submitted on 21 Jul 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Mass-loss rate and local thermodynamic state of the KELT-9 b thermosphere from the hydrogen Balmer series*

A. Wyttenbach<fnref n="1" />,<fnref n="2" />,*, P. Mollière<fnref n="1" />,<fnref n="3" />, D. Ehrenreich<fnref n="4" />, H. M. Cegla<fnref n="4" />,****, V. Bourrier<fnref n="4" />, C. Lovis<fnref n="4" />, L. Pino<fnref n="5" />, R. Allart<fnref n="4" />, J. V. Seidel<fnref n="4" />, H. J. Hoeijmakers<fnref n="4" />,<fnref n="6" />, L. D. Nielsen<fnref n="4" />, B. Lavie<fnref n="4" />, F. Pepe<fnref n="4" />, X. Bonfils<fnref n="2" />, and I. A. G. Snellen<fnref n="1" />

Received 13 December 2019 / Accepted 25 April 2020

## ABSTRACT

KELT-9 b, the hottest known exoplanet, with $T_{eq} \sim 4400$ K, is the arch type of a new planet class known as ultra-hot Jupiters. These exoplanets are presumed to have an atmosphere dominated by neutral and ionized atomic species. In particular, $H\alpha$ and $H\beta$ Balmer lines have been detected in the KELT-9 b upper atmosphere, suggesting that hydrogen is filling the planetary Roche lobe and escaping from the planet. In this work, we detected $\delta$ Scuti-type stellar pulsation (with a period $P_{puls} = 7.54 \pm 0.12$ h) and studied the Rossiter-McLaughlin effect (finding a spin-orbit angle $\lambda = -85.01^{\circ} \pm 0.23^{\circ}$) prior to focussing on the Balmer lines ($H\alpha$ to $H\zeta$) in the optical transmission spectrum of KELT-9 b. Our HARPS-N data show significant absorption for $H\alpha$ to $H\delta$. The precise line shapes of the $H\alpha$, $H\beta$, and $H\gamma$ absorptions allow us to put constraints on the thermospheric temperature. Moreover, the mass loss rate, and the excited hydrogen population of KELT-9 b are also constrained, thanks to a retrieval analysis performed with a new atmospheric model. We retrieved a thermospheric temperature of $T = 13\,200_{-720}^{+800}$ K and a mass loss rate of $\dot{M} = 10^{12.8 \pm 0.3}$ g s$^{-1}$ when the atmosphere was assumed to be in hydrodynamical expansion and in local thermodynamic equilibrium (LTE). Since the thermospheres of hot Jupiters are not expected to be in LTE, we explored atmospheric structures with non-Boltzmann equilibrium for the population of the excited hydrogen. We do not find strong statistical evidence in favor of a departure from LTE. However, our non-LTE scenario suggests that a departure from the Boltzmann equilibrium may not be sufficient to explain the retrieved low number densities of the excited hydrogen. In non-LTE, Saha equilibrium departure via photo-ionization, is also likely to be necessary to explain the data.

**Key words.** planetary systems – planets and satellites: atmospheres – planets and satellites: individual: KELT-9 b – techniques: spectroscopic – instrumentation: spectrographs – methods: observational

## 1. Introduction

Observations of exoplanet atmospheres allow us to characterize several of these remote worlds. Indeed, the measurement of planetary spectra are rich in information about these planets' physical and chemical conditions. Thus, the study of exoplanetary atmospheres are motivated by the possibility that they can provide crucial hints about their origin and evolution.

When observing the primary transit of an exoplanet as a function of the wavelength, we measure a transmission spectrum. The transmission spectrum of an exoplanet bears the signature of its atmospheric constituents. The possibility of detecting atomic or molecular species through transmission spectroscopy was theorized early on (Seager & Sasselov 2000; Brown 2001). Since giant close-in planets (Mayor & Queloz 1995), with hot and extended atmospheres, are easy targets to probe in transmission, early observations of hot Jupiter atmospheres were

achieved with sodium detections (Charbonneau et al. 2002; Redfield et al. 2008; Snellen et al. 2008). Recent results suggest that high-resolution (typically $\lambda/\Delta\lambda \sim 10^5$) transmission spectra of hot Jupiters probe their thermospheres, a key region of the upper atmosphere, in order to understand irradiated planets and allow us to achieve precise temperature and wind measurements because the lines (e.g., the Na I D doublet, the K I D doublet, the H I Balmer lines $H\alpha, H\beta, \gamma$, the He I ($^2$S) triplet) are resolved and studies benefit from large signals (e.g., Wyttenbach et al. 2015, 2017; Casasayas-Barris et al. 2017, 2018, 2019; Khalafinejad et al. 2017, 2018; Seidel et al. 2019; Yan & Henning 2018; Cauley et al. 2019; Allart et al. 2018, 2019; Nortmann et al. 2018; Salz et al. 2018; Alonso-Floriano et al. 2019; Keles et al. 2019; Chen et al. 2020). Another landmark achievement obtained in the optical domain is the detection, with cross-correlation, of iron and other metals in the atmosphere of several ultra-hot Jupiters (Hoeijmakers et al. 2018, 2019; Borsa et al. 2019; Bourrier et al. 2020a; Cabot et al. 2020; Gibson et al. 2020; Nugroho et al. 2020; Stangret et al. 2020; Ehrenreich et al. 2020).

Ultra-hot Jupiters are the hottest and the most irradiated hot Jupiters (with a planetary equilibrium temperature above 2000–2500 K, Fig. 1). Well-known members of that category

<footnote n="1">¹ Leiden Observatory, Leiden University, Postbus 9513, 2300 RA Leiden, The Netherlands
e-mail: aurelien.wyttenbach@univ-grenoble-alpes.fr</footnote>

<footnote n="2">² Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="3">³ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="4">⁴ Geneva Observatory, University of Geneva, Ch. des Maillettes 51, 1290 Versoix, Switzerland</footnote>

<footnote n="5">⁵ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands</footnote>

<footnote n="6">⁶ University of Bern, Center for Space and Habitability, Sidlerstrasse 5, 3012 Bern, Switzerland</footnote>

<footnote>* Based on observations made at TNG (Telescopio Nazionale Galileo) telescope with the HARPS-N spectrograph under program A35DDT4 and OPTICON 2018A/038.</footnote>

<footnote>** Fellow of the Swiss National Science Foundation (SNSF).</footnote>

<footnote>*** CHEOPS Fellow, SNSF NCCR-PlanetS.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Distribution of planetary masses as a function of the total power received by the planet from the incident bolometric stellar radiation, expressed in units of power received on Earth at the top of the atmosphere ($P_{\oplus}$ = 174 PW = $1.74 \times 10^{24}$ erg s$^{-1}$). The marker size scales with the planet density, and the color with the planetary equilibrium temperature. Physical properties of exoplanets were extracted from the Extrasolar Planets Encyclopedia (exoplanet.eu) in December 2019. The lower left corner represents the evaporation desert. KELT-9 (HD 195689), encircled in black, stand out of all other planetary population in term of its received power and equilibrium temperature.</figcaption></figure>

include KELT-9 b (Gaudí et al. 2017, the hottest exoplanet discovered so far around a main-sequence star, with $T_{eq} \ge 4000$ K, see Table 1), the newly discovered MASCARA-4 b/bRing-1 b (Dorval et al. 2020), as well as targets often observed due to their bright host stars, namely, WASP-189 b, WASP-33 b, MASCARA-1 b, WASP-12 b, WASP-103 b, WASP-18 b, WASP-121 b, KELT-20 b/MASCARA-2 b, HAT-P-7 b, WASP-76 b, etc. Hot exoplanet atmospheres are expected to be near chemical equilibrium from $T \ge 2500$ K (Miguel & Kaltenegger 2014; Lothringer et al. 2018; Kitzmann et al. 2018). Under such assumptions, planet atmospheres are expected to be cloud-free on the day-side, to have strong inverted temperature-pressure profiles, and to be nearly barren of molecules, apart from H$_2$, CO, and possible traces of H$_2$O, TiO, and OH. Because of thermal dissociation and ionization, mostly neutral and ionized atomic hydrogen and helium, along with such metals as iron, sodium, magnesium, titanium, etc., are expected to be found (Kitzmann et al. 2018; Lothringer et al. 2018; Lothringer & Barman 2019; Mansfield et al. 2020). In this type of atmosphere, the important source of H$^-$ bound-free and free-free opacities acts as a dominant continuum opacity (Arcangeli et al. 2018; Parmentier et al. 2018; Bourrier et al. 2020b). These predictions can change due to heat transport and circulation, depending on whether we observe the day or the night side, or the terminator (Parmentier et al. 2018; Bell & Cowan 2018). For example, in emission spectroscopy, these features of ultra-hot Jupiters are mixed up with the H$_2$O and TiO emission features (that may be absent due to dissociation), which make emission spectra appear as a continuum (see Parmentier et al. 2018; Arcangeli et al. 2018; Mansfield et al. 2018; Kreidberg et al. 2018). In transmission spectroscopy, the probed pressures are lower, thus increasing the presence of atomic species because there is a higher rate of thermal dissociation. These higher rates are produced even without taking disequilibrium processes into account. In particular, the effect of the H$^-$ opacity increases, playing a more crucial role in the modeled transmission spectra (Kitzmann et al. 2018).

<table><caption>**Table 1.** Adopted values for the orbital and physical parameters of the KELT-9 system.</caption><thead><tr><th>Parameter</th><th>Value <sup>(a)</sup></th></tr></thead><tbody><tr><td>Star: KELT-9 (HD 195689)</td><td></td></tr><tr><td>V mag</td><td>7.55</td></tr><tr><td>Spectral type</td><td>B9.5-A0 v</td></tr><tr><td>M<sub>★</sub></td><td>1.978 ± 0.023 M<sub>☉</sub><sup>(b)</sup></td></tr><tr><td>R<sub>★</sub></td><td>2.178 ± 0.011 R<sub>☉</sub><sup>(b)</sup></td></tr><tr><td>ρ<sub>★</sub></td><td>0.2702 ± 0.0029 g cm<sup>-3</sup> <sup>(b)</sup></td></tr><tr><td>T<sub>eff</sub></td><td>9 600 ± 400 K <sup>(c)</sup></td></tr><tr><td>log g<sub>★</sub></td><td>4.058 ± 0.010 cgs <sup>(b)</sup></td></tr><tr><td>v sin i<sub>*</sub> Spectro</td><td>112.60 ± 0.04 km s<sup>-1</sup> <sup>(d)</sup></td></tr><tr><td>v sin i<sub>*</sub> RM (SB)</td><td>122.5 ± 0.6 km s<sup>-1</sup> <sup>(d)</sup></td></tr><tr><td>v sin i<sub>*</sub> RM (DR)</td><td>116.9 ± 1.8 km s<sup>-1</sup> <sup>(d)</sup></td></tr><tr><td>Transit</td><td></td></tr><tr><td>T<sub>0</sub></td><td>2 457 095.68572 ± 0.00014 BJD<sub>tdb</sub></td></tr><tr><td>P</td><td>1.4811235 ± 0.0000011 d</td></tr><tr><td>(R<sub>p</sub>/R<sub>★</sub>)<sup>2</sup></td><td>0.00677 ± 0.000072</td></tr><tr><td>ΔT<sub>1-4</sub></td><td>3.91584 ± 0.00384 h</td></tr><tr><td>b</td><td>0.177 ± 0.014 R<sub>★</sub></td></tr><tr><td>a</td><td>0.0346 ± 0.001 au</td></tr><tr><td>i</td><td>86.79° ± 0.25°</td></tr><tr><td>λ SB</td><td>-85.65 ± 0.09° <sup>(d)</sup></td></tr><tr><td>λ DR</td><td>-85.01 ± 0.23° <sup>(d)</sup></td></tr><tr><td>linear LD u<sub>1</sub></td><td>0.35 (adopted)</td></tr><tr><td>Radial velocity</td><td></td></tr><tr><td>K<sub>1</sub> = K<sub>★</sub></td><td>276 ± 79 m s<sup>-1</sup></td></tr><tr><td>K<sub>2</sub> = K<sub>p</sub></td><td>234.24 ± 0.90 km s<sup>-1</sup> <sup>(b)</sup></td></tr><tr><td>e</td><td>0 (adopted)</td></tr><tr><td>ω</td><td>90° (adopted)</td></tr><tr><td>γ</td><td>-17 740 ± 40 m s<sup>-1</sup> <sup>(b)</sup></td></tr><tr><td>Planet: KELT-9 b (HD 195689b)</td><td></td></tr><tr><td>M<sub>p</sub></td><td>2.44 ± 0.70 M<sub>J</sub> <sup>(b)</sup></td></tr><tr><td>R<sub>p</sub></td><td>1.783 ± 0.009 R<sub>J</sub> <sup>(b)</sup></td></tr><tr><td>ρ<sub>p</sub></td><td>0.53 ± 0.15 g cm<sup>-3</sup> <sup>(d)</sup></td></tr><tr><td>g<sub>surf</sub></td><td>1990 ± 600 cm s<sup>-2</sup> <sup>(d)</sup></td></tr><tr><td>T<sub>eq</sub></td><td>4360 ± 200 K <sup>(e)</sup></td></tr><tr><td>H<sub>0</sub></td><td>1400 ± 400 km <sup>(e)</sup></td></tr><tr><td>2R<sub>p</sub>H<sub>0</sub>/R<sub>★</sub><sup>2</sup></td><td>150 ± 50 ppm <sup>(e)</sup></td></tr></tbody></table>

**Notes.** <sup>(a)</sup>The parameters and errors are initially taken or computed from the extended data table 3 in Gaudí et al. (2017). <sup>(b)</sup>Updated values from Hoeijmakers et al. (2019). <sup>(c)</sup>Updated values from Borsa et al. (2019). <sup>(d)</sup>Updated values from this work. <sup>(e)</sup>For the planetary equilibrium temperature $T_{eq}$ and the lower atmosphere scale-height $H_0$, we assume an albedo $A = 0$, a redistribution factor $f = 1$ and a mean molecular weight $\mu_{atm} = 1.3$ (atomic H and He atmosphere). We also consider the planet distance to the stellar surface to be $o = a - R_★$.

Recent observations of KELT-9 b in transmission spectroscopy bring in key aspects for our understanding of ultra-hot Jupiters. The detections of Fe I, Fe II, Ti II (but not of Ti I; Hoeijmakers et al. 2018) are confirming the predicted presence of neutral and ionized species in ultra-hot Jupiters. In Cauley et al. (2019); Hoeijmakers et al. (2019); Yan et al. (2019); Turner et al. (2020), more neutral and ionized atomic species are reported, such as Na I, Ca I, Ca II, Mg I, Cr II, Sc II, and Y II. In KELT-9 b, ionized species are more abundant than neutral species, which is
<table><caption>**Table 2.** Log of HARPS-N observations.</caption><thead><tr><th></th><th>Date</th><th># Spectra <sup>(a)</sup></th><th>Exp. [s]</th><th>Airmass <sup>(b)</sup></th><th>Seeing</th><th>S/N 400 nm <sup>(c)</sup></th><th>S/N 600 nm <sup>(c)</sup></th><th>S/N Hα core <sup>(d)</sup></th><th>S/N Hβ core <sup>(d)</sup></th><th>S/N Hγ core <sup>(d)</sup></th></tr></thead><tbody><tr><td>Night 1</td><td>2017-07-31</td><td>49 (22/27)</td><td>600</td><td>1.5–1.02–1.7</td><td>–</td><td>60–131</td><td>98–189</td><td>60–120</td><td>64–128</td><td>51–105</td></tr><tr><td>Night 2</td><td>2018-07-20</td><td>46 (22/24)</td><td>600</td><td>1.7–1.02–1.4</td><td>–</td><td>27–79</td><td>70–160</td><td>44–100</td><td>64–114</td><td>41–79</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>In parentheses: the number of spectra taken during the transit and outside the transit, respectively. <sup>(b)</sup>The three figures indicate the airmass at the start, middle and end of observations, respectively. <sup>(c)</sup>The range of signal-to-noise (S/N) ratio per pixel extracted in the continuum near 400 nm and 600 nm. <sup>(d)</sup>The range of S/N ratio in the line cores of the Hα, Hβ and Hγ lines.

expected for such high temperatures. However, the absorptions measured for ionized species are not reproduced by the models (Hoeijmakers et al. 2019; Yan et al. 2019). Those measurements probe altitudes between 1.01 and 1.1 planetary radii (1–15 scale-height or $10^{-3}–10^{-6}$ bar), thus, in the lower thermosphere. A complete exploration of species absorbing at optical wavelengths in ultra-hot Jupiter is very constraining in terms of aeronomic conditions. Indeed, the additional detections of Hα by Yan & Henning (2018); Turner et al. (2020), and of Hα, Hβ by Cauley et al. (2019), also bring in important clues about the KELT-9 system because measuring absorption from specific strong opacity lines, such as the sodium Na I D doublet or the hydrogen Balmer lines, allow us to probe higher layers (≥1.2 planetary radii) in the thermosphere (Christie et al. 2013; Heng et al. 2015; Huang et al. 2017). So far Hα has only been observed in six exoplanets, four of which are classed as ultra-hot Jupiters (Jensen et al. 2012, 2018; Cauley et al. 2015, 2016, 2017a,b,c, 2019; Barnes et al. 2016; Casasayas-Barris et al. 2018, 2019; Cabot et al. 2020; Chen et al. 2020). The detection of an optically thick Hα absorption in the atmosphere of KELT-9 b, which takes place at up to ~1.6 planetary radii (at ~10 000 K in the thermosphere), provides hints that the hydrogen is filling up the planetary Roche lobe (1.95 R<sub>p</sub>) and is escaping from the planet at a rate of ~10<sup>12</sup> g s<sup>-1</sup>. This confirms that such planets undergo a process of evaporation (Sing et al. 2019). However, it is important to get a more precise idea of the mass loss rate of such objects to understand the role of the stellar irradiation and its processes of interaction with the planet atmosphere (Fossati et al. 2018; Allan & Vidotto 2019; García Muñoz & Schneider 2019).

This paper reports on the results of the Sensing Planetary Atmospheres with Differential Echelle Spectroscopy (SPADES) program (Bourrier et al. 2017b, 2018; Hoeijmakers et al. 2018, 2019). This program, performed with the HARPS-N spectrograph, is a twin of the Hot Exoplanetary Atmospheres Resolved with Transit Spectroscopy (HEARTS) program based on the use of the HARPS spectrograph (Wyttenbach et al. 2017; Seidel et al. 2019; Bourrier et al. 2020a). Throughout this work, we introduce the python CHESS code (the Characterization of Exoplanetary and Stellar Spectra code) and its sub-modules, which allow us to measure, separate, and model the stellar and exoplanet spectra.

The HARPS-N transit observations<fnref n="1" /> are described in Sect. 2. In Sect. 3, we propose a comprehensive analysis of spectroscopic transit data and present an improved method to extract transmission spectra. Furthermore, we show a stellar pulsation analysis and a Rossiter-McLaughlin analysis, and we present the Balmer

lines signatures we found in the transmission spectrum of KELT-9 b. In Sect. 4, we present a new model and retrieve fundamental atmospheric properties from the high-resolution data observations. In Sect. 5, we discuss the constraints we find on the KELT-9b atmospheric temperature, mass loss rate, and local thermodynamic properties. We present our conclusions in Sect. 6.

## 2. The HARPS-N dataset

We observed two transits of the ultra-hot Jupiter KELT-9 b around its bright A0V host star with the HARPS-N spectrograph (Cosentino et al. 2012) mounted on the 3.58 m Telescopio Nazionale Galileo (TNG) telescope in the Roque de los Muchachos observatory (La Palma, Spain)<fnref n="2" />.

We dedicated a full night of continuous observations to each transit of KELT-9 b on 31 July 2017 and 20 July 2018 (see Table 2). One fiber recorded the target flux (fiber A) and the other monitored the sky (fiber B). In total, we registered 95 spectra, each with an exposure time of 600 s. The time series cover each transit in full and cover baselines of ~2–3 h before and ~1.5 h after each transit, used to build accurate out-of-transit master spectra. The typical signal-to-noise ratio (S/N) calculated in the continuum around 400 and 600 nm is between 27–131 and 70–190, respectively (no spectra were discarded in the analysis; however, the spectra collected on the second night are affected by a failure in the ADC, making their blue part much noisier than in the first night). Finally, the CHESS.KING module (KINematic and Geometric properties of the system) handles all the system parameters (see Table 1) needed during the data analysis.

We identified 44 spectra (22 in each night) fully in transit (i.e., entirely exposed between the first and fourth contacts) and we considered the 51 remaining spectra out of transit. The HARPS-N raw frames are automatically reduced with the HARPS-N Data Reduction Software (DRS version 3.5). The DRS pipeline consists of an order-by-order optimal extraction of the two-dimensional echelle spectra. Then, for each exposure, the daily calibration set is used to flat-field, blaze correct, and wavelength-calibrate all the 69 spectral orders. Finally, the orders are merged into a one-dimensional spectrum and resampled (with flux conservation) between 387 and 690 nm with a 0.001 nm wavelength step (the wavelengths are given in the air and the reference frame is the solar system barycenter one, see Fig. 2). Our main analysis is conducted on the final HARPS-N merged spectra that have a spectral resolution of $\lambda/\Delta\lambda = 115\,000$ or $2.7\,\text{km s}^{-1}$. In the present study, we took advantage of the

<footnote n="1">¹ This data set is the same as in Hoeijmakers et al. (2018, 2019). Those previous studies focused on metal detections (species between atomic numbers 3 and 78), while our present study focuses only on the hydrogen lines, making it a complementary analysis. We note that no Helium (λ10830 Å) has been found in Nortmann et al. (2018).</footnote>

<footnote n="2">² Our observational strategy of recording several transits – to assess the reproducibility of the detected spectroscopic signatures – was made possible thanks to a Director’s Discretionary Time proposal asked for in summer 2017 (A35DDT4) and a standard OPTICON proposal in the general frame of our SPADES survey (2018A/038, PI: Ehrenreich).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Left panel: normalized out-of-transit master spectrum of KELT-9 in the 387–687 nm optical region. The stellar A0 spectral type makes the spectrum dominated by strong Balmer lines, though we notice the presence of numerous shallow metal lines, as shown in the inset. The Balmer lines ($H\alpha, \beta, \gamma, \delta, \epsilon, \zeta$), that are covered by the HARPS-N spectrograph wavelength range, are identified by the vertical dashed lines. Right panel: normalized out-of-transit master cross-correlation function (CCF) obtained with a dedicated A0 mask. The averaged stellar line has a measured projected velocity broadening of $v \sin i_* = 112.6$ km s$^{-1}$. The vertical dashed line represents the systemic velocity (−17.74 km s$^{-1}$).</figcaption></figure>

cross-correlation functions (CCFs) that are computed with the DRS. The CCFs are constructed by correlating each spectral order with a weighted stellar mask (a line list containing the main lines of the spectral type of interest, Pepe et al. 2002). Then, the data are further prepared with the CHESS. R00K module (spec-tRoscopic data tOOI Kit). For the KELT-9 observations, we use a customized A0 mask (as in Anderson et al. 2018) containing about 600 lines of typically Fe I, Fe II, and other metal species. The A0 mask is sensitive to both the stellar and planetary spectra. Since the stellar spectrum is rotationally broadened, we computed the CCFs over a 700 km s$^{-1}$ window around the systemic velocity and with an oversampling step of 0.082 km s$^{-1}$ (a tenth of the pixel size). Then, we conducted the CCF analysis on the sum of the CCF of all the 69 orders$^3$ (Fig. 2). The use of a custom A0 mask boost the stellar CCFs contrast by a factor of ~10 in comparison to the standard DRS reduction, allowing us to have a precision similar to the independent reduction presented in Borsa et al. (2019).

The correction of the telluric absorption lines in each spectrum is the purpose of the CHESS.KNIGHT module (Knowledge of the NIGHTly Telluric lines). In this study, we used the molecfit software, version 1.5.7, provided by ESO (Smette et al. 2015; Kausch et al. 2015). Following the method presented in Allart et al. (2017) and applied in Hoeijmakers et al. (2018, 2019); Seidel et al. (2019), synthetic telluric spectra were computed with the real-time atmospheric profiles measured above the TNG/HARPS-N location and then fitted to each observed spectrum prior to dividing them. The correction of the changing telluric lines is measured to be at the noise level, especially for the $(4v + \delta)$ water band in the $H\alpha$ region. We note that the $H\beta$ to $H\zeta$ regions are not affected by telluric lines.

### 3. A comprehensive analysis of spectroscopic transit data

Spectroscopic transit observations (at high-resolution) have been widely used to characterize exoplanetary systems. On the one

hand, they allow us to measure exoplanetary spin-orbit alignment, via the Rossiter-McLaughlin (RM) effect, with classical velocimetric measurements or with the tomographic and reloaded methods (Queloz et al. 2000; Collier Cameron et al. 2010; Cegla et al. 2016). On the other hand, the same data also allow us to measure exoplanetary transmission spectra (Wyttenbach et al. 2015). These two signals can be entangled between each other or with any non-homogeneous stellar surface effects, such as the center-to-limb variation (CLV) of spectral lines, or stellar activity (Czesla et al. 2015; Barnes et al. 2016; Cauley et al. 2017a). The method presented here aims to separate and to accurately measure the stellar and planetary spectra (see Fig. 3). Even though we are only taking some major effects into account, our method can be, in principle, adapted to correct for any stellar feature. An upcoming paper by Bourrier et al. (in prep.) will establish a rigorous formalism and present an accompanying pipeline (but see also Bourrier et al. 2020a; Ehrenreich et al. 2020, who already used the same method). Our method is implemented as in our previous work (Wyttenbach et al. 2015, 2017), using the CHESS modules BISHOP (Basic Insights on the Stellar pHotosphere Occulted by the Planet) and QUEEN (QUantitative Extraction of Exoplanet atmospheres in transmission), and applied to the KELT-9 HARPS-N data.

#### 3.1. On the extraction of a transmission spectrum

Because of the fluctuations of the Earth's atmospheric transmission, ground-based transit spectroscopy observations at high-resolution have to be self-normalized. The stellar spectra $f(\lambda, t)$ are divided by their respective values in a reference wavelength band $\langle\lambda_{\text{ref}}\rangle$ chosen in the continuum near the region of interest to get the normalized spectra: $\tilde{f}(\lambda, t) = f(\lambda, t)/f(\langle\lambda_{\text{ref}}\rangle, t)$. The master spectra in- and out-of-transit $F_{\text{in}}(\lambda) = \sum_{t \in \text{in}} f(\lambda, t)$ and $F_{\text{out}}(\lambda) = \sum_{t \in \text{out}} f(\lambda, t)$ are self-normalized in the same manner, giving $\tilde{F}_{\text{in}}(\lambda)$ and $\tilde{F}_{\text{out}}(\lambda)$. We now base our approach from the schematic shown in Fig. 3 and will proceed step by step.

We first treat an ideal case where there is no exoplanetary atmosphere, and we only consider the stellar rotation and the limb-darkening (and ignore other stellar effects). We are interested in extracting the Rossiter-McLaughlin effect. A possibility is to use the reloaded method (Cegla et al. 2016;

<footnote>$^3$ We discarded (by visual inspection) the order number 1, 2, 50–58, 61, 64–66, 68 and 69 because they contain no or very few stellar lines, or because they are strongly affected by interstellar medium (ISM) or telluric line absorptions.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Schematic of a spectroscopic transit observation of an exoplanetary system. The difference between an out-of-transit spectrum (left) and an in-transit spectrum (middle) gives us two essential quantities: the exoplanetary atmosphere transmission spectrum and the local stellar spectrum of the surface that is hidden behind the exoplanet (right). The temporal evolution of the stellar local spectrum will mostly depend on the projected stellar rotation (represented by the black arrow going through the star from pole to pole), the planetary orbit parameters (represented with the green arrows) and particularly the projected spin-orbit alignment angle. The temporal evolution of the exoplanet atmosphere spectrum will mostly depend on the atmosphere size, the planetary orbital projected velocity, and the local stellar spectrum that is transmitted through. Hence, it is necessary to separate and to measure accurately each component to retrieve the unknown parameters. The present figure is partially inspired by figures in Dravins et al. (2015, 2017); Dumusque et al. (2014); Cegla et al. (2016); Bourrier et al. (2017a).</figcaption></figure>

Bourrier et al. 2017a, 2018):

$$f_{\text{res}}(\lambda, t) = \tilde{F}_{\text{out}}(\lambda) - [1 - \delta(t)] \cdot \tilde{f}(\lambda, t), \quad (1)$$

where $1 - \delta(t)$ represents the value of the photometric light curve (white light curve), and $\delta(t)$ is the transit depth at time $t$. We compute the light curve with a Mandel & Agol (2002) model using the batman code (Kreidberg 2015) and the system parameters from Table 1. While the observed light-curve contains errors, these latter have a negligible effect on the normalization of the spectroscopic observations (Cegla et al. 2016; Bourrier et al. 2017a). In this ideal case, $f_{\text{res}}(\lambda, t_{\text{in}})$ represents the stellar flux from the planet-occulted photosphere regions, that is to say the local stellar spectrum $s_{\text{local}}(\lambda, t)$ emitted from the region behind the planet at time $t$. The measurement and modeling of $f_{\text{res}}(\lambda, t_{\text{in}})$ over time during the transit allow us to determine the sky-projected spin-orbit angle (Cegla et al. 2016). The Rossiter-McLaughlin measurements are mostly done with CCFs in the velocity domain (because of the boost in signal). However, we underline that this measurement is also possible in a single spectral line.

We continue by treating a second ideal case where an exoplanet with an atmosphere transits a star deprived of its rotation (and of any other stellar features, except the limb-darkening). Since the star is not rotating, we know that the shape of the local stellar spectra $s_{\text{local}}(\lambda, t)$ is not time-dependent and is equal to the out-of-transit spectrum $\tilde{F}_{\text{out}}(\lambda)$. In this case, the only difference between the shape of the in- and out-of-transit spectra is the transmission through the planet atmosphere, implying that the stellar flux is multiplied by the planet apparent size. Then, the formula proposed by Brown (2001) and implemented by considering the change in radial velocity of the planet in Wyttenbach et al. (2015), is applied in a straightforward manner:

$$-\mathcal{R}'(\lambda) = 1 - \frac{1}{N_{\text{in}}} \sum_{t \in \text{in}} \left[ \frac{[1 - \delta(t)] \cdot \tilde{f}(\lambda, t)}{\tilde{F}_{\text{out}}(\lambda)} \right]_{p} . \quad (2)$$

The $-\mathcal{R}'(\lambda)$ denotes the normalized spectrum ratio (averaged over $N_{\text{in}}$ in-transit spectra), or the averaged transmission spectrum, and is equal to $(R_p(\lambda)/R_\star)^2$. The $|p|$ expresses the spectrum

Doppler shift according to the planet radial velocity and places the observer in the planet rest frame. As for the Eq. (1), we also use the light curve normalization $1 - \delta(t)$ to put the in-transit spectra to their correct flux levels compared to the out-of-transit spectra. This correction is also necessary for the planetary atmosphere since the absorption depends on the relative stellar flux that is passing through the atmosphere (Pino et al. 2018). This means that atmospheric absorptions were slightly overestimated in previous works. The calculation of Eq. (2) remains precise enough for slowly rotating stars (Wyttenbach et al. 2017).

We then present a comprehensive case that considers the presence of the planet atmosphere and an inhomogeneous stellar surface. The local stellar spectra $s_{\text{local}}(\lambda, t)$ vary as a function of the stellar rotation and other effects (center-limb variation, convective blueshift, pulsations, gravity darkening, active regions, etc.). In transit configuration, we note that the light source transmitted through the planetary atmosphere is emitted from the region behind the atmosphere limb annulus<fnref n="4" />. We assume that the light source $s_{\text{local,atm}}(\lambda, t)$ coming from this annulus (with a surface ratio $\Sigma_{\text{atm}}$ which is several times $2R_pH/R_\star^2$) is equal in shape to the local stellar spectrum $s_{\text{local}}(\lambda, t)$ (that comes from the surface behind the planet), but differs in intensity from it by a factor $\Sigma_{\text{atm}}/\delta$, meaning that $s_{\text{local,atm}}/\Sigma_{\text{atm}} = s_{\text{local}}/\delta$. The light source $s_{\text{local,atm}}$ is absorbed by the exoplanet atmosphere and constitute $s_{\text{local,atm}} \cdot s_{\text{tr}}/\Sigma_{\text{atm}}$, $s_{\text{tr}}(\lambda, t)$ being the transmission spectrum we search to measure<fnref n="5" />. Owing to the absorption by the exoplanet atmosphere, the residuals $f_{\text{res}}(\lambda, t)$, computed with Eq. (1), are not equal to the local stellar spectrum $s_{\text{local}}(\lambda, t)$, but contain the sum of $s_{\text{local}}(\lambda, t)$ and of $s_{\text{tr}}(\lambda, t) \cdot s_{\text{local,atm}}/\Sigma_{\text{atm}}$. Hence, the spectrum ratio,

$$s_{\text{tr}}(\lambda, t) = \delta(t) \cdot \frac{f_{\text{res}}(\lambda, t) - s_{\text{local}}(\lambda, t)}{s_{\text{local}}(\lambda, t)} \quad (3)$$

is equal to the surface of the atmosphere limb (alone) times its wavelength dependent absorption, meaning that $\delta(t) + s_{\text{tr}}(\lambda, t) = (R_p(\lambda, t)/R_\star(\lambda))^2$ is the transmission spectrum. Taking into

<footnote n="4">⁴ We ignore any refraction in the exoplanetary atmosphere.</footnote>

<footnote n="5">⁵ We define the planetary absorption $A$ to be $A = s_{\text{tr}}/\Sigma_{\text{atm}}$.</footnote>
account that the shape of the out-of-transit flux $\tilde{F}_{\text{out}}(\lambda)$ is relative to the continuum reference band $\langle\lambda_{\text{ref}}\rangle$, we get the final transmission spectrum formula relative to $\langle\lambda_{\text{ref}}\rangle$:

$$ -\Re'(\lambda) = \frac{1}{N_{\text{in}}} \sum_{t \in \text{in}} \left. \frac{\delta(t) + s_{\text{tr}}(\lambda, t)}{\tilde{F}_{\text{out}}(\lambda)} \right|_{p} \\ = \frac{1}{N_{\text{in}}} \sum_{t \in \text{in}} \left. \delta(t) \cdot \frac{\tilde{F}_{\text{out}}(\lambda) - [1 - \delta(t)] \cdot \tilde{f}(\lambda, t)}{\tilde{F}_{\text{out}}(\lambda) \cdot s_{\text{local}}(\lambda, t)} \right|_{p} \quad (4) $$

Placing the observer in the planet rest frame is the last performed operation before the sum. We note that by setting $s_{\text{local}}(\lambda, t) = \delta(t) \forall \lambda$, we can simplify Eq. (4) to get Eq. (2).

We first acknowledge that the presence of the planet atmosphere, as anticipated by Snellen (2004), magnifies the RM effect since the amplitude of the RM effect is proportional to the transit depth (see also Dreizler et al. 2009). Additionally, at high-resolution, the planetary orbital velocity plays a major role. First, the change in shape of a stellar line due to the planet atmosphere happens when the planetary radial velocity is smaller than the full line broadening (which is dominated by the stellar $v \sin i_*$). Second, the line shape anomaly due to the RM effect is also altered by the planetary atmosphere when the two signals overlap in terms of velocity. Hence, the presence of a planet atmosphere biases the RM effect, in a way that is not a simple amplification, but rather by creating new anomalies that depend on the planetary radial velocity and of the spin-orbit angle. Despite Snellen (2004) and Czesla et al. (2012) analyzed specific stellar lines to find RM anomalies, these latter were not modeled with the effect of the planet orbital velocity. This effect has been precisely taken into account only recently by Borsa et al. (2019), who called it the “atmospheric Rossiter-McLaughlin” effect. This effect is particularly important for ultra-hot Jupiters who share many spectral lines with their host stars (see more details in Bourrier et al. 2020a; Ehrenreich et al. 2020).

Similarly to the RM effect being modified by the planet atmosphere, the planetary atmosphere absorption depends on the stellar flux that is passing through it. For example, when the local stellar and the atmospheric spectra overlap, the atmospheric absorption is lowered because less flux is coming from the star. This effect is automatically taken into account and corrected for in Eq. (4) thanks to the division by $s_{\text{local}}(\lambda, t)$.

This method relies on the knowledge of the local stellar profile $s_{\text{local}}(\lambda, t)$. One approach is to model $s_{\text{local}}(\lambda, t)$ by using stellar spectrum model at different limb-darkening angles (e.g., Yan et al. 2017; Yan & Henning 2018; Casasayas-Barris et al. 2019). Our approach is to empirically fit the observed local stellar profile (see Sect. 3.2.1). We emphasize the caveat, particularly during the overlap, of the uncertainty whether the local stellar profile remains constant.

## 3.2. Analysis of the local stellar photospheric CCFs

The analysis of the stellar CCFs is conducted on the CCFs constructed with the aforementioned A0 mask (similarly to Anderson et al. 2018; Borsa et al. 2019). The CCF residuals are computed with Eq. (1) in the velocity domain (see Fig. 4). They clearly show the local stellar CCFs (the CCF of the spectrum of the stellar surface that is hidden behind the exoplanet, see Sect. 3.2.2), the CCFs of some stellar pulsations (see Sect. 3.2.1) and the planetary atmosphere CCFs trace. The planet atmospheric trace is an independent confirmation of the metal (mainly Fe I, Fe II) detection of Hoeijmakers et al. (2018, 2019) and Borsa et al. (2019), and is not further discussed in this study. The main

<figure><img src="image_5.png" /><figcaption>Fig. 4. Map of the residual CCF flux during the transit of KELT-9 b. The observations of two nights are binned by 2 in phase. The 1st and 4th contacts are traced with the horizontal dashed white lines; the 0 km s⁻¹ and ±v sin i* are traced with vertical dashed white lines. Top panel: local stellar spectrum tracing the spin-orbit misalignment through the Rossiter-McLaughlin effect. The solid black line shows our best fit. The exoplanet transmission CCF, that has a slanted shape following the orbital velocity (dashed black line), traces the metal detections by Hoeijmakers et al. (2018, 2019). Bottom panel: residual CCF flux after the removal of the local stellar and of the planetary atmosphere spectra. The residuals are dominated by the stellar pulsations, whose strengths are 10 × lower than the local stellar signal (see the color scale). The pulsations are due to pressure-mode in a δ Scuti-type star.</figcaption></figure>

goal of this section is to focus on the pulsation signal, and then on the local stellar CCFs to have an updated value of the spin-orbit angle λ and to use this value as a prior in our Balmer lines analysis.

### 3.2.1. Spectroscopic detection of stellar pulsation

Prior to discussing the reloaded RM study, we first mention the spectroscopic detection of stellar pulsation. As all the different signals can be disentangled, we study the CCFs cleaned from the local stellar and the planetary atmosphere CCFs, but still keep the pulsation signals (Fig. 4). A rotationally broadened model (including convolution with the instrumental profile and with a local stellar CCFs width) are fitted to each CCF (Gray 2005; Anderson et al. 2018; Borsa et al. 2019; Dorval et al. 2020). We measure an averaged v sin i* value of 112.60 ± 0.04 km s⁻¹. We note that our CCFs (see Fig. 2) and v sin i* values are similar
<figure><img src="image_6.png" /></figure>

**Fig. 5.** Stellar pulsation signal seen from variation in the stellar $v \sin i_*$. The stellar CCF are cleaned from the local stellar and planetary atmosphere spectra, and fit with a stellar rotation profile. The pulsation period is $P_{\text{puls}} = 7.54 \pm 0.12$ h, and is consistent between the two nights (in dark blue and green, respectively) that are separated by one year. The averaged $v \sin i_*$ is reported for each night with the horizontal dashed lines, and has an averaged value of $v \sin i_* = 112.60 \pm 0.04$ km s$^{-1}$.

to those published in the results of Borsa et al. (2019). The instrumental profile (2.7 km s⁻¹) with the local stellar CCFs (see Sect. 3.2.2) have an averaged broadening of $\xi = 9.17 \pm 0.05$ km s⁻¹, making the total broadening ~113 km s⁻¹. The stellar local profiles are fully described by Gaussians. Indeed, the local stellar profile have a small broadening, despite the high disk integrated rotational broadening ($\xi \ll v \sin i_*$). This is because of the high observing cadence, the small impact parameter, and the inclined orbit (see Sect. 3.2.2). We measure a sinusoidal oscillation in the time series of the stellar $v \sin i_*$ around the averaged value. The oscillation signal is not in phase with the planetary transit. However, this signal has a period $P = 7.54 \pm 0.12$ h, in phase within the two nights separated by one year (Fig. 5). This period corresponds to the stellar pulsation signal seen in the KELT-9 TESS light curve of $P_{\text{puls}} = 7.5868 \pm 0.0009$ h (Wong et al. 2019), and is compatible with pressure modes (p-mode) in a δ Scuti-type star. This is therefore an independent spectroscopic detection of the KELT-9 δ Scuti-type pulsation. The first and second CCF moments are clearly described by a single sine function with a period $P_{\text{puls}}$ (the CCF radial velocities are the first CCF moment, and the $v \sin i_*$ are proportional to the second CCF moment), hinting to a sectoral mode oscillation (the azimuthal order of the mode *m* is equal to $\pm l$, the degree of the mode, see e.g., Balona 1986; Aerts et al. 1992, 2010). A complete analysis of the pulsation mode is out of the scope of this paper, but a follow-up study could bring complementary information to the spin-orbit analysis, in particular, for the stellar inclination.

### 3.2.2. “Reloaded” Rossiter-McLaughlin analysis

For the reloaded RM analysis, we treated the CCF residuals in a similar way to that of recent studies by Bourrier et al. (2020a); Ehrenreich et al. (2020). We first removed the signals from the stellar pulsations in the out-of-transit phases by fitting and removing as many Gaussians as there are of individual pulsation signals, allowing us to recompute a new out-of-transit master, from which we recompute again the CCF residuals. At this stage, we remove the stellar pulsation in-transit and the planet atmosphere signal with Gaussian fitting. By doing so, we

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

**Fig. 6.** Velocities of the local stellar surface hidden behind the planet during its transit. The observations of the two nights are shown in dark blue and green, respectively. The first to fourth contacts are traced with the vertical dashed lines. The overlap region between the local stellar and atmospheric signals (between phases −0.025 and 0) is shown in red and is not used for the fit. The spin-orbit misalignment is modeled with a solid body rotation (in orange), and with a differential rotation (in light blue). The differential rotation model is slightly preferred.

are left with the local stellar CCF alone ($s_{\text{local}}$; see e.g., Cegla et al. (2016); Wyttenbach et al. (2017); Bourrier et al. (2020a) for more descriptions of the local stellar CCF), from which we measure each local velocity (Fig. 6). We discard the velocities coming from the overlapping region with the atmospheric signal, because they are biased (Borsa et al. 2019). Local stellar CCFs are considered to be contaminated by the planet atmosphere if the local CCF full width at half maximum (FWHM) overlaps the planetary atmosphere CCF FWHM (10 data points between phase −0.025 and 0.000). We noticed a posteriori that discarding or keeping the overlapping region changes our results by about $1\sigma$. We also discard the outliers points at the egress, for which the CCFs (and velocities) are not detected following the S/N criteria of Cegla et al. (2016) and Bourrier et al. (2018).

We model the residual CCF velocities with the formalism of Cegla et al. (2016), that computes the brightness-weighted average rotational velocity behind the planet during each exposure. The model parameters are the spin-orbit angle, and the stellar velocity field. The latter could be described by a rigid rotation (with $v \sin i_*$ as a parameter) or a solar-like differential rotation law, which could be expected for A-type star (Reiners & Royer 2004; Balona & Abedigamba 2016). In this case, we have the three parameters $v_{\text{eq}}$, $i_*$, and $\alpha$ which describe a solar-like differential rotation law (Cegla et al. 2016). We use emcee (Foreman-Mackey et al. 2013) to explore the parameter space and to find the best models and their parameter confidence intervals (CI). For each of the two nights and the combination of the two nights, we use 200 walkers for 1500 steps and a burn-in size of 500 steps. We use the Bayesian information criterion (BIC) to compare between models. For the solid rotation case (BIC = 151.5), we found $v \sin i_* = 122.5 \pm 0.6$ km s⁻¹, and $\lambda = -85.65^\circ \pm 0.09^\circ$ (see Fig. A.1). The $v \sin i_*$ differs between the two nights and is also bigger than our spectral fit, and the one of Borsa et al. (2019). Despite this, the spin-orbit angle is consistent with Gaudi et al. (2017) and Borsa et al. (2019). The differential rotation case (BIC = 147.4) has a $\Delta$BIC = 4.1 in its favor, which is interpreted as “positive evidence” but not a confirmed detection ($\sigma \sim 2.2$). In this scenario,
we found $\alpha = -0.10^{+0.03}_{-0.04}$, $v_{eq} = 195^{+37}_{-23}$ km s$^{-1}$, $i_* = 143.2^{\circ+6.3^\circ}_{-6.0^\circ}$, and $\lambda = -85.01^{\circ \pm 0.23^{\circ}}$ (best solution, see Fig. A.2). We note that $\lambda$ is still compatible (at $3\sigma$) with the rigid rotation case, and that the corresponding $v \sin i_* = 116.9 \pm 1.8$ km s$^{-1}$ is more consistent with our spectral fit. Another differential solution, with a $\Delta$BIC = 3.8, was found with $\alpha = 0.09 \pm 0.03$, $v_{eq} = 155 \pm 18$ km s$^{-1}$, $i_* = 55.9^{\circ+11.6^\circ}_{-7.2^\circ}$, and $\lambda = -85.10^{\circ \pm 0.25^{\circ}}$, with a very similar local stellar trace to the previous solution, but with a higher $v \sin i_* = 127.7^{+2.7}_{-2.4}$ km s$^{-1}$ (see Fig. A.3). Our two differential rotation solutions give a true spin-orbit angle of $\psi = 89.6^{\circ \pm 0.7}$ and $\psi = 84.2^{\circ+0.4^\circ}_{-0.3^\circ}$, respectively. The best fit, shown in Figs. 4 and 6, is later used for our Balmer lines extraction.

The differential rotation solutions give equatorial velocities that are typical of A-type stars (Zorec & Royer 2012); furthermore an A-type star is likely to show differential rotation (Zorec et al. 2017). The different solutions have a solar-type differential rotation ($\alpha = 0.08$), and an “anti-solar”-type ($\alpha = -0.1$), respectively, meaning that in the latter case, that the high latitudes are rotating slightly faster than the equator. The equatorial velocities found imply that the stellar rotation periods are $P_{\text{rot}} = 0.7 \pm 0.1$ day (~17.1 h) and $P_{\text{rot}} = 0.56 \pm 0.09$ day (~13.6 h), respectively. We note that these periods do not correspond to the pulsation period of ~7.6 h. This means that the star is turning at ~33–56% of its theoretical critical Keplerian break-up velocity ($v_{\text{crit}} \sim 416$ km s$^{-1}$). With such an equatorial velocity, the angular rotation parameter $\omega = \Omega_{\text{rot}}/\Omega_{\text{crit}}$ is about 0.45–0.75 (Maeder 2009). This implies that gravity darkening effects become significant and measurable, for example by TESS or high quality photometry (as already done by Wong et al. 2019; Barnes et al. 2011; Ahlers et al. 2015, 2020). For example, with $\omega = 0.45–0.75$, the star becomes oblate and the equatorial radius becomes up to ~4–12% bigger than the polar radius (Georgy et al. 2014). In a gravity darkened star, the polar regions are hotter and more luminous than the equatorial regions. Because of the stellar inclination measured in the differential rotation scenario, we are mainly measuring light from the polar region of the star. Hence, at $\omega = 0.45–0.75$, the stellar luminosity is overestimated up to ~5–15%, and the effective temperature up to ~1–3%, depending on the stellar inclination toward the observer (Georgy et al. 2014). This may be the source of the differences in retrieved parameters between Gaudí et al. (2017) and Borsa et al. (2019). Thus, a follow-up study with high-resolution spectroscopy and precise photometry, which take care of gravity darkening effects, must be undertaken, and will allow us to refine our solution and to choose between different degenerate geometries (e.g., Zhou et al. 2019; Dorval et al. 2020; Ahlers et al. 2020). In particular, the gravity darkening provides a third possibility to measure the stellar inclination alongside to the Rossiter-McLaughlin effect and the stellar pulsation mode analysis. Finally, because of the polar orbit, the small impact parameter and the stellar inclination, the planet may be transiting near the pole during the transit (Wong et al. 2019). Such a transit configuration is favorable to measure star-planet interaction because of anisotropic stellar flux and winds (“gravity-darkened seasons”, Ahlers et al. 2020).

### 3.3. The transmission spectrum of KELT-9 b in the Balmer lines regions

To extract the planetary atmosphere information from the spectra, we follow the procedure that leads to Eq. (4). We assume the system parameters from Table 1. In particular, for each Balmer line we use the same white light curve from Gaudí et al. (2017)

to normalize the spectra because it is the only one available<fnref n="6" />. To correct the transmission spectra from the local stellar photospheric signal, we must know $s_{\text{local}}(\lambda, t)$ for every Balmer line. We assume that the local stellar signal in the Balmer lines follows the CCF one in terms of velocity. For each Balmer line, and for each phase, we fit a Gaussian to the local stellar spectrum, with the velocity centers following our best reloaded RM model from Sect. 3.2. The use of the differential rotation model instead of the rigid rotation model does not change our final line extraction. Indeed, the two models are similar and the velocity errors from the CCFs are much smaller than the one present in a single Balmer line. Thus, we can fix the velocities to the best reloaded RM model in our fit. The Gaussian offset values follow the light curve $\delta(t)$. We therefore only fit for the contrast and the FWHM values. In a second step, we averaged all the contrast and FWHM values, and fitted again a Gaussian (for each line, for each phase) with starting positions coming from the averaged values, and by not allowing us the fit to exceed $\pm 1\sigma$ value of the averaged contrast and FWHM. Whenever the local stellar spectrum signal is too weak for the Gaussian fit to converge (or to be significant), we assume the averaged contrast and FWHM values. We note that the stellar pulsations, despite being expected in the Balmer lines, are too weak to be detectable given the S/N of a single line. Thus, the pulsation contaminations in the Balmer lines are at the noise level. Also, we removed low-frequency trends in the residual continuum by fitting polynomials (Seidel et al. 2019). The knowledge of $s_{\text{local}}(\lambda, t)$, and the $K_p$ from Hoeijmakers et al. (2019) allows us to compute the weighted averaged transmission spectra with Eq. (4)<fnref n="7" />. We discard the spectra where the atmospheric signal overlaps with the local stellar signal ($-0.027 \le \text{Phase} \le 0.0015$), because we cannot fit $s_{\text{local}}(\lambda, t)$ with a Gaussian for these phases. Figure B.1 shows a representative example of the effect of the local stellar photospheric signal correction (RM effect correction). We measure the RM effect to have the same relative amplitude for the different Balmer lines and to be corrected down to the noise level. We compute the errors on the absorption lines by considering the systematic noise, empirically measured by the standard deviation into the continuum. We also consider the stellar line contrasts that decrease the S/N in the absorption region (see Table 2). This is particularly relevant for the blue part of the detector. Each line center takes the systemic velocity into account (e.g., 6562.408 Å for $H\alpha$); the reference passbands for the normalization of every spectrum are centered 6 Å on each side of the line center, and are 6 Å wide (e.g., [6553.408:6559.408] \cup [6565.408:6571.408] for $H\alpha$). For the different Balmer lines, our main results are:

$H\alpha$. The $H\alpha$ line is detected with a high confidence of 24.3$\sigma$, which makes the planetary trail visible for each phase (Fig. 7). Figure B.1 shows that the magnitude of the RM effect on the $H\alpha$ line is about 0.3% in the blue part of the line (around 30% of the line core contrast), and about $\le 0.05–0.1\%$ in the other parts of the line. The RM effect correction is made down to the noise level. We note that our results for $H\alpha$ are in conflict with previous literature: by ~3$\sigma$ with Yan & Henning (who used the CARMENES spectrograph 2018), by ~4.5$\sigma$ with Cauley et al. (2019, PEPSI) and by ~1$\sigma$ with Turner et al. (2020, CARMENES). The difference could come from

<footnote n="6">⁶ This has little impact since the expected relative change of the line contrasts would be less than 0.5% (from the prediction of Fig. 14 of Yan et al. 2019), that is smaller than our individual error bars.</footnote>

<footnote n="7">⁷ In practice, Eq. (4) is used with a weighted mean that takes the changing data quality into account. The weights $w$ at time $t$ are the inverse of the squared flux uncertainties, $w(t) = 1/\sigma(t)^2$.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Left: detection of the planetary Hα line (significance of 24.3σ). Top left panel: map of the residual spectra ±200 km s⁻¹ around the Hα line in the stellar rest frame (Eq. (1)). The observations of two nights are binned by 2 in phase. The planetary atmosphere Hα absorption follows the orbital velocity (solid black line). The local Hα stellar line follows the reloaded RM model (dashed black line). The first and fourth contacts are traced with the horizontal dashed white lines, the 0 km s⁻¹ and ±ν sin i* are traced with vertical dashed white lines. Bottom left panel: weighted averaged Hα transmission spectrum (corrected from the local stellar photospheric signal, Eq. (4)) in the planet rest frame (light gray), also binned by 15× in black circles (the continuum is set to the white light curve transit depth δ = 0.00677). The vertical black dashed line shows the line center taking the systemic velocity into account. We show the Gaussian fit (blue) with contrast of 0.91 ± 0.04% and FWHM of 44.3 ± 1.8 km s⁻¹. Right: same as for Hα, but for the Hβ line (detection significance of 16.1σ, see Table 3).</figcaption></figure>

our new extraction method, the different resolutions, or different
normalizations. However, we note that Casasayas-Barris et al.
(2019) reported Hα detections in KELT-20 b/MASCARA-2 b,
that are different for the HARPS-N and CARMENES instru-
ments, despite using the same data reduction. It is thus likely
that there are instrumental systematics that are not yet taken
into account (detector non-linearity, blaze and bias removal,
spectrum interpolations, moonlight pollution, etc.). A follow-up
study exploring such effects must be undertaken.

Hβ. The Hβ line and its trail are also clearly detected (16.1σ, Fig. 7). Our Hβ detection also disagrees with the result of Cauley et al. (2019, by ~2σ) and is subject to the same caution. Another explanation for the Hα and Hβ differences between studies could be intrinsic variation in the atmosphere. We also note that the line absorption is stronger in the second part of the transit than in the first part, hence the average measured line could arise from different weights given to each spectrum (due to the observation conditions).

Hγ. For Hγ the atmospheric trail is barely visible, since
the total spectrum 7.5σ detection is distributed to ~1σ in
each spectrum (Fig. 8). Hβ and Hγ have only been reported
in the atmosphere of three exoplanets: HD 189733 b, KELT-
9 b, and KELT-20 b/MASCARA-2 b (Cauley et al. 2016, 2019;
Casasayas-Barris et al. 2019).

Hδ. The Hδ detection in the KELT-9 b atmosphere (4.6σ) is the first of its kind to our knowledge (Fig. 8). Despite its significance, the line is probably slightly biased because the local stellar signal becomes less visible and is more difficult to correct for. This may be due to the lack of signal in the bluer part of the detector.

Hε. The Hε absorption (2.9σ) is just below the significance
level of 3σ. Therefore it should be treated as a hint of detection
(Fig. 9). Additional data are needed to better constrain the Hδ
and Hε line shapes.

Hζ. Hζ is not detected in our data set (Fig. 9). The most
probable reason is due to the lack of flux in the very blue part of
the HARPS-N wavelength range. Our upper limit value is hinting
that the Hζ absorption is smaller than all the other Balmer lines,
as expected.

Table 3 presents a summary of our spectral analysis of the Balmer series. From the contrasts of the KELT-9 b Balmer Hα, Hβ, and Hγ lines, we know that neutral hydrogen is present at altitudes between ~1.2 R<sub>P</sub> and ~1.5 R<sub>P</sub> (see Figs. 7 and 8). Since this is close to the planetary Roche radius (1.95 ± 0.2 R<sub>P</sub>), the strong Balmer lines absorptions are a hint of atmospheric escape (Yan & Henning 2018). We note that no line is significantly blueshifted or redshifted in the line of sight. With an average shift of 0.2 ± 0.6 km s<sup>-1</sup>, we do not infer any strong winds in the upper atmosphere. This is compatible with the high circulation drag expected in very hot atmospheres (Koll & Komacek 2018, but also see Wong et al. (2019) who found a relatively high recirculation efficiency). Finally, we note that the Hα, Hβ, and Hγ lines have a similar FWHM, with an averaged value of 44.9 ± 1.4 km s<sup>-1</sup>. For the rest of our analysis, we investigate the strong and consistent Hα, Hβ, and Hγ signals that contain most of the information.

**4. Interpretation of the Balmer lines of KELT-9 b**

The measured FWHMs of the KELT-9 b Balmer Hα, Hβ, and Hγ
lines cannot be explained only by thermal broadening. Indeed,
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Same as for Fig. 7, but for the Hγ line (**left**) and the Hδ line (**right**). The Hγ and Hδ lines are detected with a significance of 7.5σ and 4.6σ, respectively (see Table 3).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Same as for Fig. 7, but for the Hε line (**left**), and the Hζ line (**right**). The Hε and Hζ lines are not statistically detected (see Table 3).</figcaption></figure>

for hydrogen, such a FWHM would require a thermal broadening with a temperature around 40 000 K, far above any theoretical prediction for a thermosphere temperature (e.g., Salz et al. 2016). Since the temperature cannot be measured directly, a rough estimate of the atmospheric properties of KELT-9 b from the Balmer series can be made following Huang et al. (2017, Sect. 2) as done in previous studies (Yan & Henning 2018; Turner et al. 2020). The estimation of Huang et al. (2017) assumes an isothermal hydrostatic atmosphere with a mean molecular weight $\mu = 1.3$ and a temperature $T = 5000$ K. For KELT-9 b, a better assumption is an atmosphere dominated by ionized hydrogen ($\mu = 0.66$) and a thermosphere temperature of $T = 10000$ K (Fossati et al. 2018; García Muñoz & Schneider 2019). With these numbers, we find that to explain the FWHM<fnref n="8" /> of our Hα detection, one

needs an optically thick layer ($\tau \sim 40$) of excited hydrogen H I(2) with a number density of $n \sim 8.2 \times 10^3$ cm⁻³ along a typical line path of $L = 2 \sqrt{2R_p N_{sc} H} \sim 13.8 \times 10^9$ cm, where $N_{sc} \sim 6$ is the number of absorbing scale-heights (Madhusudhan et al. 2014; Madhusudhan & Redfield 2015; Huang et al. 2017). This means that at an altitude of $N_{sc}H \sim 0.3 R_p$ ($H \sim 6000$ km), which is about half of the line core contrast, the Hα line core have $\tau \gg 1$. Since the line is optically thick in most of the layers around that altitude, the line FWHM is larger than the thermal width (Huang et al. 2017). Also, this is consistent with our observations since the Hα core is observed at $\sim 1.5 R_p$, where $\tau$ becomes close to unity (Lecavelier Des Etangs et al. 2008; de Wit & Seager 2013; Heng et al. 2015; Heng & Kitzmann 2017).

Another estimate of the atmospheric conditions can be done with the Lecavelier Des Etangs et al. (2008) formalism, which again assumes an isothermal hydrostatic atmosphere. Indeed, one can measure the atmospheric scale-height H thanks to the

<footnote n="8">⁸ To compare our FWHM, we removed the additional width contribution by other effects (see Sects. 4.1 and 4.2.1).</footnote>
<table><caption>**Table 3.** Summary of the hydrogen Balmer line detections.</caption><thead><tr><th></th><th>C<br>[%]</th><th>v<br>[km s<sup>-1</sup>]</th><th>FWHM<br>[km s<sup>-1</sup>]</th><th>δ<sup>(a)</sup><br>[%]</th></tr></thead><tbody><tr><td>Hα</td><td>0.91 ± 0.04</td><td>2.1 ± 0.8</td><td>44.3 ± 1.8</td><td>0.68 ± 0.03</td></tr><tr><td>Hβ</td><td>0.68 ± 0.04</td><td>-1.7 ± 1.3</td><td>46.3 ± 2.5</td><td>0.45 ± 0.03</td></tr><tr><td>Hγ</td><td>0.44 ± 0.07</td><td>0.9 ± 2.9</td><td>43.5 ± 5.5</td><td>0.31 ± 0.04</td></tr><tr><td>Hδ</td><td>0.61 ± 0.13</td><td>-4.6 ± 2.2</td><td>25.6 ± 5.7</td><td>0.42 ± 0.09</td></tr><tr><td>Hε</td><td>0.39 ± 0.14</td><td>-15.0 ± 5.0</td><td>29.5 ± 10.0</td><td>0.26 ± 0.09<sup>(b)</sup></td></tr><tr><td>Hζ</td><td>≤0.50</td><td>—</td><td>—</td><td>≤0.34<sup>(b)</sup></td></tr></tbody></table>

**Notes.** (a) The absorption depth is averaged over the line full width (see Wyttenbach et al. 2017), takes systematic noise into account, and is used to compute the line detection significance. (b) $H\epsilon$ is not statistically detected. For $H\zeta$ we show the $3\sigma$ upper limit assuming a Gaussian line shape with a fixed $v = 0$ km s$^{-1}$ and FWHM = 30 km s$^{-1}$.

linear correlation between the line contrasts (or the altitudes $z$) and the hydrogen Balmer line oscillator strengths (see the $\log_{10} gf$ values in Sect. 4.1). Knowing that $\Delta z = H\Delta \ln(gf) = H\ln(10)\Delta \log_{10}(gf)$, the slope of the $z$-versus-$\ln(gf)$ linear fit gives us $H$. We computed that $H = 9600 \pm 1400$ km around the average altitude of absorption (1.3 $R_P$). When the gravity is evaluated at 1.3 $R_P$ (~0.59 $g_{surf}$), we have that $T/\mu \approx 13\,500 \pm 2000$ [K/u]. Assuming a $\mu$ between 0.66 and 1.26 (atmospheric composition dominated by a mixture of ionized and neutral hydrogen and helium), the upper atmosphere temperature is constrained between 7500 and 19 500 K. The lower half temperature range is more probable because the thermal ionization of hydrogen decreases $\mu^9$.

Finally, our line contrast measurements and the aforementioned assumptions lead to the same conclusions as in Yan & Henning (2018) of a Jeans mass loss rate of ~$10^{12}$ g s$^{-1}$, although hydrodynamic escape is more likely the mass loss mechanism in KELT-9 b (Fossati et al. 2018).

While this kind of estimation is useful for getting a qualitative idea of the planet aeronomic conditions, the data quality calls for a quantitative estimate. A few forward self-consistent models of the hydrogen $H\alpha$ line exist (Christie et al. 2013; Huang et al. 2017; Allan & Vidotto 2019; García Muñoz & Schneider 2019). They explore different physical questions and have different atmospheric structures, but none of them are linked to any retrieval algorithms. Retrieval approaches (on similar absorption lines) have been made by Cauley et al. (2019); Fisher & Heng (2019); Seidel et al. (2020), but without the mass loss rate being a parameter. Thus, we decided to follow our own retrieval approach with a model adapted to our data set.

## 4.1. The PAWN model

We model the planet atmosphere in transmission with the customized model PAWN (Parker Winds and Saha-BoltzmanN atmospheric model) which is included in our CHESS framework as a tool for modeling exoplanetary thermospheric lines in transmission spectroscopy. Our purpose is to constrain fundamental parameters of the thermosphere regions such as the temperature, the mass loss rate and the hydrogen number density (local thermodynamical condition), if possible. Another goal is to have a

sufficiently fast model to feed a retrieval algorithm. Hence, our modeling holds several assumptions.

We model a 1D planetary atmosphere, with a special interest in the upper atmosphere part (thermosphere). A logarithmic grid with $n_r$ radius $r$ layers is chosen between 1 and $r_{up} R_P$. For the hydrogen Balmer lines, $r_{up} = 2$ ensures the convergence of the models. We model the atmospheric structure $\rho(r)$ with a hydrostatic or a hydrodynamic structure. For the hydrostatic solution, we have the base density $\rho_0(r=1)$ as a free parameter, and consider that the pressure scale-height $H = k_B T/\mu g$ can change with the altitude according to $T(r)$, $\mu(r)$, $g(r)$ (see e.g., the $\pi\eta$ code in Ehrenreich et al. 2006; Pino et al. 2018; Seidel et al. 2020). In the hydrostatic case, we can estimate a lower limit to the mass loss rate, by computing the Jeans escape rate at the planetary Roche radius, $r_{Roche}$ (Eggleton 1983; Ridden-Harper et al. 2016), or at the exobase radius, $r_{exobase}$, if below (Lecavelier des Etangs et al. 2004; Tian et al. 2005; Salz et al. 2016; Heng 2017):

$$ \dot{M}_{\text{Jeans}} = 4\pi r_{\text{Roche}}^2 \mu(r_{\text{Roche}}) F_{\text{Jeans}}(r_{\text{Roche}}), \quad (5) $$

where $F_{\text{Jeans}}$ is the Jeans escape flux (Shizgal & Arkos 1996; Hunten et al. 1989; Chamberlain & Hunten 1987). For the hydrodynamical solution, we use the isothermal Parker wind solution following the formalism of Heng (2017)<fnref n="10" />, Oklopčić & Hirata (2018), and Lampón et al. (2020). Our free parameter for the atmospheric structure becomes the atmospheric mass loss rate:

$$ \dot{M}_{\text{PW}} = 4\pi r^2 \rho(r) v(r), \quad (6) $$

where $\rho(r)$ and $v(r)$ solve the Parker equations (Parker 1958; Lamers & Cassinelli 1999) from the sonic point ($r_s$, $v_s$, $\rho_s$, $\mu_s$) down to the lowest radius in the grid. We note that $\dot{M}_{\text{PW}}$ is constant through the atmosphere. Since the Parker wind solution is isothermal, we keep an isothermal profile for both the hydrostatic and hydrodynamic structure with the temperature $T$ as a free parameter. We note that the assumption of an isothermal profile for the upper thermosphere (typically between 1.1 and 2–3 $R_P$ for hot Jupiters) is a good approximation of aeronomy models (e.g., Lammer et al. 2003; Lecavelier des Etangs et al. 2004; Yelle 2004; García Muñoz 2007; Murray-Clay et al. 2009; Koskinen et al. 2013; Salz et al. 2016). Indeed, above the thermobase where the temperature changes dramatically, the temperature encounters a maximum and is varying much slower through several orders of magnitude in pressure than in the lower regions.

For the atmospheric constituent structure, we use a pre-computed chemical equilibrium table obtained with the chemical equilibrium code of Mollière et al. (2017). The chemical equilibrium (of a variety of neutral and singly ionized atomic, molecular, and condensed species) is computed, with a Gibbs energy minimization, into a grid of pressure ($10^{-15}-10^2$ bar) and temperature (1000–20 000 K). The mean molecular weight $\mu(\rho(r), T(r))$, the mass fraction $m_t(r)$, and volume mixing ratio $n_t(r)$ are given by the chemical grid, where $t$ is a chosen element and ionization of interest (H, He, Na, K, Ca, II, etc.). For example, in the hydrogen case, we know the ratio of ionized to neutral hydrogen $n_{HII}/n_{H\alpha}$ from reading and interpolating the chemical grid with the correct $T$ and $\rho$ from our atmospheric structure. When the Parker wind solution is in use, we link it with

$$ M^2 - \ln(M^2) - 1 = 4 \ln(x) + 4(1/x - 1), $$

<footnote n="10">¹⁰ Equation (13.14) of Heng (2017, Chap. 13) should be corrected as follows:</footnote>

<footnote>where the Mach number $M = v/c_s$, and $x = r/r_s$, with $c_s$, $r_s$ being the sound speed, and the sonic point radius, respectively.</footnote>

<footnote n="9">⁹ A composition dominated by molecular hydrogen and helium ($\mu = 2.34$) is probably excluded by an unphysical $T \ge 30\,000$ K.</footnote>
the chemical grid at the sonic point down to the lowest radius. We note that for hydrogen, the Gibbs energy minimization is equivalent to computing the Saha law with an electron number density $n_e$ coming from the ionization of the hydrogen and from all the other atomic species in the grid. By default, the chemical grid is always used, but it is also possible to assume a constant (free parameter) $\mu$ through the atmosphere. In that case the $n_{H1}/n_{H1}$ can be constant or can follow the Saha equation with $n_e$ accounting only for the electron coming from the hydrogen ionization. In any case, we are ignoring photo-ionization effects that can be important in the upper atmosphere.

Once we know the vertical distribution of our element of interest, we still need to know what is the number density of the electronic level $\ell(n)$ from which the transitions we observed are arising. In the case of H1, we are observing the Balmer series in absorption, which arise from photons being absorbed by the $n=2$ electronic state, that is to say the H1(2) state. If the atmosphere is in local thermodynamic equilibrium (LTE), the distribution of the electronic states follow the Boltzmann equation (for hydrogen):

$$ \frac{n_{H1(n)}}{n_{H1}} = \frac{2n^2}{\nu(T)} e^{E_1/k_B T \cdot (n^2-1)/n^2}, \qquad (7) $$

where $E_1 = 13.606$ eV, and $\nu(T)$ is the hydrogen partition function. The Boltzmann equilibrium is assumed by default in our code. However, because we are probing high in the thermosphere, the collisions between particles become less numerous and the electronic state distribution can depart from the Boltzmann distribution, the atmosphere is in non-LTE or NLTE (Barman et al. 2002; Fortney et al. 2003; Christie et al. 2013; Huang et al. 2017; Oklopčić & Hirata 2018; Fisher & Heng 2019; Allan & Vidotto 2019; García Muñoz & Schneider 2019; Lampón et al. 2020). For the NLTE case, we define the departure coefficients (e.g., Mihalas 1970; Salem & Brocklehurst 1979; Barman et al. 2002; Draine 2011; Hubeny & Mihalas 2014) as:

$$ \beta_{H1(n)} = \frac{n_{H1(n)}^{\text{NLTE}}}{n_{H1(n)}^{\text{LTE}}}. \qquad (8) $$

Then, we can let the different $\beta_{H1(n)}$ be free parameters (or equivalently all the $n_{H1(n)}$ are free parameters). Another possibility is to force all the $n_{H1(n)}$ to depart by the same amount from LTE, meaning that only one constant free parameter $\beta$ is used for all $n_{H1(n)}$. We note that a constant $n_{H1(n)}/n_{H1}$ through the thermosphere is an assumption justified by Huang et al. (2017) and García Muñoz & Schneider (2019).

When an element and its transitions are chosen (in our case the Hydrogen Balmer series), the opacity $\kappa$ [cm$^2$ g$^{-1}$] are computed in each layer of the atmosphere according to the line list of Kurucz (1979, 1992), to the temperature, and to the number density $n_{H1(n)}$. We do assume that the line intensity arises from photon absorption and from stimulated emission (Sharp & Burrows 2007). For the hydrogen Balmer transition from electronic level i to j, we have:

$$ \kappa_{ij} = \frac{\pi e^2}{m_e c} \frac{g_i f_{ij}}{m_H n_{H1}} \left( \frac{n_{H1(i)}}{g_i} - \frac{n_{H1(j)}}{g_j} \right) \Phi_V(v_{ij}), \qquad (9) $$

where $e$ is the elementary charge, $m_e$ is the electron mass, and $c$ the speed of light in vacuum. The $f_{ij}$, $g_i$, $g_j$ are the oscillator strength, and the statistical weights of the electronic levels due to their degeneracies, respectively. For the hydrogen Balmer series, we have $\log_{10} g_2 f_{2j} = 0.71, -0.02$, and $-0.447$ for $j=3$,

4, and 5, respectively. A Voigt profile $\Phi_V$, centered on the line transition center $v_{ij}$ (or $\lambda_{ij}$), is assumed for the line shape. The Lorentzian wing component considers the natural broadening (for the hydrogen Balmer series, the Einstein A-coefficients are $\log_{10}(A_{j2} [s^{-1}]) = 8.76, 8.78, 8.79$ for $j=3, 4$, and 5, respectively). There is no pressure broadening. The Gaussian core component takes the thermal broadening into account. We can also add velocity components due to rotation, circulation, winds in the atmosphere (see below). For this exploratory study, we decided not to consider the opacity arising from $H^-$ (the main source of background continuum opacity, Kitzmann et al. 2018), as a rough estimate shows that its opacity ($\sim 10^2$ cm$^2$ g$^{-1}$) stays smaller than the one from each Balmer line ($10^3-10^8$ cm$^2$ g$^{-1}$ in the Gaussian core). Hence, the modeled absorptions are generally robust for the line cores that we observe (Hoeijmakers et al. 2019). The spectral resolution of our model is chosen to be 0.001 Å ($\lambda/\Delta\lambda \sim 10^6$), and the wavelength coverage is 20 Å around each line.

Next, we perform the radiative transfer in a grazing geometry following the method of Ehrenreich et al. (2006); Mollière et al. (2019). For each impact parameter, the optical depth $\tau = \sum \kappa \tilde{m}_i$, where $\tilde{m}_i$ is the column mass of the element $\ell$ in [g cm$^{-2}$], is computed through the line of sight (LOS). The transmission function $\mathcal{T} = e^{-\tau}$ is computed according to a pure absorption in each LOS. The transmission spectrum is computed assuming a spherically symmetric atmosphere, weighted along the impact parameter surface annuli, and normalized to the stellar surface.

Within the previous two steps, our line profile can be broadened according to different atmospheric patterns (see Brogi et al. 2016; Seidel et al. 2020). First, we assume a solid body planetary rotation (perpendicular to the orbital plane). Since the planet rotation is axisymmetric, this breaks the spherical symmetry assumed before. In that case, we divide each hemisphere of the planet in a grid of $m_\theta$ circular sector. The sector grid is evenly separated in $\cos\theta$, where $\theta$ is the angle going from the equator to the pole back to the equator (anti-clockwise). If the rotation is included, our code becomes essentially a 2D atmospheric code. The planet rotation period $P_P$ is the parameter that defines the solid rotation. $P_P$ can be a free parameter, or more simply be fixed to the planetary orbital period by assuming that the planet is tidally locked (which is our default choice). It is interesting to note that in this configuration, each LOS has a fixed velocity along the cord, which is given by:

$$ v_{\text{rot}}^{\text{LOS}}(r, \theta) = \frac{2\pi r}{P_P} \cos(\theta). \qquad (10) $$

In the case of KELT-9b, the planetary rotation have a broadening effect of $\sim 6.1$ km s$^{-1}$ at $1 R_P$ for each hemisphere ($\sim 12.2$ km s$^{-1}$ in total).

Second, when the atmosphere is hydrodynamical (Parker winds), a spherically symmetric expansion of the atmosphere occurs with a velocity $v(r)$. This expansion changes the projected velocity of the different components of each LOS, thus creating another broadening source. This Parker wind broadening can be considered in our model, but not simultaneously to the rotation (because the symmetry is not the same, the structure would be 3- and the computational cost would become too large). The LOS velocity is given by:

$$ v_{\text{PW}}^{\text{LOS}}(r_i, r_j) = \frac{v(r_j)}{r_j} \sqrt{r_j^2 - r_i^2}. \qquad (11) $$

In practice, for KELT-9b, for most of the mass loss rates (except for unphysical large rates), the upward velocity $v(r)$ is
smaller than 0.01–1 km s⁻¹ in the probed layers (i.e., the sonic point r_s, where the velocity becomes several km s⁻¹ is far above the Roche radius). This leads to a negligible (much smaller than the instrumental resolution) parker wind broadening (though see Seidel et al. 2020, for other possible scenarios). Thus, for the rest of our study, we decided to ignore this effect in favor of the planetary rotation. We note that these two aspects have been previously explored in Seidel et al. (2020). They were able to keep a 1D atmospheric structure at the price of several approximations in the atmospheric velocity patterns (pseudo-2D/3D).

In the next step, our line profile is modified according to the observational set-up. First, the orbital blurring due to the planet changing velocity during an exposure time is computed with a box-shape convolution (of a width of ~6.8 km s⁻¹ for the KELT-9 b orbit and for the 600 s exposure time). The spectrum is then convolved with the instrument profile, which is assumed to be a Gaussian (with a width of ~2.7 km s⁻¹ for HARPS-N). Then, the spectrum is binned down to match the data sampling. Finally, the differential normalization (due to the loss of the absolute depth in ground-based observations, see Wyttenbach et al. 2015, 2017; Pino et al. 2018) is performed according to the aforementioned reference passbands of each line. The models, as for the data, are normalized to the white light curve transit depth of KELT-9b ($\delta$ = 0.00677). The models are now in a form which is comparable with the data. Our PAWN model is finally linked with the Markov chain Monte Carlo (MCMC) posterior sampling code emcee (Foreman-Mackey et al. 2013). Typically, for one line, one model without rotation is computed in $\leq 0.1$ s ($n_r$ = 100), while with rotation, one model takes ~0.5 s ($n_r$ = 50, $m_\theta$ = 14). We note that a model with $n_r$ = 50 and $m_\theta$ = 14 is $\leq 1\%$ equal to a model with $n_r$ = 250 and $m_\theta$ = 28.

In the next section, we apply our PAWN model to the Balmer lines detection in the KELT-9 b atmosphere.

## 4.2. Interpretation of the Balmer series in the KELT-9 b atmosphere with the PAWN model

In this section, we investigate the Hα, Hβ, and Hγ signals based on their line shape. All lines are fitted together in a range of 20 Å (6000 data points, notably to cover the normalization passbands). Our retrieval approach uses our PAWN model with the emcee sampler to explore the parameter space and to find the best models and their parameter confidence intervals (CI). For each different set of models, we use ten walkers for each dimension (parameter) during 2500 steps and with a burn-in size of 500 steps. We use the bayesian information criterion (BIC) to compare between models.

For every model, each line center position is a free parameter (three parameters for Hα, Hβ, and Hγ). The line center positions are sensitive to circulation winds in the upper atmosphere (Snellen et al. 2010; Brogi et al. 2016; Seidel et al. 2020). Uniform priors between ±20 km s⁻¹ are used to ensure that our prior cover the full posteriors. We noticed a posteriori that these parameters are not correlated with any of the other parameters. Also, we found that our MCMC results are confirming the line shifts and errors quoted in Table 3, thus we do not show the results for these parameters. The second set of parameters defines the atmospheric structure with the temperature T, and the base density $\rho_0(r = 1)$ or the mass loss rate $\dot{M}$ depending on if the structure is hydrostatic or hydrodynamic, respectively (for these latter, we use their $\log_{10}(\rho_0)$ values). We use a uniform prior on T between 1000 and 20 000 K. We use a log-uniform prior between −18 and −4 for the $\log_{10}(\rho_0)$ [g cm⁻³]). We use a

log-uniform prior between 8 and 17 for the $\log_{10}(\dot{M})$ [g s⁻¹]). In our base models, we use the Boltzmann equilibrium and there are no other parameters. We first explore this set-up, compare the hydrostatic and hydrodynamic models, before exploring some NLTE models.

### 4.2.1. LTE models with a hydrostatic or hydrodynamic structure

We first asses how our simplest models are able to fit the observed lines. These two set of models are hydrostatic and hydrodynamic, respectively. They have two free parameters (apart of the line centers), and assume a Boltzmann equilibrium. For the hydrostatic model, we find a solution with a thermosphere temperature of T = 13 180$^{+650}_{-580}$ K, and a base density of $\log_{10}(\rho_0)$ [g cm⁻³]) = −10.74 ± 0.05 (BIC = 5016.1, Fig. 10). The corresponding Jeans mass loss rate at the Roche radius is computed to be $\dot{M}_{\text{Jeans}} \sim 6.6 \times 10^{11}$ g s⁻¹. Interestingly, our hydrodynamical model, which has a BIC = 5017.2, is found to have a similar thermosphere temperature of T = 13 200$^{+800}_{-720}$ K, while the mass loss rate is tending towards to the higher value of $\log_{10}(\dot{M}_{\text{PW}})$ [g s⁻¹]) = 12.8 ± 0.3 (i.e. $\dot{M}_{\text{PW}} \approx 5.9 \times 10^{12}$ g s⁻¹, Fig. 11). Both models are comparable in term of statistics and none can be preferred; this is because the density profiles are not different enough at the altitudes probed. The mass loss rate is in both cases around $10^{12}$ g s⁻¹, confirming the result of Yan & Henning (2018). Also, we can safely conclude that the thermosphere temperature, if in LTE, must be near 13 000 K (in line with our first estimate using the Lecavelier Des Etangs et al. 2008 equation). Despite the high temperature, the line FWHMs ($44.9 \pm 1.4$ km s⁻¹) are not dominated by thermal broadening. Indeed, for T = 13 200 K, the hydrogen thermal broadening is ~25 km s⁻¹, while the other source of broadening are ~16, ~7, and ~2.7 km s⁻¹ for the rotational (at r = 1.3 $R_p$), orbital blurring, and instrumental broadening, respectively. This means that the first source of broadening (~33 km s⁻¹) comes from optically thick absorptions through most of the lower layers (Huang et al. 2017).

We notice that the line contrasts and FWHMs are both enhanced when the temperature or the mass loss rate are increasing. For the temperature, this is caused by the atmosphere being more extended, and by thermal broadening. For the mass loss rate, this is caused by larger densities present at higher altitudes that increase the LOS optical depth, thus boosting the broadening and the contrast. Therefore, the shape of the correlation in the hydrodynamical case (Fig. 11) is explained by the following mechanism: when the temperature increases, the number density of the hydrogen second electronic level builds up. However, as for the stellar case, the total number density of neutral hydrogen decreases because of ionization. This last effect eventually dominates at high temperature. So, to keep a constant absorption corresponding to our observations, a higher density is needed. This is achieved by increasing the mass loss rate (Eq. (6)). Thus, we get the positive correlation between T and $\dot{M}_{\text{PW}}$, which shape is given by the Boltzmann equation.

### 4.2.2. NLTE models with a hydrodynamic structure

We explore the NLTE case for the hydrodynamical atmosphere case. We first run MCMC chains by letting all the electronic level number density be free parameters. The prior on each $\log_{10}(\beta_{H1(n)})n_{H1(n)}^{\text{LTE}}/n_{H1}$ was set to be a log-uniform distribution between −40 and 0. With this test, we found that the model hardly converges, and that the parameters are not
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Correlation diagrams in the case of a hydrostatic atmosphere in the Boltzmann equilibrium.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Correlation diagrams in the case of a hydrodynamic atmosphere in the Boltzmann equilibrium.</figcaption></figure>

well-constrained (see Fig. C.1). The explanation lies in two different families of solutions. One family is distributed around the Boltzmann equilibrium, while the other is away from it, and has a lower temperature and a higher mass loss rate. We therefore further explore these two families.

For the first family, we limit our prior ranges on the $log_{10}(\beta_{H1(n)}n_{H1(n)}^{LTE}/n_{H1})$ between -9, -12, -16, -20 and 0, for n=2, 3, 4, 5, respectively. The results are shown in Fig. C.2, and the model results in a BIC = 5031.4. The increased number of parameters clearly exclude this model ($\Delta$BIC = -14.2 < -5), despite a better best fit. We still notice that the retrieved temperature lies between 12 600 and 16 600 K, while the mass loss rate lies between $10^{12}$ and $10^{16}$ g s$^{-1}$, encompassing the LTE result.

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Correlation diagrams in the case of a hydrodynamic atmosphere with a departure $\beta$ from LTE.</figcaption></figure>

For the second family of solutions in NLTE, we notice that it is mainly dominated by the departure of the second electronic level. We explore this family, by reducing the number of parameters for the number densities. Instead of having all the $log_{10}(\beta_{H1(n)}n_{H1(n)}^{LTE}/n_{H1})$ be free parameters, we decide to use a unique departure coefficient $\beta$ for all the electronic levels, that is $\beta_{H1(n)} = \beta$ for all H I(n). For this last MCMC run, a log-uniform prior for $log_{10}(\beta)$ was set between -10 and 5. The result is shown in Fig. 12. We confirm that this set-up explores a parameter space close to the second aforementioned family, where the temperature is lower and the mass loss rate is higher. In this case, the best fit is again better, and have a BIC = 5014.8. The $\Delta$BIC = 2.4 from the hydrodynamical-Boltzmann model, meaning a small positive evidence against the LTE model (and only a smaller $\Delta$BIC = 1.3 from the hydrostatic Boltzmann model). In any case, the $\Delta$BIC is not sufficient to conclude in favor of any NLTE models (see the comparison in Fig 13). Nevertheless it is interesting to note a few particularities in this scenario. First, the $\beta$ parameter is strongly correlated with the temperature and mass loss rate, giving a conservative idea of the uncertainties on these two parameters. Second, the mass loss rate is $10^{14.3_{-1.4}^{+1.2}}$ g s$^{-1}$ (between $8 \times 10^{12}$ and $3 \times 10^{15}$ g s$^{-1}$). Such a mass loss rate could have dramatic consequences (see discussion in Sect. 5). Third, the departure from LTE is $10^{-6.1_{-1.1}^{+3.3}}$ (in the $10^{-8.7}-10^{-2.2}$ range at $3\sigma$), making the LTE excluded in this framework. Interestingly, our best value is far (but still compatible at ~2.5$\sigma$) from the theoretical work by García Muñoz & Schneider (2019) who predicted a departure from Boltzmann to be at the level of ~$10^{-3}$. Our finding could hint to another possible departure from LTE than only the one from the Boltzmann equilibrium (see discussion in Sect. 5).

## 5. Discussion

With our PAWN framework, we are able to retrieve a good fit to the data (Fig. 13), and show that the thermosphere temperature (in LTE) is near 13 000 K and the mass loss rate is around $10^{12}$ g s$^{-1}$. If the hydrodynamic expansion of the thermosphere
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Hα, Hβ, and Hγ transmission spectra (light gray, and binned by 10 × in black circles) in the planet rest frame and in velocity space (the continuums are set to the white light transit depth δ = 0.00677). The best PAWN models, in the case of a hydrodynamically expanding thermosphere, and in the Boltzmann equilibrium, are shown in blue. The retrieved temperature is T = 13 200<sup>+800</sup><sub>-720</sub> K, and the mass loss rate is M&#775;<sub>PW</sub> = 10<sup>12.8±0.3</sup> g s<sup>-1</sup>. The best PAWN models, in the case of a hydrodynamically expanding thermosphere, with a LTE departure, are shown in dashed red. The retrieved temperature is T = 9600<sup>+1210</sup><sub>-1180</sub> K, the mass loss rate is M&#775;<sub>PW</sub> = 10<sup>14.3±1.4</sup> g s<sup>-1</sup>, and the common LTE departure coefficient is β = 10<sup>-6.1±1.1</sup>.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Summary of the mass loss rates and their errors found in the LTE case compared to the expected mass loss rates from the theoretical studies of Fossati et al. (2018) and García Muñoz & Schneider (2019). Our finding shows that the KELT-9 b thermosphere probably undergoes a “Balmer-driven” expansion, that leads to a high escape rate. In this case, the total mass loss through the planetary life time (500 My) represents ~0.05 M<sub>Jup</sub>, which is ~2% of the planet mass.</figcaption></figure>

is a valid assumption, then the thermosphere temperature is con-
strained to be T = 13 200<sup>+800</sup><sub>-720</sub> K, and the mass loss rate to be
M&#775;<sub>PW</sub> = 10<sup>12.8±0.3</sup> g s<sup>-1</sup>. This mass loss rate is higher than the
energy-limited escape rate of 10<sup>10</sup>–10<sup>11</sup> g s<sup>-1</sup> found by Fossati
et al. (2018, see Fig. 14). It is important to note that this study
only considers the X-FUV and NUV irradiation from the star,
which may be small in the case of a A0-type star without
chromosphere. However, García Muñoz & Schneider (2019) pre-
sented a new escape mechanism that applies to ultra-hot Jupiters
orbiting early-type stars. This mechanism shows that it is impor-
tant to also consider the NUV and optical Balmer continuum and
Balmer series irradiation from the star. In this scenario, the ther-
mosphere is first heated by Lyman photons (Lyman continuum

and especially Lyα), which increases the number density of
excited hydrogen in the high atmosphere (Lyman irradiation and
de-excitation put the population levels out of Boltzmann equi-
librium), thus increasing the Balmer opacity. At a certain point,
the excited hydrogen atoms are able to absorb significantly much
more flux from the stellar Balmer photons, which constitutes an
enormous source of energy in early-type stars. This increases
even more the thermosphere temperature, which reaches an
equilibrium temperature of ~12 000–15 000 K (3000 K more
than with only the X-FUV-NUV flux). In this case the mass
loss rate increases to 10<sup>12</sup>–10<sup>13</sup> g s<sup>-1</sup>, the atmosphere under-
goes a “Balmer-driven” escape, named so because the Balmer
energy deposition dominates over all other. This escape rate is
well in line with our findings (Fig. 14). Knowing that the age
of KELT-9 is ~300–500 Myr (Gaudi et al. 2017; Borsa et al.
2019), a mass loss of this rate over the planetary life time rep-
resents a total mass loss of ~0.03–0.05 M<sub>Jup</sub> (~1.5–2% of the
planet mass).

In the thermosphere, the excited hydrogen atoms are
expected to be in departure from LTE. The main reason being
the lower frequency of collisions at high altitudes, making the
hydrogen energy levels depart from the Boltzmann equilibrium,
typically with a β factor of 10⁻³ (García Muñoz & Schneider
2019). Despite the data being too low in S/N to show clear evi-
dences for a departure from the LTE, our NLTE model hints at
a lower temperature of T = 9600 ± 1200 K, and a higher mass
loss rate of log₁₀(Ṁ<sub>PW</sub> [g s⁻¹]) = 14.3⁺¹.²⁻¹.⁴. If true, the planet may
have already lost 1–2 Ṁ<sub>Jup</sub>, and would undergo a catastrophic
evaporation that would make the planet barely survive the next
~500–700 My, when its host star will enter the red-giant phase.
Such a mass loss rate is not yet physically motivated (though see
García Muñoz & Schneider 2019). One possibility would be to
explore the real effect of the stellar gravity and magnetic field.
Indeed, such effects could push down the sonic point to much
lower altitude, forcing the density and velocity to increase at
high altitude, hence increasing the mass loss rate (Murray-Clay
et al. 2009). Another possibility would be to investigate the effect
of the high rotation of KELT-9, which is also relatively young.
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Example of PAWN models of the Ca II H line detection by Yan et al. (2019). Our model is compatible with their detection only when assuming some LTE departure. For the Ca II, the departure is likely explained by further thermal and/or photo-ionization not yet taken in our code into account. Additional data gathering, and modeling will be necessary to confirm our interpretation.</figcaption></figure>

Indeed, both effects potentially increase the mass loss rate (Allan & Vidotto 2019).

In our NLTE model, we found that the departure coefficient $\beta$ is on the order of $10^{-6}$, while theory predicts $10^{-3}$ (only for Boltzmann departure). Thus, our much smaller $\beta$ coefficient could hint to some other possible LTE departure. Indeed, because $\beta$ is eventually normalized to the total amount of neutral hydrogen, our finding may signify that there is less neutral hydrogen than in the chemical (Saha) equilibrium by a factor of $10^3$. This implies that the LTE departure may also be caused by photo-ionization (Saha equilibrium departure, e.g., Fortney et al. 2003; Sing et al. 2008; Vidal-Madjar et al. 2011; Lavvas et al. 2014). In this case, the ratio of neutral to ionized hydrogen would be $10^{-3}$ lower. The other atoms would probably also undergo photo-ionization, making the ionized species even more dominant in the high atmosphere, increasing the importance of ionospheric phenomena (Helling et al. 2019).

Moreover, our NLTE scenario may be able to explain the strong absorption by ionized atomic species, especially the Ca II and Fe II, by Hoeijmakers et al. (2018, 2019); Cauley et al. (2019); Yan et al. (2019); Turner et al. (2020). It is interesting that Yan et al. (2019) found a different line ratio between the calcium doublet and triplet for the two ultra-hot Jupiter, WASP-33 b and KELT-9 b. Hence, this type of data is also sensing the local thermodynamical conditions, in complement to the Balmer lines. We can predict absorptions of the calcium doublet Ca II H&K, and infrared triplet (IRT) lines, by taking the two main atmospheric structures of our best fit models of the Balmer lines. To model the calcium lines, we assume the calcium abundance to be solar ($\chi_{Ca} = 6.34$, Asplund et al. 2009), and $r_{up} = 10$ (for good convergence). In our code, the calcium population is dominated by ionized Ca II at the temperatures and pressures of our models (see also Kitzmann et al. 2018; Lothringer et al. 2018; Yan et al. 2019, although for very high temperatures, the Ca III and further ionization levels are becoming important). In this case, our model is compatible with the detections of Yan et al. (2019); Turner et al. (2020) only if we assume a departure from LTE on the order of $\sim 10^{-2}$, or $\sim 10^{-4.5}$ for our two main structures, respectively (Fig. 15). Furthermore, the Ca II opacities are much larger than the Balmer opacities (by $\sim 100$ times), but

the measured absorptions have about the same contrasts, which hints at the presence of some Ca II NLTE. Because the calcium H&K lines are arising from the ground state of the singly ionized calcium, Boltzmann equilibrium departure is not likely to explain this, because the lack of collision is mostly populating the ground state level (the smaller calcium IRT absorption may be explained by non-Boltzmann distribution since this triplet arises from the first excited level). The $\beta$ departure in the Ca II case may also be due to a non-solar metallicity, but it is unlikely that the planet atmosphere metallicity is far below the solar one (KELT-9 has a solar metallicity, and planets are expected to generally have a higher metallicity than their host star). Thus, the departure from LTE may be better explained by a second ionization of calcium (through thermal and photo-ionization), which would lower the singularly ionized calcium population. The departure level for Ca II is in line with the extra departure found for hydrogen (we note that the ionization energy of the H I and Ca II have similar value of 13.6 eV and 11.9 eV, respectively). More investigations have to be done to explore if the hypothesis of chemical equilibrium in the thermosphere of (ultra) hot Jupiter is still justified.

One caveat of our framework is that our set of assumptions may be too strong when interpreting, together or not, the hydrogen and calcium absorptions. For example, the assumption of an isothermal atmosphere is probably too strong because temperature gradients and temperature inversion are predicted (Kitzmann et al. 2018; Lothringer et al. 2018; Mansfield et al. 2020). Actually, such a temperature inversion has been detected for the first time by Pino et al. (2020) in the KELT-9 b atmosphere by observing neutral iron lines in emission on the day side. Their detection shows that the temperature inversion happens from the 1 bar region up to the $10^{-3}-10^{-5}$ bar region. However, Hoeijmakers et al. (2019) infers that the maximum probable pressure in transmission is in the mbar regime, due to the $H^-$ opacity in the optical. The consequence for our model is that defining the lower boundary for the radius integration at $1 R_p$ (the measured white light radius) is justified because this radius is set by the, effectively gray, $H^-$ opacity. At no wavelength should the high-resolution observations probe smaller radii. However, because we assume that the atmosphere is isothermal, at high temperatures this can lead to pressures at $1 R_p$ that are much smaller than mbar. Forcing a fit to continue to mbar pressures increases the line core to continuum amplitude unrealistically. Such a continuation to larger pressures at fixed temperature is also not correct because the atmosphere must transition from the high thermosphere temperature to the planetary equilibrium temperature. Accounting for the decrease in atmospheric scale-height, the line core to continuum amplitude would then stay compatible with our current model. Another effect, not considered in our model, is a potential difference in the vertical distribution of the individual atmospheric gases that are determined by their respective masses, whose effect is to give birth to a heterosphere on top of the homosphere (Shizgal & Arkos 1996). For example, the presence of a heterosphere may also explain the weak calcium absorption, since these atoms would stay in the lower thermosphere. However, this effect of atomic diffusion may be counterbalanced by the hydrodynamic outflow that is able to lift heavy atomic species, that has been observed in exospheres (Vidal-Madjar et al. 2004). A last caveat is the modeling of the exosphere itself. Because the Roche radius of KELT-9 b is low ($1.95 \pm 0.2 R_p$), it is expected that hydrogen and metals escape easily. The kinetic broadening in the exosphere would change the shape of the lines, making them asymmetric (Ehrenreich et al. 2015; Bourrier et al. 2015,
2016; Allart et al. 2019). This effect has not yet been seen for the KELT-9 b Balmer lines. Thus, additional data gathering, and modeling will be necessary to get a coherent image of the KELT-9 b atmosphere consistent with all the observations.

## 6. Conclusion

We analyzed two optical high-resolution transit observations of KELT-9 b gathered with HARPS-N, with our CHESS framework. We found a spin-orbit misalignment of $\lambda = -85.01^{\circ} \pm 0.23^{\circ}$, and positive evidence in favor of a differential rotation model. We also confirmed the presence of stellar pulsations with a period of $P_{\text{puls}} = 7.54 \pm 0.12$ h. We use the reloaded RM analysis to correct the Balmer lines from the local stellar photospheric spectrum effect in transmission spectroscopy. We detected the $H\alpha$, $H\beta$, $H\gamma$, and $H\delta$ absorptions coming from the planet's upper atmosphere, with high confidence (between 24.3 and $4.6\sigma$) during both nights. The $H\alpha$, $H\beta$, $H\gamma$ lines probe altitudes up to ~1.2 to ~1.5 planetary radii; they do not show any significant velocity shifts, and share a large FWHM of $44.9 \pm 1.4$ km s$^{-1}$ that is dominated by optically thick absorptions. We presented a new model, CHESS.PAWN, to interpret the absorption in the Balmer series in exoplanet atmospheres that is linked with a retrieval algorithm. We constrained the thermosphere temperature to be $T = 13\ 200^{+800}_{-720}$ K, and the mass loss rate to be $\dot{M}_{\text{PW}} = 10^{12.8 \pm 0.3}$ g s$^{-1}$, under the assumptions of a hydrodynamic expansion and of LTE. This escape rate is fully explained by the “Balmer-driven” escape mechanism proposed by García Muñoz & Schneider (2019). NLTE effects are not yet significantly detected with this data set.

Our work shows that with the infrared metastable Helium triplet ($\lambda10830$ Å, Seager & Sasselov 2000; Oklopčić & Hirata 2018; Spake et al. 2018; Allart et al. 2018; Nortmann et al. 2018), the Balmer series is an alternative new tool for probing exoplanet evaporation, alongside with classical UV observations (Vidal-Madjar et al. 2003). Furthermore, these two probes seems to be complementary since the helium lines are mainly sensitive to G and K-type exoplanetary systems (Oklopčić 2019), while the Balmer lines are more sensitive to early-type (A-type to G-type) systems (Jensen et al. 2012, 2018; Yan & Henning 2018; Casasayas-Barris et al. 2018, 2019; Cabot et al. 2020; Chen et al. 2020). Finally, such observations will benefit from the ESPRESSO spectrograph (Pepe et al. 2013, 2014). On the one hand, ESPRESSO will increase the S/N of the absorption lines (Ehrenreich et al. 2020; Chen et al. 2020). On the other hand, in combination with near-infrared spectrographs (CARMENES, SPIRou, NIRPS), it will be possible to detect more absorption lines arising from different electronic levels of a given atomic species. For example, the combination of the Balmer and Paschen (or Brackett) series for the hydrogen, of the infrared triplet, and of the D3 line for the helium, or of the H&K doublet and of the infrared triplet for the calcium, will allow us to infer the local thermodynamic state. This will eventually allow us to put constraints on several chemical and physical aeronomic processes.

**Acknowledgements.** A.W. acknowledges the financial support of the SNSF by grant number P2GEP2_178191 and P400P2_186765. This work has been carried out within the frame of the National Centre for Competence in Research “PlanetS” supported by the Swiss National Science Foundation (SNSF). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement number 724427; project Four Aces). I.S. and P.M acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 694513 and 832428; ExoplanetBio). L.P. acknowledges that the research leading to these results has

received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 679633; Exo-Atmos). We thank the anonymous referee for the reading and the pertinent comments and questions. A.W. thanks J. Alonso-Floriano for discussions about transmission spectroscopy, C. Aerts for discussions about stellar pulsations, and A. Vidotto, A. García-Muñoz, M. Bergemann, and Y. Miguel for discussions about NLTE in hot Jupiter atmospheres.

## References

Aerts, C., de Pauw, M., & Waelkens, C. 1992, A&A, 266, 294

Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, *Asteroseismology* (Berlin: Springer)

Ahlers, J. P., Barnes, J. W., & Barnes, R. 2015, ApJ, 814, 67

Ahlers, J. P., Kruse, E., Colón, K. D., et al. 2020, ApJ, 888, 63

Allan, A., & Vidotto, A. A. 2019, MNRAS, 490, 3760

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

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

Allart, R., Bourrier, V., Lovis, C., et al. 2019, A&A, 623, A58

Alonso-Floriano, F. J., Snellen, I. A. G., Czesla, S., et al. 2019, A&A, 629, A110

Anderson, D. R., Temple, L. Y., Nielsen, L. D., et al. 2018, MNRAS, submitted [arXiv:1809.04897]

Arcangeli, J., Désert, J.-M., Line, M. R., et al. 2018, ApJ, 855, L30

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

Balona, L. A. 1986, MNRAS, 219, 111

Balona, L. A., & Abedigamba, O. P. 2016, MNRAS, 461, 497

Barman, T. S., Hauschildt, P. H., Schweitzer, A., et al. 2002, ApJ, 569, L51

Barnes, J. R., Haswell, C. A., Staab, D., & Anglada-Escudé, G. 2016, MNRAS, 462, 1012

Barnes, J. W., Linscott, E., & Shporer, A. 2011, ApJS, 197, 10

Bell, T. J., & Cowan, N. B. 2018, ApJ, 857, L20

Borsa, F., Rainer, M., Bonomo, A. S., et al. 2019, A&A, 631, A34

Bourrier, V., Ehrenreich, D., & Lecavelier des Etangs, A. 2015, A&A, 582, A65

Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., Tanaka, Y. A., & Vidotto, A. A. 2016, A&A, 591, A121

Bourrier, V., Cegla, H. M., Lovis, C., & Wyttenbach, A. 2017a, A&A, 599, A33

Bourrier, V., Ehrenreich, D., Allart, R., et al. 2017b, A&A, 602, A106

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

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

Bourrier, V., Kitzmann, D., Kuntzer, T., et al. 2020b, A&A, 637, A36

Brogi, M., de Kok, R. J., Albrecht, S., et al. 2016, ApJ, 817, 106

Brown, T. M. 2001, ApJ, 553, 1006

Cabot, S. H. C., Madhusudhan, N., Welbanks, L., Piette, A., & Gandhi, S. 2020, MNRAS, 494, 363

Casasayas-Barris, N., Palle, E., Nowak, G., et al. 2017, A&A, 608, A135

Casasayas-Barris, N., Pallé, E., Yan, F., et al. 2018, A&A, 616, A151

Casasayas-Barris, N., Pallé, E., Yan, F., et al. 2019, A&A, 628, A9

Cauley, P. W., Redfield, S., Jensen, A. G., et al. 2015, ApJ, 810, 13

Cauley, P. W., Redfield, S., Jensen, A. G., & Barman, T. 2016, AJ, 152, 20

Cauley, P. W., Redfield, S., & Jensen, A. G. 2017a, AJ, 153, 217

Cauley, P. W., Redfield, S., & Jensen, A. G. 2017b, AJ, 153, 81

Cauley, P. W., Redfield, S., & Jensen, A. G. 2017c, AJ, 153, 185

Cauley, P. W., Shkolnik, E. L., Ilyin, I., et al. 2019, AJ, 157, 69

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

Chamberlain, J. W., & Hunten, D. M. 1987, *Theory of Planetary Atmospheres*. An Introduction to Their Physics and Chemistry (Cambridge: Academic Press), 36

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

Chen, G., Casasayas-Barris, N., Pallé, E., et al. 2020, A&A, 635, A171

Christie, D., Arras, P., & Li, Z.-Y. 2013, ApJ, 772, 144

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

Cosentino, R., Lovis, C., Pepe, F., et al. 2012, Proc. SPIE, 8446, 84461V

Czesla, S., Schröter, S., Wolter, U., et al. 2012, A&A, 539, A150

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

de Wit, J., & Seager, S. 2013, Science, 342, 1473

Dorval, P., Talens, G. J. J., Otten, G. P. P. L., et al. 2020, A&A, 635, A60

Draine, B. T. 2011, *Physics of the Interstellar and Intergalactic Medium* (Princeton: Princeton University Press)

Dravins, D., Ludwig, H.-G., Dahlen, E., & Pazira, H. 2015, in Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, 18th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, 18, 853

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

Dreizler, S., Reiners, A., Homeier, D., & Noll, M. 2009, A&A, 499, 615

Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132
Eggleton, P. P. 1983, ApJ, 268, 368

Ehrenreich, D., Bourrier, V., Wheatley, P. J., et al. 2015, Nature, 522, 459

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

Ehrenreich, D., Tinetti, G., Lecavelier Des Etangs, A., Vidal-Madjar, A., & Selsis, F. 2006, A&A, 448, 379

Fisher, C., & Heng, K. 2019, ApJ, 881, 25

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

Fortney, J. J., Sudarsky, D., Hubeny, I., et al. 2003, ApJ, 589, 615

Fossati, L., Koskinen, T., Lothringer, J. D., et al. 2018, ApJ, 868, L30

García Muñoz, A. 2007, Planet. Space Sci., 55, 1426

García Muñoz, A., & Schneider, P. C. 2019, ApJ, 884, L43

Gaudi, B. S., Stassun, K. G., Collins, K. A., et al. 2017, Nature, 546, 514

Georgy, C., Granada, A., Ekström, S., et al. 2014, A&A, 566, A21

Gibson, N. P., Merritt, S., Nugroho, S. K., et al. 2020, MNRAS, 493, 2215

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

Helling, C., Gourbin, P., Woitke, P., & Parmentier, V. 2019, A&A, 626, A133

Heng, K. 2017, Exoplanetary Atmospheres: Theoretical Concepts and Foundations (Princeton: Princeton University Press)

Heng, K., & Kitzmann, D. 2017, MNRAS, 470, 2972

Heng, K., Wyttenbach, A., Lavie, B., et al. 2015, ApJ, 803, L9

Hoeijmakers, H. J., Ehrenreich, D., Heng, K., et al. 2018, Nature, 560, 453

Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165

Huang, C., Arras, P., Christie, D., & Li, Z.-Y. 2017, ApJ, 851, 150

Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton: Princeton University Press)

Hunten, D. M., Donahue, T. M., Walker, J. C. G., & Kasting, J. F. 1989, Escape of Atmospheres and Loss of Water., eds. S. K. Atreya, J. B. Pollack, & M. S. Matthews (Tucson: University of Arizona Press), 386

Jensen, A. G., Redfield, S., Endl, M., et al. 2012, ApJ, 751, 86

Jensen, A. G., Cauley, P. W., Redfield, S., Cochran, W. D., & Endl, M. 2018, AJ, 156, 154

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

Keles, E., Mallonn, M., von Essen, C., et al. 2019, MNRAS, 489, L37

Khalafinejad, S., Salz, M., Cubillos, P. E., et al. 2018, A&A, 618, A98

Khalafinejad, S., von Essen, C., Hoeijmakers, H. J., et al. 2017, A&A, 598, A131

Kitzmann, D., Heng, K., Rimmer, P. B., et al. 2018, ApJ, 863, 183

Koll, D. D. B., & Komacek, T. D. 2018, ApJ, 853, 133

Koskinen, T. T., Harris, M. J., Yelle, R. V., & Lavvas, P. 2013, Icarus, 226, 1678

Kreidberg, L. 2015, PASP, 127, 1161

Kreidberg, L., Line, M. R., Parmentier, V., et al. 2018, AJ, 156, 17

Kurucz, R. L. 1979, ApJS, 40, 1

Kurucz, R. L. 1992, IAU Symp., 149, 225

Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999, Introduction to Stellar Winds (Cambridge: Cambridge University Press)

Lammer, H., Selsis, F., Ribas, I., et al. 2003, ApJ, 598, L121

Lampón, M., López-Puertas, M., Lara, L. M., et al. 2020, A&A, 636, A13

Lavvas, P., Koskinen, T., & Yelle, R. V. 2014, ApJ, 796, 15

Lecavelier des Etangs, A., Vidal-Madjar, A., McConnell, J. C., & Hébrard, G. 2004, A&A, 418, L1

Lecavelier Des Etangs, A., Pont, F., Vidal-Madjar, A., & Sing, D. 2008, A&A, 481, L83

Lothringer, J. D., & Barman, T. 2019, ApJ, 876, 69

Lothringer, J. D., Barman, T., & Koskinen, T. 2018, ApJ, 866, 27

Madhusudhan, N., Knutson, H., Fortney, J. J., & Barman, T. 2014, Protostars and Planets VI (Tucson: University of Arizona Press), 739

Madhusudhan, N., & Redfield, S. 2015, Int. J. Astrobiol., 14, 177

Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars (Berlin: Springer)

Mandel, K., & Agol, E. 2002, ApJ, 580, L171

Mansfield, M., Bean, J. L., Line, M. R., et al. 2018, AJ, 156, 10

Mansfield, M., Bean, J. L., Stevenson, K. B., et al. 2020, ApJ, 888, L15

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

Miguel, Y., & Kaltenegger, L. 2014, ApJ, 780, 166

Mihalas, D. 1970, Stellar Atmospheres (Wales: W.H. Freeman & Co Ltd;)

Mollière, P., van Boekel, R., Bouwman, J., et al. 2017, A&A, 600, A10

Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67

Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23

Nortmann, L., Pallé, E., Salz, M., et al. 2018, Science, 362, 1388

Nugroho, S. K., Gibson, N. P., de Mooij, E. J. W., et al. 2020, MNRAS, submitted [arXiv:2003.04856]

Oklopčić, A., & Hirata, C. M. 2018, ApJ, 855, L11

Oklopčić, A. 2019, ApJ, 881, 133

Parker, E. N. 1958, ApJ, 128, 664

Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, A&A, 617, A110

Pepe, F., Mayor, M., Galland, F., et al. 2002, A&A, 388, 632

Pepe, F., Cristiani, S., Rebolo, R., et al. 2013, The Messenger, 153, 6

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

Pino, L., Ehrenreich, D., Wyttenbach, A., et al. 2018, A&A, 612, A53

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

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

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

Reiners, A., & Royer, F. 2004, A&A, 415, 325

Ridden-Harper, A. R., Snellen, I. A. G., Keller, C. U., et al. 2016, A&A, 593, A129

Salem, M., & Brocklehurst, M. 1979, ApJS, 39, 633

Salz, M., Czesla, S., Schneider, P. C., & Schmitt, J. H. M. M. 2016, A&A, 586, A75

Salz, M., Czesla, S., Schneider, P. C., et al. 2018, A&A, 620, A97

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

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

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

Sharp, C. M., & Burrows, A. 2007, ApJS, 168, 140

Shizgal, B. D., & Arkos, G. G. 1996, Rev. Geophys., 34, 483

Sing, D. K., Vidal-Madjar, A., Lecavelier des Etangs, A., et al. 2008, ApJ, 686, 667

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

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

Snellen, I. A. G. 2004, MNRAS, 353, L1

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

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

Spake, J. J., Sing, D. K., Evans, T. M., et al. 2018, Nature, 557, 68

Stangret, M., Casasayas-Barris, N., Pallé, E., et al. 2020, A&A, 638, A26

Tian, F., Toon, O. B., Pavlov, A. A., & De Sterck, H. 2005, ApJ, 621, 1049

Turner, J. D., de Mooij, E. J. W., Jayawardhana, R., et al. 2020, ApJ, 888, L13

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

Vidal-Madjar, A., Désert, J.-M., Lecavelier des Etangs, A., et al. 2004, ApJ, 604, L69

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

Wong, I., Shporer, A., Morris, B. M., et al. 2019, AJ, submitted [arXiv:1910.01607]

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

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

Yan, F., & Henning, T. 2018, Nat. Astron., 2, 714

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

Yan, F., Casasayas-Barris, N., Molaverdikhani, K., et al. 2019, A&A, 632, A69

Yelle, R. V. 2004, Icarus, 170, 167

Zhou, G., Bakos, G. Á., Bayliss, D., et al. 2019, AJ, 157, 31

Zorec, J., & Royer, F. 2012, A&A, 537, A120

Zorec, J., Rieutord, M., Espinosa Lara, F., et al. 2017, A&A, 606, A32
**Appendix A: Rossiter-McLaughlin MCMC results**

Here, we show the correlation diagrams for the reloaded
Rossiter-McLaughlin study obtained with the emcee algorithm.

<figure><img src="image_17.png" /><figcaption>**Fig. A.1.** Correlation diagrams in the case of a stellar solid body rotation.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. A.2.** Correlation diagrams for the case of a stellar differential rotation, with an “anti-solar” law.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. A.3.** Correlation diagrams for the case of a stellar differential rotation, with a solar-like law.</figcaption></figure>

**Appendix B: Rossiter-McLaughlin correction for transmission spectroscopy**

Here, we show the effect of the Rossiter-McLaughlin correction
on the Hα line presented in Sect. 3.3.

<figure><img src="image_20.png" /><figcaption>**Fig. B.1.** Hα transmission spectrum in the planet rest frame (gray line and black circles), uncorrected (top) and corrected (bottom) from the Rossiter-McLaughlin effect. The vertical black dashed line shows the line center taking the systemic velocity into account. We show the Gaussian fit (blue) to the corrected line in both sub-panels for the purposes of comparison.</figcaption></figure>
**Appendix C: Transmission spectroscopy MCMC results**

Here, we show complementary correlation diagrams for the Balmer lines NLTE study presented in Sect. 4.2.2.

<figure><img src="image_21.png" /><figcaption>**Fig. C.1.** Correlation diagrams in the case of a complete NLTE atmosphere (each neutral hydrogen electronic level is a separated free parameter). Two families of solution are visible, one being compatible and not so far from the Boltzmann equilibrium (shown in red dashed lines), the other being significantly away from the Boltzmann equilibrium. These two families are studied in the main text.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. C.2.** Correlation diagrams in the case of a complete NLTE atmosphere (each neutral hydrogen electronic level is a separated free parameter). Exploration of the high temperature, near Boltzmann equilibrium (red dashed lines) family solution (see main text).</figcaption></figure>