© The Authors 2024

1 Introduction

The wide sample of more than 5000 detected exoplanets includes a broad selection of different alien worlds spanning from cold super-Earths to extremely hot Jupiter-sized planets. This latter category hosts planets with equilibrium temperatures above 2000 K, which thanks to their extended atmospheres are optimal targets for chemical characterisation. Among the ultra-hot Jupiters (UHJs), we locate KELT-9 b, the hottest planet discovered to date (Gaudi et al. 2017). The exoplanet is in a nearly polar (λ = 85.78 deg), short-period (1.48 days) orbit at a separation of about 0.03 AU from its host star. Also known as HD 195689, the B9.5–A0 host star features an effective temperature in the range of 10000 K, a radius of R*=2.36 R_⊙, and a mass (M_*=2.52 M_⊙) more than twice that of the Sun (Gaudi et al. 2017). The tremendous stellar irradiation experienced by KELT-9 b leads its day-side equilibrium temperature to reach about 4600 K (Gaudi et al. 2017), making this planet hotter than most stars. Hence, its scorching temperature induces the atmosphere to inflate substantially, making this planet an ideal target for transmission spectroscopy.

In the context of exoplanet characterisation, high-resolution spectroscopy is one of the most powerful methods of investigating the atomic and molecular atmospheric make-up. Compared to their space-based counterparts, ground-based spectrographs are able to achieve a higher resolving power, such as the ℛ~115 000 of HARPS-N (Cosentino et al. 2012). At such resolutions, individual lines can be resolved, breaking the degeneracy induced by broad spectral features in low-resolution spectra (Birkby 2018).

A wide variety of metallic species, both ionised and neutral, have been identified in the atmosphere of KELT-9 b. The cross-correlation technique led to the detection of Fe I, Fe II, Ti II (Hoeijmakers et al. 2018), Na I, Cr II, Sc II and Y II signatures (Hoeijmakers et al. 2019), Ca I, Cr I, Ni I, Sr II, and Tb II at the 5σ level, and Ti I, V I, and Ba II above the 3σ level (Borsato et al. 2023). Furthermore, line studies of its primary transit unveiled the presence of O I (Borsa et al. 2021b), Mg I and Fe II (Cauley et al. 2019), Hα (Yan & Henning 2018), ionised calcium (i.e. H&K doublet and near-infrared triplet) (Yan et al. 2019; Turner et al. 2020), the hydrogen Balmer series up to Hδ (Yan & Henning 2018; Cauley et al. 2019; Turner et al. 2020; Wyttenbach et al. 2020), and the Paschen β line (Sánchez-López et al. 2022).

In this study, we combine six primary transits of KELT-9 b observed with the HARPS-N instrument mounted on the Telescopio Nazionale Galileo (TNG). By targeting the hydrogen Balmer line series, we aim to identify and characterise all the single H I lines up to Hζ, after a careful correction of the strong Doppler shadow effect generated by the geometry of the transit. The same framework adopted for the Balmer series has been used to analyse several metal lines looking for both neutral and ionised species. Finally, we aim to compare the observed features with published theoretical transmission spectra of KELT-9b computed accounting for non-local thermodynamic equilibrium (NLTE) effects that have been shown to reproduce well past observations of the Hα and Hβ lines (Fossati et al. 2021).

The paper is organised as follows. In Section 2, we describe the observational data. Section 3 illustrates the methodology employed to analyse the data, including the telluric correction, the Doppler shadow removal, and the transit spectrum construction. In Section 4, we report the results we obtained via the single-line analysis technique, while we compare our transmission spectra with models computed accounting for NLTE effects in Section 5. We discuss the implications of our study in Section 6, before giving a short summary highlighting the main conclusions in Section 7.

2 Observations

We analysed a total of six KELT-9 b primary transits, observed from July 2017 to September 2018 with HARPS-N at the TNG. Four of these transits were acquired in the context of the long-term observing programme at the TNG telescope ‘GAPS2: the origin of planetary systems’ – awarded to the Italian Global Architecture of Planetary System (GAPS) Collaboration (P.I Micela; see Guilluy et al. 2022), while the two remaining ones were collected by other programmes (P.I. Ehrenreich). HARPS-N, being the northern-hemisphere’s twin of ESO’s HARPS/3.6 m telescope, is a fibre-fed cross-dispersed echelle spectrograph that covers the 3830–6900 Å spectral range at an average resolution of ℛ ~ 115 000. For each night of observation, the exposure time per spectrum was set to either 600 s or 300 s, leading to an average signal-to-noise ratio (S/N), averaged over all the spectral orders, of 126 and 88, respectively. In addition to the spectra taken in transit, a number of out-of-transit observations were also acquired. These act as a baseline for the stellar flux and allow one to derive the master-out spectrum, which includes only the stellar contribution (see Section 3). A complete overview of the observations, including the date, programme number, PI, number of in- and out-of-transit spectra, average S/N, and airmass range, is given in Table 1. We also show the S/N as a function of wavelength in Fig. A.1¹. The night reports for the different nights do not give any important observing conditions to remark with a seeing below 1″ expect for Night 4. For the analysis of Hζ, we had to discard eight spectra (out of 202 in-transit spectra used) for which the S/N was too low in the bluest part of the spectrum.

Borsa et al. (2021b) discussed the negligibility of line broadening due to the exposure time (i.e. 4.7, 3.5, and 1.3 km/s on average for exposure times of 400, 300, and 111 s, respectively) for CARMENES data. In our case, combining all six nights, we obtained a mean exposure time of 433s, weighting the intransit exposure times for the mean S/N of each of our 202 in-transit spectra. We therefore included the line broadening in our analysis with the NLTE models, as is discussed in Section 5.

3 Transmission spectra extraction

We employed the 3.7 version of the HARPS-N Data Reduction Software (DRS; Pepe et al. 2002) to do an initial reduction of the HARPS-N raw data, which includes a correction for the blaze function that takes care of both the instrumental blaze and the Earth atmospheric change thanks to the Atmospheric Dispersion Corrector. We operated with S1D spectra that were created by the DRS, starting from the S2D images. During the merging procedures, the DRS takes care of the overlaps between the 69 HARPS-N orders. At this stage, the data contain the stellar signal, the planetary signal, and the telluric contamination, all given in the solar system barycentric reference frame. The wavelength information is given in air. For each transit, we calculated the transmission spectra of the single lines following the procedure described in Wyttenbach et al. (2015), and hence comparing out-of-transit and in-transit spectra.

