© The Authors 2025

1 Introduction

Ultra-hot Jupiters (UHJs) are close-in gas giant planets with equilibrium temperatures of T_eq ≳ 2000 K. The high temperature favours metal ionisation in the planetary atmosphere and leads to the production of a significant amount of H⁻, which dominates the continuum of the planetary transmission and emission spectra (e.g. Arcangeli et al. 2018; Parmentier et al. 2018). Thanks to the relatively large atmospheric pressure scale heights driven by the high equilibrium temperatures and the brightness of their host stars, UHJs have become prime targets for atmospheric observations conducted at both transmission and emission geometry. This has fostered a significant modelling effort whose aim is to interpret the observations.

Several UHJs orbit around relatively bright F- and A-type stars, which has driven the community to study their characteristics in greater detail (e.g. Ahlers et al. 2020; Kama et al. 2023). Among the properties of planet-hosting stars, the high-energy emission – X-ray and extreme ultraviolet (EUV; together XUV) – is one of the most important because it has a profound impact on the planetary atmospheric properties and evolution (e.g. Lammer et al. 2003; Yelle 2004; Owen & Wu 2017; Kubyshkina et al. 2018). Observations have suggested that intermediate-mass stars cooler than about 8250 K have strong XUV emission, possibly significantly stronger than solar, while the XUV emission of hotter stars can be approximated as black-body emission and is thus relatively low (e.g. Fossati et al. 2018; Günther et al. 2022). Therefore, we expect that UHJs orbiting intermediate-mass stars cooler than about 8250 K have a very hot and extended upper atmosphere with hydrogen ionisation driven by the incident XUV stellar irradiation, while the atmospheres of planets orbiting hotter intermediate-mass stars are mostly controlled by the stellar ultraviolet (UV) and optical incident radiation (e.g. Fossati et al. 2018, 2021, 2023; Lothringer & Barman 2019; García Muñoz & Schneider 2019).

Among the currently known UHJs, KELT-9b has the highest equilibrium temperature (T_eq ≈ 4000 K; Gaudi et al. 2017; Jones et al. 2022), and thus it is one of the most studied UHJs to date. KELT-9b has been observed from the ground using different high-resolution spectrographs to study the emission and transmission spectra at optical wavelengths. Both hydrogen and a range of neutral and ionised metal species have been detected in the planetary atmosphere through cross-correlation and/or single line profile analyses (e.g. Yan & Henning 2018; Hoeijmakers et al. 2018; Turner et al. 2020; Wyttenbach et al. 2020; Pino et al. 2020; Sánchez-López et al. 2022; Borsato et al. 2023; D’Arpa et al. 2024; Stangret et al. 2024). Eclipse and phase curve observations of this planet have enabled constraints to be placed on the albedo, the lower atmospheric temperature, and atmospheric dynamics (e.g. Hooton et al. 2018; Wong et al. 2020; Jones et al. 2022).

Atmospheric modelling of UHJs, further supported by observations, have indicated that the atmospheres of these planets are characterised by a strong temperature inversion, typically occurring between the millibar and microbar level (e.g. Lothringer et al. 2018; Kitzmann et al. 2018; Gandhi & Madhusudhan 2019; Lothringer & Barman 2019; Borsa et al. 2022; Yan et al. 2022; Fossati et al. 2020, 2021, 2023). Several studies have aimed to identify the primary cause of the temperature inversion. Hydrogen Balmer line heating has been proposed as a key agent responsible for the temperature inversion (García Muñoz & Schneider 2019). Several studies also searched for molecules (e.g. TiO, VO, FeH) that are believed to cause temperature inversions but have had little to no success (e.g. Nugroho et al. 2020; Kesseli et al. 2020; Johnson et al. 2023). Different models of UHJ atmospheres predict that the temperature inversion is mostly driven by metal-line absorption of stellar UV radiation, more specifically of the near-UV (NUV) radiation, which is particularly strong for intermediate-mass stars (e.g. Lothringer et al. 2018; Lothringer & Barman 2019; Fossati et al. 2021, 2023). This is supported by numerous detections of neutral and singly ionised species in the atmospheres of UHJs.

Fossati et al. (2021,2023) showed that non-local thermodynamical equilibrium (NLTE) effects enhance the temperature inversion driven by metal-line absorption of stellar UV radiation. As a matter of fact, they have shown that FeII and MgII are, respectively, the species mostly responsible for heating and cooling the planetary atmosphere and that NLTE leads to significant FeII over-population and under-population of MgII, increasing the magnitude of the inversion. This occurs because the NLTE-driven over-population of excited FeII significantly increases the absorption of stellar NUV radiation, further increasing the heating rate. Therefore, accurate characterisation of the atmosphere is only possible if NLTE effects are accounted for. This means that transmission spectra computed while accounting for NLTE effects have been the only reliable forward models capable of fitting the observed hydrogen Balmer line profiles of the UHJs KELT-9b and MASCARA-2b/KELT-20b (see also Stangret et al. 2024; D’Arpa et al. 2024).

WASP-178b is a UHJ that orbits an intermediate-mass star hotter than 8250 K, though there are contrasting results in the literature on the actual value of the stellar effective temperature (Hellier et al. 2019; Rodríguez Martínez et al. 2020; Fouesneau et al. 2022; Lothringer et al. 2025). Lothringer et al. (2022) presented optical and NUV low-resolution (R ∼ 70) observations conducted with the WFC3/UVIS instrument on board the Hubble Space Telescope (HST). They used the NUV observations to retrieve the temperature-pressure (TP) profile while assuming local thermodynamical equilibrium (LTE), obtaining an inversion characterised by a temperature in the upper atmosphere that is about 3000 K higher than the temperature in the lower atmosphere. Furthermore, they interpreted the strong NUV absorption as being indicative of the presence of SiO, which is considered to be the main Si-bearing species in high-temperature environments.

Cont et al. (2024) reported the results obtained from highresolution emission spectroscopy observations carried out with CRIRES⁺ at the Very Large Telescope (VLT). The emission spectrum shows clear indications of CO and H₂O emission, indicative of a temperature inversion. Cont et al. (2024) used retrievals (based on LTE models), including optical eclipse measurements from TESS¹ (Ricker et al. 2015) and CHEOPS² (Benz et al. 2021), to constrain the inversion, obtaining an upper limit on the lower atmospheric (∼1 bar) temperature of about 2000 K and a lower limit on the middle atmospheric (∼10⁻⁶ bar) temperature of about 3200 K. The retrievals also constrain the atmospheric metallicity to be solar, or possibly slightly super-solar, with a solar C/O ratio. The presence of CO and H₂O in the atmosphere of WASP-178b was further confirmed by Lothringer et al. (2025) on the basis of infrared transmission spectroscopy observations conducted with JWST. In particular, they ran retrievals, assuming LTE, on the HST (optical and NUV) + JWST transmission spectra and ultimately obtained a solar C/O ratio as the best fit, in agreement with Cont et al. (2024), though about a third of the retrieval scenarios reported very low C/O ratios, which have been excluded by Cont et al. (2024). They also did not detect SiO in the JWST spectrum, but the covered infrared SiO band is significantly weaker than the NUV one, indicating that there is no discrepancy between the NUV detection by Lothringer et al. (2022) and non-detection in the JWST data.

High-resolution transmission spectroscopy observations collected with ESPRESSO (Pepe et al. 2021) at the VLT led to the detection of Hα, Hβ, and Na I D lines with the tentative detection of one of the Mg I b lines (Damasceno et al. 2024). The transmission spectrum has also been analysed using cross-correlation, leading to the detection of Mg I, Fe I, and Fe II (Damasceno et al. 2024). Furthermore, from the analysis of Balmer lines, Damasceno et al. (2024) found significant broadening of the Hα and Hβ lines with widths of 39.6 ± 2.1 km s⁻¹ and 27.6 ± 4.6 km s⁻¹, respectively, that they interpret as a signature of atmospheric escape.

WASP-178b is an ideal candidate for NLTE atmospheric modeling. It provides an opportunity to further explore the impact of NLTE effects in the atmospheres of UHJs and start comparing results obtained for different planets. In addition, it allows us to compare a synthetic NLTE transmission spectrum with NUV observations. This enables us to test if NLTE effects could be an alternative explanation to SiO in reproducing the NUV observations.

This paper is organised as follows. Section 2 presents the methodology and results of a stellar spectroscopic analysis aimed at independently measuring the stellar effective temperature to constrain the stellar XUV emission. Section 3 describes the planetary atmospheric modelling and transmission spectra calculation, the results of which are presented in Section 4. We discuss our interpretation of the observations in Section 5, focusing on the properties of the planetary upper atmosphere and inference of the mass-loss rate (Section 5.1), on the comparison of the obtained NLTE temperature-pressure profile with those of KELT-9b and MASCARA-2b/KELT-20b (Section 5.2), and on the comparison of our model with the NUV (Section 5.3) and optical (Section 5.4) observations. Finally, in Section 6 we draw the conclusions of this work.

2 Stellar atmospheric analysis

An accurate value for the stellar effective temperature (T_eff) plays a critical role in models of the planetary atmosphere. This is because stars cooler than about 8250 K are believed to have strong XUV emission, similarly to Sun-like stars, while for hotter stars the high-energy emission is several orders of magnitude smaller (e.g. Fossati et al. 2018; Günther et al. 2022). For WASP-178, there are somewhat contrasting T_eff values reported in the literature (e.g. ${8640}_{-240}^{+500}$ K, Rodríguez Martínez et al. 2020; 9047 K, Fouesneau et al. 2022; 9360 ± 150 K, Hellier et al. 2019; ${9200}_{-170}^{+200}$ K, Lothringer et al. 2025), and thus we decided to independently conduct a full stellar atmospheric parameter retrieval and abundance analysis, to also test the classification of the star as a metallic-line (Am) star, as reported by Hellier et al. (2019). To this end, we analysed an archival spectrum of the star collected with the ESPRESSO high-resolution spectrograph on March 31, 2022. The spectrum covers the 3770−7900 Å wavelength range at a resolution of R ≈ 140 000 with a signal-to-noise ratio per pixel of ≈ 620 at λ ≈ 5300 Å.

We computed stellar atmosphere models employing the LLMODELS stellar atmosphere code (Shulyak et al. 2004) assuming local thermodynamical equilibrium (LTE), planeparallel geometry, and the approach of Canuto & Mazzitelli (1991, 1992) to model convection. We used the VALD database (Piskunov et al. 1995; Ryabchikova et al. 2015) as a source of atomic line parameters for opacity calculations. We derived the atmospheric parameters from the ESPRESSO spectrum employing a mixture of equivalent widths, calculated through either direct integration or line profile fitting using BINMAG (Kochukhov 2018), and spectral synthesis, through the SYNTH3 spectral synthesis code (Kochukhov 2007). We converted the equivalent widths in LTE abundances using a modified version (Tsymbal 1996) of the WIDTH9 code (Kurucz 1993). We followed an iterative procedure in which every time any of T_eff, surface gravity (log g), microturbulence velocity (v_mic), or abundances changed during the iteration process, we recalculated a new model with the implementation of the last measured quantities. This ensures that the stellar model structure is consistent with the derived abundances (Kochukhov et al. 2009; Shulyak et al. 2009; Fossati et al. 2009, 2010, 2011b).

We estimated T_eff by imposing excitation equilibrium considering TiII, CrI, CrII, FeI, FeII, and NiI, while we estimated log g by imposing ionisation equilibrium for Mg, Si, Ca, Cr, Mn, Fe, and Ni. Finally, we measured v_mic by imposing equilibrium between line abundances and equivalent widths considering CI, CaI, TiII, CrI, CrII, FeI, FeII, and NiI. The use of multiple elements for the determination of the stellar atmospheric parameters increases the reliability of the results compared to an analysis based on iron only (e.g. Ryabchikova et al. 2009). Because of uncertainties in the spectral normalisation, we could not use the hydrogen Balmer lines to further constrain T_eff and log g, but we used them to independently check the final values. In the end, we obtained T_eff = 9100 ± 100 K, log g = 4.1 ± 0.1, and v_mic = 2.5 ± 0.2 km s⁻¹ (Table 1). The T_eff value is in excellent agreement with that of Fouesneau et al. (2022) and lies within 1σ of those of Rodríguez Martínez et al. (2020), Hellier et al. (2019), and Lothringer et al. (2025). Our derived log g value is lower than those present in the literature, but still in agreement within 1σ, except for those derived by Hellier et al. (2019) and Lothringer et al. (2025) that are higher by about 2σ. The derived T_eff value indicates that the star has a low XUV emission.

