Probing the extent of WASP-52 b's atmosphere. High-resolution observations and 3D modeling insights

TL;DR

This study combines high-resolution spectroscopy and 3D hydrodynamic modeling to investigate the atmospheric escape of WASP-52 b, revealing a transitional escape regime with no strong evidence of extended helium absorption or anisotropic outflow.

Contribution

It introduces a comprehensive approach using broad phase high-resolution spectra and advanced 3D modeling to characterize atmospheric escape in WASP-52 b, a regime not extensively studied before.

Findings

No significant helium absorption detected beyond transit.

Model suggests a possible escape structure below detection limits.

WASP-52 b exhibits a transitional atmospheric escape regime.

Abstract

WASP-52 b is an inflated hot Jupiter with a large Roche lobe filling fraction, positioned in the hot Neptune desert. Previous in-transit observations of the helium triplet at 10833 A have reported a range of excess absorption values (1.5%-5.5%) and a lack of net blueshift relative to the planet's rest frame, distinguishing it from other escaping atmospheres. This study investigates the extent and morphology of material escaping from WASP-52 b, assessing whether its outflow resembles a stream-like structure, as suggested for HAT-P-67 b and HAT-P-32 b. We obtained high-resolution spectra with CRIRES+ and CARMENES, covering a broader orbital phase range ($\varphi \approx \pm0.1, \pm0.2, 0.5$) than previous studies. By analyzing the He I 10833 A line as a tracer of escape, we search for extended absorption beyond transit. Additionally, we explore possible outflow morphologies with three-dimensional (3D) hydrodynamic simulations, coupled with an improved radiative transfer approach, assessing the He I 10833 A triplet. The helium line shows no significant evidence of planetary material at the orbital phases observed in this work, though 3D modeling suggests such a structure could exist below observational detection limits. We conclude that the atmospheric outflow of WASP-52 b can be characterized by an intermediate hydrodynamic escape parameter, placing it in a transitional regime between cold outflows forming a stream-like morphology and hot outflows forming a tail. Additionally, the absence of a detectable in-transit blueshift in the helium line rules out a strong day-to-nightside anisotropy scenario.

Tables

Table 1. Table 1: System parameters of WASP-52 b from Hébrard et al. ( 2013 ) .

Parameter		Value
Orbital period	$P_{\mathrm{orb}}$ [d]	$1.7497798\pm 0.0000012$
Transit midpoint	$T_c$ [d]	$2455793.68143\pm 0.00009$
Transit duration	$T_{14}$ [h]	$1.8096\pm 0.0120$
Planet mass	$M_p$ $\left[M_{\mathrm{jup}}\right]$	$0.46\pm 0.02$
Planet radius	$R_p$ $\left[R_{\mathrm{jup}}\right]$	$1.27\pm 0.03$
Eccentricity	$e$	0
Inclination	$i$ [deg]	$85.35\pm 0.20$
Semi-major axis	$a$ [au]	$0.0272\pm 0.0003$
Star mass	$M_{\ast }$ [$M_{☉}$]	$0.87\pm 0.03$
Star radius	$R_{\ast }$ [$R_{☉}$]	$0.79\pm 0.02$
Rotational period	$P_{\mathrm{rot},\ast }$ [d]	$11.8\pm 3.3$
Stellar Age	[Gyr]	${0.4}_{-0.2}^{+0.3\left(a\right)}$

Table 2. Table 2: Overview of observation nights with CRIRES + .

Night	Date	$\overline{\phi }$	# Spectra
1	30/09/2023	$-0.17$	4
		$-0.14$	4
2	08/10/2023	$+0.42$	4
3	13/10/2023	$+0.27$	4

Table 3. Table 3: Overview of observation nights with CARMENES.

Night	Date	$\overline{\phi }$	# Spectra
			VIS	NIR
1	05/09/2023	$+0.50$	9	3
2	08/09/2023	$+0.21$	2	0
3	09/09/2023	$-0.22$	4	4
4	20/09/2023	$+0.05$	4	(4)∗
5	02/10/2023	$-0.11$	4	4

Table 4. Table 4: Input parameters of the 3D simulation models.

model	${\lambda }_p$	$f_{\mathrm{pres}}$	${\rho }_p$	${\rho }_{\ast }$	$c_{s,\mathrm{sub}}$	$c_{s,\mathrm{anti}}$
		[%]	[g cm-3]	[g cm-3]	[km s-1]	[km s-1]
ISO	5	100	$5.42\times {10}^{-16}$	$2.87\times {10}^{-14}$	11.3	11.3
ANISO	5	10	$1.81\times {10}^{-15}$	$2.87\times {10}^{-14}$	11.3	1.1
TORUS	5	100	$5.42\times {10}^{-16}$	$2.87\times {10}^{-15}$	11.3	11.3
TAIL	5	100	$5.42\times {10}^{-16}$	$7.22\times {10}^{-14}$	11.3	11.3

Table 5. Table 5: Density scaling factors and time-averaged mass-loss rates from 3D simulation models.

model	$S$	$⟨{\dot{m}}_p⟩$	$⟨{\dot{m}}_{\ast }⟩$
		[g s-1]	[g s-1]
ISO	1.40	$2.03\times {10}^{11}$	$9.98\times {10}^{11}$
ANISO	0.86	$1.73\times {10}^{11}$	$5.92\times {10}^{11}$
TORUS	1.20	$1.36\times {10}^{11}$	$6.85\times {10}^{10}$
TAIL	1.40	$1.74\times {10}^{11}$	$2.79\times {10}^{12}$

Equations

$$P_p={\rho }_p\frac{G{M}_p}{\gamma {\lambda }_p{R}_p}\text{leavevmode},$$

$$T=\gamma \frac{\mu m_H}{k_{\mathrm{B}}}\frac{P}{\rho },$$

$${\tau }_{\nu }=\int n_s{\sigma }_0\Phi (\Delta \nu )\text{leavevmode}\mathrm{d}s.$$

$$\mathrm{EW}=f\cdot \Delta \lambda =f\cdot \frac{{\lambda }_0}{R},$$

Full text

1 Introduction

WASP-52 b is a transiting exoplanet first identified in 2013 through the Wide Angle Search for Planets (WASP) survey (Hébrard et al., 2013). This inflated hot Jupiter orbits a K-dwarf star of young to middle age (system parameters are listed in Table 1). WASP-52 b is notable for its large size relative to its Roche Lobe, with $R_{p}/R_{\rm RL}\approx 0.58$, suggesting significant tidal interactions with its host star. In terms of surface gravity and Hill radius, the planet is similar to HAT-P-32 b (see Fig. 9 in MacLeod et al., 2024). This makes WASP-52 b an intriguing candidate for investigating atmospheric escape in the regime of extended, stream-like outflows, similar to those observed in planets like HAT-P-67 b and HAT-P-32 b (Gully-Santiago et al., 2024; Zhang et al., 2023a).

WASP-52 b has been the subject of numerous spectroscopic studies aiming to characterize its atmosphere (Kirk et al., 2016; Louden et al., 2017; Mancini et al., 2017; Chen et al., 2017; Alam et al., 2018; Bruno et al., 2018; May et al., 2018; Kirk et al., 2022; Allart et al., 2023; Cannock et al., 2024; Fournier-Tondreau et al., 2024). Evidence of atmospheric escape has been obtained through detections in the helium triplet at 10833 Å. Using ultra-narrowband photometry, Vissapragada et al. (2020) detected helium absorption of $0.29\%\pm 0.13\%$, which was later confirmed by Keck/NIRSPEC observations with a spectral resolution of R = 25 000. These observations revealed an excess absorption of up to 3.4%, with a velocity shift relative to the planet’s rest frame consistent with zero (Kirk et al., 2022).

