Strong Fractionation of Deuterium and Helium in Sub-Neptune Atmospheres along the Radius Valley

We simulate atmospheric fractionation in escaping planetary atmospheres using IsoFATE, a new open-source numerical model. We expand the parameter space studied previously to planets with tenuous atmospheres that exhibit the greatest helium and deuterium enhancement. We simulate the effects of extreme-ultraviolet-driven photoevaporation and core-powered mass loss on deuterium–hydrogen and helium–hydrogen fractionation of sub-Neptune atmospheres around G, K, and M stars. Our simulations predict prominent populations of deuterium- and helium-enhanced planets along the upper edge of the radius valley with mean equilibrium temperatures of ≈370 K and as low as 150 K across stellar types. We find that fractionation is mechanism dependent, so constraining He/H and D/H abundances in sub-Neptune atmospheres offers a unique strategy to investigate the origin of the radius valley around low-mass stars. Fractionation is also strongly dependent on retained atmospheric mass, offering a proxy for planetary surface pressure as well as a way to distinguish between desiccated enveloped terrestrials and water worlds. Deuterium-enhanced planets tend to be helium dominated and CH4 depleted, providing a promising strategy to observe HDO in the 3.7 μm window. We present a list of promising targets for observational follow-up.

Extreme ultraviolet (EUV)-driven photoevaporation and core-powered mass loss are two important proposed mechanisms for driving hydrodynamic escape of planetary atmospheres (Sekiya et al. 1980a;Watson et al. 1981;Ginzburg et al. 2018).Photoevaporation results from atmospheric absorption of high-energy radiation originating from the host star and subsequent hydrodynamic outflow in the form of a Parker wind, the bulk of which occurs in the first several hundred Myr postformation for most stellar types (Parker 1965).Corepowered mass loss is instead driven by a combination of host star bolometric luminosity and remnant planet formation energy and drives hydrodynamic outflow over Gyr timescales.Both mechanisms may have profound impacts on atmospheric composition in part due to mass fractionation of atmospheric constituents, including iso-topes (Sekiya et al. 1980b;Watson et al. 1981;Hunten et al. 1987;Zahnle et al. 1990;Yung & DeMore 2000).
Atmospheric fractionation is of particular interest for sub-Neptunes (i.e.R p ≲ 4.0R ⊕ ), the most common planetary class in the galaxy.Model predictions suggest that sub-Neptunes may become enriched in helium through escape (Hu et al. 2015;Malsky & Rogers 2020;Malsky et al. 2023).Transmission spectroscopy observations of escaping helium tails with supersolar H:He ratios on sub-Neptunes support these predictions (Zhang et al. 2023;Orell-Miquel et al. 2023).Deuteriumprotium (D/H) fractionation is of great interest as a tracer of planet formation and atmospheric evolution, and is potential observable with JWST and highdispersion ground-based instruments (Lincowski et al. 2019;Mollière & Snellen 2019;Morley et al. 2019;Kofman & Villanueva 2019;Gu & Chen 2023).Constraining D/H can also help assess planetary chemistry and habitability, since hydrogen loss is linked to both water loss and planetary oxidation (Genda & Ikoma 2008;Wordsworth et al. 2018;Wordsworth & Kreidberg 2022).
Here we explore the differences in deuterium and helium fractionation of sub-Neptune atmospheres between thermally-driven escape mechanisms -EUV-driven photoevaporation and core-powered mass loss -using our newly developed numerical model, IsoFATE: Isotopic Fractionation of ATmospheric Escape.Building on previous calculations in Hu et al. (2015), Malsky & Rogers (2020), Malsky et al. (2023), and Gu & Chen (2023), we expand the analysis across a wider parameter space, introduce core-powered mass loss effects, and simulate much higher levels of fractionation (D/H mole ratios exceeding Venusian D/H ≈ 1,000x Solar in helium-dominated atmospheres).We find that the two escape mechanisms predict helium and deuterium enhancement in distinct regions of planet mass-radiusinstellation space, offering a unique means to identify which mechanism dominates the evolution of sub-Neptune atmospheres.We describe the model in Sections 2 and 3 and present our results in Section 4. We discuss potential applicability to observational surveys in Section 5, and state our conclusions in Section 6.

ATMOSPHERIC ESCAPE MODEL
IsoFATE simulates mass fractionation via two atmospheric escape mechanisms: EUV-driven hydrodynamic escape and core-powered mass loss, described in the following sections.In both cases, planetary radius evolves temporally as a result of atmospheric loss and thermal contraction.The overall mass loss rate is governed by Equations 2, 3, 4, and 5, defined in the following sections.Radius evolution due to thermal contraction on  1).The dashed lines show the integrated flux density for 10 nm < λ < 130 nm from the semi-empirical HAZMAT spectra (Peacock et al. 2020).These spectra are computed for the lower quartile, median, and upper quartile EUV flux density samples of early M stars at five ages: 10 Myr, 45 Myr, 120 Myr, 650 Myr, and 5 Gyr.We chose FEUV,0 = 3500 W m −2 for our model, approximately equal to 0.1% of the bolometric stellar surface flux.
the other hand is governed by thermal evolution models presented in Lopez & Fortney (2014).This model component is derived from hydrostatic structure models for rocky planets enveloped in Solar-composition H/He atmospheres, assuming one third of the mass is in an iron core and two thirds is in a silicate mantle, which we assume too.

