Transit Spectroscopy of K2-33b with Subaru/IRD: Spin-Orbit Alignment and Tentative Atmospheric Helium

Exoplanets in their infancy are ideal targets to probe the formation and evolution history of planetary systems, including the planet migration and atmospheric evolution and dissipation. In this paper, we present spectroscopic observations and analyses of two planetary transits of K2-33b, which is known to be one of the youngest transiting planets (age $\approx 8-11$ Myr) around a pre-main-sequence M-type star. Analysing K2-33's near-infrared spectra obtained by the IRD instrument on Subaru, we investigate the spin-orbit angle and transit-induced excess absorption for K2-33b. We attempt both classical modelling of the Rossiter-McLaughlin (RM) effect and Doppler-shadow analyses for the measurements of the projected stellar obliquity, finding a low angle of $\lambda=-6_{-58}^{+61}$ deg (for RM analysis) and $\lambda=-10_{-24}^{+22}$ deg (for Doppler-shadow analysis). In the modelling of the RM effect, we allow the planet-to-star radius ratio to float freely to take into account the possible smaller radius in the near infrared, but the constraint we obtain ($R_p/R_s=0.037_{-0.017}^{+0.013}$) is inconclusive due to the low radial-velocity precision. Comparison spectra of K2-33 of the 1083 nm triplet of metastable ortho-He I obtained in and out of the 2021 transit reveal excess absorption that could be due to an escaping He-rich atmosphere. Under certain conditions on planet mass and stellar XUV emission, the implied escape rate is sufficient to remove an Earth-mass H/He in $\sim$1 Gyr, transforming this object from a Neptune to a super-Earth.


INTRODUCTION
Planets in their infancy are invaluable to the study of planet formation and evolution.In particular, Neptune-to Jupiter-size planets younger than ∼ 100 million years (Myr) can provide important clues to their evolution, as pronounced changes to their orbits (e.g., scattering) and atmospheres have been predicted to take place on this timescale (Winn & Fabrycky 2015;Owen 2019).Since the K2 mission (Howell et al. 2014), dozens of young transiting exoplanets have been identified, some of which are amenable to orbital and atmospheric characterization by transit spectroscopy (e.g., David et al. 2019;Newton et al. 2019;Plavchan et al. 2020;Thao et al. 2023).
Specifically, measurement of projected stellar obliquity with respect to the orbit plane (denoted by ) tests scenarios for the evolution of planets onto close-in orbits, e.g.migration within a primordial disk ★ E-mail: hd17156b@gmail.com(TH) (producing a low obliquity) vs. planet-planet scattering or a secular resonance such as Kozai-Lidov followed by tidal circularization (producing a high obliquity).Measurement of obliquities in young systems is especially valuable since the realignment of the stellar spin and planetary orbit is expected to occur on a Gyr timescale (e.g., Winn et al. 2010;Albrecht et al. 2012), after which the different scenarios are indistinguishable.Moreover, while migration within a disk must obviously have been completed while the disk was extant (≲10 Myr), both planet-planet scattering and the Kozai-Lidov mechanism can operate after the disk is gone, on timescales of ≳ 1 − 10 Myr (e.g., Fabrycky & Tremaine 2007;Nagasawa & Ida 2011;Martin et al. 2016).
The technique of transmission spectroscopy can probe the atmospheres of transiting planets if the atmospheric scale height is sufficiently large, e.g. by a high planet equilibrium temperatures and/or low molecular weight atmosphere, and especially if the atmospheres are escaping in a wind.Exceptionally young, close-in planets that are heavily irradiated by their host stars may have extended, escaping H/He-rich atmospheres.Atmospheric escape driven by stellar highenergy irradiation over ∼ 100 Myr is thought to be responsible for the formation of the gap between the radius distributions of super-Earths and sub-Neptunes (Owen & Wu 2017;Fulton et al. 2017) as well as the "Neptune" desert close to stars, although there are other mechanisms proposed to explain those observed properties of exoplanets including the "core-powered mass-loss" (e.g., Ginzburg et al. 2018;Gupta et al. 2022).While the allowed transitions of H I (e.g., Lyman and Balmer series) are either at vacuum wavelengths not accessible from the ground or are not excited at typical wind temperature, the strong 1083.3nm triplet of metastable He I from its ortho-electronic state has been used to detect or constrain the extended atmospheres or winds of several young gas-rich exoplanets (Vissapragada et al. 2021;Gaidos et al. 2022;Zhang et al. 2022bZhang et al. ,a, 2023;;Gaidos et al. 2023).
K2-33b, super-Neptune-sized planet in a 5.42−day orbit, is the youngest transiting exoplanet known to date (Mann et al. 2016;David et al. 2016).Its host is a pre-main-sequence low-mass star, residing in the 8-11 Myr Upper-Scorpius star forming region (Mann et al. 2016).Despite the star's youth, there is no evidence for a circumstellar disk.One aspect of interest is the planet's apparent radius: the planet size was measured as   = 5.04 +0.34  −0.37  ⊕ based on transit observations with K2 (≈ 0.6 m) and the ground-based MEarth telescope array (0.84m) (Mann et al. 2016), but a recent study using transit observations at 1.1-1.7 m with the Wide Field Planetary Camera (WFPC3) on the Hubble Space Telescope and 3.6 and 4.5 m with the IRAC camera on the Spitzer telescope indicate a 30% smaller radius than derived from transit observations optical wavelengths (Thao et al. 2023), possibly due to photochemical haze in an extended atmosphere, or star spots.Ohno et al. (2022) proposed an alternative explanation involving wavelength-dependent scattering in a circumplanetary dust ring, a scenario which can better explain the transmission spectrum of the planet.Further spectroscopic characterisation of K2-33b would help us solve the puzzle of the radius discrepancy, as well as better understand the origin and evolution of this infant planet.
Here, we present high-resolution near-infrared spectroscopy of transits of K2-33b using the InfraRed Doppler (IRD) instrument (Tamura et al. 2012;Kotani et al. 2018) on the Subaru telescope, with the aim of investigating the projected stellar obliquity and the extended helium atmosphere of the planet.The large radius of K2-33b for its insolation level (Mann et al. 2016;Hirano et al. 2018) makes the planet a particularly appropriate target for such a study, and IRD has also been successfully used to search for extended atmospheres around other young planets (Gaidos et al. 2020a,b;Hirano et al. 2020c;Gaidos et al. 2022).

