© ESO 2020

1 Introduction

More than twenty years of exoplanet research has greatly broadened and deepened our understanding of planets. In particular, the detection of close-in giant planets has opened the possibility of characterising planetary atmospheres (e.g. Seager & Sasselov 2000; Brown 2001; Charbonneau et al. 2002; Snellen et al. 2010; Brogi & Line 2019). Atmospheric characterisation has become one of the central aspects of exoplanet science such that current (e.g. TESS) and future (e.g. PLATO) planet detection facilities are aimed at finding planets orbiting stars that are bright enough to enable transmission and/or emission spectroscopy (e.g. Rauer et al. 2014; Ricker et al. 2015).

The atmospheric characterisation of observations conducted so far has led to the detection of a wide range of atomic and molecular species in exoplanetary atmospheres (e.g. Vidal-Madjar et al. 2004; Redfield et al. 2008; Fossati et al. 2010; Deming et al. 2013; Barman et al. 2015; Kreidberg et al. 2015; Line & Parmentier 2016; Sing et al. 2016; Casasayas-Barris et al. 2017; Evans et al. 2017; Lendl et al. 2017; Line et al. 2017; Birkby et al. 2017; Wyttenbach et al. 2017; Brogi et al. 2018; Jensen et al. 2018; Nikolov et al. 2018; Nortmann et al. 2018; Tsiaras et al. 2018; Wakeford et al. 2018; Hoeijmakers et al. 2019; von Essen et al. 2019). Furthermore, the community has devised sophisticated forward models and retrieval techniques to extract relevant physical information (e.g. temperature structure, abundance profiles) from the observations (e.g. Irwin et al. 2008; Madhusudhan & Seager 2009, 2010; Howe & Burrows 2012; Lee et al. 2012; Line et al. 2013; Waldmann et al. 2015a,b; Lavie et al. 2017; MacDonald & Madhusudhan 2017; Malik et al. 2017; Mollière et al. 2015, 2019; Shulyak et al. 2019; Barstow et al. 2013, 2017, 2020). The main general result is that exoplanetary atmospherespresent a large variety of physical and chemical properties. The picture is further complicated by the fact that most planets hostaerosols which hide spectral features and/or sequester elements producing detectable features (e.g. Sing et al. 2016; Parmentier et al. 2016; Gibson et al. 2017; Wakeford et al. 2017; Benneke et al. 2019). Therefore, for some hot Jupiters, a thorough observational atmospheric characterisation may be limited to the upper atmospheric layers.

However, ultra-hot Jupiters, gas giant exoplanets with equilibrium temperatures exceeding 2000 K and typically orbiting intermediate-mass stars (A- and early F-type), are not subject to this limitation, because their high atmospheric temperature and intense irradiation ensure that the molecules composing aerosols are dissociated, at least on the day-side (e.g. Parmentier et al. 2018; Kitzmann et al. 2018; Lothringer et al. 2018). About a dozen ultra-hot Jupiters have been detected to date and the most studied one is KELT-9b, also known as HD 195689 b, which is the hottest planet orbiting a non-degenerate star among those found so far (Gaudi et al. 2017). Both day- and night-side temperatures of KELT-9b have been measured through phase curve observations at optical (with TESS) and infrared (with Spitzer at 4.5 μm) wavelengths. From TESS photometry, Wong et al. (2020) obtained a night-side temperature of 3020 ± 90 K and a day-side temperatureof 4570 ± 90 K, while from Spitzer data Mansfield et al. (2020) derived a night-side temperature of 2556 K and a day-side temperatureof 4566 K. Both hydrogen and several metals have been detected in the atmosphere of KELT-9b through ground-based high-resolution transmission spectroscopy (e.g. Yan & Henning 2018; Yan et al. 2019; Hoeijmakers et al. 2018, 2019; Borsa et al. 2019; Cauley et al. 2019; Turner et al. 2020; Pino et al. 2020; Wyttenbach et al. 2020). The detection of metals, in particular ofions such as FeII, indicates that the atmospheric temperature in the region where lines of these elements form is above 4000 K (Hoeijmakers et al. 2019).

Self-consistent atmospheric modelling of KELT-9b conducted with the PHOENIX stellar and planetary model atmosphere code indicates that the temperature profile should be inverted, reaching temperatures of the order of 6300 K in the upper atmosphere (Lothringer et al. 2018; Fossati et al. 2018). Pino et al. (2020) detected such a thermal inversion from secondary eclipse high-resolution spectroscopic observations. García Muñoz & Schneider (2019) modelled the upper atmosphere of KELT-9b assuming a pure hydrogen composition accounting for stellar irradiation, with a particular focus on the heating caused by the near-ultraviolet photons, and accounting for non-local thermodynamic equilibrium (NLTE) effects. They obtained a temperature profile reaching a maximum of about 15 000 K at a pressure of 10⁻⁹ bar, but mentioned that the implementation of cooling by FeII would decrease the temperature by 1000–2000 K, presumably more with the inclusion of additional metals. They also showed that in their model NLTE effects play a significant role in shaping the properties of the planetary upper atmosphere.

The Hα hydrogen Balmer line has been the first one detected in the atmosphere of KELT-9b, shortly followed by some of the higher-order hydrogen Balmer lines (e.g. Yan & Henning 2018; Cauley et al. 2019; Turner et al. 2020; Wyttenbach et al. 2020). The detection of the Hα line revealed that the planet hosts an extended, hot hydrogen envelope, while the large transit depth gave rise to the idea that Hα could be used to directly probe the extended, escaping planetary upper atmosphere, hence constrain mass loss (Yan & Henning 2018; Fossati et al. 2018; García Muñoz & Schneider 2019). Turner et al. (2020) show that the Hα line does not probe the planetary atmosphere beyond the Roche lobe; however this does not prevent one from constraining the atmospheric mass-loss rate (Yan & Henning 2018; Fossati et al. 2018; Wyttenbach et al. 2020). As a matter of fact, the detection and modelling of the hydrogen Balmer lines enable one to constrain the atmospheric temperature between about the mbar and the nbar pressure level (García Muñoz & Schneider 2019; Wyttenbach et al. 2020), hence the energetics driving the escape.

García Muñoz & Schneider (2019), Turner et al. (2020), and Wyttenbach et al. (2020) employed the observed Hα line profile of KELT-9b to extract information on the planetary atmospheric properties. García Muñoz & Schneider (2019) and Turner et al. (2020) showed that the Hα synthetic transmission spectrum computed accounting for NLTE effects and with the system parameters obtained by Gaudi et al. (2017) was significantly weaker than the observed one. García Muñoz & Schneider (2019) concluded that this was possibly the result of an overestimation of the planetary mass, particularly because at that time the mass measurement was affected by a large uncertainty (2.88 ± 0.84 M_Jup; Gaudi et al. 2017). However, additional and more accurate radial velocity measurements confirmed the planetary mass derived by Gaudi et al. (2017), further providing a significantly smaller uncertainty (2.88 ± 0.35 M_Jup; Borsa et al. 2019). In order to fit the observed Hα line profile, Turner et al. (2020) employed a smaller planetary mass by about 30%, but they suggested that the line profile may be also reproduced by increasing the temperature compared to that provided by PHOENIX, instead of decreasing the planetary mass. Wyttenbach et al. (2020) employed a retrieval-like technique to reconstruct the atmospheric temperature, mass-loss rate, and density of excited hydrogen by fitting the observed Hα, Hβ, and Hγ lines. They assumed an isothermal profile and a constant density ratio of excited hydrogen to total hydrogen, and did not consider photoionisation. By keeping the hydrogen density profile equal to that given by the Boltzmann equation, hence assuming local thermodynamic equilibrium (LTE), they obtained an atmospheric temperature of 13 200 K. Instead, by leaving the density of excited hydrogen as a free parameter, they obtained an atmospheric temperature of 9600 ± 1200 K and a density of excited hydrogen relative to the total amount of hydrogen of about 10⁻¹¹. However, there are significant degeneracies among the free parameters.

In this work, we construct a large number of atmospheric temperature-pressure (TP) profiles, ranging from about 1 to 10⁻¹¹ bar, and employ the Cloudy NLTE radiative transfer code (Ferland et al. 2017) to constrain the TP profile, or family of TP profiles, that best reproduce the available transmission spectroscopy observations of the Hα and Hβ lines. Observationally constraining the planetary atmospheric TP profile is an important step towards the characterisation of the prototype ultra-hot Jupiter and it significantly advances models of the planetary upper atmosphere, further helping to place ultra-hot Jupiter escape on more solid observational ground. The knowledge, even if approximate, of the atmospheric TP profile would also enable one to produce more accurate synthetic transmission spectra to be employed for detecting metals in the planetary atmosphere, for example, to refine the inference of the local thermodynamic state and to better constrain physico-chemical processes in the atmospheres of ultra-hot Jupiters.

This paper is organised as follows. In Sect. 2, we present the considered observational material and compare the available Hα and Hβ line profiles witheach other. Section 3 describes how we computed the synthetic TP profiles, while Sect. 4 presents the modelling scheme we employed to calculate the synthetic transmission spectra from the TP profiles. Section 5 presents the results that we discuss in Sect. 6. We gather the conclusions of this work in Sect. 7.