EUV-driven hydrodynamic escape
Our stellar radiation model assumes EUV (10 nm < λ < 130 nm) radiation from the host star is the primary driver of mass loss for the EUV-driven hydrodynamic escape mechanism.We ignore X-ray irradiation since approximately 80-95% of the high energy flux is in the EUV (Peacock et al. 2020).This also allows us to remain consistent with available semi-empirical spectra, as discussed below (Loyd et al. 2016;Fontenla et al. 2016;Peacock et al. 2019aPeacock et al. ,b, 2020;;Richey-Yowell et al. 2023).We assume a H/He composition with solar metallicity (mean molecular mass = 2.35 u) in the bulk atmosphere and we assume that hydrogen in the escaping wind is atomic due to photodissociation of H 2 catalyzed by water vapor (Liang et al. 2003;Yelle 2004;Moses et al. 2011;Hu et al. 2012).
To model the temporal evolution of the incident EUV flux, F EUV , for planets around M stars, we constructed a power law function as in Cherubim et al. (2023), Cloutier et al. (2023), andCadieux et al. (2023).Our model is informed by semi-empirical spectra generated by the HAZMAT team for populations of M1 stars over a range of ages from 10 Myr to 5 Gyr (Figure 1; Peacock et al. 2020): where t sat is the saturation time marking the transition from constant EUV flux to power law decay.Note that F EUV is related to stellar EUV luminosity L EUV by: F EUV = L EUV /4πa 2 , where a is the orbital semimajor axis.We found t sat = 500 Myr to be consistent with the HAZMAT spectra and β = −1.23 to be consistent with the HAZMAT spectra as well as previously reported Solar data (Ribas et al. 2005).We set F EUV,0 equal to 0.1% of the incident planetary flux, F p , in our model simulations to remain consistent with HAZMAT data.For planets around K dwarfs, we adopt t sat = 200 Myr, consistent with high energy spectra for stars of various ages in the HAZMAT survey Richey-Yowell et al. (2023).For planets around G dwarfs, we adopt t sat = 100 Myr (Ribas et al. 2005;Garcés et al. 2011).For planets around both K and G dwarfs, we adopt F EUV,0 = 10 −3.5 F p (M star /M ⊙ ) based on Jackson et al. (2012) and commonly used in photoevaporation studies (e.g.Owen & Wu 2017;Lopez & Rice 2018;Rogers & Owen 2021).
In the energy-limited escape regime, EUV-driven escape generates a mass flux [kg m −2 s −1 ], which we compute as a function of F EUV , an efficiency factor ϵ, and the planetary gravitational potential V pot = GM p /R 2 p , where M p and R p are planetary mass and radius respectively: ϵ encapsulates several heat transfer processes, ultimately representing the fraction of incident radiation that drives escape.We chose a value of ϵ = 0.15, consistent with previous photoevaporation studies and model estimates (Watson et al. 1981;Owen & Wu 2013;Shematovich et al. 2014;Schaefer et al. 2016;Kubyshkina et al. 2018).As discussed later in Section 4.2, various cooling mechanisms may lead to ϵ values below 0.15, one of which we model, as shown in the next paragraph.It may be that radiative cooling becomes more important for higher metallicity atmospheres as fractionation takes place, which should allow more planets to retain fractionated atmospheres.Simulations with ϵ = 0.05 -0.10 showed a slight shift in the radius valley, and an accompanying shift in the planet radius-flux space, but still a prominent population of fractionated atmospheres.We leave more detailed modeling of escape efficiency for future work.
A reduction in escape efficiency is expected due to radiative cooling, especially for closer-in planets whose winds are thermostatted at ≈ 10,000 K at the τ = 1 surface so that additional EUV flux does not increase the escape rate (Murray-Clay et al. 2009).In this case, atmospheric gas is sufficiently ionized such that hydrogen recombination and emission of Lyα photons becomes the dominant cooling mechanism.Planets limited by this cooling mechanism rather than by F EUV are in the "radiative/recombinationlimited" escape regime (Lopez & Rice 2018;Lopez 2017;Wordsworth et al. 2018).Mass flux [kg m −2 s −1 ] in the radiative/recombination-limited escape regime is computed as: where c s = 2k B 10 4 /µ H is the sound speed at the sonic point [m/s], g = GM p /R 2 p is the gravitational field strength at the base of the flow [m/s 2 ], R p is taken as the minimum of the Bondi radius, Hill radius, and Lopez & Fortney (2014) radius prescription [m], h is the Planck constant, ν 0 = 4.835 x 10 15 is the EUV ionizing radiation frequency [Hz; ≈60 nm/ 20 eV], α rec = 2.7 x 10 −13 (T eq /10 4 ) −0.9 /10 6 is the case B recombination coefficient for hydrogen [m 3 /atom/s], k B is the Boltzmann constant, T eq is the planetary equilibrium temperature assuming zero albedo, R s = GM p /(2c 2 s ) is the sonic point radius [m], and M p is the planet mass [kg].For all EUV-driven hydrodynamic escape simulations performed with IsoFATE, the mass flux is set by the minimum of ϕ EUV (Equation 2) and ϕ RR (Equation 3).

