Stringent Upper Bounds on Atmospheric Mass Loss from Three Neptune-Sized Planets in the TOI-4010 System

TL;DR

This study uses helium transit observations to set stringent upper limits on atmospheric mass loss rates for three Neptune-sized planets in the TOI-4010 system, challenging existing hydrodynamic escape models.

Contribution

It provides the first upper bounds on atmospheric escape for multiple planets in a single system, highlighting discrepancies with standard hydrodynamic models.

Findings

Upper bounds on helium absorption are below 1.3%.

Mass loss rate limits are around 10^10.5 g/s.

Results challenge predictions from one-dimensional hydrodynamic models.

Abstract

Photoevaporative models predict that the lower edge of the Neptune desert is sculpted by atmospheric mass loss. However, the stellar high energy fluxes that power hydrodynamic escape and set predicted mass loss rates can be uncertain by multiple orders of magnitude. These uncertainties can be bypassed by studying mass loss for planets within the same system, as they have effectively undergone scaled versions of the same irradiation history. The TOI-4010 system is an ideal test case for mass loss models, as it contains three Neptune-sized planets with planet b located in the `Neptune desert', planet c in the `Neptune ridge', and planet d in the `Neptune savanna'. Using Keck/NIRSPEC, we measured the metastable helium transit depths of all three planets in order to search for evidence of atmospheric escape. We place upper bounds on the excess helium absorption of 1.23\%, 0.81\%, and 0.87\% at 95\% confidence for TOI-4010~b, c and d respectively. We fit our transmission spectra with Parker wind models and find that this corresponds to 95th-percentile upper limits of $10^{10.17}$g~s$^{-1}$, $10^{10.53}$g~s$^{-1}$, and $10^{10.50}$g~s$^{-1}$ on the mass loss rates of TOI-4010~b, c, and d respectively. Our non-detections are inconsistent with expectations from one-dimensional hydrodynamic models for solar composition atmospheres. We consider potential reductions in signal from a decreased host star XUV luminosity, planetary magnetic fields, enhanced atmospheric metallicities, and fractionation, and explore the implications of our measurements for the past evaporation histories of all three planets.

Full text

††facilities: ADS, NASA Exoplanet Archive, Keck II/NIRSPEC††software: AIOLOS (Schulik & Booth, 2023), arviz (Kumar et al., 2019), astropy (Astropy Collaboration et al., 2013, 2018, 2022), emcee (Foreman-Mackey et al., 2013), matplotlib (Hunter, 2007), numpy (Harris et al., 2020), photoevolver (Fernández Fernández et al., 2023), p-winds (Dos Santos et al., 2022), pyTPCI (Rosener et al., 2025), scipy (Virtanen et al., 2020)

\savesymbol

tablenum

\restoresymbolSIXtablenum

1 Introduction

Planets on close-in orbits can lose significant portions of their atmospheres to photoevaporative mass loss (e.g., Owen, 2019). It is thought that atmospheric mass loss may carve out part of the ‘Neptune desert’ (see Figure 1), which is defined by a lack of sub-Neptune to Neptune-sized planets on close-in orbits (Szabó & Kiss, 2011; Beaugé & Nesvorný, 2013; Lundkvist et al., 2016). Both theoretical (Owen & Lai, 2018; Ionov et al., 2018) and observational work (Vissapragada et al., 2022; Guilluy et al., 2023) suggest that the higher gravities of planets near the upper edge of the Neptune desert make them effectively invulnerable to photoevaporation. As a result, it has been proposed that the shape of the upper edge is created by tidal disruption of planets that came too close to the star during high eccentricity migration (Owen & Lai, 2018). The same theoretical models predict that the lower boundary of the desert is likely sculpted by atmospheric escape, and planets in this region are thought to have either formed in situ or migrated inward via interactions with the protoplanetary gas disk (Batygin et al., 2016; Bailey & Batygin, 2018; Ida & Lin, 2008).

Our current story for the lower edge of the Neptune desert is based on predicted mass loss rates from population-level photoevaporative models, which can be uncertain by multiple orders of magnitude. These uncertainties come from poorly-constrained system properties that control mass loss, including the stellar Extreme Ultra-Violet (EUV) and X-ray (together XUV) radiation flux (e.g., Tu et al., 2015), stellar wind properties (Carolan et al., 2020; Mitani et al., 2022; MacLeod & Oklopčić, 2022; Rumenskikh et al., 2023), the presence or absence of strong planetary magnetic fields (Schreyer et al., 2023), and planetary atmospheric compositions (e.g., Zhang et al., 2022b). In energy-limited mass loss models, the total amount of mass lost is directly tied to the EUV and X-ray output of the host star. Most of this mass loss is thought to occur within the first $\sim$100 Myrs - 1 Gyr of the star’s lifetime, when its UV and X-ray output is highest (Güdel, 2004; King et al., 2018). However, King & Wheatley (2021) found that significant EUV irradiation of exoplanet atmospheres can occur on Gyr timescales, indicating that EUV-driven mass loss is not necessarily dominated by the high-energy activity early in the star’s lifetime, and can instead be important on Gyr timescales. Since different stars will have different time-integrated high energy fluxes, it can be difficult to compare mass loss histories for planets in different systems. However, in a multi-planet system all planets experienced the same irradiation history scaled by orbital distance, effectively removing this as a source of uncertainty, although the overall uncertainty on the system-wide scaling of XUV fluxes remains (Owen & Campos Estrada, 2020).