IRD SPECTROSCOPY FOR K2-33
We obtained spectra of K2-33 during transits of "b" using IRD, a fiber-fed spectrograph developed for precise RV measurements in the near infrared (970 − 1730 nm) with the spectral resolution of  ≈ 70, 000 (Kotani et al. 2018).During the science exposures, we also injected light from the laser-frequency comb (LFC) into IRD as a simultaneous wavelength reference.Typical integration time was 900-1500 sec for each frame.We obtained 17 and 19 spectra on UT 2020 June 8 and UT 2021 April 4, respectively, covering the full transit but only limited out-of-transit baseline on both nights.Since the first two frames on UT 2020 June 8 were severely affected by poor weather (i.e., counts less than 100 e − at 1000 nm), we discarded them from the subsequent analyses.
We reduced the raw IRD frames with IRAF (Tody 1993) and our custom code, and extracted wavelength-calibrated, one-dimensional spectra for the star and LFC separately.K2-33's reduced spectra typically had a signal-to-noise (S/N) ratio of 30-45 per pixel at 1000 nm.We analysed those reduced spectra and computed relative RVs for individual frames using IRD's standard pipeline for precise RV measurements; we generated a telluric-free template spectrum for K2-33 using multiple observed spectra, to which the relative RVs were measured by the forward modelling technique (Hirano et al. 2020a).Instrumental drifts were corrected by taking into account the variations in the instantaneous instrumental profile of the spectrograph.For each frame, RVs were calculated for spectral segments each covering ≈ 2 nm, and we derived the final RV for each frame by taking the weighted mean of the segment RVs.Those relative RVs and their errors are summarized in Table 1.

MEASUREMENTS OF THE STELLAR OBLIQUITY
In this section, we present analyses of IRD data to measure the stellar obliquity.To model the spectroscopic transits of K2-33b, we first derived the central transit times (  ) of our transit observations by IRD.Since an updated ephemeris was not given in the follow-up paper by  Thao et al. (2023), we combined the ephemeris in the discovery paper (Mann et al. 2016) with the transit times by Spitzer observations (  = 2457701.57374±0.00091;A. Mann, private communications) and obtained an updated period of  = 5.424871 ± 0.000010 days for K2-33b.Using this revised period, we estimated the predicted transit times of IRD observations as   = 2459008.9676± 0.0026 and 2459329.0350± 0.0032 for the 2020 and 2021 transits, respectively.When deriving the system parameters below, we take into account this uncertainty in   .We note that no evidence of TTVs for K2-33b was reported in any of the past photometric data sets (A.Mann, private communications).