Core-powered mass loss
Our core-powered mass loss model mirrors that of EUV-driven hydrodynamic escape except that mass loss is powered by the intrinsic planetary luminosity L p , which results from left over heat of formation rather than by stellar EUV radiation.We follow the prescription outlined in Gupta & Schlichting (2020), which assumes the escaping atmosphere is isothermal at T eq .In the energy-limited escape regime, the mass flux [kg m −2 s −1 ] expression for the core-powered mass loss mechanism is: where ϵ is a heat transfer efficiency factor and V pot is the gravitational potential, as in Equation 2. Note that our formulation is equivalent to Equation 9in Gupta & Schlichting (2020) with the exception of the added heating efficiency term ϵ, and L p is equivalent to their Equation 6.This formulation corresponds to an absolute upper limit on escape rate, assuming all cooling luminosity goes into driving escape.
In the "Bondi-limited" escape regime, on the other hand, the escape rate is determined by the thermal velocity of the escaping particles at the Bondi radius.For this case, we again follow the prescription outlined in Gupta & Schlichting (2020) (Equation 10): where ρ rcb = 1.0 kg/m 3 is the density at the radiativeconvective boundary and R rcb is the height of the radiative-convective boundary, calculated as R core +R env following Lopez & Fortney (2014) (i.e.tropopause height).Equation 5implicitly assumes that the atmosphere is isothermal in the escaping region.For all corepowered mass loss simulations performed with IsoFATE, the mass flux is set by the minimum of ϕ L (Equation 4) and ϕ B (Equation 5).

ATMOSPHERIC MASS FRACTIONATION
For a hydrogen-dominated planetary atmosphere undergoing hydrodynamic escape, the primary escaping species is hydrogen.The remaining composition of the escaping wind depends on the upward flux of minor species present below the exobase.If the total escape flux is low, H escapes alone.When the escape flux exceeds a critical value, heavier species are dragged along with the escaping flow.Hence, the loss of heavier species is governed by the competition between upward drag by escaping H and downward diffusion by gravity.
Our model simulates this process by numerically integrating via the forward Euler method several differential equations of the general form where N i is the total number of moles of species i, Φ i is the number flux of species i [particles m −2 s −1 ] defined in Equations 7, 8, and 11, and A is the planetary surface area [m 2 ].
We consider diffusive mass fractionation of escaping atmospheres composed of hydrogen and helium present in Solar abundances (Lodders 2003).We consider two cases: helium-hydrogen fractionation (He/H) and deuterium-hydrogen fractionation (D/H).For He/H fractionation, we assume a binary mixture where H is the light species and He is the heavy species.We do not expect deuterium to impact this system since it is never the dominant species in the atmosphere.However, for D/H fractionation, we model a ternary mixture of H, D, and He.
For a binary mixture, we follow the prescription for diffusive fractionation derived in Wordsworth et al. (2018) (see also Zahnle & Kasting 1986;Zahnle et al. 1990;Hu et al. 2015).The number flux [particles m −2 s −1 ] of a light species Φ 1 and a heavy species Φ 2 can be calculated as a function of their molar concentrations in the bulk atmosphere (x i = N i /N tot ), their particle masses m i , and the total mass flux (from Equation 2or 4; simplified here as ϕ): where m = m 1 x 1 +m 2 x 2 is the mean particle mass of the escaping flow and b 1,2 is the binary diffusion coefficient for species 1 and 2 and H i is the effective scale height of species i at the base of the escaping region.We use b H,D = 7.183 × 10 19 T 0.728 , b H,He = 1.04 × 10 20 T 0.732 , and b He,D = 5.087 × 10 19 T 0.728 (Mason & Marrero 1970;Genda & Ikoma 2008).The quantity ϕ c is the critical mass flux required for the light species to drag the heavier species along with the flow.It is defined as: This results from setting the two definitions of Φ 2 in Equation 8 equal to each other, setting ϕ = ϕ c , and using the scale height definition H i = k B T /m i g where g is gravitational field strength, k B is the Boltzmann constant, and T is temperature.
Helium enrichment can have an important effect on deuterium fractionation, so we model a ternary mixture of H, D, and He for all D/H fractionation simulations.In this case, we model He/H fractionation as usual with Equations 7 and 8 where species 1 is H and species 2 is He.We derive our D escape flux from an analytical expression for an arbitrary number of species escaping in a subsonic wind from Zahnle et al. (1990): where f j = N j /N 1 is the mixing ratio of species j, N 1 is moles of the primary escaping species (i.e., H), r is radius, r 0 is the planet radius, M is the planet mass, m is atomic mass, Φ is escape flux as previously defined, and b i,j is the binary diffusion coefficient for species i and j.
After setting the total number of species to 3, expanding out sums and solving for Φ 3 , we arrive at an expression for the D (species 3) escape flux in a H/He background gas: where ).An equivalent expression with a more detailed derivation is reported by Gu & Chen (2023).
Figure 2 shows a comparison between our numerical model and an analytic solution (derived in Appendix A) for a representative simulated planet with M p = 12 M ⊕ , orbital period P = 10 days, and initial atmospheric mass fraction f atm,0 = 1.5%.Unlike our numerical results presented in the following sections, planet radius and escape flux were held constant in this case for simplicity.These results demonstrate that planetary atmospheres experience the most mass fractionation -be it D/H or He/H -toward the end of the escape process, i.e. when nearly the entire atmosphere is lost.Hence, we expect planets that are able to retain a thin gaseous envelope to exhibit the greatest He-and D-enhancement.However, we find that He-and D-enhancement is not solely a function of final atmospheric mass; the escape history is important too.If the escape flux is too high (ϕ >> ϕ c ), lighter species are dragged along with heavier species and fractionation is minimal.Given the canonical interpretation of the radius valley as separating enveloped terrestrial planets from bare rocky cores, we therefore expect planets along the upper radius valley at sufficient orbital separations to exhibit the most fractionation.We validate this prediction in the following section.