The lack of a net blueshift is unusual compared to most exoplanets with detected helium signatures, where the line centroid typically shows a net blueshift corresponding to velocities of several km s*-1*(e.g. Nortmann et al., 2018; Spake et al., 2018; Allart et al., 2018; Palle et al., 2020; Kirk et al., 2020; Orell-Miquel et al., 2023; Zhang et al., 2023a; Masson et al., 2024). This blueshift is generally interpreted as evidence of day-to-nightside winds in the planet’s upper atmosphere, driven by pressure gradients between the strongly irradiated dayside and the cooler, shadowed nightside (e.g. Nail et al., 2024). Interestingly, HAT-P-67 b and HAT-P-32 b also show no significant velocity shift during transit, which aligns with the 3D hydrodynamic model predictions of Nail et al. (2025), that a cold stream-forming outflow should have an approximately zero velocity shift. This further motivates our investigation into whether WASP-52 b has a similar outflow structure.

Additionally, Allart et al. (2023) observed the helium triplet using the high-resolution spectropolarimeter SPIRou (R = 70 000) on the Canada–France–Hawaii Telescope (CFHT), placing an upper limit of 1.69% on the in-transit absorption. Furthermore, JWST NIRISS/SOSS observations (R = 700) by Fournier-Tondreau et al. (2024) detected helium excess absorption with a significance of 7.3$\sigma$, corresponding to an estimated absorption of 5.5 $\pm$ 0.9% in high-resolution observations. Additionally, they reported a tentative detection of post-transit absorption at the 2.9$\sigma$ level. Notably, all previously mentioned transit observations only covered orbital phases within a limited range of $-0.06$ to $0.07$ at most.

WASP-52 is an active star (Bruno et al., 2020; Rosich et al., 2020; Salisbury et al., 2021), with a $\log R^{\prime}_{HK}$ index of approximately -4.4 (Hébrard et al., 2013). It shows frequent star-spot crossings during planetary transits (e.g. Mancini et al., 2017; Louden et al., 2017; May et al., 2018; Bruno et al., 2018, 2020; Fournier-Tondreau et al., 2024). Mancini et al. (2017) accounted for star-spot occultation events during transit observations to derive a more precise sky-projected orbital obliquity of $\lambda_{\rm ob}=3.8^{\circ}\pm 8.4^{\circ}$. This value remains broadly consistent with the earlier estimate from the Rossiter-McLaughlin effect, $\lambda_{\rm ob}=24^{\circ}{}^{+17}_{-9}$, reported by Hébrard et al. (2013).

WASP-52 b is located in a unique position within the hot Neptune desert, residing in the high-mass range with a short orbital period (Owen & Lai, 2018; Vissapragada et al., 2020). This location makes it an intriguing case for studying hot Jupiter formation and evolution. Owen & Lai (2018) suggested that WASP-52 b may have reached its current position through a mechanism other than high-eccentricity migration, or that its inflation occurred only after the planet had fully circularized.

WASP-52 b is a compelling candidate for studying the escape of atmospheric material due to its location in the hot Neptunian desert, its significant Roche lobe filling, and the observed absence of a blueshift in helium absorption. We have conducted observations using the CARMENES spectrograph at the Calar Alto Observatory and the CRIRES+spectrograph at the Very Large Telescope (VLT) in the visible (VIS) and near-infrared (NIR). These observations cover a larger phase range of the planet’s orbit than usual to constrain the full extent of escaping planetary material. By combining these observations with 3D hydrodynamic models, we aim to constrain the atmospheric outflow morphology. Specifically, we investigate whether WASP-52 b has a widely extended outflow similar to those observed in HAT-P-67 b and HAT-P-32 b (Gully-Santiago et al., 2024; Zhang et al., 2023a), which may be associated with a stream morphology (Nail et al., 2025; MacLeod et al., 2024). In addition, we introduce a new and improved approach for performing radiative transfer analysis on our 3D simulations. The new method is based on the use of the Cloudy code (Ferland et al., 1998). Cloudy’s extensive database of reaction rates and cross-sections makes this method more appropriate for simulating low-temperature outflows that produce stream-like morphologies, as indicated by our previous work (Nail et al., 2025). Additionally, Cloudy calculations include dozens of atomic and ionic species, which can potentially support investigations of additional tracers of atmospheric escape (e.g. Linssen & Oklopčić, 2023) in the context of large-scale 3D simulations.

This study is organized as follows: In Section 2, we describe our observing strategy, the data obtained for this study, the data reduction process, and we present our observational results from the helium line. In Section 3, we describe the 3D simulation setup with Athena++, the new framework developed for post-processing to generate synthetic spectra from these models, and we present the 3D modeling results and compare them to the observational helium data. This is followed by a discussion in Section 4 and a brief summary in Section 5.

2 He i 10833 Å observations

We obtained high-resolution spectra of WASP-52 in the visible and near-infrared, covering multiple orbital phases outside of transit, to investigate the extent of atmospheric escape of planet WASP-52 b. The metastable helium line He i 10833 Å is a commonly used tracer for atmospheric escape in low-density regions (e.g. Oklopčić & Hirata, 2018; Nortmann et al., 2018; Allart et al., 2023; Ballabio & Owen, 2025), making it an ideal starting point for our analysis. Additional potential spectral tracers of atmospheric escape, particularly those detectable with a high-resolution spectrograph in the visible wavelength range, may also provide valuable insights into the composition and dynamics of the outflow. However, exploring such tracers is beyond the scope of this study; here, we focus exclusively on helium.

Figure 1 illustrates our observing strategy. Unlike traditional transit observations, where a time series of spectra is collected throughout the transit with a baseline of a few hours, we targeted individual phases independently to detect any excess absorption compared to spectra taken near eclipse, when the absorption signal from the planetary outflow is expected to be minimal. We did not gather in-transit spectra, as these observations have already been conducted in previous studies. Instead, we focused solely on the pre- and post-transit phases at $\varphi\approx\pm 0.1$ and $\varphi\approx\pm 0.2$, allowing for a margin of $\pm 0.05$ to accommodate observational constraints.

2.1 CRIRES+data

We obtained spectra of WASP-52 using the high-resolution spectrograph CRIRES+ (Dorn et al., 2023), mounted on the Very Large Telescope (VLT). Utilizing a $0.2^{\prime\prime}$ slit, we achieved a nominal spectral resolution of approximately R = $100\leavevmode\nobreak\ 000$. To effectively subtract the sky background and reduce instrumental effects, we applied an ABBA nodding pattern, alternating the target between two distinct slit positions (A and B).

The observing nights and according orbital phases are listed in Table 2. On night 1, phases $-0.2$ and $-0.1$ were observed consecutively, yielding a total of 8 spectra (4 x AB nodding). On nights 2 and 3, two exposures were taken at each nodding position, resulting in 4 spectra per night (2 x AB nodding). Unfortunately, phase 0.1 could not be observed due to poor weather conditions. Each spectrum had an integration time of 900 seconds and was taken in the Y-band using the Y1029 wavelength setting.