2 Observed transmission spectra

Transmission spectra of KELT-9b covering at least one of the hydrogen Balmer lines, typically Hα, have been published by Yan & Henning (2018), Cauley et al. (2019), Turner et al. (2020), and Wyttenbach et al. (2020). Both Yan & Henning (2018) and Turner et al. (2020) obtained the Hα transmission spectrum using the CARMENES high-resolution échelle spectrograph (R ≈ 95 000) attached to the 3.5 m telescope of the Calar Alto Observatory. The transmission spectrum of Yan & Henning (2018) was obtained combining transits observed during two nights in 2017 (August and September), while that of Turner et al. (2020) was obtained from one transit observation recorded in 2018 (June). Cauley et al. (2019) obtained the Hα and Hβ transmission spectra of KELT-9b following one transit observation collected in 2018 (July) using the PEPSI high-resolution échelle spectrograph (R ≈ 50 000) attached to the Large Binocular Telescope. We employ here the transmission spectrum obtained considering only the first third of the transit observed by Cauley et al. (2019), because the rest of the transit appears to be characterised by variability. Wyttenbach et al. (2020) obtained transmission spectra of the Hα to Hδ lines following two transit observations conducted in July 2017 and July 2018 with the HARPS-N high-resolution spectrograph (R ≈ 115 000) attached to the 3.6 m Telescopio Nazionale Galileo.

Yan & Henning (2018) and Cauley et al. (2019) normalised the already blaze-corrected spectra using linear fits to adjacent continuum points, finally stitching the spectral orders together. They further corrected telluric contamination employing a theoretical telluric spectrum, finally dividing each spectrum by the average out-of-transit spectrum and then shifting all spectra to the planet rest frame before co-adding the in-transit residual spectra. Turner et al. (2020) applied SYSREM (Tamuz et al. 2005) to remove telluric and stellar lines and other systematics from the blaze-corrected spectra, which have then shifted to the planet rest frame before co-adding. Wyttenbach et al. (2020) followed spectral reduction and analysis procedures similar to those of Yan & Henning (2018) and Cauley et al. (2019), but they extracted the transmission spectrum computing the difference between in- and out-of-transit spectra to correct for the Rossiter-McLaughlin (RM) effect, and only in a second step, the residuals have been divided by the average out-of-transit spectrum. All authors corrected the data for the RM effect, while Yan & Henning (2018) and Cauley et al. (2019) further corrected the data for the centre-to-limb variation (CLV) effect (see, for example Yan et al. 2017).

Figure 1 compares the Hα (right panel) and Hβ (left panel) transmission spectra published by Yan & Henning (2018), Cauley et al. (2019), Turner et al. (2020), and Wyttenbach et al. (2020). To ease the comparison, in the making of Fig. 1 we removed the small wavelength shifts present among the different line profiles, which are probably caused by the use of different systemic radial velocity values. Because of differences in the number of observed transits or in the collecting area of the telescopes employed to obtain the data, the noise level in the transmission spectra of Yan & Henning (2018), Cauley et al. (2019), and Wyttenbach et al. (2020), as indicated by the error bars, is significantly smaller than that of Turner et al. (2020).

The measurements of the full width at half maximum (FWHM) given by the Gaussian fits indicate that there are significant differences among the observed line widths (see also Wyttenbach et al. 2020, for a similar discussion). As indicated in Table 1, the FWHM obtained from the Hα profile of Yan & Henning (2018) is about 1.5 times larger than those obtained from the profiles of Cauley et al. (2019) and Turner et al. (2020), while that of Wyttenbach et al. (2020) lies in between those. Also for the Hβ line, the FWHM derived from the profile of Wyttenbach et al. (2020) is significantly larger than that obtained from the profile of Cauley et al. (2019). The Hα and Hβ profiles of Wyttenbach et al. (2020) are significantly shallower than those of the other authors. These differences are not due to differences in the spectral resolution of the instruments or data, but may stem out of differences in the data reduction and analysis procedures, although for example Cauley et al. (2019) accounts for the CLV effect, while Turner et al. (2020) does not, but the two profiles are comparable. Differences may, however, rise also from the treatment of the Doppler shadow (Borsa et al. 2019; Wyttenbach et al. 2020) or from the fact that the different authors considered significantly different systemic velocities. Such differences could also arise as a result of instrumental systematics, which are different for each instrument and for which they have not been properly corrected. Finally, it is also possible that the observed differences are of an astrophysical nature and due to intrinsic variations of the hydrodynamically escaping atmosphere. Indeed, Cauley et al. (2019) report on the possible presence of transit depth variations along a single transit. Table 1 also shows that the equivalent widths derived from the profiles of Wyttenbach et al. (2020) and Cauley et al. (2019) are comparable, which is because the former are shallower and broader, while the latter are deeper and narrower.

For the comparison of the synthetic spectra with the observations, we decided to consider only the transmission spectra obtained by Cauley et al. (2019) and Wyttenbach et al. (2020). This is because Yan & Henning (2018) and Turner et al. (2020) provided only the Hα line, the profile of Turner et al. (2020) is much noisier than the others, and the profile of Yan & Henning (2018) appears to be significantly different (deeper and broader) from the others. To account for the differences in the Hα and Hβ line profiles and to increase the signal-to-noise of the transmission spectra, we aligned the profiles of Cauley et al. (2019) and Wyttenbach et al. (2020) employing the line centres obtained from the Gaussian fits and computed the mean profiles using a weighted average. Therefore, in the end, the comparison with synthetic spectra is based on the profiles obtained by Cauley et al. (2019) and Wyttenbach et al. (2020), and on the weighted average profiles obtained considering the results of Cauley et al. (2019) and Wyttenbach et al. (2020).

Fig. 1 Transmission spectra of the hydrogen Balmer line profiles (Hβ, left; Hα, right) published by Yan & Henning (2018, green), Cauley et al. (2019, black), Turner et al. (2020, blue), and Wyttenbach et al. (2020, orange). The magenta line is the weighted average profile obtained considering just the profiles of Cauley et al. (2019) and Wyttenbach et al. (2020). Wavelengths are in vacuum to match those of the synthetic transmission spectra (see Sect. 4). The data have been rebinned by a factor of five for visualisation purposes. The dashed lines are Gaussian fits to the data. The horizontal dotted line at one sets the continuum level to guide the eye.

3 Temperature-pressure profiles

Lothringer et al. (2018) employed the PHOENIX stellar and planetary atmosphere code to compute TP profiles of ultra-hot Jupiters showing that they present strong temperature inversions in the middle part of the atmosphere, at pressures between 10⁻¹ and 10⁻⁵ bar, caused by absorption of the intense UV and optical stellar radiation (see also García Muñoz & Schneider 2019). Fossati et al. (2018) used a PHOENIX TP profile, specifically computed for KELT-9b, to estimate the planetary mass-loss rate. Also this profile shows a rather steep temperature inversion, with the increase occurring between 10⁻¹ and 10⁻⁶ bar (Fig. 2). The presence of an inverted temperature profile has been observationally confirmed by Pino et al. (2020).

Therefore, we took inspiration from the presence of the inversion around the mbar level to draw a large number of empirical TP profiles, which are then used as a basis to compute the synthetic transmission spectra (Sect. 4) to be compared to the observations (Sect. 5). We computed 126 empirical TP profiles using a modified version of the Three-channel Eddington Approximation (TCEA) temperature profile model by Guillot (2010), in the form used by Line et al. (2013). In this approximation, the temperature profile (T) as a function of the atmospheric pressure (p) is parametrised as (1)

where t_int is the planetary internal temperature, (2) (3)

and (4)

The parameters of the model that we vary to modify the shape of the TP profile are κ, γ, β, and s. The parameter κ, which is related to the Planck thermal infrared opacity, shifts the TP profile in pressure. The parameter γ, which is related to the visible-to-thermal stream Planck mean opacity ratio, controls the difference in temperature between the top and bottom of the atmosphere. The parameter β, which is related to various atmospheric temperature distribution effects (e.g. albedo, emissivity, day-night redistribution), shifts the TP profile in temperature. The parameter s was not present in the original formalism of Guillot (2010) and Line et al. (2013) and we introduce it to control the slope of the temperature gradient in the atmosphere. To reduce the number of free parameters, we set t_int equal to zero, hence Eq. (1) becomes (5)

In Eq. (2), g_p is the planetary surface gravity computed considering the measured planetary mass and transit radius and, in Eq. (3), E₂ (γτ) is the second order exponential integral function of γτ. In Eq. (4), R_s is the stellar radius, a is the planetary orbital separation, and T_eff is the stellar effective temperature. Table 2 lists the relevant system parameters adopted in this work.

