High-resolution Transmission Spectroscopy of Ultrahot Jupiter WASP–33b with NEID

We report an attempt to detect molecular and atomic species in the atmosphere of the ultrahot Jupiter WASP-33b using the high-resolution echelle spectrograph NEID with a wavelength coverage of 380–930 nm. By analyzing the transmission spectrum of WASP-33b using the line-by-line technique and the cross-correlation technique, we confirm previous detection of Hα, Hβ, Hγ, and Ca ii infrared triplets. We find no evidence for a significant day-to-night wind in WASP-33b, taking into account the effects of stellar pulsations using a relatively novel Gaussian process method and poorly constrained systemic velocity measurements. We also detect the previously reported pretransit absorption signal, which may be a pulsation mode induced by the planet. Combined with previous CARMENES and HARPS-N observations, we report the nondetection of TiO, Ti i, and V i in the transmission spectrum, while they were already detected in the dayside atmosphere of WASP-33b. This implies a difference in the chemical compositions and abundances between the dayside and terminator atmospheres of WASP-33b and certainly requires further improvements in the sensitivity of the detection methods.


INTRODUCTION
Ultra-hot Jupiters (UHJs) are gas giant planets with high dayside temperatures (T day ≥ 2200 K, Parmentier et al. 2018) that orbit very close to their host stars.The characterization of the atmospheres of UHJs has been a rapidly growing field in exoplanet science over the past five years.Theoretical studies of UHJ atmospheres suggest that the atmospheres may be dominated by atoms and ions rather than molecules due to thermal dissociation and ionization (Lothringer et al. 2018;Parmentier et al. 2018;Arcangeli et al. 2018).In addition, because their atmospheres are in an extreme state, some extreme phenomena may occur, such as atmospheric evaporation and escape (e.g., Yan & Henning 2018;Sing et al. 2019).Hydrodynamic atmospheric escape driven by intense stellar irradiation (extreme-ultraviolet and/or nearultraviolet) would lead to significant mass loss, which would Corresponding author: Guo Chen guochen@pmo.ac.cn affect the long-term evolution of the planet (Fossati et al. 2018;García Muñoz & Schneider 2019).And UHJs are strongly irradiated and tidally locked by their host stars, which can lead to significant differences in the chemical composition and abundances of their atmospheres between the day and night sides (Arcangeli et al. 2018;Bell & Cowan 2018;Helling et al. 2019;Parmentier et al. 2018).
High-resolution transmission spectroscopy has been widely used to characterize the atmospheres of UHJs, in particular to reveal their chemical composition.For example, Yan & Henning (2018) found the excess Hα absorption in KELT-9b (Gaudi et al. 2017), revealing an extended hydrogen atmosphere.Hoeijmakers et al. (2018) further reported the presence of neutral and singly ionized atomic iron (Fe i and Fe ii) and singly ionized atomic titanium (Ti ii) in its atmosphere.Prinoth et al. (2022) presented the first unambiguous detection of TiO on an exoplanet's terminator, which plays a crucial role in the formation of thermal inversions (Hubeny et al. 2003;Fortney et al. 2008), based on transmission spectroscopy of WASP-189b.Deibert et al. (2021) and Casasayas-Barris et al. (2021) found the presence of ionized calcium (Ca ii) at the terminator of WASP-76b by studying the Ca ii infrared triplet (∼850 nm).Neutral Fe i has also been detected on this planet (Ehrenreich et al. 2020).The atmosphere of MASCARA-2b (also known as KELT-20b;Talens et al. 2018;Lund et al. 2017) has been studied by Casasayas-Barris et al. (2018, 2019), and the hydrogen Balmer series (Hα, Hβ, and Hγ) are detected along with Ca ii, Fe ii, and Na i. Hydrogen Balmer and Ca ii lines can be used as probes to study atmospheric escape and the dynamics of the upper atmosphere.WASP-33b (Collier Cameron et al. 2010) is a UHJ that has attracted considerable attention in recent years, orbiting a δ-Scuti A-type star with a period of 1.22 days.Its host star is very bright (V ∼ 8 mag), making it an ideal target for groundbased high-resolution observations.The dayside atmosphere of WASP-33b has been well studied.Several studies have found evidence of a temperature inversion in its dayside atmosphere (e.g., Haynes et al. 2015;Nugroho et al. 2017;Cont et al. 2021).Species such as Si i, Fe i, Ti i, Ti ii, V i, OH, CO, H 2 O, and TiO have been detected in the dayside atmosphere by high-resolution emission spectroscopy (Nugroho et al. 2020(Nugroho et al. , 2021;;Serindag et al. 2021;Cont et al. 2021Cont et al. , 2022a,b;,b;Herman et al. 2022;van Sluijs et al. 2022;Finnerty et al. 2023).In contrast, studies of the terminator atmosphere have lagged behind due to contamination from stellar pulsations.Its host star has a complex pulsation mode, with certain modes possibly being excited by the planet (Herrero et al. 2011;Borsa et al. 2021a).Nevertheless, the hydrogen Balmer lines and Ca ii have been identified using highresolution transmission spectroscopy (Yan et al. 2021;Borsa et al. 2021a;Cauley et al. 2021;Yan et al. 2019).In addition, von Essen et al. (2019) claimed the detection of aluminum oxide (AlO) using a low-resolution optical transmission spectrum at the terminator of WASP-33b.
To date, many ground-based high-resolution spectrographs have been used for transmission spectroscopy of the hot Jupiter and ultra-hot Jupiter atmospheres.In this paper, we present a new study of the atmosphere of WASP-33b by transmission spectroscopy using the stable high-resolution Doppler spectrograph NEID (Schwab et al. 2016), which was designed to detect exoplanets with extreme precision radial velocity and newly commissioned in late 2019.With its high spectral resolution and stability, NEID has the potential to unravel the atmospheric composition of exoplanets.In addition, we also combine previous observational data in order to obtain more robust results.
This paper is organized as follows.We describe the observations and data reduction in Sect. 2 and 3. We summarize our methods, including the line-by-line analysis and the cross-correlation analysis, in Sect. 4. We present and discuss the resulting transmission spectra in Sect. 5.The conclusions are given in Sect.6.