We used the CRIRES+ pipeline, cr2res, provided by ESO111https://www.eso.org/sci/software/pipelines/index.html#pipelines_table to reduce the data. The raw calibration frames, acquired as part of the observatory’s daily routine, were processed through the standard calibration cascade outlined in the CRIRES+ Pipeline User Manual. This step generated the processed calibration files, including the bad pixel mask, normalized flat field, and the trace-wave file containing the wavelength solution and spectral order positions. Next, we reduced the raw science frames with the cr2res_obs_nodding recipe. This step involved applying the processed calibration files to the science frames, subtracting the B frames from the A frames in each AB (or BA) nodding sequence, and performing optimal extraction to produce one-dimensional science spectra. We then applied the wavelength solution to the extracted spectra.

To enhance the signal-to-noise ratio, we combined two consecutive spectra from each nodding position. The resulting spectrum was then corrected for telluric absorption lines using Molecfit (Smette et al., 2015; Kausch et al., 2015). Figure 2 shows the air mass evolution alongside the mean signal-to-noise ratio (S/N) near the He i 10833 Å line in the telluric-corrected spectra.

2.2 CARMENES data

Similarly, we collected spectra of WASP-52 using the CARMENES spectrograph (Quirrenbach et al., 2018), mounted on the 3.5 m telescope at the Calar Alto Observatory. CARMENES provides simultaneous coverage of both the visual ($5500-9600$ Å) and near-infrared ($9600-17200$ Å) regions, with nominal spectral resolutions of R=$94\leavevmode\nobreak\ 600$ and $80\leavevmode\nobreak\ 400$, respectively. During the observation campaign, each spectrum was obtained with an exposure time of 900 seconds. The observing nights are listed in Table 3.

We processed the raw spectra using the CARMENES pipeline CARACAL (Zechmeister et al., 2014; Caballero et al., 2016), which performs standard data reduction steps, including dark subtraction, flat fielding, wavelength calibration, and spectrum extraction. The pipeline also provides noise estimates for each data point. To enhance data quality, wavelength points with low signal-to-noise ratios (S/N) were masked. We used Molecfit to performed the telluric correction. Figure 3 shows the air mass and S/N near the helium line in the near-infrared, along with the S/N near the H i 6565 Å line for the reduced spectra.

2.3 Searching for planetary excess absorption in the helium line core

Since our observing strategy does not permit the creation of a high signal-to-noise master-out-of-transit spectrum, typically constructed by co-adding multiple out-of-transit spectra, we instead focus on line core equivalent widths (EW) of the stellar He i 10833 Å line. Studies by Gully-Santiago et al. (2024) and Zhang et al. (2023b) show that changes in the stellar equivalent width caused by planetary excess absorption are clearly detectable in their respective light curves (see their Fig. 9 and Fig. 1).

In Figure 4, we present the stellar metastable helium triplet as observed with CARMENES and CRIRES+. CRIRES+ at the VLT uses AB nodding, a technique that effectively eliminates telluric emission features. This method alternates between two positions (A and B) on the sky, capturing both the target spectrum and telluric emission lines. By subtracting the two exposures (A – B), the telluric emission, which is uniform across the slit, cancels out, leaving a cleaner target spectrum. In contrast, CARMENES data still contains telluric emission lines, as it does not use this technique. Hence, while the CARMENES spectra shown (top panel) were corrected for telluric absorption lines with Molecfit, telluric emission lines remain. We chose not to model and subtract the telluric OH emission lines to avoid potential inaccuracies. Instead, we carefully scheduled the observations to ensure that the OH emission lines do not overlap with the core of the helium main component. As such, we focus exclusively on the line core ($\lambda=10832.90-10833.62\,$Å) of the main component, highlighted by the solid vertical lines. We chose this wavelength range in an effort to minimize interference from OH emission lines, though some effect may still remain. Furthermore, Gully-Santiago et al. (2024) and Zhang et al. (2023b), in their observational studies of extended atmospheric escape, found that the planetary absorption signal remains approximately stationary in the stellar rest frame (see their Fig. 10), exhibiting little correlation with the planet’s orbital motion. Based on these findings, we are confident in restricting our analysis to the stellar line core, which is centered in the stellar rest frame. In Appendix A, we present an additional analysis of the full helium line using only the CRIRES+ data, which allows full line coverage due to its corrections for both telluric absorption and emission lines. The qualitative results align with those obtained from the line core analysis presented below.

In the bottom panel of Figure 4, we present the core equivalent width derived from the helium spectra of both data sets, integrating the flux exclusively over the wavelength range corresponding to the helium line core, as indicated.

Furthermore, we use the results by Kirk et al. (2022) (Keck/NIRSPEC) as a reference for planetary excess helium absorption. Using WebPlotDigitizer222https://github.com/automeris-io/WebPlotDigitizer, we extracted the helium line profile of WASP-52 b presented in their work (see their Fig. 7). Since the variability of the stellar line is accounted for by dividing by the master-out-of-transit spectrum, the derived excess absorption is independent of stellar variability. From the extracted line profile, we measured the equivalent width of the helium line core within the same wavelength ranges used for our observations, resulting in an excess absorption core EW of $(23\pm 2)$ mÅ. The uncertainty in this measurement reflects an estimated error margin, as the direct error derived from the data points’ standard deviation is only $1.14$ mÅ. However, due to reliance on WebPlotDigitizer for data extraction and the sensitivity of the EW to data binning, we adopted a more conservative error estimate. We interpreted the pre- and post-transit absorption signals reported by Kirk et al. (2022) to be consistent with zero within the estimated uncertainties. For the out-of-transit phases, we assumed the same noise level and uncertainties as during mid-transit, using the same error value of 2 mÅ. In the bottom panel of Figure 4, we added the results from Kirk et al. (2022) to the baseline of our observations, which itself carries some uncertainty. Consequently, the error bars of the Keck observations represent the combined uncertainties from both the absorption measurement and the standard deviation of the baseline.

From our observations, the EW light curve of WASP-52 b shows no evidence of excess absorption outside the transit, with the measured data points clustering closely around a consistent EW value. Therefore, we interpret all our data points as reflecting a baseline of stellar absorption, with no apparent contribution from the planet. Notably, there is one outlier with lower absorption at phase $-0.11$, which corresponds to the CARMENES spectrum with the strongest OH emission line. We interpret this as the likely cause of the deviation.

To account for the different uncertainties in each data point, we calculated the weighted mean and weighted standard deviation of the equivalent widths using inverse-variance weighting. Each data point was assigned a weight proportional to the inverse square of its associated uncertainty, ensuring that points with smaller errors contributed more significantly to the results. The weighted mean was computed as the sum of the product of the weights and the equivalent widths, divided by the total weight. The weighted standard deviation was then determined by measuring the weighted spread of the data around the weighted mean. We derive a stellar helium line core EW of 194.5 mÅ, with a weighted standard deviation of 2.5 mÅ. Our data sensitivity allows us to detect changes in the equivalent width (EW) of the line core as small as 7.5 mÅ with 3$\sigma$ confidence. This level of precision is sufficient to detect the in-transit excess EW reported by Kirk et al. (2022).

Unfortunately, due to bad weather conditions and technical difficulties, we lack post-transit data at the phase of 0.1 in the near-infrared, which would have been instrumental in exploring potential trailing gas from the planet. Based on the absence of any pre-transit excess absorption, we rule out the presence of a significantly extended leading helium stream like that observed in HAT-P-67 b and HAT-P-32 b. Our observations constrain the extent of the absorbing gas to $31\leavevmode\nobreak\ R_{p}$ ($\varphi=-0.11$) in the leading direction and $73\leavevmode\nobreak\ R_{p}$ ($\varphi=0.26$) in the trailing direction.

3 3D simulations with Athena++

