The PEPSI Exoplanet Transit Survey (PETS). II. A Deep Search for Thermal Inversion Agents in KELT-20 b/MASCARA-2 b with Emission and Transmission Spectroscopy^*

After receiving the reduced and normalized PEPSI spectra, more steps are necessary in order to detect the atmospheric signal. Our overall procedure is inspired by that of Nugroho et al. (2017), although with significant differences. We first run the molecfit package (Kausch et al. 2015; Smette et al. 2015) on the red arm of each PEPSI spectrum to model out the telluric lines, as is described in more detail in Section 3.3; no such treatment is necessary for the largely telluric-free blue arm. We median-stack the telluric-corrected spectra and subtract this median spectrum from each of the time-series spectra in order to remove the stellar lines and the time-invariant component of the telluric lines. We then run the SYSREM algorithm (Tamuz et al. 2005) on these residuals, as described in Section 3.3. Next, we cross-correlate a model spectrum generated as described in Section 3.2 with each spectrum. The cross-correlation algorithm is adapted from that in the BANZAI-NRES package (McCully et al. 2022). These procedures, and the tests that we conducted to show that this is the optimal method, are detailed in Section 3.3. We show the spectra and the results of molecfit in Figure 2.

Zoom In Zoom Out Reset image size

Finally, for a grid of possible radial velocity (RV) semi-amplitude K_P values for the planet, we shift each of the time-series cross-correlation functions (CCFs) to remove the expected Doppler shift from that K_P at the observation time; that is, we shift the spectra into the planetary rest frame. We assumed the ephemeris from Lund et al. (2017). We then stack the resulting shifted CCFs. For the transmission data we only combine the CCFs derived from the in-transit spectra. In this step we combine the CCFs derived from the blue and red arm data, and data taken on different nights. The exceptions are for FeH, where there are no lines present in the blue arm so we use only the red arm, and CaH, where we use only the blue arm due to the lack of red lines. We combine the CCFs using a weighted sum, where the weights are the 95th percentile S/N from that arm multiplied by the total equivalent width in the model emission or transmission lines within that bandpass.

We search for a CCF peak near the expected values of K_P and v_sys in this space; we calculated K _P given the system parameters from Lund et al. (2017), and v_sys was measured by the PEPSI pipeline. These parameters are listed in Table 2. We measure the S/N of the peak (if any) present at the location expected from the known parameters of the planet. We estimate the significance of any peaks using the standard deviation of the points in the shifted and combined CCFs with ∣v∣ > 100 km s⁻¹ from the expected location of the planetary signal.

We generate model spectra using the petitRADTRANS package (Molliere et al. 2019). We assumed the stellar and planetary parameters from Lund et al. (2017), which are listed in Table 2. We assumed a P-T profile of the form given by Equation (29) of Guillot (2010) as implemented in petitRADTRANS. We tested this profile with a simple grid search for the parameters κ_IR and γ, along with the P-T profiles retrieved for KELT-20 b by Yan et al. (2022b), Fu et al. (2022), Borsa et al. (2022), and Kasper et al. (2022). We computed a model Fe i spectrum for each profile assuming a volume mixing ratio (VMR) of 5.4 × 10⁻⁵, cross-correlated it with the data as described below, and logged the S/N of the resulting detection. Our best model is a Guillot profile with γ = 30 and κ_IR = 0.04, which resulted in a 16.94σ detection. We adopt this P-T profile for the remainder of this work. None of the published retrieved profiles outperformed our model, the best being that from Borsa et al. (2022) with 16.43σ. We show these profiles in Figure 3, and list the results of the P-T profile comparison in Table 3. We note that we used only an approximation to the P-T profile retrieved by Fu et al. (2022), as they used a more sophisticated profile tied to an atmospheric model. Instead, we approximated their profile using a simple model which is isothermal above and below the inversion, with a transition between the two regimes which is linear in log space, similar to those of Borsa et al. (2022) and Yan et al. (2022b). We also note that we only tested the best-fit P-T profile from each of the cited works and did not attempt to account for the quoted uncertainties on the profile parameters. For simplicity we use the same P-T profile for both the emission and transmission spectra; although in reality the P-T profile is likely to be different on the dayside versus the terminator, the dependence of the transmission spectrum on the details of the P-T profile is weak.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. System and Atmospheric Model Parameters

