Atmospheric characterization of the ultra-hot Jupiter WASP-33b Detection of Ti and V emission lines and retrieval of a broadened line proﬁle

Ultra-hot Jupiters are highly irradiated gas giant exoplanets on close-in orbits around their host stars. The dayside atmospheres of these objects strongly emit thermal radiation due to their elevated temperatures, making them prime targets for characterization by emission spectroscopy. We analyzed high-resolution spectra from CARMENES, HARPS-N


Introduction
In-depth characterization of exoplanet atmospheres is an emerging field in astronomy. Most known exoplanets have no analogs in our Solar System and therefore challenge our understanding of their formation and evolution. In recent years, observations and theoretical work have increased our knowledge about the subclass of ultra-hot Jupiters (UHJs; Parmentier et al. 2018), which are gas giant planets with the highest equilibrium temperatures identified to date (T eq ≥ 2200 K). Located on close-in orbits, they are expected to be tidally locked to their host stars with highly irradiated daysides. The extreme irradiation regime in combination with a synchronous rotation causes a strong temperature contrast between the permanent day-and nightsides. Therefore, the atmospheric composition is significantly different between the two planetary hemispheres (e.g., Arcangeli et al. 2018;Bell & Cowan 2018;Komacek & Tan 2018;Tan & Komacek 2019;Molaverdikhani et al. 2020). UHJ daysides are likely dominated by atomic species with a high degree of ionization, while most of the molecules are predicted to dissociate (Kitzmann et al. 2018;Lothringer et al. 2018). In contrast, the planetary nightsides harbor a large variety of molecules and can be covered by clouds with complex molecular chemistry (Helling et al. 2019;Gao et al. 2020;Gao & Powell 2021).
Ultra-hot Jupiters are prime targets for atmospheric characterization by observing their thermal emission signal, given the elevated dayside temperatures of these objects. Various techniques have been used to study the emission signals from UHJ atmospheres, such as space-based photometric light curves (e.g., Zhang et al. 2018;von Essen et al. 2020;Deline et al. 2022), low-spectral-resolution observations (e.g., Arcangeli et al. 2018;Changeat & Edwards 2021), and ground-based high-resolution spectroscopy. In particular, emission spectroscopy has proven to be a powerful tool for identifying thermal inversions in the atmospheres of UHJs (e.g., Haynes et al. 2015;Nugroho et al. 2017;Kreidberg et al. 2018;Yan et al. 2020;Kasper et al. 2021). In these atmospheres, strong absorption of visible and ultraviolet radiation from the host star causes the temperature to increase with altitude. The presence of TiO and VO was initially proposed to cause this heating mechanism in upper planetary atmospheres (Hubeny et al. 2003;Fortney et al. 2008). However, theoretical work and more recent observations at high spectral resolution suggest that atomic metal species may also be fundamental for maintaining thermal inversion layers (e.g., Arcangeli et al. 2018;Lothringer et al. 2018). Neutral Fe is the species most commonly detected in thermal inversion layers, and is therefore thought to play an important role in the temperature structure of UHJ atmo-  Lehmann et al. (2015) with parameters from Kovács et al. (2013) , (b) Maciejewski et al. (2018) , (c) Nugroho et al. (2017) , (d) Kovács et al. (2013) , (e) Collier Cameron et al. (2010) , (f) Johnson et al. (2015) -other studies (e.g., Collier Cameron et al. 2010) report larger uncertainties; hence, we have also performed our calculations with a rot sin i * that deviates by ± 10 km s −1 , but could not find any significant differences in the results.
spheres (Pino et al. 2020;Nugroho et al. 2020a;Yan et al. 2020Yan et al. , 2021Kasper et al. 2021;Cont et al. 2021;Herman et al. 2022). Another significant result of UHJ high-spectral-resolution observations is the frequent nondetection of Ti, V, and their oxides, although their signature should be present in the planetary spectra when assuming equilibrium chemistry (e.g., Merritt et al. 2020;Hoeijmakers et al. 2020b;Tabernero et al. 2021;Yan et al. 2022b). This result suggests that depletion mechanisms may limit the presence of these species in the upper atmospheres of UHJs. Theoretical work has explored different mechanisms that may remove Ti-and V-bearing species from the atmospheres of hot giant exoplanets. For example, Spiegel et al. (2009) proposed the depletion of Ti and V in the upper part of planetary atmospheres via gravitational settling of TiO and VO. Cold trapping of Ti-and V-bearing molecules on the planetary nightsides has been suggested as another possible depletion scenario (Parmentier et al. 2013). However, these theoretical studies are mostly limited to the temperature range below that of UHJs. Further theoretical and observational work is therefore needed to better understand the circumstances under which Ti-and Vbearing species are depleted in UHJ atmospheres and how this impacts the planetary energy balance.
Several UHJs have been characterized using high-resolution Doppler spectroscopy. This technique uses the Doppler-shift of a planet's orbital motion relative to the telluric and stellar lines to identify its spectral signature (e.g., Snellen et al. 2010;Brogi et al. 2012). However, constraining the atmospheric parameters from the observed spectra remains a difficult task. For example, the elemental abundances and ratios or the atmospheric temperature structure cannot be accurately determined using highresolution Doppler spectroscopy alone. To overcome this dif- Notes. (a) Total number of CARMENES spectra was 105; due to their poor quality ten spectra were removed from the time series of the VIS channel and 17 spectra were removed from the time series of the NIR channel.
ficulty, atmospheric retrieval frameworks have been developed to fit high-resolution spectroscopy observations with parameterized model spectra. These retrievals result in statistical estimates of the physical parameters in the exoplanet atmosphere. Most of these retrievals rely on a likelihood function that compares the observed data to a forward model. A Bayesian estimator is then used to optimize the forward model parameters (e.g., Brogi et al. 2017;Brogi & Line 2019;Yan et al. 2020;Gibson et al. 2020). The use of alternative techniques for retrievals, such as supervised machine learning, has also been explored recently (Fisher et al. 2020).
In this study, we investigate the dayside emission spectrum of the UHJ WASP-33b. We report the first detection of atomic Ti and V via emission spectroscopy. Also, we present the physical parameters of the planetary atmosphere obtained from a Bayesian retrieval. WASP-33b orbits an A5-type star undergoing δ Scuti pulsations. Its bright host star (V ∼ 8 mag), high equilibrium temperature (T eq ∼ 2700 K), and short orbital period (P orb ∼ 1.22 d) make WASP-33b a prime target for emission spectroscopy. Different studies have shown that the planet possesses a temperature inversion in its dayside hemisphere (e.g., Haynes et al. 2015;Nugroho et al. 2017;Cont et al. 2021). To date, the hydrogen Balmer lines, Si i, Ca ii, Fe i, OH, CO, H 2 O, and TiO have been identified via transmission or emission spectroscopy in the planetary atmosphere (Nugroho et al. 2017(Nugroho et al. , 2020aYan et al. 2019Yan et al. , 2021Cont et al. 2021Cont et al. , 2022Herman et al. 2022;van Sluijs et al. 2022). The WASP-33 system parameters used in this work are summarized in Table 1.
This work is organized as follows. We describe the observations and our data reduction routine in Sects. 2 and 3. The technique to search for individual species in the planetary atmosphere and the resulting detections are presented in Sect. 4. Our retrieval method and the inferred atmospheric parameters are described in Sect. 5. The conclusions of our work are given in Sect. 6.

Observations
The thermal phase curve of WASP-33b was observed at high spectral resolution over a total of eight visits. We observed the planetary emission spectrum on 15 November 2017 with CARMENES at the 3.5 m Calar Alto telescope (Quirrenbach et al. 2014(Quirrenbach et al. , 2020. Both the visible (VIS) and near-infrared (NIR) wavelength channels of CARMENES were used. We consider the CARMENES channels independently as two different instruments. We also observed the spectrum of WASP-33b during the two nights of 15 October 2020 and 7 November 2020 with HARPS-N at the Telescopio Nazionale Galileo (Mayor et al. 2003;Cosentino et al. 2012). Another five observations were obtained between September 2013 and November 2014 by Herman et al. (2020) with ESPaDOnS at the Canada-France-Hawai'i telescope (Donati 2003). The combined set of nine observation datasets (two datasets from the CARMENES VIS and NIR observations, one dataset from each of the other seven observations) has already been used in previous studies to investigate the emission signature of TiO, Fe i, and Si i in the dayside atmosphere of WASP-33b (Herman et al. 2020(Herman et al. , 2022Cont et al. 2021Cont et al. , 2022. Details of the observations and the main characteristics of the used spectrographs are summarized in Table 2. Prior to the data reduction, we discarded a number of spectra due to poor data quality. All the discarded spectra belong to the CARMENES observation night on 15 November 2017. For seven spectra, the star was not aligned with the telescope fiber, and for three spectra the flux drastically dropped due to a passing cloud. Another seven spectra were removed due to the insufficient flux level in the NIR channel caused by the low elevation of the target towards the end of the observation night (airmass > 2). Consequently, we removed a total of 10 VIS channel spectra and 17 NIR channel spectra from the CARMENES observation. As no data quality issues were identified for the HARPS-N and ES-PaDOnS observations, we included all spectra from the two instruments in the further analysis.

Data reduction
3.1. Pre-processing the spectra The extraction of the one-dimensional spectra from the raw frames was performed with the data reduction pipelines of the respective instruments. We extracted the order-by-order spectra of CARMENES using the instrument pipeline caracal v2.20 (Zechmeister et al. 2014;Caballero et al. 2016) and obtained the order-merged spectra of HARPS-N by running the Data Reduction Software (Cosentino et al. 2014). The order-by-order ESPaDOnS spectra were extracted by the observatory using the Upena pipeline, which relies on the data reduction package Libre-ESpRIT (Donati et al. 1997). Except for the echelle orders 45-43 in the CARMENES NIR channel, we included the entire wavelength range of all instruments in our analysis. The three excluded CARMENES NIR orders coincide with the telluric absorption band of water at 1.4 µm and therefore suffer from an insufficient flux level.
Each of the nine observation datasets was reduced separately. To normalize all the spectra to the same continuum level, we applied a polynomial fit to the individual spectra and divided them by the fit function. We used a second-order polynomial for the order-by-order spectra of CARMENES and ESPaDOnS and a seventh-order polynomial for the order-merged spectra of HARPS-N. Outliers were removed by applying a 5σ clip to the time evolution of each pixel. Wavelength bins with the absorption level larger than 80 % of the spectral continuum were masked. We also masked the strong sky emission lines, which are present in the CARMENES NIR channel.

Removal of telluric and stellar features
We applied the detrending algorithm SYSREM to the spectral matrix of each observation to correct for the contribution of telluric and stellar lines. The algorithm was originally proposed for the removal of systematic effects from large sets of photometric light curves (Tamuz et al. 2005). SYSREM models the common systematic features of the light curves by iteratively fitting their trends as a function of time. Subsequently, the modeled systematics are subtracted from the data. When applied to the search for exoplanets, each wavelength bin of the spectral matrix is considered as an independent light curve. This procedure results in the socalled residual spectral matrix, which is the spectral matrix after removal of the systematics.
The input data for SYSREM are the spectral matrix and its uncertainties. For the CARMENES and ESPaDOns data, we used the propagated uncertainties from the instrument pipelines. As the HARPS-N pipeline does not compute uncertainties, we estimated them according to the procedure described by Yan et al. (2020). This method consists of running SYSREM in a first step with five iterations and uniform uncertainty values. The resulting residual matrix is dominated by noise. We then calculated the average noise row by row and column by column. Finally, the uncertainty values of each data point are calculated as the mean of the respective row and column average noise.
We followed the approach from Gibson et al. (2020) and ran SYSREM in flux space instead of magnitude space (Tamuz et al. 2005). This method allows the relative strength of the planetary spectral lines to be preserved during the correction for the telluric and stellar contamination. First, we used the algorithm in the standard way to calculate a model of the systematic features. This model corresponds to the sum of the models from each SYSREM iteration. We then divided the model of the systematic features out of the original spectral matrix. Also, the uncertainties were divided by the model. We ran the algorithm over ten consecutive iterations. This results in a residual spectral matrix for each iteration. For the CARMENES and ESPaDOnS observations, we created the order-merged residual spectral matrix by combining the orderby-order SYSREM reduced spectral matrices. This step was not necessary for the HARPS-N data, because the spectra were already provided by the instrument pipeline in an order-merged format. An overview of the data reduction including SYSREM is provided in Fig. 1.

Detection of the planetary emission lines
We used the cross-correlation technique (e.g., Snellen et al. 2010;Brogi et al. 2012;Sánchez-López et al. 2019;Prinoth et al. 2022) to extract the weak emission signature of WASP-33b from the noise-dominated residual spectra. This method combines numerous spectral lines and translates them into a single peak via calculating the cross-correlation function (CCF) between the residual spectra and model spectra as described below. We searched for the emission lines of the metal species Ti i, Ti ii, V i, and V ii. Also, we aimed to confirm the presence of the hydroxyl radical (OH), which was recently detected by Nugroho et al. (2021). Moreover, we reassessed the detections of Fe i and Si i in previous studies (Cont et al. 2021(Cont et al. , 2022 and investigated the presence of Fe ii and Si ii in the planetary atmosphere.

Model spectra
We modeled a planetary atmosphere with 61 layers equispaced in log pressure from 1 to 10 −6 bar. As the atmospheric temperature profile of WASP-33b is not yet known in detail, we adopted the T -p profile of WASP-189b, which was measured by Yan et al. (2020) via the Fe i emission signature (Fig. 2). We believe 10 11 10 10 10 9 10 8 10 7 10 6 10 5 10 4 VMR and T -p profile (right panel) used to generate the model spectra. We assumed equilibrium chemistry and solar elemental abundances for computing the VMRs. As the T -p pattern is not known in detail, we assumed the two-point profile that was retrieved for WASP-189b by Yan et al. (2020).
that using the T -p profile of WASP-189b is a reasonable approximation because the planet has similar properties to WASP-33b. This temperature profile has also been used successfully in previous work to characterize the atmosphere of WASP-33b (Cont et al. 2021(Cont et al. , 2022. The T -p profile is parametrized by a low pressure point (T 1 , p 1 ) and a high pressure point (T 2 , p 2 ). An isothermal atmosphere is assumed at pressures below p 1 or higher than p 2 . Between these two isothermal layers, we assumed a temperature that changes linearly as a function of log p. We used easyCHEM (Mollière et al. 2017) to calculate the mean molecular weight and the volume mixing ratios (VMRs) of the investigated chemical species (Fig. 2). For this purpose, we assumed equilibrium chemistry and a solar elemental abundance. An opacity grid of each species was computed for modeling the planetary spectra. The metal opacities were calculated from the Kurucz line list (Kurucz 2018) and the OH opacities were obtained from the MoLLIST line database (Brooke et al. 2016;Yousefi et al. 2018;Bernath 2020). Eventually, we ran the radiative transfer code petitRADTRANS (Mollière et al. 2019) to compute the model spectrum of each species. As the reduced spectra were continuum normalized, a normalization of the model spectra was also required. First, we calculated the planet-to-star flux ratio of the model spectra by dividing with the blackbody spectrum of the host star. The result was then normalized to the planetary continuum. As a last step, we convolved each model spectrum with the instrument profiles, obtaining the final model spectra for cross-correlation. The normalized model spectra of all investigated species are illustrated in Figs. 3 and 4.

Cross-correlation method
The cross-correlation method was applied to each species independently. We Doppler-shifted the model spectrum from -520 km s −1 to +520 km s −1 with steps of 1 km s −1 . Each ordermerged residual spectrum was multiplied with the shifted model spectrum and weighted by the uncertainties. This resulted in the weighted CCF, defined as The data points of the residual spectra weighted by the inverse of the squared uncertainties are denoted by r i , and m i is the model spectrum as a function of the Doppler-shift velocity . By combining the CCFs from the different datasets in a two-dimensional array, we obtained the final CCF map. In this step, we included the CCFs of those instruments that cover the wavelength ranges with significant emission features of the species in consideration. In contrast, we did not include the data of the instruments for which the model spectra predicted insignificant spectral features. The instruments and wavelength ranges included in the analysis and the corresponding chemical species are listed in Table 3. WASP-33 is a δ Scuti variable star, showing time-dependent luminosity variations and a variable stellar line profile (Herrero et al. 2011). Given that SYSREM corrects only efficiently for stationary features, the variable stellar lines are not entirely removed by the algorithm from the observed spectra. The presence of residual stellar lines of the same species as the one under investigation can lead to artifacts in the CCF map. In this case, the affected radial velocity (RV) range in the CCF map depends on the stellar rotation velocity and lies between ± rot sin i * in the stellar rest frame (i.e., ±87 km s −1 for WASP-33). We find these artifacts to be present in the CCF maps of Ti ii, V ii, Fe i, Fe ii, Si i, and Si ii, and therefore masked the affected RV range in the CCF maps of these species. For Ti i, V i, and OH, we find no residual stellar lines present because there are no artifacts in their CCF maps, nor is a substantial change in the detections of these species with or without masking achieved. No masking is therefore applied to the CCF maps of these three species. In Fig. 5 we show two example CCF maps, one with artifacts, the other without.
We assumed a circular orbit for shifting the CCF map of each species to the planetary rest frame. For this purpose, the CCFs were Doppler-shifted with a planetary RV of with K p the orbital velocity semi-amplitude, sys the velocity of the WASP-33 system, bary the barycentric correction, ∆ the velocity deviation from the planetary rest frame, and φ the orbital phase. We repeated the alignment over a range of different K p values. For each alignment, the CCF map was collapsed into a Article number, page 5 of 22  one-dimensional CCF by calculating the mean value over all outof-eclipse CCFs. The one-dimensional CCFs were then combined in a two-dimensional array. This array was normalized by its standard deviation under exclusion of the region around the strongest signal peak. In this way, we obtained a signal-tonoise map (S/N map) of the investigated chemical species. If the spectral signature is present in the data, the S/N map shows a significant peak at the location of the expected K p and close to ∆ =0 km s −1 .

Cross-correlation results and discussion
We find significant detections (> 4σ) of Ti i, V i, OH, Fe i, and Si i that correspond to the Doppler-velocity of WASP-33b's orbital motion. A clear peak is seen for all SYSREM iterations greater than one in the S/N maps at the location of the expected orbital parameters (Fig. 6). Figure 3 shows the S/N maps and the corresponding model spectra used for cross-correlation. We summarize the S/N of the detections and the measured orbital parameters in Table 3. This is the first detection of Ti i and V i in the atmosphere of WASP-33b. Previously, these two species were detected exclusively in transmission, making this their first detection in an    (Cont et al. 2021(Cont et al. , 2022 and are shown here for completeness. We note that the orbital parameters found for the individual species differ slightly, but within the one-sigma range. Detecting the emission signature from multiple species also unambiguously proves the existence of a temperature inversion layer in the planetary atmosphere, which is in agreement with the predictions from theoretical work (Hubeny et al. 2003;Fortney et al. 2008;Lothringer & Barman 2019). Figure 4 illustrates the cross-correlation results for the ionic species Ti ii, V ii, Fe ii, and Si ii. Emission lines of Ti ii, Fe ii, and, to a lesser extent, V ii are expected to be present near the blue end of the observed wavelength range. Only in the specific case of Ti ii and when considering the data from HARPS-N alone do we find tentative evidence for emission lines consistent with the planetary rest frame. Given the strong detection of Ti i and the higher spectral resolution of HARPS-N in compari- son to the other instruments used, we consider the tentative Ti ii peak as a real signal. On the other hand, no significant detections of V ii and Fe ii are achieved. We ran an injection-recovery test to investigate the nondetections of the two species. To this end, we Doppler-shifted the convolved model spectra used to attempt the detection of V ii and Fe ii with the reversed K p of WASP-33b ). Subsequently, we injected the shifted model spectra into the raw data and performed the preprocessing and cross-correlation analysis as described in Sects. 3 and 4.2. Doppler-shifting the injected signal with the reversed K p value of -231 km s −1 prevents interference with potentially undetected V ii and Fe ii signals from the planetary atmosphere. No recovery could be achieved for V ii, and the injected Fe ii model spectrum resulted in a negligible peak in the S/N map. We show the results of the injection-recovery test in Fig. 7. Given the poor recovery of the injected V ii and Fe ii signals, we conclude that our nondetections of these species are likely due to the relatively small number of prominent emission lines rather than resulting from their absence in the planetary atmosphere. Our nondetection of the Si ii signal was expected because the emission lines are very weak over the considered wavelength range. In summary, obtaining a cross-correlation signal of the metal ions studied is more challenging than detecting the corresponding neutral species because of the smaller number and strength of their emission lines, which are restricted to a relatively narrow wavelength range.
Our results show that atomic Ti and V are not rained out, cold trapped, or otherwise depleted in the hot atmosphere of WASP-33b. Moreover, the detections of Ti i, V i, and the tentative signal of Ti ii are consistent with the identification of these species in the atmospheres of a number of other strongly irradiated exoplanets (e.g., Hoeijmakers et al. 2018;Stangret et al. 2022;Prinoth et al. 2022;Bello-Arufe et al. 2022;Kesseli et al. 2022). We therefore conclude that the presence of significant Ti and V concentrations is favored at elevated temperatures. The injected V ii signal cannot be detected, and that of Fe ii causes a negligible peak in the S/N map. We conclude that the two species are not detectable even if they are present in the planetary atmosphere, given that the injected signals are poorly recovered.
Because emission spectroscopy preferentially probes spectral lines emerging from the hot planetary dayside, we propose this method to be particularly suitable for the search for refractory species such as Ti and V in exoplanet atmospheres. However, some observations suggest that the presence or absence of refractory species may not be determined by the atmospheric temperature alone. For example, Ti i was detected in the planetary atmosphere of HD 149026b (Ishizuka et al. 2021), the temperature of which is significantly below that of UHJs. On the other hand, the signature of the same species was not identified in the very hot atmosphere of KELT-20b/MASCARA-2b (Yan et al. 2022b). These observations suggests that physical parameters other than the temperature may also be important for the presence of Ti and V in the atmospheres of gas giant exoplanets.
We further applied the cross-correlation method to the datasets with a model spectrum that includes all the detected species (i.e., Ti i, V i, OH, Fe i, Si i, and Ti ii). The S/N map and the model spectrum are shown in Fig. 8. We find that using the model spectrum of all the species together results in an overall detection strength of S/N = 8.5 after eight SYSREM iterations (Fig. 6). A comparison with the S/N values from the individ-ual species in Table 3 shows that the spectral signature of Fe i is driving the detection. This is a reasonable finding, because the Fe i lines are expected to dominate the planetary emission spectrum over the entire wavelength range of our observations. The retrieved orbital parameters result in K p and ∆ values of 225.0 +3.0 −2.5 km s −1 and 0 +3 −2 km s −1 , respectively. Our K p is slightly lower than the expected value of 231 ± 3 km s −1 (Kovács et al. 2013;Lehmann et al. 2015), but is consistent with the results of previous studies (Nugroho et al. 2020a;Cont et al. 2021Cont et al. , 2022. The retrieved ∆ value is consistent with the orbital motion of the planet. This indicates that the planet may not have a strong dayside to nightside wind at the altitudes probed by these emission lines, which would cause a deviation of the detection peak of the order of a few km s −1 from the planetary rest frame. Figure 9 shows two versions of the CCF map, one aligned with the detected K p of 225 km s −1 , the other with the expected value of 231 km s −1 . The emission signal aligned with the detected K p value appears as a vertical trail with zero offset from the planetary rest frame. However, aligning the CCFs with the expected K p results in a tilted trail that is blue-and redshifted by a few km s −1 before and after secondary eclipse, respectively. Assuming that the expected K p reflects the true orbital velocity of the planet, we propose that these Doppler shifts are caused by the fast rotation velocity of WASP-33b of ∼ 7 km s −1 and the possible presence of super-rotation. In this scenario, the signature of the dayside atmosphere undergoes a spectral blueshift before eclipse and a redshift after eclipse (Fig. 9).

Retrieval of the atmospheric properties
We used the retrieval method developed by Yan et al. (2020) to constrain the properties of WASP-33b's atmosphere. This method has already been successfully used to determine the thermal structure in the atmosphere of two other UHJs, WASP-189b and KELT-20b/MASCARA-2b (Yan et al. 2020(Yan et al. , 2022bBorsa et al. 2022). In our implementation, we further developed this method for use with data from multiple instruments with different wavelength coverage.

Retrieval method
First, we calculated an individual master residual spectrum for each instrument. We used Eq. 2 to shift the residual spectra to the planetary rest frame. For aligning the spectra, we used the values of K p and ∆ measured for each instrument individually. The measurements were performed using the spectral model that consists of all the detected species. The corresponding S/N maps are reported in Fig. 10; the orbital parameters and SYSREM iterations used are summarized in Table 4. To obtain the master residual spectrum, we computed the average of the shifted residual spectra, weighted by the squared S/N of each exposure frame. As each spectrum corresponds to a different orbital phase value, the final result of our retrieval reflects the average atmospheric conditions over the observed phase interval. The master residual spectrum is still dominated by noise and contains the continuum normalized planet-to-star flux ratio.
In a second step, we defined the model spectrum for fitting with the master residual spectra from the four instruments. To this end, we used the radiative transfer code petitRADTRANS (Mollière et al. 2019). Our opacity grid covers a temperature range up to 25 000 K for all metal species. In the specific case of OH, the available partition function limited the opacity calculations to 5000 K. All OH opacities above this temperature were   for the line broadening. Following our description in Sect. 4.1, the model spectrum was further converted to the planet-to-star flux ratio and convolved with the instrument profile. The model spectrum was then interpolated to the wavelength solution of the master residual spectrum of each instrument. Finally, we fitted the spectral model to the master residual spectra. For each instrument, we used a standard Gaussian log likelihood function with R i the data points of the residual spectrum, σ i their uncertainties, β a scaling term to correct for possible overestimation or underestimation of the uncertainties, and m i the spectral model. We co-added the functions of the different instruments to get the combined log likelihood function of all the data. This results in an independent noise scaling factor for each instrument. For evaluating the combined likelihood function and estimating the fit parameters, we ran the Markov Chain Monte Carlo (MCMC) method from the emcee software package (Foreman-Mackey et al. 2013). Our retrieval includes the following free parameters: the temperature profile T 1 , p 1 , T 2 , p 2 ; a noise scaling parameter for each instrument β CV , β CN , β H , β E 1 ; the overall metallicity [M/H]; the average broadening velocity broad . For 1 Each instrument has an independent noise scaling factor. We use the following abbreviations: CARMENES VIS (CV), CARMENES NIR (CN), HARPS-N (H), and ESPaDOnS (E). each free parameter, we used 24 walkers with 4000 steps in the sampling.

Retrieval including all the detected species
We included the emission lines of all the detected chemical species from Sect. 4 into the retrieval (i.e., Ti i, V i, OH, Fe i, Si i, and Ti ii). Only the spectra with a planetary Doppler velocity outside the RV interval ± rot sin i * were used. In this way, we can be sure that the residual stellar lines do not overlap with the planetary spectral signature. The excluded spectra correspond to ∼10% of the cumulative exposure time. Table 5 summarizes the best-fit parameters resulting from the retrieval. The posterior distributions of the best-fit parameters are well constrained and are shown in Fig. A.1. Also, the retrieved noise scaling terms have values close to one, indicating an appropriate error estimation. Figure 11 illustrates the retrieved T -p profile. A more exhaustive presentation is provided in Fig. B.1. We find that the inversion layer extends over the pressure range 10 −5.1 bar to 10 −3.1 bar, with temperatures of 3981 +213 −108 K and 3424 +107 −111 K in the upper and lower planetary atmosphere, respectively. The thermal inversion layer is weaker when compared to the retrieval results of similar UHJs (Yan et al. 2020(Yan et al. , 2022b. This implies that the emission lines of WASP-33b are relatively shallow. We investigated whether or not the presence of H − opacity could explain the relatively low intensity of the emission lines, but found that the species has a negligible impact on our results (Fig. 12). To our knowledge, we performed the first retrieval of the thermal profile of WASP-33b's atmosphere using high-resolution spectroscopy. Similar work was recently carried out by van Sluijs et al. (2022), who maximized the S/N detection strengths of the CO emission lines to study the atmosphere of the same planet. We retrieve temperatures that are consistent with their results at low atmospheric pressures, but are able to get tighter constraints on the thermal profile at higher pressures.
The retrieval is mainly driven by the spectral signature of the neutral chemical species, which we assume are mostly ionized in the upper atmospheric layers due to the elevated temperatures. Nevertheless, we suggest that the emission lines of neutral species emerging from these layers are strong enough to allow the determination of the temperature T 1 at the upper limit of the thermal inversion. This suggestion is motivated by the fact that the posterior distribution of T 1 is well constrained in Fig. A.1 profile computed using all the spectra outside the RV interval ± rot sin i * . We also measured the temperature profiles by running the retrieval on two subsets of the spectral time series. The T -p profile inferred from the spectra far from secondary eclipse is indicated in pink, while that inferred from the spectra close to eclipse is in light blue. The right panel shows the orbital phase intervals used to compute the different T -p profiles (orbital phases in the RV range ± rot sin i * in transparent dark blue).
contrast, undetectable emission lines would cause a flat pattern toward higher temperatures in the posterior distribution, which is inconsistent with our results. Our retrieval constrains the atmospheric metallicity to [M/H] = 1.49 +0.82 −0.76 dex, which corresponds to a supersolar elemental abundance in the upper planetary atmosphere. This is in line with the results of van Sluijs et al. (2022), who obtain their strongest CO detection at abundances of ∼ 1 dex in the atmosphere of WASP-33b. Moreover, our result is consistent with previous studies that measured atmospheric abundances greater than solar for a number of hot giant exoplanets (e.g., Madhusudhan et al. 2014;Sedaghati et al. 2017;Nikolov et al. 2018). A correlation of the inferred [M/H] with the pressure parameters p 1 and p 2 can be recognized as diagonal distribution patterns in Fig. A.1. This is consistent with the expected degeneracy between [M/H] and the atmospheric temperature profile. The VMRs of the investigated species are computed by assuming that the elemental abundances all vary with the same metallicity value. The reason for this approximation is that spectra with very high S/N would be needed to successfully determine the abundances of individual species, a task that may be addressed in future studies.
We point out that our forward model approximates the VMRs of the different species by assuming equilibrium chemistry. Fossati et al. (2021) on the other hand suggested that nonlocal thermodynamic equilibrium (NLTE) may play an important role in the upper atmospheres of UHJs. NLTE is expected to alter the population levels and thus the VMRs of different chemical species. For example, NLTE is predicted to strongly affect the population levels of Fe by lowering the VMR of Fe i in favor of Fe ii. The alteration of the population levels has important consequences for the amplitude of the planetary emission spectrum, the strength of atmospheric absorption of incoming stellar light, and therefore on the atmospheric T -p profile. Given the great importance of Fe in the thermal inversion layers of UHJ atmospheres, future consideration of NLTE effects will increase the reliability of atmospheric retrievals.
We find that the spectral emission lines are significantly broadened. The retrieved Gaussian broadening profile has a standard deviation of broad = 1.9 ± 0.3 km s −1 , corresponding to a FWHM of 4.5 ± 0.7 km s −1 . The thermal and pressure broadening information is already included in the forward model via the opacities of the radiative transfer calculation. Therefore, broad is expected to account only for the broadening effects from atmospheric dynamics and the rotation of the planet. In particular, rotational broadening is supposed to affect the spectral line width, given the high rotational velocity of ∼ 7 km s −1 at WASP-33b's equator when assuming tidal locking. The line broadening from our retrieval is lower than that expected from a planetary sphere where the flux is emitted homogeneously over its entire surface, which would yield a FWHM of ∼ 9 km s −1 . We calculated this rotational broadening value with the PyAstronomy software package (Czesla et al. 2019). The retrieved line broadening therefore suggests that most of the contribution to WASP-33b's emission spectrum may originate from the hottest region of the planetary atmosphere, which is located close to the substellar point.
We tested whether our method is capable of detecting the presence of spectral line broadening. To this end, we applied our retrieval method to synthetic data. We first took the model spectrum without additional broadening as shown in Fig. 8. We then broadened the spectrum with broad = 2 km s −1 . Sequentially, we injected random white noise with different uncertainty levels to the model spectrum and performed the retrieval. We found that the correct broadening velocity can be retrieved even with the noise level of ten times the uncertainties of our observations (Fig. A.2). We therefore conclude that the quality of the observational data used is good enough to allow an appropriate determination of the spectral line broadening with our retrieval method.  Notes. In total, we ran six different retrievals. Three retrievals were performed with the spectral lines of all the detected species (i.e., Ti i, V i, OH, Fe i, Si i, and Ti ii). We used (i) all the spectra, (ii) a subset of spectra at orbital phases close to the secondary eclipse, and (iii) a subset of spectra far from the secondary eclipse. Another three retrievals were performed by including all the spectra, but only with the emission lines of (iv) Fe i, (v) Ti i + Ti ii combined, and (vi) V i, respectively. We conducted the retrievals by using uniform priors with the boundaries as follows.

Phase resolved retrieval
The observer's line of sight aligns with different geographical longitudes of the planet as a function of the orbital phase. Consequently, different regions of the planetary atmosphere are expected to contribute to the observed signal during our observations. Performing the retrieval on subsets of the spectral time series that correspond to different orbital phase intervals will allow us to study the physical conditions at different longitudes in the planetary atmosphere.
To perform a phase-resolved retrieval, we subdivided the spectral time series into two subsets for calculating the master residual spectra. One subset corresponds to the spectra close to the secondary eclipse when the dayside hemisphere faces the observer. The other subset consists of the spectra far from the secondary eclipse, and therefore contains the information from regions of both the planetary day-and nightsides. The orbital phase coverage of the two subsets is illustrated in Fig. 11. To avoid any interference with residual stellar lines, only the spectra at orbital phases outside the RV interval ± rot sin i * were used.
The posterior distributions of the fit parameters are shown in Figs. A.3 and A.4. Figure 11 compares the T -p curves obtained from the two subsets with the thermal profile obtained from all the spectra. A more detailed overview on the T -p profiles is also shown in Fig. B.1. Depending on the pressure level, the thermal profiles calculated from the two spectral subsets differ by about 300 K to 700 K, and the profile determined from all the spectra lies between the two. We find that the spectra close to the eclipse deliver a hotter temperature profile in comparison to those far from the eclipse. As the retrieved temperature profile corresponds to the average T -p curve of the visible hemisphere, this is consistent with the expectation of a temperature gradient away from the substellar point. Close to the secondary eclipse, mainly the dayside is facing the observer, leading to a hot average temperature profile. However, for spectra far from the secondary eclipse, the terminator region and a significant fraction of the cooler nightside are aligned with the observer's line of sight, resulting in lower atmospheric temperatures being measured.
The values of [M/H] and broad obtained for the two spectral subsets are consistent. We note a trend that spectra far from the secondary eclipse are less affected by line broadening than spectra near the eclipse. This is a reasonable result because far from the secondary eclipse, a significant fraction of the planetary disk that faces the observer is not illuminated by the host star. The nonilluminated planetary atmosphere is not expected to significantly contribute to the emission signal, decreasing the impact of the rotational line broadening. The dependence of the line broadening from the orbital phase can also be recognized when considering the CCF trail in Fig. 8, showing an increased width of the CCF toward the secondary eclipse. However, we note that the retrieved broad values should be compared with caution because of the relatively large uncertainties. All results of the phase resolved retrieval are listed in Table 5.

Retrieval of individual species
To compare the temperature profiles probed by individual chemical species, we ran retrievals for Ti, V, and Fe. These are the species with the most prominent detections in Sect. 4.2. Retrieving the T -p profile for the individual species also allows us to test the consistency of the results of our method. We set [M/H] to the previously determined value of 1.49 dex (cf. Sect. 5.2.1), because the metallicity was poorly constrained when leaving it as a free parameter. For the Ti retrieval, we included the opacities of both Ti i and Ti ii. As we were not able to constrain the presence of V ii and Fe ii by cross-correlation, only V i and Fe i were included in the other two retrievals, respectively. We used the spectra with a planetary RV outside ± rot sin i * for the retrievals including Ti ii and Fe i, which are affected by residual stellar lines. For V i, we included all the out-of-eclipse spectra because we did not find any significant residuals of stellar V i lines. The posterior distributions of the fit parameters are shown in Figs. A.5,A.6,and A.7. Figure 13 compares the retrieved Tp profiles and a detailed overview is provided in Fig. B.2. All results are summarized in Table 5. For Fe i, the parameters are well constrained and almost perfectly match the results from the retrieval that uses all the species in Sect. 5.2.1. We attribute this similarity to the larger number of Fe i lines when compared to that of all the other species. This also suggests that the retrieval that includes all the species is mainly driven by the Fe i lines. Also, the T -p profile from Ti i and Ti ii is in agreement with the retrieval result of all the species. The retrieval yielded a well constrained set of fit parameters, albeit at a somewhat lower precision than for Fe i. This is probably caused by the lower number and decreased strength of the Ti i and Ti ii spectral lines in comparison to the Fe i emission signature. The fact that retrievals for independent species achieve consistent results demonstrates the reliability of our method and gives confidence to the calculated parameter values.
The posterior distributions of the V i retrieval provide a significantly poorer constraint of the atmospheric parameters. This results in extended uncertainties of the T -p profile in the lower and upper planetary atmosphere. In view of the relatively low detection significance and the scarce number of features in the V i model spectrum in Fig. 3, this agrees with our expectations. While the temperature profile inferred from V i is close to that derived from Ti i/Ti ii and Fe i in the upper planetary atmosphere, we find a deviation at pressures higher than 10 −3 bar. Apart from the large uncertainties, the interpretation of this difference is unfortunately difficult, because both physical and methoddependent effects can cause the observed discrepancy. One possible scenario could be that the individual species are confined to different regions in the planetary atmosphere, and each has its own temperature profile. Figure 13 shows the chemical equilibrium VMRs of Ti i, Ti ii, V i, V ii, Fe i, and Fe ii as a function of temperature. We consider a pressure level of 10 −4 bar, which corresponds to the location of the thermal inversion layer. The VMRs of Ti i and Ti ii show a similar pattern when compared to that of V i and V ii, respectively. The colder temperatures probed by V i at higher pressure levels could therefore be explained by a lower abundance of the species than assumed by our method. In this case, the hottest region of the planetary atmosphere would be depleted by V i due to ionization and the hypothesized low abundance. Therefore, the V i signature would encode mainly the information from the coolest regions of the planetary dayside, leading to a more moderate T -p profile. The nominal abundance of Ti i and the inclusion of Ti ii would instead allow information from the entire dayside atmosphere to be considered, resulting in a hotter T -p profile. Alternatively, the use of oversimplified modeling of the VMRs could prevent an accurate estimation of the temperature profile for individual species. If the equilibrium chemistry assumption is not met, this could result in a poorly constrained temperature curve due to the degeneracy between temperature and VMR. Also, an oversimplification of the T -p model used and the presence of horizontal variations of the temperature profile from the substellar point to the nightside could explain the discrepancy measured. Spectra with increased S/N, the inclusion of nonequilibrium chemistry, and the use of a more comprehensive T -p model will enable us to better understand the information encoded in the spectral signature of individual species.

Conclusions
We used observations from the high-resolution spectrographs CARMENES, HARPS-N, and ESPaDOnS to analyze the thermal emission spectrum of the UHJ WASP-33b. A joint analysis of the data from the three instruments allowed us to investigate the planetary emission spectrum over an extended wavelength range from the near-ultraviolet to the NIR (3700-17 100 Å). We applied the cross-correlation technique to extract the faint spectral signature of the planetary atmosphere from the observations. This analysis led to the first detection of the emission signature of Ti i and V i in an exoplanet atmosphere, as these two species had previously been detected exclusively by transmission spectroscopy. Also, we detected a tentative emission signal of Ti ii. These detections are an important finding, given the frequently observed depletion of Ti-and V-bearing species in the atmospheres of UHJs. Moreover, we confirmed the presence of OH, Fe i, and Si i detected by previous studies. No significant signature from the ionic species V ii, Fe ii, or Si ii could be found in our spectral time series. The identification of spectral emission lines unambiguously proves the presence of a thermal inversion layer in the dayside atmosphere of WASP-33b, which is in line with theoretical work on highly irradiated planetary atmospheres.
We conducted a retrieval for the atmospheric T -p profile, the elemental abundances, and the spectral line broadening. The retrieval was performed using the data from all three instruments together. For this purpose, we forward modeled the emission lines of all the detected species via the radiative transfer code petitRADTRANS and assumed equilibrium chemistry. Compared to other UHJs (e.g., WASP-189b, KELT-20b/MASCARA-2b), our retrieval results in a relatively weak thermal inversion that extends roughly from 3400 K to 4000 K at pressures near 10 −3 bar and 10 −5 bar. We determined supersolar elemental abundances around 1.5 dex in the upper planetary atmosphere and found a spectral line broadening with a Gaussian FWHM of about 4.5 km s −1 . By running the retrieval on two distinct subsets of the spectral time series, we obtained temperature profiles that differ by 300 K to 700 K. This confirms the expectation that the temperature is higher on the planetary dayside than on the nightside. We also performed the retrieval for different chemical species individually. The temperature profiles from Ti i/Ti ii and Fe i are in good agreement with the overall result from all the species. However, the T -p profile of the V i signature slightly deviates from that of the other species. We suggest that a V idepleted planetary atmosphere could explain the measured discrepancy between the T -p profiles.
Our work shows that high-resolution emission spectroscopy offers the possibility to study the physical conditions in UHJ atmospheres in great detail. Further progress on atmospheric retrievals will be achieved by deploying more sophisticated spectral forward models, increasing the number of species included, and expanding the orbital phase coverage of UHJ observations.