Modelling of the RM Effect
The RM effect has been routinely used to measure the projected stellar obliquity , which is manifested as an anomalous RV variation during transits (e.g., Queloz et al. 2000;Winn et al. 2005;Ohta et al. 2009).K2-33b is a challenging target for this type of observations due to its shallow transit and faintness of the host star ( ≈ 15.7,  = 11.1).Besides, K2-33b is reported to exhibit a chromatic variation of the transit depth (or   /  ), possibly originating from a geometrically thick, hazy atmosphere (Thao et al. 2023) or planetary ring (Ohno et al. 2022).The uncertainty in   /  leads to a systematic error in  when it is combined with the uncertainty in  sin  as well as the impact parameter .
In order to break the degeneracy between  sin  and   /  and better model the RV anomaly due to the RM effect, we first estimated  sin  from the IRD spectrum.Following the procedure in Hirano et al. (2020c), we computed the mean cross-correlation function (CCF) between an observed IRD spectrum of K2-33 and a template spectrum of a similar spectral type (i.e., GJ 436, an M3 star).CCF was computed for each spectral order of an IRD spectrum, and they were averaged to obtain the mean CCF.The black solid line in Figure 1 indicates the resulting CCF for K2-33 (on UT 2021 April 24).To compare it with model CCFs, we generated a number of mock IRD spectra with differing  sin  by convolving the observed IRD spectrum of GJ 436 ( sin  < 0.5 km, Bourrier et al. 2018) with the rotation plus macroturbulence broadening kernel (Hirano et al. 2011), in which we assumed the macroburbulent velocity of  = 1 km s −1 .We then computed the CCF for each mock spectrum as for the observed spectrum of K2-33.Using this set of mock CCFs with various  sin  values, the observed CCF was fitted by interpolation to estimate  sin .We obtained  sin  = 8.62 ± 0.04 km s −1 , and the best-fit CCF model is overplotted in Figure 1 by the red dashed line.The uncertainty in  sin  (0.04 km s −1 ) is small, which only accounts for the statistical error in the fit and is likely underestimated principally due to uncertainty in the macroturbulence velocity.In what follows, we introduce a systematic error of 0.4 km s −1 for  sin  based on the difference between our  sin  value and that derived in Mann et al. (2016).
We analysed the observed RVs from the two transit nights.Since the overall slopes of the RV baseline differ significantly between the two transit nights (Figure 2), we adopted two different baseline parameters, introduced as an instantaneous Keplerian RV semiamplitude  and RV offset parameter  for each night, which empirically represent the RV baselines on the transit nights.
We employed the RM velocity anomaly model by Hirano et al. (2011), which was derived in cases that observed RVs are extracted by the forward-modelling (template-matching) technique as for IRD's  (Mann et al. 2016).The last prior on  14 was introduced in place of a prior on , which was reported to have highly asymmetric errors (i.e.,  = 0.16 +0.19  −0.11 ).Implementing the MCMC analysis (Hirano et al. 2016) in which the two transit data sets are simultaneously modeled and fitted, we derived the system parameters for K2-33b.The result of this analysis is presented in Table 2.The stellar obliquity  was found to be −6 +61 −58 degrees (the 68 % confidence interval), which is compatible with spinorbit alignment, but alignment/misalignment is inconclusive from this RM analysis due to the large uncertainty.The low   /  value of 0.037 +0.013  −0.017 is more consistent with the HST result (  /  ≈ 0.0367 for the HST white light curve at 1.088−1.68m) than those reported by optical transit observations (  /  ≈ 0.04735 for K2 and 0.0489 for MEarth, Thao et al. 2023), although its large uncertainty makes it inconclusive, too.
To check for the dependency of the derived parameters ( in particular) on the choice of prior distributions that we imposed in the fit, we repeated the MCMC analyses after relaxing or removing the Gaussian prior for  sin  or   .We found that while changing the Gaussian width of the prior to three times the original width did not significantly change  for both  sin  and   , completely removing the Gaussian prior (i.e., imposing a uniform prior) on   led to a very poor constraint on the obliquity ( = 3 +115 −124 degrees).Assuming a uniform prior for  sin  resulted in an apparently consistent result ( = −12 +70 −58 degrees), but we found a strong degeneracy in the fitting parameters ( sin  and   /  in particular), and almost no meaningful constraints were obtained for  sin  and   /  .

