Atmosphere Escape Inferred from Modeling the Hα Transmission Spectrum of WASP-121b

We used the hydrodynamic model (Yan & Guo 2019) to simulate the atmospheric structure of WASP-121b and obtained the atmospheric temperature, velocity, and particle number densities. The planetary and stellar parameters are based on the observations (Delrez et al. 2016; Cabot et al. 2020). The equilibrium temperature (T_eq) is 2361 K, which is also the temperature at the bottom boundary in our model. The high temperature hints at the dissociation of H₂. The chemical composition of WASP-121b is assumed to be the same as that of WASP-121, which is calculated by the solar abundance that is modified by [Fe/H] = 0.13 (Delrez et al. 2016). The integrated flux in the XUV band is an important input in the simulations. Due to the lack of the stellar XUV observations, we used the age–luminosity relation (Sanz-Forcada et al. 2011) to calculate the F_XUV received by the planet. WASP-121b is about 1.5 Gigayear (Gyr), so the F_XUV is about 37387 erg cm⁻² s⁻¹ at the orbital distance. In our calculation the value is divided by a factor of 4, which accounts for the uniform redistribution of the stellar radiation energy around the planet. Finally, the XUV spectral energy distribution (SED) is obtained by the XSPEC-APEC software (Arnaud 1996). In the simulations the pressure at the bottom boundary is 1 μbar. The upper boundary is 7.6 R_p, which covers the radius of the host star. In this Letter, the above input values are called the fiducial ones. In general, the escaping models assume that the photons of Lyα can freely escape from the atmosphere (Murray-Clay et al. 2009). In the process of resonant scattering, the number of scattering that a Lyα photon takes to escape the atmosphere is comparable to its line center optical depth τ. In the atmosphere above ∼1.1–1.2 R_p where Lyα cooling is most efficient, τ ≪ 1/p_abs, where p_abs is the Lyα photon destroy probability per scattering (Huang et al. 2017). Thus, the Lyα cooling is included in the simulations. Furthermore, the stellar tidal force is also considered in the model.

The Hα absorption is caused by the hydrogen atoms in the first excited state (n = 2, where n is the principal quantum number). Because of the coupling of spin and orbital angular momentum, this state is split to 2s and 2p substate. In the upper thermosphere, both the collisions among particles and the radiation process affect the population of H(2). Thus, we used a non-local thermal equilibrium scheme to calculate the hydrogen populations based on the hydrodynamic results. Assuming that the atmosphere is in a stationary state, the production rates of H(2) are equal to their loss rates. We find that the Lyα mean intensity (${\bar{J}}_{\mathrm{Ly}\alpha }$) is dominant in determining the number density of H(2p) (as shown by Equation (1)), which is consistent with that of Christie et al. (2013) and Huang et al. (2017). Therefore, it can be approximately expressed as

$$\begin{eqnarray}&&{n}_{2{\rm{p}}}\approx \displaystyle \frac{{B}_{1{\rm{s}}\to 2{\rm{p}}}{\bar{J}}_{\mathrm{Ly}\alpha }}{{A}_{2{\rm{p}}\to 1{\rm{s}}}}{n}_{1{\rm{s}}}\end{eqnarray} \tag{ 1 }$$

where n(1s) and n(2p) are the number densities of H(1s) and H(2p), B_1s→2p and A_2p→1s are the Einstein coefficients (Rybicki & Lightman 2004). The number density of H(2s) is solved in the meantime. The source of H(2s) is mainly from l-mixing; and the sink of H(2s) is mostly due to the l-mixing (<2.5 R_p) and photoionization (>2.5 R_p). Therefore, it can be approximately expressed as

$$\begin{eqnarray}&&{n}_{2{\rm{s}}}\approx \displaystyle \frac{{C}_{2{\rm{p}}\to 2{\rm{s}}({\rm{p}})}{n}_{2{\rm{p}}}{n}_{{\rm{p}}}+{C}_{2{\rm{p}}\to 2{\rm{s}}({\rm{e}})}{n}_{2{\rm{p}}}{n}_{{\rm{e}}}}{{C}_{2{\rm{s}}\to 2{\rm{p}}({\rm{p}})}{n}_{{\rm{p}}}+{C}_{2{\rm{s}}\to 2{\rm{p}}({\rm{e}})}{n}_{{\rm{e}}}+{{\rm{\Gamma }}}_{2{\rm{s}}}}\end{eqnarray} \tag{ 2 }$$