The first step of the analysis consists of correcting the observed spectra for telluric contamination. To do so, we employed Molecfit² (Smette et al. 2015; Kausch et al. 2015), a specialised ESO tool designed to handle the correction of telluric atmospheric lines in astronomical spectra through the use of a line-by-line radiative transfer model. Although the software was specifically developed to correct data obtained with ESO instruments, in principle it is able to correct spectra from non-ESO ground-based spectrographs as well, and recently ESO released an additional experimental support for HARPS-N. Telluric correction is particularly important in the red part of the HARPS-N spectra, where the Earth atmosphere’s H₂O and O₂ absorption features become relevant. In Fig. 1, we show an example of telluric removal in the vicinity of the Hα region. Molecfit returned telluric-corrected spectra in the solar barycentric reference frame. Therefore, we shifted the spectra in the stellar reference frame by correcting for the star’s systemic radial velocity, υ_sys. The telluric-corrected spectra were then normalised to a common continuum level in a narrow wavelength range around each absorption line of interest. The normalisation around a small region should also remove all the possible wavelength gradients arising from changes in the overall transmission of the atmosphere and other slow changes in the spectrograph.

For each night, we computed the transmission spectra in the stellar reference frame as the ratio between each spectrum and the weighted average out-of-transit spectrum (master-out). By stacking all spectra sorted by phase, we obtained a 2D map (referred to as tomography) such as the one shown in Fig. 2. From this map, we can clearly see how the joint contribution of the Rossiter–McLaughlin effect (RME) (Rossiter 1924; McLaughlin 1924) and the centre-to-limb variations (CLVs), which arise from the stellar rotation and the limb darkening, respectively, is a major contribution that needs to be modelled and removed. The RME and CLV compete with the planetary atmospheric signal to shape what we observe in the 2D map, such as Fig. 2, according to the geometry of the transit.

We corrected for the RME+CLV effects by adopting an analogous method to that described by Yan et al. (2017) and Casasayas-Barris et al. (2017), and already applied in Guilluy et al. (2024), which involves the use of stellar models. We used ATLAS9 stellar models (Kurucz 1992, 2005, 2014, 2017) to compute the disc-integrated stellar model considering the system parameters listed in Table 2. We computed the spectra for the case of a non-rotating star as well as a rotating star. We used PyLightCurve (Tsiaras et al. 2016) to calculate the planet’s path during the transit to evaluate the part of the stellar disc obscured by the planet at each phase (hereafter, obscured region). We then built a 0.01 R_* × 0.01 R_* pixel size grid to approximate the stellar disc and for each phase we performed the following operations:

• for each pixel belonging to the obscured region, we computed the intensity according to the limb darkening angle and the radial velocity shift with respect to the stellar rotation axis;

• we computed the spectrum of the obscured region using the ATLAS9 models, summing the spectrum of each pixel properly shifted by its own radial velocity.

Instead of using different model spectra at different limb darkening angles, μ, and interpolating to compute the right intensity for each pixel, we used an analytical approach. We computed the limb darkening coefficients for our stellar model using the ExoTethys (Morello et al. 2020) function Sail, adopting a quadratic limb darkening law: $$\frac{I_{\lambda }(\mu )}{I_{\lambda }(1)}=1-c_{1,\lambda }(1-\mu )-c_{2,\lambda }(1-\mu {)}^2,$$(1)

where:

• λ indicates a specific spectral bin or passband;

• μ = cos θ, θ being the angle between the line of sight and the normal to the stellar surface;

• I_λ(μ) is the stellar intensity profile and I_λ(1) is the intensity at the centre of the disc;

• c_1,λ and β_2,λ are the limb darkening coefficients.

We note that the different choice of limb darkening law could be a possible cause of discrepancy with other analyses, since the RME+CLV modelling strongly relies on the part of the stellar disc occulted. However, since we do not have any spatially resolved spectrum of the star, it is not possible to state exactly which limb darkening law would better mimic the real observations.

The integrated flux, F_λ, as a function of the limb darkening angle, μ, can be expressed as $$F_{\lambda }=2\pi {\int }_0^1{I}_{\lambda }(\mu )\mu d\mu .$$(2)

By combining Equations (1) and (2) and solving the integral analytically, we obtained I_λ(1). At this point, knowing I_λ(1) and the limb darkening coefficients, we were able to compute the flux for each μ (and hence for each pixel) using Eq. (1). The integrated flux, F_λ, corresponds to the non-rotating model spectrum. The rotational broadening is given by the sum of the different contributions of each pixel properly shifted by their own radial velocity. The observed master-out is broadened by the stellar rotation, and hence we used the rotationally broadened model as a comparison to the master-out to normalise the non-rotating model continuum. We then subtracted the modelled obscured region from the master-out spectra and divided all spectra by this quantity. The average continuum for the RME+CLV models was set to one before applying the correction. An example of the planet absorption for the Balmer series spectra in the stellar reference frame before and after the correction for RME and CLV is shown in Fig. 2. From the comparison between the top and bottom panels of Fig. 2, we can see that, by following this framework, we were able to remove the RME+CLV contaminations, while not changing the planetary signal. This is particular easy to observe in a target such as KELT-9 b with its peculiar transit, but may be more challenging when the two effects overlap, such as in the case of HAT-P-67 b (Sicilia et al. 2024). We have not considered the potential distortion of line profiles caused by gravity darkening. However, it is noteworthy that such distortions are rarely considered in the literature (Borsa et al. 2021c), and are anticipated to be negligible for this target, as is indicated by Cauley & Ahlers (2022), who state that the gravity darkening is expected to be less relevant with respect to the other effects considered (such as RME and CLV) with the present level of achievable precision. Furthermore, we did not include in our models broadening due to the exposure time, since the polar orbit of KELT-9 b causes the planet to obscure the same portion of the star in terms of velocity. In fact, during the 600s exposure the difference in the velocity occulted is below 1 km/s. Once the spectra in the stellar reference frame had been corrected, we shifted all the in-transit spectra to the planet reference frame by correcting for its radial velocity profile in a circular orbit scenario, $$RV(\mathrm{\Phi })=K_{\mathrm{p}}\sin 2\pi \mathrm{\Phi },$$(3)

where Φ is the orbital phase and K_p is the radial velocity of the planet, computed using the parameters listed in Table 2. The final transmission spectrum of each line was calculated as the error-weighted average of all the full in-transit spectra (T₂-T₃) in the planet reference frame. The uncertainties were originally computed as the square root of the DRS-corrected flux and then propagated during the entire process.

Fig. 1 Example of telluric correction on a KELT-9 spectrum observed by HARPS-N on July 31 2017. After normalisation, the DRS-processed data (in red) was telluric-corrected with Molecfit (green line) using the telluric model (in blue), here shifted upwards for clarity.