MODEL SIMULATIONS
We are interested in understanding the thermallydriven atmospheric escape histories for sub-Neptunes around G, K, and M stars and how their final atmospheric compositions depend on their initial conditions.To explore these questions, we ran a suite of atmospheric evolution models via Monte Carlo simulations over a broad parameter space.For each trial, 500,000 samples were randomly drawn from log uniform grids of initial M p between 1 -20 M ⊕ , initial f atm between 0.1 -50%, and P between 1 -300 days.Atmospheric mass loss was initiated at 1 Myr.Orbital migration effects were ignored so P was held constant.Model simulations were halted at 5 Gyr.We assume Solar composition for initial D/H and He/H values according to Lodders (2003), though in reality these values are expected to deviate for different systems.An average trial of 500,000 simulated planets runs for approximately 36 hours, affording rapid exploration of a broad parameter space to yield population-level predictions.We focus on the results for planets around low-mass stars, given their greater observability.We also choose to simulate escape with EUV-driven photoevaporation alone, core-powered mass loss alone, and the combination of the two.There is a great body of empirical evidence (Section 1) to suggest that these mechanisms may be sculpting the radius distribution around G, K, and M stars.More recent work suggests that each may operate independently or in tandem for sub-Neptunes (Owen & Schlichting 2023).
In order to test our prediction of fractionation along the upper radius valley and to compare our model to previous studies, we determined the location of the radius valley in our simulations by fitting a linear model to logR p vs. logP following the approach taken by Van Eylen et al. (2018).This approach uses support vector machines to determine the hyperplane of maximum separation between planets above and below the radius valley, which is a line in this case.For this purpose, we classify planets below the valley as "super-Earths," defined as having lost their entire atmospheres, while planets above the valley are classified as "sub-Neptunes" and have retained a H/He envelope.The line of separation maximizes the distance to points of these two predefined classes of simulated data.To obtain the model fits, we employed the Python support vector classification routine in the scikit machine learning package scikit-learn.This method requires choice of a penalty parameter C, which sets a tolerance for misclassification, with high values allowing lower misclassification.We tested values between C = 0.1 -100 and determined that C = 100 provided the best visual fit to the data and best agreed with previously reported values.Larger values of C resulted in equally good fits.