Doppler-Shadow Analysis
Unfortunately, our measurements of the RM effect did not provide a good constraint on the stellar obliquity, likely owing to the large RV uncertainties and degeneracy in the system parameters; the transit impact parameter  of K2-33b is relatively low, in which case the classical RM modeling is known to severely suffer from the param-eter degeneracy (e.g., Albrecht et al. 2011).We thus attempted a secondary analysis to constrain the obliquity based on the Dopplershadow (or more often referred as "Doppler tomography") technique (e.g., Collier Cameron et al. 2010;Johnson et al. 2017).
The Doppler-shadow technique has already been applied to IRD data to measure  (e.g., Hirano et al. 2020b,c;Gaidos et al. 2020bGaidos et al. , 2022)), and we follow the same approach here.Briefly, we computed CCFs for individual frames using GJ 436's spectrum as a CCF template as described in Section 3.1, and derived the mean out-of-transit CCF after correcting for the barycentric motion of Earth.Since K2-33 exhibits very high stellar activity and two transits (2020 and 2021) were observed by IRD, we combined nightly frames on each transit night individually to focus on the relative line-profile variations within the night.We then subtracted the mean CCF from the individual CCFs for each night to extract "residual" CCF variations against time1 .For the 2020 data set, we discarded the first frame from the analysis due to the very low S/N ratio of the original spectrum.The top panels of Figure 3 display those residual CCF maps for the 2020 (panel (a)) and 2021 (panel (b)) transits.The three vertical dotted lines represent the CCF line-profile center (middle) and its ±8.6 km s −1 , which approximately corresponds to  sin  of K2-33.
The residual CCF maps are noisy and the planet shadow is hardly recognisable.We thus opted to combine the two transits to mitigate the impact of low S/N ratios for each transit.Panel (c) of Figure 3 presents the combined (by weighted mean) residual CCFs using both 2020 and 2021 data sets.The map still suffers from systematic patterns, but it suggests a planet shadow around the blueshifted edge of the line profile ( ≈ −15 km s −1 ) at ingress and that moves redward until the egress ( ≈ 2.5 km s −1 ).
To model the observed CCF residuals and estimate the system parameters, we synthesized mock IRD spectra that mimic in-transit spectra during a transit of K2-33b.The details of the model calculation are given in Hirano et al. (2020c); we created a total of 2142 mock spectra, changing  sin  (from 4.0 to 12.0 km s −1 ) and the twodimensional position of the planet on the stellar disc.Given the good agreement between our estimate of   /  from the classical RM modelling (Section 3.1) and that obtained by HST observations in the near-infrared (1.088 − 1.68 m), we fixed   /  at 0.367 (Thao et al. 2023) generating mock transit spectra.All mock IRD spectra were then subject to the CCF analysis using GJ 436 as a template.With all these CCFs, one can produce a planet shadow model in the residual CCF map for a given set of system parameters ( sin , , /  , etc) by interpolations.
Following Hirano et al. (2020c), we performed an MCMC analysis to estimate the obliquity  as well as the other system parameters for K2-33b.In doing so, the time sequence of the combined residual CCF (panel (c) of Figure 3) was fitted by the above CCF model, and we derived the posterior distributions for the fitting parameters:  sin , , /  , , and   .As in Section 3.1, we imposed Gaussian priors on  sin , /  ,  14 , and   .The other transit parameters, in particular   /  , ,  1 , and  2 , were held fixed to the literature and theoretical values.The result of the fit is also given in Table 2, and the residual CCF model with the best-fit parameters is depicted in panel (d) of Figure 3.The spin-orbit angle  is now constrained with a better precision than that in the classical RM modelling, likely reflecting the discernible signal of the planet shadow.Its estimation of  = −10 +22 −24 degrees implies that K2-33b's orbit is almost aligned with the equatorial plane of the host star.