Fig. 2 Tomographies in the stellar reference frame for the hydrogen Balmer lines. Top panels: raw tomographies before the RME+CLV correction. Top row, left to right: Hα, Hβ, and Hγ; bottom row, left to right: Hδ, Hϵ, and Hζ. The planetary atmospheric absorption can be distinguished as a darker region following the expected planetary radial velocity profile, and the RME as a brighter feature. We also show the planetary radial velocity profile during out-of-transit with a magenta line. The contact points are represented with dashed and continuous white lines. Middle panels: model tomography computed in the manner explained in Section 3. Bottom panels: corrected tomography obtained by dividing the raw one by the model one.

4 Single-line analysis results

For each analysed line, we extracted the transmission spectrum, as is described in Section 3, and then performed a Gaussian fit returning the absorption depth, the full width at half maximum (FWHM), and the velocity shift with respect to the line wavelength corresponding to a zero radial velocity. We used the Python package Scipy CurveFit³, which employs non-linear least squares to fit a function, setting 100 km/s as the velocity range. For some lines, we adopted a smaller wavelength range (e.g. 80 km/s) due to the closeness with other strong lines. We left the continuum offset as a free parameter and we used the retrieved values to normalise the spectra to one, and hence we do not report the continuum values in the best-fit tables. Our definition of significance is based on the retrieved value of the Gaussian fit and it is defined as the ratio between the absorption depth value and its error.

Our procedure of modelling and removing of RME + CLV strongly depends on the size of the obscured region of the stellar disc, which, at the first order, we consider equal to the planetary radius. To include the atmospheric extension, we translated the line absorption depths, returned by the Gaussian fit, into planetary radii and used these values to repeat the RME + CLV removal with an increased size that accounts for both the planetary radius and the atmospheric extension. We repeated the loop, retrieving a new depth expressed in planetary radii, until the threshold convergence (0.001 R_p) was reached. Adopting this procedure, as opposed to Yan et al. (2017), we are including both the RME and the effective size of the planet occulting the stellar disc by using a wavelength-dependent line profile derived from the planet that accounts for the velocity change in the planetary atmosphere’s signal. Therefore, our correction, similar to other ones adopted for several instruments and targets (Casasayas-Barris et al. 2017; Yan et al. 2019; Borsa et al. 2021a; Stangret et al. 2022; Fossati et al. 2023), manages to deal with both the planet’s atmospheric lines and also the RM signal in the overlapping regions.

4.1 Hydrogen Balmer series

We focussed on the hydrogen Balmer series by analysing all the lines observable in the HARPS-N wavelength range, spanning from Hα to Hζ. In Fig. 3, we show the transmission spectra around each Balmer line. We significantly detect all the lines observable in the available wavelength range with the first significant detection of Hϵ and Hζ. The parameters of the Gaussian best fit are listed in Table 3.

4.1.1 List of Balmer lines detected

• Hα: we detect Hα with a significance of ~60σ. Hα has already been detected in this target by Yan & Henning (2018) and Turner et al. (2020) with CARMENES, and by Cauley et al. (2019) using PEPSI. The Balmer series has been analysed by Wyttenbach et al. (2020) using Night 1 and Night 3 data, detecting Hα, Hβ, Hγ, and Hδ, but not Hϵ (tentative detection) and Hζ (nondetection). As is shown in Fig. A.2, our results are in good agreement (<3σ) with Turner et al. (2020) and Wyttenbach et al. (2020), but in disagreement (>5σ) with Cauley et al. (2019) and Yan & Henning (2018).

• Hβ: we detect Hβ with a significance of ~31σ and, as for Hα, our result agrees (2σ) with Wyttenbach et al. (2020), but it is at odds (7σ) with Cauley et al. (2019).

• Hγ and Hδ: these lines have been detected in the atmospheres of just a few exoplanets (HD 189733 b, KELT-9 b, and KELT-20b/MASCARA-2b; Cauley et al. 2016; Cauley et al. 2019; Casasayas-Barris et al. 2019), with Wyttenbach et al. (2020) being the only study to have detected both lines in the atmosphere of KELT-9 b. We detect Hγ and Hδ with a significance of 25σ and 12.2σ, respectively, further supporting Wyttenbach et al. (2020) detection, especially for Hδ.

• Hϵ: we significantly detect Hϵ for the first time in the atmosphere of an exoplanet with a significance of 6.8σ. Hϵ, already marginally detected by Wyttenbach et al. (2020), falls in a region in which two other strong lines are present; namely, Ca II H at 3968.47 Å and Fe I at 3969.25 Å. The Ca II H line is discussed in Section 4.2.3; it is more than 100 km/s away from the core of Hϵ, and hence it does not affect our detection. The Fe I line lies between the Ca II H lines and Hϵ, and we find that is shifted by 8.67 km/s, similar to the other Fe lines (see Section 4.2.8). The region including Ca II H, Fe I, and Hϵ is shown in Fig. A.3 with the corresponding Gaussian fits.

• Hζ: we detect Hζ for the first time in the atmosphere of an exoplanet, with a significance of 5.7σ. We find that it is deeper than the upper limit measured by Wyttenbach et al. (2020), who expected the absorption to be the smallest among the Balmer lines in the HARPS-N range. Instead, we find Hζ to be the second deepest line after Hα, though with a very large uncertainty (almost three to ten times larger than for the other lines). The discrepancy is most likely due to the low S/N in the bluest region covered by HARPS-N that forced us to discard eight spectra because of their low S/N. The low signal significantly affects both the uncertainties and the normalisation process, because of the difficulty of finding a pseudo-continuum region around the line core, which may have led to overestimating the line depth. We also investigated the possibility that the larger depth is due to the presence of other lines in the region; namely, Fe I at 3888.51 Å and Ca I at 3889.10 Å lines. We analysed these lines adopting the procedure described in Section 3 and found that the depth and FWHM are similar to those obtained analysing Hζ, while the velocity shifts in both cases correspond to the exact position of Hζ. Furthermore, the analysis of the Ca I and Fe I lines presented in Sections 4.2.3 and 4.2.8 reveals that the detected lines of those species are shallower compared to those obtained for Ca II at 3889.10 Å and Fe I at 3888.51 Å.

Fig. 3 Balmer series transmission spectra: all the figures have the same range on the y axis to underline the differences between the single lines. The dashed red lines represent the Gaussian best fit, the results of which are listed in Table 3, while the black dots are for a 20× binning. Hζ presents a wider scatter compared to the other lines due to the low S/N in the bluest part of the HARPS-N range. The purple lines show the NLTE broadened models discussed in Section 5. The small gradient of the NLTE model for Hϵ is due to the vicinity with the Ca H line at 3968.47 Å.

4.1.2 Comparison with previous works and night-to-night variability