To date, there have been four published sets of mass loss measurements for multiple planets in the same system, including GJ 9827 b and d (Carleo et al., 2021), HD 63433 b and c (Zhang et al., 2022a), TOI-776 b and c Loyd et al. (2025), and TOI-836 b and c (Zhang et al., 2025). All of the planets in these systems have radii less than 3 R⊕, and therefore host relatively low-mass (a few percent of the total planet mass at most) envelopes with a wide range of possible compositions. Atmospheric escape observations in these systems have typically focused on whether or not the planet can retain its primordial atmosphere. These studies have found signatures of atmospheric escape from lower density mini-Neptunes ($1.5-3$ R⊕ Zhang et al., 2022a; Loyd et al., 2025; Zhang et al., 2025), and non-detections for super-Earths with higher bulk densities (radii $\lesssim 1.7$ R*⊕* Carleo et al., 2021; Zhang et al., 2022a, 2025). These observations provide additional support for the hypothesis that atmospheric escape can strip the primordial hydrogen/helium envelopes of mini-Neptunes, leaving behind their rocky cores (e.g., Lopez & Fortney, 2013; Owen & Wu, 2017).

In this study, we instead aim to use multi-planet systems to explore the origin of the Neptune desert. The TOI-4010 system (Kunimoto et al., 2023) provides an ideal test case for this purpose, as this system consists of three Neptune-sized planets orbiting a K dwarf host. All three planets are located near the lower edge of the Neptune desert (see Figure 1), with TOI-4010 b (3.02 $R_{\oplus}$, $P=1.35$ days) located inside the desert, TOI-4010 c (5.93 $R_{\oplus}$, $P=5.41$ days) just outside the desert in the ‘Neptune ridge’, and TOI-4010 d (6.18 $R_{\oplus}$, $P=14.71$ days) in the middle of the ‘Neptune savanna’ (for definitions of these regions see Castro-González et al., 2024, and also Figure 1). All three planets transit their stellar host, providing us with a unique opportunity to measure and compare mass loss rates from three planets located at different distances within the same system. The closely packed orbits of the three planets indicate that they cannot have migrated inward via high eccentricity channels, and are instead more analogous to the sub-Neptune-sized multi-planet systems along the lower edge of the desert. The densities of these three planets (mass uncertainties are $<$12%) show increasing envelope mass fractions with increasing distance (Kunimoto et al., 2023), and it is therefore particularly interesting to explore whether this gradient is primordial or could instead have been sculpted by mass loss processes.

We can characterize the present-day properties of the atmospheric outflows from the TOI-4010 planets by measuring the absorption signature in the metastable helium feature (He∗) when each planet passes in front of its host star (Spake et al., 2018; Oklopčić & Hirata, 2018). Measuring He∗ outflow signals allows us to constrain the potential effects of varying atmospheric composition (Salz et al., 2016; Crossfield et al., 2020; Dragomir et al., 2020; Zhang et al., 2022b), the presence or absence of a strong planetary magnetic field (Trammell et al., 2011; Owen & Adams, 2014; Schreyer et al., 2023), and stellar winds (Christie et al., 2016; Vidotto & Cleary, 2020; Carolan, 2021; Wang & Dai, 2021) on the observed outflow properties of the three planets. These mechanisms should result in distinct outflow morphologies, velocities, and overall mass loss rates (e.g., Wang & Dai, 2021; MacLeod & Oklopčić, 2022; Schreyer et al., 2023).

In this study, we use the Near Infrared Spectrograph (NIRSPEC) on Keck II Observatory to obtain He∗ transit observations of all three Neptune-sized planets in the TOI-4010 system. These observations allow us to expand our studies of atmospheric escape in multi-planet systems to larger (Neptune-sized) planets, which have massive hydrogen- and helium-dominated gas envelopes, for the first time. TOI-4010 is also the first three-planet system with mass loss constraints for all three planets. In Section 2, we detail our He∗ transit observations. In Section 3, we report the measured He∗ excess absorption signals from our Keck/NIRSPEC observations and in Section 4 we model the corresponding mass loss rates for each planet. Lastly, in Section 5 we discuss the effect of varying stellar high-energy flux and atmospheric metallicity on the excess absorption measurements.

2 Observations and Data Reduction

Using Keck/NIRSPEC, we observed full transits of TOI-4010 b on October 13, 2023 UT, TOI-4010 c on January 6, 2024 UT, and TOI-4010 d on December 31, 2023 UT. All observations were obtained using the high-resolution mode and $Y$-band filter. Observations of TOI-4010 b and TOI-4010 d were obtained using the [math]$\times$ [math] slit, with an exposure time of $300\text{\,}\mathrm{s}$ for TOI-4010 b and an increased exposure time of $600\text{\,}\mathrm{s}$ for TOI-4010 d due to poor weather conditions which resulted in an increased PSF width and correspondingly lower flux through the slit. The seeing was even worse during our observations of TOI-4010 c, and as a result we used a wider [math]$\times$ [math] slit with an exposure time of $600\text{\,}\mathrm{s}$ in order to maximize flux. For all observations we used an ABBA nodding pattern to subtract the background.

We analyzed the data using the pipeline described in Zhang et al. (2022b); Saidel et al. (2025). Briefly, we identified bad pixels and produced a master flat, which we used to construct A-B subtracted images for each AB pair. We then used optimal extraction to obtain the 1D spectrum and corresponding uncertainties, and wavelength calibrated this spectrum using a PHOENIX model spectrum (Husser et al., 2013) and model tellurics. Finally, we removed tellurics using molecfit (Smette et al., 2015) and interpolated all of the spectra onto a common wavelength grid.