The TP profiles have been constructed considering the already available constraints. For setting the range of temperatures spanned by the models in the lower atmosphere, we considered the day- and night-side temperatures measured from phase curve observations(Mansfield et al. 2020; Wong et al. 2020). We further considered that the TP profile should present a temperature rise (Lothringer et al. 2018; Pino et al. 2020) and that the temperature in the PHOENIX TP profile at pressures below about 10 mbar may be underestimated (Turner et al. 2020). At the bottom of the atmosphere, around the 1 bar level, the empirical TP profiles havetemperatures ranging between about 3500 and 5500 K (i.e. about 1000 K above and below the measured day-side temperature), while at the top of the atmosphere the TP profiles have temperatures ranging between about 6000 and 11 000 K. We did not consider higher upper atmospheric temperatures because it would be unlikely for the planet to have an upper atmosphere hotter than the stellar photosphere (see Table 2), which is the source of heating (see also Mitani et al. 2020). Giant planets in close orbit to late-type stars present an upper atmospherictemperature hotter than the stellar photosphere, because the heating source is the radiation emitted by the stellar chromosphere and transition region, which are significantly hotter than the stellar photosphere. However, this is not the case for KELT-9b because the star does not possess a chromosphere and transition region (Fossati et al. 2018). Therefore, in the absence of an atmospheric heating mechanism other than photospheric stellar irradiation, an atmospheric temperature significantly higher than that of the stellar photosphere would be unlikely.

The 126 empirical TP profiles are divided into families having four different minimum temperatures at the bottom of the atmosphere (i.e. about 3500, 4000, 4750, and 5500 K) and six different maximum temperatures at the top of the atmosphere (i.e. about 6000, 7000, 8000, 9000, 10 000, and 11 000 K). In addition to the temperature minima and maxima, we varied also the pressure level at which the temperature rises and the slope of the temperature increase. We computed the TP profiles by setting κ equal to 2.5 or 3.1, γ ranging between 0.28 and 2.32 (81 different values), β ranging between 0.98 and 1.62 (34 different values), and s equal to 1, 2.5, or 4. For each TP profile, we set the maximum pressure (i.e. bottom of the atmosphere) at 2 bar and the minimum pressure (i.e. top of the atmosphere) at 8 × 10⁻¹² bar, and divided the atmosphere into 29 layers equally spaced in log p (i.e. steps of 0.408 in logp). This number of layers is a compromise between describing the TP profiles with enough accuracy and computation time. In this respect, for a few TP profiles, we run additional models with a larger number of layers (up to 35) obtaining the same results. This range of pressures is wide enough to fully contain the atmospheric formation region of the hydrogen Balmer lines (Turner et al. 2020) and ensures that the atmosphere is completely transparent to optical light at the top of the atmosphere and completely opaque to optical light at the bottom of the atmosphere. Figure 2 shows all TP profiles, also comparing them with that obtained with PHOENIX for KELT-9b, while Table A.1 lists the TP profiles and the parameters employed to obtain them. For reference, among our set of TP profiles, that computed by PHOENIX for KELT-9b is best reproduced by the TP profile number 025, which has κ equal to 3.1, γ equal to 0.99, β equal to 1.03, and s equal to 2.5. As mentioned above, the majority of the TP profiles is on average hotter than that computed by PHOENIX, because this is what was suggested by the analyses of García Muñoz & Schneider (2019), Turner et al. (2020), and Wyttenbach et al. (2020).

Fig. 2 Empirical TP profiles used as input for the spectral synthesis calculations. Black, red, blue, and green lines identify the models having a lower temperature, at the bottom of the atmosphere, close to 3500, 4000, 4750, and 5500 K, respectively. Solid lines are for profiles computed setting κ equal to 2.5, while dashed lines are for profiles computed setting κ equal to 3.1 The thick orange line shows the TP profile obtained for KELT-9b using PHOENIX.

4 Synthetic transmission spectra

To compute the synthetic Hα and Hβ transmission spectra on the basis of the TP profiles shown in Fig. 2, we employed the Cloudy spectral synthesis code following the scheme described in Turner et al. (2020) and Young et al. (2020). Cloudy is a spectral synthesis code designed to simulate physical conditions within an astrophysical plasma, predicting the emitted or absorbed spectrum, further accounting for NLTE effects (Ferland et al. 2017), which is important for correctly computing the population of atoms in excited states, such as those leading to the formation of the hydrogen Balmer lines (e.g. García Muñoz & Schneider 2019; Young et al. 2020). Cloudy further accounts for hydrogen and metal collisional excitation and ionisation, as well as photexcitation and photoionisation.

Since our empirical TP profiles are one-dimensional, we assume a spherically symmetric planetary atmosphere varying only with altitude. Cloudy is capable to consider all elements up to Zn and, to speed up the calculations,we took into account just H, He, C, O, Na, Mg, Si, K, Ca, Ti, and Fe, which are abundant and provide a significant contribution to the electron density through their ionisation, hence also to the mean molecular weight and continuum level. The only molecules we took into account are those of hydrogen, namely H₂, H, and H (see Young et al. 2020, for more details). Other molecules play a minor role in determining the atmospheric physical properties (Lothringer & Barman 2019), particularly around the formation region of the hydrogen Balmer lines. We took the atmospheric abundance profile of the considered elements from the PHOENIX model, which assumes solar metallicity, and let Cloudy compute dissociation,ionisation, and excitation.

In the atmospheric temperature profiles, the radial dimension is given by the pressure, but the Cloudy calculations require a physical length. Therefore, we computed the radius of each layer i using the pressure scale height (6)

where k_B is the Boltzmann constant, T_i is the temperature of layer i, μ_i is the mean molecular weight of layer i, and g_i is the planetary gravity at layer i. We then obtained the planetary radius corresponding to each layer by computing (7)

where p_i is the atmospheric pressure at layer i and p₀ is the reference pressure at the observed planetary radius R_p.

Because of the short planetary orbital distance and rather high stellar mass, the Roche lobe of KELT-9b is smaller than that of most hot Jupiters and it plays a significant role in shaping the atmospheric properties (e.g. Fossati et al. 2018). For this reason, when computing g_i in Eq. (6), we accounted for the shape and location of the Roche lobe in the direction perpendicular to the sub-stellar point, that is to say the one that was probed by transmission spectroscopy.

We mapped each one-dimensional (1D) TP profile onto concentric circles and then calculated the lengths through successive layers of atmosphere along line-of-sight transmission chords (see Turner et al. 2020; Young et al. 2020). These lengths, along with the atmospheric properties of their respective layers, have been stacked and entered into Cloudy as the line of sight transmission medium. We then computed separate transmission spectra with Cloudy for each layer, at a spectral resolution of R = 100 000 and without adding any turbulent velocity. Finally, we computed the total transmission spectrum of the planet by adding up the single layer spectra, weighted by their relative area projected on the stellar disc and accounting for the planetary impact parameter (b). When doing this last operation, we did not account for limb darkening effects because the observations of Yan & Henning (2018), Cauley et al. (2019), and Wyttenbach et al. (2020), which are those with the higher signal-to-noise ratio and hence more relevant for the analysis, are already corrected for the CLV effect. Further details on the algorithm employed to generate transmission spectra with Cloudy are given by Young et al. (2020). Because of the duration of Cloudy calculations, it is not possible to carry out this work accounting for the three dimensional geometry of the atmosphere, but we discuss the impact of the 1D assumption in Sect. 6.7. Overall, the broadening terms considered in computing the synthetic profiles are natural, temperature, Stark, and instrumental broadening. We run a few models testing the relative impact of the different broadening terms (excluding instrumental broadening) obtaining that thermal broadening is the most relevant, while Stark broadening is the least important. We do not consider rotational broadening, which would however have a small impact on the line profiles, compared in particular to thermal broadening.

Following what was indicated by the PHOENIX model, we first took p₀ to be equal to 0.01 bar (Turner et al. 2020). However, because the transmission spectra are normalised to the continuum, the choice of the reference pressure affects the line strength of the synthetic spectra. Furthermore, because each synthetic spectrum is computed employing a different TP profile, the reference pressure depends on the TP profile. To identify the reference pressure to employ to compute each transmission spectrum, we applied an iterative procedure. We first computed one transmission spectrum for each TP profile assuming p₀ = 0.01 bar. Then, for each transmission spectrum, we set p₀ equal to the pressure value of the layer for which the average continuum opacity calculated by Cloudy in the region covered by the Hα and Hβ lines is closest to 2/3 (the continuum opacity at the wavelengths covered by the Hα and Hβ lines is due to H⁻). Finally, we recomputed a new transmission spectrum with the updated reference pressure and followed this procedure again, until convergence. We reached convergence for all profiles following at most four iterations. The Cloudy code is not parallel, therefore we parallelised the computation of the different transmission rays. Employing 126 CPUs, split into two clusters to avoid memory overflow, we computed the transmission spectra, including the convergence to obtain the continuum level, within about 2.5 weeks.

Figure 3 shows, as an example, the continuum optical depth as a function of wavelength for each atmospheric layer computed with Cloudy in transmission geometry at the last iteration considering the TP model number 047. The continuum optical depth around the position of the Hα and Hβ lines is closest to 2/3 at the level corresponding to a pressure of 10 μbar.

5 Results

Before comparing the observed and synthetic line profiles, we normalised the synthetic spectra by fitting first order polynomials (one for the Hα and one for the Hβ line) to continuum points adjacent to and on both sides of each line. We also corrected the observations for wavelength shifts employing the line centres obtained from the Gaussian fits.