Helium Enhancement
Following Hu et al. (2015), Malsky & Rogers (2020), and Malsky et al. (2023), we first explored helium enhancement of sub-Neptune atmospheres through hydrodynamic escape.Figure 3 shows a set of simulations for planets around an early M dwarf experiencing EUVdriven photoevaporation.As expected, significant helium enhancement is seen for planets along the upper edge of the radius valley.For planets with high fractionation, there is a negative correlation between orbital period and planetary mass, which is further highlighted in Figure 4. Figure 3 also shows reasonable agreement of the radius valley location with Lopez & Rice (2018) and Owen & Wu (2017).Note that only the radius valley slope is reported in these works and the vertical offset is approximated by eye here.
In order to maximize fractionation for a given planetary atmosphere, the condition ϕ ≈ ϕ c must be met so that a sufficient quantity of the lighter species is driven away and a sufficient inventory of the heavier species is retained.For ϕ >> ϕ c , the heavier species fails to diffuse downward through the escaping flow rapidly enough and is dragged along with the lighter escaping species.Hotter, more irradiated planets closer to their host stars must possess stronger gravitational fields and hence be more massive in order to undergo fractionation and avoid excessive escape fluxes.This critical planetary mass decreases with increasing orbital period, hence the negative correlation seen in Figures 3 and 4. The correlation is stronger for EUV-driven photoevaporation since it is a stronger driver of escape than corepowered mass loss (hereafter "CPML"), leading to more massive fractionated planets as seen in Figure 4.
The top panel of Figure 5 shows helium molar concentrations for EUV-driven photoevaporation (left), CPML (middle), and the combined mechanisms (right) (Compare to Malsky & Rogers (2020) Figure 5 and Malsky et al. ( 2023) Figure 1).Qualitatively, our results resemble those of Malsky et al., with significant helium enrichment observed along the upper edge of the radius valley.However our results differ in that the area near the radius valley is more populated with fractionated planets, whereas the results of Malksy et al. show significant gaps in the distribution near the radius valley, presumably due to the reported failure of the MESA code's equations of state for planets with tenuous atmospheres.Another difference is that the most fractionated planets appear at greater orbital distances in our results while they appear at intermediate orbital distances in the results of Malsky et al.Finally, while Malsky et al. report planets with atmospheres ≳ 80% helium by mass, our simulations yield planets with greater overall helium enhancement, reaching ≈ 100%.
Figure 5 was generated by binning data from simulations represented in Figure 3.The mean helium molar concentration is computed for each grid cell only for planets that have retained an atmosphere at the end of the simulation.We observe closer-in helium-enhanced planets (P ≲ 10 days) have greater planetary radii in the EUV-driven photoevaporation simulations compared to those in the CPML simulations.This follows from Figure 4 which shows that these planets are more massive and possess greater final atmospheric mass fractions in the EUV-driven case compared to the CPML case.This is also reflected in the steeper radius valley reported in previous photoevaporation simulations  2023) report similar findings for EUV-driven photoevaporative mass loss, with planets in their simulations achieving He mass fractions ≥ 0.40 after 2.5 Gyr and ≥ 0.80 after 10 Gyr.The different parameter spaces occupied by helium-enhanced planets observed in this study offers a novel approach to investigate whether EUV-driven photoevaporation or CPML dominates sub-Neptune atmosphere evolution around low-mass stars.We discuss observational prospects to this end in Section 5.