Despite this correction, we found that the observations for all three planets were affected to varying degrees by residual telluric contamination (see Fig. 2). Our observations of TOI-4010 b contain an OH telluric feature redward of the He∗ line, which is reasonably well-separated in wavelength space. We have masked out this line ($1.083\,401.083\,48\text{\,}\mathrm{\SIUnitSymbolMicro m}$) as shown in the first panel in Figure 2. The location of the OH line shifts between observations, likely due to changes in location of the stellar helium line between observations due to Earth’s changing barycentric velocity, and we find that during our observations of TOI-4010 c and TOI-4010 d, the OH line now overlapped with the right side of the He∗ line. In order to ensure that the residual variations in this OH line do not affect our He∗ measurement, we have masked out the contaminated red wing of the line core ($1.083\,321.083\,39\text{\,}\mathrm{\SIUnitSymbolMicro m}$), as shown in Figure 2 and focused our analysis only on the blue wing of the He∗ line.

Next, we isolated the time-varying component of the spectroscopic timeseries by taking the log of all telluric-corrected spectra on the common wavelength grid and subtracting the average of every row and column. The resulting timeseries of residual spectra contain the fractional flux variations, where each pixel corresponds to a specific wavelength and time. We then constructed our in-transit excess absorption spectrum by taking all columns (corresponding to wavelength) of the residuals image and subtracting off the average out-of-transit residual (defined as the residuals preceding first contact and following fourth contact). Because the residuals are in log space, this is equivalent to dividing by the geometric mean. We removed continuum variations by masking out strong lines in the spectrum and fitting a third-order polynomial to the spectrum, then subtracting it from the residuals. We then computed the radial velocity of the planet relative to the star using each epoch’s barycentric Julian date and shifted the spectrum onto the planetary rest frame. The resulting excess absorption spectrum for each planet is shown in Figure 3.

3 Results

All three planets display no discernible excess He∗ absorption. To set an upper limit to the possible He∗ absorption for each planet, we take the mean excess absorption within a $0.75\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ bin centered on the main helium triplet peak ($10\,833.327\text{\,}\mathrm{\SIUnitSymbolAngstrom}$). Following e.g. Vissapragada et al. (2024), we estimate $\sigma$ as the rms scatter of the transmission spectrum binned to $0.75\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. For TOI-4010 b, we rule out absorption to $<$ 1.23% at 95% confidence, for TOI-4010 c and TOI-4010 d we rule out absorption to $<$ 0.81% and $<$ 0.87% respectively at 95% confidence.

4 Mass Loss Modeling

4.1 Joint Constraints on Outflow Temperature and Mass Loss Rate from Retrievals

Using our non-detections of He∗, we are able to place an upper limit on the mass loss rate for all three planets. We calculated this upper limit assuming an approximately solar ratio of hydrogen to helium (90% hydrogen, 10% helium by number, e.g., Aschwanden, 2014) in a one-dimensional isothermal Parker wind model. This is the most commonly utilized model in the literature (e.g., Oklopčić & Hirata, 2018; Dos Santos et al., 2022; Linssen et al., 2024), making it easier to compare our results to those of other published studies. We used p-winds (Dos Santos et al., 2022) to calculate the expected helium excess absorption spectrum as a function of the assumed mass loss rate ($\dot{M}$) and outflow temperature ($T_{0}$). When calculating these models, we adopted the planetary and stellar parameters from Kunimoto et al. (2023). TOI-4010 does not have a measured high-energy stellar spectrum, and we therefore used the MUSCLES spectrum of HD 85512 (Youngblood et al., 2017) as a proxy as the star is a close match in spectral type (both K stars), effective temperature (HD 85512: $4715\text{\,}\mathrm{K}$; TOI-4010: $4960\text{\,}\mathrm{K}$), surface gravity (HD 85512: $10^{4.39}$ cm/s2; TOI-4010: $10^{4.54}$ cm/s2; Youngblood et al., 2017; Kunimoto et al., 2023), stellar rotation period (HD 88512: $\approx$ 44-58 days; TOI-4010: 37.7 days; Lovis et al., 2011; Pepe et al., 2011; Laliotis et al., 2023; Kunimoto et al., 2023), and log $R^{\prime}_{\mathrm{HK}}$ (HD 85512: -4.98; TOI-4010: -4.90; Costes et al., 2021; Kunimoto et al., 2023) to TOI-4010.

We allowed $\dot{M}$ and $T_{0}$ to vary as free parameters and fit them to our data for each planet in order to obtain corresponding upper bounds on their outflow properties in $\dot{M}$-$T_{0}$ space. We defined prior ranges of $\mathcal{U}(6,12)$ for log $\dot{M}$ in $\text{\,}\mathrm{g}\,\mathrm{s}^{-1}$ and $\mathcal{U}(3900,15000)$ for $T_{0}$ in Kelvin, based on predicted isothermal outflow temperatures for sub-Neptune and Neptune-sized planets (e.g., Salz et al., 2016; Linssen et al., 2022) and used emcee to perform an MCMC analysis in order to retrieve the posterior probability distributions for log $\dot{M}$ and $T_{0}$. We verify that convergence is achieved by first confirming that the number of steps in our chain is greater than fifty times the autocorrelation length, indicating that we are running our chains for enough steps to generate independent samples. We confirm convergence by verifying that the Gelman-Rubin statistic $\hat{R}$ is less than 1.01 for all sampled parameters. The resulting constraints on the outflow properties of each planet as well as the predicted excess absorptions corresponding to a 50$\times$50 grid sampling our $\dot{M}$-$T_{0}$ prior range are shown in Figure 4.