We estimated the projected rotational velocity (v sin i) and macroturbulence velocity (v_mac) values in two independent ways, namely (1) through fitting synthetic spectra to the observed spectrum of unblended and weakly blended lines and (2) following the procedure of Murphy et al. (2016). The latter involves the comparison of least-squares deconvolution (LSD; Kochukhov et al. 2010) profiles of synthetic spectra broadened with different v sin i and v_mac values with the LSD profile of the observed spectrum to look for the best fit. Employing both techniques, we obtained compatible results and adopt those derived through LSD fitting: v sin i = 8.5 ± 1.5 km s⁻¹, which is in agreement with Hellier et al. (2019) and Damasceno et al. (2024), and v_mac = 5.0 ± 2.0 km s⁻¹ (Table 1). The tension between the v sin i values given by Hellier et al. (2019) and Damasceno et al. (2024) is likely due to the fact that, having ignored v_mac broadening, the uncertainty reported by Hellier et al. (2019) is underestimated³. The higher v sin i value reported by Rodríguez Martínez et al. (2020, ${12.28}_{-0.82}^{+0.78}$ km s⁻¹) might be the consequence of a too low effective temperature and/or of the fact that they did not consider macroturbulence broadening.

The obtained abundances are listed in Table A.1 and presented in Figure 1. The abundance uncertainties account for the line-to-line abundance scatter (i.e. the standard deviation from the average abundance) and the uncertainty due to the error bars on T_eff, log g, and v_mic (see Section 4.1 of Fossati et al. 2010). We obtained an [Fe/H] abundance of 0.23 ± 0.15 dex, in excellent agreement with that obtained by Hellier et al. (2019), while it is likely that the sub-solar iron abundance obtained by Rodríguez Martínez et al. (2020, −0.06 ± 0.34 dex) is driven by the too low T_eff value. The ESPRESSO spectrum shows also the presence of weak HeI features and their strength matches that of a synthetic spectrum computed considering the finally adopted atmospheric parameters and solar He abundance. We anticipate here that in the case of WASP-178, the obtained stellar metallicity is proper only for the stellar photosphere and not for the entire star, and therefore the derived abundance pattern should not be taken as a reference for stellar evolution or comparisons with the planetary atmospheric abundance pattern (see below).

We independently estimated the stellar parameters and abundances by spectral fitting using the ZEEMAN spectral synthesis code (Landstreet 1988; Wade et al. 2001; Folsom et al. 2012). To estimate the stellar parameters, the spectrum was divided into smaller wavelength regions and then fit iteratively. The scatter in the parameter estimates between different spectral regions is reported as the uncertainty. The stellar parameters estimated from this method are T_eff = 9080 ± 100 K, log g = 4.0 ± 0.1, v sin i = 9.6 ± 0.1 km s⁻¹ and v_mic = 2.5 ± 0.3 km s⁻¹. The stellar abundances are estimated by fitting individual lines of each element, and the line-by-line scatter is reported as the uncertainty. In the case of elements with single lines (e.g. Sr, Al and He), the uncertainty is estimated by varying T_eff by 1σ. Finally, we found stellar parameters and abundances within 1σ of those obtained with the previous method.

The abundance pattern is characterised by an underabundance of C and Sc and by an overabundance of the Fe-peak elements, of Sr, Y, Zr, Ba, and of the rare-Earth elements, which is typical of Am stars (e.g. Fossati et al. 2007, 2008). This conclusion is further supported by the somewhat enhanced microturbulence velocity compared to that of stars of similar temperature, which is also typical of Am stars (Landstreet 1998; Landstreet et al. 2009). Chemical peculiarities in A-type stars develop as a result of diffusion processes that alter only the surface (i.e. photospheric) composition (e.g. Michaud 1970). Therefore, the measured photospheric abundance pattern is not representative of that of the whole star.

Rodríguez Martínez et al. (2020) found that the planetary orbital plane is almost perpendicular to the projected stellar equator, which led them to conclude that the star is likely seen pole-on as a consequence of A-type stars being typically fast rotators. Pagano et al. (2024) derived a possible stellar rotation period of about 3 days that, together with the stellar radius and v sin i value, leads to a stellar rotation rate of about 26 km s⁻¹ and inclination angle of about 17 degrees, supporting the almost pole-on geometry. This inferred stellar rotation rate is, however, significantly smaller than the typical rotational velocity of A-type stars (>100 km s⁻¹; e.g. Abt & Morrell 1995; Royer et al. 2007; Sun & Chiappini 2024) and is in line with the Am chemical peculiarities that develop as a result of diffusion processes in stars rotating slower than ∼90 km s⁻¹ (e.g. Charbonneau & Michaud 1991; Fossati et al. 2008). Given the rather uncertain nature of the stellar rotation period measurement, the strongest constrain on the stellar rotation comes from the Am nature of the star that implies a stellar rotational velocity lower than 90 km s⁻¹, leaving large uncertainties on the stellar inclination angle. This rotational velocity is too low to lead to a significant pole-to-equator temperature difference, in agreement with the non-detection of gravity darkening (e.g. von Zeipel 1924; Espinosa Lara & Rieutord 2011) in the TESS and CHEOPS light curves (Pagano et al. 2024).

Fig. 1 Ion abundance relative to the Sun (Asplund et al. 2009) of the WASP-178 atmosphere. On the y-axis, N corresponds to the number of atoms of a given ionic species indicated on the x-axis and legend, while Ntot is the total number of atoms. Black open squares show the abundance of the neutral elements, red open diamonds are for singly ionised elements, and blue asterisks are for doubly ionised elements.

3 Planetary atmospheric modelling and synthetic transmission spectrum

To simulate the structure of the planetary atmosphere, we followed the same procedure as in Fossati et al. (2021) and Fossati et al. (2023), which we briefly describe below. We computed the TP profile of the planetary atmosphere, employing the HELIOS code (Malik et al. 2017, 2019) at pressures higher than 0.3 mbar, and the CLOUDY NLTE radiative transfer code (version 17.03; Ferland et al. 1998, 2013, 2017), through the CLOUDY for Exoplanets (CfE) interface (Fossati et al. 2021; Young et al. 2024), at lower pressures. We applied this separation because HELIOS does not account for NLTE effects that are relevant in the middle and upper atmosphere, while CLOUDY, which considers NLTE effects, is unreliable at densities greater than 10¹⁵ cm⁻³ (see Ferland et al. 2017, for more details). The TP profiles computed by HELIOS and CLOUDY were then joined together to obtain single TP profiles under two different conditions (i.e. LTE HELIOS + LTE CLOUDY and LTE HELIOS + NLTE CLOUDY), which are used to derive the atmospheric chemical composition and transmission spectra.

For the computation of the TP profile with HELIOS, we considered additional opacities not present in the public version of the code (Fossati et al. 2021) and included cross sections of all molecules available from the DACE⁴ database. Most of the corresponding line lists were provided by the ExoMol⁵ project (Tennyson et al. 2016), with a few compiled from the databases HITRAN⁶ (HBr, HCl, O₂, O₃; Gordon et al. 2022) and HITEMP⁷ (N₂O, NO₂; Rothman et al. 2010; Hargreaves et al. 2019). The atomic line opacity was extended to include Si I-II, Ca I-II, and Ti I-II, which we pre-tabulated using the HELIOS-K⁸ package (Grimm & Heng 2015) and the original line lists produced by R. Kurucz⁹ (Kurucz 2018). Our calculations of the HELIOS models covered the 10²−10⁻⁹ bar pressure range divided into 120 layers and assumed a heat redistribution parameter (f), which accounts for the day-to-night side heat redistribution efficiency, equal to 0.25. In this way, the TP profile in the lower atmosphere is comparable to that obtained by Lothringer et al. (2022), who also considered full heat redistribution.

For the planetary atmospheric modelling and transmission spectra calculation, we considered the system parameters listed in Table 2, which lie within 1σ of those of Damasceno et al. (2024), and solar atmospheric composition (Lodders 2003). Since we obtained a T_eff value comparable to that given by Hellier et al. (2019), we decided to employ their system parameters so that our results can be more directly compared to those present in the literature. Nevertheless, to double-check the system parameters are still consistent with the adoption of our stellar effective temperature, we re-derived the stellar and planetary mass and radius. To this end, we extracted the stellar mass and radius from interpolating across PARSEC v1.2S (Marigo et al. 2017) isochrones and evolutionary tracks (Bonfanti et al. 2015, 2016) using the updated T_eff, G-band magnitude, parallax (Gaia Collaboration 2016, 2023), and metallicity as input. Since the metallicity derived from the spectroscopic analysis does not represent that of the whole star, but it is proper only of the photosphere, we chose to adopt a more moderate supersolar metallicity of +0.1 dex, with an uncertainty of 0.15 dex. We then derived the planetary mass and radius using the stellar reflex velocity, orbital inclination, orbital velocity, and transit depth given by Hellier et al. (2019), further assuming no orbital eccentricity. Finally, we obtained M_s = 1.96 ± 0.08 M_⊙, R_s = 1.59 ± 0.04 R_⊙, M_p = 1.61 ± 0.11 M_J, R_p = 1.72 ± 0.05 R_J, which lie within 1σ of the values given by Hellier et al. (2019, see Table 2) and Lothringer et al. (2025), and lead to a planetary reflex velocity fully consistent with that measured by Cont et al. (2024). For the HELIOS and CLOUDY calculations, we employed the synthetic stellar spectral energy distribution (SED) computed with LLMODELS for the atmospheric parameters reported in Section 2, and thus we did not add any extra high-energy emission to the LLMODELS photospheric fluxes (Fossati et al. 2018; Günther et al. 2022).

The 10¹⁵ cm⁻³ upper density limit for CLOUDY lies at a pressure of about 0.3 mbar, while the continuum lies at a pressure of about 60 mbar. This is the pressure at which the HELIOS model gives an optical depth of 2/3 around 7500 Å, which is the wavelength corresponding to the peak efficiency of the EulerCAM photometer (Lendl et al. 2012) used to measure the planetary radius from the transit light curve (Hellier et al. 2019). For this reason, we run CfE setting the planetary transit radius at the reference pressure (p₀) of 60 mbar, although this choice has no impact on the final TP profile (Fossati et al. 2021, 2023). To save computation time, we run Cloudy considering all elements up to Zn and only hydrogen molecules (i.e. H₂, ${\mathrm{H}}_2^{+},{\mathrm{H}}_3^{+}$, HD) because for planets as hot as WASP-178b, the inclusion of all molecules present in the CLOUDY database does not affect the resulting middle and upper atmospheric temperature (Fossati et al. 2021, 2023). In particular, this means that the CLOUDY calculations do not consider SiO. We mimicked atmospheric heat redistribution in the computation of the CLOUDY TP profiles by scaling the input stellar SED multiplying it by a factor f₁ ≤ 1, looking for the one leading to the TP profile best matching the one given by HELIOS around the 10⁻⁴ bar level. We finally adopted f₁ = 1.0, but we note that the value employed for f₁ has no significant (>200 K) impact on the TP profile at pressures lower than about 10⁻⁵ bar (Fossati et al. 2021, 2023). This is likely because the heating rate is relatively high for all values of f₁ exceeding 0.25, but also substantially offset by thermostatic radiative cooling rates.