3.1 Simulation setup and models

We simulate a planet in a circular orbit around its host star, incorporating thermal winds from both the planet and the star, using the 3D Eulerian (magneto)hydrodynamic code Athena++333https://github.com/PrincetonUniversity/athena, version 2021 (Stone et al., 2020). For these hydrodynamic simulations, we use the setup previously developed in MacLeod & Oklopčić (2022) and extended to account for planetary day-night anisotropies in Nail et al. (2024). In the following, we outline the most important features of our simulations and the modifications made for this study, building on the methods described in Section 2 of those works. In this approach, we remain agnostic about the detailed physical mechanisms driving the winds. This approach allows us to flexibly explore a broader parameter space. We apply a boundary condition that heats the base of the outflow, and approximate the gas temperature using a nearly isothermal equation of state. This simplification enables us to focus on the large-scale interactions between the planetary and stellar winds.

Within the Athena++ simulations, the equations for mass, momentum, and energy conservation are solved for an inviscid gas. We use the ideal gas law with an adiabatic index of $\gamma=1.0001$, which approximates isothermal behavior along adiabats. Nonetheless, temperatures at the points where planetary winds originate can vary, enabling us to account for temperature differences between the planet’s dayside and nightside, as outlined in Nail et al. (2024). The simulation includes the gravitational influences of the planet and the star, while neglecting the outflow’s back-reaction on the planetary orbit and the effect of He i 10833 Å radiation pressure on gas dynamics.

Hydrodynamic winds are introduced by parameterizing the conditions at the bases of the stellar and planetary winds. The stellar surface pressure is determined by the constant surface density $\rho_{*}$ and the hydrodynamic escape parameter $\lambda_{*}$, which represents the ratio of gravitational potential energy per particle, $GM_{*}mR_{*}^{-1}$, to thermal energy, $k_{\rm B}T(R_{*})$, where $m$ is the mass of the escaping particle. For all models, we set $\lambda_{*}=15$. Similarly, the pressure at the planetary surface is given by

[TABLE]

where $\lambda_{p}=\frac{GM_{p}}{c_{s}^{2}R_{p}}$ denotes the planet’s hydrodynamic escape parameter. In models incorporating day-night anisotropy, surface pressure may vary across the planetary surface (Nail et al., 2024).

We use a spherical polar mesh centered on the star in a co-rotating frame aligned with the planet’s orbit. The computational domain spans from the star’s surface out to a distance of 0.37 au in all directions. The planet is positioned along the negative x-axis at a distance of $a=0.027$ au from the star. Both the star and the planet rotate with the planet’s orbital frequency, $\Omega=\Omega_{\text{orb}}=G(M_{p}+M_{*})a^{-3}$, with their angular momenta oriented in the positive z-direction.

The base mesh consists of $9\times 6\times 12$ mesh blocks, each containing $16^{3}$ zones that are logarithmically spaced in the radial ($r$) direction and evenly spaced in the angular ($\theta$ and $\varphi$) directions, ensuring nearly cubic zone shapes. Near the poles, we mitigate extreme aspect ratios by averaging conserved quantities, effectively reducing the number of zones in the $\varphi$ direction.

To accurately resolve the atmospheric scale height near the planet, we used a mesh scaled to the planetary radius, setting the static mesh refinement (SMR) level to $N_{\mathrm{SMR}}=3$ ensuring increased resolution within a radius of 12 $R_{p}$. We analyze snapshots of the simulation only after a quasi-steady state had been reached after five orbits, defined as the point where the planetary and stellar mass-loss rates stabilized. These rates are calculated by tracking changes in the planet’s and star’s mass over time while accounting for mass flux across the simulation boundaries. This approach ensures mass conservation by incorporating ejected, accreted, and redistributed mass within the system.

We investigate four models to explore the influence of stellar and planetary wind properties on the outflow morphology of WASP-52 b (see Tab. 4):

• $\bullet$

ISO (Isotropic planetary wind): A spherically symmetric planetary outflow, interacting with a moderate stellar wind ($\dot{m}_{*}/\dot{m}_{p}\approx 5$).

• $\bullet$

ANISO (Anisotropic planetary wind): Same as ISO, but accounting for a day-night pressure gradient following the approach in Nail et al. (2024), with the nightside pressure reduced to 10% of the dayside value.

• $\bullet$

TORUS (Reduced stellar wind): Reduced stellar mass-loss rate relative to model ISO ($\dot{m}_{*}/\dot{m}_{p}\approx 0.5$).

• $\bullet$

TAIL (Enhanced stellar wind): Enhanced stellar mass-loss rate relative to the ISO model ($\dot{m}_{*}/\dot{m}_{p}\approx 16$).

3.2 Post-processing of 3D simulations

To post-process the Athena++ simulation snapshots, we introduce a new framework that connects Athena++ snapshot outputs (3D pressure and density profiles) with Cloudy444gitlab.nublado.org/cloudy/cloudy v23.01 (Ferland et al., 1998, 2017; Chatzikos et al., 2023; Gunasekera et al., 2023), a non-local thermodynamic equilibrium (NLTE) plasma simulation code.

Using atmospheric density and temperature profiles from Athena++ and an incoming stellar spectral energy distribution (SED), Cloudy determines the electron temperature, ionization levels, and excitation states of various atomic species. Although Cloudy usually computes the electron temperature, it can also accept a predefined temperature profile, as is done here using data from the Athena++ simulation. In reality, photoionization by stellar EUV/XUV photons can locally heat the gas above the parameterized temperatures, particularly in low-density regions. However, we chose the Athena++ temperature profiles to maintain consistency with the hydrodynamic evolution of the wind. Our approach offers a computationally efficient compromise, combining realistic 3D dynamics with detailed 1D radiative transfer, though it neglects lateral radiative effects and self-consistent temperature evolution.

To simulate transmission spectra originating from orbital phases similar to the observations (see Fig. 1), one-dimensional rays are generated for the corresponding observing angles along the line-of-sight between the observer and the star. The code first calculates a direction vector based on the azimuthal and polar angles, then determines the number of rays required to cover the stellar surface based on the observing angle and position of the planet. To ensure higher resolution near the planet, rays are allocated more densely in that region, with their weighting in the spectra being adjusted according to their relative area, following the approach of the previous radiative transfer analysis described in MacLeod & Oklopčić (2022). The radial extent of each ray is set by user-defined inner and outer limits (here from -20 to 40 $R_{p}$), as well as the mesh grid containing the ray. To assign values to the zone centers, we use nearest-neighbor interpolation from the simulation grid, ensuring that any point within a zone’s volume receives the corresponding value of that zone.

Density and pressure data from these rays, as well as the derived temperature

[TABLE]

assuming a mean molecular weight of $\mu=0.6$, were extracted from Athena++ and provided to Cloudy.

Since Cloudy performs radiative transfer calculations, it does produce a spectrum that is transmitted through each ray. However, since it is a hydrostatic code, this spectrum does not include spectral line broadening due to the radial velocities. Therefore, instead of using Cloudy’s own predicted spectrum, we instead extract from it the number densities of each energy level of each different element, and use these to calculate the transit spectrum in post-processing. In this way, we are able to include Doppler shifts and broadening of the spectral lines based on the velocity structure derived from Athena++.

Using the number densities of the energy levels for each ray as calculated by Cloudy, assuming solar composition, we performed radiative transfer calculations for the metastable helium triplet. For the radiative transfer, spectral line data for helium were sourced from the NIST database555www.nist.gov/pml/atomic-spectra-database, such as wavelength, energy levels, oscillator strengths, and Einstein coefficients.