We find a 95th-percentile upper limit on log $\dot{M}$ of $10^{10.17}$g s*-1* for TOI-4010 b, $10^{10.53}$g s*-1* for TOI-4010 c and $10^{10.50}$g s*-1* for TOI-4010 d. This is likely an over-estimate, as mass loss rates near these upper bounds would require unphysically high outflow temperatures (see discussion in §4.2). Nonetheless, if we assume a constant mass loss rate equal to the stated upper limit for each planet, we estimate atmospheric lifetimes ($M_{p}/\dot{M}$) of $>$ $110\text{\,}\mathrm{Gyr}$ for all three planets (we note that this is a lower limit on the planet’s lifetime as the star’s XUV flux should decrease with age). We conclude that even with this relatively conservative limit, photoevaporation will not significantly alter their atmospheric masses over the main sequence lifetime of the host star.

4.2 Comparison to Predictions from Forward Models

Next, we consider whether or not these upper bounds are consistent with predictions from self-consistent forward models for solar metallicity atmospheres. We calculated the predicted excess absorptions for each planet using pyTPCI (Rosener et al., 2025), a self-consistent one-dimensional radiative-hydrodynamics code to model atmospheric outflows. Following the methodology in Zhang et al. (2022b, 2025), we set a lower atmospheric boundary condition at each planet’s radii, defined a particle number density of $1\text{\times}{10}^{14}\text{\,}\mathrm{c}\mathrm{m}^{-3}$, and computed the corresponding atmospheric pressure using each planet’s equilibrium temperature. As in Zhang et al. (2025), we did not include molecules but did include atomic species with solar abundances greater than $10^{-5}$ (C, N, O, Ne, Mg, Si, S, and Fe), along with hydrogen and helium. We used each planet’s mass and semi-major axis ($M_{p}$, $a$), and the same stellar spectrum of HD 85512 (Youngblood et al., 2017) as a proxy for the unknown stellar spectrum of TOI-4010. We then modeled the predicted excess absorption signals and corresponding mass loss rates for each planet. In order to facilitate comparisons with the isothermal outflow models in our MCMC retrievals, we identified the maximum value of the radially varying outflow temperature and used this value when comparing the retrieved versus predicted outflow properties in Fig. 4.

We find that the predicted excess absorption signals for planets c and d should have been readily detectable in our observations, with a somewhat weaker predicted detection significance for planet b (see Fig. 3) likely due in part to its much smaller radius and higher model-predicted outflow temperature ($T_{0}$ TOI-4010 b: $8900\text{\,}\mathrm{K}$; $T_{0}$ TOI-4010 c: $5600\text{\,}\mathrm{K}$; $T_{0}$ TOI-4010 d: $5500\text{\,}\mathrm{K}$). The predicted outflow temperatures for all three planets lie near the lower end of the range considered in our retrievals. When we compare the predicted outflow properties in $\dot{M}$-$T_{0}$ space to the constraints from our p-winds retrievals, we find that our upper bounds on the mass loss rates for planets b and c are approximately an order of magnitude lower than the predicted mass loss rates from the pyTPCI models ($\dot{M}$ for TOI-4010 b: $10^{11.15}$g s*-1*; $\dot{M}$ for TOI-4010 c: $10^{10.98}$g s*-1*) with the same outflow temperatures (see Fig. 4). For planet d, our upper bound is consistent with the predicted mass loss rate of $10^{9.96}$g s*-1*.

5 Discussion

There are several factors that could potentially explain our non-detections, including a lower than expected XUV luminosity for the host star, confinement and/or suppression of the outflow by stellar winds or planetary magnetic fields, an enhanced atmospheric metallicity, or fractionation in the outflow. We investigate each of these explanations in turn in the following sections as well as the past mass loss histories of these planets.

5.1 Decreased Stellar XUV Flux

Photoevaporative outflows are powered by high-energy radiation from the star. If the true XUV flux of the host star was lower than that of our chosen stellar proxy (HD 85512), we would over-predict the mass loss rates and He∗ excess absorption signals for all three planets (we do note that stellar activity indicators such as the stellar rotation period and log $R^{\prime}_{\mathrm{HK}}$ are in good agreement for HD 85512 and TOI-4010 as discussed in Section 4.1, and the XUV flux of HD 85512 is in good agreement with the predicted present-day XUV flux of TOI-4010, see Section 5.5). LTT 9779 b is a useful case in point, as XMM-Newton observations revealed that the X-ray luminosity of LTT 9779 is 15$\times$ lower than expected for its age, drastically reducing its predicted mass loss rate (Fernández Fernández et al., 2024). Changes in the XUV flux can also result in a net change in the fractional abundance of He∗ in the outflow, as He∗ is populated via the ionization and subsequent recombination of helium atoms (Oklopčić & Hirata, 2018). This can additionally alter the predicted magnitude of the absorption signal.

We investigate the effect of a lower XUV flux on the predicted excess absorption signals by reducing the XUV flux from the proxy star by a factor of five, ten and fifteen relative to our original model. We find that decreasing the stellar XUV luminosity for planets b and c results in an initially steep decrease in the outflow temperature, which only decreases marginally with further XUV flux reductions (see Figure 4). Both planets see a similar, although less severe effect in the outflow rates, which undergo relatively modest reductions with decreasing XUV flux. These changes in outflow temperature and mass loss rate results in a modest increase in the excess absorption spectra of TOI-4010 b and minimal reductions in the absorption strength for TOI-4010 c (see Figure 3).