Deuterium Enhancement
We next explore the atmospheric fractionation of deuterium and hydrogen (i.e., protium) as a result of hydrodynamic escape.Previous work by Gu & Chen (2023) report modest D/H fractionation for sub-Neptunes around Solar-type stars (max N D /N H ≈ 1.7x Solar).Like the He fractionation studies of Malsky & Rogers (2020) and Malsky et al. (2023), Gu & Chen (2023) utilized the MESA code and encountered anomalies in the data for small planet radii.This shortcoming prevented the smooth transition to a bare rocky core as a planet lost its atmosphere.As we expect from analytic solutions (e.g. Figure 2) and as we see from our simulations, planetary atmospheres are significantly fractionated in the final stages of losing their atmospheres.The loss process is imprinted on the planet only if the planet manages to retain a gaseous envelope.Hence, previous work that employed the MESA code to compute equations of state was unable to probe the parameter space in which planets become deuterium-enhanced.
Here we build on these previous studies and explore a wider parameter space, to provide a more complete picture of deuterium fractionation.The bottom panel of   2023)).Note that the D/H mole fraction is capped at 100x Solar in Figure 5 and many planets retain atmospheres with essentially no hydrogen remaining, only deuterium.Though qualitatively similar to past work, our results demonstrate that deuterium enhancement was previously dramatically underestimated.We observe significant fractionation along the upper radius valley for all three mechanism combinations in our simulated data.Not surprisingly, we observe deuteriumenriched planets in the same parameter space as that of helium-enriched planets.In fact, we find that 100% of deuterium-enriched planets are also helium enriched.
Our results suggest that atmospheric D/H may serve as a useful diagnostic of surface pressure for sub-Neptunes with thin hydrogen/helium envelopes.Figure 6 shows deuterium concentration as a function of final atmospheric mass fraction, f ret , for planets around an M dwarf experiencing EUV photoevaporation.We observe a sharp transition between f ret = 10 −4 -10 −3.5 where planets essentially lose all of their hydrogen, and only trace deuterium in a thin, helium-dominated atmosphere remains (surface pressure ≲ 1,000 bar).Differ-ences in results for other stellar types were negligible, though the curve is shifted to lower f ret values by about 0.5 dex for the CPML mechanism.As we have demonstrated, the most tenuous atmospheres tend to be the most fractionated.Hence, atmospheric D/H can serve as a proxy for the upper limit of surface pressure.Misener & Schlichting (2021) propose that the core-powered mass loss mechanism may shut down when the cooling timescale equals the loss timescale, leaving behind a tenuous atmosphere.Other factors such as molecular line cooling may decrease escape efficiency, similarly resulting in retention of tenuous atmospheres (Nakayama et al. 2022).Misener & Schlichting (2021) report analytical estimates of f ret ∈ [10 −4 , 10 −8 ] which overlap with values corresponding to elevated D/H in our simulations (see their Figure 5).The f ret values of 10 −3 , 10 −4 , and 10 −5 correspond to mean surface pressures, P surface , of approximately 10,000, 700, and 60 bar respectively in our simulations, as seen in Figure 6.These values are calculated as the mean P surface for all simulated planets within f ret ±0.5 dex, where P surface = GM 2 p f ret /4πR 4 core .Since D-enriched atmospheres tend to be He-dominated, a similar relationship is observed for f ret and x He suggesting He/H can be used as a proxy for P surface too.
While constraining atmospheric D/H presents an observational challenge, depletion of hydrogen in deuterium-enhanced, helium-dominated atmospheres may confer an advantage.Mollière & Snellen (2019) show that HDO may be detectable in sub-Neptune atmospheres with D/H as low as the VSMOW value, so long as methane (CH 4 ) is depleted to leave the 3.7 µm window open.Fortuitously, deuterium-enhanced planets are typically helium-dominated and are expected to be depleted in CH 4 .In a solar composition atmosphere ≲ 700 K, the dominant carbon-bearing species is CH 4 (Fortney et al. 2020).For an evolved atmosphere, CH 4 is depleted in favor of CO and CO 2 when x H < x C + x O (Hu et al. 2015), a condition commonly achieved in our simulations using He as a proxy for heavier species.Therefore, a residual helium atmosphere will have most of its carbon inventory in CO and CO 2 rather than CH 4 , making HDO more observable and D/H measurement more feasible.
Our calculations assume that the planet forms with a rocky core and a solar-composition atmosphere.Water world exoplanets, which may have significant fractions of their mass as H 2 O, would likely exhibit much less D/H fractionation due to rapid H and D exchange between H 2 O reservoirs and atmospheric H 2 , unless they also lost all of their H 2 O inventory (Genda & Ikoma 2008).Hence D/H observations should also provide a way to distinguish water worlds from fractionated sub-Neptunes.To identify planet candidates with helium-and deuterium-enriched atmospheres, we queried the NASA Exoplanet Archive on August 24, 2023 for sub-Neptunes in the relevant parameter space.We filtered the sub-Neptunes for those with mass and radius measurement precision of at least 55% and 8% respectively and massradius profiles inconsistent with purely rocky compositions.We then filtered our simulated data for those planets with helium molar concentration of x He ≥ 0.145 (equivalent to Malsky et al. (2023) cut off at mass fraction = 0.4) and deuterium mole fraction of N D /N H ≥ 100x Solar, motivated by deuterium isotopologue observability predictions with JWST and high resolution ground-based spectroscopy for terrestrial exoplanets (Lincowski et al. 2019;Mollière & Snellen 2019).Using our simulated data, we created a three-dimensional grid of planet mass, planet radius, and planet incident flux.The mean planet flux was calculated for all fractionated planets in each mass-radius grid cell (Figure 7).Finally, we used the RegularGridInterpolator routine in the SciPy Python library to interpolate the flux at the mass-radius location of our planet candidates.We have demonstrated that atmospheric fractionation is largely a function of planet mass, planet radius, and stellar irradiation history.Of these parameters, irradiation history is the most uncertain for the population of exoplanets with measured masses and radii.Our analysis ignores orbital migration, short-term stellar activity, planet albedo, and other factors affecting irradiation history.Hence we consider a range of flux values in which our planet candidates could fall to be considered an observational prospect.Following from Malsky et al. (2023), we compared the interpolated flux to the flux calculated for each planet.If the interpolated flux fell within 0.1x and 10x the calculated flux, then we included the planet in our list.The results are presented in Figure 7 and Table 1.The rightmost columns in Table 1 show the mechanisms for which He/H and D/H fractionation are predicted along with an estimate of the probability of fractionation for each target, P frac .This probability is the averaged number of fractionated planets divided by the total number of planets in each grid cell occupied by the target and its error bars in M p -R p -F p space (Figure 7).The transmission spectroscopy metric (TSM) was calculated following the prescription outlined in Kempton et al. (2018).