We compared the observed and synthetic profiles employing the χ² and the line equivalent widths as diagnostics. For the former, χ² has been computed employing the non-rebinned observed spectra and by interpolating the synthetic spectra on the sampling of each observation. For the latter, we compared the observed and synthetic equivalent widths and also the ratio of the Hα-to-Hβ equivalent widths. Both χ² and equivalent widths have been evaluated considering a range around the line centre of 2 Å for Hα and of 1.4 Å for Hβ. In the case of the Hα line there are 141 degrees of freedom, while in the case of the Hβ line there are 133 degrees of freedom.

We consider both χ² and equivalent widths, instead of just χ², to reduce the importance of reproducing line shape over line strength. Reproducing line profiles implies fitting both strength and shape of thelines, including line widths and shape of the wings, which is significantly more challenging than reproducing line strengths alone, that is equivalent widths. The widths and wings shape of spectral lines are influenced by several factors (e.g. turbulence) and physical processes (e.g. collisions), which are less likely to be correctly reproduced by the modelling. Instead, line strengths are more influenced by the general physical properties of the gas, such as temperature and densities, which are those we primarily aim to characterise.

Figure 4 shows the results obtained from the χ² analysis and from the comparison of the equivalent widths. The top-left panel of Fig. 4 indicates that there is a family of TP profiles that both minimises χ², reaching values of χ² between 100 and 400 (i.e. reduced χ² – – between about 1 and 3; see top-right panel of Fig. 4) for both lines, and matches well the measured equivalent widths. However, the bottom panels of Fig. 4 show that this family of models does not fit equally well the measured Hα-to-Hβ equivalent width ratio.

In general, because of the smaller error bars, the χ² obtained from the analysis of the weighted average profiles is larger than that obtained from the analysis of the single profiles. We also find that the profiles of Wyttenbach et al. (2020) are a better match to the synthetic lines. This is because the Hα and Hβ lines obtained by Wyttenbach et al. (2020) are shallower and broader than those of Cauley et al. (2019). We come back to this point later in this Section.

We identified the family of TP profiles best fitting the observations on the basis of the considered diagnostics by extracting those fulfilling the conditions (8)

where and are the values computed for the weighted average Hα and Hβ line profiles, respectively,|Δ EQW_Hα,WA| and |Δ EQW_Hβ,WA| are the absolute values of the difference between observed and synthetic equivalent widths considering the weighted average Hα and Hβ line profiles, respectively,and |Δ EQW_ratio,WA| is the absolute value of the difference between observed and modelled Hα-to-Hβ equivalent width ratio computed considering the weighted average line profiles. Under these conditions, the best fitting models correspond to those of the TP profiles number 030, 047, 054, 101, 108, and 125 that are shown in Fig. 5. We further looked for the three best fitting models by setting more stringent conditions, namely (9)

obtaining that the best fitting TP profiles are 047, 108, and 125, which are highlighted in Fig. 5 by thicker lines. The values given in the conditions listed in Eqs. (8) and (9) are arbitrary. This is because we just aim at identifying the general properties of the family of TP profiles leading to best fit the observations. In the next section, we then focus as an example on just one of those TP profiles, which has been chosen on the basis of considerations independent from the fitting of the hydrogen Balmer lines.

Figure 5 clearly shows that the observations are best reproduced by models based on TP profiles having an upper atmospheric temperature around 10 000–11 000 K, namely about 4000 K hotter than what predicted by the PHOENIX model and comparable to the stellar effective temperature. In general, synthetic lines computed on the basis of TP profiles with larger γ values are a better fit to the observations. Furthermore, the family of TP profiles best fitting the observations presents the temperature inversion in the middle of the atmosphere, roughly between 10 mbar and 1 μbar, indicating that the best fitting TP profiles are neither those with a very steep inversion at the bottom of the atmosphere, nor those with the inversion high up in the atmosphere. Figure 5 also suggests that the observed Hα and Hβ line profiles do not depend of the temperature of the lower atmosphere (i.e. ≳0.1 bar). However, the planetary temperature measured through phase curve observations can be at aid in further identifying the best fitting TP profile. As a matter of fact, at the photosphere (≈10 mbar level), the TP profiles number 047 and 125 are significantly hotter than the measured planetary day-side temperature, making TP model number 108 the most likely of the three (see also Lothringer & Barman 2020). Therefore, the TP profile number 108 is the one we employ to represent the family of TP profiles best fitting the observations and that we thoroughly discuss in the next section.

We performed this same analysis again, but considering the profiles of Cauley et al. (2019) and Wyttenbach et al. (2020) separately. We obtained that the result on the temperature of the upper atmosphere is robust as is does not depend on the chosen dataset, while the result on the shape and position of the temperature inversion depends slightly on the considered dataset, with the profiles obtained by Cauley et al. (2019) leading to a steeper temperature inversion located deeper in the atmosphere compared to what we obtained from the analysis of the profiles given by Wyttenbach et al. (2020).

Figure 6 shows the observed Hα and Hβ line profiles in comparison with the synthetic lines obtained employing the TP profiles number 108 (thick red solid line in Fig. 5), 025, which is the TP profile most resembling the PHOENIX one, and 114, which is the TP profile leading to the strongest Balmer lines. The best fitting profiles are slightly broader and shallower than the weighted averaged Hα and Hβ lines. The line depths of the best fitting synthetic profiles match well those of the profiles given by Wyttenbach et al. (2020), but are broader, while they are shallower and even broader than those of Cauley et al. (2019).

In agreement with Turner et al. (2020), the synthetic profile obtained employing a TP structure comparable to that given by PHOENIX leads to significantly weaker lines. The synthetic spectra presenting the strongest Hα and Hβ absorption features are only slightly deeper than those of Cauley et al. (2019), but are significantly broader, which is why both equivalent width and analyses do not favour these models.

In general, as shown in Fig. 7, the synthetic profiles that best fit the observed line depths, overestimate the observed widths. We remark that the broadening terms considered in computing the synthetic profiles are natural, temperature, Stark, and instrumental broadening, and that Stark broadening has a negligible impact on the total line broadening. It may be possible that the data analysis procedure, particularly the spectral normalisation, removes the signal coming from the far line wings, which the synthetic spectra suggest being not negligible, particularly for Hα, hence artificially reducing the width of the observed lines.

Fig. 3 Continuum optical depth as a function of wavelength in the optical range for each atmospheric layer in transmission geometry considering the TP model number 047. The values shown on the left-hand side give the logarithm of the pressure of the corresponding layer. For visualisation purposes, the pressure level is shown for every second layer starting from the top of the atmosphere and stop where the continuum optical depths of the layers start overlapping. To guide the eye, the red horizontal line is at 2/3, while the vertical blue dotted lines indicate the position of the Hα and Hβ lines.

Fig. 4 Top left: χ2 vs. difference between observed and modelled equivalent width (EQW) considering the profiles obtained by Cauley et al. (2019, black), Wyttenbach et al. (2020, orange), and for the weighted average (purple). Asterisks and rhombi are for the Hα and Hβ profiles, respectively. Bottom left: χ2 obtained from the analysis of the Hα line profiles vs. difference between observed and modelled Hα-to-Hβ equivalent width ratio. The error bars at the top-left (top panel) and bottom-left (bottom panel) corners indicate the size of the equivalent width and equivalent width ratio uncertainties. Right: same as left, but with on the x-axis. The dashed lines at zero and one are to guide the eye. In all panels, the position of the results obtained with the TP profiles number 025 and 108 is marked by red and blue circles, respectively.

Fig. 5 Same as Fig. 2, but for the TP profiles best fitting the observed Hα and Hβ lines according to the conditions listed in Eq. (8). The thicker lines indicate the three TP profiles fulfilling the stricterconditions listed in Eq. (9). The hatched area shows the main formation region of the Hα and Hβ lines according tothe three best fitting models.

6 Discussion

We analyse here in detail the results. We focus on the TP model profiles number 025, which is the closest to the one obtained with PHOENIX, and number 108, which, of the three models best fitting the planetary Hα and Hβ line profiles,is the one with the temperature of the lower atmosphere closer to what obtained from phase curve observations (Wong et al. 2020; Mansfield et al. 2020; Lothringer & Barman 2020).

6.1 Atmospheric composition

Figure 8 shows the details of the atmospheric hydrogen composition, namely the densities of neutral hydrogen (HI), protons (HII), H⁻, molecular hydrogen (H₂), H, H, and electrons (e⁻) with respect to the total hydrogen density, as a function of pressure for the TP models 025 and 108. For these calculations, we employed a geometry in which the atmosphere is illuminated from the top.

For the cooler TP profile, the atmosphere is dominated by neutral hydrogen at pressures higher than about 10⁻⁹ bar, while below that pressure protons are the most abundant species. For the hotter TP profile, instead, this threshold moves at a higher pressure of about 10⁻⁶ bar.