HE I TRANSMISSION SPECTROSCOPY
We constructed summed, normalised spectra of K2-33 inside and outside of each of the two transits of "b" , and calculated difference spectra that should contain any absorption by the planet's extended atmosphere/wind (Fig. 4).During both events the He I triplet is affected by the weaker, bluer of two neighboring OH lines (vertical red lines in upper panels of Fig. 4).The quality of the difference spectrum from the 2020 transit (lower panel of Fig. 4a) is poor, likely due to the low signal-to-noise of the available out-of-transit spectra, and there is no indication of excess absorption.However, we identify possible absorption during the 2021 transit (lower panel of Fig. 4b).
To quantify any transit-associated absorption and evaluate its sig-nificance, we fit a two-parameter model He I triplet line profile to unaffected regions of the spectrum (black points in the lower panels of Fig. 4).The profile has a variable amplitude (expressed as an equivalent width, EW) and is broadened by the instrument resolution and gas temperature (10 4 K), plus variable Doppler broadening that represents the collective motion of the planet and its wind.We determined 95% confidence ratios by fitting 1000 Monte Carlo representations of the actual data with added random Gaussian noise.The best-fit profiles for the 2020 and 2021 transits are shown as the heavy blue curves in the lower panels of Fig. 4 and have EWs of -29 mÅ and +63 mÅ, respectively.The 95% confidence range are plotted as the light blue shaded areas in Fig. 4 and are from -63 to +8 mÅ and and 55-71 mÅ, respectively.Thus while there is no significant feature in the 2020 transit, the anomalous absorption during the 2021 transit appears statistically significant.However, these calculations do not include systematic (time-correlated) and "red" (wavelength correlated) noise and thus the significance of the 2021 detection should be interpreted with caution.We further searched for anomalous absorption by performing fits of the same profile to each individual spectrum differenced by the mean out-of-transit spectrum, in these cases fixing the Doppler broadening term to the best-fit value (39 km s −1 ) obtained from the fit to the mean spectra of the 2021 transit.Uncertainties in each EW are calculated using the same Monte Carlo scheme describe above.The He I EW time series (Fig. 5) show no significant variation during the 2020 transit ( 2 = 31.0for 13 d.o.f.) but contains a clear transitlike trend corresponding to the predicted transit interval during the 2021 transit ( 2 = 160.1 for 17 d.o.f.).
The 2021 He I time-series contains a large positive excursion near transit center.We speculate that this is the result of a flare on this active star: strong flare-associated He I emission has been observed on other active M dwarfs (Schmidt et al. 2012;Hintz et al. 2020;Kanodia et al. 2022;Howard et al. 2023) and the decay timescale (∼ 30 min) is characteristic of flares.Conclusive evidence for a flare comes from an analysis of the 1.2822 m Paschen  line of H I in the same spectra: this reveals simultaneous, transient emission superposed on the photospheric absorption line (Figs.6 and 7).The 11-12 km sec −1 red-shift is too large to be produced by stellar rotation ( eq ≈ 8 km s −1 ) but might be produced by electron beaming (Kowalski et al. 2022).Even more curious: in contrast to observations of emission in the Paschen lines during other flares (Kanodia et al. 2022;Fuhrmeister et al. 2023, e.g.,), this line is narrow and only marginally resolved (FWHM ≈ 0.03 nm).We re-fit the in-transit spectrum excluding the two epochs most affected by the flare; the best-fit EW was unchanged (63 mÅ), but the 95% confidence range shrank to 58-68 mÅ.
We used the best-fit and upper limit values of the EW from the 2021 transit to constrain possible mass loss from the planet, based on a model of an isothermal, spherically symmetric Parker wind with a solar H/He composition as described in Gaidos et al. (2020a).A key input to the model is the X-ray/EUV emission from the host star, which is responsible for formation of the triplet He I state via ionization of He I, but for which essentially nothing is known observationally.In lieu of data, we adopted the panchromatic spectrum of GJ 729, an active, rapidly-rotating (2.85 day period) M3.5 dwarf (compared to the M3 spectral type of K2-33) from the Mega-MUSCLES survey planet-hosting K and M dwarfs (Wilson et al. 2021).2Fluxes in specific X-ray and UV bands were adjusted using the empirical relations between these bands and Lyman- emission developed by Linsky et al. (2014).The Lyman- emission was in turn estimated from the star's X-ray luminosity and the empirical relation of Linsky et al. (2013).In turn the X-ray luminosity was estimated using the star's Rossby number and the relation of Wright et al. (2018).(K2-33 does not have a corresponding source in the Second ROSAT All Sky Survey source catalog (Boller et al. 2016), but the flux at Earth based on the estimated X-ray luminosity is only 2 × 10 −13 ergs s −1 cm −2 , comparable to the detection limit of the survey.)The Rossby number /  was calculated to be 0.10 (in the saturated regime) using the rotation period  from Klein & Donati (2020) and a convective turnover time   of 65 days calculated using the stellar luminosity of Thao et al. (2023) and scaling by the squareroot of the luminosity (Pizzolato et al. 2003) relative to a solar value of 25 days.We assumed two values for the planet mass (which sets the atmospheric scale height and gravitational potential well but is a poorly constrained parameter of the wind model) of 25 and 5 M ⊕ .The first is motivated by the empirical relation between mass and radius of   ∼  2  , while the second invokes an "extreme" scenario of a super-Earth rocky core surrounded by about an Earth-mass of H and He.
In Fig. 8, calculated equivalent widths in mÅ are plotted as a function of mass-loss rate (between 0.0003-0.3M ⊕ Myr −1 , bracketing the estimated energy-limited escape rates, see below) and wind temperature (3000-30000K).The blue regions are the 95% confidence range for anomalous absorption in the 2021 transit (Fig. 4b).A lower planetary mass permits higher escape rates, but for a given escape rate the He I column density along the line of sight is lower due to slower wind speeds.The total incident X-ray plus EUV and Ly- flux (averaged over the planet) is predicted to be 1.5 × 10 4 ergs s −1 cm −2 , and the corresponding XUV energy-limited escape rates (for 100% efficiency) are 0.0047 and 0.023 M ⊕ Myr −1 for 25 and 5 M ⊕ (vertical red dashed lines in Fig. 8).For the 5 M ⊕ case, the available XUV emission is insufficient to induce a wind of sufficient density to produce the observed anomaly. 3In the 25 M ⊕ case, and XUV-driven wind can plausibly produce the anomaly but only if the wind temperature is low and the escape efficiency is high, or the stellar XUV luminosity is significantly higher than predicted.