where C_2p↔2s(p) and C_2p↔2s(e) are the protons' and electrons' collisional transition rates (Seaton 1955) between 2p and 2s, n_p and n_e are the number densities of protons and electrons, and Γ_2s is the photoionization rate of H(2s). Although n_2p and n_2s can be estimated by Equations (1) and (2) in the simulations we solved the equation of rate equilibrium for H(2p) and H(2s) by including the reactions of radiative excitation and de-excitation, collision excitation and de-excitation, photoionization and recombination, etc. The equations of rate equilibrium are the same to that of Huang et al. (2017).

Because there is currently no available Lyα profile of the host star, in our model we took the Lyα flux of ζ Dor (Linsky et al. 2013) to replace that of WASP-121. ζ Dor is a F7V star whose spectral type is similar to that of WASP-121 (F6V). The Lyα profile of ζ Dor is referenced from Wood et al. (2005). Based on this, we constructed a double Gaussian profile with the full width at half maximum (FWHM) = 0.7 Å, and the two centers being 1215.5 and 1215.9 Å (Lyα₁, see the solid black line in Figure 1(a)). The Lyα integrated flux is about 71,900 erg cm⁻² s⁻¹. The ${\bar{J}}_{\mathrm{Ly}\alpha }$ is then calculated by the Lyα resonant scattering method of Huang et al. (2017). Because the stellar Lyα profile is an important physical input in calculating H(2p) population but cannot be measured precisely, we investigated another Lyα profile (Lyα₂, see the dashed red line in Figure 1(a)) with the same integrated flux and two centers but with a different FWHM = 0.45 Å. The Lyα intensity is higher around the line center of Lyα₂ profile, which leads to higher Voigt line profile weighted mean intensity (see Figure 1(b) for ${\bar{J}}_{\mathrm{Ly}\alpha }$ of Lyα₁ and Lyα₂). The Hα absorption (see Figure 1(c); for more details, see the Method section below) is slightly deeper for Lyα₂ case compared to that of Lyα₁. The difference of absorption at the Hα line center is less than 0.2% and it is indistinguishable at the line wings. Our results show that the FWHM of Lyα input profile in the range of 0.45–0.7 Å has minor influence on the final results. Therefore, we use Lyα₁ in our models.

Zoom In Zoom Out Reset image size

Finally, the ionization of H(2s) and H(2p) is an important process that is caused by the photons in the wavelength range of 912–3646 Å. The spectrum in this wavelength range is taken from the stellar atmosphere model of Castelli & Kurucz (2003). The photoionization cross-sections of H(2s) and H(2p) are cited from TOPbase of The Opacity Project (Cunto & Mendoza 1992; Cunto et al. 1993).

After we obtained the hydrogen populations, we simulated the Hα radiative transfer as the stellar Hα line travels along the ray path in the planetary atmosphere during transit (Yan & Guo 2019). The transmission spectrum is defined by Equation (3) as a function of wavelength,

$$\begin{eqnarray}&&\mathrm{TS}(\lambda )=\displaystyle \frac{{F}_{\mathrm{IT}}}{{F}_{\mathrm{OT}}}(\lambda )-1.0\end{eqnarray} \tag{ 3 }$$

and the excess absorption (not including the planet itself) is " − TS(λ)," where F_IT and F_OT are the in-transit and out-of-transit flux. Hα absorption is a bound-bound transition whose line center is at 6562.8 Å (in air). The calculation of cross section of Hα absorption is similar to that of Lyα (Yan & Guo 2019), and the oscillator strength is 0.64108 at 6562.8 Å taken from NIST Atomic Spectra Database.⁵

We define the model with the fiducial inputs as our fiducial model. The mass-loss rate of WASP-121b is 1.28 × 10¹² g s⁻¹ for the fiducial model. The number density of hydrogen atoms is obtained from the hydrodynamic simulation. The atmospheric structures are showed in Figures 2(a)–(c). As we can see, the highest temperature is higher than 10,000 K. The velocity becomes supersonic when the altitude is larger than 1.7 R_p, and reaches 100 km s⁻¹ beyond 7 R_p. By solving a detailed equation of statistical equilibrium, we obtained the number densities of H(2s) and H(2p) (see Figure 2(a)). The number density of H(2s) plus H(2p) is about a few 10⁻⁷ times H(1s). The ratio of H(2p) to H(2s) changes significantly with the increase of r/R_p, mainly due to the large photoionization of the n = 2 state hydrogen, which especially affects H(2s). The optical depth of Hα line center is shown in Figure 2(d), which is larger than unity when r/R_p is less than 1.1.

Zoom In Zoom Out Reset image size