OBSERVATIONAL PROSPECTS
The parameter space populated by simulated, fractionated planets is visibly different between mechanisms.That of the CPML simulations particularly stands out.As seen in Figure 4, planets undergoing CPML alone have lower masses than their counterparts undergoing EUV photoevaporation.These planets also have lower f ret values (≈ 0.02%).A narrow distribution of f ret values will produce a narrow distribution of planetary radii.Hence we see a narrower band of fractionated planets in M p -R p space for the CPML scenario.We see wider bands for planets undergoing EUV photoevaporation because the additional parameter of variable EUV flux over time results in a wider distribution of allowed planet masses and initial atmospheric mass fractions that result in residual fractionated atmospheres.We also note that in the EUV scenarios, the maximum planetary flux is lower than that in the CPML scenario for a given planetary mass because atmospheres of closein planets cannot survive the intense irradiation which results in complete stripping of the atmosphere.
Since the M p -R p -F p parameter space is different for each mechanism combination, the planets predicted to have fractionated atmospheres differ by mechanism as well (See Table 1).Planets around K dwarfs are particularly interesting since they are especially amenable to helium detection via the 1083 nm line (Oklopčić 2019).Three targets predicted to have helium atmospheres stand out among this group: Kepler-80 e, TOI-836 b, and TOI-178 c.These three targets are predicted to have helium enhancement in the CPML scenario, but not when EUV photoevaporation is included.These planets offer a unique means to test which thermallydriven escape mechanism dominates the sculpting of planets around K dwarfs.Of these targets, TOI-836 b and TOI-178 c have been allocated time on JWST for atmospheric characterization (Batalha et al. 2021;Hooton et al. 2021).Kepler-80 e has a TSM of 6, and is likely unobservable with current facilities.
The remaining targets around K dwarfs and those predicted to be fractionated around M dwarfs offer a means to test if thermally-driven escape mechanisms work to sculpt the radius valley versus alternative hypotheses such as gas-depleted formation (Lee et al. 2014;Lopez & Rice 2018;Lee & Connors 2021).Recent work has demonstrated that the radius valley for planets around low mass stars may result from a unique channel of planet formation entirely different than the thermallydriven escape processes evidenced to be sculpting the radius valley of planets around Sun-like stars (Cloutier & Menou 2020;Cherubim et al. 2023).Here we propose a new strategy for investigating this open question through constraining He/H and D/H ratios of planets around low mass stars.
6. CONCLUSION 1. Helium-dominated, deuterium-rich atmospheres are predicted for sub-Neptunes along the upper edge of the radius valley for planets around G, K, and M stars.This novel class of planets spans a large range of T eq ∈ [150 K, 890 K], with a mean of 370 K, largely independent of stellar type, placing some planets in their habitable zones.Closer-in planets require a greater mass to retain fractionated atmospheres in the face of enhanced escape.This is especially true for EUV photoevaporation.
We show that EUV-driven photoevaporation fractionates planets in a wider planet/atmospheric mass distribution relative to CPML, endowing them with a wider distribution of final atmospheric mass fraction (f ret,EUV ≈ 0.06%, f ret,CPML ≈ 0.03%).More massive planets tend to become fractionated via EUV-driven photoevaporation compared to CPML.2023) that preclude resolution of atmospheric composition for planets with the most tenuous atmospheres, which are the ones that undergo the most fractionation.Hence we expand the parameter space for helium and deuterium enhancement predictions and demonstrate that fractionation was previously underestimated.

Our numerical model
4. Constraining atmospheric He/H and D/H abundance for sub-Neptunes offers a unique means of investigating competing mechanisms that may explain the origin of the radius valley for low mass stars.Observational studies are needed to test existing model predictions.We leave more detailed modeling of processes that may confound the fractionation/escape relationship, such as cometary delivery and interior-atmosphere exchange, for future studies.
5. Helium-dominated atmospheres are expected to be depleted in CH 4 and instead possess CO 2 (Hu et al. 2015).This leaves the 3.We present an analytic solution for D/H in a ternary mixture of H, D, and He to better understand the sharp transition to enhanced D/H in helium-dominated atmospheres shown in Figure 6.For simplicity, we allow only D and H to escape from a static He background.The diffusion-limited escape flux for species i in a He-dominated atmosphere is where x i = N i /(N He + N i ) is molar concentration and x i ≈ N i /N He for a minor species i. b is the binary diffusion coefficient for He and species i. Combining with Equation A1, we get Expanding the scale height terms gives us where m p is the proton mass and m He ≈ 4m p .This simplifies to for H, where τ f is the fractionation timescale: The total number of moles of He in the atmosphere is Substituting this into τ f , we get which simplifies to: Integrating Equation B13 and a similar expression for D, d log N D /dt = −2/τ f , the number of moles of H and D at time t is To solve for D/H at time t, we simply divide these expressions: The result is a classic expression for Rayleigh fractionation.

Figure 1 .
Figure 1.The solid red line shows the stellar surface EUV flux evolution model used in our photoevaporation simulations for planets around M stars (Section 2.1, Equation1).The dashed lines show the integrated flux density for 10 nm < λ < 130 nm from the semi-empirical HAZMAT spectra(Peacock et al. 2020).These spectra are computed for the lower quartile, median, and upper quartile EUV flux density samples of early M stars at five ages: 10 Myr, 45 Myr, 120 Myr, 650 Myr, and 5 Gyr.We chose FEUV,0 = 3500 W m −2 for our model, approximately equal to 0.1% of the bolometric stellar surface flux.

Figure 2 .
Figure 2. Comparison of our numerical results with an analytic solution.The top panel shows the time evolution of total D and H inventory in moles.The analytic solution for NH is represented by the blue dashed line: A is the planetary surface area, ΦH is the H escape flux and t is time.The bottom panel shows the time evolution of the mole fraction of D/H.The analytic solution for ND/NH is represented by the dashed orange line: x3 is D mole fraction and τ = NH,0/(AΦH) is the atmospheric loss timescale.The derivation is shown in Appendix A, Equation A9.