At the bottom of the atmosphere, at pressures higher than about 10⁻³ bar, the atmosphere is dominated by HI and H₂, with the latter decreasing rapidly with decreasing pressure. The third most abundant hydrogen species in the lower atmosphere is H⁻, but HII becomes quickly more abundant, at pressures of about 10⁻² bar and >1 bar for the cooler and hotter TP models, respectively. This is mostly due to thermal ionisation, because hydrogen photoionisation occurs at higher altitudes and it is believed to be small due to the shape of the stellar spectral energy distribution (Fossati et al. 2018). Despite this, at pressures higher than about 10⁻⁸ bar, the profiles obtained considering the hotter TP model indicate that H⁻ is on average about 10 times more abundant than what obtained considering the cooler TP model. This has an effect on the continuum level at optical and near infrared wavelengths that in ultra-hot Jupiters is controlled mostly by H⁻ and HI bound-free opacities (e.g. Arcangeli et al. 2018; Parmentier et al. 2018).

Figure 9 shows the mixing ratio for some of the species observationally most relevant as a function of pressure for the TP models number 025 and 108. Because of the assumption of a solar composition, HeI is the second most abundant species throughout most of the atmosphere. The distribution of Ca atoms is similar in the two cases, with CaI dominating below the ~0.1 bar level, CaII dominating up to the ~10⁻⁵ bar level, above which CaIII is the dominant Ca species. For the cooler TP model, Na is mostly neutral at pressures higher than 0.1 bar, while it is singly ionised throughout the whole atmosphere for the hotter TP model.

Because of their similar ionisation potentials, Mg and Fe present comparable behaviours, with FeI and MgI being the dominant species at pressures higher than 10–100 mbar, while FeII and MgII are more abundant at lower pressures. Interestingly, for the hotter TP model, at pressures below 10⁻⁶ bar there is a mixture of almost equally abundant FeII and FeIII. It may be therefore possible to further constrain the TP profile in KELT-9b by looking for the presence of FeIII lines in the planetary transmission spectrum. Some FeIII resonance lines are convenientlylocated in the near-ultraviolet (3200–3300 Å) and in the blue part of the optical spectrum (4000–5000 Å), hence where the stellar flux is high, enabling one to obtain high-quality spectra with a range of facilities such as HST, CUTE (Fleming et al. 2018), and high-resolution ground-based spectrographs.

Fig. 6 Top: same as Fig. 1, but the dashed black, red, and blue lines are the hydrogen Balmer line profiles obtained employing the TP profiles number 108, 025, and 114, respectively. The TP profile number 025 is the most similar to the PHOENIX self-consistent TP profile, the TP profile number 108 is the one leading to best fit the observed hydrogen Balmer lines, while the TP profile number 114 is the one leading to generate the strongest Balmer lines within our sample. The black, vertical dashed lines enclose the wavelength ranges considered to compute χ2 and the equivalent widths. Bottom: TP profiles corresponding to the 108 (black), 025 (red), and 114 (blue) model numbers.

Fig. 7 Difference between observed and synthetic line FWHM as a function of the difference between observed and synthetic line depth for the profiles obtained by Cauley et al. (2019, black), Wyttenbach et al. (2020, orange), and for the weighted average profile (purple). Asterisks and rhombi are for the Hα and Hβ lines, respectively.The dashed horizontal and vertical lines at 0 are to guide the eye. The models best fitting the observed line depths overestimate the lines FWHM, and vice versa. The error bars at the bottom-right corner indicate the size of the uncertainties on the measured line depth and FWHM. The position of the results obtained with the 025 and 108 TP profiles is marked by red and blue circles, respectively.

Fig. 8 Density relative to the total density of hydrogen for neutral hydrogen (HI; black solid), protons (HII; red), molecular hydrogen (H2; green), H (blue), H (violet), H− (orange), and electrons (e−; black dashed) computed for the TP model number 025 (top) and 108 (bottom).

Fig. 9 Mixing ratios for hydrogen (black), H2 (black-dotted), He (red), Ca (blue), Fe (green), Mg (violet), Na (orange), and electrons (brown) as a function of atmosphericpressure. Neutral (XI), singly ionised (XII), and doubly ionised (XIII) species are marked by solid, dashed, and dash-dotted lines, respectively. Top panel: TP model number 025. Bottom panel: TP model number 108.

6.2 Comparison with previous results

We compare here Cloudy mixing ratios obtained for the cooler TP model, which is comparable to that computed with PHOENIX, with the results of Kitzmann et al. (2018) and Lothringer et al. (2018). However, one has to keep in mind that our runs do not include molecules, except for H-bearing molecules, while Kitzmann et al. (2018) and Lothringer et al. (2018) consider a range of molecules. Kitzmann et al. (2018) employed an inverted TP profile with the upper atmospheric temperatureabout 1500 K cooler than that of our TP model number 025, and 1700 K cooler than that computed with PHOENIX. Furthermore, our calculations, as well as those of Kitzmann et al. (2018), account for both thermal ionisation and photoionisation, while the results of Lothringer et al. (2018) were obtained considering only thermal ionisation.

The Cloudy simulation indicates that at a pressure of about 2 bar, our upper pressure boundary, the mixing ratios of HI is higher than that of H₂ implying that the twomixing ratios become equal at even higher pressures and in particular higher than those obtained by both Lothringer et al. (2018, ~10⁻² bar) and Kitzmann et al. (2018, ~2 bar). The Fe mixing ratios computed by Cloudy are very similar to those obtained by Kitzmann et al. (2018), who also found that FeII becomes the dominant species at pressures lower than 10 mbar. We obtain a good match also for the Na mixing ratio with what presented by Lothringer et al. (2018). A significant difference is found instead for the electron mixing ratio in the upper part of the atmosphere, where Cloudy gives a ~ 100 times larger electron density than that given by Lothringer et al. (2018). We ascribe this to the fact that Cloudy considers photoionisation, while PHOENIX does not (Lothringer et al. 2018). There is also a difference in the H⁻ abundance between what we obtained with Cloudy and what presented by Lothringer et al. (2018). While their H⁻ mixing ratio decreases almost monotonically with decreasing pressures, our H⁻ mixing ratio increases up to 10⁻⁴–10⁻⁵ bar to then decrease monotonically at lower pressures (see Fig. 8). Since at high pressures the electron density can be considered to be consistent between the two computations, we ascribe this difference to different implementations of the physics controlling H⁻. It will be important in the future to perform a detailed comparison, because of the importance of H⁻ in controlling the planetary atmospheric continuum, and hence the predicted strength of spectral lines.

6.3 Relevance of NLTE effects

We took advantage of the NLTE radiative transfer capabilities of Cloudy to study the presence and impact of NLTE effects in the modelling of the Hα and Hβ lines. To this end, we ran Cloudy for the TP models 025 and 108 employing a geometry in which the atmosphere is illuminated from the top, extracting then the NLTE hydrogen level populations and the Hα and Hβ optical depth as a function of atmospheric pressure. We further used the TP models to extract the ratio between the densities of hydrogen in the n = 2 and n = 1 levels employing the Boltzmann equation, hence assuming LTE. Figure 10 shows the results of this analysis.

A comparison of the top- and bottom-left panels of Fig. 10 indicates that the hotter TP model indeed leads to a larger density of excited hydrogen atoms by a factor of 10³ –10⁴ compared to what obtained from the cooler TP model, hence stronger Balmer lines better fitting the observations. Figure 10 also shows that, except for the case of the hotter TP profile and at pressures below ~10⁻⁶ bar, there seemsto be a general good match between the LTE and NLTE n = 2 and n = 1 hydrogen level populations, possibly indicating that the LTE approximation might be valid for computing the hydrogen Balmer lines. However, the profiles obtained from TP model number 025 indicate that at pressures corresponding to the main Hα and Hβ line formation region (i.e. between 10⁻⁵ and 10⁻⁸ bar; right) the n = 2 LTE level population is almost 10 times smaller than that computed accounting for NLTE effects, hence leading to LTE underestimating the strength of the lines. For the hotter TP model, the LTE and NLTE level populations at pressures below ~10⁻⁶ bar diverge significantly, leading to LTE overestimating the strength of the lines. We also note that, considering an isothermic profile with a temperaturesimilar to that of the upper atmosphere of TP model number 108, Wyttenbach et al. (2020) found an even larger difference (about 10⁶) in the LTE population of excited hydrogen compared to NLTE, but their result is based on constant temperature and density ratio profiles that are further strongly correlated with the mass-loss rate, which hamper comparisons between the two results.

By construction, Cloudy is unable to perform calculations of line profiles assuming LTE for the n = 1 and n = 2 hydrogen levels. For this reason, we computed the Hα transmission line profiles assuming LTE approximation as follows. First, we computed the LTE number densities of excited (i.e. n = 2) hydrogen using the Boltzmann formula, further considering the lowering of the ionisation potential due to Debye screening by electrons and ions in the calculation of the partition function. Next, we assumed the absorption profile to be described by a Voigt function and considered temperature (Doppler), natural, and pressure (Stark and van der Waals) broadening. Because of the relatively low gas density in the region of the atmosphere probed by transmission spectroscopy, pressure broadening is weak, in agreement with our Cloudy simulations and the results of Wyttenbach et al. (2020). We additionally cross checked our calculations by using a more elaborate treatment of the Hα profile employing pressure broadening tables after Lemke (1997) obtaining a good match between the two approaches. This confirms that a Voigt profile is a good approximation for interpreting transmission observations of the Hα line. All relevant numerical routines were taken from the LLMODELS stellar model atmosphere code (Shulyak et al. 2004). To calculate the LTE transmission profiles we employed the τ-REx (Tau Retrieval for Exoplanets) forward model (Waldmann et al. 2015a). We used an extended version of the original package, which now includes additional continuum opacity sources relevant for the high atmospheric temperatures of KELT-9b, as described in Shulyak et al. (2020).