We compared the simulated transmission spectrum with the observation of Borsa et al. (2021, hereafter B21).⁶ In their work, they reported the observations in 1-UT and 4-UT mode, the observation time of which is 2019 January 6 and 2018 November 30. The Hα absorption depth was about 1.4% and 1.7% for 1-UT and 4-UT mode, with a blueshift of 4.64 and 3.9 km s⁻¹, respectively. B20S is obtained by shifting the observed transmission spectrum of Borsa et al. (2020) toward the red side by the corresponding blueshift for 1-UT and 4-UT.

To investigate the contribution of different atmospheric regions to the Hα absorption, we calculated the Hα absorption for different altitudes of the atmosphere until the altitude reaches 7.6 R_p. Figure 3(a) shows the results of the fiducial model, in which the gray dots with error bars are B20S. The light-gray and dark-gray points are for 1-UT and 4-UT, respectively. The different lines represent the absorption of Hα produced in different atmospheric altitudes. Our simulations show that different altitudes of the atmosphere contribute differently to the final transmission spectrum of Hα. The increase of the altitude leads to a deeper absorption. From 3 R_p to 4 R_p, the absorption at Hα line center only increases by 0.05%. The Hα absorption of WASP-121b by the atmosphere above 4 R_p is negligible because the H(2) are sparse. In addition, we compared the equivalent width (EW) of the model transmission spectrum with that of the observations in the passbands 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, and 3.6 Å. We found that the EW will decrease with the increase of the passbands when the passbands are larger than 1.5 Å for 1-UT, due to some minus values of the absorption at the line wings. We also found that the EW calculated by our model is lower than that of the 4-UT if the passband is greater than 1.5 Å. Thus, we only analyzed the EW in passbands 0.75, 1.0, and 1.5 Å. Figures 3(b)–(d) show the equivalent width calculated in the three passbands as a function of the atmospheric altitude, respectively. The black line represents the fiducial model, in which the sonic point is marked by the purple cross. The orange and cyan horizontal solid lines represent the mean EW of 1-UT and 4-UT, respectively. The corresponding dashed lines are the upper (+1σ) and lower (−1σ) limit of the observations. One can see that the EW of the fiducial model increases with the increase of the atmospheric altitude. For 1-UT, an atmosphere higher than 1.7 Rp can fit the observation in the three passbands. For 4-UT, a supersonic atmosphere beyond the Roche lobe (1.86 R_p) is required to fit the EW in 0.75 and 1.0 Å. For example, an atmosphere at least up to 6 R_p can fit the lower limit of the observation in passband 1.0 Å. The EW of the model in passband 1.5 Å, however, cannot reach the lower limit of the 4-UT observation, because our model cannot reproduce the relatively high absorption at the Hα line wings.

Zoom In Zoom Out Reset image size

The results above show that the absorption of Hα close to the line center can be well fitted by our model and the contribution of supersonic regions cannot be neglected. However, the strong absorption in the wings of Hα for 4-UT should be investigated in the future. Furthermore, the 1D simulation cannot reproduce the blueshift found by B21. A blueshift or redshift has also been found in other absorption lines in exoplanetary atmosphere (Allart et al. 2018; Casasayas-Barris et al. 2018; Bourrier et al. 2020; Cabot et al. 2020; Gibson et al. 2020), which can be led by the atmospheric winds or the circulations from the day-side to night-side (Guo 2013).