The optical depth $\tau$ at each frequency $\nu$ was calculated along the path length d$s$ of the rays:

[TABLE]

Equation 3 takes into account the species number density $n_{s}$, the line cross-section $\sigma_{0}$, and the Voigt $\Phi(\Delta v)$ profile, which includes a Lorentzian component due to natural line broadening and a Gaussian component that depends on the local gas temperature $T$. Additionally, the wavelength shift $\Delta\nu$ reflects the gas velocity in the inertial frame, projected along the line of sight. The extinction of each ray is calculated by $1-\exp(-\tau_{\nu})$. The rays for a given observing angle are weighted according to the quadratic limb darkening law, and their extinction contributions are summed to compute the total transmission spectrum. We tested the convergence of the spectra by varying the total number of rays used in the radiative transfer calculation, finding that approximately 1000 rays are sufficient to achieve reliable results (see Appendix B).

A spectral energy distribution (SED) for WASP-52 is unavailable, thus we adopted the SED of a comparable star, HD 40307, which has a spectral type of K2.5 and $T_{\rm eff}=4880$ K, provided by the MUSCLES Treasury Survey (France et al., 2016), and scaled it to the orbital distance of WASP-52 b. The SED serves as a close approximation to the reported properties of WASP-52 b’s host star, which has a spectral type of K2.0 and $T_{\rm eff}=5000$ K (Hébrard et al., 2013). The age of HD 40307 is estimated to be 2.6 Gyr (Turnbull, 2015), while Hébrard et al. (2013) estimated WASP-52 to be 0.4 Gyr. However, WASP-52 is an active star, and its actual age remains uncertain, leading to potential differences in their XUV fluxes. If WASP-52 is younger and more active, its XUV flux would be higher; conversely, if it is older, as suggested by Mancini et al. (2017), its XUV flux would be lower. These variations could affect the absolute populations of metastable helium in the gas to some extent.

To assess the performance of the new radiative transfer framework, we compared its results to those obtained with the previous code introduced by MacLeod & Oklopčić (2022), which focuses solely on the He i 10833 Å line, presented in Appendix C. We found that the main difference between the codes is the inclusion of doubly ionized helium in Cloudy. The production of this species significantly reduces the metastable helium number density, as it is primarily produced through the recombination of singly ionized helium. In Cloudy, less singly ionized helium is available due to the new He2+ pathway. The previous code uses fixed values of recombination rates and collision cross-sections that are well suited for modeling planetary outflows of temperatures around $10^{4}$ K, while Cloudy incorporates more accurate values across a broader temperature range. Therefore, the new radiative transfer framework is better suited for calculating abundances in large-scale 3D simulations of atmospheric escape, especially those aimed at studying stream-like outflows, as they might include a wide range of densities and temperatures a factor of a few lower than $10^{4}$ K. In Figure 9, we compare the helium and hydrogen species number densities along rays generated by the different codes. The largest differences are observed in low-density environments, where a significant population of He2+ can build up due to the low probability of recombination.

Additionally, we assess the physical timescales relevant to the population and observability of helium in the metastable triplet state. Based on Cloudy output, we estimate the total photoionization timescale of $\sim 1$ hour for neutral helium. The radiative lifetime of the metastable state is about 2.2 hours (Drake, 1971). The hydrodynamic escape timescale, estimated as the time required for gas to flow from the planetary surface to the Hill radius ($r_{\rm H}\approx 2.5R_{p}$), is approx. 5.7 hours, assuming a sound speed of 11 km s*-1*. Other depopulation channels, such as collisional de-excitation by electrons, are included in our post-processing but are much slower. Importantly, ionizations and recombinations occur continuously throughout the gas, and the metastable state reflects a statistical population of many atoms rather than the history of any single one. These timescale comparisons demonstrate that while ionization equilibrium is not strictly assumed, the metastable population is maintained under quasi-steady-state conditions captured by our non-equilibrium post-processing approach.

We also compared temperature profiles along a representative ray near the planet from both Athena++ and Cloudy. Athena++ assumes a nearly isothermal wind, while Cloudy’s predicted temperature varies with the density profile. Although their temperature structures differ on large scales, the resulting spectra agree to within  1%, with Cloudy showing slightly stronger absorption due to marginally higher temperatures near the planet. Achieving self-consistency between photoionization heating and dynamics would require iterative coupling, as implemented in the 1D sunbather code (Linssen et al., 2024), which links a Parker wind solver with Cloudy. Notably, sunbather yields base wind temperatures ($\approx 8500$ K) similar to Athena++. We adopt the Athena++ temperature profiles to maintain consistency with the hydrodynamic structure.

3.3 Simulation results

We modeled the atmospheric escape of WASP-52 b hydrodynamically in 3D using Athena++ as described in Section 3.1. As a first step, we wanted to identify a simulation setup that reproduces the mid-transit helium absorption observations reported by Kirk et al. (2022). Once this was established, we extended the analysis to explore additional out-of-transit phases, providing insights into the potential morphology of escaping tails or streams, analogous to our observational strategy (see Fig. 1). To generate the synthetic spectra based on these 3D simulation snapshots, we used the radiative transfer framework described in Section 3.2.

In our observational analysis (Sect. 2), we focused on the line core of the metastable helium triplet, extracting the equivalent width of the line core (see Fig. 4, bottom panel). To compare our models to observations, we calibrated the synthetic spectra generated from the 3D hydrodynamic simulations to match the mid-transit line core EW observed by Kirk et al. (2022).

To achieve this calibration, we rescaled the densities and pressures of the simulation snapshots, using a scaling factor $S$. The hydrodynamic equations that govern the outflow are scale-invariant, which means that we can adjust the pressure and density after the simulation without affecting the overall morphology of the wind-wind interactions. Crucially, this approach preserves the ratio of mass-loss rates between the planetary and stellar winds, ensuring the relative dynamics remain consistent. This rescaling allows us to match the simulated line core EWs to the observed mid-transit values reported by Kirk et al. (2022), without changing the fundamental properties of the outflow. The scaling factors we applied, along with the resulting time-averaged mass-loss rates $\langle\dot{m}\rangle$ for the star and planet, are summarized in Table 5.

As a starting point, we explored planetary hydrodynamic escape parameters $\lambda_{p}$ between 3 and 9. For a direct comparison with the observations of Kirk et al. (2022), we convolved the mid-transit synthetic spectra with a Gaussian to match the Keck/NIRSPEC resolution ($R=25\leavevmode\nobreak\ 000$). The results are shown in Appendix D. A low escape parameter ($\lambda_{p}=3$), corresponding to a hot planetary wind of T\approx 15\leavevmode\nobreak\ 700\leavevmode\nobreak\ K, produces a highly thermally broadened line profile that is inconsistent with the observed line width. Conversely, a high escape parameter ($\lambda_{p}=9$), representing a cold planetary wind of $T\approx 5200$ K, results in narrow line profiles that do not match the observed broadening. Unlike our previous study (Nail et al., 2025), we do not find kinematic broadening of the helium line for our WASP-52 b simulations. This broadening typically occurs when the line forms in a high-density outflow, a scenario inconsistent with the modest excess helium absorption observed for WASP-52 b. These findings rule out both low and high $\lambda_{p}$ scenarios.

The best agreement with the observed line width is achieved for $\lambda_{p}=5$, corresponding to a planetary wind temperature of approximately T\approx 9400\leavevmode\nobreak\K. Building on this result, we explored four models: one assuming isotropic winds for both the planet and the star (ISO), and three variations involving adjustments to the stellar wind strength (TORUS, TAIL) and the inclusion of a day-night anisotropy in the planetary wind (ANISO).