Finally, we used atmospheric structure (i.e. TP and abundance profiles) both in LTE and NLTE to compute the LTE and NLTE transmission spectra using CLOUDY ranging from the far-ultraviolet to the near-infrared at a spectral resolution of 100 000. To this end, we followed the procedure described by Young et al. (2020) and Fossati et al. (2020), in which the atmosphere is divided into 100 layers equally spaced in log p. We computed the NLTE transmission spectrum employing the NLTE TP profile and enabling NLTE in the CLOUDY calculations. Similarly, we derived the LTE transmission spectrum using the LTE TP profile and imposing the LTE assumption for the CLOUDY transmission spectrum calculation. Within CLOUDY, the difference between LTE and NLTE computations lies in the fact that in LTE level populations are computed using the Boltzmann equation with the exception of the first two levels of hydrogen, which are always computed in NLTE. Furthermore, by default CLOUDY include photoionisation, also in LTE. Instead, HELIOS does not account for photoionisation, and thus thermal ionisation and level populations were computed using the Saha and Boltzmann equations, respectively.

The TP profile is for the sub-stellar point, but the computation of transmission spectra implies a different geometry (i.e. from emission geometry to transmission geometry), which then calls for a new calibration of the reference pressure, p₀, that corresponds to the reference radius, R₀, that is, the measured transit radius, R_p. We performed this calibration following Fossati et al. (2023), which gives more accurate results compared to that followed by Fossati et al. (2021). In short, we computed transmission spectra at different p₀ values in the 0.001−0.1 bar range, then convolved each synthetic transmission spectrum with the EulerCAM efficiency curve to obtain R_p/R_s as a function of p₀. The final transmission spectrum is that computed with the p₀ value leading to best match the observed R_p/R_S of 0.111 (Hellier et al. 2019). For both LTE and NLTE theoretical transmission spectra, we find a best-fitting reference pressure value of 0.002 bar.

4 Results

4.1 Atmospheric temperature-pressure structure

We joined the helios and Cloudy TP profiles at a pressure of about 4 × 10⁻⁴ bar resampling the final theoretical TP profile over 200 layers equally spaced in log p ranging from 10 bar to 5 × 10⁻¹² bar. This range is wide enough to cover the formation region of most UV, optical, and infrared lines considered for the computation of the theoretical transmission spectrum. Figure 2 shows the composite theoretical TP profile in comparison to the original helios and Cloudy theoretical profiles. The TP profile is very similar to that obtained for MASCARA-2b/KELT-20b using the same modelling scheme (Fossati et al. 2023): isothermal at pressures higher than about 10 mbar and lower than about 10⁻⁸ bar, with a linear temperature rise from about 2200 K to about 8100 K in between.

The comparison of the TP profiles obtained in LTE and NLTE (Figure 2) indicates that the LTE assumption appears to be appropriate at pressures higher than about 10⁻⁴ bar, but it underestimates the temperature by about 3000 K at lower pressures. This difference in upper atmospheric temperature between LTE and NLTE TP profiles is similar to that obtained for MASCARA-2b/KELT-20b (T_eq ≈ 2260 K; Talens et al. 2018), which has an equilibrium temperature comparable to that of WASP-178b (Fossati et al. 2023). Interestingly, this temperature difference in the upper atmosphere is about 1000 K larger than that obtained for KELT-9b (T_eq ≈ 4050 K; Gaudi et al. 2017), which has instead a significantly higher equilibrium temperature (see Section 5.2, for a more thorough comparison of the TP profiles of the three planets).

Fig. 2 Theoretical atmospheric structure obtained for WASP-178b. Top: HELIOS (dashed orange line) and CLOUDY (NLTE; dashed cyan line) TP profiles. The solid black line shows the composite TP profile. The horizontal dotted black line indicates the location of the continuum according to the HELIOS model. The horizontal dash-dotted red line gives the location of CLOUDY’s upper-density limit of 1015 cm−3. The dashed dark green line shows the CLOUDY TP profile computed assuming LTE. Bottom: Pressure (black; left y-axis) and temperature (red; right y-axis) composite theoretical profiles as a function of the planetary polar radius. For context, the L 1 point lies at about 4.00 Rp, while the Roche radius in the limb direction lies at about 2.67 Rp.

4.1.1 Heating and cooling

Following Fossati et al. (2021) and Fossati et al. (2023), we identify the elements primarily responsible for the difference between the LTE and NLTE TP profiles by computing CLOUDY NLTE TP profiles excluding one of the elements at a time. Similarly to previous results, we find that Fe dominates heating and Mg dominates cooling in the middle and upper atmosphere (Figure 3). Removing Fe leads to a 1000−2500 K cooler temperature in the middle and upper atmosphere, while excluding Mg leads to an about 300 K hotter TP profile, particularly in the upper atmosphere, at pressures lower than the μ bar. The other elements, instead, contribute less than 50 K to the atmospheric heating or cooling, with the exception of Ca for which we find that removing it from the calculation of the TP profile leads to an about 200 K hotter temperature around the 10 μ bar level.

Figure 3 suggests that metals play a pivotal role in the heating and cooling of the middle and upper atmosphere. To gather more insights, we extracted from the CLOUDY run the three main heating and cooling agents that we show in Figure 4. We find that metal line absorption (particularly of Fe; Figure 3) is the main heating mechanism, with Hi photoionisation playing a minor role, while numerous species contribute to the cooling of the middle and upper atmosphere, with Mg being the dominant coolant in our steady state model. We further explore the general impact of photoionisation on the atmospheric heating by running CLOUDY without photoionisation, which leads to an upper atmospheric temperature that is just a few hundred Kelvin cooler than that obtained with photoionisation, while the temperature remains the same in the middle atmosphere (see Figure 3). This indicates that photoionisation heating does not play a major role in the atmospheric thermal balance (see also the magenta line in the top panel of Figure 4).

To better understand the relative roles of Mg and Fe in modifying the TP profile, we took the output obtained after the last NLTE CfE iteration and used it as input for a further CLOUDY run, but we considered only FeI (i.e. FeI is not allowed to ionise) and FeII (i.e. FeII is not allowed to ionise or recombine), fixing the FeI or FeII density profile to that obtained when accounting for all elements and ions. We then applied the same procedure for Mg. Figure 5 shows the temperature profiles obtained in all these cases. As previously found for KELT-9b and MASCARA-2b/KELT-20b, FeII dominates the heating in the middle and upper atmosphere and FeI contributes to the cooling in the middle atmosphere, as shown by the fact that without FeI we obtain a higher temperature by a few hundred Kelvin at pressures higher than about 10⁻⁸ bar. As also shown above, removing Mg increases the atmospheric temperature across the entire middle and upper atmosphere. Considering just Mgi being present in the planetary atmosphere leads to a TP profile that is similar to that obtained without Mg, which suggests that Mg ions are responsible for the cooling at altitudes corresponding to pressures lower than 0.1 mbar. Considering only Mgir leads to a slightly cooler TP profile within the 10⁻⁴ and 10⁻⁸ bar pressure range, which is where Mg cooling is the strongest (see Figure 4), while at higher altitudes the TP profile is almost identical to that obtained considering all elements and ions. This indicates that both MgI and MgII contribute to the cooling in the 10⁻⁴−10⁻⁸ bar pressure range, with MgI being more important that MgII, while MgII is the dominant Mg cooling species at lower pressures.

Removing Fe from the atmosphere decreases the temperature significantly more than removing photoionisation, bringing it close to that obtained in LTE. Therefore, the temperature difference between the LTE and NLTE TP profiles is driven by the differences in atomic level populations of Fe, particularly of FeII, driven by NLTE effects, which increase the level population of long-lived states presenting transitions mainly in the NUV band, where the stellar emission is strong. Instead, the fact that by removing Mg cooling the temperature increases by only 300 K indicates that other cooling mechanisms replace it and are nearly as efficient as Mg cooling in the new steady state. The bottom panel of Figure 4 shows that this other cooling agent is hydrogen. Therefore, in general atomic line cooling is important and it is hard to heat the atmosphere above 9000 K without removing all major cooling agents.

Fig. 3 Comparison between the NLTE CLOUDY TP profiles obtained considering all elements up to Zn (solid black line) and excluding Fe (solid red line), Mg (solid blue line), Ca (solid green line that is almost overlapping with the black line), or photoionisation (yellow). The dashed orange line shows the HELIOS (i.e. LTE) TP profile for reference.

Fig. 4 Heating and cooling contributions in the middle and upper atmosphere (i.e. at pressures lower than CLOUDY’s density limit) of WASP-178b. Top: Contribution to the total heating as a function of pressure. At each pressure bin, the plot shows the three most important heating processes. The main heating processes are hydrogen photoionisation (red; photoionisation of HI lying in the ground state), metal line absorption (blue), H− absorption (green), photoionisation of hydrogenic species (magenta), collisions with H2 (brown), and Compton heating (i.e. electron absorption; yellow). Bottom: Contribution to the total cooling as a function of pressure. At each pressure bin, the plot shows the three most important cooling agents. The main cooling agents are Mg (red), Ca (blue), C (green), ${\mathrm{H}}_2^{+}$ (magenta), H (brown), S (orange), HD (i.e. hydrogen deuteride; yellow), Na (grey), and Si (light blue). Metal line absorption (primarily from FeII) and Mg line emission are the primary heating and cooling mechanisms, respectively.

Fig. 5 Temperature-pressure profiles obtained from isolating the impacts of Fe (left) and Mg (right) ions. The black line is for when all elements are accounted for (same as in Figures 2 and 3), while the other lines were obtained considering that in the planetary atmosphere, Fe or Mg exist only in the form of a neutral species (blue) or a singly ionised species (green) or is absent (red).

Fig. 6 Density (relative to the total density of hydrogen)-pressure profiles for neutral hydrogen (HI; solid black line), protons (HII; red), molecular hydrogen (H2.; dark green), ${\mathrm{H}}_2^{+}$ (blue), ${\mathrm{H}}_3^{+}$ (magenta), H− (orange), and electrons (e−; dashed black).

4.2 Atmospheric abundance profiles

To obtain a homogeneous synthetic physical and chemical structure throughout the entire atmosphere, we passed the composite TP profile to CLOUDY (see Fossati et al. 2021, for more details). We used this final CLOUDY model to extract the abundance profiles of the most important species. Figure 6 shows the CLOUDY density profiles with respect to the total hydrogen density of the main hydrogen-bearing species, that is neutral hydrogen (HI), protons (HII), H⁻, molecular hydrogen (H₂), ${\mathrm{H}}_2^{+}$ and ${\mathrm{H}}_3^{+}$, plus electrons (e⁻). The abundance profiles are very similar to those previously obtained for MASCARA-2b/KELT-20b. Neutral hydrogen is the most abundant species in the middle and upper atmosphere (<5 × 10⁻³ bar), with ionised hydrogen taking over just at the very top of the atmosphere, at pressures lower than 10⁻¹¹ bar. Instead, molecular hydrogen dominates in the lower atmosphere, but its abundance decreases rapidly at pressure lower than 5 × 10⁻³ bar as a result of thermal dissociation. As previously found for other UHJs (e.g. Arcangeli et al. 2018; Fossati et al. 2021, 2023), H⁻ is relatively abundant with its density first increasing with decreasing pressure up to the 1 μbar level and then decreasing at lower pressures.

Figure 7 shows the mixing ratio as a function of pressure given by CLOUDY for some of the most abundant elements. Oxygen remains neutral up to the top of the atmosphere, where ionisation becomes dominant around 10⁻¹¹ bar, while CII becomes the dominant carbon species above the 10⁻⁷ bar level. In the lower atmosphere (>10 mbar), Ca, Na, and K are the main contributors to free electrons, while around the mbar level Fe, Mg, and Si ionisation contribute significantly to the electron population, which then starts following closely the proton abundance at lower pressures as a result of the increasing importance of hydrogen thermal ionisation. Among the species shown in Figure 7, Ca is the only element for which the second ionised species become dominant within the simulated pressure range.

Fig. 7 Atmospheric mixing ratios obtained for metals. Top: mixing ratios for hydrogen (solid black line), H2 (dotted black), He (red), C (blue), O (dark green), Na (magenta), K (orange), and electrons (bright green). Bottom: same as the top panel but for Mg (red), Si (blue), Ca (dark green), and Fe (magenta). The hydrogen, H2, and e− mixing ratios are shown in both panels for reference. Neutral (XI), singly ionised (XII), and doubly ionised (XIII; present just in the bottom panel) species are shown as solid, dashed, and dash-dotted lines, respectively.