A detailed discussion about the differences in the Hα and Hβ transmission spectra observed by Cauley et al. (2019), Yan & Henning (2018), Turner et al. (2020) and Wyttenbach et al. (2020) is given by Fossati et al. (2020). Their conclusions stress that the variations could stem from distinct instruments and resolving powers, but more likely they arise from the slightly differing methodologies employed to extract the planetary signal. The major contributions to the discrepancies seem to be the normalisation, the removal of the RME+CLV effects, and the systemic velocity adopted in the different studies, with the latter being responsible for the velocity shifts (see Fossati et al. 2020, for more details). We can extend the conclusions that Fossati et al. (2020) drew for the other works to our results, which fit well in the scenario, considering the slightly different approaches used in the analyses. In Fig. A.2, we show that the Hα results are in good agreement (<3σ) with those of Turner et al. (2020) and Wyttenbach et al. (2020), but at odds with those of Cauley et al. (2019) and Yan & Henning (2018), while Hβ is in good agreement with Wyttenbach et al. (2020), but not with Cauley et al. (2019).

Another plausible explanation for the discrepancy may be the intrinsic presence of in-transit variations along a single transit, as was pointed out by Cauley et al. (2019). Thanks to our sample of six nights, we were able to explore the presence of variability using our method by analysing Hα, Hβ, and Hγ individually for each night. Fig. A.2 and Table A.1 show that the individual nights’ best-fit results scatter around the value obtained by analysing all the nights together. None of the parameters seem to show a trend that depends on the different exposure time used in the first three nights (600s) and in the last three (300s). Focussing on the Hα absorption depths, we see an intrinsic variability of the order of ~30% in our dataset, with depths spanning over a range (0.78 ± 0.05–1.18 ± 0.03) that includes all the results from other works. Since our datasets have two common transit observations (July 31 2017 and July 20 2018), we deepened the comparison with Wyttenbach et al. (2020). Fig. A.5 shows that: (i) our results are generally in good agreement (<3σ) with those of Wyttenbach et al. (2020), except for the velocity shift of Hϵ, which has been previously only tentatively detected; and (ii) the best-fit parameters considering two nights are closer to the ones that we obtained with six nights than the results obtained analysing the single nights separately, hinting that increasing the global S/N plays a major role in the final results.

We conclude that the discrepancies may generally arise from the joint contribution of several factors, such as different instrumentation and datasets, or they originate in the slightly different methods used to extract the planetary signal, as has already been discussed by Fossati et al. (2020), or from the different K_p used. However, the analysis of six nights using the same instrument and the same framework hints at an intrinsic variability that covers the range of the different values retrieved in the different works. The possible causes of this variability are still unclear, but one cannot exclude that it originates in the planetary atmosphere. A similar intrinsic variability has already been noted by Fossati et al. (2023) in analysing HARPS-N data of UHJ KELT-20 b, despite, in their case, the amplitude of the variation not being significant.

4.2 Metal lines

We analysed several metal species by employing the same method used for the Balmer series. We used the NIST Atomic Spectra Database Lines Data⁴ to aid identifying the strongest lines in the HARPS-N range for each species, focussing on both neutral and ionised species. We do not detect any He (listed here despite not being a metal), Li, Be, K, Sr, Y, Ni, and Ba lines. We stress that these species have not been detected on KELT-9 b in previous studies (Hoeijmakers et al. 2018, 2019; Yan & Henning 2018; Cauley et al. 2019; Borsato et al. 2023; Ridden-Harper et al. 2023) neither via single-line analysis nor by using cross-correlation with templates. Instead, we detect Na, Mg, Ca, Sc, Ti, V, Cr, and Fe lines with a significance greater than 3σ. In the following sections, we list our results for each species in order of ascending atomic number. The best-fit parameters are listed in Table A.2.

4.2.1 Sodium

The neutral sodium doublet lines, named Na D1 and Na D2 in this study, are among the most studied and detected lines in single-line analyses due to their capability to probe the upper layers of the atmosphere. The sodium doublet has already been detected in KELT-9 b by Langeveld et al. (2022) using HARPS-N data collected during Night 1 and Night 3. Their results are in excellent agreement with ours, since they retrieve a mean depth of 0.16 ± 0.03%, while we retrieve 0.15 ± 0.02% and 0.16 ± 0.01% for Na D1 and Na D2, respectively. The mean velocity shift they retrieve is −4.1 ± 2.9, while ours is −4.59 ± 0.84 and −6.68 ± 1.05 for Na D1 and Na D2, respectively. The sodium doublet will also be analysed, using the same dataset employed in this work, in a dedicated study (Sicilia et al. in prep.) that aims to perform a population study on several GAPS targets.

4.2.2 Magnesium

We identify a significant absorption on each of the neutral Magnesium b triplet lines (~5167–5183 Å), as has previously been reported by Cauley et al. (2019), plus the Mg I line at 5528.40 Å. Our Mg I triplet line depths are more than 3σ away from the results reported by Cauley et al. (2019), except for the Mg I 5167.23 Å line, which is consistent within 1σ. Instead, our FWHM results are more in agreement (1σ difference) with Cauley et al. (2019), the first line at 5167.3216 Å being 2σ away. The line shifts have remarkably large error bars, similar to Cauley et al. (2019). The discrepancies in the fitted values between our study and Cauley et al. (2019) are probably due to the latter employing lower-resolution spectra from PEPSI (ℛ~50000) compared to HARPS-N. The Mg I lines at 5167.23 and 5172.70 Å are located in the vicinity of the Fe II 5169.0282 Å line, as is highlighted by Hoeijmakers et al. (2019) in their Fig. 8. To avoid any potential contamination from this iron line, we restricted the velocity range around the centre of the reddest line to ± 80 km/s. This excludes the Fe II line, which does not appear in either transmission spectrum (see Fig. A.7). For the first time, we detect an additional Mg I line at 5528.40 Å with a 4.5σ significance.

4.2.3 Calcium

We extracted transmission spectra for five Ca I and two Ca II lines. The best Gaussian fit parameters are listed in Table A.2, while the transmission spectra are shown in Fig. A.8. We detect for the first time five individual lines of neutral calcium with a significance larger than 3σ and with absorption depths of 0.04–0.08% that translate into planetary radii of 1.03–1.06 R_p. The region around these lines is not affected by other species observed in the atmosphere of KELT-9 b, supporting their identification. This is the first detection of individual Ca I lines in the atmosphere of KELT-9 b, though Ca I has been detected in the atmosphere of KELT-9 b by Borsato et al. (2023) through using the cross-correlation technique on HARPS-N Night 1 and Night 3 data along with two nights obtained with CARMENES.