Figure 5 presents the scaled snapshots of the four simulation models, captured after five orbits, by which time a quasi-steady state had been reached. In the ISO model, the planetary wind extends up to an orbital phase of approximately $-0.2$, suggesting that it could be probed during pre-transit observations. Near the planet, the outflow shows characteristics intermediate between a bubble and stream morphology, as described by MacLeod et al. (2024). A cavity of unshocked gas is still present close to the planet, but it is stretched along the orbit by tidal forces.

In the ANISO model, the cavity on the planet’s night and trailing sides disappears due to the reduced gas pressure on the nightside, which is insufficient to counteract the stellar wind. Consequently, the gas in these regions becomes completely shocked. On the planet’s dayside, however, the ram pressure of the planetary outflow is high enough to enable the formation of an extended cavity that is accelerated along the leading direction of the planet’s orbit.

The TORUS model represents a scenario where the stellar wind is too weak to confine the planetary outflow. As a result, the planetary wind escapes freely and is not significantly deflected or compressed. The gas lingers in the planet’s orbital path, forming a torus-like structure around the star. In contrast, the TAIL model shows a scenario where the planetary wind is highly confined by a strong stellar wind, resulting in the formation of a dense, trailing tail and the truncation of the leading arm.

3.4 Synthetic helium results

Using the snapshots described above, we generated synthetic helium spectra for different observing angles. The resulting line profiles, shown alongside the observational mid-transit data from Kirk et al. (2022) in Figure 11, allow for a direct comparison.

In this section, however, we focus on the helium line core equivalent width (EW), following the approach used in the observational data analysis in Section 2. To do so, we integrate over the core of the spectral line within the same wavelength range specified in the observational study. The selected orbital phases for reference include $\pm 0.2$, $\pm 0.1$, and 0. Additionally, we examine ingress and egress phases at $\varphi=\pm 0.02$, as well as $\varphi=\pm 0.04$, corresponding to the immediate pre- and post-transit points, just before the first contact point ($t_{1}$) and after the fourth contact point ($t_{4}$). Although our observations did not cover these additional phases, analyzing them helps constrain the density gradients near the planet predicted by different models and their impact on the expected EW light curves. The integrated EW results for the four models are shown in Figure 6. We convolved the synthetic spectra with the instrument resolution of Keck/NIRSPEC (R = 25 000) to match the simulation density scaling with the observational results. Therefore, the synthetic line-core EW presented in Figure 6 are slightly higher than the mid-transit observational data point, as the figure displays the results of the unconvolved spectra. This figure is intended as a comparison to the observations conducted in this study, and given the high resolving power of the spectrographs we used, the impact on the line profile is negligible.

The EW light curve shows that for all models, no significant out-of-transit absorption in the He i 10833 Å line core would be apparent in our observations. This means that even if there is a torus of gas lingering in the orbit of the planet, as shown in the TORUS model, no excess absorption is expected to be detected, considering the sensitivity of our observations. Therefore, we cannot definitively rule out any of the scenarios presented here. In all models, except the ANISO model, enhanced helium absorption could be detected during the transit contact point $t_{4}$ in the trailing direction. For the ANISO model, the night-side pressure is too low to form a tail, and the velocity shift appears too high (see Fig. 11). This suggests that a strong day-night anisotropy is unlikely for WASP-52 b, which may explain the lack of a blueshift in the He i 10833 Å line during transit.

4 Discussion

4.1 Outflow morphology of WASP-52b

Our He i 10833 Å line observations suggest that the atmospheric outflow of WASP-52 b is not detectably extended. We constrain its extent to $<73\leavevmode\nobreak\ R_{p}$ ($\varphi=0.26$) in the trailing and to $<31\leavevmode\nobreak\ R_{p}$ ($\varphi=-0.11$) in the leading direction. Comparisons with 3D models suggest that the observed in-transit absorption line width is best explained by thermal broadening at a temperature of T\approx 9400\leavevmode\nobreak\K, corresponding to a moderate hydrodynamic escape parameter of $\lambda_{p}=5$.

In our previous work (Nail et al., 2025), which explores the outflow morphology of HAT-P-67 b and HAT-P-32 b, we examined the kinematic broadening of the helium line in a stream-like morphology, where the outflow diverges into two cold streams, corresponding to a high $\lambda_{p}\approx 9$ value and a low temperature of $T\approx 4700$ K at the dayside of the planet. This results in a broad and irregular helium line profile, with kinematic broadening dominating over thermal broadening. However, the in-transit spectra of WASP-52 b do not show kinematic broadening, suggesting that the outflow likely does not follow a stream-like morphology.

Certain limitations have affected the scope of our analysis. One challenge was the lack of NIR observations at phase 0.1 due to a readout error in CARMENES and unfavorable weather conditions during the CRIRES+observations. This phase would have been particularly valuable for probing the presence of a trailing tail of escaping gas, as suggested by JWST/NIRISS observations from Fournier-Tondreau et al. (2024). They report a tentative helium detection at 2.9$\sigma$ during egress, which could be attributed to a trailing tail of escaping gas. Although their observations were taken at a much lower resolution (R $\approx$ 700), making direct comparison and translation to the high-resolution spectra presented in this work challenging, their results are broadly consistent with the TAIL, TORUS, and ISO models presented in this study. For a direct comparison, see Figure 7 in Appendix A. Specifically, both egress and ingress observations with JWST show enhanced helium absorption, with egress absorption being slightly stronger. The ANISO model, however, can be ruled out for two reasons: first, the light curve does not show enhanced absorption during egress (see Fig. 6), and second, the mid-transit velocity shift of the helium line appears too high (-4.5 km s*-1*, see App. D), whereas Kirk et al. (2022) reported a shift of $0.00\pm 1.19$ km s*-1*.

Although WASP-52 b and HAT-P-32 b have similar Hill radii and surface gravities, leading to comparable Rossby numbers, important for distinguishing between stream and bubble morphologies as predicted by MacLeod et al. (2024) (see Fig. 9 therein), the helium outflow from WASP-52 b appears significantly less extended. While the results from Nail et al. (2025) suggest that HAT-P-32 b has a stream-like outflow (high $\lambda_{p}$), this raises the question of what differentiates these two planets in terms of atmospheric escape, or what causes the variation in their escape parameters. One possible explanation lies in differences in optical depth, which may be influenced by a higher planetary mass-loss rate as well as variations in the population of metastable helium. Alternatively, assuming similar tidal and mass-loss conditions, as well as composition, differences in outflow temperature could also play a role. Under this assumption, the outflow temperature of WASP-52 b seems to be higher (lower $\lambda_{p}$) than that of HAT-P-32 b.

One possibility is that a different wind-driving mechanism is at play. According to Owen & Schlichting (2024), the Bondi radius (the radius at which bolometrically heated gas becomes unbound) lies below the altitude of the XUV-heated region if core-powered mass-loss is dominant. The altitude of this XUV-heated region, depends strongly on atmospheric composition and the stellar XUV flux. Given that WASP-52 is highly active, it may produce a stronger XUV flux than HAT-P-32, leading to a deeper XUV penetration depth and, as a result, a hotter, photoevaporation-driven outflow. Although our understanding of heat redistribution and gas circulation at very low pressures in the upper atmosphere is limited, one could speculate that heat at the wind launching altitude might be distributed efficiently. If heat redistribution were less effective, we would expect a strong day-to-nightside temperature contrast in the photoevaporation scenario. However, our models show that a pronounced day-night anisotropy (ANISO) seems unlikely, suggesting that other factors, such as atmospheric composition or thermal structure, may play a more significant role in shaping the outflow.