4.3 Synthetic transmission spectrum

To enable comparisons with observations, we computed the synthetic NLTE and LTE transmission spectra ranging from the far-UV to the near-infrared (i.e. from 912 Å to 2.85 μm). To this end, we followed the procedure described by Young et al. (2020) and Fossati et al. (2020), which consists in first mapping the TP profile onto concentric circles and calculating the lengths through successive layers of atmosphere along line-of-sight transmission chords. These lengths, along with the atmospheric properties of their respective layers, are then stacked and entered into CLOUDY as the line of sight transmission medium. The final transmission spectrum of the planet is then computed by adding up the single layer spectra, weighted by their relative area projected on the stellar disc and accounting for the planetary impact parameter.

Figure 8 shows the LTE and NLTE synthetic transmission spectra covering the entire wavelength range, while similar plots of specific wavelength ranges can be found in Appendix (Figures A.2, A.3, A.4, A.5, A.6, and A.7). For completeness and easier comparison with observations, we show in Figure A.8 the transit depth difference between the theoretical LTE and NLTE transmission spectra.

Figure 8 indicates that, in general, the LTE transmission spectrum tends to underestimate absorption, mostly due to the fact that the LTE model is generally cooler than the NLTE model, though there are some typically weak features showing the opposite effect. The obtained behaviour is similar to that previously found for MASCARA-2b/KELT-20b and KELT-9b (Fossati et al. 2021, 2023), but the amplitude of the difference between the LTE and NLTE transmission spectra is significantly larger in the case of WASP-178b. For example, in the case of the Hα line the impact of NLTE effects on the transmission spectrum is of about 15% for KELT-9b and about 11% for MASCARA-2b/KELT-20b, while we find it to be about 32% for WASP-178b, that is more than twice that found for KELT-9b. This is likely due to the more inflated nature of the atmosphere of WASP-178b resulting from the differences in planetary mass. Figure 8 also shows that the core of the strongest UV absorption lines predicted by the NLTE model lies beyond the Roche lobe. Given that our simulations assume one-dimensional geometry, the actual strength of these lines should be taken with caution as significant differences can be expected when compared with observations.

As in previous analyses, we find the strongest deviation from LTE in the UV band, with deviations around the 10−40%, except for strong resonance lines (e.g. Lyα, MgII H&K) where the deviation reaches 50%. In the optical and near-infrared bands, we find deviations from LTE typically below 10%, except for specific lines, such as the hydrogen lines (Balmer and Paschen) for which we find deviations larger than 20%. Similarly to the case of MASCARA-2b/KELT-20b, the absorption of the OI triplet at about 7780 Å, which shows a prominent deviation from LTE, appears to be weaker than that of the Hα line, while for KELT-9b the two lines had comparable strengths (Fossati et al. 2021, 2023; Borsa et al. 2021). From an observational point of view, this result suggests that it will require data of significantly higher quality than that of the data analysed by Borsa et al. (2021) to detect the OI triplet in the atmosphere of both MASCARA-2b/KELT-20b and WASP-178b.

Also in the case of WASP-178b, the difference between the LTE and NLTE transmission spectra, particularly in the UV band, is partly due to the difference in underlying TP profiles, in addition to the different physical assumptions in the radiative transfer calculations. Most of the spectral lines lying in the UV belong to ionised species, which are more abundant in the NLTE model as a consequence of the higher temperature of the NLTE TP profile compared to the LTE TP profile, particularly in the line forming region. Furthermore, the hotter temperature of the NLTE model leads to a higher pressure scale height compared to the LTE model, which also affects the strength of the lines in the transmission spectrum.

Computationally, it is not yet possible to run full retrievals that account for NLTE effects; however, it might be interesting to explore the possibility of deriving abundances and/or abundance profiles on the basis of a fixed TP profile computed using forward modelling that accounts for NLTE effects. To this end, we used the NLTE TP profile as input to ClOUDY to compute a transmission spectrum assuming LTE. Figure 9 shows a comparison between the NUV and optical transmission spectrum computed with partial NLTE (in red) and the one computed with full NLTE (in black). In the optical, the difference is typically small and below 10%; however, we remark that CLOUDY in LTE still assumes NLTE level populations for the first two hydrogen states, which is most likely why there is no difference between the two Hα lines. In the NUV, the difference increases significantly up to around 20%. Nevertheless, a comparison between Figure 8 and 9 shows that the partial NLTE transmission spectrum is on average closer than the LTE one to the full NLTE transmission spectrum. This indicates that retrievals based on the partial NLTE approach might lead to more reliable results than retrievals fully based on the LTE assumption.

Fig. 8 Comparison between the theoretical LTE (red) and NLTE (black) transmission spectra computed at a resolution of 100 000 in the UV and optical range (top) and in the near-infrared (bottom). Within each plot, the bottom panel shows the relative difference between LTE and NLTE (in percent). In the top panel, the blue horizontal dashed line indicates the location of the Roche radius in the planetary limb direction (i.e. perpendicular to the star-planet connecting line).

Fig. 9 Comparison between theoretical transmission spectra computed with full NLTE (A; black; same as the black line in Figure 8), namely the NLTE TP profile and accounting for NLTE effects for the computation of the transmission spectrum (TS), and with partial NLTE (B, red), namely NLTE TP profile and assuming LTE for the computation of the transmission spectrum, at a resolution of 100 000 in the NUV and optical range. The bottom panel shows the relative difference between the two spectra in percent. In the top panel, the blue horizontal dashed line indicates the location of the Roche radius in the planetary limb direction.

Fig. 10 Line formation analysis in transmission geometry. Top: relative peak absorption for Hα (black), the MgIb line at ∼5183 Å (red), the Cal line at ∼6439 Å (blue), the Cair K line (green), the FeI line at ∼4071 Å (magenta), and the FeII line at ∼5169 Å (orange). Bottom: relative absorption in a 100 Å spectral bin centered at the MgI resonance line at ∼2850 Å (black), the MgII h&k resonance lines at ∼2800 Å (red), the FeI NUV absorption band at ∼2720 Å (blue), the FeII resonance line at ∼2600 Å (green), and the FeII NUV absorption band at ∼2370 Å (magenta) as a function of pressure. In both panels, the dashed line shows the TP profile (top x-axis). The line formation region corresponds to the pressure values in which the relative absorption ranges between zero and one.

4.4 Line formation region

We used the Cloudy outputs obtained for each atmospheric layer in the synthetic transmission spectra calculations to estimate the line formation region in transmission geometry for a few of the strongest lines detected in the optical transmission spectra of UHJs: Hα, the reddest of the MgI b lines at ∼5183 Å, the CaI line at ∼6439 Å, the Cair K line at ∼3933 Å, the Fer line at ∼4071 Å, and the FeII line at ∼5169 Å. In particular, we started from the top of the atmosphere (at 5 × 10⁻¹² bar) looking at the fraction of light absorbed through increasingly thick atmospheric layers at the centre of each line. In this scheme, the top boundary of the formation region of a given line (i.e. at low pressure) is set where, with increasing pressure of the lower considered layer, the line strength in the transmitted spectrum starts deviating from zero, while the bottom boundary of the line formation region (i.e. at high pressure) is set at the pressure for which the line strength stops increasing with increasing pressure.

The top panel of Figure 10 shows the relative peak absorption of the considered lines as a function of pressure, and thus the line formation region corresponds to the pressure range in which the relative line peak absorption lies between zero and one. Among those considered, Hα is the line forming higher up in the atmosphere, in the 1−10 nbar range, while the CaI line at ∼6439 Å and the FeI line at ∼4071 Å are those forming at the higher pressures. In general, the lines considered here have a rather narrow formation region, spanning 1–2 orders of magnitude in pressure, except for the FeII line at ∼5169 Å that has a line formation region spanning about five orders of magnitude in pressure (10⁻⁵−10⁻¹⁰ bar). The line formation region shown in Figure 10 reflects the abundance profile of each ion.

Figure 2 and the top panel of Figure 10 indicate that in the optical band lines form at a temperature higher than 5000 K and that the Hα line forms in the hottest part of the atmosphere, where the temperature is higher than 8000 K. Similarly to the cases of KELT-9b and MASCARA-2b/KELT-20b, the complete lack of Lyα emission from the host star implies that the excitation of HI to the n = 2 level, which leads to the formation of the Balmer lines, has to occur mostly thermally. Furthermore, because of the very low EUV emission of the host star, photoionisation and subsequent recombination rates to the n = 2 level are low compared to thermal excitation (Fossati et al. 2023). Figure 10 also shows the pressure range probed by the lines considered here and thus the atmospheric layers constrained by observations through both forward modelling and retrieval analyses.

We carried out the same line formation analysis described above, but computing Cloudy models at a much lower spectral resolution of 100. In particular, we computed the average relative absorption in a 100 Å range centered on the strongest spectral features in the NUV band, namely the MgI resonance line at ∼2850 Å, the MgII h&k resonance lines at ∼2800 Å, the Fer NUV absorption band at ∼2720 Å, the FeII resonance line at ∼2600 Å, and the Fell NUV absorption band at ∼2370 Å. These are the approximate spectral resolution and wavelength bin width covered by the WFC3/UVIS observations presented by Lothringer et al. (2022). The bottom panel of Figure 10 shows the average relative absorption as a function of pressure, indicating that at the resolution of WFC3/UVIS the observations probe mainly the ≈ 10⁻⁵−10⁻¹ bar pressure range.

Figure A.1 shows a comparison of the LTE and NLTE TP profiles obtained joining helios and Cloudy with those presented by Lothringer et al. (2022), Cont et al. (2024), and Lothringer et al. (2025) assuming LTE. All results support the presence of a temperature inversion, as also evinced by the observations, but the magnitude of the inversion is significantly larger when accounting for NLTE effects. In the lower atmosphere (>1 mbar), all TP profiles are comparable (particularly when considering that at high pressures the temperature of Cont et al. 2024 is a lower limit), which is partly driven by the fact that the HELIOS calculations have been set up to produce a profile comparable to that of PHOENIX. The main difference lies in the middle and upper atmosphere, where the profiles of Lothringer et al. (2022) and Cont et al. (2024), as well as our LTE profile, predict a temperature that is a few thousand kelvins cooler than our NLTE profile. In the 10⁻⁴−10⁻⁵ bar range, the temperature retrieved by Lothringer et al. (2025) is particularly low compared to all other profiles, although they obtained it including in the retrieval the observed NUV transmission spectrum.

In the case of the TP profiles obtained through forward modelling, the difference between our result and previous results can be ascribed mainly to the impact of accounting for NLTE effects. Instead, in the case of the retrieved profiles the difference with our NLTE profile is due not only to the different assumptions in the radiative transfer modelling (i.e. LTE vs NLTE), but also to the fact that the observations probe mainly the lower part of the atmosphere, at pressures higher than 10⁻⁵ bar.

Fig. 11 Atmospheric sound speed (black; left y-axis) and Jeans escape parameter (red; right y-axis) profiles as a function of pressure computed in the direction of the sub-stellar point (i.e. the Roche lobe lies at the L1 point).

5 Discussion

5.1 Upper atmosphere and mass-loss rate