Because LLMODELS is a module structured code, we turned some of its modules into numerical libraries and loaded them into τ-REx to compute continuum opacity additional to those described in Shulyak et al. (2020). In particular, we calculated opacities due to HI, He, and various metals, which contribute a small, but non-negligible, ~8% to the total transmission depth at the wavelengths around the Hα line.

The number densities of neutral and ionised species were taken from the Cloudy calculations for the respective model structures. This ensures that the background opacity is as similar as possible between the τ-REx and Cloudy computations, enabling us to isolate and study the effect of NLTE on the Hα transmission spectra.

Figure 11 presents the comparison of the LTE and NLTE Hα transmission spectra computed considering the TP model number 025 and 108. As expected from the comparison of the level populations, the LTE Hα profiles obtained considering the TP model number 025 is comparable to that obtained considering NLTE effects. Instead, with the hotter model, the LTE Hα transmission spectrum is stronger than the one obtained considering NLTE, particularly for the line wings, though the line core appears to be a good fit, in a way similar to our model producing the strongest lines (TP model number 114; see Fig. 6). This difference is caused by the difference in the LTE and NLTE level populations at pressures lower than 10⁻⁶ bar.

In agreement with García Muñoz & Schneider (2019), our simulations indicate that NLTE effects are extremely important for correctlymodelling the hydrogen Balmer lines in the transmission spectrum of KELT-9b. Although Wyttenbach et al. (2020) do not find strong evidence in favour of NLTE, they find a departure from LTE of the order of 10⁶ and assumed a constant density ratio of excited hydrogen to total hydrogen throughout the atmosphere, which our simulations clearly indicate being not realistic, particularly for the hotter models (i.e. temperatures similar to those obtained by Wyttenbach et al. 2020). Indeed, contrary to Wyttenbach et al. (2020), who considered only a simple scaling for the hydrogen n = 2 level population, we carried out realistic and self-consistent NLTE calculations taking into account the relevant radiative and collisional transitions between bound and free states in the hydrogen atom.

In an attempt to more meaningfully compare our LTE calculations with those of Wyttenbach et al. (2020), we computed an LTE Hα profile assuming an isothermal temperature at T = 13 000 K, which is theLTE temperature profile obtained by Wyttenbach et al. (2020) best fitting the observed Hα line, and the chemicalabundances obtained from a Cloudy run employing the same isothermal profile. However, we find that the depth of our LTE Hα profile constructed under these assumptions is about twice that of Wyttenbach et al. (2020) and that it would require to decrease thedensity of excited hydrogen by a factor of ten to be able to fit the observations. It is possible that this further discrepancy is caused by differences in the underlying chemical composition, that is Cloudy vs. equilibrium chemistry without photoionisation. It is also possible that meaningful comparisons of our results with those of Wyttenbach et al. (2020) are in general not possible, because they considered the planetary mass-loss rate as a free parameter, further obtaining that it strongly correlates with the atmospheric temperature and hydrogen density profiles.

Fig. 10 Top left: ratio of the n = 2 and n = 1 level populations for neutral hydrogen (bottom x-axis) as a function of atmospheric pressure (left y-axis) and planetary radius (right y-axis) computedwith Cloudy in NLTE (solid line) and using the Boltzmann equation (i.e. LTE; dashed line) considering the TP model number 025, shown by the red line (red top x-axis). Top right: line optical depth as a function of atmospheric pressure for Hα (solid line) and Hβ (dashed line) computed with Cloudy considering the TP model number 025. Both hydrogen lines form mostly between 10−5 and 10−8 bar. Bottom: same as top panels, but for the TP model number 108.

Fig. 11 Comparison between observed and synthetic Hα line profiles computed in LTE (red dashed line) and NLTE (black dashed line) considering the TP model numbers 025 (left) and 108 (right). The colour scheme of the observed line profiles is the same as in Fig. 1.

6.4 Impact of stellar effective temperature

So far, we attempted to fit the observed hydrogen Balmer lines by modifying the TP profile, obtaining that it should be significantly different than that obtained self-consistently with PHOENIX. However, there are other ways of modifying the line profiles without necessarily changing significantly the TP profile, for example by decreasing the mean molecular weight, hence increasing the pressure scale height, and/or by increasing the temperature of the host star, hence increasing the atmospheric energy available for exciting hydrogen. Therefore, we discuss whether increasing the effective temperature of the planetary atmosphere (this section) or decreasing the metallicity of the planetary atmosphere (next section) would enable one to fit the observed hydrogen Balmer lines employing the self-consistently computed TP profile. Furthermore, looking for the possible effect of a different stellar spectral energy distribution is important at the light of the gravity darkening effect recently detected for KELT-9, with the temperature ranging from 10 200 K at the poles and 9400 K at the equator, and with the stellar inclination angle tilted by almost 40 degrees (Ahlers et al. 2020).

As shown by Lothringer & Barman (2019), the shape of the input stellar spectral energy distribution has an impact on the atmospheric properties of ultra-hot Jupiters (see also Fossati et al. 2018). Therefore, we checked whether a variation in the temperature of the star compatible with what given by stellar spectroscopic analyses leads to significant changes in the strength of the Hα and Hβ line profiles.

To this end, we calculated a further transmission spectrum with Cloudy considering the TP model number 025 and a stellar spectral energy distribution computed with an effective temperature of 10 200 K, hence 600 K hotter than what employed for the analysis presented in Sect. 5. We remark that the stellar effective temperature we adopted in this work is that given by Borsa et al. (2019), while the hotter one we chose for studying the impact of stellar effective temperature on the transmission spectrum is that given by Gaudi et al. (2017).

Figure 12 shows the comparison between the Hα and Hβ line profiles obtained for the two stellar spectral energy distributions and considering the TP model number 025. The hotter star leads to a slight increase in the level population of excited hydrogen and therefore to stronger lines, but the increase is small compared to what would be necessary to fit the observed profiles. This indicates that fitting the observed line profiles requires a planetary atmosphere with a temperature as high as ~10 000–11 000 K, even when considering a hotter star. This result further indicates that the fact that the planet crosses regions of the star with different temperatures due to the fast stellar rotation (Ahlers et al. 2020) has no detectable effect on the transmission spectra.

Fig. 12 Comparison between the Hβ (left) and Hα (right) line profiles computed in NLTE with Cloudy employing an incident stellar flux corresponding to a star with an effective temperature of 9600 K (black; Borsa et al. 2019) and 10 200 K (red; Gaudi et al. 2017). To ease the comparison, the range of the y-axis is the same as that of the top panels of Fig. 6.

6.5 Planetary atmospheric metallicity

Following Eq. (6), it is possible to increase the pressure scale height and hence the size of the absorption features, without exceedingly increasing the temperature, by decreasing the mean molecular weight. Therefore, we computed a synthetic transmission spectrum covering the Hα and Hβ line profiles, employing the TP profile number 025 and a metallicity of 0.1 × solar.

Figure 13 shows the comparison between the Hα and Hβ line profiles obtained for the TP model number 025 and considering a solar-like and sub-solar planetary atmosphericmetallicity. Decreasing the metallicity only very slightly increased the depth of the absorption lines furtherconfirming that fitting the observations requires a planetary atmosphere with a temperature of ~10 000–11 000 K.

Fig. 13 Comparison between the Hβ (left panel) and Hα (right panel) transmission profiles obtained considering the TP profile number 025 and a solar metallicity (Z = 1; black line) or a sub-solar metallicity (Z = 0.1; red line). To ease the comparison, the range of the y-axis is the same as that of the top panels of Fig. 6.

6.6 Possible origin and consequences of the additional atmospheric heating

Our results indicate that the atmosphere of KELT-9b, particularly at pressures lower than ~10⁻³ bar, should be significantly hotter than that estimated by self-consistent calculations, for example with PHOENIX. However, the difficulty of simultaneously fitting the strength and the width of the Hα and Hβ line profiles suggests that either our results slightly overestimate the upper atmospheric temperature or the analysis of the observations, possibly the spectral normalisation, partially removes the outer edges of the line wings. We find that the upper atmospheric temperature of the family of TP profiles leading to best fit the observations is around 10 000–11 000 K that is comparable to that of the stellar photosphere, which is the main heating source. This temperature also reminds of the upper atmospheric temperature of classical hot Jupiters orbiting late-type stars (e.g. Yelle 2004). However, there are significant differences in the heating source and mechanisms between classical hot Jupiters and ultra-hot Jupiters like KELT-9b orbiting early A-type stars (Fossati et al. 2018).