OBSERVATIONS
We observed one transit of WASP-33b on October 16, 2021 using the NEID (Schwab et al. 2016) spectrograph on the WIYN 3.5m Telescope at Kitt Peak Observatory1 .NEID covers the wavelength range from 380 to 930 nm at a resolving power of R ∼ 110, 000 or ∼2.7 km s −1 (Halverson et al. 2016).The exposure time was set to 300 s.We obtained 65 spectra, but discarded four due to low signal-to-noise ratio (S/N).The average S/N of the remaining 61 spectra is around 98 at 5500 Å.The raw data were reduced using the NEID Data Reduction Pipeline2 (NEID-DRP).After data reduction of the raw frames, we obtained the order-by-order spectra.These spectra are at the vacuum wavelength and in the terrestrial rest frame.For the convenience of the study, we converted the wavelengths to air wavelengths.Although the NEID wavelength range covers down to 380 nm, the S/N of the bluer orders is too low to be used for transmission spectroscopy.
In addition, we collected archival data for four transits of WASP-33b from the Calar Alto Archive3 and the Italian center for Astronomical Archive4 .These were used by Yan et al. (2021) to detect the hydrogen Balmer lines.Two of them were observed on January 5 and January 16, 2017, with the CARMENES (Quirrenbach et al. 2018) spectrograph installed on the 3.5 m telescope at the Calar Alto Observatory.The CARMENES visual channel has a wavelength coverage of 520-960 nm with a high spectral resolution (R ∼ 94, 600).The CARMENES data were reduced using the CARACAL pipeline (Caballero et al. 2016).Similar to NEID, we also converted the wavelengths of the order-byorder one-dimensional spectra to air wavelengths.The first 19 spectra with shorter exposure times (<120 s) were discarded for the first night.Weather conditions were ideal on the first night and overcast on the second night.
The remaining two transits were observed on October 17 and November 8, 2018, with the HARPS-North (HARPS-N; Cosentino et al. 2012) spectrograph on the Telescopio Nazionale Galileo (TNG).This spectrograph has a high resolution of R ∼ 115, 000 and covers the wavelength range from 383 to 690 nm.The raw data were processed using the HARPS-N pipeline (HARPS-N Data Reduction Software, version 3.7).We used the order-merged onedimensional spectra with a wavelength step of 0.01 Å from the pipeline.Since the raw wavelength step is oversampled, we re-binned the spectrum so that the wavelength step becomes 0.03 Å (as described in Yan et al. 2021).Unlike NEID and CARMENES, the HARPS-N spectra were at air wavelength and corrected for the barycentric Earth radial velocity (BERV).In terms of data quality, the second night was better than the first because the observed flux on the first night (October 17, 2018) dropped significantly when the telescope was pointed near the zenith (Yan et al. 2021).
The observation log summary, including our NEID observation and four publicly available observations, is presented in Table 1.It should be noted that the CARMENES and HARPS-N observations were used only for the crosscorrelation search and not for the line-by-line analysis.

DATA REDUCTION
For each NEID spectral order, we first performed a blaze correction using the blaze function provided by the NEID pipeline.We then combined all the spectral orders from each order-by-order spectrum and resampled them to a uniform grid with a wavelength step of 0.02Å to obtain an ordermerged one-dimensional spectrum.To correct for outliers that are unlikely to be caused by the observed source, we flagged the pixels that were greater than ten standard deviations from the flux values of the smoothed spectrum produced by a median filter with a window size of 9 pixels.The flux values of these pixels were replaced by the median of the ten surrounding pixels (Langeveld et al. 2021).We normalized each spectrum and corrected for the continuum variation following Chen et al. (2020).Specifically, we first averaged all the out-of-transit spectra to create a master-out spectrum as a preliminary reference spectrum.We then fitted the ratio of each individual spectrum to the reference spectrum with a fourth-order polynomial function.Finally, we divided each spectrum by the best-fit polynomial function to correct for the continuum variation.After this data reduction, we obtained the final one-dimensional normalized spectra for the following analysis.
Similar to the data reduction for NEID, we obtained the order-merged normalized CARMENES spectra.We did not apply any blaze correction to these spectra as this was already done by the CARACAL pipeline.In addition, for HARPS-N, we do not need the merging operation since the HARPS-N spectrum (s1d) provided by the pipeline is one-dimensional, but the other procedures are consistent with NEID.

METHODS
In this section, we apply two methods commonly used in transmission spectroscopy to search for atmospheric species, the line-by-line analysis and the cross-correlation analysis.

Creation of transmission spectral matrix
We first used version 1.2.0 of Molecfit (Smette et al. 2015;Kausch et al. 2015) to perform the telluric correc-tion.Molecfit is an ESO tool for correcting telluric features in ground-based spectra.Molecfit produces a synthetic model of the Earth's atmospheric absorption for each spectrum based on a line-by-line radiative transfer mode (LBLRTM).Unlike CARMENES and NEID, the HARPS-N spectra are given in the Solar system barycentric rest frame.Therefore, we first shifted each HARPS-N spectrum to the terrestrial rest frame using the BERV values.The parameters used to build the telluric spectrum are similar to those of Allart et al. (2017), but are slightly adjusted to account for instrumental differences.In addition, we carefully selected the wavelength regions containing only the separated unsaturated strong lines of telluric H 2 O or O 2 to fit for the Earth's atmospheric absorption.
After the telluric correction, we shifted all the spectra to the stellar rest frame by considering the BERV values and the stellar systemic velocity ( sys ).We then computed a new master-out spectrum by combining all the continuum variation and telluric corrected out-of-transit spectra that have been corrected for continuum variation and telluric using the square of the S/N as the weight.We divided each individual spectrum by the master-out spectrum to remove the stellar lines.The spectrum in which the stellar lines are removed, called the residual spectrum, contains planetary signals and noise.We applied a Saviztky-Golay filter using savgol filter from the scipy package with a window size of 401 to remove low frequency variations on the continuum of the residual spectrum, which could be caused by stellar pulsation, blaze variation correction, or instrumental stability.Finally, we arranged all the residual spectra by phase to construct a transmission spectral matrix.