CLOUDY is a hydrostatic code, but part of the modelled atmosphere lies beyond the Roche lobe (Figure 8), which calls for checking the validity of the hydrostatic approximation throughout the considered pressure range. Figure 11 shows the profiles of the Jeans escape parameter and of the sound speed as a function of pressure. For the computation of the Jeans escape parameter, we used Equation (7) of Fossati et al. (2017, see also Volkov et al. 2011), setting the bulk flow velocity equal to zero. In this way, the Jeans escape parameter accounts for the gravitational potential difference between every point in the atmosphere and the Roche lobe, where the Jeans escape parameter goes to zero. The strong metal lines in the optical band form roughly in the 10⁻⁴−10⁻⁸ bar range, where the Jeans escape parameter is larger than five, and thus the atmosphere is hydrostatic. This is also the case of the atmospheric regions probed by the WFC3/UVIS and JWST observations. However, the Hα line forms in the 1–10 nbar range, where the Jeans escape parameter is between two and five, and thus the bulk outwards velocity starts becoming non-negligible, though the velocity remains subsonic. Therefore, in the case of WASP-178b, in contrast to KELT-9b and MASCARA-2b/KELT-20b, high-resolution transmission spectroscopy observations of the H I Balmer lines probe the upper atmosphere and can thus provide information on the ongoing hydrodynamic escape and eventually constrain the mass-loss rate. Our NLTE model supports the suggestion of Damasceno et al. (2024) that the significant broadening they found characterising the hydrogen Balmer line profiles can be related to atmospheric escape. However, the velocities Damasceno et al. (2024) inferred from the line broadening (39.6 ± 2.1 km s⁻¹. for Hα and 27.6 ± 4.6 km s⁻¹ for Hβ) are significantly higher than the sound speed in the line formation region, indicating that only part of the broadening might be related to the motion of the escaping material, while other processes (e.g. thermal and non-thermal broadening and winds) must contribute to it.

We estimated the mass-loss rate from WASP-178b as $$\dot{M}\approx 4\pi x_R{z}_R\frac{{\rho }_R}{2\sqrt{\pi }}\sqrt{\frac{2k{T}_R}{m_R}},$$(1)

where x_R is the radial distance to the L1 point, z_R is the radial distance to pole, ρ_R is the mass density at the Roche lobe boundary, T_R is the temperature at the Roche lobe boundary, and m_R is the mean molecular weight at the Roche lobe boundary. Detailed simulations of atmospheric escape show that this approach is reasonably accurate, to a factor of a few, as long as the exobase is below the Roche lobe boundary (Koskinen et al. 2022). Based on our model of the atmosphere, we obtain a mass-loss rate of 4.5 × 10¹⁰ g s⁻¹, which is comparable to that obtained for the archetypal hot Jupiter HD 209458b. This mass-loss rate is likely high enough to affect the structure of the atmosphere above the heating peak, by causing temperatures to decrease with altitude above the ∼1 nbar level. As we indicate above, this means that self-consistent models of escape and thermal structure are required to simulate absorption signatures of strong atomic lines probing lower pressures at high radial distances from the planet.

Fig. 12 Comparison between the theoretical NLTE TP profiles of WASP-178b (black; this work), MASCARA-2b/KELT-20b (red; from Fossati et al. 2023), and KELT-9b (blue; from Fossati et al. 2021). For context, the legend also gives the effective temperature of the host stars from Section 2 (for WASP-178b), Talens et al. (2018, for MASCARA-2b/KELT-20b), and Borsa et al. (2019, for KELT-9b).

5.2 WASP-178b vs MASCARA-2b/KELT-20b vs KELT-9b

The HELIOS and CLOUDY models have been used to construct the theoretical TP profiles, accounting for NLTE effects in the middle and upper atmosphere, of three UHJs: WASP-178b, MASCARA-2b/KELT-20b, and KELT-9b. Figure 12 shows a direct comparison between the theoretical TP profiles of the three planets¹⁰.

As a result of the hotter host star, shorter orbital separation, and therefore a higher planetary equilibrium temperature, the lower atmosphere of KELT-9b is significantly hotter than that of the other two planets, which have similar equilibrium temperatures of about 2400 K (Lund et al. 2017; Talens et al. 2018; Hellier et al. 2019). All three planets have a significant temperature inversion in the middle atmosphere and the three TP profiles get closer to each other with decreasing pressure. The temperature inversion is more significant for the two cooler planets, and in particular for WASP-178b, which shows the strongest inversion with a temperature difference of about 6000 K between the top and bottom of the atmosphere. Remarkably, at the top of the atmosphere, the temperature of WASP-178b is comparable to that of KELT-9b, while that of MASCARA-2b/KELT-20b is about 500 K lower.

The characteristics of the TP profiles in the upper atmosphere do not appear to be directly related to either the temperature of the host star or the orbital separation. Furthermore, the host stars are very similar in terms of high-energy emission as all have low XUV emission. One of the main parameters distinguishing the three planets is the bulk density, which impacts the atmospheric pressure scale height and thus the penetration depth of the stellar radiation in the planetary atmosphere and in turn the shape of the TP profile. As a matter of fact, among the three planets WASP-178b has the lowest bulk density (≈0.35 g cm⁻³; Hellier et al. 2019), while MASCARA-2b/KELT-20b has the highest bulk density (≈0.71 g cm⁻³; Talens et al. 2018), with KELT-9b lying roughly in between (≈0.49 g cm⁻³.; Borsa et al. 2019). Furthermore, WASP-178b and KELT-9b have Roche radii that are significantly closer to the planetary radius compared to MASCARA-2b/KELT-20b. At the same time, atomic line cooling could act as a thermostat such that independently of the energy input, the heating rate is offset by a cooling rate that increases with temperature. Therefore, a parameter space study could lead to the explanation of the origin of the behaviour shown in Figure 12.

5.3 Comparison with HST observations

Lothringer et al. (2022) presented NUV transmission spectroscopy observations performed at low-resolution (R ≈ 70) with the WFC3/UVIS instrument on board HST, which they interpreted using both forward and retrieval modelling based on the PHOENIX code assuming LTE. They fitted the data by adding SiO opacity to the model comprising Mg and Fe, among other species, though the data could have also been fitted by increasing the atmospheric metallicity. In the end, they preferred the first option. This choice was based on a NUV medium-resolution (R ≈ 30 000) transmission spectrum collected with the STIS spectrograph on board HST, which, despite the higher spectral resolution, did not show extra-absorption at the position of strong Fe and Mg resonance lines compared to the UVIS spectrum.

Below, we present a comparison of the LTE and NLTE synthetic transmission spectra with the HST observations. We extracted the low-resolution WFC3/UVIS transmission spectrum from Lothringer et al. (2022), but we could not do the same for the STIS mid-resolution spectrum because Lothringer et al. (2022) did not provide it in machine readable format and their extended data Figure 8 shows just short portions (≈90 Å wide) around the MgII h&k lines and the strongest FeII NUV resonance lines. Therefore, to enable a comparison of the synthetic transmission spectra with the STIS observations, we downloaded the data from the HST archive and analysed them independently. The details of the STIS data analysis are given in Appendix B.

5.3.1 WFC3/UVIS

As shown in Figure A.1, by considering NLTE effects, we calculate a significantly higher temperature in the middle and upper atmosphere compared to those obtained by Lothringer et al. (2022) and Lothringer et al. (2025) using forward modelling or retrievals. This leads to larger pressure scale heights, and thus a larger extent of the atmosphere, which, together with the stronger absorption features, might allow us to reproduce the UVIS data without the need of SiO absorption or super-solar abundances.

The top panel of Figure 13 shows the comparison between the theoretical transmission spectra we obtained in LTE and NLTE with the UVIS data. The theoretical transmission spectra have been obtained by averaging the high-resolution spectra shown in Figure 8 for each of the bins used to obtain the UVIS transmission spectrum. We find that in the NUV, the NLTE transmission spectrum provides an excellent fit to the data with a reduced ${\chi }_v^2$ of 0.95 in the 2000–3000 Å range, with 11 degrees of freedom, while the LTE transmission spectrum lies significantly below the observations returning a higher reduced ${\chi }_v^2$ of 3.77. Indeed, accounting for NLTE effects increases NUV absorption compared to what we obtained assuming LTE and this enables us to fit the UVIS observation without adding SiO absorption or increasing metallicity.

However, in the optical band, the UVIS transmission spectrum lies slightly below the models, indicative of the presence of a small offset between the UVIS data and the broad-band radius ratio. Given the magnitude of the residual structured red noise present in the data, this offset might not be significant. However, to be conservative, we shifted the LTE and NLTE models to match the UVIS data at wavelengths longer than 4000 Å (bottom panel of Figure 13). This leads to a worse fit to the data in the NUV band as indicated by the reduced ${\chi }_v^2$ values of 2.59 and 6.69 obtained for the NLTE and LTE models, respectively. Therefore, taking into account the possible offset, we find that to match the UVIS data one might require additional absorption across the entire NUV band. This could be achieved indeed by including some SiO absorption, or by a slight increase in the metallicity, which would agree with the results of Cont et al. (2024) and Damasceno et al. (2024). In fact, for example, a slight increase in Fe abundance would not only lead to an increase in the Fe opacity, and thus NUV absorption, but also an increase in the atmospheric temperature in the middle and upper atmosphere, which would further favour Fe thermal ionisation and increase the pressure scale height, with the net effect of increased NUV absorption. Also, we cannot rule out the possibility that the middle and upper atmosphere is even hotter than predicted by the NLTE model, which would further strengthen NUV absorption. In fact, low-resolution NUV observations of the UHJ WASP-189b conducted with the Colorado Ultraviolet Transit Experiment (CUTE) CubeSat suggest that the temperature in the upper atmosphere may be significantly higher than that predicted by state-of-the-art hydrodynamic models (Sreejith et al. 2023). Finally, it is important to note that our model is one-dimensional while three-dimensional effects are likely to be important particularly at the terminator region. This might impact the TP profile, and in turn the shape of the transmission spectrum.

Fig. 13 Comparison of the NLTE (black) and LTE (blue) theoretical transmission spectra with the WFC3/UVIS observations published by Lothringer et al. (2022, red). The horizontal error bars shown for the red points correspond to the size of the wavelength bins. The top plot is without scaling of the models, while in the bottom plot the NLTE and LTE models have been scaled to match the average of the UVIS observations at wavelengths longer than 4050 Å, marked by the vertical grey dashed line. In the top plot, the grey horizontal dashed line indicates the broad-band Rp/Rs obtained from the EulerCAM optical transits (Hellier et al. 2019). In each plot, the bottom panel shows the residuals between the observation and each of the models, with a green dashed line at zero for reference.

5.3.2 STIS

The top panel of Figure 14 shows the comparison of the LTE and NLTE synthetic transmission spectra with the mid-resolution (≈9 Å binning) STIS transmission spectrum obtained using the ‘multiple jitter parameter’ analysis method to remove systematic noise, as detailed in Appendix B. The data present a significant scatter around the expected transit depth and the shape of the transmission spectrum does not follow what is predicted by any of the two models. Therefore, to check the quality of the STIS data and systematic removal we compare the STIS transmission spectrum with the UVIS one at the same resolution. To this end, we re-extracted the STIS transmission spectrum from the data, but using the same binning as that of the UVIS data (≈74 Å binning). The bottom panel of Figure 14 shows the comparison of the LTE and NLTE synthetic transmission spectra with the UVIS and STIS transmission spectra, all at the same binning. The lowresolution STIS transmission spectrum does not match the UVIS transmission spectrum. Notably, the STIS transmission spectrum we extracted from the data following the procedure described in Appendix B matches well what presented by Lothringer et al. (2022, see Figure B.4). Furthermore, Appendix B shows that the results obtained from the analysis of the STIS data is model-dependent. Therefore, we conclude that the STIS dataset of WASP-178b may not be reliably used to gather information about the atmosphere of WASP-178b. The residual noise present in the STIS data might result from slit losses caused by the use of a too narrow slit. Therefore, we recommend using a wider slit for future STIS transmission spectroscopy observations¹¹.

Lothringer et al. (2022) fitted the UVIS transmission spectrum using an LTE model and required additional opacity to be able to fit the deeper NUV transit depths shown by the data. They identified in SiO the source of this opacity and used the lack of strong absorption at the position of the MgII h&k and of FeII resonance lines in the STIS data to exclude the scenario in which hydrodynamic escape is responsible for the increased NUV absorption. However, our results indicate that SiO absorption is not required to fit the data. Nevertheless, even if the shape of the UVIS transmission spectrum in the NUV is the result of NLTE effects and not of SiO absorption, our findings support the conclusion of Lothringer et al. (2022) that Fe and Mg have not rained out.

