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

KELT-9 b, the hottest known exoplanet with $T\sim4400$ K, is the archetype 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=7.54\pm0.12$ h) and studied the Rossiter-McLaughlin effect (finding a spin-orbit angle $\lambda=-85.01^{\deg}\pm0.23^{\deg}$) 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=13200^{+800}_{-720}$ K and a mass loss rate of $\dot{M}=10^{12.8\pm0.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.


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 Based on observations made at TNG (Telescopio Nazionale Galileo) telescope with the HARPS-N spectrograph under program A35DDT4 and OPTICON 2018A/038.
Ultra-hot Jupiters are the hottest and the most irradiated hot Jupiters (with a planetary equilibrium temperature above 2,000-2,500 K, Fig. 1). Well-known members of that category include . 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.
KELT-9 b (Gaudi et al. 2017, the hottest exoplanet discovered so far around a main-sequence star, with T eq 4, 000 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 2, 500 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 dayside, 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. 2019). These predictions can change due to heat transport and circulation, depending on whether we observe the day or the night side, or the terminator 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). 2.178 ± 0.011 R † ρ 0.2702 ± 0.0029 g cm −3 † T eff 9 600 ± 400 K § log g 4.058 ± 0.010 cgs † v sin i * Spectro 112.60 ± 0.04 km s −1 ‡ v sin i * RM (SB) 122.5 ± 0.6 km s −1 ‡ v sin i * RM (DR) 116.9 ± 1.8 km s −1 ‡ Notes.
( * ) The parameters and errors are initially taken or computed from the extended data table 3 in Gaudi et al. (2017). ( †) Updated values from Hoeijmakers et al. (2019). ( ‡) Updated values from this work. ( §) Updated values from Borsa et al. (2019). ( ¶) 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 µ 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. Notes.
( * ) In parentheses: the number of spectra taken during the transit and outside the transit, respectively. ( †) The three figures indicate the airmass at the start, middle and end of observations, respectively. ( ‡) The range of signal-to-noise (S /N) ratio per pixel extracted in the continuum near 400 nm and 600 nm. ( §) The range of S /N ratio in the line cores of the Hα, Hβ and Hγ lines. (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 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 to 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(Jensen et al. , 2018Cauley et al. 2015Cauley et al. , 2016Cauley et al. , 2017aCauley et al. ,b,c, 2019Barnes et al. 2016;Casasayas-Barris et al. 2018, 2019Cabot 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 p ) and is escaping from the planet at a rate of ∼10 12 g s −1 . 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. 2017bHoeijmakers et al. 2018Hoeijmakers et al. , 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. 2020). 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 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 Sec 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.

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) 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, gen lines, making it a complementary analysis. We note that no Helium (λ10830Å) has been found in Nortmann et al. (2018). 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).
Article number, page 3 of 21 the orders are merged into a one-dimensional spectrum and resampled (with flux conservation) between 387 nm 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 λ/∆λ = 115, 000 or 2.7 km s −1 . In the present study, we took advantage of the 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.ROOK module (spectRoscopic data tOOl 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 Kausch et al. 2015). Following the method presented in Allart et al. (2017) and applied in Hoeijmakers et al. (2018Hoeijmakers et al. ( , 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

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 . 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. will establish a rigorous formalism and present an accompanying pipeline (but see also Bourrier et al. 2020;Ehrenreich et al. 2020, who already used the same method). Our method is implemented as in on our previous work (Wyttenbach et al. , 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.

On the extraction of a transmission spectrum
Because of the fluctuations of the Earth's atmospheric transmission, ground-based transit spectroscopy observations at highresolution have to be self-normalized. The stellar spectra f (λ, t) are divided by their respective values in a reference wavelength 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 The master spectra in-and out-of-transit F in (λ) = t∈in f (λ, t) and F out (λ) = t∈out f (λ, t) are self-normalized in the same manner, givingF in (λ) andF out (λ). 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;Bourrier et al. 2017aBourrier et al. , 2018: where 1 − δ(t) represents the value of the photometric light curve (white light curve), and δ(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 res (λ, t in ) represents the stellar flux from the planet-occulted photosphere regions, that is to say the local stellar spectrum s local (λ, t) emitted from the region behind the planet at time t. The measurement and modeling of f res (λ, t 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 local (λ, t) is not time-dependent and is equal to the out-of-transit spectrumF out (λ). 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: The −R (λ) denotes the normalized spectrum ratio (averaged over N in in-transit spectra), or the averaged transmission spectrum, and is equal to (R p (λ)/R ) 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 − δ(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 local (λ, 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 4 . We assume that the light source s local,atm (λ, t) coming from this annulus (with a surface ratio Σ atm which is several times 2R P H/R 2 ) is equal in shape to the local stellar spectrum s local (λ, t) (that comes from the surface behind the planet), but differs in intensity from it by a factor Σ atm /δ, meaning that s local,atm /Σ atm = s local /δ. The light source s local,atm is absorbed by the exoplanet atmosphere and constitute s local,atm · s tr /Σ atm , s tr (λ, t) being the transmission spectrum we search to measure 5 . Owing to the absorption by the exoplanet atmosphere, the residuals f res (λ, t), computed with Eq. 1, are not equal to the local stellar spectrum s local (λ, t), but contain the sum of s local (λ, t) and of s tr (λ, t) · s local,atm /Σ atm . Hence, the spectrum ratio, is equal to the surface of the atmosphere limb (alone) times its wavelength dependent absorption, meaning that δ(t) + s tr (λ, t) = (R p (λ, t)/R (λ)) 2 is the transmission spectrum. Taking into account that the shape of the out-of-transit fluxF out (λ) is relative to the continuum reference band λ ref , we get the final transmission spectrum formula relative to λ ref : Placing the observer in the planet rest frame is the last performed operation before the sum. We note that by setting s local (λ, t) = δ(t) ∀λ, 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 highresolution, 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. 2020;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 local (λ, t).
This method relies on the knowledge of the local stellar profile s local (λ, t). One approach is to model s local (λ, 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 Sec. 3.2.1). We emphasize the caveat, particularly during the overlap, of the uncertainty whether the local stellar profile remains constant. 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.

Analysis of the local stellar photospheric CCFs
The analysis of the stellar CCFs is conducted on the CCFs constructed with the aforementioned A0 mask ( 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 puls = 7.54±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 ± 0.04 km s −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 −1 . We note that our CCFs (see Fig. 2) and v sin i * values are similar to those published in the results of Borsa et al. (2019). The instrumental profile (2.7 km s −1 ) with the local stellar CCFs (see Sec. 3.2.2) have an averaged broadening of ξ = 9.17 ± 0.05 km s −1 , making the total broadening ∼113 km s −1 . 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 (ξ v sin i * ). This is because of the high observing cadence, the small impact parameter, and the inclined orbit (see Sec. 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 ± 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 puls = 7.5868 ± 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 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 ±l, the degree of the mode, see e.g., Balona 1986;Aerts et al. 1992Aerts et al. , 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.  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.

"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. (2020); 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 are left with the local stellar CCF alone (s local ; see e.g., Cegla et al. (2016); Wyttenbach et al. (2017); Bourrier et al. (2020) 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 σ. 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 eq , i * , and α which describe a solarlike 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 burnin 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 ± 0.6 km s −1 , and λ = −85.65 • ± 0.09 • (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 ∆BIC = 4.1 in its favor, which is interpreted as "positive evidence" but not a confirmed detection (σ∼2.2). In this scenario, we found α = −0.10 +0.03 . We note that λ is still compatible (at 3 σ) with the rigid rotation case, and that the corresponding v sin i * = 116.9 ± 1.8 km s −1 is more consistent with our spectral fit. Another differential solution, with a ∆BIC = 3.8, was found with α = 0.09 ± 0.03, v eq = 155 ± 18 km s −1 , i * = 55.9 •+11.6 • −7.2 • , and λ = −85.10 • ± 0.25 • , 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 ψ = 89.6 • ± 0.7 and ψ = 84.2 •+0.4 • −0.3 • , 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 (α = 0.08), and an "anti-solar"-type (α = −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 rot = 0.7 ± 0.1 day (∼17.1 h) and P rot = 0.56 ± 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 crit ∼ 416 km s −1 ). With such an equatorial velocity, the angular rotation parameter ω = Ω rot /Ω 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. 2015Ahlers et al. , 2020. For example, with ω = 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 ω = 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 Gaudi 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 configura-tion is favorable to measure star-planet interaction because of anisotropic stellar flux and winds ("gravity-darkened seasons", Ahlers et al. 2020).

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 Gaudi et al. (2017) to normalize the spectra because it is the only one available 6 . To correct the transmission spectra from the local stellar photospheric signal, we must know s local (λ, 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 Sec. 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 δ(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 ±1σ 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 local (λ, t), and the K p from Hoeijmakers et al. (2019) allows us to compute the weighted averaged transmission spectra with Eq. 4 7 . We discard the spectra where the atmospheric signal overlaps with the local stellar signal (-0.027≤Phase≤0.0015), because we cannot fit s local (λ, 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.  Right: Same as for Hα, but for the Hβ line (detection significance of 16.1 σ, see Table 3). 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).

Hγ Hδ
-Hα: The Hα line is detected with a high confidence of 24.3 σ, which makes the planetary trail visible for each phase (Fig.7).  Table 3). ferent 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 instruments, 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 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.

Notes.
( * ) 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. ( †) H is not statistically detected. For Hζ we show the 3σ upper limit assuming a Gaussian line shape with a fixed v = 0 km s −1 and FWHM = 30 km s −1 . 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 P and ∼ 1.5 R P (see Fig 7 and 8). Since this is close to the planetary Roche radius (1.95 ± 0.2 R P ), 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 −1 , 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 −1 . For the rest of our analysis, we investigate the strong and consistent Hα, Hβ, and Hγ signals that contain most of the information.

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, 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, section 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 µ=1.3 and a temperature T =5000 K. For KELT-9 b, a better assumption is an atmosphere dominated by ionized hydrogen (µ=0.66) and a thermosphere temperature of T =10 000 K (Fossati et al. 2018;García Muñoz & Schneider 2019). With these numbers, we find that to explain the FWHM 8 of our Hα detection, one needs an optically thick layer (τ∼40) of excited hydrogen H i(2) with a number density of n∼8.2 × 10 3 cm −3 along a typical line path of L=2 √ 2R P N sc H∼13.8×10 9 cm, where N sc ∼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 ∼ 0.3 R P (H ∼ 6 000 km), which is about half of the line core contrast, the Hα line core have τ 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 ∼ 1.5R P , where τ becomes close to unity 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 linear correlation between the line contrasts (or the altitudes z) and the hydrogen Balmer line oscillator strengths (see the log 10 g f values in Sec. 4.1). Knowing that ∆z = H∆ ln(g f ) = H ln(10)∆ log 10 (g f ), the slope of the z-versus-ln(g f ) linear fit gives us H. We computed that H = 9 600 ± 1 400 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/µ 13 500 ± 2 000 [K/u]. Assuming a µ 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 7 500 and 19 500 K. The lower half temperature range is more probable because the thermal ionization of hydrogen decreases µ 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α 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 8 To compare our FWHM, we removed the additional width contribution by other effects (see Sec. 4.1 and 4.2.1). 9 A composition dominated by molecular hydrogen and helium (µ = 2.34) is probably excluded by an unphysical T 30 000 K. 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.

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 1-D 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 ρ(r) with a hydrostatic or a hydrodynamic structure. For the hydrostatic solution, we have the base density ρ 0 (r = 1) as a free parameter, and consider that the pressure scale-height H = k B T/µg can change with the altitude according to T (r), µ(r), g(r) (see e.g., the π η 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 Tian et al. 2005;Salz et al. 2016;Heng 2017): where F 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) 10 , Oklopčić & Hirata (2018), and Lampón et al. (2020). Our free parameter for the atmospheric structure becomes the atmospheric mass loss rate: where ρ(r) and v(r) solve the Parker equations (Parker 1958;Lamers & Cassinelli 1999) from the sonic point (r s , v s , ρ s , µ s ) down to the lowest radius in the grid. We note thatṀ 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 precomputed 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-20000 K). The mean molecular weight µ(ρ(r), T (r)), the mass fraction m ι (r), and volume mixing ratio n ι (r) are given by the chemical grid, where ι is a chosen element and ionization of interest (H i, He i, Na i, K i, Ca ii, etc.). For example, in the hydrogen case, we know the ratio of ionized to neutral hydrogen n H ii /n H i from reading and interpolating the chemical grid with the correct T and ρ from our atmospheric structure. When the Parker wind solution is in use, we link it with 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) µ through the atmosphere. In that case the n H ii /n H i 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 ι(n) from which the transitions we observed are arising. In the case of H i, 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 H i(2) state. If the atmosphere is in local thermodynamic equilibrium (LTE), the distribution of the electronic states follow the Boltzmann equation (for hydrogen): where E 1 =13.606 eV, and υ(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 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: Then, we can let the different β H i(n) be free parameters (or equivalently all the n H i(n) are free parameters). Another possibility is to force all the n H i(n) to depart by the same amount from LTE, meaning that only one constant free parameter β is used for all n H i(n) . We note that a constant n H i(n) /n H i through the thermosphere is an assumption justified by Huang et al. (2017) andGarcía Muñoz &Schneider (2019).
When an element and its transitions are chosen (in our case the Hydrogen Balmer series), the opacity κ [cm 2 g −1 ] are computed in each layer of the atmosphere according to the line list of Kurucz (1979Kurucz ( , 1992, to the temperature, and to the number density n H i(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: where e is the elementary charge, m e is the electron mass, and c the speed of light in vacuum. The f i j , 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 2 j =0.71, -0.02, and -0.447 for j=3, 4, and 5, respectively. A Voigt profile Φ V , centered on the line transition center ν i j (or λ i j ), 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 (∼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 Å (λ/∆λ ∼ 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 τ = κm ι , wherem ι is the column mass of the element ι in [g cm −2 ], is computed through the line of sight (LOS). The transmission function T = e −τ 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 θ circular sector. The sector grid is evenly separated in cos θ, where θ 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 2-D 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: In the case of KELT-9 b, the planetary rotation have a broadening effect of ∼ 6.1 km s −1 at 1 R P for each hemisphere (∼ 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-D and the computational cost would become too large). The LOS velocity is given by: 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 −1 in the probed layers (i.e., the sonic point r s , where the velocity becomes several km s −1 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 −1 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 a be a Gaussian (with a width of ∼2.7 km s −1 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. 2015Wyttenbach et al. , 2017Pino 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 (δ=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 0.1 s (n r =100), while with rotation, one model takes ∼0.5 s (n r =50, m θ =14). We note that a model with n r =50 and m θ =14 is 1% equal to a model with n r =250 and m θ =28.
In the next section, we apply our PAWN model to the Balmer lines detection in the KELT-9 b atmosphere.

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 −1 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 ρ 0 (r = 1) or the mass loss rateṀ depending on if the structure is hydrostatic or hydrodynamic, respectively (for these latter, we use their log 10 values). We use a uniform prior on T between 1000 and 20000 K. We use a log-uniform prior between -18 and -4 for the log 10 (ρ 0 [g cm −3 ]). We use a log-uniform prior between 8 and 17 for the log 10 (Ṁ [g s −1 ]). 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.

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 (ρ 0 [g cm −3 ]) = −10.74 ± 0.05 (BIC=5016.1, Fig. 10). The corresponding Jeans mass loss rate at the Roche radius is computed to beṀ Jeans ∼ 6.6 × 10 11 g s −1 . 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 (Ṁ PW [g s −1 ]) = 12.8 ± 0.3 (i.e.Ṁ PW 5.9 × 10 12 g s −1 , 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 −1 , 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 ± 1.4 km s −1 ) are not dominated by thermal broadening. Indeed, for T = 13 200 K, the hydrogen thermal broadening is ∼ 25 km s −1 , while the other source of broadening are ∼ 16 km s −1 , ∼ 7 km s −1 , and ∼ 2.7 km s −1 for the rotational (at r = 1.3R P ), orbital blurring, and instrumental broadening, respectively. This means that the first source of broadening (∼ 33 km s −1 ) 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 increas- ing. 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Ṁ PW , which shape is given by the Boltzmann equation.

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 (β H i(n) n LT E H i(n) /n H i ) 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 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 (β H i(n) n LT E H i(n) /n H i ) 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 , 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. 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 (β H i(n) n LT E H i(n) /n H i ) be free parameters, we decide to use a unique departure coefficient β for all the electronic levels, that is β H i(n) = β for all H i(n). For this last MCMC run, a log-uniform prior for log 10 (β) 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 Altitude R /R p Velocity [km/s] LTE non-LTE 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 +800 −720 K, and the mass loss rate iṡ M PW = 10 12.8±0.3 g s −1 . 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 = 9 600 +1210 −1180 K, the mass loss rate isṀ PW = 10 14.3 +1.2 −1.4 g s −1 , and the common LTE departure coefficient is β = 10 −6.1 +1.3 −1.1 .
fit is again better, and have a BIC=5014.8. The ∆BIC = 2.4 from the hydrodynamical-Boltzmann model, meaning a small positive evidence against the LTE model (and only a smaller ∆BIC = 1.3 from the hydrostatic Boltzmann model). In any case, the ∆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 β 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.2 −1.4 g s −1 (between 8 × 10 12 and 3 × 10 15 g s −1 ). Such a mass loss rate could have dramatic consequences (see discussion in Sec. 5). Third, the departure from LTE is 10 −6.1 +1.3 −1.1 (in the 10 −8.7 -10 −2.2 range at 3σ), making the LTE excluded in this framework. Interestingly, our best value is far (but still compatible at ∼2.5 σ) 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 Sec. 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 is a valid assumption, then the thermosphere temperature is constrained to be T = 13 200 +800 −720 K, and the mass loss rate to beṀ PW = 10 12.8±0.3 g s −1 . This mass loss rate is higher than the energy-limited escape rate of 10 10 -10 11 g s −1 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) presented a new escape mechanism that applies to ultra-hot Jupiters orbiting early-type stars. This mechanism shows that it is im-portant to also consider the NUV and optical Balmer continuum and Balmer series irradiation from the star. In this scenario, the thermosphere 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 equilibrium), 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 (3 000 K more than with only the X-FUV-NUV flux). In this case the mass loss rate increases to 10 12 -10 13 g s −1 , the atmosphere undergoes 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 My (Gaudi et al. 2017;Borsa et al. 2019), a mass loss of this rate over the planetary life time represents a total mass loss of ∼0.03-0.05 M Jup (∼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 −3 (García Muñoz & Schneider 2019). Despite the data being too low in S/N to show clear evidences for a departure from the LTE, our NLTE model hints at a lower temperature of T = 9 600 ± 1200 K, and a higher mass loss rate of log 10 (Ṁ PW [g s −1 ]) = 14.3 +1.2 −1.4 . If true, the planet may have already lost 1-2 M Jup , 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  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 Jup , which is ∼2% of the planet mass.
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. Indeed, both effects potentially increase the mass loss rate (Allan & Vidotto 2019). In our NLTE model, we found that the departure coefficient β is on the order of 10 −6 , while theory predicts 10 −3 (only for Boltzmann departure). Thus, our much smaller β coefficient could hint to some other possible LTE departure. Indeed, because β 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 phenomenons (Helling et al. 2019).
Moreover, our NLTE scenario may be able to explain the strong absorption by ionized atomic species, especially the Ca ii  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 (χ 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;   (2020) only if we assume a departure from LTE on the order of ∼10 −2 , or ∼10 −4.5 for our two main structures, respectively (Fig 15). Furthermore, the Ca ii opacities are much larger than the Balmer opacities (by ∼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 β 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 scaleheight, 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 . A last caveat is the modeling of the exosphere itself. Because the Roche radius of KELT-9 b is low (1.95 ± 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 Bourrier et al. 2015Bourrier et al. , 2016Allart 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.

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 λ = −85.01 • ± 0.23 • , and positive evidence in favor of a differential rotation model. We also confirmed the presence of stellar pulsations with a period of P puls = 7.54 ± 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α, Hβ, Hγ, and Hδ absorptions coming from the planet's upper atmosphere, with high confidence (between 24.3 and 4.6 σ) during both nights. The Hα, Hβ, Hγ 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 ± 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 bė M PW = 10 12.8±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 (λ10830Å, 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 Gtype) systems (Jensen et al. 2012(Jensen et al. , 2018Yan & Henning 2018;Casasayas-Barris et al. 2018, 2019Cabot et al. 2020;Chen et al. 2020). Finally, such observations will benefit from the ESPRESSO spectrograph (Pepe et al. 2013(Pepe et al. , 2014. On the one hand, ESPRESSO will increase the S/N of the absorption lines 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.