Parameter	Value	Source
RP (R⊕)	19.51	Lund et al. (2017)
MP (M⊕)	1072	Lund et al. (2017)
Teq (K)	2262	Lund et al. (2017)
P (d)	3.4741085	Lund et al. (2017)
T0 (BJD_TDB)	2457485.74965	Lund et al. (2017)
R⋆ (R⊙)	1.565	Lund et al. (2017)
M⋆ (M⊙)	1.76	Lund et al. (2017)
Teff (K)	8720	Lund et al. (2017)
[Fe/H]	−0.29	Lund et al. (2017)
K P (km s−1)	169 ± 6	Lund et al. (2017)
v sys (km s−1)	−26.0	⋯
$v\sin {i}_{P}$ (km s−1)	2.5	Lund et al. (2017)
κIR	0.04	⋯
γ	30	⋯
P0 (bar)	1	⋯
${X}_{{{\rm{H}}}_{2}}$	0.7496	⋯
XHe	0.2504	⋯
VMR (H−)	1 × 10−9	⋯

Download table as:  ASCIITypeset image

We generate spectra over the wavelength range 3850–7500 Å, covering the full PEPSI CD III+V bandpass, with an even spacing of 0.01 Å. We set up a model atmosphere with 130 grid points spanning pressures from 10² to 10⁻¹⁰ bar.

We include continuum opacity due to H⁻ bound-free absorption, H₂-H₂ and H₂-He collisions, and Rayleigh scattering due to H₂ and He. We assume a cloud-free atmosphere. We conservatively assumed an H⁻ VMR of 1 × 10⁻⁹, similar to that found by Arcangeli et al. (2018) for WASP-18 b, a planet several hundred Kelvin hotter than KELT-20 b. This is the major source of continuum opacity in the optical in our models, and is degenerate with the VMR ratio limits that we set using our transmission spectra (Section 3.4; Section 4.2). Including this level of H⁻ absorption worsens the transmission VMR limits by approximately two orders of magnitude over no H⁻ absorption; however, this has no significant effect on the emission results. We also note that we assume a reference pressure level P₀ of 1 bar; changing P₀ to 1 mbar has a similar effect upon the VMR limits to including the noted H⁻ absorption. We conservatively assume the high H⁻ abundance in order to account for the unknown P₀ level and H⁻ abundance. The single line absorber in question is the only line opacity source we include in our models.

After computation of the model spectra, we convolve the models with a rotational broadening kernel and an instrumental broadening kernel. For the emission spectra we assume an analytic form for the rotational broadening kernel from Equation (18.14) of Gray (2005). For the transmission spectra we instead assume the form from Equation (15) of Gandhi et al. (2022), which results in a double-peaked profile from the limbs of the planet. We further assume that KELT-20 b is tidally locked; given the planetary parameters from Lund et al. (2017), this implies an equatorial velocity of 2.5 km s⁻¹. For the instrumental broadening, we assume a Gaussian line spread function with a width corresponding to the PEPSI resolving power of R = 130,000 (Strassmeier et al.2018).

The treatment of the emission versus transmission spectra is slightly different. For the emission spectra, after calculating the planetary flux using petitRADTRANS, we generate a blackbody spectrum corresponding to the stellar T_eff and use this to estimate the planet-to-star flux ratio F_P /F_⋆. It is this quantity that we use as as our model spectrum for the cross-correlation.

For the transmission spectra, petitRADTRANS calculates the planetary radius as a function of wavelength R_P (λ). As our spectra have been normalized by the data reduction pipeline and are not flux-calibrated, we do not have any information in our data on the planetary radius, only the differential change in radius as a function of wavelength. We therefore estimate the transmission spectrum as:

$$\begin{eqnarray}&&F(\lambda )=1-\displaystyle \frac{{R}_{P}{\left(\lambda \right)}^{2}-{R}_{P0}^{2}}{{R}_{\star }^{2}}\end{eqnarray} \tag{ 1 }$$

where R_P0 is the broadband planetary radius, which we adopt following Lund et al. (2017).