Fig. 14 Comparison between observed and synthetic transmission spectra in the region covered by the STIS data. Top: mid-resolution (≈9 Å binning) STIS transmission spectrum is shown with red points. Bottom: low-resolution (≈74 Å binning) STIS and UVIS transmission spectra are shown with the red circles and green triangles, respectively. In both panels, the blue and black lines show the LTE and NLTE synthetic transmission spectra, respectively, while the grey horizontal dashed line indicates the broad-band Rp/Rs obtained from the EulerCAM optical transits (Hellier et al. 2019). The horizontal error bars shown for the red and green points correspond to the size of the wavelength bins.

Fig. 15 Comparison between synthetic and observed hydrogen Balmer line profiles. Top: NLTE (black) and LTE (red) theoretical transmission spectra compared to the observed Hα line profile (blue dots) presented by Damasceno et al. (2024). The green line is a Gaussian fit to the observed line profiles. The black and red dashed lines correspond, respectively, to the NLTE and LTE line profiles convolved with a Gaussian profile with a width such that the height of the lines equals the height of the Gaussian fit to the observations. Middle: same as top panel but for the Hβ line. Bottom: comparison of the observed Hα line (blue dots) with NLTE synthetic spectra computed with the system parameters listed in Table 2 (black; same as in the top panel), with Mp = 2.0 MJ (red) and 2.5 MJ (magenta), and with Teff = 8600 K (orange).

5.4 Comparison with optical observations

Figure 15 shows the comparison of the synthetic NLTE and LTE transmission spectra with the high-resolution observations of the Hα and Hβ lines presented by Damasceno et al. (2024). The NLTE synthetic line profiles are significantly stronger than the observations. Instead, the line profiles obtained assuming LTE appear to be a better match to the data, particularly following the inclusion of some additional broadening.

Previous works have shown that for the UHJs KELT-9b and MASCARA-2b/KELT-20b, the NLTE synthetic transmission spectra provide a significantly better match to the observed H I Balmer lines compared to the LTE models (Fossati et al. 2021, 2023). In terms of system properties, WASP-178b shares several similarities with KELT-9b and MASCARA-2b/KELT20b, and thus the results shown in Figure 15 are unexpected and call for a deeper investigation.

We first attempted to increase planetary mass or decrease stellar T_eff as both decrease the strength of the Hα line (see e.g. Fossati et al. 2023). The bottom panel of Figure 15 shows a comparison between the observed Hα line and NLTE synthetic spectra computed for M_p = 2.0 and 2.5 M_J, and for T_eff = 8600 K. As expected, a higher planetary mass decreases the strength of the Hα line, but fitting the observation would require a mass higher than 2.5 M_J that is inconsistent at more than 8σ with the measured stellar and planetary reflex motions (Hellier et al. 2019; Cont et al. 2024). Similarly, considering a 500 K cooler input stellar SED decreases the strength of the Hα line, but fitting the observation would require reducing T_eff to values inconsistent with the observed spectrum.

Having discarded uncertainties on the basic system parameters as origin of the discrepancy, we turn our attention to the planetary atmosphere as the weaker observed Balmer lines compared to the NLTE model suggest a cooler TP profile, particularly around the line formation region. Considering that by dropping the LTE assumption the NLTE model can be deemed being physically more accurate than the LTE one, the explanation for this mismatch can be sought by looking for planetary atmospheric properties or physical processes unaccounted for in the model that would lead to a decrease in the atmospheric temperature.

As can be inferred from Figures 3, 4, and 5, the Mg and/or Fe abundance can significantly impact the shape of the TP profile. Therefore, a lower middle and upper atmospheric temperature can be obtained by decreasing the Fe abundance, which decreases the heating, and/or by increasing the Mg abundance, which increases the cooling. A low Fe abundance would be in contrast with the clear detection of both FeI and Fe II by Damasceno et al. (2024). However, a high Mg abundance would be in line with the detection of Mg I by (Damasceno et al. 2024). Nevertheless, our NLTE transmission spectrum computed assuming solar abundances fits the UVIS observations well. Therefore, varying the atmospheric abundance pattern might not significantly improve the fit of the NLTE synthetic spectrum to the available NUV and optical observations. In relation to the planetary abundance pattern, we remark that in the case of the WASP-178 system the stellar photospheric abundances cannot be taken as reference for inferring planetary properties. This is because in Am stars (and generally in all intermediate-mass chemically peculiar stars), the peculiar abundance pattern is a characteristic of the stellar atmospheric layers, while the bulk stellar abundance pattern is different and can hardly be accurately constrained for field stars (see e.g. Fossati et al. 2008, 2011a).

The most notable difference between WASP-178b and the other two UHJs is the lower planetary mass of WASP-178b and thus the lower planetary gravity. As already mentioned in Section 5.2, the lower gravity moves towards higher pressures the level at which hydrodynamic escape (not accounted for in our modelling) becomes important. This might significantly impact the TP profile and the strength of the absorption lines, especially of the H I Balmer lines that form in the 1−10 nbar range, that is the pressure level at which hydrodynamic effects are expected to become important (see Section 5.1). In particular, with the inclusion of hydrodynamic effects, the temperature in the upper atmosphere would decrease as a consequence of adiabatic cooling, leading to a decrease in the strength of the Balmer lines. However, the inclusion of hydrodynamics in the modelling might also lead to strengthen the lines that (partly) form in the upper atmosphere as a result of advection. Therefore, although the presence of hydrodynamic effects in the upper atmosphere might explain the discrepancy between the NLTE synthetic spectrum and the observed Balmer line profiles, dedicated hydrodynamic modelling that properly accounts for metal line heating and cooling processes as well as advection (e.g. Huang et al. 2023; Kubyshkina et al. 2024) is needed to confirm or reject this hypothesis. Furthermore, dedicated hydrodynamic modelling would clarify whether the significant line broadening measured by Damasceno et al. (2024) is indeed related to atmospheric escape or not.

6 Conclusions

We have performed forward NLTE atmospheric modelling of the UHJ WASP-178b (Hellier et al. 2019). To identify the most adequate stellar SED to use as input for the model, we conducted a stellar spectroscopic analysis from which we derived an effective temperature of 9100 ± 100 K and confirmed that the stellar atmosphere is characterised by Am chemical peculiarities. Therefore, the star does not present a strong XUV emission (i.e. only photospheric emission). Meanwhile, the metallicity, and more in general the abundance pattern, cannot be taken as reference for the system because Am peculiarities are the result of chemical diffusion processes localised to the stellar atmosphere.

We conducted planetary atmospheric modelling combining TP profiles computed with HELIOS (in LTE) for the lower atmosphere and ClOUDY (in both LTE or NLTE) for the middle and upper atmosphere, following the same approach previously used to model the atmospheres of the UHJs KELT-9b and MASCARA-2b/KELT-20b (Fossati et al. 2021, 2023). Ultimately, we obtained an isothermal TP profile at pressures higher than 10 mbar and lower than 10⁻⁸ bar, with a roughly linear rise from ∼2200 K to ∼8100 K in between. A detailed analysis of the atmospheric thermal balance shows that accounting for NLTE effects significantly increases the middle and upper atmospheric temperature. This is due to a combination of increased heating due to the overpopulation of long-lived FeII levels presenting transitions mainly in the NUV band, where the stellar emission is strong, and of decreased cooling due to the underpopulation of MgI and MgII levels that dominate radiative cooling. We also identified that hydrogen line cooling, which is the second most important coolant after Mg, acts as a thermostat because it does not allow the atmospheric temperature to go above 9000 K independently of the Mg level underpopulation and abundance. This shows that NLTE effects, and particularly the resulting increased Fe line heating and decreased Mg line cooling, cannot be ignored when modelling hot planetary atmospheres.

Despite the hot upper atmosphere and the rather low planetary mass, the lack of strong stellar XUV emission implies that the atmosphere remains mostly hydrostatic, up to pressures of ∼1 nbar. We used the NLTE model to constrain the planetary mass-loss rate, for which we obtained a value of 4.5 × 10¹⁰ g s⁻¹, which is comparable to that of classical hot Jupiters such as HD 209458b. We computed the planetary synthetic transmission spectrum covering from the far-UV to the near-infrared, obtaining that a number of strong UV lines extend to pressures lower than 1 nbar and even beyond the Roche lobe. Therefore, more sophisticated models that also account for hydrodynamics (e.g. Huang et al. 2023; Kubyshkina et al. 2024) are required to more adequately constrain the mass-loss rate and model the strongest lines, particularly at high spectral resolution.

We compared the synthetic transmission spectra computed in LTE and NLTE with the available observations at optical (Hα and Hβ, from Damasceno et al. 2024; high spectral resolution) and NUV (HST STIS and WFC3/UVIS from Lothringer et al. 2022; medium and low spectral resolution, respectively) wavelengths. We find that the NLTE transmission spectrum overestimates the observed Hα and Hβ absorption, while the LTE model is more in line with the data. This is surprising given that the NLTE models of the hydrogen Balmer lines computed for UHJs KELT-9b and MASCARA-2b/KELT-20b match the observations remarkably well (and much better than the LTE models). We could not identify the origin of this discrepancy, and we suggest that because of the hydrodynamic nature of the upper atmosphere where the Hα and Hβ lines form, particularly at high spectral resolution, future attempts to reproduce the observations should account for hydrodynamics in addition to NLTE effects. We find that the NLTE transmission spectrum is an excellent match to the low resolution UVIS observation, particularly in the NUV band, without the need to consider SiO absorption or super-solar metallicity, as required when assuming LTE instead (Lothringer et al. 2022, 2025). This is mainly the outcome of the hotter NLTE TP profile compared to the LTE one. We performed an independent analysis of the STIS observations aimed at minimising the impact of systematic noise, but we obtained that red noise is too strong, and the resulting transmission spectrum is not consistent with the UVIS one. Furthermore, the resulting transmission spectrum is significantly dependent on the applied data analysis method. Therefore, the STIS spectrum cannot be used to support the detection of SiO absorption presented by Lothringer et al. (2022), and thus SiO may have condensed in the atmosphere of WASP-178b (Moran et al. 2024; Lothringer et al. 2025), but Fe- and Mg-bearing silicates did not. Furthermore, our results indicate that SiO absorption might not be present at all in the NUV transmission spectrum of WASP178b, which would indicate a significantly lower Si abundance compared to that derived by Lothringer et al. (2025) with consequences on inferences with respect to comparisons with other UHJs and possible formation scenarios. Therefore, accounting for NLTE effects in the computation of the TP profile and transmission spectrum might be critical for making the best use of the available observations to advance our understanding of the hottest planets.

Acknowledgements

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. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. LF and SB gratefully acknowledge funding from the European Union’s Horizon Europe research and innovation programme under grant agreement No. 101079231 (EXOHOST), and from UK Research and Innovation (UKRI) under the UK government’s Horizon Europe funding guarantee (grant number 10051045). AGS gratefully acknowledge funding from the Austrian Science Fund (FWF) grant J4596-N. D.S. acknowledges financial support from the project PID2021-126365NB-C21(MCI/AEI/FEDER, UE) and from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033. M.E.Y. acknowledges support from the UK Science and Technology Facilities Council (STFC, grant number ST/X00094X/1). We thank the anonymous referee for the useful comments that led to improve the paper.

Appendix A Additional figures and tables

Fig. A.1 Comparison between the composite TP profiles computed in NLTE (black) and LTE (green) with those obtained by Lothringer et al. (2022, LTE) employing forward modelling (red) and retrieval (blue), Cont et al. (2024) employing retrieval (yellow), and Lothringer & Barman (2019) employing retrieval (magenta).

Fig. A.2 Same as Figure 8 but for the 1100−1550 Å band.

Fig. A.3 Same as Figure 8 but for the 1500−2350 Å band.

Fig. A.4 Same as Figure 8 but for the 2300−3050 Å band.

Fig. A.5 Same as Figure 8 but for the 3000–4050 Å band.

Fig. A.6 Same as Figure 8 but for the 4000−6100 Å band.

Fig. A.7 Same as Figure 8 but for the 6000−11000 Å band.

Fig. A.8 Transit depth difference between the NLTE and LTE transmission spectra shown in Figure 8. The top plot covers the UV and optical range, while the bottom plot covers the infrared band.