Correction of RM and CLV effects
The residual spectra still contain the contamination from variations in the stellar line profile caused by additional systematic effects during the transit.Two important effects are the Rossiter-McLaughlin (RM) effect (Rossiter 1924;McLaughlin 1924;Queloz et al. 2000;Louden & Wheatley 2015) and the center-to-limb variation (CLV) effect (e.g., Yan et al. 2015Yan et al. , 2017;;Czesla et al. 2015;Cegla et al. 2016;Casasayas-Barris et al. 2018;Yan & Henning 2018;Chen et al. 2020).A strong RM effect is expected in the transmission spectrum of WASP-33b because its host star is a rapidly rotating A-type star.Correction for the RM effect is therefore critical to the transmission spectroscopy study of WASP-33b.
In addition, WASP-33b's orbit is subject to nodal precession, and the orbit change rates have been measured using Doppler tomography (Johnson et al. 2015;Iorio 2016;Watanabe et al. 2020).Measuring the orbital inclination (i) and the projected spin-orbit misalignment angle (λ) at the epoch of each observation is difficult under the effects of the stellar pulsation and is beyond the scope of this study.Therefore, based on the epoch difference between the reference epoch ( 2014) and each observation epoch, we adopted the orbital change rates measured by Johnson et al. (2015) to calculate the orbital inclination and spin-orbit misalignment angle at the epoch of each observation.This strategy has been used in previous studies such as Yan et al. (2019Yan et al. ( , 2021)).The stellar and planetary parameters of the WASP-33 system are listed in Table 2.
To correct for the CLV and RM effects, we followed the method described in Casasayas- Barris et al. (2019) to construct the CLV and RM model.In this method, the stellar spectra at 21 different limb-darkening angles (µ = cos θ) are computed with the Spectroscopy Made Easy (SME) tool (Valenti & Piskunov 1996;Piskunov & Valenti 2017), using the VALD3 line lists (Ryabchikova et al. 2015) and the Kurucz ATLAS9 models (Castelli & Kurucz 2003).The stellar spectrum calculation assumes local thermodynamic equilibrium (LTE) and solar abundance.However, we found that the synthetic stellar spectrum is slightly different from the observational data.In this case, we fitted the abundances of some crucial elements (e.g., H for the analysis of Hα, Hβ, and Hγ lines) using the SME tool instead of directly assuming the solar abundance.We roughly assumed that the size of the planetary disk was equal to the apparent radius of the planet, and did not account for the increase in radius due to the planetary atmosphere.Using this method, we obtained a combined CLV and RM model and normalized it by creating a master-out CLV+RM model, similar to the creation of the master-out spectrum in Sect.4.1.1.
For the case of WASP-33b, we assumed a solar abundance as input for each species and an isothermal atmosphere at 2700 K, which is close to the equilibrium temperature of Effective temperature Projected rotation velocity Planetary parameters Planetary surface gravity log g p [cgs] 3.46 +0.08 Orbital semi-major axis a [R ⋆ ] 3.69 ± 0.05 (a)   Orbital period 1.219870897 (a)   Transit epoch (BJD) T 0 [d] 2454163.22449 (a)   Transit depth 2017 January 5 2017 January 16 References-(a) Yan et al. ( 2019). (b) Kovács et al. (2013). (c)  Adopted from Yan et al. (2019).The changes of the orbital parameters between the two HARPS-N or two CARMENES observations are negligible, as the dates are very close. (d) Predicted value using parameters in Johnson et al. (2015).
WASP-33b.The planetary radius and surface gravity we used are listed in  2019), the atmosphere becomes opaque between 1 and 10 mbar due to H − absorption in ultra-hot Jupiters.Therefore, according to the nature of WASP-33b, we set the continuum level to 10 mbar.Before the spectral models were used in the crosscorrelation analysis, they were convolved with a Gaussian function to match the instrument resolution.The spectral models were normalized with a Gaussian high-pass filter of 12 points to remove the continuum level.With these procedures, the spectral models have values around one and are ready for use in the cross-correlation analysis.

Cross-correlation
For efficiency and simplicity, we used the SYSREM algorithm (Tamuz et al. 2005), which is widely used in crosscorrelation analysis, to reduce the telluric and stellar contamination and obtain the observed transmission spectrum.The cross-correlation analysis is different from the line-byline analysis, but the two methods produce consistent results.
We first constructed the spectral matrix by collecting the one-dimensional normalized spectra of each observation after data reduction (see Sect. 3), and then ran SYSREM on the spectral matrix.In one iteration, SYSREM captures the systematics by iteratively fitting each wavelength bin in the spectral matrix, resulting from variations in telluric absorption, stellar lines, and instrumental effects.SYSREM then removes these systematics from the spectral matrix.Through multiple iterations, SYSREM can effectively remove the telluric and stellar lines and obtain the residual transmission spectral matrix containing the planetary features and noise.We implemented SYSREM following Gibson et al. (2020) by dividing the spectral matrix by the systematics rather than subtracting the systematics from the spectral matrix.We ran the algorithm through ten iterations and selected the residual spectral matrix of the iteration with the highest S/N as the algorithm output.This matrix was then used in the crosscorrelation calculations.
We implemented the cross-correlation analysis (Snellen et al. 2010) for each species independently.The spectral models were shifted from −250 to +250 km s −1 with a step of 1 km s −1 .At each shift, the weighted residual spectrum was multiplied by the shifted spectral model to obtain a weighted cross-correlation function (CCF).The CCF is given by where r i is the residual spectrum, m i is the spectral model shifted by velocity , and σ i is the error at wavelength point i.
We calculated the variance of each wavelength bin (i.e., along the column) and the variance of each frame (i.e., along the row) in the residual spectral matrix, and then took the square root of their sum as the error.For each residual spectrum, we obtained a CCF.We then stacked all the CCFs to generate the CCF map for each observation night and spectral model.
To detect planetary signals, we used a grid of assumed planetary orbital velocity semi-amplitudes (K p ) to convert the CCF map from the terrestrial rest frame to the planetary rest frame.Assuming the planet has a circular orbit, the radial velocity (RV) of the planet p is defined as where sys is the systemic velocity, bary is the barycentric Earth radial velocity (i.e., BERV), ∆ is the radial velocity deviation from zero value, and ϕ is the orbital phase.In general, ∆ can be closely related to the dynamical properties of the planetary atmosphere (e.g., wind).For each value of the assumed K p , we generated a one-dimensional CCF by averaging the CCF map of the in-transit portion, which had been shifted to the planetary rest frame by that value.We then created a K p map by stacking the one-dimensional CCFs of the different K p values ranging from 0 to 300 km s −1 with a step of 1 km s −1 .The K p map was expressed in the form of S/N, with the CCF values being normalized by the standard deviation calculated from +250 to +300 km s −1 in K p and +200 to +250 km s −1 in ∆ .