DISCUSSION
Our measurements suggest a low stellar obliquity for K2-33.Our constraint on  (−10 +22 −24 deg) can be combined with the stellar inclination, which is the obliquity projected onto the line-of-sight, to constrain the unprojected obliquity.The stellar inclination is estimated via  sin  and the equatorial rotation velocity derived from the rotation period and radius of the star.Mann et al. (2016) constrained the stellar inclination angle of K2-33 to   > 63 • at 1 ; here we revised this constraint based on the updated rotation velocity of  sin  = 8.6±0.4 km s −1 .Adopting the rotation period of  rot = 6.29 ± 0.17 days (Mann et al. 2016) and making use of the Bayesian inference scheme described by Masuda & Winn (2020), we computed the posterior distribution of cos   .This yielded cos   = 0.20 +0.17 −0.06 , corresponding to   = 78 +8 −9 deg.Thus, both the sky-projected and line-of-sight stellar obliquities of K2-33 are low.Together with  = −10 +22 −24 degree (Section 3.2) and the orbital inclination of   = 89.1 +0.6  −1.1 deg (Mann et al. 2016) and using the geometric equation for the stellar obliquity (e.g., Eq. ( 9) of Fabrycky & Winn 2009), we further derived the 3-dimensional (unprojected) stellar obliquity as  = 23 +16 −11 deg.The low stellar obliquity of the K2-33 system is in line with the results of other obliquity measurements for very young transiting planets; To this point, stellar obliquities for six transiting systems younger than 100 Myr (AU Mic, HIP 67522, V1298 Tau, DS Tuc A, TOI-942, and K2-33) have been reported (Hirano et al. 2020c;Palle et al. 2020;Martioli et al. 2020;Addison et al. 2021;Heitzmann et al. 2021;Gaidos et al. 2022;Johnson et al. 2022;Zhou et al. 2020;Montet et al. 2020;Wirth et al. 2021), and all measurements indicated low obliquities that are compatible with spin-orbit alignment.In other words, no evidence for a primordial spin-orbit misalignment was found for newborn planets.Since the timescale of obliquity damping by tides is generally longer than ∼ 1 Gyr (Albrecht et al. 2012(Albrecht et al. , 2022)), those planets (< 1 Gyr) were born with a primordially aligned orbit with the stellar equatorial plane, evidence that they formed at their present location or migrated inwards along the original disc plane (e.g., Lin et al. 1996;Ida & Lin 2008;Baruteau et al. 2014) without chaotic orbital disturbance.This is in marked contrast to the distribution of  for older counterparts (≳ 1 Gyr), a significant fraction of which is known to exhibit spin-orbit misalignment (Fig. 9) (also see e.g., Albrecht et al. 2022).However, one should notice that many of these cases are isolated hot Jupiters while young planets with obliquity measurements are generally smaller in size or have longer orbital periods; those measurements are possibly probing different populations of close-in planets.Note that when we extend the age range to 100 − 1000 Myr, we find a fraction of planets having highly misaligned orbits (Fig. 9), but most of those are very hot jovian planets orbiting early-type stars with  eff > 6250 K, a boundary that divides the mostly aligned and misaligned systems (Winn et al. 2010;Albrecht et al. 2012).Further observations for planets with differing ages and sizes would help us unveil the age evolution of the obliquity distribution.
Although our 2021 observations of the K2-33 system show transitassociated excess absorption in the 1.083 m of metastable ortho-He I consistent with an extended, escaping atmosphere, the significance of this detection awaits a better understanding of systematic effects, particularly changes in the stellar chromosphere line on hours-long timescales (see also Gaidos et al. 2022).Although we lack sufficient baseline observations to definitely rule out a stellar source for the variation, the absence of significant variation during the 2020 transit (Fig. 5a) is encouraging.As another assessment, we used other spectra of K2-33 taken outside of transit for RV monitoring (Hirano et al. in prep), to generate residuals between pairs of spectra obtained at least 30 min and at most 5 hours apart, and performed fits of the same He I triplet profile to retrieve spurious EWs.We found that, with the exception of spectra from a single run in 2019 that appears to be affected by systematics, the scatter in best-fit EWs fell below that determined from the 2021 transit.
The lack of a detectable He I anomaly in 2020 then suggests that the planetary wind varies on year-long timescales, e.g., with the changing level of magnetic activity on and XUV irradiation from the star.Magnetic cycling of a few years period has been observed on very cool dwarfs stars, including T Tauri stars (e.g., Klein et al. 2021;Finociety et al. 2023;Lin et al. 2023).Ultimately, there is no substitute for additional transit observations with adequate out-of-transit baseline to resolve the issue.Higher signal-to-noise observations that can avoid or better correct contamination by tellurics, especially the intense OH lines, are also needed.Finally, observations and analysis of K2-33 in X-rays, e.g. with eROSITA and XMM-Newton, are needed to better estimate the high-energy irradiance of its planet. 4If the atmosphere of K2-33b is escaping at a rate of 3 × 10 −3  ⊕ Myr −1 (Fig. 8), after several hundred Myr most or all of a primordial H/He envelope would be lost, leaving this sub-Neptune as a future super-Earth.The reduced radius in the near infrared revealed by Thao et al. (2023) was not confirmed for our RV data due to the large RV uncertainty and/or degeneracy in the fitting parameters.To see what extent we can constrain   /  in case of a perfect spin-orbit alignment, which is supported by the Doppler-shadow analysis (Sec.3.2), we repeated the MCMC fit to the RV data assuming  = 0 degree (while  sin  was allowed to vary with the same Gaussian prior).As a result, we found   /  = 0.035 +0.009 −0.012 , which gives a slightly tighter constraint (i.e., deviated from the K2 value by > 1 ), but the smaller radius in the IRD wavelength range is still inconclusive.
To further discuss the atmospheric escape/evolution of K2-33b as well as the origin of the smaller near-infrared radius, constraining its planet mass would be an urgent task.Near-infrared observations such as by IRD have a marked advantage over visible RV measurements in that activity-induced RV variations are known to be mitigated by virtue of improved spots' contrasts at longer wavelengths (e.g., Crockett et al. 2012;Miyakawa et al. 2021).However, the detection of a large flare in the IRD data during a mere 10 hr of monitoring suggests that such events are frequent, and less energetic events may be even more frequent.Flares are a source of astrophysical noise or "jitter" in RV measurements (e.g., Reiners 2009;Pietrow et al. 2023) and indeed, we see our RV measurements contemporaneous with the flare (see an RV excursion around   in Fig. 2b).Flare may be an impediment to precise infrared RVs for young, active stars such as the host of K2-33b, but their identification using "fingerprint" lines such as Paschen  offers a potential way forward.