Appendix B HST STIS data analysis

Before proceeding with the re-analysis of the STIS observations, we derived precise transit ephemerides by analysing the available TESS photometry, which has been obtained during sectors 11 (2019), 38 (2021), and 65 (2023). We modelled all transit events present in the data using ALLESFITTER (Günther & Daylan 2021), while simultaneously detrending the light curves with splines. We computed the posterior distributions of the parameters of interest using the Dynamic Nested Sampling inference algorithm (see e.g. Feroz & Hobson 2008; Feroz et al. 2019) that is integrated in ALLESFITTER via the DYNESTY package (Speagle 2020). Finally, we obtained an orbital period $P={3.34483532}_{-7.7\times {10}^{-7}}^{+7.6\times {10}^{-7}}$ days and a reference mid-transit time T_0,BJD = 2459348.73196±0.00013 days, which agree within 1σ to the values adopted by Lothringer et al. (2022).

The STIS data are affected by significant systematics. We tackled their removal using two techniques, both based on polynomial detrending. The first technique is described in Cubillos et al. (2020), Cubillos et al. (2023), and Sreejith et al. (2023) and builds up on the jitter detrending scheme introduced by Sing et al. (2019). We corrected for systematic effects the order-byorder white light curves using a function S_λ(t, ϕ, j_i) based on time t, telescope phase ϕ, and the jitter parameter that exhibited the highest degree of correlation with the data j_i. The transit model T_λ(t) has been implemented following Mandel & Agol (2002); the transit ephemerides were derived from the TESS data as described above and the limb-darkening coefficients were calculated through the open-source routines of Espinoza & Jordán (2015). In detail, the model F_λ(t) applied to the raw light curves is $$F_{\lambda }(t)=T_{\lambda }(t)S_{\lambda }(t,\varphi ,j_i).$$(B.1)

This approach, which we call ‘one jitter parameter’ method, assumes that within each spectral order the systematics are wavelength independent. It turns out that the optimal model of the instrumental systematics takes the form $$\begin{array}{rl}& S_{\lambda }(t,\varphi ,j_i)=1+a_0(t-t_0)+a_1{(t-t_0)}^2+a_2{(t-t_0)}^3\\ & +b_0(\varphi -{\varphi }_0)+b_1{(\varphi -{\varphi }_0)}^2+b_2{(\varphi -{\varphi }_0)}^3+b_3{(\varphi -{\varphi }_0)}^4\\ & +c_0(j_i-⟨j_i⟩)\end{array}$$(B.2)

where a_k, b_k, and c_k are the polynomial coefficients of the fit, and t₀ and ϕ₀ are the reference values for time and phase, set as midtransit time t₀ = T₀ and the HST telescope mid-phase ϕ₀ = 0.2, respectively. The term 〈j_i〉 denotes the mean value of the jitter parameter j_i along the visit. The best fit in this case corresponds to the jitter parameter LOS_Zenith.

To select the optimal systematic model, we used the Bayesian Information Criterion (BIC) and Akaike Information Criterion (AIC). Once the optimal systematic model has been identified, we extracted the transmission spectra using the ‘dividewhite’ spectral analysis method (e.g. Cubillos et al. 2020, 2023; Sreejith et al. 2023), which corrects instrumental systematics by dividing white-light curves by the best-fit transit model, assuming weak wavelength dependence. Systematics-corrected light curves are constructed by propagating uncertainties through all steps and dividing the data into wavelength bins. These wavelength-specific light curves are modelled with a Mandel & Agol (2002) transit model, fixing the orbital parameters and the limb-darkening coefficients. The free parameters for each fit are transit depth, mid-transit time, and visit-specific out-of-transit flux. Remarkably, the BIC criterion led us using a 3rd order polynomial in time (BIC = 3313.4), while for example when using a 1st order polynomial in time we obtained a higher BIC of 3637.5. This is further addressed later in this section.

However, this approach led to a transmission spectrum that was incompatible with both the UVIS data and the broad-band EulerCAM photometry (Figure B.1), indicating that the model did not adequately capture the underlying systematics influencing the data. Therefore, additional factors not accounted for in our initial model (Equation B.2) contribute significantly to the observations. To better understand the root cause of these discrepancies, we calculated the correlation between the observed data and a wide range of parameters, including time and the HST jitter parameters available through the data distributed by the HST archive. The result of the correlation analysis are summarised in Figure B.2, which identifies Termang, LOS_Zenith, Latitude, Longitude, MagV2, MagV3 as the parameters with the highest correlation. As we found the correlation to be uniform across the covered wavelength range, we chose to use the same correlation parameters across all orders. Notably, these parameters show a degree of correlation with the data stronger than that shown by time.

Fig. B.1 Comparison of the NLTE (black) theoretical transmission spectrum with the HST WFC3/UVIS observations published by Lothringer et al. (2022, red) and the STIS transmission spectrum (green) re-analysed employing a detrending based on one jitter parameter and 3rd order polynomial in time (see text), binned to the same binning of the UVIS spectrum. The horizontal error bars shown for the red and green points correspond to the size of the wavelength bins. The grey horizontal dashed line indicates the broad-band Rp/Rs obtained from EulerCAM (Hellier et al. 2019).

To account for these findings, the second technique we employed to analyse the data, which we call ‘multiple jitter parameter’ method, enables one to account for several parameters simultaneously against which to decorrelate. Therefore, we considered the parameters showing the strongest degree of correlation mentioned above as input to the decorrelation model based on the MCMCI code¹² (Bonfanti & Gillon 2020) that we briefly describe below. The MCMCI is a Fortran-based fitting code that integrates a Markov Chain Monte Carlo (MCMC) algorithm with an isochrone placement algorithm to fully characterise exoplanetary systems and supports analysis of light curves and radial velocity (RV) time series by modelling transits, occultations, and the Rossiter-McLaughlin effect. It samples the parameter space using the Metropolis-Hastings algorithm and Gibbs sampling, to infer system properties such as planetary radius, mass, orbital eccentricity, and host star characteristics. The code is capable of correcting for observational trends and systematics in the data. Bayesian priors are incorporated and the code allows for flexible handling of the orbital eccentricity to reduce biases. The code selects optimal models based on the BIC, enabling precise characterisation of single and multi-planet systems. This tool can take up to six parameters (J_i) at the same time against which to decorrelate and employs the following function of time t to model the transit light curve $$F_{\lambda }(t)=T_{\lambda }(t)S_{\lambda }(t,J_1,J_2,J_3,J_4,J_x,J_y),$$(B.3)

where T_λ(t) is the Mandel & Agol (2002) transit model. The S_λ() term is a polynomial model (with order not greater than 4) of the instrumental systematics taking the form $$\begin{array}{rl}{S}_{\lambda }=& \sum _{i=1}^{n_t}{a}_i\cdot t^i+\sum _{i=1}^{n_{J_1}}{b}_i\cdot J_1{}^i+\sum _{i=1}^{n_{J_2}}{c}_i\cdot J_2{}^i+\sum _{i=1}^{n_{J_3}}{d}_i\cdot J_3{}^i+\\ & \sum _{i=1}^{n_{J_4}}{e}_i\cdot J_4{}^i+\sum _{i=0}^{n_{J_{xy}}}\sum _{j=0}^{n_{J_{xy}}-i}{f}_{ij}\cdot J_x^i\cdot J_y^j,\end{array}$$(B.4)

where a_i, b_i, c_i, d_i, and e_i are the polynomial coefficients of the fit for time and the jitter parameters (J₁−J₄; i.e. TermAng, LOS_Zenith, MagV2, MagV3), while the terms n_t and n_Ji indicate the order of the polynomials used for each parameter. The f_ij are the polynomial coefficients associated with the two jitter parameters J_x and J_y, which can have mutual correlation (in our case Latitude and Longitude).

Fig. B.2 Heat map showing the correlation of jitter parameters and time with flux across each spectral order. The x-axis represents the spectral order, while the y-axis lists the parameters. The colour indicates the strength of the correlation, with red and blue representing positive and negative correlations, respectively, as indicated by the colour bar on the right hand side.

We generated order-specific white light curves and conducted model fitting using MCMCI for various polynomial combinations of the jitter parameters. To ensure consistency across spectral orders, the same polynomial configuration was applied for all orders during the analysis. The optimal model (i.e. the sequence of polynomial orders n_t, n_J1, n_J2, n_J3, n_J4, and n_Jxy) was identified by evaluating the BIC, with the combination (3rd order for J₁ and J₂, 1st order for J₄, and 4th order for J_xy) yielding the lowest BIC selected as the best-fit model (BIC = 5667.84 as compared to 68570.7 for an all-zero polynomial baseline). Based on this optimal solution, we constructed a systematics model using the best-fit parameters to account for instrumental and observational trends. Then, we applied the ‘divide-white’ approach, as described earlier, to extract the transmission spectrum for the wavelength bins of interest.

Figure B.3 shows example light curves for three spectral orders located around the beginning, middle, and end of the wavelength range covered by the STIS observations. In particular, this plot compares the raw light curves (left column), which clearly show the presence of strong systematic noise, with the light curves obtained using the ‘one jitter parameter’ (high and low polynomial order in time; second and third columns, respectively) and ‘multiple jitter parameter’ (last column) analysis methods described above. The light curves obtained using the ‘one jitter parameter’ analysis method return either an unrealistic transmission spectrum (i.e. for the high polynomial order in time; see Figure B.1) or light curves that appear to be still affected by residual systematic noise (i.e. for the low polynomial order in time), in addition to be disfavoured by the BIC. Instead, the light curves obtained applying the ‘multiple jitter parameter’ analysis method present significantly less systematic noise. This demonstrates the higher quality of the ‘multiple jitter parameter’ analysis method with respect to the ‘one jitter parameter’ analysis method.

Fig. B.3 Light curves extracted from the STIS data for three spectral orders located around the beginning (21st spectral order; 2326.7–2365.8 Å; top), middle (11th spectral order; 2629.1-2673.3 Å; middle), and end (2nd spectral order; 2977.4–3027.3 Å; bottom) of the range covered by the data. Within each plot, the left panel is for the raw light curve, the two middle panels are for the data analysed using the ‘one jitter parameter’ (high and low polynomial order in time) analysis method, and the right panel is for the data analysed using the ‘multiple jitter parameter’ analysis method.

Figure B.4 compares the STIS transmission spectra of WASP-178b obtained using the different systematics correction procedures described above, with that presented by Lothringer et al. (2022, priv. comm. from the authors). Figure B.4 further shows the re-analysed STIS transmission spectra extracted binning at the resolution of UVIS and the observed UVIS transmission spectrum. This plot clearly shows that using just one jitter parameter with the higher polynomial order in time introduces a bias in the transmission spectrum resulting in much deeper transits, possibly due to over-fitting of the data (see also Figure B.1). Artificially removing the correlation with time or limiting it to a first order polynomial, either considering one or multiple jitter parameters results in a transmission spectrum comparable to that presented by Lothringer et al. (2022) and more in line with what expected following the broad-band transit radius. Remarkably, the transmission spectrum obtained arbitrarily using a low polynomial order in time within the ‘one jitter parameter’ analysis method matches well the UVIS transmission spectrum, though this specific analysis method is disfavoured by the BIC (compared to the one obtained considering a high polynomial order in time) and by the presence of residual systematic noise in the light curves.

Following the comparisons shown in Figure B.4 and the higher degree of sophistication of the second technique in removing systematic noise (Figure B.3), we consider the transmission spectrum obtained using the ‘multiple jitter parameter’ analysis method as the final one. However, the analysis and comparisons presented here clearly indicate that the results obtained from this data set are significantly model-dependent.