Species detection using the line-by-line technique
We present here the results obtained by analyzing individual lines in different wavelength regions.Using the NEID data, we confirm previous detections of the Balmer series of H (Hα, Hβ, and Hγ) and Ca ii.We also try to analyze the spectral lines of several atoms and ions expected to be present in WASP-33b's atmosphere.However, since almost all individual lines of atoms and ions are affected by stellar pulsations, we cannot distinguish between planetary atmospheric absorption signals and stellar pulsation signals.Therefore, we could not report any additional new significant detections of atoms and ions in WASP-33b using line-by-line analysis.While in principle a large number of additional transit observations or pulsation modeling with simultaneous (spectro-)photometric monitoring can help to eliminate the effects of pulsations, we decided to treat stellar pulsations as correlated noise in our analysis, as will be detailed in this section.

Hα transmission spectrum
The first to third rows of Figure 1 shows the transmission spectral matrix and the stellar CLV and RM effects around the Hα line.The planetary Hα absorption is clearly visible in the phase-resolved transmission spectra observed with NEID after correction for the CLV and RM effects.The velocity shifts of the absorption are fully consistent with the orbital motion of the planet during the transit.We removed the CLV and RM effects from the transmission spectral matrix, and shifted the residual spectra to the planetary rest frame using the expected K p = 231 km s −1 calculated with the planetary orbital parameters.We then averaged all the in-transit spectra with the square of the S/N as weight to derive the final one-dimensional transmission spectrum.The phase-folded Hα transmission spectrum is shown in the fourth row of Figure 1.We fit a Gaussian profile to the one-dimensional Hα transmission spectrum, assuming that the planetary absorption signal has a Gaussian profile controlled by the line contrast (h), the velocity offset of the center of the observed line (V center ), and the full width at half maximum (FWHM).
The best-fit parameters derived by this classical white-noise based method are presented in Table 3.
The center of the Hα line measured with NEID has a velocity offset of 1.9 ± 0.7 km s −1 .Velocity offsets inferred by transmission spectroscopy could be due to planetary winds (Wyttenbach et al. 2015;Louden & Wheatley 2015;Brogi et al. 2016).Based on HARPS-N measurements, Borsa et al. (2021a) reported a significant blueshift of −8.2 ± 1.4 km s −1 for Hα and attributed it to the presence of winds in the atmosphere of WASP-33b.However, Yan et al. (2021) found no evidence for winds (velocity offset of 1.2 ± 0.9 km s −1 ) by combining HARPS-N and CARMENES observations.Our results are more consistent with the latter.In the case of WASP-33b, these contradictory results are not surprising because of the correlated noise introduced by stellar pulsations.The pulsations would temporarily alter the stellar line profile and make it difficult to remove the stellar lines using the master-out spectrum, leaving strong residuals on top of the planetary absorption.Such correlated noise could introduce spurious velocity offsets when fitting a Gaussian function to the expected planetary absorption line.
Meanwhile, because its host star is a rapidly rotating Atype star, it is extremely difficult to accurately measure the stellar systemic RV.Without a precise RV value, we are cautious about interpreting the measured velocity offset as the winds in the planetary atmosphere.We conclude that our NEID data do not show a significant velocity shift, but we cannot exclude the possibility that winds exist at the terminator of WASP-33b.
The FWHM of the Hα line is measured to be 30.5 ± 1.7 km s −1 , which is slightly lower than the results of Borsa et al. (2021a) and Yan et al. (2021).The difference could be due to instrumental effects, different data reduction procedures (e.g., removal of telluric and stellar lines, correction of CLV and RM effects, and normalization), and stellar pulsations.Given the visible variations of stellar pulsations on different nights, stellar pulsations could be the more likely origin.Such a large FWHM is the result of a combination of mechanisms, including rotational broadening, thermal broadening, and dynamical processes, suggesting that excited hydrogen atoms may be present in the upper atmosphere and likely undergo hydrodynamic escape.Note-* The results of the classical method of CARMENES and HARPS-N are obtained from Yan et al. (2019Yan et al. ( , 2021)).
The measured line contrast of the Hα line is 0.91 ± 0.04 %.Following Chen et al. (2020), we derived the effective planetary radius (R eff ) defined by: R 2 eff /R 2 p = (δ + h)/δ, where δ = (R p /R ⋆ ) 2 is obtained from Table 2 and h is the measured line contrast.For the Hα line, the effective planetary radius is 1.302 +0.012 −0.012 R p , which is below the effective Roche radius (1.71 +0.08 −0.07 R p for WASP-33b).Therefore, the observed Hα absorption signal originates from the gravitationally bound atmospheric envelope within the Roche lobe of the planet.Consequently, we will not be able to observe out-of-transit absorption like Nortmann et al. (2018) and Ehrenreich et al. (2015) when planetary atmospheric materials do not escape beyond the Roche lobe.We note that our measured line contrast is consistent with Yan et al. (2021) (0.99 ± 0.05 %), but is almost twice the value obtained by Borsa et al. (2021a) (0.54 ± 0.04 %) and significantly lower than a third measurement (1.68 ± 0.02 %) reported by Cauley et al. (2021) using the PEPSI spectrograph mounted on the Large Binocular Telescope.As with the FWHM, we speculate that the discrepancies in the measured line contrast are caused by different levels of stellar pulsations.
In order to obtain proper parameter estimates for the line profile, we decided to adopt Gaussian process (GP) regression in the line profile modeling.GP can capture the correlation between data points based on a covariance matrix, where data points that are close to each other in input space are highly correlated.Since its first implementation in lowresolution transmission spectroscopy by Gibson et al. (2012), GP has been widely used in time-series exoplanet data to account for correlated noise, and has recently been proposed to treat telluric and planetary signals for high-resolution transmission spectroscopy (Meech et al. 2022).Here we used GP to account for the correlated systematics caused by stellar pulsations.
We assumed that the one-dimensional transmission spectrum (R) follows a multivariate Gaussian distribution: where M is the mean function with the parameter vector θ, C is the covariance matrix with the parameter vector φ, and λ is the wavelengths.The elements in C are determined by a kernel function k τ i j ; φ , where τ i j ≡ λ i − λ j is the wavelength interval between two data points.Here we assumed that the mean function M is the Gaussian function to model a pure planetary absorption profile.We implemented GP with the python code celerite (Foreman-Mackey et al. 2017).The kernel function consists of a 3/2-order Matérn kernel and a diagonal jitter ("white noise") term.Finally, we used the Markov chain Monte Carlo (MCMC) sampling algorithm implemented in emcee (Foreman-Mackey et al. 2013) to find the best-fit parameters and their uncertainties.The best-fit parameters and adopted priors are presented in Table 3 and Table 4, respectively.The blue line in the fourth row of Figure 1 shows the bestfit model derived from the GP method, where the GP does a good job of accounting for the variation caused by stellar pulsations.The GP method derived an FWHM of 42.8 +3.1 −3.2 km s −1 , a line contrast h of 1.24 +0.20  −0.19 % (R eff = 1.396 +0.054 −0.053 R p ), and a velocity offset of 1.1 +1.6 −1.6 km s −1 .Compared to those derived from the classical white-noise based method, the GP-derived parameters are more robust with larger error bars after accounting for the correlated noise.
The purple and cyan lines in the fifth row of Figure 1 show the hydrostatic models of Hα with 2,710 K (i.e., the equilibrium temperature of WASP-33b) and 10,000 K, respectively.The planetary rotation broadening is considered assuming tidal locking using the method of Boucher et al. (2023).Each hydrostatic model is calculated by petitRADTRANS assuming an isothermal temperature structure, chemical equilibrium, and other settings consistent with Sect.4.2.1.The observed Hα excess absorption is significantly stronger than the equilibrium-temperature model prediction and slightly stronger than the 10,000 K model prediction.We note that the maximum temperature allowed by our adopted H line opacity is 6,000 K, which might lead to the underestimation.On the other hand, the discrepancy could be due to the lack of consideration of non-local thermodynamic equilibrium (NLTE) effects and/or hydrodynamic escape (e.g., García Muñoz & Schneider 2019;Fossati et al. 2021Fossati et al. , 2023;;Huang et al. 2023).
We also applied the GP method to the Hα transmission spectra derived in Yan et al. (2021).As shown in Table 3, the results obtained by the classical method and the GP method are very consistent for the Hα transmission spectra of Yan et al. (2021).This is probably because the line core of the planetary absorption profile is less affected by stellar pulsations in the corresponding observations with different phasing.
We note that there is an absorption feature moving at a velocity semi-amplitude of ∼460 km s −1 (approaching twice the expected K p of the planet) in the pre-transit phase, as shown at the third panel of Figure 1.When shifted to the planetary rest frame, as shown in Figure 2, the pre-transit absorption feature is located at red-shifted velocities of more than 100 km s −1 .Similar pre-transit absorption features have been reported by Borsa et al. (2021a) and Cauley et al. (2021).Borsa et al. (2021a) found that the pre-transit absorption features did not cancel out like other pulsations when averaging different transits, and suggested that they could be stellar pulsation mode excited by the planet.The evidence of tidally perturbed stellar pulsation modes has also been found in the light curve of WASP-33 observed by the Transiting Exoplanet Survey Satellite, which were attributed to the misalignment of the stellar rotational axis and the planetary orbit (Kálmán et al. 2022).