We also detect the Ca II H&K doublet, finding an absorption depth of 0.67±0.05% and 0.54±0.06%, which translates into planetary radii of 1.42 R_p and 1.33 R_p, respectively. Ca II has been detected by Borsato et al. (2023) and Turner et al. (2020) as well. The Ca II H&K lines have been identified by Yan et al. (2019) using both cross-correlation and single-line analysis. Yan et al. (2019) reports a combined H&K line depth of 0.78±0.04 %, corresponding to 1.47±0.02 R_p, in agreement with our measurement of the Ca II H line. We already discussed the fact that the Ca II H line is blended with the Hϵ and Fe I lines, as is shown in Fig. A.3. However, by choosing a narrower range around the line, we were able to avoid the contamination and perform a Gaussian fit.

4.2.4 Scandium

Singly ionised scandium displays a few strong lines in the HARPS-N range, among which we could identify five (see Fig. A.9 and Table A.2) ranging between 1.06 R_p and 1.15 R_p. Despite ionised scandium also being identified through cross-correlation by both Hoeijmakers et al. (2019) and Borsato et al. (2023) with a high significance (>5σ and > 10σ, respectively), this is the first time that Sc II individual lines have been observed. We also searched for neutral scandium signatures, but failed to identify any of them, which supports the non-detection by Hoeijmakers et al. (2019) obtained via cross-correlation.

4.2.5 Titanium

Starting from the NIST database, we searched for the strongest Ti lines in the wavelength range that we have available. We confidently detected ten Ti II lines with a significance greater than 3σ (Table A.2). We confirm the previous three Ti II detections made by Cauley et al. (2019) and uncover eleven new ones. Overall, we find an average line depth of 0.14±0.03% and an average FWHM of 17.34±4.61 km/s for the identified Ti II lines, which are comparable to previous results. Individual values for the detection significance, absorption depth, R_p, FWHM, and v_sys are reported in Table A.2 and shown in Fig. A.10 and Fig. 4. To mitigate blending effects and avoid false detections, we excluded from the analysis some Ti II lines that were too close to other lines of different atomic species.

We could not identify Ti I, despite it being tentatively detected by Borsato et al. (2023), but with a lower significance (~3σ) even with cross-correlation, hinting at the weakness of the signal. Ti II has been previously detected via cross-correlation by Hoeijmakers et al. (2018, 2019) and Borsato et al. (2023) with a significance of ~18σ, ~25σ, and ~17σ, respectively, supporting our single-line detections. We notice a discrepancy between our resulting mean line shift from the zero point and those obtained by Hoeijmakers et al. (2019) and Borsato et al. (2023). We recover a mean velocity shift of −5.19±1.45 km/s, while they reported systemic velocity values in the −18 to −20 km/s range. By subtracting the stellar systemic velocity that we adopted from their result (−17.74 km/s), we obtain υ_sys ~ 0–3 km/s, and hence a difference in retrieved values, which may arise from the different methods and datasets used.

4.2.6 Vanadium

Neither V I nor V II were found by Hoeijmakers et al. (2019), but recently V I was tentatively detected by Borsato et al. (2023) at ~3.5σ. We do not detect any V II single lines, but report the possible detection of one V I line at 4379.23 Å (4.3σ), as is shown in the tomography and in the transmission spectrum in Fig. A.11. In a 2 Å range around the line of interest, we do not find possible contaminating species, supporting the detection of V I in the atmosphere of KELT-9 b.

4.2.7 Chromium

Ionised chromium was successfully observed in previous studies in the atmosphere of KELT-9 b via cross-correlation (Hoeijmakers et al. 2019; Borsato et al. 2023) with a high statistical significance (>7σ). Neutral chromium was also detected at more than 5σ by Borsato et al. (2023), while Hoeijmakers et al. (2019) only report a tentative detection. Our single-line analysis reveals for the first time four Cr II lines with a significance larger than 3σ (see Fig. A.12 and Table A.2). Despite being detected by Borsato et al. (2023) via cross-correlation, we did not detect any neutral chromium line due to the weakness of the features.

Fig. 4 Best-fit parameters for each Ti II line. The red and dashed blue lines represent, respectively, the mean and its uncertainty.

4.2.8 Iron

Iron is the chemical species that displays the most lines in the optical range, and thus it has always been an ideal target for cross-correlation analysis. We focussed mainly on Fe I and Fe II, which have been detected both via cross-correlation (Hoeijmakers et al. 2018, 2019; Borsato et al. 2023) and singleline analysis (Cauley et al. 2019), detecting 10 Fe I and 25 Fe II lines. We also conducted an investigation into the presence of lines corresponding to Fe iii, Fe iv, and Fe v; however, no detections were made in these cases. From Figs. 5 and 6, we can see that Fe I lines are on average shallower and more narrow than Fe II lines. The velocity shift values of the Fe I lines are in agreement among them, except for a few isolated points (which are still within 2σ). Instead, Fe II lines present a larger scatter in velocity shifts. The Fe I and Fe II spectra are shown in Figs. A.13 and A.14, while the best-fit parameters are shown in Figs. 5 and 6.

5 Comparison with non-local thermodynamic equilibrium models

The intense ultraviolet (UV) radiation that the planet experiences as a result of the host star’s high temperature and close orbital separation leads to significant deviations from the local themodynamical equilibrium (LTE; Fossati et al. 2021). These authors highlighted how the Hα and Hβ transmission spectra are in an excellent match with their models once NLTE effects are taken into account in the computation of both the temperature-pressure structure and the transmission spectrum. Hence, the observations can be used to further constrain the theoretical models, and thus the physical properties of the exoplanet atmosphere. Furthermore, as Fossati et al. (2021) mentioned, the NLTE synthetic transmission spectrum can be used to guide future observations aiming to detect features in the observed transmission spectrum, as has been done by Borsa et al. (2021b). The NLTE transmission spectra considered in this work and published by Fossati et al. (2021) were computed using the Cloudy for Exoplanets (CfE) interface, which enables one to use the Cloudy general-purpose NLTE radiative transfer code to model the atmospheric structure of middle and upper planetary atmospheres and generate transmission spectra. All necessary information about CfE, Cloudy, and the NLTE models of KELT-9b can be found in Fossati et al. (2021).

As for the results in Section 4, we divided our comparison with NLTE models into two sections, one for the hydrogen Balmer series and one for the metal lines. Fossati et al. (2021) and Fossati et al. (2020) already compared their NLTE models with those obtained by other works detecting Balmer lines mentioned in Section 4.1. For the analysis of the detected metal lines, we used the NLTE models to support our detections, checking the presence of features in the NLTE spectra in the same region in which we find an absorption signal in the observed transmission spectrum. To further verify that an absorption feature in the NLTE spectrum was due to a specific species, we also used NLTE models computed accounting only for hydrogen and the specific species considered. Hereafter, we refer to the latter as ‘individual NLTE models’, while we call global NLTE models the NLTE transmission spectrum computed accounting for all species.