CONCLUSIONS
We obtained high-resolution infrared ( -band) spectroscopy of the 8-11 Myr-old low-mass star K2-33 during two transits of its closein 5 ⊕ planet, one of the youngest such objects known.By modelling the RV variations during the transits, we found hints of low projected stellar obliquity for K2-33 ( = −6 +61 −58 deg) and smaller radius of K2-33b in the near-infrared band (  /  = 0.037 +0.013 −0.017 ) than in the optical, but those findings are inconclusive due to the lower RV precision.We also constrained the obliquity by analysing the "shadow" of the planet in cross-correlation profiles with wavelength to better constrain the obliquity ( = −10 +22 −24 deg).Combining the latter constraint with the orbital and stellar inclination angles revealed by the transit and  sin  measurements, we found a low de-projected spin-orbit angle of  = 23 +16 −11 deg for K2-33b.We also compared the mean spectra inside and outside of transits of K2-33b to search for excess absorption in the 1083 nm triplet of metastable ortho-He I.While the 2020 data contain no evidence for anomalous absorption, the 2021 transit shows possible absorption consistent with an extended atmosphere, although we cannot definitely rule out an astrophysical signal from the star.Our modeling of the signal with a spherical model of a planetary wind suggests mass loss of at most a few Earth-mass of solar-composition H/He over a Gyr, potentially driven by XUV irradiation, and dependent on the (as yet determined) planet mass.