We list the species considered and the sources of the opacity data used in Table 4, and show example spectra in Figure 4. We chose to search for the "classical" inversion agents TiO and VO (Fortney et al. 2008), and the metal hydrides FeH and CaH. petitRADTRANS does not currently have high-resolution opacity tables for several other species of interest, such as AlO, CaO, NaH, or MgH (Gandhi & Madhusudhan 2019). In the interests of simplicity, for TiO we only use the single isotope ⁴⁸TiO, as it makes up the majority of Ti atoms, and ignore the other stable isotopes. Likewise, for FeH the opacity table only considers the most common isotopologue, but unlike for TiO there are no petitRADTRANS opacity tables available for the other isotopologues.

In order to detect the minute planetary signal, it is necessary to correct for systematic sources of error, most notably the presence of telluric absorption lines in the spectra, but also systematics due to detector imperfections, etc. Most recent studies of exoplanetary atmospheres at high resolution use one of two methods to correct systematics: SYSREM (e.g., Nugroho et al. 2020; Cont et al. 2021a; Merritt et al. 2021) or molecfit (e.g., Cabot et al. 2021; Stangret et al. 2021; Kesseli et al. 2022). The SYSREM algorithm (Tamuz et al. 2005) is a principle component analysis (PCA)-based algorithm that removes common linear or time-varying systematics. The molecfit package (Kausch et al. 2015; Smette et al. 2015) fits and removes the telluric lines in the spectrum. After implementing both methods, we performed tests to evaluate which method works best for our PEPSI data.

In order to implement SYSREM, we modified the source code from the PySysRem package ¹⁵ to handle the PEPSI data. We added code to propagate uncertainties formally through the SYSREM algorithm. For all runs, we used the airmass of the exposure as the initial guess for the starting values of the first systematic removed. For the red arm, we split the spectra into regions with telluric contamination versus no telluric contamination, with a threshold of >1% absorption over multiple nearby lines, varying the number of systematics removed for the former regions and removing one systematic from the latter. Our tests indicated that the recovered S/N did not change significantly if we removed two versus one systematic(s) from the telluric-free regions.

Although molecfit was designed to work on data from ESO facilities, it has also been used for data from other spectrographs (e.g., Cabot et al. 2021). We only used molecfit to correct the PEPSI red arm data, as there are no significant telluric lines in the blue arm. We fit three narrow regions spanning the red arm (6290–6320 Å, 6470–6520 Å, and 7340–7410 Å) containing moderately strong telluric lines, and included H₂O and O₂ in our models. While the telluric model removes weak lines adequately, it does not fit the strong lines to any good precision as they are saturated, and we consequently mask out regions of the spectra containing poorly corrected lines in our further analysis. The telluric correction is sufficiently good for weak lines that any remaining residuals can be removed well by SYSREM (see below).

For all our tests, we inject a model emission spectrum of TiO, VO, FeH, or CaH into the PEPSI data. We constructed the model spectra as described in Section 3.2, using VMRs chosen to give a moderate but not overwhelming detection significance. We passed these spectra through SYSREM and then attempted to recover the signal using the cross-correlation methodology described in Section 3.1. We evaluated different numbers of SYSREM systematics removed in order to find the optimal number. For the red arm, we evaluated using SYSREM alone, molecfit alone, or molecfit followed by SYSREM.

We show the results of these tests in Figure 5. In all cases for the red arm, using molecfit alone outperforms using SYSREM alone. In many cases, the combination of molecfit and SYSREM outperforms using either alone. For most of the data sets, the optimal number of systematics removed is small (1–4). In only two of the eighteen data sets considered (arms × nights) is the optimal number greater than five.

Zoom In Zoom Out Reset image size

We summarize the optimal number of systematics removed that we use for the remainder of the paper in Table 5.

For the blue arm, the optimal number of SYSREM iterations is generally ≤5, except for one TIO data set which attained a maximum S/N after 14 iterations.

We did not inject the FeH signal in the blue arm because it does not have any significant lines in that region of the spectrum, and similarly we did not inject the CaH signal into the red arm.