Hβ and Hγ transmission spectra
The NEID spectrograph also covers the Hβ and Hγ lines.The best-fit parameters of the Gaussian function using the classical white-noise based method and the GP method are shown in Table 3.The Hβ and Hγ transmission spectral matrices and one-dimensional transmission spectra are shown −0.038 R p , and 1.163 +0.048 −0.054 R p , respectively.The velocity shifts of the Hβ and Hγ lines are consistent with planetary orbital motion, but their absorption features are weaker than the Hα line.Moreover, their absorption depths are significantly stronger than hydrostatic model predictions at the equilibrium temperature.The case of the hydrogen Balmer series has also been detected in KELT-9b (Cauley et al. 2019;Wyttenbach et al. 2020) and KELT-20b (Casasayas-Barris et al. 2019).

Ca ii transmission spectra
Although the Ca ii H&K lines (K line 3933.661Å, and H line 3968.467Å at air wavelength) induced by resonant transitions from the Ca ii ground state are covered by the wavelength region of the NEID spectrograph, the S/N of these lines is too low to extract the absorption signals from the planetary atmospheric layer.Therefore, we only analyze the  3 shows the best-fit parameters of the Gaussian function based on the classical method and the GP method, respectively.The mean velocity offset of the Ca ii lines is 0.51 +1.36  −1.32 km s −1 for the GP method.This result is consistent with the Balmer lines when stellar pulsations are treated as correlated noise by the GP, indicating that no strong day-to night-side winds are observed.The mean FWHM of the Ca ii IRT lines (29.7 +3.6 −3.3 km s −1 ) is consistent with previous results.We translate the line contrast of Ca ii into an effective planetary radius of R eff = 1.143 +0.033 −0.034 R p for Ca ii IRT (λ8498), R eff = 1.264 +0.045 −0.046 R p for Ca ii IRT (λ8542), and R eff = 1.273 +0.036 −0.037 R p for Ca ii IRT (λ8662), respectively.In addition, similar to the Balmer lines, the Ca ii absorption is also significantly stronger than predicted by the equilibrium temperature hydrostatic model.This is consistent with Yan et al. (2019).