We suggest that the new radiative transfer framework introduced in this work may help resolve some of these uncertainties, as it allows us to calculate the number densities of multiple species simultaneously in a 3D outflow. This capability opens up the possibility of identifying potential tracers of atmospheric escape and the composition of the outflow. Different species probe different line-forming regions: lower altitudes, where the wind is launched, and higher altitudes, where the wind interacts with the stellar wind. Helium, for instance, serves as a good probe of the upper atmosphere due to atmospheric fractionation, while heavier species could be probing the lower atmosphere.

An important consideration in the study of WASP-52 b is the distinction between the physical extent of its outflow and its detectable extent. While simulations indicate the presence of an extended atmosphere, the observed signal depends on the outflow’s geometry, optical depth, and level populations. The densities in the extended outflow may be too low to effectively shield ionizing radiation, potentially contributing to the production of doubly ionized helium. Additionally, recombination is sensitive to gas density; in a low-density environment, He2+ ions encounter electrons less frequently, reducing the likelihood of recombination into He+. Consequently, while the material is present, its transparency in the observed helium line could explain the weak or absent signal during out-of-transit observations. This highlights the importance of considering not just the presence of an outflow, but also its radiative properties. Alternative tracers that are more sensitive to conditions in lower-density regions, such as ionized species that remain detectable even at lower column densities, may offer additional insight into the structure and composition of the outflow.

4.2 Challenges from stellar activity in detecting extended atmospheric features

Stellar activity remains a potential complication when studying WASP-52 b as a test case for atmospheric escape. Previous observations of the planet have shown that stellar activity can interfere with the interpretation of atmospheric escape tracers (e.g. Bruno et al., 2020; Chen et al., 2020; Fournier-Tondreau et al., 2024), as atomic lines observed in the planetary atmosphere are also chromospherically sensitive lines. Cauley et al. (2018) demonstrated that heterogeneities on the stellar surface can introduce spurious in-transit absorption signals in high-resolution transmission spectra. Their simulations showed that unocculted active regions can produce apparent absorption depths of up to 0.3% in H$\alpha$ and 0.2% in the Na i D lines. By contrast, they found that the He i 10833 Å line is largely unaffected by stellar activity; in many cases, activity tends to dilute rather than mimic the absorption signal. This, of course, strongly depends on the location of the active regions relative to the transit chord. If active regions are predominantly unocculted, they can weaken the contrast between the stellar and planetary signals, whereas occulted active regions may still introduce subtle distortions. Nonetheless, compared to other chromospheric lines like H i 6565 Å , the He i 10833 Å triplet seems to be one of the most reliable tracers for atmospheric escape in active stars. This further motivated our focus on this line in our study.

Given the potential impact of stellar activity on transmission spectra, one approach to disentangle stellar and planetary contributions is to monitor both H$\alpha$ and He i simultaneously. Investigating a potential anticorrelation between the two lines could indicate stellar activity (Guilluy et al., 2020). This highlights another advantage of our new radiative transfer framework, which can help quantify the planetary contribution of atmospheric escape tracers that also probe stellar activity. However, a limitation of the present study is the lack of spectroscopic coverage during the immediate pre- and post-transit phases. These orbital phases are particularly valuable for identifying extended or asymmetric features, such as trailing tails or leading arms of escaping material. In an ideal future experiment, these phases would be sampled more densely, with higher cadence observations in both the near-infrared and visible wavelength ranges. Obtaining at least five or more data points outside of transit, would improve our ability to track variability in the line profiles and to test whether H i 6565 Å and He i 10833 Å signals are anti-correlated. This would help to separate stellar and planetary contributions in the tracer lines.

To better constrain the baseline variability of the stellar lines, a long-term monitoring campaign during the secondary eclipse, or at other out-of-transit phases, would also be valuable. Such monitoring, carried out over multiple nights, could characterize the intrinsic variability of stellar chromospheric lines and provide the context needed to confidently interpret in-transit absorption features.

5 Summary

In this study, we present high-resolution observations of the helium triplet at 10833 Å using the spectrographs CRIRES+and CARMENES, covering an extended range of orbital phases ($\varphi\approx\pm 0.1,\pm 0.2,0.5$) to investigate the extent of the planet’s helium outflow. To explore possible outflow morphologies, we performed 3D simulations with Athena++, testing four distinct scenarios and comparing them to the observational data. Additionally, we used an improved radiative transfer post-processing approach, extending our previous methods to generate synthetic spectra, which could potentially be used to model other atmospheric escape tracer lines. Our results on the atmospheric escape of WASP-52 b include:

• $\bullet$

Our helium observations constrain the extent of metastable helium in the orbit of WASP-52 b, ranging from a minimum at phase $-0.11$ ($31\leavevmode\nobreak\ R_{p}$) to a maximum at phase $0.26$ ($73\leavevmode\nobreak\ R_{p}$) in the trailing direction.

• $\bullet$

We consider a stream-like morphology, as suggested for HAT-P-67 b and HAT-P-32 b, to be unlikely for WASP-52 b. Instead, we propose that its outflow morphology likely represents an intermediate regime between bubble and stream structures. This conclusion is supported by a hydrodynamic escape parameter of approximately $\lambda=5$, corresponding to an outflow temperaure of T\approx 9400\leavevmode\nobreak\K, which is required for our 3D simulations to reproduce in-transit observations by Kirk et al. (2022).

• $\bullet$

Assuming a fixed hydrodynamic escape parameter, the extent and shape of the planetary wind, whether it forms a torus or a trailing tail, depend on the stellar mass-loss rate; however, from our modeling perspective, the two scenarios appear indistinguishable, as the metastable helium densities in a potentially extended structure are too low to be detected with our current observations.

• $\bullet$

A high day-night anisotropy scenario is unlikely, as it would produce a blueshift in the helium line during transit, which is not observed by Kirk et al. (2022). Furthermore, such a scenario would result in higher helium absorption during ingress than during egress, which contradicts the low-resolution observations with JWST by Fournier-Tondreau et al. (2024), where egress absorption is observed to be slightly stronger.

• $\bullet$

The new radiative transfer framework presented here, which uses Cloudy to calculate the number densities of multiple species in the 3D simulations of wind-wind interaction between the star and planet, has the potential to model multiple atmospheric escape tracer lines and provide more robust results in low-temperature and low-density regimes. When combined with observations, it offers great potential for accessing the planetary wind structure in terms of composition and temperature.

6 Data availability

The observational data, simulation snapshots, and Jupyter notebook for reproducing the figures666https://doi.org/10.5281/zenodo.16687147, along with the new radiative transfer framework and tutorial notebooks for running Athena++ simulations and Cloudy post-processing777https://doi.org/10.5281/zenodo.16687969, are all publicly available on Zenodo.

Acknowledgements.

Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under program ID 112.25M6.001. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101004719. We thank SURF (www.surf.nl) for the support in using the National Supercomputer Snellius. A. Oklopčić gratefully acknowledges support from the Dutch Research Council NWO Veni grant.

Appendix A Full helium line analysis