Overall, the optimal strategy that we have used as is follows. For TiO and VO, we use both molecfit and SYSREM. For FeH, we use only molecfit, while for CaH, we use only SYSREM, since its lines are only present in the blue arm. For each species we adopt a number of SYSREM iterations approximating the average of the best number for the three data sets, which we list in Table 5.

In cases where we obtain a non-detection, it is useful to attempt to quantify these non-detections. Do the data constrain the non-detection to an interesting level, or are the data simply not high enough quality to place useful limits? We attempt to address this issue through injecting a model spectrum into our data before treating the data with SYSREM, and then recovering it with our cross-correlation process described in Section 3.1. For each species, we inject models with VMRs spaced 0.25 dex apart, i.e., VMR = [1.0, 1.8, 3.2, 5.6] × 10^x, and find the smallest VMR for which we can recover the signal at >4σ (see Kesseli et al. 2020). We note that these results are contingent upon the choice of P-T profile, as the strength of the inversion affects the strength of the emission lines, and thus the inferred VMR limit. However, Yan et al. (2022b), Fu et al. (2022), Borsa et al. (2022), and Kasper et al. (2022) all independently retrieved broadly similar P-T profiles using independent methodology and data, giving us some confidence in the broad accuracy of our P-T profile. Additionally, the presence of additional or stronger continuum opacity sources than we have assumed could affect the limits. These results should therefore only be taken as approximate under the assumption of the P-T profile and other model atmospheric parameters.

One other potential issue is with the inaccuracy of line lists, which can cause systematic reductions in signals strong enough to result in a non-detection when a detection would otherwise be expected (e.g., de Regt et al. 2022). In order to attempt to account for this issue, for species with multiple opacity tables derived from different line lists available (i.e., TiO and VO), we generate separate model spectra using the two line lists, inject one into the data, and attempt to recover it using the other line list. We have found that this worsens the detectable VMR by approximately one order of magnitude for TiO. For VO, on the other hand, this worsens the limits by two orders of magnitude for the emission spectra and the transmission spectra are not detectable at any VMR when using a different template due to their lower S/N. For completeness we quote the VO limits using both the same and differing line lists.

In order to perform a final validation of our methodology and data, we attempt to reproduce the detections of emission from Fe i from Yan et al. (2022b), Borsa et al. (2022), and Kasper et al. (2022). We estimated the VMRs for this model using a FastChem (Stock et al. 2018) equilibrium chemical model as implemented in PyFastChem. We assumed a solar abundance mixture, and the same P-T profile and pressure grid that we used to construct our model spectra. We assumed the maximum VMR present in the FastChem models, 5.4 × 10⁻⁵ for Fe i. We show the model spectra in the left panels of Figure 4, and the CCFs shifted and stacked according to the expected planetary orbit in Figure 6.

Zoom In Zoom Out Reset image size

We were able to recover a signal from Fe i near the expected (v, K_P ) values at a significance of 16.9σ. Yan et al. (2022b) obtained a 7.7σ detection, and Borsa et al. (2022) obtained a 7.1σ detection. It is unsurprising that we obtain a stronger detection, due to the much larger aperture of LBT (2 × 8 m) versus the smaller telescopes used by those works (3.5 m at Calar Alto Observatory, and 3.58 m at the Telescopio Nazionale Galileo), although our bandpass is smaller than that of these instruments, partially offsetting this advantage.

We searched for TiO, VO, FeH, and CaH emission and transmission in our spectra of KELT-20 b, and did not make any detections with a signficance of >4σ. Although the TiO transmission data do contain a bump near the correct (v, K_P ) peaking at 2.94σ (Figure 8), due to the presence of other similarly strong peaks elsewhere in the (v, K_P ) space we do not consider this to be a robust detection. We show the data and minimum recovered peaks in the Appendix.

Instead, we set limits upon the VMRs of these species as described in Section 3.4. These limits range from 3.2 × 10⁻⁷ for FeH to 1 × 10⁻⁹ for TiO for the emission data, and from 5.6 × 10⁻⁷ for FeH to 3.2 × 10⁻⁹ for TiO for the transmission data. We list the full results in Table 6, and show the results graphically in Figure 7.

Zoom In Zoom Out Reset image size