The upper atmospheres of planets orbiting stars later than spectral type A3-A4 are heated by absorption of the stellar high-energy emission, namely X-ray and EUV (XUV) radiation. The temperature reached in the upper atmosphere of these planets is typically of the order of 10⁴ K, which is also of the same order of magnitude as that of the plasma in the stellar chromosphere and transition region emitting the majority of the XUV radiation absorbed by the planet and causing the heating. Planets orbiting early A-type stars, instead, are believed not to possess a chromosphere and transition region, hence not to emit significant XUV radiation¹ (Fossati et al. 2018). Indeed, the atmospheres of ultra-hot Jupiters are heated primarily by the less energetic UV and optical stellar photons. Given that the stellar photospheric temperature is of the order of 10 000 K, from a thermodynamic point of view it would be possible for the planetary atmosphere to reach such values. For the same reason, it is then unlikely for the planetary atmosphere to become much hotter, unless there is significant atmospheric heating due to other mechanisms, such as ohmic dissipation or magnetohydrodynamic waves (e.g. Thorngren & Fortney 2018). This conclusion is also supported by detailed hydrodynamic simulations of the upper atmosphere of KELT-9b conducted by Mitani et al. (2020).

It is, however, interesting to notice that the upper atmospheric temperature of the TP profile we consider as reference for leading to best fit the observations lies at the upper limit of the explored temperature range. Therefore, we generated 30 more TP profiles with a lower atmospheric temperature close to that of the TP profile number 108 and upper atmospheric temperatures ranging between 11 000 and 16 000 K. We further produced synthetic transmission spectra on the basis of these profiles and compared them with the observations in the same way as described in Sects. 4 and 5. We obtained that four of these hotter models satisfied the conditions set in Eq. (8) and that three models satisfied the more restrictive conditions set in Eq. (9). However, all these models display an upper atmospheric temperature in excess of 12 000 K, which, given the stellar effective temperature, we consider not being realistic for describing the TP profile of KELT-9b. Furthermore, none of these models leads to a better fit to the data with respect to the TP profile number 108 and also have the problem of being shallower and broader than the observations.

It is therefore possible that self-consistent models computing the planetary TP structure underestimate atmospheric heating and/or overestimate cooling. García Muñoz & Schneider (2019) showed that, because of the shape of the stellar spectral energy distribution and of the large density of excited hydrogen in the planetary atmosphere, heating through absorption of the stellar UV and optical radiation plays an important role.

In Sect. 6.1, we showed that the planetary atmospheric composition we obtained from Cloudy differs from those of Kitzmann et al. (2018) and Lothringer et al. (2018) also in the electron density, which we find being significantly higher, particularly in the upper atmosphere. Indeed, the Cloudy runs suggest a higher density of ionised metals and a higher metal ionisation degree compared to previous models, particularly in the upper atmosphericlayers. This may be due to the fact that, in contrast to PHOENIX, Cloudy includes metal photoionisation. Given that metal absorption and ionisation are among the main heating mechanisms in ultra-hot Jupiters, it is possible that self-consistent calculations of the atmospheric TP profile accounting for photoionisation would lead to metal ionisation properties similar to those found with Cloudy, hence to increased upper atmospheric temperatures.

An upper atmosphere at a temperature of 10 000–11 000 K would also explain the strong CrII and FeII absorption features obtained by Hoeijmakers et al. (2019) from the analysis of optical high-resolution transmission spectroscopy observations. The observed features were surprisingly large compared to the model, which was produced considering an isothermal temperature profile at 4000 K, while a higher atmospheric temperature leads to a stronger Fe and Cr ionisation, hence to larger features of both species. Therefore, it would be important to re-analyse these observations considering a model computed on the basis of a hotter TP profile (e.g. the TP model number 108), which would lead to increase the detection level of several ions and possibly also to the detection of additional elements that would constrain the shape of the atmospheric TP profile (e.g. FeIII).

6.7 Impact of the model assumptions

By employing Cloudy, we were able to lift the LTE approximation for modelling lines rising from excited hydrogen levels and therefore run all computations accounting for NLTE effects for all considered elements. Despite this, the synthetic hydrogen Balmer line profiles did not fit well enough the observations and therefore we discuss here the modelling assumptions that may possibly play a role in this outcome, focusing in particular on the assumption of a 1D geometry and of a hydrostatic atmosphere.

Because of the 1D approximation, we assumed that the day and night side of the planetary atmosphere through which the stellar light travels have the same properties. This is clearly not the case, though global circulation models of ultra-hot Jupiters suggest that there is a relatively small temperature (i.e. 1000–1500 K) and compositional difference between the day and the night side (e.g. Bell & Cowan 2018; Tan & Komacek 2019) and phase curve observations of KELT-9b confirm this prediction (Wong et al. 2020; Mansfield et al. 2020), hence alleviating the effect of this assumption on our results. The 1D assumption likely leads one to overestimate the number of excited hydrogen atoms, but also to overestimate hydrogen ionisation. These two effects counteract each other in terms of their impact on the strength of the hydrogen Balmer lines and future simulations accounting for the three-dimensional nature of the planetary atmosphere are therefore necessary to assess how the 1D assumption affects synthetic transmission spectra of KELT-9b.

The assumption of a hydrostatic atmosphere may have two important effects on our results: (1) the atmospheric layers characterised by strong upward acceleration (i.e. in a hydrodynamic state) have a larger pressure scale height, because the upward acceleration reduces locally the planetary gravity; (2) a hydrodynamic atmosphere is also characterised by a significant upward transport of gas from the bottom layers, which would for example lead to increase the density of (excited) neutral hydrogen atoms higher up in the atmosphere. Both these effects would lead to an increase in the strength of the Hα and Hβ lines, without the need of increasing the atmospheric temperature. Fossati et al. (2018), García Muñoz & Schneider (2019), and Mitani et al. (2020) showed that the upper atmosphere of KELT-9b is probably expanding hydrodynamically.

Within Cloudy, to be able to fix the temperature profile, we had to exclude from the calculations all non-hydrogen-based molecules, which may have affected the lower atmosphere. However, our results indicate that the lower atmosphere has a negligible impact on the shape and strength of the hydrogen Balmer lines. Furthermore, the overall hotter atmospheric temperature would lead to an even lower molecular abundance, decreasing the impact these would have on the results.

7 Conclusion

We constructed 126 empirical atmospheric temperature-pressure profiles for the ultra-hot Jupiter KELT-9b varying the lower and upper atmospheric temperatures, and the location and gradient of the temperature rise, further considering the available observational constraints (Wong et al. 2020; Mansfield et al. 2020; Pino et al. 2020; Turner et al. 2020). We then employed the Cloudy NLTE radiative transfer code in a one-dimensional geometry to produce transmission spectra of the Hα and Hβ lines on the basis of the constructed TP profiles, comparing them with the available observations.

We found that the family of TP profiles leading to best fit the observations is characterised by an upper atmospheric temperature of 10 000–11 000 K, hence about 4000 K hotter than predicted by the PHOENIX model, and by a temperature rise starting at pressures higher than 10⁻⁴ bar, which is around the highest pressure level probed by the Hα and Hβ lines. For the TPprofile leading to best fit the observations, we compared the n = 2 excited hydrogen level population and the Hα and Hβ line profiles computed in LTE and NLTE obtaining that the LTE approximation leads one to overestimate the level population by several orders of magnitude, and hence to overestimate the strength of the lines. Furthermore, the Cloudy simulations also clearly indicate that other commonly considered assumptions, such as constant temperature and/or density ratio profiles, are not realistic and would lead to misinterpret the data.

We further analysed the atmospheric chemical composition obtained with Cloudy employing the TP profile leading to best fit the Hα and Hβ lines. The main results of this analysis are that the high upper atmospheric temperature is likely caused by metal photoionisation and that at pressures lower than about 10⁻⁶ bar, FeIII is almost as abundant as FeII. The latter finding indicates that FeIII may be also detectable in the transmission spectra of KELT-9b. Its identification in the observed transmission spectra would be a confirmation that the upper atmospheric temperature of KELT-9b is significantly higher than previously thought. More in general, the detection and modelling of features rising from different ions of the same element and from different electronic levels of the same ion is important to constrain the physical properties of the planetary atmosphere, including the local thermodynamic state, which is a key piece of information for the further development of planetary atmosphere models. However, the detection of such metal features should be based on synthetic spectra computed with adequate TP profiles and accounting for NLTE effects that become relevant when looking at lines rising from excited states, as is the case of the hydrogen Balmer lines.

The family of TP profiles leading to best fit the observed hydrogen Balmer lines, however, leads still to profiles that are shallower and broader than the observations, though the line strength (i.e. equivalent width) is comparable to the measured one. There are a range of possible explanations for this result connected to both data analysis and modelling assumptions. The observed transmission spectra covering the hydrogen Balmer lines of KELT-9b have been obtained with different instruments and employing slightly different data analysis methods, which may have introduced systematic differences among the results. It is also not possible to exclude that the data analysis procedures have removed part of the far line wings, leading to narrower spectral lines. Therefore, we suggest that a uniform analysis of all available spectroscopic data covering the primary transit of KELT-9b would be beneficial for future theoretical studies aiming at fitting the observations.