For planet d, the outflow rate and temperature only marginally decrease with decreasing XUV flux, however the excess absorption decreases significantly. This suggests that TOI-4010 d is likely in the energy-limited regime of mass loss, as a decrease in the energy deposition results in a decrease in the outflow rate and consequently the observed excess absorption (e.g., Murray-Clay et al., 2009; Owen & Alvarez, 2016; Lampón et al., 2021). However, for TOI-4010 d, even the lowest excess absorption (2.4%) is still readily detectable with Keck. Each planets’ predicted mass loss rates fall outside the distribution of mass loss rates and outflow temperatures allowed in our retrievals, and we conclude that a lower stellar XUV flux is unlikely to explain our non-detections.

5.2 Stellar Winds and Planetary Magnetic Fields

Strong stellar winds and planetary magnetic fields also play an important role in setting the observability of an outflow. Strong stellar winds can redirect an outflow into a cometary-like tail (e.g., Wang & Dai, 2021; MacLeod & Oklopčić, 2022) and in the most severe cases, drastically reduce the outflow to a fraction of its Hill sphere (Carolan et al., 2020). Published studies show that the effect of the stellar wind on the He∗ signal can vary. MacLeod & Oklopčić (2022) found that a strong stellar wind can lead to an enhancement of the He∗ line depth and extended post-egress absorption, as the confinement of the outflow creates an extended tail of high-density material. Alternatively, studies such as Wang & Dai (2021), find that outflow models with a stellar wind can at most lead to a reduction in the He∗ excess depth by a factor of two. Such a reduction would not reproduce the He∗ non-detections observed for the TOI-4010 planets, consistent with the similarly-massed GJ 436 system, in which strong stellar winds were unable to reproduce their He∗ non-detection (Rumenskikh et al., 2023). We therefore conclude that strong stellar winds do not provide a plausible explanation for our non-detections.

Alternatively, strong planetary magnetic fields can reduce the He∗ transit depth by forcing the outflow to follow magnetic field lines. Schreyer et al. (2023) showed that an increase in planetary magnetic field strengths from a $0.3\text{\,}\mathrm{G}$ to $10\text{\,}\mathrm{G}$ field can reduce the excess absorption at most by a factor of three, although we note that the excess absorption does not decrease monotonically with magnetic field strength. If all three planets have surface magnetic field strengths similar to Jupiter, such a reduction in signal size is plausible. However, if all three planets have magnetic field strengths representative of Uranus or Neptune, which are $<5\%$ the surface planetary magnetic field strength of Jupiter (e.g., Bagenal, 2013), the corresponding reduction in excess absorption strength is likely to be much less than a factor of three. In the event of a strong Jupiter-like planetary magnetic field, a reduction of this magnitude could only produce a non-detection for TOI-4010 b. Both TOI-4010 c and TOI-4010 d have much larger predicted excess absorption signals that should still be readily detectable even with significant magnetic confinement.

5.3 Enhanced Atmospheric Metallicity

Up to this point we have assumed a solar metallicity for the planetary outflows in our models. However, increasing the atmospheric metallicity to values greater than $100\times$ solar can drastically reduce the observed He∗ signature and overall mass loss rate (Zhang et al., 2022b; Linssen et al., 2024; Zhang et al., 2025). As the atmospheric metallicity increases, heating in the atmosphere becomes dominated by the photoionization of metals and metal lines become one of the dominant sources of cooling. This efficient cooling by metal lines reduces the outflow temperature and corresponding mass loss rate. The decreased mass loss rate also results in increased ionization of the outflow, which in turn increases the population of doubly ionized helium. This reduces the fractional He∗ population in the outflow, as the primary pathway for producing this state is the recombination of singly ionized helium with an electron. With a reduced population of singly ionized helium, the population of He∗ is also reduced (Zhang et al., 2022b, 2025). The number density of free electrons available for recombination is also limited at high metallicities. Assuming metal particles are not completely ionized (which is not necessarily true high in the atmosphere, but He∗ is formed quite close to the planet), a metal particle typically creates fewer free electrons per unit mass than hydrogen and helium (assuming a non-ionized metal which is not always the case in the upper atmosphere). At high metallicities, the increased metal content will therefore result in a reduction in the free electron density, limiting the amount of He∗ that can be produced via radiative recombination (Linssen et al., 2024).

We quantify the effect of increased atmospheric metallicity on the observed He∗ signal for the TOI-4010 system by using the same pyTPCI framework from Section 4.2 to model outflows with atmospheric metallicities of 10$\times$, 30$\times$, and 100$\times$ solar. Although we also attempted to model outflows with 200$\times$ solar metallicities, these models failed to converge at the onset of advection for all three planets. Following Zhang et al. (2025), we set the abundance of the metals (N, O, Ne, Mg, Si, S, Fe, C) equal to the solar abundance times the atmospheric metallicity. The resulting He∗ signals are convolved with the resolution of Keck/NIRSPEC and shown in Figure 5.

For all planets, increasing the atmospheric metallicity from solar to 10$\times$ and 30$\times$ solar results in an increase in the predicted He∗ signal. However, as the atmospheric metallicity is increased from 30$\times$ to 100$\times$ solar, the He∗ signal drops below solar predictions. This reduction is sufficient to explain our non-detection for TOI-4010 b, but still results in readily detectable He∗ signals for both TOI-4010 c (excess absorption: 4.9%) and TOI-4010 d (excess absorption: 4.4%). Zhang et al. (2025) found that an increase in atmospheric metallicity from 100$\times$ to 200$\times$ solar can result in an order of magnitude reduction in the predicted He∗ excess absorption, suggesting that the non-detections for planets c and d might be explained by a modest additional increase in atmospheric metallicity. TOI-4010 has a metallicity of [Fe/H] = $0.37\pm 0.07$ dex, making it one of the most metal-rich hosts of sub-Saturn planets (Kunimoto et al., 2023). It is therefore reasonable to expect that the atmospheric metallicities of the Neptune-sized planets in this system might be higher than those of Neptunes orbiting stars with lower metallicities.