Then, we carried out an analysis similar to Borsa et al. (2021b), whereby we used different microturbulence (v_mic; velocity of gas on a scale smaller than the pressure scale height, implemented by adding it in quadrature to the thermal velocity) and macroturbulence (v_mac; velocity of gas on a scale larger than the pressure scale height) velocity values to account for any additional line broadening. We created a grid of models with different v_mic values ranging between 1 and 14 km/s (in steps of 1 km/s), and for each model we convolved the transmission spectrum with: (i) the rotational profile of the annulus representing the atmosphere (as done by Brogi et al. 2016), (ii) a step function to include the broadening due to exposure time (with a size of 4.7 km/s, retrieved from the in-transit mean exposure time of 404 s) and, (iii) a Gaussian kernel with an FWHM equal to a v_mac broadening between 1 and 25 km/s, in steps of 1 km/s. We then selected the v_mic-v_mac pair that minimises the χ². To this end, for each spectral feature, we shifted the NLTE synthetic spectrum to match the centres of the observed lines. This was necessary because, when generating spectra, the Cloudy code (used to generate the NLTE spectra; Ferland et al. 2017) does not locate spectral lines at the expected wavelength, but at the centre of the spectral bin in which the line falls, which thus depends on the spectral resolution used to compute the synthetic spectra. This is why the NLTE synthetic spectrum cannot be used to derive line shifts from the observations. We remark that this coarse placement of spectral lines impacts exclusively the computation of the output spectra and not the NLTE radiative transfer, which is instead computed considering the actual line wavelengths. Before the comparisons, all the models are normalised by dividing for a fitted continuum in a small range around the line. The best-fitting NLTE synthetic spectra for each line are shown along with the observed transmission spectra and Gaussian fits from Fig. 3 to Fig. A.14. The best v_mic-v_mac values along with the ${\chi }_{\text{red }}^2$ values are listed in Table 3. We stress that the ${\chi }_{\text{red }}^2$ values that we report in Table 3, Table A.2, and Table A.3 may be underestimated due to the search for v_mic-v_mac models that minimise ${\chi }_{\text{red }}^2$.

Fig. 5 Same as Fig. 4, but for Fe I.

Fig. 6 Same as Fig. 4, but for Fe II.

5.1 Hydrogen Balmer series

Fossati et al. (2021) already presented a comparison between their NLTE and LTE synthetic spectra for Hα and Hβ and the transmission spectra presented by Yan & Henning (2018), Turner et al. (2020), Cauley et al. (2019) and Wyttenbach et al. (2020), noticing a good χ² agreement, and that the NLTE synthetic transmission spectrum fits the observations significantly better than the LTE model. We extended their comparison for Hα and Hβ to our results including the other Balmer lines detected. The ${\chi }_{\text{red }}^2$ reported in Table 3 suggests that NLTE models are in good agreement with our results: the Hα and ${\chi }_{\text{red }}^2$ values are compatible with the ones listed in Fossati et al. (2021) and all the values are below 1 except for Hζ. This discrepancy may be explained by the low S/N of the observations in that spectral region that gives rise to normalisation problems. However, despite the low ${\chi }_{\text{red }}^2$, the NLTE model predicts the existence of the Hζ line, supporting our detection.

5.2 Metals

For the analysis of the metal lines, it is important to consider that the NLTE models have been computed using solar abundances and that the synthetic transmission spectra have been computed using the sub-stellar temperature-pressure profile, which probably overestimates the temperature in the terminator region probed by the observations (Fossati et al. 2021). These two effects may lead to discrepancies in the comparison with the observations in terms of the depth and shape of the lines.

• Na: the Na lines are both predicted by the NLTE models, which, however, seem to expect deeper lines, possibly due to a different abundance.

• Mg: the Mg I triplet b lines detected are predicted by the individual Mg I NLTE model, which also confirms the presence of a Mg I line at 5528.40 Å, corroborating our detections. We stress also that the best v_mic-v_mac pair for this single line is in agreement with the ones found for the two reddest lines of the triplet.

• Ca: the synthetic NLTE spectra predict the existence of both the Ca I lines that we detected, although these two lines are just above the 3σ threshold for detection. The Ca II H line shows a larger broadening than the Ca II K line, as was already found by the Gaussian fit, with a v_mic of 8 km/s and 4 km/s, respectively. – Sc: the comparison with NLTE models and the use of individual species is particularly important for species like scandium that have never been detected via single-line analysis and that could be confused for other lines. In our case, the global NLTE model predicts the existence of all detected Sc II lines, but individual models show that the 5657 Å and 5684 Å lines might be blended with other weaker lines, respectively a close Sc II and Na I, with the latter line not detected in the atmosphere of KELT-9 b.

• Ti: none of the detected lines are predicted by the individual Ti II NLTE models, while few of them may originate in the contamination of other species predicted by the global NLTE model. For the lines for which the global NLTE model does not predict the existence of any feature due to other species, we checked the presence of other lines of species detected in the VALD3 database (Piskunov et al. 1995), which is more complete compared to NIST. We conclude that the detected features may be due to ionised titanium, which has been detected via single-line analysis and cross-correlation by Hoeijmakers et al. (2018, 2019) and Borsato et al. (2023).

• V: neither the global nor the individual NLTE models predict the presence of the V I line at 4379.23 Å that we tentatively detected. The line is present in the NIST database and the reason why the NLTE model does not predict it may be related to the solar abundance adopted by the model or to the assumption on the temperature profile in the terminator region.

• Cr: as for scandium, the comparison with the NLTE model is particularly useful for the identification and detection of chromium lines; the individual Cr II synthetic spectrum predicts the existence of all detected lines showing also a blend with weaker lines (at 4558.64 Å and 4618.8 Å) that seems not to affect our detections. Also, in this case, the individual NLTE model does not forecast the presence of any of the Cr I features, supporting our non-detection.

• Fe: the synthetic NLTE models predict the presence of the detected Fe I lines, except for the line at 4824.17 Å previously detected by Cauley et al. (2019), which is instead due to Cr II. The large FWHM obtained for the line at 5328.03 Å is due to blending with two Fe I lines, while the analysis of the 5269.53 Å line shows the presence of a strong nearby Ca I line. Overall, we report a good match between our Fe II detections and the theoretical NLTE model. We remark that the Fe II lines at 5316.6 Å and 4351.76 Å might be affected by blending.

Fig. 7 Height distribution of chemical species detected expressed in planetary radii as a function of atmospheric temperature. H I and Ca II lie on top, but still well below the Roche (1.95 Rp). The other metals are distributed at lower altitudes. The dashed black line represents the temperature profile presented by Fossati et al. (2021) accounting for NLTE effects. The full lines representing each species height distribution are shifted with respect to the dashed black line for clarity.