Fig. B.4 Comparison of STIS transmission spectra of WASP-178b obtained following different data analysis procedures. The black empty diamonds show the transmission spectrum obtained using the ‘multiple jitter parameter’ method (‘Multiple JP’). The blue and green empty diamonds show, respectively, the results obtained from the ‘one jitter parameter’ analysis methods considering a lower (‘1 JP lower order time’) and higher (‘1 JP higher order time’) polynomial order in time, with the latter being favoured by the BIC. For each of these, the filled squares show the transmission spectra obtained binning at UVIS resolution. For comparison, the red empty diamonds show the transmission spectrum presented by Lothringer et al. (2022, priv. comm. from the authors). The grey circles show the UVIS transmission spectrum from Lothringer et al. (2022). The horizontal grey dashed line indicates the broad-band Rp/Rs obtained from EulerCAM (Hellier et al. 2019).

References

All Tables

All Figures

Fig. 1 Ion abundance relative to the Sun (Asplund et al. 2009) of the WASP-178 atmosphere. On the y-axis, N corresponds to the number of atoms of a given ionic species indicated on the x-axis and legend, while Ntot is the total number of atoms. Black open squares show the abundance of the neutral elements, red open diamonds are for singly ionised elements, and blue asterisks are for doubly ionised elements.

In the text

Fig. 2 Theoretical atmospheric structure obtained for WASP-178b. Top: HELIOS (dashed orange line) and CLOUDY (NLTE; dashed cyan line) TP profiles. The solid black line shows the composite TP profile. The horizontal dotted black line indicates the location of the continuum according to the HELIOS model. The horizontal dash-dotted red line gives the location of CLOUDY’s upper-density limit of 1015 cm−3. The dashed dark green line shows the CLOUDY TP profile computed assuming LTE. Bottom: Pressure (black; left y-axis) and temperature (red; right y-axis) composite theoretical profiles as a function of the planetary polar radius. For context, the L 1 point lies at about 4.00 Rp, while the Roche radius in the limb direction lies at about 2.67 Rp.

In the text

	Fig. 3 Comparison between the NLTE CLOUDY TP profiles obtained considering all elements up to Zn (solid black line) and excluding Fe (solid red line), Mg (solid blue line), Ca (solid green line that is almost overlapping with the black line), or photoionisation (yellow). The dashed orange line shows the HELIOS (i.e. LTE) TP profile for reference.
In the text

Fig. 4 Heating and cooling contributions in the middle and upper atmosphere (i.e. at pressures lower than CLOUDY’s density limit) of WASP-178b. Top: Contribution to the total heating as a function of pressure. At each pressure bin, the plot shows the three most important heating processes. The main heating processes are hydrogen photoionisation (red; photoionisation of HI lying in the ground state), metal line absorption (blue), H− absorption (green), photoionisation of hydrogenic species (magenta), collisions with H2 (brown), and Compton heating (i.e. electron absorption; yellow). Bottom: Contribution to the total cooling as a function of pressure. At each pressure bin, the plot shows the three most important cooling agents. The main cooling agents are Mg (red), Ca (blue), C (green), ${\mathrm{H}}_2^{+}$ (magenta), H (brown), S (orange), HD (i.e. hydrogen deuteride; yellow), Na (grey), and Si (light blue). Metal line absorption (primarily from FeII) and Mg line emission are the primary heating and cooling mechanisms, respectively.

In the text

	Fig. 5 Temperature-pressure profiles obtained from isolating the impacts of Fe (left) and Mg (right) ions. The black line is for when all elements are accounted for (same as in Figures 2 and 3), while the other lines were obtained considering that in the planetary atmosphere, Fe or Mg exist only in the form of a neutral species (blue) or a singly ionised species (green) or is absent (red).
In the text

	Fig. 6 Density (relative to the total density of hydrogen)-pressure profiles for neutral hydrogen (HI; solid black line), protons (HII; red), molecular hydrogen (H2.; dark green), ${\mathrm{H}}_2^{+}$ (blue), ${\mathrm{H}}_3^{+}$ (magenta), H− (orange), and electrons (e−; dashed black).
In the text

Fig. 7 Atmospheric mixing ratios obtained for metals. Top: mixing ratios for hydrogen (solid black line), H2 (dotted black), He (red), C (blue), O (dark green), Na (magenta), K (orange), and electrons (bright green). Bottom: same as the top panel but for Mg (red), Si (blue), Ca (dark green), and Fe (magenta). The hydrogen, H2, and e− mixing ratios are shown in both panels for reference. Neutral (XI), singly ionised (XII), and doubly ionised (XIII; present just in the bottom panel) species are shown as solid, dashed, and dash-dotted lines, respectively.

In the text

Fig. 8 Comparison between the theoretical LTE (red) and NLTE (black) transmission spectra computed at a resolution of 100 000 in the UV and optical range (top) and in the near-infrared (bottom). Within each plot, the bottom panel shows the relative difference between LTE and NLTE (in percent). In the top panel, the blue horizontal dashed line indicates the location of the Roche radius in the planetary limb direction (i.e. perpendicular to the star-planet connecting line).

In the text

Fig. 9 Comparison between theoretical transmission spectra computed with full NLTE (A; black; same as the black line in Figure 8), namely the NLTE TP profile and accounting for NLTE effects for the computation of the transmission spectrum (TS), and with partial NLTE (B, red), namely NLTE TP profile and assuming LTE for the computation of the transmission spectrum, at a resolution of 100 000 in the NUV and optical range. The bottom panel shows the relative difference between the two spectra in percent. In the top panel, the blue horizontal dashed line indicates the location of the Roche radius in the planetary limb direction.

In the text

Fig. 10 Line formation analysis in transmission geometry. Top: relative peak absorption for Hα (black), the MgIb line at ∼5183 Å (red), the Cal line at ∼6439 Å (blue), the Cair K line (green), the FeI line at ∼4071 Å (magenta), and the FeII line at ∼5169 Å (orange). Bottom: relative absorption in a 100 Å spectral bin centered at the MgI resonance line at ∼2850 Å (black), the MgII h&k resonance lines at ∼2800 Å (red), the FeI NUV absorption band at ∼2720 Å (blue), the FeII resonance line at ∼2600 Å (green), and the FeII NUV absorption band at ∼2370 Å (magenta) as a function of pressure. In both panels, the dashed line shows the TP profile (top x-axis). The line formation region corresponds to the pressure values in which the relative absorption ranges between zero and one.

In the text

	Fig. 11 Atmospheric sound speed (black; left y-axis) and Jeans escape parameter (red; right y-axis) profiles as a function of pressure computed in the direction of the sub-stellar point (i.e. the Roche lobe lies at the L1 point).
In the text

	Fig. 12 Comparison between the theoretical NLTE TP profiles of WASP-178b (black; this work), MASCARA-2b/KELT-20b (red; from Fossati et al. 2023), and KELT-9b (blue; from Fossati et al. 2021). For context, the legend also gives the effective temperature of the host stars from Section 2 (for WASP-178b), Talens et al. (2018, for MASCARA-2b/KELT-20b), and Borsa et al. (2019, for KELT-9b).
In the text

Fig. 13 Comparison of the NLTE (black) and LTE (blue) theoretical transmission spectra with the WFC3/UVIS observations published by Lothringer et al. (2022, red). The horizontal error bars shown for the red points correspond to the size of the wavelength bins. The top plot is without scaling of the models, while in the bottom plot the NLTE and LTE models have been scaled to match the average of the UVIS observations at wavelengths longer than 4050 Å, marked by the vertical grey dashed line. In the top plot, the grey horizontal dashed line indicates the broad-band Rp/Rs obtained from the EulerCAM optical transits (Hellier et al. 2019). In each plot, the bottom panel shows the residuals between the observation and each of the models, with a green dashed line at zero for reference.

In the text

Fig. 14 Comparison between observed and synthetic transmission spectra in the region covered by the STIS data. Top: mid-resolution (≈9 Å binning) STIS transmission spectrum is shown with red points. Bottom: low-resolution (≈74 Å binning) STIS and UVIS transmission spectra are shown with the red circles and green triangles, respectively. In both panels, the blue and black lines show the LTE and NLTE synthetic transmission spectra, respectively, while the grey horizontal dashed line indicates the broad-band Rp/Rs obtained from the EulerCAM optical transits (Hellier et al. 2019). The horizontal error bars shown for the red and green points correspond to the size of the wavelength bins.

In the text

Fig. 15 Comparison between synthetic and observed hydrogen Balmer line profiles. Top: NLTE (black) and LTE (red) theoretical transmission spectra compared to the observed Hα line profile (blue dots) presented by Damasceno et al. (2024). The green line is a Gaussian fit to the observed line profiles. The black and red dashed lines correspond, respectively, to the NLTE and LTE line profiles convolved with a Gaussian profile with a width such that the height of the lines equals the height of the Gaussian fit to the observations. Middle: same as top panel but for the Hβ line. Bottom: comparison of the observed Hα line (blue dots) with NLTE synthetic spectra computed with the system parameters listed in Table 2 (black; same as in the top panel), with Mp = 2.0 MJ (red) and 2.5 MJ (magenta), and with Teff = 8600 K (orange).

In the text

	Fig. A.1 Comparison between the composite TP profiles computed in NLTE (black) and LTE (green) with those obtained by Lothringer et al. (2022, LTE) employing forward modelling (red) and retrieval (blue), Cont et al. (2024) employing retrieval (yellow), and Lothringer & Barman (2019) employing retrieval (magenta).
In the text

	Fig. A.2 Same as Figure 8 but for the 1100−1550 Å band.
In the text

	Fig. A.3 Same as Figure 8 but for the 1500−2350 Å band.
In the text

	Fig. A.4 Same as Figure 8 but for the 2300−3050 Å band.
In the text

	Fig. A.5 Same as Figure 8 but for the 3000–4050 Å band.
In the text

	Fig. A.6 Same as Figure 8 but for the 4000−6100 Å band.
In the text

	Fig. A.7 Same as Figure 8 but for the 6000−11000 Å band.
In the text

	Fig. A.8 Transit depth difference between the NLTE and LTE transmission spectra shown in Figure 8. The top plot covers the UV and optical range, while the bottom plot covers the infrared band.
In the text

Fig. B.1 Comparison of the NLTE (black) theoretical transmission spectrum with the HST WFC3/UVIS observations published by Lothringer et al. (2022, red) and the STIS transmission spectrum (green) re-analysed employing a detrending based on one jitter parameter and 3rd order polynomial in time (see text), binned to the same binning of the UVIS spectrum. The horizontal error bars shown for the red and green points correspond to the size of the wavelength bins. The grey horizontal dashed line indicates the broad-band Rp/Rs obtained from EulerCAM (Hellier et al. 2019).

In the text

	Fig. B.2 Heat map showing the correlation of jitter parameters and time with flux across each spectral order. The x-axis represents the spectral order, while the y-axis lists the parameters. The colour indicates the strength of the correlation, with red and blue representing positive and negative correlations, respectively, as indicated by the colour bar on the right hand side.
In the text

Fig. B.3 Light curves extracted from the STIS data for three spectral orders located around the beginning (21st spectral order; 2326.7–2365.8 Å; top), middle (11th spectral order; 2629.1-2673.3 Å; middle), and end (2nd spectral order; 2977.4–3027.3 Å; bottom) of the range covered by the data. Within each plot, the left panel is for the raw light curve, the two middle panels are for the data analysed using the ‘one jitter parameter’ (high and low polynomial order in time) analysis method, and the right panel is for the data analysed using the ‘multiple jitter parameter’ analysis method.

In the text

Fig. B.4 Comparison of STIS transmission spectra of WASP-178b obtained following different data analysis procedures. The black empty diamonds show the transmission spectrum obtained using the ‘multiple jitter parameter’ method (‘Multiple JP’). The blue and green empty diamonds show, respectively, the results obtained from the ‘one jitter parameter’ analysis methods considering a lower (‘1 JP lower order time’) and higher (‘1 JP higher order time’) polynomial order in time, with the latter being favoured by the BIC. For each of these, the filled squares show the transmission spectra obtained binning at UVIS resolution. For comparison, the red empty diamonds show the transmission spectrum presented by Lothringer et al. (2022, priv. comm. from the authors). The grey circles show the UVIS transmission spectrum from Lothringer et al. (2022). The horizontal grey dashed line indicates the broad-band Rp/Rs obtained from EulerCAM (Hellier et al. 2019).

In the text