Figure 1 .
Figure 1.Normalized CCFs between K2-33's IRD spectrum and the template spectrum of GJ 436.The black solid line represents the observed CCF while the red dashed line corresponds to the best-fit theoretical model with  sin  = 8.62 km s −1 .

Figure 2 .
Figure 2. RV variations by Subaru/IRD for the 2020 (panel (a)) and 2021 (panel (b)) transits.Panel (c) plots the two transit observations after subtracting the RV baseline for each dataset.In panels (a)-(c), the best-fit theoretical models are overplotted by the black solid lines.The RV residuals around the best-fit model are shown in panel (d).

Figure 3 .
Figure 3. Residual CCF contours after subtracting the mean out-of-transit CCFs for the 2020 (panel (a)) and 2021 (panel (b)) datasets.The color contour represents the fractional variation of the residual CCF.The horizontal dashed line show the transit center and times of ingress and egress.The vertical dotted lines correspond to the stellar-line center and approximate line edges (± sin  from the line center).The combined CCF contour using the two datasets is depicted in panel (c), and its best-fit theoretical model is presented in panel (d).

Figure 4 .Figure 5 .
Figure 4. IRD spectra of K2-33 in the vicinity of the 1083.3nm He I triplet obtained during the transit opportunities on 8 June 2020 (left) and 24 April 2021(right).The upper panels show the in-transit (cyan) and out-of-transit (orange) spectra with significant stellar lines (including the He I singlet and unresolved doublet) labeled, as well as telluric OH emission (red lines) and H 2 O absorption (blue inverted triangles).The lower panels contain the difference spectra (black lines) along with the best-fit He I line (blue curves) and the 99% upper limit (magenta lines), with the EWs reported in the legend.

Figure 6 .Figure 7 .
Figure 6.Time-series spectra of K2-33 around the 1.282 m Paschen  line of H I during the 2021 transit of 'b', showing the photospheric absorption line and the red-shifted emission line from a flare.The spectra are temporally ordered from bottom to top and the dashed red line is the best-fit Voigt line profile used as a baseline to compute emission strength.

Figure 8 .
Figure 8. Predicted absorption in the He I from a spherical, isothermal Parker wind model of an escaping H/He atmosphere around K2-33b, assuming two different planet masses (25 and 5 ⊕ ).Contours of EW are in mÅ.The blue region is the absorption inferred from observations of the 2021 transit (Fig. 4b).Note the horizontal axes are different between the left and right panels.Vertical dashed red lines mark the estimated XUV energy-limited escape rate assuming 100% efficiency.

Figure 9 .
Figure 9. Observed projected obliquity || as a function of system age (< 500 Myr).Planets around hot stars ( eff > 6250 K) are plotted by red triangle, while those around cool stars ( eff < 6250 K) are shown in blue.K2-33b (black square) is the youngest planet with a measured obliquity.The right panel shows the histograms of || for systems older than 500 Myr.Those for hot and cool stars are shown in red and blue, respectively.Data downloaded from the NASA Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu/index.html).

Table 2 .
(Thao et al. 2023)nalysisRVs.To take into account the chromatic dependence of the transit depth(Thao et al. 2023), we allowed   /  to float freely, although a single   /  value was fitted in the whole observing band covered by IRD ( , , and ).The fitting parameters in our RM model are ,  sin ,   , the scaled semi-major axis /  , the transit impact parameter , and   /  in addition to the baseline parameters for each transit data set described above.The limb-darkening parameters  1 and  2 for the quadratic law are required in the RM model.
Claret et al. (2013)precision of the individual RV data points, we fixed those parameters at  1 = 0.09 and  2 = 0.35 based on the theoretical values for a  eff = 3500 K, log  = 4.0 star byClaret et al. (2013)for the -band.To ensure the convergence of the fit, we imposed Gaussian priors of  sin  = 8.6 ± 0.4 km s −1 ,   − 2454833 = 4496.0350± 0.0032, /  = 10.29 ± 0.39, and the full transit duration of  14 = 4.08 ± 0.07 hours, based on the above discussion and literature values