6 Discussion

The large number of single lines detected allows us to better constrain the behaviour of KELT-9 b’s atmosphere. In the analyses discussed below, we excluded the lines with large FWHMs that the comparison with NLTE models proved to be blended with other lines. By coupling the depth of the line with the altitude at which it forms, we were able to construct a map of the species detected as a function of altitude. Our observations allowed us to probe the atmosphere up to 1.5 R_p, which is still well below the Roche lobe (1.95R_p). In the upper atmospheric layers, ~1.25 R_p, neutral hydrogen and Ca II H&K form, while ionised iron forms below, extending down to 1.1–1.0 R_p, where it mixes with neutral Fe. Just below 1.2 R_p we find Ti II, while between 1.10 and 1.04 we find scandium and Cr II. Neutral calcium is the species that forms deeper in the atmosphere, just above 1.03 R_p. Using the temperature–radius profile from Fossati et al. (2021), we can couple the heights that we retrieved from the Gaussian fits with the modelled atmospheric temperature obtained accounting for NLTE effects, as is shown in Fig. 7. We find that, in the upper atmosphere, H I Balmer lines form at a temperature higher than 7500 K, while Ca II lies above 8200 K. Slightly below 1.3 Rp, we find Fe II ranging from 5000 K to 8000 K. Ti II and Fe I show a common range in height corresponding to a temperature of 4800–6600 K. Sc II and Mg I lie between 5000 K and 6000 K, while Ca I and Cr II, which lie lower in the atmosphere, have temperatures below 5000 K.

The NLTE model’s predicted number densities for each neutral and ionised species support our detections. As a matter of fact, the general trend is that ionised species have higher abundances compared to the neutral ones (see Fig. 8), corroborating the detection of Sc II, Ti II, and Cr II. We are not able to find ionised magnesium due to the lack of strong lines in the HARPS-N wavelength range, while for Ca and Fe we detect both the neutral and the ionised lines, probably because of their higher abundance compared to the elements for which we find only ionised species.

Fig. 9 shows the possible relations among our Gaussian best-fit parameters. The distribution in height has already been discussed before, but from the histogram we can clearly see how most of the lines detected are below 1.15 planetary radii, while only a few lines of Fe II, Ca II and the Balmer lines lie above. Overall, we find all lines to be blueshifted, with the distribution peaking around −6 km/s and ranging from 0 km/s to −11 km/s, except for two points above 0 km/s. We do not find a clear correlation between the velocity shift and the height distribution, suggesting the homogeneous motion of the atmosphere in the probed layers. The systematic blueshift is expected by the presence of night-side to day-side winds due to the fact that we are observing the planet terminator. General circulation models for HJ predict wind speeds of typically 2–4 km/s, which is slightly smaller than the 6 km/s that we retrieve. The FWHM of the lines detected returns a quasi-symmetric distribution peaked at ~20 km/s, with a tail due to the Ca II and H Balmer lines, which have a different physical origin since their broadening is mostly due to collisional effects rather than rotational ones.

We focussed on the possible correlation between the FWHM and the height distribution, comparing the line widths and the height distribution in Fig. 9 in a similar way to that done by Borsa et al. (2021a). We excluded the H I and Ca II lines for the reason explained above. Assuming the planet and its atmosphere rotate to be a rigid body with a rotation period equal to 6.6 km/s (tidally locked), we obtained the velocity-radius profile shown by the purple line in Fig. 9. Despite most of our points being at lower radii, our results follow the tidally locked profile, supporting the expectation of tidally locked rotation. In our analysis, we retrieved the FWHM from a Gaussian fit. Although this does not take into consideration the fact that during the transit only an annulus is visible, and that another profile would be required (Brogi et al. 2016), we decided to fit a Gaussian profile as we aim for a qualitative comparison between the statistics inferred from our results and the expected rotational profile shown in purple in Fig. 9.

The impact of NLTE effects in the atmosphere of KELT-9 b has already been discussed by various studies (Fossati et al. 2020; Fossati et al. 2021), and the inclusion of NLTE effects in the modelling scheme has been found to be necessary to reproduce the observations. The key impact of NLTE effects in the planetary atmosphere is the significant increase in temperature in the middle and upper atmosphere. The overpopulation of Fe II drives heating, significantly increasing the absorption of stellar near-UV radiation, where the host star’s spectral energy distribution peaks, further increasing the heating rate. This process leads to the strong temperature inversion that we see in Fig. 7.

We compared our line fitting results to NLTE models, finding that the line profiles, including the line depths, are well described by the NLTE synthetic spectra. Furthermore, NLTE models have played a key role in identifying and confirming the single-line detections presented here, hinting at the blending of different lines, which allowed us to explain why few lines show large FWHM values. We computed the mean v_mic and v_mac for each species, obtaining the results listed in Table 4.

Borsa et al. (2021b) found v_mic = 3.0 ± 0.7 km/s and v_mac = 13 ± 5 km/s, analysing the oxygen triplet at about 7770 Å with CARMENES data. The mean values of v_mic and v_mac that we retrieve including all the species except for H I, Ca II, and Ti II agree with their results. We investigated the possible correlation between v_mic and v_mac and the height distribution, shown in Fig. 10, and we did not find any significant correlation. We stress that we consider v_mic and v_mac to be fudge parameters needed to explain the observed extra broadening, but they still do not have a well-defined physical meaning. The v_mic values, considering their uncertainties, are compatible with zero.

Fig. 8 Number densities predicted by the NLTE model of Fossati et al. (2021) for the detected atoms in neutral and singly ionised form. The dashed lines indicate the height probed by the observations. Generally, the ionised species have higher number densities compared to the neutral ones.

7 Conclusions

In summary, we analysed six transit observations of the UHJ KELT-9 b taken with the HARPS-N spectrograph mounted at the TNG. The elevated temperature of the planet results in the significant extension of the atmosphere and leads to an increase in the S/N of the data, making it the perfect target for transmission spectroscopy. We employed the single-line analysis technique to search for a variety of metal lines across different species, ranging from the lightest elements, such as hydrogen, to the heaviest ones, such as iron. Thanks to the strong atmospheric signal due to the combination of six observational nights, we were able to identify 70 individual lines belonging to seven different species, Na I, Ca I, Ca II, Fe I, Fe II, Mg I, Ti II, Sc II, and Cr II. This is more than any other study utilising this technique has ever found. Our approach allowed us to confirm previous detections and make new ones, such as the discovery of Hϵ and Hζ, which had only been deemed tentative up to now. Hζ displays a larger depth than expected and we stress that our results might be influenced by normalisation problems arising from the low S/N in the bluest part of the wavelength range covered by HARPS-N. We also detected Cr II and Sc II single lines for the first time in the atmosphere of an exoplanet, along with many new lines from already-detected species. Among the chemical species that we looked for, we could not find He, Li, Be, K, Sr, Y, Ni, Ba, and Mn single lines. Even though the atmospheric signal is strong, with an average S/N of ~100, one possible reason for these non-detections is the weakness of those lines, which might not be sufficient for a single-line analysis study.