5.4 H/He Fractionation

Mass-dependent fractionation of the outflow could also act to suppress the He∗ transit absorption signal. In this process heavier components, such as helium, decouple from the hydrogen in the uppermost layers of the atmosphere. As a consequence, the outflow can be preferentially depleted in helium, while the atmosphere becomes more enriched in helium over time (Zahnle & Kasting, 1986; Hunten et al., 1987). This process has occurred in the solar system planets (Lammer et al., 2008) and has likely influenced the transition of early Earth from a hydrogen-rich to a metal-rich atmosphere (Zahnle et al., 2019). It is therefore reasonable to investigate its potential impact on the TOI-4010 planets. Indeed, fractionation has already been invoked to predict the existence of He-rich sub-Neptune-sized exoplanets (Malsky et al., 2023; Cherubim et al., 2024). The methodologies used by the latter authors differ from the simulations we have presented up to this point, which assumed no fractionation.

The exact quantity of helium that is dragged into the outflow generally depends on the temperature profile, the mass-flux of hydrogen, the ionization state of helium, and the ‘injection radius’ of helium relative to the wind launching radius. Schulik & Owen (2024) have shown that for gas giant exoplanets, outflows can be significantly fractionated at semi-major axis distances greater than $\sim$0.05 au for K stars. TOI-4010 is also a K star and the outer two planets have semi-major axes within this range ($a=0.06$ and 0.11 au, respectively; Kunimoto et al., 2023), suggesting that fractionation may play a role in the observed outflows.

In order to estimate the expected depletion in helium for the TOI-4010 planets, we perform simple 1D simulations of this problem using the AIOLOS code (Schulik & Booth, 2023) including the dissociation of a molecular $\mathrm{H}_{2}$ layer and molecular cooling from the tables by Malygin et al. (2014). For our simulations, we set the lower boundary radius of the simulation to the $10\rm~mbar$ level, while injecting helium at a homopause pressure of $1\rm~nbar$. The latter value corresponds to the maximum altitude from the solar system planets (Atreya et al., 1991; Catling & Kasting, 2017). Overall, this setup allows us to compute the maximum expected helium abundance in a simulation in which the hydrogen is sourced from a molecular layer. In these simulations we measure the “reduction factor”, which describes how reduced the helium abundance is in the outflow. In particular, this reduction factor measures the ratio of the mass-loss rates at the sonic point ($\dot{M_{\mathrm{He}}}/\dot{M_{\mathrm{H}}}$, where H and He represent the sum of neutral and ionized hydrogen and helium respectively) normalized to the composition of the atmospheric base ($\rho_{\mathrm{He}}/\rho_{\mathrm{H}}$), which should be approximately representative of the bulk mass of the atmosphere. Normalizing the mass-loss rates results in a reduction factor that measures the composition of the outflow relative to the composition of the atmospheric base (approximately the composition of the bulk atmosphere). Reduction factors equal to one indicate that the helium composition in the outflow is not reduced (fractionation is not occurring) and therefore the compositions of the flow and atmospheric base are equal. For reduction factors less than one, fractionation is occurring, with the severity of fractionation increasing with a decreasing reduction factor (in other words more fractionation reduces the helium present in the outflow).

Our simulation results predict heavily fractionated outflows for planets c and d, with reduction factors of $10^{-2}$ for TOI-4010 c and $10^{-6}$ for TOI-4010 d. We additionally find that molecular cooling effectively shuts down the outflow for planet d, as shown in Figure 6. This effect is not seen in our previous pyTPCI simulations, as those models do not include molecular line cooling. TOI-4010 b is a special case, as its outflow lies in the photon-limited regime (Owen & Alvarez, 2016), where the stellar heating distributes itself over the entire radial extent of the flow. This leads to a rapid mass loss rate ($\sim 20$ km s*-1* at 10 $R_{p}$) and near-perfect coupling of helium into the outflow, resulting in a reduction factor of $1$. However, the helium in the outflow has insufficient time to ionize and remains mostly neutral, with the ion fraction reaching only $5\%$ at 10 $R_{p}$, reducing the strength of the predicted He∗ absorption signal during transit as presented in Fig. 6.

These simulations result in nearly complete reductions in the predicted excess absorption signals for TOI-4010 b, c, and d relative to the solar metallicity pyTPCI model predictions. These reductions are consistent with our non-detections, and we conclude that fractionation and/or molecular cooling can indeed suppress the detection of He∗ signals for this system. We note that the dominant source of cooling in the AIOLOS models is the molecular cooling and Ly$\alpha$ cooling for the atomic region, which differs from the dominant source of cooling (bremsstrahlung or free-free radiation cooling) in the pyTPCI models. These differences explain the difference in the predicted mass loss rates, as cooling mechanisms and the mean molecular weight control the outflow rate, which in turn determines the ionic state of the elements in the outflow. More detailed models that include both fractionation and tunable atmospheric metallicity are needed in order to study the competing effects of these variables on outflow strengths for these planets.

5.5 Constraints on Past Evaporation Histories

In the preceding subsections, we have evaluated possible mechanisms that would suppress signatures of present-day mass loss. In this subsection, we quantify what our non-detections of present-day mass loss mean for the past mass loss histories of all three planets. In order to study the past evaporation of these planets and place constraints on their internal structures, we performed simulations of their past evaporation histories taking into account the XUV emission history of their host star. To do so, we adopted the method described in Fernández Fernández et al. (2023), which combines three models: (1) a description of the XUV emission history of the host star, (2) a model of the internal structure of the planet (including the thermal evolution of its gaseous envelope), and (3) a model of atmospheric escape to compute the mass loss rate from the planet at each point in time.