Figure 3 .
Figure 3. EUV-driven photoevaporation simulations for a random subset of n = 1 × 10 4 planets around an early M dwarf.Helium molar concentration is displayed as a function of orbital period and final planetary radius after 5 Gyr.Planets that lost their entire atmospheres are circled in red and labeled "rocky".Marker size depicts planet mass.Planets along the upper edge of the radius valley undergo the most fractionation.compared to that in CPML simulations (Lopez & Rice 2018; Gupta & Schlichting 2020).Malsky et al. (2023) report similar findings for EUV-driven photoevaporative mass loss, with planets in their simulations achieving He mass fractions ≥ 0.40 after 2.5 Gyr and ≥ 0.80 after 10 Gyr.The different parameter spaces occupied by helium-enhanced planets observed in this study offers a novel approach to investigate whether EUV-driven photoevaporation or CPML dominates sub-Neptune atmosphere evolution around low-mass stars.We discuss observational prospects to this end in Section 5.

Figure 4 .
Figure 4. Planet mass and final atmospheric mass fraction, fret, as a function of orbital period for the subset of simulated planets with helium-enriched atmospheres (xHe ≥ 0.145) for both loss mechanisms.The shaded areas represent the standard deviation.Mass is negatively correlated with orbital period for both mechanisms.Highly fractionated planets retain larger atmospheres and are generally more massive in the EUV-driven photoevaporation simulations.

Figure 5 .
Figure 5. Helium molar concentration (top panel) and deuterium mole fraction (bottom panel) as a function of orbital period and planet radius for EUV-driven photoevaporation (left) core-powered mass loss (middle) and the combined mechanisms (right).These plots are derived from simulations with n = 5 x 10 5 planets around M dwarfs, like that shown in Figure 3.The mean He (top) or D (bottom) abundance of all sub-Neptunes in each grid cell that have retained atmospheres after 5 Gyr is shown.The dashed green line represents the calculated radius valley as described in Section 4. Significant helium and deuterium enhancement is observed along the radius valley in all cases.

Figure 5
Figure 5 shows deuterium mole fraction as a function of orbital period and planet radius (Compare to Figure 2 in Gu & Chen (2023)).Note that the D/H mole fraction is capped at 100x Solar in Figure5and many planets retain atmospheres with essentially no hydrogen remaining, only deuterium.Though qualitatively similar to past work, our results demonstrate that deuterium enhancement was previously dramatically underestimated.We observe significant fractionation along the upper radius valley for all three mechanism combinations in our simulated data.Not surprisingly, we observe deuteriumenriched planets in the same parameter space as that of

Figure 6 .
Figure 6.Molar concentrations of deuterium relative to H (top) and He (bottom) as a function of retained atmospheric mass fraction, fret, for EUV-driven photoevaporation.The upper axis corresponds to the mean surface pressure for each fret bin.Simulated data correspond to the final state for 5 × 10 5 simulations of planets around an early M dwarf evolved for 5 Gyr.The red line shows an analytic solution for molar concentration of D as a function of fret and Teq (derived in Appendix B).The mean Teq = 500 K was used in this case.Venus and Earth dashed lines show deuterium concentrations corresponding to Venusian atmospheric D/H and the Vienna Standard Mean Ocean Water (VSMOW) values respectively.The "Gu & Chen" dashed line shows the maximum D/H value reported by Gu & Chen (2023) in their photoevaporation simulations.The bottom panel shows that most planets contain atmospheric deuterium in the hundreds of ppm.

Figure 7 .
Figure 7. Masses and radii for the subset of simulated planets with helium-enriched atmospheres (xHe ≥ 0.145; top panel) and deuterium-enriched planets (ND/NH ≥ 100x Solar; bottom panel) overlaid with planets around M dwarfs.Colors represent the mean incident planetary flux for all simulated planets in each grid cell.We include planet candidates whose calculated flux is within ± 1 dex of the interpolated flux to account for uncertainty in instellation history.
IsoFATE overcomes challenges presented by the MESA model used in previous studies by Malsky & Rogers (2020), Malsky et al. (2023), and Gu & Chen ( 7 µm spectral window open to detect HDO, a novel observing strategy potentially feasible for even modest D/H values with high-resolution ground-based spectroscopy (Mollière & Snellen 2019).Software: scikit-learn (Pedregosa et al. 2011), SciPy (Virtanen et al. 2020) B. ANALYTIC SOLUTION OF D/H IN A HELIUM-DOMINATED ATMOSPHERE

Table 1 .
Fractionated planet candidates predicted from Mp-Rp-Fp interpolation Planet Stellar type Mp [M⊕] Rp [R⊕] P [days] J-band mag TSM He/H Mechanisms (P frac ) D/H Mechanisms (P frac )