Single-line analysis in transmission has proven to be a powerful tool to support and confirm detections made with cross-correlation studies, while delivering important physical information about the width and depth of the lines, which in turn provides clues about the turbulence and the stratification of the atmosphere, in addition to information on the physical and chemical structure of the atmosphere. By juxtaposing our absorption lines with NLTE models, we could verify the presence of detected features and distinguish regions in which line blending may be occurring. Most of the lines analysed experience a velocity blueshift, corroborating the existence of night-side to day-side winds in the atmosphere, especially for the lines forming between 1.07 and 1.20 planetary radii. The rotational velocity and the height distribution generally agree with the hypothesis of a tidally locked rigid rotating body. In conclusion, our study marks a significant contribution to the application of the single-line analysis technique: we yielded an unprecedented number of individual lines, detecting many of them for the very first time in an exoplanetary atmosphere. The huge number of lines detected allowed us to further corroborate the presence of winds and to clarify the atmospheric stratification of KELT-9b.

Fig. 9 Analysis of Gaussian best-fit parameters of all detected lines. Top panel: histograms of height, velocity, and FWHM. Middle panel: relations between velocity, height, and FWHM. Bottom panel: rotational velocity against the height in the atmosphere (in planetary radii); the dashed purple line represents the profile expected for the tidally locked scenario, while the black line represents the planetary radius.

Fig. 10 Distribution of vmic and vmic as a function of height, for the different species detected. We also report, as a comparison, the values retrieved by Borsa et al. (2021b) with CARMENES data.

Data availability

An extended version of the manuscript with an Appendix including Tables A.1–A.3 and Figures A.1–A.14 can be found on Zenodo https://zenodo.org/records/13348484.

Acknowledgements

The authors acknowledge financial contribution from PRIN INAF 2019 and from the European Union – Next Generation EU RRF M4C2 1.1 PRIN MUR 2022 project 2022CERJ49 (ESPLORA), PRIN MUR 2022 project PM4JLH (“Know your little neighbours: characterising low-mass stars and planets in the Solar neighbourhood”) and PRIN MUR 2022 (project no. 2022J7ZFRA, EXO-CASH). The authors acknowledge the support of the ASI-INAF agreement 2021-5-HH.O and the project “INAF-Astrofisica Fondamentale GAPS2”. AS acknowledges financial contribution from the University College London MAPS PGR Travel Grant 2021/22, the European Research Council under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 758892, ExoAI) and the European Union’s Horizon 2020 COMPET program (grant agreement no. 776403, ExoplANETS A). Furthermore, AS is supported by a UKRI-STFC CDT-DIS studentship under project no. 541466. GMa acknowledges support from CHEOPS ASI-INAF agreement no. 2019-29-HH.0. Part of the research activities described in this paper were carried out with contribution of the Next Generation EU funds within the National Recovery and Resilience Plan (PNRR), Mission 4 – Education and Research, Component 2 – From Research to Business (M4C2), Investment Line 3.1 – Strengthening and creation of Research Infrastructures, Project IR0000034 – “STILES – Strengthening the Italian Leadership in ELT and SKA”. This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. We sincerely thank the Referee for their comments stressing some important factors and analyses that we did not consider initially and that we think improved the quality of the manuscript. Author contributions. MCD led the data analysis for this project with contributions from AS and FB. LF provided and contributed to the analysis of NLTE atmospheric models. MCD wrote the manuscript along with AS, FB, LF and GM. All authors discussed the results and commented on the draft.

References

All Tables

All Figures

	Fig. 1 Example of telluric correction on a KELT-9 spectrum observed by HARPS-N on July 31 2017. After normalisation, the DRS-processed data (in red) was telluric-corrected with Molecfit (green line) using the telluric model (in blue), here shifted upwards for clarity.
In the text

Fig. 2 Tomographies in the stellar reference frame for the hydrogen Balmer lines. Top panels: raw tomographies before the RME+CLV correction. Top row, left to right: Hα, Hβ, and Hγ; bottom row, left to right: Hδ, Hϵ, and Hζ. The planetary atmospheric absorption can be distinguished as a darker region following the expected planetary radial velocity profile, and the RME as a brighter feature. We also show the planetary radial velocity profile during out-of-transit with a magenta line. The contact points are represented with dashed and continuous white lines. Middle panels: model tomography computed in the manner explained in Section 3. Bottom panels: corrected tomography obtained by dividing the raw one by the model one.

In the text

Fig. 3 Balmer series transmission spectra: all the figures have the same range on the y axis to underline the differences between the single lines. The dashed red lines represent the Gaussian best fit, the results of which are listed in Table 3, while the black dots are for a 20× binning. Hζ presents a wider scatter compared to the other lines due to the low S/N in the bluest part of the HARPS-N range. The purple lines show the NLTE broadened models discussed in Section 5. The small gradient of the NLTE model for Hϵ is due to the vicinity with the Ca H line at 3968.47 Å.

In the text

	Fig. 4 Best-fit parameters for each Ti II line. The red and dashed blue lines represent, respectively, the mean and its uncertainty.
In the text

	Fig. 5 Same as Fig. 4, but for Fe I.
In the text

	Fig. 6 Same as Fig. 4, but for Fe II.
In the text

Fig. 7 Height distribution of chemical species detected expressed in planetary radii as a function of atmospheric temperature. H I and Ca II lie on top, but still well below the Roche (1.95 Rp). The other metals are distributed at lower altitudes. The dashed black line represents the temperature profile presented by Fossati et al. (2021) accounting for NLTE effects. The full lines representing each species height distribution are shifted with respect to the dashed black line for clarity.

In the text

	Fig. 8 Number densities predicted by the NLTE model of Fossati et al. (2021) for the detected atoms in neutral and singly ionised form. The dashed lines indicate the height probed by the observations. Generally, the ionised species have higher number densities compared to the neutral ones.
In the text

Fig. 9 Analysis of Gaussian best-fit parameters of all detected lines. Top panel: histograms of height, velocity, and FWHM. Middle panel: relations between velocity, height, and FWHM. Bottom panel: rotational velocity against the height in the atmosphere (in planetary radii); the dashed purple line represents the profile expected for the tidally locked scenario, while the black line represents the planetary radius.

In the text

	Fig. 10 Distribution of vmic and vmic as a function of height, for the different species detected. We also report, as a comparison, the values retrieved by Borsa et al. (2021b) with CARMENES data.
In the text