Figure 7 shows the equivalent width (EW) measurements derived from a full-line integration over the helium triplet feature ($\lambda=10831.556-10834.300$ Å), using only the CRIRES+ data. The CRIRES+ spectra enable such an analysis, as they provide complete coverage of the line profile without contamination from telluric emission lines. For consistency, we applied the same integration procedure to the four simulation models presented in this study (see Tab. 4). The CRIRES+ EW measurements cluster around a nearly constant value of $190.6\pm 2.1$ m$\AA$, which we interpret as the contribution from the stellar line alone. To compare the observational and simulated results, we added this stellar baseline offset to the model EWs.

We also include data points extracted from the JWST/NIRISS observations published by Fournier-Tondreau et al. (2024), which offer higher time resolution but considerably lower spectral resolution ($R\approx 700$) compared to the KECK/NIRSPEC data presented by Kirk et al. (2022) (see Fig. 4 and 6). To enable a direct comparison with the high-resolution analysis presented in this work, we digitized the flux values (reported in ppm) from their Figure 11 using WebPlotDigitizer888https://github.com/automeris-io/WebPlotDigitizer, and converted them into fractional excess absorption. To derive the equivalent width (EW), we first converted the absorption values from parts per million (ppm) to unitless fractional values by dividing by $10^{6}$. Assuming a Gaussian instrumental profile, the equivalent width was then computed using the following relation:

[TABLE]

where $f$ is the fractional excess absorption, $\lambda_{0}$ is the central wavelength of the helium feature, and $R$ is the spectral resolution of the instrument. The quantity $\Delta\lambda=\lambda_{0}/R$ corresponds to the effective width of the resolution element or kernel area. For the JWST/NIRISS data, we adopted $\lambda_{0}=10830$ Å  which is the central wavelength reported in their analysis, and $R=700$. The uncertainty in the flux measurements is approximately 1000 ppm, which corresponds to an equivalent width error of 15.5 mÅ. For clarity of presentation, we omit these error bars in Figure 7.

One advantage of using the full wavelength range of the line, is that it captures planetary signals that are shifted due to the planet’s motion. This approach is also more robust to data sensitivity, as it integrates over a broader wavelength range. By using the higher time-resolution JWST data for comparison, we obtain a better sampling of the ingress and egress phases, providing valuable insights into the immediate extent of the atmosphere. The JWST data suggest that egress and post-transit absorption are slightly enhanced compared to ingress and pre-transit, which could indicate a trailing, escaping atmosphere. However, this conclusion is tentative and requires further investigation.

Appendix B Convergence testing of the new radiative transfer framework

To assess the numerical convergence of the new radiative transfer framework, described in Section 3.2, we tested its sensitivity to the number of rays used in the spectrum calculation. Specifically, we computed the helium absorption spectrum for the TORUS model while varying the total number of rays. In this test, we distinguish between regions of differing density: a high-density region, probed during mid-transit observations at phase $\varphi=0$, and lower density region, probed for example during pre-transit observations, here at phase $\varphi=-0.024$.

Figure 8 shows the resulting helium absorption spectra computed with varying numbers of rays. Spectra calculated with approximately 100 to 300 rays show noticeable fluctuations, particularly in low-density regions. However, the spectra converge rapidly as the ray count increases, with stable results achieved around 900 rays. While a slight difference in the spectrum is observed between 1000 and 3000 rays, the additional computational cost does not yield a significant improvement in accuracy. This indicates that using approximately 1000 rays per model snapshot provides sufficient convergence for the cases examined in this study.

Appendix C Comparison of radiative transfer methods

The previous radiative transfer (RT) code, first introduced in MacLeod & Oklopčić (2022) and further used in Nail et al. (2025) and Nail et al. (2025), was built for high-temperature regimes of $10^{4}$ K. However, planetary winds formed under lower hydrodynamic escape parameters, such as those seen in stream-morphology, have significantly lower temperatures on the order of $10^{3}$ K. Therefore, we were motivated to improve our RT method by incorporating the more sophisticated NLTE plasma simulation code Cloudy (see Sect. 3.2), that gives access to a more complex chemical network and that reacts more appropriately to different temperature regimes.

To further evaluate the differences between the previous and new RT methods, we compared the spectra and species number densities from a snapshot of the TORUS model, as presented in Section 3. The comparison is shown in Figure 9, which illustrates the helium spectrum resulting from the radiative decay of the metastable helium state (top panel), along with the number densities of hydrogen and helium species (middle and bottom panel) for both RT methods. While Cloudy provides access to an even broader range of species, we focus here on those included in the previous code.

The dotted lines in Figure 9 indicate results obtained using the new RT code implemented in this study, while the solid lines represent results from the previous RT code. The middle panel displays the number density profiles of a ray that closely crosses the planet in the $yz$-plane (the plane of observation during the planet’s transit in front of the star) at $\mathrm{YZ}(0.16|0.02)R_{*}$, whereas the bottom panel illustrates the number densities for a ray that passes farther from the planet at $\mathrm{YZ}(0.86|-0.28)R_{*}$, both in a mid-transit configuration $(\varphi=0)$.

While hydrogen abundances are similar in both codes, the helium species show notable differences. In the new RT code, neutral helium abundances are lower, particularly in low-density regions far from the planet. For regions very close to the planet (within a few planetary radii) the number density profiles of both codes are nearly identical, but differences become more pronounced farther out, as seen in the bottom panel. This discrepancy is due to the formation of He2+ in the new code, which the previous code does not account for. At distances from the planet greater than $\sim 5\leavevmode\nobreak\ R_{p}$, most helium in the new RT model transitions to He2+, reducing the availability of He*+* for recombination into the metastable He(23S) state responsible for the spectral lines at 10833 Å. For a detailed description of helium transitions, we refer to Fig. 1 in Oklopčić & Hirata (2018).

To further illustrate this point, Figure 10 compares the helium spectra produced by the new and previous RT codes. In this test, we include only simulation cells that exceed a given optical depth threshold in the metastable helium component ($\Delta\tau$), gradually relaxing the threshold to include more of the low-density material. As expected, both codes yield nearly identical spectra when only high-opacity $\Delta\tau>10^{-2.5}$ (i.e., high-density, close-in) regions are considered. However, as lower-opacity cells, located farther from the planet, are included, the spectra begin to diverge. This behavior supports the idea that the main difference between the two codes arises from the treatment of doubly ionized helium.

The absence of He2+ formation in the previous RT code is acceptable for 1D modeling, which is centered on the planet, but poses a limitation in more complex 3D outflows, in which many light rays do not come close (within a few planet radii) to the planet. This comparison emphasizes the importance of including detailed atomic processes, like He2+ formation, to accurately capture low-density outflow regions.

Furthermore, we tested whether the inclusion of metals in Cloudy affects the spectra. However, excluding metals yielded the same spectra as including them, allowing us to narrow the focus to the formation of doubly ionized helium.

Appendix D Comparison of 3D model spectra with in-transit observations

As previously described, we tested various values of $\lambda_{p}$ for the simulation, which dictates the outflow temperature of the planetary wind (see Eq. 2), and consequently, the width of the line due to thermal broadening (see top panel in Figure 11). We chose $\lambda_{p}=5$ as the optimal value, which corresponds to an outflow temperature of $T\approx 9400$ K.

The bottom panel of Figure 11 presents the results for the four models presented in Section 3. All models align with the observational data, extracted from Fig. 7 of the publication using WebPlotDigitizer; however, the ANISO model shows a blueshift of $-4.5$ km s*-1*, which does not match the observational data, which was reported as $0.00\pm 1.19$ km s*-1*. Nevertheless, we wanted to explore this scenario to investigate how it would affect the symmetry of the outflow.