First Detection of Hydroxyl Radical Emission from an Exoplanet Atmosphere: High-dispersion Characterization of WASP-33b Using Subaru/IRD

High-resolution spectroscopy is one of the most successful methods to characterize exoplanet atmospheres, especially those of hot Jupiters. The resolved planetary lines are disentangled from the telluric and stellar lines due to the planetary orbital motion allowing us to unambiguously detect atomic/molecular signatures in the atmosphere of exoplanets (e.g., Snellen et al. 2010). By comparing hundreds/thousands of unique absorption/emission lines to model templates through cross-correlation, we can highly boost the signal and constrain the chemical abundances, the planetary rotation, the projected equatorial wind, and even the temperature–pressure (T-P) profile of the atmosphere for emission spectroscopy data.

With the equilibrium temperature (T_eq) similar to M-dwarfs, similar prominent atomic/molecular opacity sources (Fe i, Fe ii, and Na i, as well as Ti ii, TiO, VO, AlO, FeH, CO, H₂O, and OH) are expected to be found in the atmosphere of hot Jupiters (T_eq < 2200 K) and ultra-hot Jupiters (T_eq > 2200 K). Many of these optical opacity sources have been detected (e.g., Nugroho et al. 2017; Hoeijmakers et al. 2019; von Essen et al. 2019; Yan et al. 2019, 2021; Pino et al. 2020). In the near-infrared, however, the most frequently detected species using high-resolution spectroscopy are CO and H₂O (Snellen et al. 2010; Birkby et al. 2013; De Kok et al. 2013; Lockwood et al. 2014; Wang et al. 2018; Cabot et al. 2019; Webb et al. 2020). The opacities in this wavelength region play an important role as a coolant in the atmosphere, and constraining their abundance would allow us to study the climate of a unique planetary population and estimate the C/O ratios thus inferring planetary formation history (Öberg et al. 2011). We, therefore, analyzed the dayside spectrum of WASP-33b, one of the hottest ultra-hot Jupiters (T_day > 3100 K, e.g., De Mooij et al. 2013) orbiting a fast-rotating δ-Scuti A5-type star (Collier Cameron et al. 2010), to search for molecular signatures in the near-infrared.

In this Letter, we present the first detection of high-resolution OH emission and the evidence of H₂O emission in the dayside spectra of WASP-33b. In Section 2, we describe the observations and data reduction. We then describe our modeling of the planetary emission spectrum in Section 3, and in Section 4, we detail our methodology in validating the accuracy of the OH line lists and searching the signal of OH and H₂O in the atmosphere of the planet. Finally, in Section 5, we present and discuss our findings and their implications for the planetary atmosphere.

We observed WASP-33 in the second half of the night of 2020 September 30 using the InfraRed Doppler instrument (IRD; R ≈ 70,000; λ ≈ 0.97–1.75 μm; Tamura et al. 2012; Kotani et al. 2018) on the Subaru 8.2 m telescope (PID: S20B-008, PI: S.K. Nugroho). We continuously observed the target with an exposure time of 300 s per frame without the laser frequency comb in natural guide star mode. The weather during the observations was not always stable; therefore, we were only able to obtain 33 exposures covering the orbital phase of WASP-33b from ≈0.597 to 0.700 (there were some gaps close to the middle of the observations due to clouds). We converted the Julian Date UTC to Barycentric Julian Date in Barycentric Dynamical Time (BJD_TDB) using the online calculator from Eastman et al. (2010) then calculated the orbital phase using the transit epoch taken from Johnson et al. (2015).

The data were reduced following Hirano et al. (2020) resulting in 70 spectral orders ranging from ≈9260 to 17419 Å with an average signal-to-noise ratio (S/N) of 140. We fitted the continuum of the spectrum with the highest average S/N using the continuum task in IRAF ²⁰ and divided it out from the data. Any possible blaze function variations were then corrected following the procedure in Nugroho et al. (2020b). Then, the sky emission lines, bad pixels, and regions were visually identified and masked. Additionally, we also masked any pixels that had a flux less than 10% of the continuum. In total, we masked 13.9% of the total number of pixels of the data. Finally, the spectra of each spectral order were aligned into two-dimensional arrays with wavelength along one axis and orbital phase along with the other. We estimated the uncertainty of each pixel by taking the outer product of the standard deviation of each wavelength and exposure bin, then normalized by the standard deviation of the whole array.

To check if there is any wavelength shift during the observation, the data were cross-correlated with the Doppler-shifted telluric templates produced using the Cerro Paranal Sky Model (Noll et al. 2012; Jones et al. 2013) over a velocity range of −50 to 50 km s⁻¹ in 0.01 km s⁻¹ steps. We found no significant shift (<0.05 km s⁻¹) compared to the precision that we need for this analysis; therefore, we did not attempt to correct for this.