Species detection using the cross-correlation technique
We used the cross-correlation technique to search for some species that are expected to be present at the terminator of WASP-33b, in particular TiO, Ti i, V i.These have been detected in the dayside atmosphere of WASP-33b.In general, materials from the planet's dayside atmosphere can be transported through the terminator to the nightside by certain dynamical processes.Thus, if species are detected in the dayside atmosphere, we are possibly to find them in the terminator and nightside atmosphere.In addition, we looked for some other atoms and ions that are expected to be present in ultra-hot Jupiter.Unfortunately, these atoms and ions are affected by the pulsations of WASP-33, and we have diffi- culty distinguishing between pulsation signals and planetary atmospheric absorption signals.

Presence of Ca ii confirmed in CCF
In Sect.5.1.3,we confirmed the presence of Ca ii in the atmosphere of WASP-33b using the line-by-line technique with the NEID data.Here we combined the five observations (see Table 1) to obtain the CCF and K p -maps for Ca ii.To reduce the CLV+RM effects, we applied the CLV+RM model to correct each residual spectrum before cross-correlation.The corrected maps are shown in Figure 5.
For WASP-33b, the Ca ii signal is very strong, and we can observe it directly from the CCF map (see the top panel of Figure 5).The K p map is shown in the middle panel of Figure 5.This map shows a strong cross-correlation signal at the expected values with an S/N of 14.2.The bottom panel of Figure 5 shows the CCFs at the K p = 225 km s −1 , corresponding to the S/N peak.The peak of the CCFs located at ∆ ∼ 0 km s −1 indicates the absence of wind in the atmosphere of WASP-33b.
We derived a K p value of 225 +11 −21 km s −1 using the Ca ii absorption lines, including the H&K lines and the IRT lines.The derived K p value using planetary Ca ii absorption is in agreement with the expected K p value ( 231 −2.5 km s −1 using the planetary emission lines (including the lines of Ti i, V i, OH, Fe i, Si i, and Ti ii) in WASP-33b and speculated that the discrepancy between the derived value and the expected value may be due to the super-rotation of WASP-33b.
In addition, the pulsations of WASP-33b's host star also affect the constraint on the K p value.This is because pulsa- whelmed by noise, making it difficult to detect species in the planetary atmosphere.For Ca ii, the strong cross-correlation signal is located at the planetary velocity and appears only during transit, it is unlikely to be induced by pulsation.However, some species with weak features in WASP-33b's atmosphere or with high abundances in WASP-33b's host star are difficult to detect at the terminator using transmission spectroscopy, especially for metal atoms and ions, as will be discussed in the following text.

Non-detection of TiO
The detection of TiO is not expected to be confounded by stellar pulsations.Due to the extremely high temperature of the host star's atmosphere, molecules such as TiO could not be present.Therefore, we can search for TiO by applying the cross-correlation technique in transmission spectroscopy.
We combined the five observations to search for the presence of TiO at the terminator of WASP-33b.The crosscorrelation results of TiO are shown in the first row of Fig-ure 6.There are no planetary features as obvious as Ca ii in the CCF map of TiO, and an S/N peak is not located at the expected value (K p = 231 km s −1 , ∆ = 0 km s −1 ).We conclude that the peak is likely due to noise and that there is no evidence for the presence of TiO at the terminator of WASP-33b in the data we used.
We performed injection-recovery tests for TiO.We injected models with different volume mixing ratios (VMRs) of TiO (see Sect. 4.2.1)into the data, and reduced the data after injection and performed the cross-correlation analyses according to Sect.4.2.2.We injected models with the VMR of TiO from 10 −11 to 10 −5 spaced by 1 dex.We also considered some cloudy scenarios (i.e., P cloud = 10 −2 , 10 −3 , and 10 −4 bar).In general, the S/N > 5 at the expected values in the K p map is used as a threshold for robust detection.We used FastChem (Stock et al. 2018) to calculate the VMRs of TiO from the equilibrium chemical models with the planet's equilibrium temperature and different metallicities.Figure 7 shows the results of the injection-recovery tests for TiO.We recovered the injection signal at the minimum VMR of TiO of 10 −10 assuming cloud free.However, we did not find TiO in the real data, which means that the true value of the VMR of TiO at the terminator of WASP-33b is even lower and is much lower than expected from the chemical equilibrium calculations.Even in the presence of high-altitude clouds, the upper limit derived from the injection-recovery tests is lower than the VMR of TiO expected for WASP-33b with its high metallicity (2−15× solar, Finnerty et al. 2023).Overall, the possibility of TiO not being detected at the terminator owing to the presence of high-altitude clouds and/or the influence of H − (corresponding to P cloud = 10 −2 bar) can be roughly ruled out.However, TiO has been detected in the thermal emission spectrum of WASP-33b (Nugroho et al. 2017;Cont et al. 2021), there should be certain mechanisms to allow its presence on the dayside but absence at the terminator.
The lack or low abundance of TiO at the terminator is still an unsolved puzzle.Several mechanisms have been proposed to explain the absence of TiO at the terminator.The colder temperature on the nightside could cause condensation or settling of TiO in the deeper atmosphere, while the strong daynight winds on ultra-hot Jupiters may prevent efficient transport of TiO to the dayside (Hubeny et al. 2003;Spiegel et al. 2009;Parmentier et al. 2013Parmentier et al. , 2016)).TiO could be strongly photo-dissociated due to the high UV flux of the host star, especially during the active cycle of the star (Knutson et al. 2010).TiO could also be strongly thermally dissociated on the dayside for ultra-hot Jupiters, and the detection of TiO at the terminator would depend on its recombination timescale relative to the atmospheric circulation timescale (Lothringer et al. 2018;Parmentier et al. 2018).The detection of TiO on the dayside of WASP-33b and the absence of TiO at the terminator in our results seems to indicate that TiO is condensed at the terminator (and the nightside), but is transported to the dayside by the atmospheric circulation and converted to the gas phase for detection.Further theoretical modelling and observations are needed to confirm this.