To model the past XUV emission of TOI-4010, we adopted the rotational and XUV stellar evolution models by Johnstone et al. (2021). Rotational evolution models serve as a proxy for the XUV emission history as spin period and stellar X-ray activity are linked via the rotation-activity relation, where faster rotators are more X-ray active (e.g. Pizzolato et al., 2003; Wright et al., 2011). We find that the measured age and spin period of TOI-4010 (age: $6.1\pm 3.1$ Gyr; spin period: 37.7 days; Kunimoto et al., 2023) is consistent with the Johnstone et al. (2021) rotation models, which predict spin periods between 23 and 51 days at ages between 3 and 9 Gyr. However, we cannot determine whether the star was a fast or slow rotator when it was younger, and thus whether its earlier X-ray emission was bright or faint. This is because the diversity of spin periods that young stars exhibit is erased at later ages as stars spin down and converge into a single mass-dependent relation (Johnstone et al., 2015). For that reason, we considered both bright and faint XUV emission histories for the host star in our simulations of the planets’ evaporation histories – these correspond to the 5th and 95th percentile XUV emission models from Johnstone et al. (2021), which span the $2\sigma$ spread in rotation periods observed in stellar populations.

Next, to model the internal structures of the planets, we assumed a two-layer structure with a solid core surrounded by a H/He solar metallicity envelope. We note that our constraints on present-day mass loss rates are lower than the mass loss rates we predict for solar metallicity envelopes, indicating the atmospheric composition might deviate from solar (see Section 5.3 & 5.4). However, we use the assumption of solar metallicity envelopes as an upper bound on the past evaporation histories of these three planets. For the solid core, we adopted the empirical mass-radius relation for rocky cores by Otegi et al. (2019), and for the gaseous atmosphere we adopted the envelope structure model by Chen & Rogers (2016) based on MESA simulations. Moreover, we adopted the mass loss model by Kubyshkina et al. (2018), based on hydrodynamic simulations of escaping atmospheres; their model accounts for both X-ray and EUV driven escape as well as multiple regimes of mass loss, including energy-limited and recombination-limited photoevaporation as well as core-powered mass loss.

Finally, we adopted the photoevolver code (Fernández Fernández et al., 2023) to combine these models and simulate the evaporation past of the three planets. We used the RK45 algorithm as the integration method with a variable time step no larger than 1 Myr. Additionally, we explored a diversity of evaporation histories for each planet in order to account for the uncertainties in their radii, masses, and ages, as well as the unknown XUV emission history of the host star.

For TOI-4010 c and d, we evolved their current states backwards in time from the present day to an age of 10 Myr, corresponding to the estimated time when the protoplanetary disk has dispersed and boil-off processes have stopped (Rogers et al., 2024). Figure 7 shows the results of our simulations, which suggest that both TOI-4010 c and d have lost relatively little of their overall H/He envelope mass to atmospheric escape. Figure 7 also shows that the initial radii of both planets were 3 R⊕ larger than their present-day values at an age of 10 Myr; their radii decrease over time as a result of radiative cooling (Chen & Rogers, 2016).

For TOI-4010 b we considered initial gas envelopes with mass fractions of 1%, 2%, and 5%, and evolved them from the age of 10 Myr to the present day. We do not consider larger starting envelopes, as they would have been stripped by boil-off processes during the dispersal of the protoplanetary disk (e.g. Rogers et al., 2024). Our results, shown in Figure 8, suggest that all scenarios result in the complete loss of this initial gas envelope within the first 500 Myr, much earlier than the current age of the planet. Given that TOI-4010 b has a present-day bulk density consistent with a H/He envelope (Kunimoto et al., 2023), it is clear that its time-integrated mass loss rate must be significantly lower than predicted by these models. This is consistent with our previous conclusion that either a high atmospheric metallicity (at least 100$\times$ solar) or low He∗ abundance are needed to explain our non-detection of a He∗ absorption signal during the transit of TOI-4010 b. The former would result in a globally reduced mass loss rate, making it easier for TOI-4010 b to retain a hydrogen-rich envelope.

Using our model framework we also estimate present-day mass loss rates of all three planets. Using the present-day radii and masses for the TOI-4010 planets, along with the mass-loss model of Kubyshkina et al. (2018) and adopting a predicted XUV flux of $(3.0\pm 1.0)\times 10^{28}$ erg s${}^{-}1$ (uncertainties propagated from errors in stellar age and rotation period, we note that this estimate agrees with the XUV flux of HD 85512: $1.2\times 10^{28}$ erg s${}^{-}1$) from the rotation-activity models of Johnstone et al. (2021), we calculate present-day mass-loss rates (assuming the upper atmosphere is pure H/He) of $10^{11.00\pm 0.16}$ g s*-1*, $10^{10.65\pm 0.15}$ g s*-1* and $10^{9.73\pm 0.12}$ g s*-1* for TOI-4010 b, c and d respectively. The estimated present-day mass loss rate for TOI-4010 c is slightly higher than our 95% upper bound ($10^{10.53}$ g s*-1*) from Section 4.1, but still within $1\sigma$, whereas our estimated mass loss rate for TOI-4010 d is consistent with our upper bound ($10^{10.50}$ g s*-1*). However, the present-day mass loss rate for TOI-4010 b is much higher than our estimated 95% upper bound ($10^{10.17}$ g s*-1*). This discrepancy indicates that there must be some method of mass loss suppression not accounted for in these solar metallicity H/He models that could reduce the predicted mass-loss to a value consistent with our 95% upper bound and non-detection, such as an enhanced atmospheric metallicity (see Section 5.3) or a reduced He∗ abundance, although we note that the latter would only suppress the detectable signal and not the outflow itself (see Section 5.4).