Before searching for any exoplanetary atomic or molecular signal using cross-correlation, we removed the telluric and stellar lines using a detrending algorithm, SysRem (Tamuz et al. 2005), which has been successfully adopted for high-resolution Doppler spectroscopy (e.g., Birkby et al. 2013). SysRem fits the systematic trend in the wavelength bin direction, which might be due to variation in airmass, water vapor column level, and other factors. Following Gibson et al. (2020), we run SysRem directly in flux for each spectral order independently. For each iteration, we summed the best-fit SysRem model and divided out from the data and the uncertainty array to propagate the error. Finally, any outliers more than five times the standard deviation of the residual were masked. A step-by-step overview of these procedures is shown in Figure 1. As in the previous analyses using SysRem (Gibson et al. 2020; Merritt et al. 2020; Nugroho et al. 2020a, 2020b; Yan et al. 2020), instead of determining the optimal SysRem iteration of each order, we used the same number of iteration for all orders. The results are shown in Section 5.

Zoom In Zoom Out Reset image size

The planetary emission spectrum template was created by assuming 70 atmospheric layers evenly spaced in log pressure from 10² to 10⁻⁸ bar of 1D plane-parallel hydro-static atmosphere, and a planetary-mass and radius of 3.266 M_J and 1.679 R_J, respectively (Kovács et al. 2013). We adopted a thermally inverted T-P profile used in Nugroho et al. (2020a) which was calculated using the equation in Guillot (2010) assuming the visible mean opacity is twice the infrared mean opacity (0.01 cm² g⁻¹, e.g., dominated by H⁻ opacity), an internal temperature of 100 K and T_eq of 3100 K (assuming uniform dayside only reradiation).

We produced five planetary spectrum models with single molecular opacity for H₂O and OH, with three and two different line lists, respectively. The cross-sections of molecular species were computed using HELIOS-K (Grimm & Heng 2015) at a resolution of 0.01 cm⁻¹ from 9100 to 17800 Å assuming a Voigt line profile taking into account natural and thermal broadening only and a line wing cutoff of 100 cm⁻¹. For H₂O, we used the line list database of POKAZATEL (Polyansky et al. 2018), HITEMP 2010 (Rothman et al. 2010), and BT2 (Barber et al. 2006); for OH, we used the updated line-list database of HITEMP (the updated HITEMP 2020; Rothman et al. 2010; Brooke et al. 2016; Yousefi et al. 2018; Noll et al. 2020) and MoLLIST (Brooke et al. 2016; Yousefi et al. 2018; Bernath 2020). For continuum opacity, we included the bound–free and free–free absorption of H⁻ using the equation from John (1988), and collision-induced absorption (CIA) of H₂–H₂ (Abel et al. 2011) and H₂–He (Abel et al. 2012).

We used FastChem (Stock et al. 2018) to estimate the abundances of chemical species, and the mean molecular weight of each atmospheric layer assuming chemical equilibrium and solar C/O. We then produced the emission spectrum by solving the Schwarzchild equation following Nugroho et al. (2017, 2020a). We divided the resulting spectra by the flux of the star assuming a blackbody spectrum, R_⋆ of 1.509 R_⊙, and T_eff of 7400 K then convolved with a Gaussian kernel to the spectral resolution of IRD. ²¹ Finally, we subtracted the planetary continuum from each model, which was determined by the continuum opacity of CIA and H⁻, estimated using a minimum filter with a window of 55 Å. The final result is the line contrast relative to the stellar continuum profile (see Figure 2).

Zoom In Zoom Out Reset image size

4.1. Accuracy of the Position of the Lines in the OH Line Lists

4.2. Searching for OH and H₂O Signatures

Even after removing the telluric and stellar lines, the planetary signal is still expected to be buried under the noise. As there are many resolved lines of H₂O and OH in the IRD wavelength range, we combined them to boost the planetary signal by cross-correlating the residual of each SysRem iterations with the Doppler-shifted planetary spectrum templates from −500 to +500 km s⁻¹ in 1 km s⁻¹ steps following:

$$\begin{eqnarray}&&\mathrm{CCF}(v)=\displaystyle \sum _{i}\displaystyle \frac{{f}_{i}{m}_{i}(v)}{{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 1 }$$

where f_i is the mean-subtracted data, m_i is the mean-subtracted spectrum model Doppler-shifted to a radial velocity of v, and ${\sigma }_{i}^{2}$ is the variance at ith wavelength bin. We performed this for each spectral order and summed them excluding the spectral order that has significant telluric removal residuals (spectral order with central wavelength of 13529.12, 14317.24, and 14457.64 Å).

We calculated an orbital velocity−systemic velocity (K_p–v_sys) map by shifting the cross-correlation functions (CCFs) to the planetary rest-frame over a range of K_p, from 0 to +300 km s⁻¹, and v_sys, from −125 to +125 km s⁻¹, both in 0.2 km s⁻¹ steps using linear interpolation and summed over time. The radial velocity of the planet at a given orbital phase (ϕ), RV_p(ϕ), assuming the planet has a circular orbit is

$$\begin{eqnarray}&&{\mathrm{RV}}_{{\rm{p}}}(\phi )={K}_{{\rm{p}}}\sin (2\pi \phi )+{v}_{\mathrm{sys}}+{v}_{\mathrm{bary}},\end{eqnarray} \tag{ 2 }$$

where v_bary is the barycentric correction, ²² and ϕ is the orbital phase of the planet.

From Collier Cameron et al. (2010), Nugroho et al. (2017, 2020a), and Yan et al. (2019), the planet signal is expected at K_p of ≈230 km s⁻¹ and v_sys of ≈−3 km s⁻¹. We computed the S/N by dividing the K_p–v_sys map by its standard deviation calculated by avoiding the area ±25 km s⁻¹ from the expected planet signal.

Next, we converted the cross-correlation map to a likelihood map (${ \mathcal L }$) using the β-optimized likelihood function following (Gibson et al. 2020, see also Brogi & Line 2019):

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{N}{2}\mathrm{ln}\left[\displaystyle \frac{1}{N}\left(\displaystyle \sum \displaystyle \frac{{f}_{i}^{2}}{{\sigma }_{i}^{2}}+{\alpha }^{2}\displaystyle \sum \displaystyle \frac{{m}_{i}^{2}}{{\sigma }_{i}^{2}}-2\alpha \mathrm{CCF}\right)\right],\end{eqnarray} \tag{ 3 }$$

where α is the scale factor of the model and N is the total number of pixels. A three-dimensional likelihood or posterior data cube (assuming uniform priors) was then produced with a range of α from 0.01 to 1.50 in 0.01 steps by subtracting the global maximum value from the cube and calculating the exponential, this normalized the likelihood to 1. We then marginalized it by summing the maps over parameters to get the best-fit parameters and uncertainties. We estimated the significance of detection by dividing the median value of the conditional distribution of α at the best-fit value of K_p and v_sys by its uncertainty.

5.1. OH Emission in the Dayside of WASP-33b

We detected the OH emission signature at an S/N of 5.4 and a significance of 5.5σ at K_p of ${230.9}_{-7.4}^{+6.9}$ km s⁻¹ and v_sys of −0.3${}_{-5.6}^{+5.3}$ km s⁻¹ (see Figures 3(b) and 4) consistent with previous results although with larger uncertainties due to narrower orbital phase coverage (e.g., Nugroho et al. 2017, 2020a; Yan et al. 2019). From Figure 3(a), the planet signal appears as a bright stripe shown by the white arrows. The strength of the signal varied with time, which might be due to the unstable weather during the observation that potentially affects the telluric removal using SysRem. The α is constrained to 0.47 ± 0.09, which means that we have overestimated the strength of the signal. This could be due to the effect of SysRem, which might have eroded and altered the observed exoplanet signal. Furthermore, the inhomogeneity of the dayside of WASP-33b or the overestimation of the T-P profile and/or the OH abundance in the modeling (photo-dissociation and other possible overlapping opacity sources) could also be the cause. We leave a more detailed analysis to future works.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For both line lists, the detected signals become prominent after two SysRem iterations and at their highest S/N after three iterations before getting weaker with more iterations. The cross-correlation maps, the S/N and/or detection significance, and the constraint on alpha for both templates are the same; therefore, we show the result for OH HITEMP 2020 only. The only difference of the result using the two line lists is in the velocity constraint, which differs by ≈0.1 km s⁻¹. This is expected as the two line lists are based on the same data, although the OH HITEMP 2020 line list was just recently updated based on the observed high-resolution OH telluric lines (Noll et al. 2020).

Compared to other ultra-hot Jupiters, the atmosphere of WASP-33b's atmosphere is more difficult to characterize with low-resolution spectroscopy/photometry due to the δ-Scuti pulsations of its host star. These pulsations can also affect high-spectral resolution searches for atomic species that occur both in the stellar photosphere and in the planet's atmosphere, such as Fe i (Nugroho et al. 2020a; Herman et al. 2021) where a region of ∼$\pm {v}_{\mathrm{rot}\star }\,\sin \,i$ has to be excluded from analysis, as the stellar pulsations overlap with the signal from the planet. Since OH is not present in the stellar atmosphere of an A-star, this does not pose a problem for the results presented in this paper. In addition, the lack of signal at K_p of 0 km s⁻¹ or RV of 0 km s⁻¹ indicates that there is no contamination from telluric OH emission. Furthermore, as the trail of the signal appeared only at the expected planetary velocity (see Figures 3(a) and (b)), we are confident that the detected OH emission signature is originating from the exoplanet.

5.2. Marginal Detection of Weak H₂O Emission?

While the signature of OH has been detected in the atmosphere of Earth, the Saturn magnetosphere, Venus, and Mars (Meinel 1950; Shemansky et al. 1993; Piccioni et al. 2008; Todd Clancy et al. 2013), this is the first time that its signature has been detected in the atmosphere of an exoplanet. Along with O, OH is one of the most important radical species that drive atmospheric chemistry. For a hot Jupiter like HD 209458b, OH is mainly produced from the photolysis of H₂O by the stellar UV (Liang et al. 2003). However, for a much hotter planet like WASP-33b, the thermochemical reaction is expected to be the dominant source of OH as the atmosphere is closer to thermochemical equilibrium (Visscher et al. 2006). Our result, which only marginally detected weak emission of H₂O, indicates that most of the H₂O in the upper atmosphere is thermally dissociated consistent with the theoretical predictions (Parmentier et al. 2018). Thus, OH is expected to be one of the most dominant O-bearing molecules along with CO and should be considered when analyzing the emission spectrum of ultra-hot Jupiters and searched for in other planetary atmosphere.

Through an injection test, if H₂O were present in the dayside of WASP-33b at the abundance and temperature retrieved by Haynes et al. (2015), we would detect it at high significance. As low-resolution spectroscopy probes a relatively deeper atmospheric layer than high-resolution spectroscopy, the retrieved parameters might not provide a reliable measurement for the upper atmosphere. Moreover, when there are overlapping unresolved features from multiple species, the retrieved parameters would be incorrect if the model does not consider all of the possible chemical species. Thus, combining low-resolution spectroscopy and high-resolution spectroscopy would be required to get a more accurate and precise characterization of the exoplanet atmosphere (Brogi & Line 2019; Gandhi et al. 2019). Finally, this work demonstrates the capability of IRD in characterizing the atmosphere of an exoplanet and its potential to complement the space-borne facilities (e.g., Hubble Space Telescope, James Webb Space Telescope).

We are extremely grateful to the anonymous referee for constructive and insightful comments that greatly improved the quality of this Letter. This work is based on data collected at the Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. Our data reductions benefited from (PyRAF and) PyFITS, which are the products of the Space Telescope Science Institute, which is operated by AURA for NASA. The M-dwarf spectra are based on data collected at the Subaru Telescope and obtained from the SMOKA, which is operated by the Astronomy Data Center, National Astronomical Observatory of Japan. H.K. is supported by a Grant-in-Aid from JSPS (Japan Society for the Promotion of Science), Nos. JP18H04577, JP18H01247, and JP20H00170. This work was also supported by the JSPS Core-to-Core Program "Planet²" and SATELLITE Research from Astrobiology Center (AB022006). N.P.G. gratefully acknowledges support from the Science Foundation Ireland and the Royal Society in the form of a University Research Fellowship. T.H. acknowledges support from JSPS KAKENHI grant No. 19K14783. Y.K. is supported by Special Postdoctoral Researcher Program at RIKEN. C.A.W. would like to acknowledge support from UK Science Technology and Facility Council grant ST/P000312/1. M.T. would like to acknowledge support from MEXT/JSPS KAKENHI grant Nos. 18H05442, 15H02063, and 22000005. M.I. acknowledges support from JSPS KAKENHI grant No. 19J11805. J.L.B. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 805445. M.B. acknowledges support from the UK Science and Technology Facilities Council (STFC) research grant ST/S000631/1. We are also grateful to the developers of the Numpy, Scipy, Matplotlib, Jupyter Notebook, and Astropy packages, which were used extensively in this work (Hunter 2007; Astropy Collaboration et al. 2013; Kluyver et al. 2016; Price-Whelan et al. 2018; Virtanen et al. 2020).

The result of the accuracy test that we performed for OH line lists can be seen in Figure 5. We also searched for H₂O signal with two other line lists, HITEMP 2010 and BT2 and the results are shown in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size