On the modelling side, having lifted the LTE approximation, the two main assumptions affecting the results are those of a 1D geometry and of a hydrostatic atmosphere. The impact of the former on the synthetic spectra is unclear, because of the competing effects (lower hydrogen ionisation and lower hydrogen excitation) that the cooler nightside temperature would have on the transmission spectra. Instead, the implementation of a model accounting for hydrodynamic effects may change the shape of the line profiles precisely in the desired direction, namely making them deeper without the need of further increasing the temperature, which would further broaden the line wings. Therefore, we suggest that future modelling and fitting of the hydrogen Balmer lines of KELT-9b shall attempt to account for hydrodynamic effects in the upper atmosphere, further to accounting for NLTE effects.

Therefore, future observational works should primarily attempt to understand the origin of the differences detected in the currently available observations, particularly by collecting higher-quality data, if possible, and by analysing them in a homogeneous way. However, unfortunately, given the significant computing time currently necessary for producing synthetic spectra accounting for NLTE effects, at least for the near future, fitting the observations will be possible only through forward modelling.

An upper atmospheric temperature of 10 000–11 000 K has strong implications for the planetary atmospheric mass loss. Indeed, the extra energy available as a result of such a high atmospheric temperature would lead to a significantly higher mass-loss rate compared to the 10¹¹ g s⁻¹ obtained on the basis of the PHOENIX temperature-pressure profile (Fossati et al. 2018). A significantly higher mass-loss rate may have important consequences for the atmospheric evolution of the planet, but may also indicate that the mass-loss rates expected also for other ultra-hot Jupiters, particularly those orbiting hotter stars, might need to be reconsidered. Indeed, the analyses of Yan & Henning (2018) and Wyttenbach et al. (2020) suggested for KELT-9b mass-loss rates as high as 10¹² g s⁻¹. Interestingly, the ultra-hot Jupiters detected so far appear to be on average more massive than the classical hot Jupiters of the same radius, which is not a bias effect because all these planets have been detected with the transit method and have similar radii. Therefore, planets with masses smaller than those typical of ultra-hot Jupiters may be difficult to detect because of their very small radius as a consequence of the strong atmospheric escape. Future work shouldtherefore also aim to constrain and estimate planetary mass-loss rates for ultra-hot Jupiters as a function of planetary parameters, to identify whether there could be a signature of mass loss in the observed planet population.

Acknowledgements

M.E.Y. acknowledges funding from the ÖAW-Innovationsfonds IF_2017_03. D.S. and M.R. acknowledge the support of the DFG priority program SPP-1992 “Exploring the Diversity of Extrasolar Planets” (DFG PR 36 24 602/41). A.G.S. and L.F. acknowledge financial support from the Austrian Forschungsförderungsgesellschaft FFG project CONTROL P865968. L.M.L. acknowledges the financial support from the State Agency for Research of the Spanish MCI through the “Center of Excellence Severo Ochoa” award SEV-2017-0709, and from the research project PGC2018-099425-B-I00. A.W. acknowledges the financialsupport of the SNSF by grant number P400P2_186765. The authors thank the anonymous referee for the insightful comments.

Appendix A Empirical TP profiles and their parameters

References

All Tables

All Figures

Fig. 1 Transmission spectra of the hydrogen Balmer line profiles (Hβ, left; Hα, right) published by Yan & Henning (2018, green), Cauley et al. (2019, black), Turner et al. (2020, blue), and Wyttenbach et al. (2020, orange). The magenta line is the weighted average profile obtained considering just the profiles of Cauley et al. (2019) and Wyttenbach et al. (2020). Wavelengths are in vacuum to match those of the synthetic transmission spectra (see Sect. 4). The data have been rebinned by a factor of five for visualisation purposes. The dashed lines are Gaussian fits to the data. The horizontal dotted line at one sets the continuum level to guide the eye.

In the text

Fig. 2 Empirical TP profiles used as input for the spectral synthesis calculations. Black, red, blue, and green lines identify the models having a lower temperature, at the bottom of the atmosphere, close to 3500, 4000, 4750, and 5500 K, respectively. Solid lines are for profiles computed setting κ equal to 2.5, while dashed lines are for profiles computed setting κ equal to 3.1 The thick orange line shows the TP profile obtained for KELT-9b using PHOENIX.

In the text

Fig. 3 Continuum optical depth as a function of wavelength in the optical range for each atmospheric layer in transmission geometry considering the TP model number 047. The values shown on the left-hand side give the logarithm of the pressure of the corresponding layer. For visualisation purposes, the pressure level is shown for every second layer starting from the top of the atmosphere and stop where the continuum optical depths of the layers start overlapping. To guide the eye, the red horizontal line is at 2/3, while the vertical blue dotted lines indicate the position of the Hα and Hβ lines.

In the text

Fig. 4 Top left: χ2 vs. difference between observed and modelled equivalent width (EQW) considering the profiles obtained by Cauley et al. (2019, black), Wyttenbach et al. (2020, orange), and for the weighted average (purple). Asterisks and rhombi are for the Hα and Hβ profiles, respectively. Bottom left: χ2 obtained from the analysis of the Hα line profiles vs. difference between observed and modelled Hα-to-Hβ equivalent width ratio. The error bars at the top-left (top panel) and bottom-left (bottom panel) corners indicate the size of the equivalent width and equivalent width ratio uncertainties. Right: same as left, but with on the x-axis. The dashed lines at zero and one are to guide the eye. In all panels, the position of the results obtained with the TP profiles number 025 and 108 is marked by red and blue circles, respectively.

In the text

	Fig. 5 Same as Fig. 2, but for the TP profiles best fitting the observed Hα and Hβ lines according to the conditions listed in Eq. (8). The thicker lines indicate the three TP profiles fulfilling the stricterconditions listed in Eq. (9). The hatched area shows the main formation region of the Hα and Hβ lines according tothe three best fitting models.
In the text

Fig. 6 Top: same as Fig. 1, but the dashed black, red, and blue lines are the hydrogen Balmer line profiles obtained employing the TP profiles number 108, 025, and 114, respectively. The TP profile number 025 is the most similar to the PHOENIX self-consistent TP profile, the TP profile number 108 is the one leading to best fit the observed hydrogen Balmer lines, while the TP profile number 114 is the one leading to generate the strongest Balmer lines within our sample. The black, vertical dashed lines enclose the wavelength ranges considered to compute χ2 and the equivalent widths. Bottom: TP profiles corresponding to the 108 (black), 025 (red), and 114 (blue) model numbers.

In the text

Fig. 7 Difference between observed and synthetic line FWHM as a function of the difference between observed and synthetic line depth for the profiles obtained by Cauley et al. (2019, black), Wyttenbach et al. (2020, orange), and for the weighted average profile (purple). Asterisks and rhombi are for the Hα and Hβ lines, respectively.The dashed horizontal and vertical lines at 0 are to guide the eye. The models best fitting the observed line depths overestimate the lines FWHM, and vice versa. The error bars at the bottom-right corner indicate the size of the uncertainties on the measured line depth and FWHM. The position of the results obtained with the 025 and 108 TP profiles is marked by red and blue circles, respectively.

In the text

	Fig. 8 Density relative to the total density of hydrogen for neutral hydrogen (HI; black solid), protons (HII; red), molecular hydrogen (H2; green), H (blue), H (violet), H− (orange), and electrons (e−; black dashed) computed for the TP model number 025 (top) and 108 (bottom).
In the text

	Fig. 9 Mixing ratios for hydrogen (black), H2 (black-dotted), He (red), Ca (blue), Fe (green), Mg (violet), Na (orange), and electrons (brown) as a function of atmosphericpressure. Neutral (XI), singly ionised (XII), and doubly ionised (XIII) species are marked by solid, dashed, and dash-dotted lines, respectively. Top panel: TP model number 025. Bottom panel: TP model number 108.
In the text

Fig. 10 Top left: ratio of the n = 2 and n = 1 level populations for neutral hydrogen (bottom x-axis) as a function of atmospheric pressure (left y-axis) and planetary radius (right y-axis) computedwith Cloudy in NLTE (solid line) and using the Boltzmann equation (i.e. LTE; dashed line) considering the TP model number 025, shown by the red line (red top x-axis). Top right: line optical depth as a function of atmospheric pressure for Hα (solid line) and Hβ (dashed line) computed with Cloudy considering the TP model number 025. Both hydrogen lines form mostly between 10−5 and 10−8 bar. Bottom: same as top panels, but for the TP model number 108.

In the text

	Fig. 11 Comparison between observed and synthetic Hα line profiles computed in LTE (red dashed line) and NLTE (black dashed line) considering the TP model numbers 025 (left) and 108 (right). The colour scheme of the observed line profiles is the same as in Fig. 1.
In the text

	Fig. 12 Comparison between the Hβ (left) and Hα (right) line profiles computed in NLTE with Cloudy employing an incident stellar flux corresponding to a star with an effective temperature of 9600 K (black; Borsa et al. 2019) and 10 200 K (red; Gaudi et al. 2017). To ease the comparison, the range of the y-axis is the same as that of the top panels of Fig. 6.
In the text

	Fig. 13 Comparison between the Hβ (left panel) and Hα (right panel) transmission profiles obtained considering the TP profile number 025 and a solar metallicity (Z = 1; black line) or a sub-solar metallicity (Z = 0.1; red line). To ease the comparison, the range of the y-axis is the same as that of the top panels of Fig. 6.
In the text