We derive the pressure levels for the VMR limits shown in Figure 7 from the contribution functions calculated by petitRADTRANS; we calculate the pressure of peak contribution and show as the error bars the range of pressures over which 68% of the contribution takes place. Counterintuitively, in our models the transmission spectra probe deeper into the atmosphere than the emission spectra. This is because, in our models, the lack of other optical opacity sources besides the absorber in question allows the transmission rays to penetrate deep into the atmosphere, below the temperature inversion. The emission lines, on the other hand, by definition are formed only in the inversion. Although the nominal transmission and emission models are consistent with each other, they are unlikely to be completely realistic. We also calculated spectra with an artificially increased VMR of H⁻ to 1 × 10⁻⁹, similar to the maximum VMR found by Arcangeli et al. (2018) for WASP-18 b, a planet ∼600 K hotter than KELT-20 b. Including this artificially increased H⁻ absorption does not significantly affect the emission results, but for the transmission results the limits worsen by 1–2 orders of magnitude, and the contribution moves by ∼1 dex lower in pressure.

In order to interpret the results, we used a FastChem (Stock et al. 2018) equilibrium chemical model as implemented in PyFastChem. We constructed two models, one assuming a solar abundance mixture and the other 10 × solar metallicity, with the same P-T profile and pressure grid that we used to construct our model spectra. We assumed that quenching occurs at a pressure of 1 bar, such that the VMRs of our species of interest are constant at P < 1 bar (see Zahnle & Marley 2014; Marley & Robinson 2015; Fortney et al. 2020). FastChem does not include photoionization or photodissociation, which may be important processes in the upper atmosphere of KELT-20 b. Photodissociation should serve to depress the concentrations of molecules including inversion agents, while photoionization should increase the concentrations of ionized species and H⁻. These effects are counteracted by the replenishing process of vertical mixing when we assume a quench pressure at 1 bar. Although these assumptions regarding quenching and the model's disregard of irradiation effects are likely to be only approximately correct, this nonetheless serves as a useful tool to help interpret our limits.

We show these predicted VMRs from FastChem as vertical lines in Figure 7. Our VMR limits for TiO and CaH lie several orders of magnitude below the expectations for those species at both metallicities. Our limits for FeH lie above the expected concentrations for all metallicities and are thus not constraining. Most of the strong FeH lines lie redward of the PEPSI CD V bandpass; new observations with PEPSI CD VI (which covers 7419–9067 Å; Strassmeier et al. 2015) would be required to better constrain FeH with PEPSI. For VO, our limits are constraining only if the ExoMol line list (McKemmish et al. 2016) is accurate; otherwise, we cannot claim a constraining limit (see de Regt et al. 2022). We again note that our VMR limits are only approximate, and in detail subject to our assumptions regarding the P-T profile, continuum opacity from H⁻, and reference pressure level. Nonetheless, the large mismatch between the expected VMRs and our limits for TiO and CaH suggest that these species are depleted in the atmosphere of KELT-20 b and cannot be responsible for the temperature inversion, while VO may also be depleted if the ExoMol line list is sufficiently accurate.

This could be due to cold trapping of these species by rainout of condensed species on the cold planetary nightside (e.g., Spiegel et al. 2009; Parmentier et al. 2013), which was previously suggested to occur for KELT-20 b by Nugroho et al. (2020). If these species can condense out into particles that are large enough to settle deep into the atmosphere despite any convection on the nightside, then these constituents will be depleted from the upper layers of the atmosphere that we can probe with transmission or emission spectroscopy. Neither atomic Ti nor V have been detected in the transmission spectrum of KELT-20 b (Nugroho et al. 2020). Since neither the expected atomic nor molecular Ti- and V-bearing species are seen in the atmosphere of KELT-20 b, it is plausible that rainout has removed these molecules from the upper atmosphere. Rainout, however, cannot explain the lack of detections of CaH, as Nugroho et al. (2020) and Casasayas-Barris et al. (2019) detected Ca in their transmission spectra, indicating that this species does not rain out. Alternately, molecular species could be photodissociated by the intense ultraviolet radiation from the early-type host star KELT-20. Tsai et al. (2021), however, found that for the even hotter UHJ WASP-33 b that photochemistry should dominate only at high altitudes (P < 10⁻⁴ bar).