5.2.3.
Results of Ti i, V i, and other atomic and molecular species Cont et al. (2022b) reported the first detection of the emission signature of Ti i and V i in a dayside of an exoplanet atmosphere and found no residual stellar pulsation features present in the CCF maps for Ti i and V i in WASP-33b.The absence of stellar pulsation is conducive to the detection of Ti i and V i.In our case, we tried to combine five observations to detect Ti i and V i in the transmission spectrum.The results are shown in the second and third rows of Figure 6.
Although the K p map of Ti i shows a suspicious signal (K p = 247 +7 −34 km s −1 , ∆ = 3 +2 −2 km s −1 ) at an S/N of 6.9, which is close to the expected value, there is no significant cross-correlation signal originating from the planetary atmosphere during the transit in the CCF map of Ti i.The suspicious signal could also be due to a weak stellar pulsation pattern in the CCF map.Therefore, we are cautious in concluding that Ti i is present at the terminator of WASP-33b.For V i,there is no evidence for its presence at the terminator of WASP-33b from the data we used, in agreement with the results of Borsa et al. (2021a), who found no evidence for V i using other HARPS-N observations.A large number of Ti i and V i lines are covered in the wavelength range by the data we used.Thus, the nondetection of Ti i and V i is unlikely to be caused by a low number of spectral lines.The significant detection of Ti i and V i on the dayside, but not on the terminator, could be caused by the difference in the atmospheric temperature between the dayside and the terminator.Moreover, some observations suggest that temperature is not the only factor, as other physical parameters also influence the presence of refractory species such as Ti and V (Cont et al. 2022b).
For completeness, we also searched for other atomic and molecular species using the cross-correlation technique, including Fe i, Fe ii, Ca i, Mg i, Cr i, Y i, VO, AlO, and FeH.The fourth row of Figure 6 illustrates the CCF results of Fe i as a typical example of stellar pulsation.The CCF results for the other atomic and molecular species are shown in Figure 8. Almost all CCF maps of the atomic species are dominated by the stellar pulsation patterns and the RM effect.In the stellar rest frame, the affected velocity range is determined by the rotational velocity of the host star, about ±87 km s −1 for WASP-33.This range just overlaps with the radial velocity of the planet during transit.The CCF maps show that although the pulsation is different from the orbital velocity semi-amplitude of the planetary motion, the strong crosscorrelation signal of the pulsation can completely obscure the planetary signal.As a result, the K p maps are strongly dominated by stellar pulsations.
In contrast, the CCF maps of the molecular species are free of the stellar pulsation patterns and the RM effect because they are not present in the high-temperature atmosphere of the host star.As for TiO, we find no evidence for VO, FeH, and AlO.However, a spurious signal (K p = 208 +14 −15 km s −1 , ∆ = −4 +3 −2 km s −1 ) at an S/N of 2.89 in the K p map of AlO cannot be ruled out.Further high-resolution observations are required to confirm the inference of AlO based on the lowresolution transmission spectrum (von Essen et al. 2019).

CONCLUSIONS
We have observed a single transit of the ultra-hot Jupiter WASP-33b with the NEID spectrograph.With the newly acquired NEID data, we confirmed the previous detections of the Balmer series of hydrogen (Hα, Hβ, and Hγ) and the Ca ii infrared triplet using the line-by-line technique.We introduced the Gaussian process method to account for the correlated systematics caused by stellar pulsations when fitting the line profiles of the Balmer lines and the Ca ii triplet.As a result, we robustly derived a line contrast of 1.24 +0.20  −0.19 % for Hα, 0.62 +0.11  −0.12 % for Hβ, 0.46 +0.15 −0.16 % for Hγ, 0.40 +0.10 −0.10 % for Ca ii λ8498, 0.78 +0.15  −0.15 % for Ca ii λ8542, 0.81 +0.12 −0.12 % for Ca ii λ8662 from the NEID data.We did not find a significant net Doppler shift traced by the Balmer lines or the Ca ii infrared triplet, indicating no evidence for a day-to night-side wind in WASP-33b.
Yan & Henning (2018) detected 1.15±0.05% of the Hα absorption in the atmosphere of KELT-9b, which orbits an A0-type star.The absorption corresponds to a hydrogen atmosphere about 1.64 times the radius of the planet, close to the Roche lobe (1.91 +0.22  −0.26 R p ). Casasayas- Barris et al. (2018Barris et al. ( , 2019) ) claimed that the detection of Balmer series of H and Ca ii in the ultra-hot Jupiter KELT-20b, orbiting an A2-type star, implies an extended atmosphere produced by irradiation far from the Roche lobe of this planet.Fossati et al. (2023) argued that KELT-9b and KELT-20b host hydrostatic atmospheres, and the Balmer series of H and Ca ii lines do not probe the exosphere.This is because their host stars are early A-type and do not have strong X-ray and EUV emission to turn a hydrostatic atmosphere into hydrodynamics.Given an A5-type for WASP-33, the upper atmosphere of WASP-33b could be strongly heated to the hydrodynamic regime (Fossati et al. 2018).In addition, Fossati et al. (2023) shows that it is necessary to consider NLTE effects in the atmospheric modelling of the Balmer series of H and Ca ii.Overall, the absorption of Ca ii and/or Balmer series of H can be linked to the state of the planet's upper atmosphere, but is highly model dependent and requires consideration of hydrodynamic and/or NLTE effects.More detailed atmospheric modelling and comparative studies between planets with different spectral types of host stars are needed.
Using the cross-correlation technique and combining the NEID data with the archival data (HARPS-N and CARMENES), we further confirmed that the Ca ii absorption is from the planet, since it is consistent with the orbital motion of the planet.However, we found no evidence for the species expected to be present at the terminator, such as TiO, Ti i, and V i.These have already been detected in the dayside atmosphere and are not affected by stellar pulsation.The discrepancy between the species inferred from the dayside and terminator atmospheres can be attributed to factors such as differences in abundance, temperature, data reduction procedure, and systematics, which require further highquality observations to reconcile.In addition, we looked for other atomic species and found them difficult to detect in the transmission spectrum because they are strongly affected by stellar pulsation.
So far, a number of ultra-hot Jupiters have been found.Due to the extreme temperature difference between the dayside and nightside atmospheres, the chemical composition and physical properties of day and night are expected to be significantly different (Parmentier et al. 2018), requiring more phase-resolved, high-precision observations.For WASP-33b, its dayside has been well studied, but its terminator has been poorly studied due to the impact of stellar pul-sation.A detailed modeling of the stellar pulsation may be necessary for future transmission spectroscopy of WASP-33b to robustly recover the signals of planetary origin.