Because the XUV integrated flux was calculated by the age–luminosity relation of Sanz-Forcada et al. (2011), there could be some uncertainties due to the uncertainty of the stellar age (Delrez et al. 2016). Here we investigated to what extent the F_XUV will influence the Hα transmission spectrum. We adopted 0.5, 0.75, 1.5, 2.0, 3.0, and 4.0 times the value of fiducial F_XUV (F₀), and the mass-loss rates are 8.20 × 10¹¹ g s⁻¹, 1.07 × 10¹² g s⁻¹, 1.58 × 10¹² g s⁻¹, 1.84 × 10¹² g s⁻¹, 2.29 × 10¹² g s⁻¹ and 2.68 × 10¹² g s⁻¹ accordingly. The energy-limited theory (Lammer et al. 2003; Erkaev et al. 2007) proposes that the mass-loss rates of the atmosphere are proportional to the F_XUV. However, our results found that they do not increase linearly with the increase of the F_XUV because the escape of Lyα photons takes away a portion of the heat from the atmosphere. The Hα absorption increases moderately with F_XUV (see Figure 4(a), the transmission spectra are calculated within 7.6 R_p). For instance, an increase of a factor of 8 in F_XUV only increases the absorption at line center by ∼0.5%. Two reasons can explain this phenomenon. For higher F_XUV, the corresponding temperature of the atmosphere becomes higher at the bottom of the atmosphere but will drop dramatically with the increase of atmospheric altitude. The high temperature occurs in the relatively high pressure and is close to the bottom of the atmosphere, so that the Lyα photons will spend a longer time to be scattered out of the atmosphere. Therefore, the ${\bar{J}}_{\mathrm{Ly}\alpha }$ will be more intense and then will excite more hydrogen atoms into the n = 2 state. However, the atoms can be ionized easily owing to higher F_XUV. The combined effect is to increase the Hα absorption slightly. It is clear from Figure 4(a) that a higher F_XUV is needed to fit the 4-UT data, while a relatively lower F_XUV is preferable for the data of 1-UT. This can be verified by Figure 4(b), which shows the χ² as a function of F_XUV. The χ² is calculated in passband 1.5 and 3.6 Å for 1-UT and 4-UT, respectively. The minimum χ² for 1-UT and 4-UT occurs at F_XUV = F₀ ∗ 0.5 and F_XUV = F₀ ∗ 1.5, respectively. B21 reported that log${R}_{\mathrm{HK}}^{{\prime} }$ = −4.87 ± 0.01 and log${R}_{\mathrm{HK}}^{{\prime} }$ = −4.81 ± 0.01 for the 1-UT and 4-UT transits, respectively, and proposed that the star was more active at 4-UT transit. Therefore, the different F_XUV levels for 1-UT and 4-UT may reflect the different stellar activities of WASP-121 in the two observations. A higher F_XUV of 4-UT could probably attributes to its higher activity level compared to that of 1-UT.

Zoom In Zoom Out Reset image size

We also simulated the Hα transmission spectra for the cases of F_XUV = F₀ ∗ 0.5 and F_XUV = F₀ ∗ 1.5 as a function of different altitudes. Figure 5(a) shows the case of F_XUV = F₀ ∗ 0.5. The absorption caused by the atmosphere below the sonic point (1.7 R_p) is shallower compared with 1-UT of B20S. Figure 5(b) shows the equivalent width as a function of different altitudes. It shows when WASP-121b receives 0.5 times the fiducial XUV irradiation of the host star, the EW produced by the atmosphere lower than 2.5 R_p cannot reach the lower limit of the observation. Above 2.5 R_p, the more the atmosphere expands, the better the fit to the observation of 1-UT. Figures 5(c) and (d) are the same to Figures 5(a) and (b), respectively, but for the case of F_XUV = F₀ ∗ 1.5 and in comparison with 4-UT of B20S. The EW of the model cannot reach the lower limit of the EW of 4-UT, because our model cannot reproduce the relatively high absorption at the Hα line wings.

Zoom In Zoom Out Reset image size

XUV SEDs can also influence the photoionization process in the atmosphere (Guo & Ben-Jaffel 2016). According to Owen & Jackson (2012), X-ray can solely drive hydrodynamic escape of planetary atmosphere. Here we study the effect of different SEDs on the transmission spectra by introducing a modified spectral index β_m, defined as F(1–100 Å)/F(1–912 Å), where F(1–100 Å) is the integrated flux in the band of X-ray and F(1–912 Å) is the integrated flux of the whole XUV band. For the fiducial model, the value is 0.1475. We tested the cases of 0.03 (almost no X-ray but all EUV), 0.103, 0.221, and 0.5 (half X-ray and half EUV, which could be not real according to the evolution of XUV radiation of late-type stars, and it is for model experiment), while the XUV integrated flux F_XUV = F₀. Figure 4(c) shows the transmission spectra (calculated within 7.6 R_p). One can see that a larger X-ray proportion will lead to deeper Hα absorption. The main reason is that more hydrogen atoms will be retained instead of being ionized, due to the lower photoionization cross section (inversely proportional to the cube of the frequency) in X-ray band in comparison with that in EUV. To compare with the observations, Figure 4(d) shows the χ² as a function of β_m. For 1-UT, the minimum χ² appears at β_m = 0.003. For 4-UT, the χ² decreases with the increase of β_m. A higher β_m is required to fit the absorption at Hα line wings for 4-UT. However, the variation of χ² is less than 1 for β_m from 0.103 to 0.5. In addition, for the case of β_m = 0.5 the absorption at Hα line center exceeds the upper limit (+1σ) of the observation. Thus, a high β_m cannot be applicable for 4-UT although the value of the χ² is smaller. Thus, we suggest that β_m should be confined to a relatively low level in order to fit the observation at different times, which is consistent with the evolution of XUV radiation of late-type stars (Sanz-Forcada et al. 2011).