6 Conclusion

In this work, we use Keck/NIRSPEC to search for evidence of He∗ absorption in the transmission spectra of all three planets in the TOI-4010 system, and find that none have detectable He∗ absorption signals. We place an upper limit on the possible He∗ absorption for each planet by averaging the excess absorption within $0.75\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ of the main helium triplet peak, which allows us to rule out absorption to 1.23%, 0.81%, and 0.87% for TOI-4010 b, c, and d respectively at 95% confidence. Using our transmission spectra and Parker wind models, we obtain upper limits on the possible mass loss rates for each planet of: $10^{10.17}$g s*-1*, $10^{10.53}$g s*-1*, and $10^{10.50}$g s*-1* for TOI-4010 b, c, and d respectively.

One-dimensional outflow models assuming a solar composition predict that all three planets in the TOI-4010 system should have easily detectable He∗ excess absorption. We attempt to explain our non-detections by considering four potential factors that could weaken or suppress the outflow: a decreased host star XUV luminosity, stellar winds or planetary magnetic fields, enhanced planetary atmospheric metallicities, and fractionation. We find that a decreased host star XUV luminosity, as well as strong stellar winds and planetary magnetic fields, are not individually sufficient to reduce the 1D outflow predictions to a level that is consistent with our non-detections. However, it is possible that a combination of these factors could explain our non-detections. This could be more rigorously quantified in future studies by developing mass loss models that incorporate tunable stellar wind strengths, XUV luminosities and planetary magnetic field strengths.

For TOI-4010 b, an atmospheric metallicity of 100$\times$ solar reduces the predicted He∗ signal to a level that is consistent with our non-detection. However, our upper limits for TOI-4010 c and d appear to require much larger metallicity enhancements ($\simeq 200\times$ solar) to suppress the outflow signals. Given the high stellar metallicity of TOI-4010, it is possible that the atmospheric metallicities of these planets could be enhanced relative to those of other Neptune-sized planets. This can be tested by upcoming JWST observations of the TOI-4010 system, which should be able to constrain the atmospheric metallicities of all three planets.

Alternatively, solar metallicity models including helium fractionation in the outflow and/or molecular cooling and other radiative effects can reproduce our non-detections for all three planets. It is also possible that both atmospheric metallicity and these other effects work in tandem to reduce the expected outflow strength on the TOI-4010 planets. Studies of hydrodynamic escape at non-solar metallicities show that molecular and atomic cooling can suppress the mass-loss rate (Yoshida & Kuramoto, 2020; Yoshida et al., 2022, 2024; Nakayama et al., 2022), with this suppression increasing with increasing atmospheric metallicity (Zhang et al., 2022b; Linssen et al., 2024; Zhang et al., 2025). This mass loss suppression will slow down the outflow, rendering it more prone to fractionation. Therefore, it is plausible that the combined effects of atmospheric metallicity and fractionation could also suppress the outflows on the TOI-4010 planets. Future mass loss models that incorporate the effects of tunable atmospheric metallicity and fractionation are therefore paramount to understanding the complex outflow properties and behavior on Neptune and sub-Neptune-sized worlds. We note that while our calculations show that it is plausible that fractionation suppresses the visible He∗ signal for the TOI-4010 system, it is currently unclear whether this applies to the population of He∗ detections and non-detections as a whole, see e.g., the compendia by Dos Santos (2023) and Orell-Miquel et al. (2024). Those authors particularly note a strong dependency on stellar type for the non-detections, with planets hosted by M stars having the largest frequency of non-detections. Future studies should explore whether fractionation can be applied universally to explain this trend.

We also evaluate the past evaporation histories for all three planets and find that TOI-4010 c and d have lost relatively little of their initial H/He envelope while TOI-4010 b should have been completely stripped of it’s envelope within the first $\sim$500 Myr. However, TOI-4010 b’s current bulk density is consistent with a puffy, H/He envelope, indicating a source of outflow suppression that must persist over the planet’s lifetime. Our modeling efforts indicate that such an outflow suppression can only be achieved by high atmospheric metallicities (100$\times$ solar), as a low He∗ abundance would only suppress the helium outflow signal and not the overall mass loss rate. The presence of a significant gas envelope on TOI-4010 b, as well as the low predicted percentage of past envelope loss on TOI-4010 c and d, place limits on the degree to which these planets may have been affected by photoevaporative mass loss. Our observations indicate that the higher envelope mass fraction of planet d relative to planet c is likely primordial, but our constraints for planet b are weaker and allow for a correspondingly wide range of initial envelope mass fractions.

acknowledgments

This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This research also made use of nep-des (available in https://github.com/castro-gzlz/nep-des). This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE‐1745301. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

We thank the Keck Observatory telescope and support astronomers, with special thanks to Greg Doppmann, Percy Gomez, Anthony Connors, Arina Rostopchina, and John Pelletier. We thank Justin Moore for his significant help in debugging software issues.

Co-author M. Zhang thanks the Heising-Simons Foundation for funding his 51 Pegasi b postdoctoral fellowship. Co-author J. Fernández Fernández acknowledges studentship support from the UK Science and Technology Facilities Council (STFC) and co-author P.J. Wheatley acknowledges STFC support under consolidated grants ST/T000406/1 and ST/X001121/1.