Figure 1 .
Figure 1.Transmission spectral matrix, CLV+RM model, and onedimensional transmission spectrum for the Hα line observed with NEID.First row: the transmission spectral matrix before correction for the CLV and RM effects in the stellar rest frame.The horizontal black dashed lines indicate the first and fourth contacts of the transit (i.e., the beginning and end of the transit).The slanted green dashed lines mark the radial velocity shift due to the orbital motion of the planet.The radial velocity regions around −158 km s −1 (i.e.∼ 6559.35Å) are masked due to invalid values (NaN).Second row: the model for the CLV and RM effects.Third row: the observed transmission spectral matrix after correction for the CLV and RM effects.The color bar is indexed to the residual flux.Fourth row: phase-folded one-dimensional transmission spectrum.The CLV and RM effects have been corrected.The y-axis shows the relative residual flux minus one in %.The gray line shows the unbinned observed transmission spectrum.The black circles with error bars show the spectrum binned every ten points (∼0.2 Å).The blue line shows the best-fit GP model.Fifth row: transmission spectrum after removal of systematic noise using the GP method.The red line shows the best-fit Gaussian line profile.The purple and cyan lines are hydrostatic models of Hα with temperatures of 2,710 K and 10,000 K, respectively, including planetary rotation.

Figure 2 .
Figure 2. Transmission spectral matrix of the Hα absorption in the planetary rest frame.The horizontal white lines indicate the first and fourth contacts of the transit.The vertical black dashed lines represent the values of the planetary radial velocity at the first and fourth contacts of the transit.The slanted cyan dashed lines show the radial velocity trace of the planet moved to bracket the pre-transit absorption feature.The color bar indicates the value of residual flux.

Figure 3 .
Figure 3.The transmission spectral matrices and one-dimensional transmission spectra for the Hβ and Hγ lines observed by NEID.Each column corresponds to an individual line.The first row to the fifth row are similar to the forms of Figure 1. three near-infrared triplet (IRT) lines with the NEID data in this section.Similar to Sect.5.1.2,the transmission spectral matrix and one-dimensional transmission spectrum for each line of the Ca ii IRT (8498.018Å, 8542.089Å, and 8662.140Å at air wavelength, labeled as Ca ii IRT λ8498, Ca ii IRT λ8542, and Ca ii IRT λ8662, respectively) are shown in Figure 4 from left to right.Table3shows the best-fit parameters of the Gaussian function based on the classical method and the GP method, respectively.The mean velocity offset of the Ca ii lines is 0.51+1.36−1.32 km s −1 for the GP method.This result is consistent with the Balmer lines when stellar pulsations are treated as correlated noise by the GP, indicating that no strong day-to night-side winds are observed.The mean FWHM of the Ca ii IRT lines (29.7 +3.6 −3.3 km s −1 ) is consistent with previous results.We translate the line contrast of Ca ii into an effective planetary radius of R eff = 1.143 +0.033 −0.034 R p for Ca ii IRT (λ8498), R eff = 1.264 +0.045 −0.046 R p for Ca ii IRT (λ8542), and R eff = 1.273 +0.036 −0.037 R p for Ca ii IRT (λ8662), respectively.In

Figure 4 .
Figure 4.As in Figure 3, but for the Ca ii IRT lines, including Ca ii IRT (λ8498), Ca ii IRT (λ8542), and Ca ii IRT (λ8662).The name of each line is given in the title of each top panel.
+3 −3 km s −1 , taken from Yan et al. 2019) calculated using Kepler's third law.Yan et al. (2019) reported a similar K p (= 224 km s −1 ) value using the Ca ii lines at the terminator of WASP-33b.Cont et al. (2022b) obtained a K p of 225 +3.0

Figure 5 .
Figure 5. Cross-correlation functions of Ca ii averaged over the five nights of observation.First row: the CCF map of Ca ii after correction for CLV and RM effects.The horizontal blue dashed lines represent the first and fourth contacts of the transit.The purple dashed line shows the orbital motion of the planet.Second row: K p map.The magenta pentagram marks the S/N peak.The expected values are indicated by the dashed lines.The color bar represents the value of the S/N.Third row: the one-dimensional cross-correlation function (the pink dashed line) corresponding to the K p value at S/N peak.The vertical gray dashed line indicates ∆ = 0.

Figure 6 .Figure 7 .
Figure 6.Cross-correlation functions of TiO, Ti i, V i, and Fe i without correction for CLV and RM effects.Each row of panels represents one species.The format of each panel from left to right corresponds to each row from top to bottom of Figure 5, respectively.The name of each species is given in each top-left panel.

Table 4 .
Adopted priors for the free parameters in the GP method.the theoretical wavelength of each individual line.

Table 1 .
Summary of the observation logs.

Table 2 .
Parameters of the WASP-33 system.

Table 3 .
The best-fit parameters from the classical and GP methods fit to the one-dimensional transmission spectra.