# Deuterium Escape on Photoevaporating Sub-Neptunes

Pin-Gao Gu  Institute of Astronomy & Astrophysics, Academia Sinica, Taipei
10617, Taiwan  Howard Chen  Department of Aerospace, Physics, and Space
Sciences, Florida Institute of Technology, Melbourne, FL 32901, USA  Planetary
Environments Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD
20771, USA

###### Abstract

We investigate the evolution of the deuterium-to-hydrogen (D/H) mass ratio
driven by EUV photoevaporation of hydrogen-rich atmospheres of close-in sub-
Neptunes around solar-type stars. For the first time, the diffusion-limited
approach in conjunction with energy-limited photoevaporation is considered in
evaluating deuterium escape from evolving exoplanet H/He envelopes. We find
that the planets with smaller initial gas envelopes and thus smaller sizes can
lead to weaker atmospheric escape, which facilitates hydrogen-deuterium
fractionation. Specifically, in our grid of simulations with low envelope mass
fraction less than 0.005, a low-mass sub-Neptune (4-$5M_{\oplus}$) at about
0.25-0.4 au or a high-mass sub-Neptune (10-$15M_{\oplus}$) at about 0.1-0.25
au can increase the D/H values by greater than 20% over 7.5 Gyr. Akin to the
helium-enhanced envelopes of sub-Neptunes due to photoevaporating escape, the
planets along the upper boundary of the radius valley are the best targets to
detect high D/H ratios. The ratio can rise by a factor of $\lesssim$ 1.65
within 7.5 Gyrs in our grid of evolutionary calculations. The D/H ratio is
expected to be higher in thinner envelopes as long as the planets do not
become bare rocky cores.

##  1 Introduction

Deuterium has long been one of the most studied isotopes in terms of its
abundance relative to hydrogen in a variety of astronomical environments,
despite being a trace element in the Universe. Deuterium was largely produced
in the Big Bang. The primordial deuterium-to-hydrogen (D/H) ratio inferred
from observations of the high-redshift intergalactic medium and cosmic
microwave background radiation constrained the baryon abundance in the Big
Bang nucleosynthesis and cosmological parameters (e.g., Tytler et al., 1996;
Planck Collabration, 2016). The D/H ratio in the present local interstellar
medium was inferred through Lyman-$\alpha$ line observations (e.g., Linsky et
al., 1995), which are lower than the protosolar value 4.55 Gyr ago implied
from solar wind measurements (e.g., Geiss & Gloeckler, 1998). This decrease in
D/H over time is expected because deuterium has been destroyed by
nucleosynthesis in stars (e.g., Lellouch et al., 2001). By contrast, the D/H
values can be considerably high in prestellar cores and protoplanetary disks
where the temperature is low enough to form deuterated molecules, which in
turn enrich the icy dust mantle through grain-surface chemistry (e.g., Cleeves
et al., 2014; Ceccarelli et al., 2014). Furthermore, within the Solar System,
the D/H ratios were measured or observed in chondrites, comets, planetary
atmospheres, and Earth’s ocean water to investigate how the different ratios
are possibly evolved through various fractionation processes (e.g., Morley et
al., 2019; Atreya et al., 2020; Piani et al., 2020; Nomura et al., 2023).

In the context of formation of giant planets in the Solar System, the D/H
ratio of Jupiter and Saturn is comparable to the protosolar value $\approx
2\times 10^{-5}$ (e.g., Geiss & Gloeckler, 1998; Pierel et al., 2017), roughly
consistent with formation through massive gas accretion from the solar nebula.
On the other hand, the D/H ratio of Uranus and Neptune is higher than the
protosolar value by a factor of $\sim 2.5$ (Feuchtgruber et al., 2013). The
theoretical interpretation is unclear (e.g., Guillot & Gautier, 2015). It
could be a natural consequence of icy giant planets, which accrete the nebular
gas and a significant amount of D‐enriched ice during their formation (e.g.,
Watson, 1974; Geiss & Reeves, 1981; Lecluse et al., 1996; Drouart et al.,
1999; Hersant et al., 2001).

For terrestrial planets, atmospheric escape has been posited to explain their
current D/H ratios. The D/H ratio in Venus’s atmosphere is much higher than
that of Earth’s ocean ($\approx 1.56\times 10^{-4}$) by about a factor of 100.
Previous work has shown that the rapid evaporation of the putative early ocean
via runaway greenhouse effect can possibly induce more loss of hydrogen to
space than the heavier isotope D, leading to the high D/H ratio on Venus
(Donahue & Pollack, 1983; Kulikov et al., 2006). The scenario of the runaway
greenhouse effect has been theorized to delineate the inner boundary of the
circumstellar habitable zone, where liquid water can exist on the crust of a
temperate planet in a range of orbital distances (Kasting et al., 1993;
Kopparapu et al., 2013; Kasting et al., 2015). In the cases of Earth and Mars,
their current D/H ratios are similar to those of carbonaceous chondrites and
Oort cloud comets, respectively. Assuming that these ratios have remained
constant for billions of years, the ratios can be attributed to the late
accretion of the carbonaceous asteroids and comets that were delivered from
the outer Solar System at or near the end of the stage of planet formation
(e.g., Matsui & Abe, 1986; Gomes et al., 2005; Drake, 2005; Raymond & Izidoro,
2017). However, atmospheric escape could further sculpt the D/H ratio after
the completion of planet formation through the accretion. Analogous to Venus,
early Mars may also have surface water as the major reservoir of hydrogen
(e.g., Ramirez et al., 2014; Batalha et al., 2015). As a result, the loss of
water and atmosphere from Mars’ surface could bring the D/H ratio to 5–7 times
larger than the terrestrial value (Krasnopolsky, 2015; Villanueva et al.,
2015). Based on formation modeling or geochemical evidence, the early Earth
could gravitationally accrete hydrogen gas from the remanent solar nebula to
form its primordial atmosphere (Ikoma & Genda, 2006; Marty, 2012; Lee &
Chiang, 2016). The subsequent loss of this H-rich atmosphere would raise the
protosolar D/H ratio to the current value, followed by deuterium exchange
between hydrogen gas and water vapor during the ocean formation (Genda &
Ikoma, 2008). The ingassing process of the H-rich atmosphere to Earth’s
interior was also posited to contribute to some of the water content in
Earth’s core during the impact phase of growing planetary embryos (Wu et al.,
2018)111An early steam atmosphere could also have been lost via impact erosion
by planetesimals during planetary accretion (e.g., Chen & Jacobson 2022),
which could lead to degassing of mantle hydrogen and removal from the
atmospheric reservoir..

Beyond the solar system, atmospheric escape plays a major role in the
demographic, evolutionary, and observational narratives of exoplanets. The
atmospheric loss from close-in exoplanets of Neptune to Jupiter masses has
been observed during planet transits (e.g., Vidal-Madjar et al., 2003;
Ehrenreich et al., 2015; Spack et al., 2018). The transit depths indicated by
the broad profiles of Lyman-$\alpha$ and He triple lines are greatly enhanced,
implying that the H-rich atmosphere can reach the Roche lobe with a velocity
beyond the escape velocity of these close-in planets (e.g., Owen, 2019, for a
recent review). Additionally, the Kepler transit survey (Borucki et al., 2010)
in conjunction with follow-up spectroscopic observations and parallax
measurements has revealed the so-called “radius valley”, i.e., the low
occurrence rate for super-Earths of planetary radii $\sim 1.8R_{\oplus}$
within the orbital period of $\sim 100$ days around FGK dwarfs (Fulton et al.,
2017; Fulton & Petigura, 2018). The radius valley was also identified by
asteroseismic measurements for a sample of bright planet-hosting stars, of
which the stellar radii were more accurately determined (Van Eylen et al.,
2018). The bimodal distribution of the planet radius separated by this radius
valley is consistent with the photoevaporation and core-powered mass-loss
models (e.g., Owen & Wu, 2013; Jin et al., 2014; Owen & Wu, 2017; Ginzburg et
al., 2018; Gupta & Schlichting, 2019). In these atmosphere loss models, a sub-
Neptune consists of a hydrogen-rich envelope and a rocky core. When the
envelope mass fraction is about a few percent of the total planet mass, the
mass-loss timescale becomes large enough for a sub-Neptune to survive (e.g.,
Chen & Rogers, 2016; Owen & Wu, 2017; Ginzburg et al., 2018).222 More
specifically, the mass-loss timescale needs to be larger than the radiative
cooling timescale for a contracting envelope in the core-powered mass-loss
model (Ginzburg et al., 2018; Gupta & Schlichting, 2019). Therefore, the rocky
planets larger than the radius valley can retain the H-rich atmosphere against
atmospheric escape over billions of years, whereas the rocky planets smaller
than the radius valley, i.e., super-Earths, almost lose their primordial
atmosphere and become almost bare.

Alternatively, the Kepler small planets larger than the radius valley could be
massive icy planets with a thin H atmosphere, which could have been partially
sculpted by atmospheric escape as well during planetary evolution (Zeng et
al., 2019; Venturini et al., 2020). Given the composition degeneracy,
Mordirrousta-Galian et al. (2021) subsequently modeled a wide range of
interior compositions and envelope mass fractions of extreme-ultraviolet (EUV)
photoevaporative planets compared to the observed bimodal size distribution.
The authors reached a similar conclusion to previous models for the
compositions of super-Earths and sub-Neptunes, with a more stringent
constraint on the initial mass distribution of these planets. Furthermore, the
atmospheric photoevaporation was applied to the post-evolution of the final
assembled planets through giant impacts in the attempt to simultaneously model
the similarity of planet size and orbital spacing in multiple planetary
systems (a.k.a. peas-in-a-pod, see a review by Weiss et al. (2023)), as well
as the radius valley for the Kepler small planets (Matsumoto et al., 2021).

As the loss of the H-rich atmosphere of terrestrial planets can produce
element fractionation in the Solar System, a similar effect is expected to
happen in exoplanets. To explain the lack of CH4 and dominance of CO in the
atmosphere of GJ 436b, Hu et al. (2015) proposed that the atmospheric loss
from warm Neptune- and sub-Neptune-sized exoplanets can strongly induce
diffusive separation of hydrogen and helium, leading to a preferential escape
of hydrogen and thereby yielding a helium-dominated atmosphere. Based on the
extension to the MESA333The Modules for Experiments in Stellar Astrophysics
(Paxton et al., 2011, 2013, 2015, 2018, 2019); also refer to
https://docs.mesastar.org. module that Chen & Rogers (2016) initially
developed, Malsky & Rogers (2020) revisited the problem by considering the
coevolution of the envelope photoevaporation and planet radius. The
compositional simulation coupled with thermal structure enabled the authors to
show that GJ 436b is too large to possess a He-rich atmosphere around a rocky
core. Nevertheless, the fractionation effect proposed by Hu et al. (2015) can
still be possible in the atmosphere of highly irradiated sub-Neptunes with low
envelope fractions. The predictive consequence is that the atmospheres of sub-
Neptunes along the edge of the radius valley may be helium-enhanced (Malsky et
al., 2023). If this is the case, however, then other light elements and
isotopes could also be subjected to dramatic modulations by atmospheric loss.

Motivated by the immense interest in the H-D fractionation for the solar-
system planets (i.e., formation, atmospheric evolution, and water delivery),
potential observations of deuterium fractionation for exoplanets have also
been posited. This would be achieved by detecting isotopic molecules in
exoplanet atmospheres with ongoing and upcoming facilities, such as the James
Webb Space Telescope (JWST) and the Extremely Large Telescope (Kofman &
Villanueva, 2019; Lincowski et al., 2019; Molliére & Snellen, 2019; Morley et
al., 2019). Despite these recent observational efforts and propositions, there
is as yet no theoretical study to model H-D fractionation processes for
exoplanets with significant H/He envelopes. In this Letter, we examine the
evolution of the D/H ratio for close-in sub-Neptunes due to atmospheric escape
based upon the aforementioned modeling for the H-He fractionation driven by
the photoevaporation of a primordial H-rich envelope.

##  2 Fractionation equations including deuterium escape

Using MESA version 12778, we simulate the coupled thermal and compositional
evolution of photoevaporating sub-Neptunes. The code we employed has been
substantially modified based on the module set up by Malsky et al. (2023);
i.e., we add the diffusive escape of deuterium in addition to the H-He
fractionation scheme. We consider a fiducial model for a solar-type star of
surface temperature 6000 K, a planet Bond albedo equal to 0.2, and the same
initial entropy of the planetary gas envelope as that specified by Malsky &
Rogers (2020). Malsky & Rogers (2020) found that the enhancement of helium
fraction due to the preferential hydrogen loss can shape the mass-radius
relation of the sub-Neptune-mass planet population. The effect is more
prominent for low-mass ($M_{p}\lesssim 10M_{\oplus}$) highly irradiated
planets with an initial envelope mass fraction $f_{env}$ below 1.0% (see
Figure 3 of Malsky & Rogers (2020)). Because deuterium’s mass lies between the
masses of H and He, it is expected that H and D fractionation due to
photoevaporation can also occur on close-in sub-Neptunes.

We begin the calculation of deuterium escape with the photoevaporative wind
described by Equation(A1) in Appendix A. After the inclusion of D in
Equation(A1, Equations (A2) and (A3) become444It can be argued that the energy
$f_{r}\Phi_{EL}$ is deposited differently among H, D, and He.

$\displaystyle\Phi\approx\Phi_{H}+\Phi_{D}+\Phi_{He}=4\pi
R_{h}^{2}(\phi_{H}m_{H}+\phi_{D}m_{D}+\phi_{He}m_{He}),$ (1)
$\displaystyle{\phi_{He}\over X_{He}}\approx{\phi_{H}\over
X_{H}}-\phi_{DL,He}+\left({1\over b_{D,He}}-{1\over
b_{D,H}}\right)b^{\prime}_{H,He}\phi_{D},$ (2)

where the escape rate $\Phi$ is given by Equation(B1) in Appendix B, $\phi$ is
the escaping number flux evaluated at the homopause radius $R_{h}$, $m$ is the
atomic mass of a species, $X$ is the mixing ratio of a species, $b$ is the
binary diffusion coefficient between two species, $b^{\prime}_{H,He}$ is
$b_{H,He}$ taking into account the H ionization, and $\phi_{DL,He}$ is the
diffusion-limited escape flux for He. Here $R_{h}$ is determined by the
location where the eddy diffusivity $K_{zz}$ equals the binary diffusivity
$\mathcal{D}_{H,He}$. Detailed descriptions of these quantities can be found
in Appendix B. In deriving Equation(2), we have assumed
$X_{D}\ll(b_{D,He}/b_{He,H})X_{H}$, $X_{D}\ll(b_{D,H}/b_{He,H})X_{H}$, and
$X_{H}+X_{He}\approx 1$. Because $1/b_{D,He}-1/b_{D,H}>1$ (see below), the
outflow of D may reduce the mass fractionation between H and He due solely to
the gravity (i.e., diffusion-limited) effect given by $\phi_{DL,He}$.
Furthermore, in Equation(2), the term due to $\phi_{D}$ is likely much smaller
than $\phi_{H}/X_{H}$ and could be ignored. In this case, Equation(2) can just
be approximated to Equation(A3).

Similarly, the number flux of D from Equation(A1) is given by

${\phi_{D}\over X_{D}}\approx{\phi_{H}-{GM_{p}(m_{D}-m_{H})\over
kTr_{0}^{2}}b_{D,H}+{b_{D,H}\over b_{He,H}}\left({X_{He}\over
X_{H}}-{\phi_{He}\over\phi_{H}}\right)\phi_{H}+{b_{D,H}\over
b_{D,He}}\phi_{He}\over X_{H}+{b_{D,H}\over b_{D,He}}X_{He}},$ (3)

where $G$ is the gravitational constant, $M_{p}$ is the planetary mass, $T$ is
the temperature, $k$ is the Boltzmann constant, $r_{0}$ is the homopause
radius for evaluating $\phi_{H}$, $\phi_{D}$, and $\phi_{He}$, and $X_{D}\ll
X_{H}$ has been applied to simplify the equation. If $\phi_{He}\neq 0$, the
above equation along with Equation(2) can be rewritten as

${\phi_{D}\over
X_{D}}\approx{\phi_{H}-\phi_{DL,D}+\alpha_{2}\phi_{DL,He}X_{He}+\alpha_{3}\phi_{He}\over
X_{H}+\alpha_{3}X_{He}},$ (4)

where

$\phi_{DL,D}\equiv{GM_{p}(m_{D}-m_{H})b^{\prime}_{D,H}\over
kTr_{0}^{2}},\qquad\alpha_{2}\equiv{b^{\prime}_{D,H}\over
b^{\prime}_{He,H}},\qquad\alpha_{3}\equiv{b^{\prime}_{D,H}\over b_{D,He}}.$
(5)

Note that in the above expressions, some of the $b$ coefficients have been
replaced with $b^{\prime}$ to take into account the H ionization (see below).
The four terms on the right-hand side of Equation(4) can be physically
interpreted as hydrogen drag (first term), helium drag (fourth term), and the
diffusion relative to hydrogen due to gravity (second term), which is
mitigated by the diffusion between hydrogen and helium due to gravity (third
term). It is evident from Equation(4) that
$\phi_{D}/X_{D}\rightarrow\phi_{H}/X_{H}$ as $b^{\prime}_{D,H}\rightarrow 0$
or as both $\phi_{DL,D}$ and $\phi_{DL,He}\rightarrow 0$ (i.e., no H-D
fractionation when D and H are strongly coupled). In addition, when $X_{He}=0$
and thus $X_{H}\approx 1$, Equation(4) becomes the same as Equation(A3) with
He replaced by D, as expected. For simplicity, we assume $r_{0}$ to be the
homosphere for hydrogen and helium $R_{h}$; hence, $T$ equals the homopause
temperature $T_{h}$. This assumption will be discussed near the end of the
paper.

Therefore, Equations(1), (2), and (3) (or eq(4) if $\Phi_{He}$ is nonzero) can
be solved for $\phi_{H}$, $\phi_{He}$, and $\phi_{D}$ evolving with time as
the envelope escapes at the rate given by $\Phi$. Since deuterium is a trace
element with $X_{D}\sim 10^{-5}$ to $10^{-4}$, it is expected that $\Phi_{H}$
and $\Phi_{He}$ would not change noticeably due to the presence of D. Hence,
it is an excellent approximation to use the same expressions for H and He mass
loss as those derived by Hu et al. (2015). The mass-loss rates of H, He, and D
then read

If $\Phi\leq\Phi_{crit,He}\equiv\phi_{DL,He}X_{H}m_{H}4\pi R_{h}^{2}$,

$\displaystyle\Phi_{H}$ $\displaystyle\approx$ $\displaystyle\Phi,$ (6)
$\displaystyle\Phi_{He}$ $\displaystyle\approx$ $\displaystyle 0,$ (7)
$\displaystyle\Phi_{D}$ $\displaystyle\approx$ $\displaystyle
X_{D}m_{D}{(1+\alpha_{2}X_{He}/X_{H})\Phi_{H}/m_{H}-\phi_{DL,D}4\pi
R_{h}^{2}\over X_{H}+\alpha_{3}X_{He}},$ (9) $\displaystyle{\rm however\ if\
}\Phi\leq\Phi_{crit,D}\equiv{1\over 1+\alpha_{2}X_{He}/X_{H}}\phi_{DL,D}4\pi
R_{h}^{2}m_{H},\ \Phi_{D}=0.$

If $\Phi>\Phi_{crit,He}$

$\displaystyle\Phi_{H}$ $\displaystyle\approx$ $\displaystyle{\Phi
m_{H}X_{H}+\phi_{DL,He}m_{H}m_{He}X_{H}X_{He}4\pi R_{h}^{2}\over
m_{H}X_{H}+m_{He}X_{He}},$ (10) $\displaystyle\Phi_{He}$
$\displaystyle\approx$ $\displaystyle{\Phi
m_{He}X_{He}-\phi_{DL,He}m_{H}m_{He}X_{H}X_{He}4\pi R_{h}^{2}\over
m_{H}X_{H}+m_{He}X_{He}},$ (11) $\displaystyle\Phi_{D}$ $\displaystyle\approx$
$\displaystyle X_{D}m_{D}{\Phi_{H}/m_{H}-\phi_{DL,D}4\pi
R_{h}^{2}+\alpha_{2}\phi_{DL,He}X_{He}4\pi
R_{h}^{2}+\alpha_{3}\Phi_{He}/m_{He}\over X_{H}+\alpha_{3}X_{He}}.$ (12)

Note that if we ignore the He species (i.e., $X_{He}=\Phi_{He}=0$),
Equations(9) and (12) are identical and have the same form as Equation(11)
with He replaced by D. In this case, without He, $\Phi_{D}=0$ when
$\Phi\leq\phi_{DL,D}4\pi R_{h}^{2}m_{H}$, analogous to the criterion for the
He escape with He replaced by D and $X_{H}\approx 1$.

We adopt the following values for diffusion coefficients: $b_{H,He}=1.04\times
10^{18}T_{h}^{0.732}$ cm-1 s-1 (Hu et al., 2015) and $b_{D,H}=7.183\times
10^{17}(T_{h})^{0.728}$ cm-1 s-1 (Genda & Ikoma, 2008). $b_{D,He}$ is crudely
estimated from $b_{D,H}$ using the square root of the reduced mass (Kasting &
Pollack, 1983; Genda & Ikoma, 2008), i.e.,
$b_{D,He}/b_{D,H}\approx\sqrt{\left({m_{D}+m_{He}\over
m_{D}m_{He}}\right)/\left({m_{D}+m_{H}\over m_{D}m_{H}}\right)}\approx$ 0.7,
which validates the aforementioned relation $1/b_{D,He}-1/b_{D,H}>1$ in
Equation(2). The correction from ionized species can be made following Hu et
al. (2015). Specifically, the binary diffusion coefficient $b_{D,H}$ corrected
for H ionization, denoted by $b^{\prime}_{D,H}$, is given by the relation
$kT/b^{\prime}_{D,H}=(1-x)kT/b_{D,H}+xm_{H^{+}}\nu_{D,H^{+}}/n_{D}$, where $x$
is the ionization fraction of H and the momentum transfer collisional
frequency $\nu_{i,H^{+}}=n_{i}\langle\sigma
v\rangle_{i,H^{+}}m_{i}/(m_{i}+m_{H^{+}})$ for the neutral species denoted by
$i$. We use $\nu_{He,H^{+}}/n_{He}=10.6\times 10^{-10}$ cm3/s (Schunk & Nagy,
1980; Hu et al., 2015) to estimate
$\nu_{D,H^{+}}\approx\sqrt{m_{H^{+}}m_{D}/(m_{H^{+}}+m_{D})\over
m_{H^{+}}m_{He}/(m_{H^{+}}+m_{He})}\nu_{He,H^{+}}(n_{D}/n_{He})$. In the MESA
module by Malsky et al. (2023), $x$ is determined by the Saha equation
ignoring photoionization. These values are expected to be insensitive to the H
ionization fraction $x$ (Hu et al., 2015). In this study, we ignore deuterium
ionization and its effect on the binary diffusion coefficients for simplicity,
as its momentum transfer collision frequencies with hydrogen and helium are
unknown.

Following a similar procedure for evolving $X_{H}$, $X_{He}$ and $X_{Z_{i}}$
by Malsky & Rogers (2020), we can evolve the abundance of D as follows:

$X_{D,n}={M_{env,n-1}X_{D,n-1}-(\Phi_{D}dt)\over
M_{env,n-1}-(\Phi_{He}+\Phi_{H}+\Phi_{D})dt},$ (13)

where $n$ is the index for the time step. We consider the same elemental
species as those in Malsky & Rogers (2020)555Eight elemental species are
considered in Malsky & Rogers (2020): 1H, 3He, 4He, 12C, 14N, 16O, 20Ne, and
24Mg. and add deuterium to the reaction network of MESA. The 3He escape is not
modeled because its effect on the deuterium escape is expected to be
negligible due to its extremely low abundance compared to H and He.

##  3 Results

We find that stellar EUV-induced hydrodynamic escape can elevate H-D
fractionation in the hydrogen-dominated envelopes of our modeled Kepler
planets; in many scenarios, the elevation in D/H is greater than 20% over 7.5
Gyr. Unless otherwise stated, we follow Malsky & Rogers (2020) and adopt the
fiducial values for the free parameters for the EUV photoevaporation
efficiency ($\eta=0.1$) and homopause conditions ($K_{zz}=10^{9}$ cm2/s and
$T_{h}=10000$ K) in this study. Changing the values of these free parameters
has a minor influence on our results for D/H (see Appendix C).

###  3.1 Typical evolution

Figure 1: Fractionation evolution of the envelope for 5 and 10 $M_{\oplus}$ at
$d=0.2$ au with various initial mass fractions of the envelope $f_{evn}$. Left
panels: mass-loss evolution of hydrogen ($X_{H}/X_{H,0}-1$, solid), deuterium
($X_{D}/X_{D,0}-1$, dashed), and helium ($X_{He}/X_{He,0}-1$, dotted). Middle
panels: Y evolution. Right panels: D/H evolution. The fractionation between H,
D, and He becomes more substantial at later times due to the decrease in the
mass-loss rate, leading to the increases in Y and D/H with time. The trend is
more prominent for the planet with a smaller initial mass fraction of the gas
envelope $f_{env}$ due to the even lower mass-loss rate. Notably, the mass
loss of He or D can almost stop at later times, as indicated by the flat
curves in the left panels.

The inclusion of deuterium escape into the compositional evolution of super-
Earths and sub-Neptunes allows us to simultaneously track the variation of
hydrogen, helium, and deuterium with time. Example evolution tracks of a 5 and
10 $M_{\oplus}$ planet with different $f_{env}$ at $d=0.2$ au experiencing
photoevaporation-driven fractionation are shown in Figure 1. The left panels
plot the deviation from the initial abundances, the middle panels present the
helium abundance Y (i.e., $X_{He}$), and the right panels display the D/H
ratio. As illustrated in the left panels of Figure 1, while the helium loss
fraction $X_{He}/X_{He,0}-1$ is clearly less than the hydrogen loss fraction
$X_{H}/X_{H,0}-1$ (Malsky & Rogers, 2020; Malsky et al., 2023), the deuterium
loss fraction $X_{D}/X_{D,0}-1$ is also fairly distinguishable; it is less
than the hydrogen but more than the helium loss fraction, as expected from
their mass differences. The abundance deviations of H, D, and He from their
initial values become more substantial at later times. This change is
expected. As the gas envelope contracts over time and the EUV flux from the
star also decreases substantially after 1 Gyr, the mass-loss rate $\Phi$
decreases with time according to Equation(B1).

Planets with lower $f_{env}$ experience the most pronounced change in their H,
D, and He abundance compared to their initial values. A planet with lower
$f_{env}$ has smaller cross-sectional radius and thus a lower mass-loss rate
during the evolution. Notably, the flat curves for $f_{env}=0.005$ and 0.01 at
the late stage in the bottom left panel of Figure 1 indicate that at the age
of $\sim$ a few Gyrs, the escape of deuterium and helium halts with the escape
of H in the atmosphere of the planet with 10 $M_{\oplus}$. This is because the
total mass-loss rate $\Phi$ starts to become smaller than both the critical
values for the helium loss rate $\Phi_{crit,He}$ and the deuterium loss rate
$\Phi_{crit,D}$ at this point, as described by Equations (6)-(9). In
comparison, $\Phi$ has not reached its critical value for deuterium in the
case of 5 $M_{\oplus}$ within 7.5 Gyr, primarily due to the higher mass-loss
rate from the weaker gravitational potential of a less massive planet.

As shown in the left panels of Figure 1, hydrogen can be lost by more (less)
than 20% for $f_{evn}\leq 0.01$ of a planet with 5 $M_{\oplus}$ (10
$M_{\oplus}$), and deuterium can be lost by more (less) than 10% for
$f_{evn}\leq 0.01$ of a planet with 5 $M_{\oplus}$ (10 $M_{\oplus}$).
Consequently, for a planet of 5 (10) $M_{\oplus}$ with $f_{evn}=0.005$-0.01,
the helium abundance Y increases to $\approx 0.26$-0.29 (0.25-0.26) from the
initial value 0.24 and the D/H ratio increases to $\approx 2.12$-2.3
(2.08-2.17) $\times 10^{-5}$ from the initial ratio 2$\times 10^{-5}$ within
7.5 Gyr, as illustrated in the middle and right panels of Figure 1. Although
all of our models start from $2\times 10^{-5}$ for the D/H ratio of the
protosolar value, it is worth mentioning that the deuterium abundance $X_{D}$
should almost scale as its initial value for a given planetary mass and
$f_{env}$. This is because deuterium is a trace element; thus, its variability
hardly changes the gas envelope mass during evolution. This can also be
realized from Equation(13) by noticing $\Phi_{D}\ll\Phi_{H}+\Phi_{He}$ and
$\Phi_{D}\propto X_{D,n-1}$. Therefore, for a planet of 5 (10) $M_{\oplus}$
with $f_{env}=0.005$-0.01, the D/H ratios increase by a factor of $\approx$
6-15% (4-8.5%) within 7.5 Gyr, regardless of the initial deuterium abundance.

Figure 2: Co-evolution of Y (top panels), D/H ratio (bottom panels), and
$R_{p}$ (color coding) for a grid of planetary models: $M_{p}=$4, 5, 10, 15
$M_{\oplus}$ and $f_{env}=$0.001, 0.003, 0.005, 0.1, 0.03, 0.05, 0.1 at each
orbital distance given by $d=$0.1, 0.15, 0.2, 0.25, 0.3, 0.35, and 0.4 au.
Some data with anomalous evolution of planet radius are removed (see text).
Similar to the evolution of helium-enhanced planets (top panels), the D/H
ratios increase with time due to atmospheric escape and have larger values for
the photoevaporating sub-Neptunes along the upper edge of the radius valley
(i.e., $R_{p}\approx 1.8-2R_{\oplus}$).

###  3.2 A grid of simulations

After understanding the basic outcomes of the H-D-He fractionation driven by
EUV photoevaporation, we perform a grid of simulations in a parameter space
covering the planet’s mass $=$(4, 5, 10, 15) $M_{\oplus}$, the orbital
distance $d=$(0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4) au, and $f_{env}=$(0.001,
0.003, 0.005, 0.1, 0.03, 0.05, 0.1). We find that in the presence of
fractionation, either a large escape rate (i.e., small $d$) or a thin envelope
(i.e., small $f_{env}$) of a low-mass planet somehow causes the atmosphere to
contract abruptly and thus unphysically due to the evolution to a wrong
equation of state in the MESA code. The data with the anomalous evolution of
planet radius are removed. Consequently, there are no planets with bare rocky
cores in our results. Nevertheless, the fate of these planets with anomalous
size evolution could be either the planets with a very thin H-He atmosphere or
bare cores. The parameter space and time frame for which the bare cores could
be produced in the simulations will be discussed in Section 4.1.

We find that planetary H/He envelopes are more helium-abundant toward
$R_{p}\approx 1.8$-$2R_{\oplus}$, particularly at later times. This can be
seen in Figure 2, in which we show the results of the helium mass fraction Y
and D/H ratio in the gas envelope of the model planets, along with the planet
radius $R_{p}$, at $t=2.5$, 5, and 7.5 Gyr. This figure resembles Figure 5 of
Malsky & Rogers (2020) and Figure 1 of Malsky et al. (2023) for the results of
the Y evolution from their grids of planetary models. However, because we are
not able to simulate a few cases with small $f_{evn}$, the Y values larger
than 0.4 in the photoevaporating envelope, as presented in Malsky et al.
(2023), are not produced during the evolution from our grid of models.
Nevertheless, Figure 2 shows that the photoevaporating planetary envelopes are
more helium-abundant toward $R_{p}\approx 1.8$-$2R_{\oplus}$ at later times,
consistent with the result of Malsky et al. (2023) that helium-enhanced
planets lie along the upper edge of the radius valley where $R_{p}\approx
1.8$-$2R_{\oplus}$.

Our calculated D/H ratios evolve in a similar manner to Y. This can be seen by
comparing the top and bottom panels of Figures 2, in which the gas envelopes
of the planets with larger Y generally exhibit higher D/H ratios toward the
upper edge of the radius valley at later times. It is also expected from the
typical evolutions shown in Figure 1. Specifically, the D/H ratio of the
planets evolving to $R_{p}\approx 1.8$-$2R_{\oplus}$ can rise to $\approx
2.5$, 3.0, and 3.3 $\times 10^{-5}$, or increase by a factor of 1.25, 1.5, and
1.65, at the ages of 2.5, 5, and 7.5 Gyrs, respectively. As explained in the
preceding paragraph, our grid of models limits the Y values because some cases
with small values of $f_{env}$ are removed. Hence, the D/H ratio of the
planets with a thinner atmosphere along the upper edge of the radius valley is
expected higher than those shown in the bottom panels of Figure 2, as long as
these planets do not evolve to bare rocky cores.

Figure 3: The mass-radius-distance (top panels) and the envelope-radius-
distance (bottom panels) parameter spaces for the deuterium-enhanced
atmosphere in our grid of simulations shown in Figure 2. In the parameter
space from $(R_{p},\ d)\approx(2.4R_{\oplus},\ 0.1$ au) to (1.7$R_{\oplus}$,
0.4 au), the data points increase with time, meaning more deuterium-enhanced
sub-Neptunes along the upper edge of the radius valley during the
photoevaporation evolution. As illustrated in the bottom panels, the planets
with initial low envelope fractions, i.e., $f_{env}=0.001$ and 0.003, can
become deuterium-enhanced, which is expected from Figure 1. The top panels
show that more massive planets, i.e., $M_{p}=10$ and 15$M_{\oplus}$, with low
$f_{env}$, can become deuterium-enhanced as the orbital distance decreases.

Our simulations result in more deuterium-enriched planets along the upper
boundary of the radius valley, particularly those with lower $R_{p}$. These
results are shown in Figure 2. By including $M_{p}$ and $f_{env}$, we identify
the parameter space for which greater deuterium enhancement occurs in the
planetary envelopes. Without the loss of generality, we define the D/H values
$\geq 2.4\times 10^{-5}$ to be deuterium-enriched, which corresponds to an
increase of 20% and more in the ratio. The results are shown in Figure 3 for
the $R_{p}$-$d$ parameter plane with color-coded values of $M_{p}$ (top
panels) and $f_{env}$ (bottom panels). As expected, the data points for
deuterium-enhanced planets lie along the lower boundary of the simulated
results presented in Figure 2. Furthermore, the bottom panels of Figure 3 show
that these planets initially possess thin atmospheres, i.e., $f_{env}=0.001$
and 0.003, as expected from the typical evolution illustrated in Figure 1. In
other words, the sub-Neptunes with a higher $f_{env}$ become less deuterium-
enhanced, which populate the upper part of the data points in Figure 2 (i.e.,
larger $R_{p}$). It is evident from Figure 3 that the deuterium-enhanced
planets populate more widely in the parameter space of $f_{env}$-$R_{p}$-$d$
from 2.5 to 7.5 Gyr. This trend is expected; namely, more deuterium-enriched
sub-Neptunes appear along the upper boundary of the radius valley during the
photoevaporation evolution. Moreover, the top panels of Figure 3 show that
more massive sub-Neptunes with low $f_{env}$ can become deuterium-enhanced as
the orbital distance decreases from 0.4 to 0.1 au. It arises because although
the incident EUV flux is stronger at smaller $d$, the deeper gravitational
potential of massive planets of low $f_{env}$ (and thus small $R_{p}$) can
reduce the mass-loss rate and hence facilitate the deuterium fractionation in
the photoevaporating atmosphere. In summary, in our grid of simulations with a
low envelope mass fraction of less than 0.005, a low-mass sub-Neptune
(4-$5M_{\oplus}$) at $\approx$0.25-0.4 au or a high-mass sub-Neptune
(10-$15M_{\oplus}$) at $\approx$0.1-0.25 au can increase the D/H values by 20%
or more over 7.5 Gyr. Note that the parameter spaces in each time frame shown
in Figure 3 are expected to be broader because some sub-Neptunes with small
$f_{env}$ at small $d$ are disregarded in the analysis due to an unphysical
radius anomaly in the simulations.

##  4 Discussion

###  4.1  Emergence of bare rocky cores and surviving sub-Neptunes

In our grid of simulations, the planets with rapid mass loss (i.e., lower $d$)
and thin initial atmospheres (i.e., lower $f_{env}$) can go through a rapid
contraction due to the incorrect equation of state in their evolution. As
described in Section 3.2, we expect that the planets in some of our false
simulations would lead to higher D/H ratios if they do not evolve to bare
cores. It is difficult to further identify which planets would become bare in
the simulations with incorrect evolution tracks. Nevertheless, we make an
attempt to explore the possible parameter space and time frame for the
emergence of bare cores in this subsection. We note that a subset of these
planets with anomalous radius evolution loses almost all the envelope and
stops evolving, while the rest of them can continue to evolve up to $t=7.5$
Gyr at the end of the simulations.

Despite the numerical limitation for the simulations with fractionations, no
issues associated with radius anomalies are found in the same grid of MESA
simulations in the absence of H-D-He fractionation. We find that the fates of
the planets, bare rocky cores or surviving sub-Neptunes, are almost identical
between the cases with and without fractionations, even though the simulations
involving fractionations with radius anomaly are incorrect. Given the great
similarities, we may use the results without fractionations to provide clues
about the parameter space and time frame for the appearance of bare cores in
the grid of simulations with fractionations.

Figure 4 presents the three-dimensional parameter spaces for simulated results
starting from the same initial conditions as those shown in Figure 2 but with
evolution in the absence of fractionations. The sub-Neptunes shown by gray
dots are almost identical to those presented in Figure 2 in terms of $R_{p}$
and $d$. We only show the results with $R_{p}\leq 2.5R_{\oplus}$ for clear
visualization of the parameter spaces and time frames for the emergence of
bare cores, which are illustrated by color-coded crosses. Because no sub-
Neptunes evolve to bare cores until $t\approx 1.8$ Gyr in the simulations,
Figure 4 presents the results in the time frames of 2.5, 5, and 7.5 Gyr, as
shown in Figure 2. It is clear that a radius valley appears between most bare
cores and sub-Neptunes (i.e., $1.6R_{\oplus}\lesssim R_{p}\lesssim
2R_{\oplus}$) and that the radius of the valley decreases with the orbital
period ($d=0.1$-0.4 au corresponds to an orbital period of about 10-90 days).
This location of the radius valley is generally consistent with the valley in
the photoevaporation model by Owen & Wu (2017, see their Figure 5). A more
detailed comparison is outside the scope of this study, as the radius valley
also depends on photoevaporation models and the core mass distribution, which
are not considered here (also see the next subsection for more
discussion).666For instance, Owen & Wu (2017) considered a Rayleigh
distribution for the core mass, which peaks around 3 $M_{\oplus}$ with a
standard deviation of 3 $M_{\oplus}$.

From Figure 4, it is clear that more sub-Neptunes evolve to bare cores over
gigayear timescales. Specifically, the $M_{p}$-$R_{p}$-$d$ and
$f_{env}$-$R_{p}$-$d$ parameter spaces of bare cores expand to larger $M_{p}$
(4-$10M_{\oplus}$), $R_{p}$ (1.5-$1.8R_{\oplus}$), $f_{env}$ (0.001-0.003),
and $d$ (0.1-0.25 au) from 2.5 to 7.5 Gyr. The bare cores emerge a bit late in
comparison with the timescale for the high EUV flux of young Sun-like stars,
which is about $10^{8}$ years (see Appendix B). It arises due to the
additional saturation factor $f_{r}$ in our photoevaporation model compared to
the energy-limited models.

On the other hand, the sub-Neptunes shown by green dots almost correspond to
those removed from Figure 2 due to radius anomaly but with complete evolution
over 7.5 Gyr. These planets populate the parameter spaces between those for
bare cores (crosses) and those for sub-Neptunes that correspond to the planets
shown in Figure 2 (gray dots). We suggest that the atmospheres of these
surviving sub-Neptunes would potentially yield higher D/H values than those
shown in the bottom panels of Figure 2.

Figure 4: Mass-radius-distance (top panels) and envelope-radius-distance
(bottom panels) parameter spaces and time frames for the appearance of bare
rocky cores (denoted by color-coded crosses) in our grid of simulations
without the H-D-He fractionation. The bare cores of the same $d$ and $R_{p}$
are plotted side by side to distinguish them. Despite no fractionations, the
sub-Neptunes shown by gray dots are identical to those in Figure 2. The sub-
Neptunes in the $d$-$R_{d}$ parameter space shown by green dots survive up to
$t=7.5$ Gyr in the grid of simulations without fractionations, whereas they
are disregarded from Figure 2 due to radius anomaly in the grid of simulations
with fractionations. The parameter spaces of bare cores become wider with
time, meaning more sub-Neptunes evolving to bare cores over time.

###  4.2 Photoevaporation models

In this study, we include deuterium in the coupled thermal/mass-
loss/compositional evolution of sub-Neptunes. However, our approach is similar
to Malsky & Rogers (2020) and Malsky et al. (2023) in that an energy-limited
approach for the photoevaporation is adopted, with the mass-loss rate
suppressed by a direct simulation Monte Carlo (DSMC) factor of $f_{r}$
(Johnson et al., 2013; Hu et al., 2015, also see Appendix B). The
photoevaporation efficiency $\eta$ in Equation(B1) depends on the radiative
cooling of the winds and is often taken as a constant of $\sim$1.5-4 (Johnson
et al., 2013). In the energy-limited approach ignoring kinetic effects, $\eta$
can be modeled as a power law of the escape velocity at the planet photosphere
(e.g., Owen & Wu, 2017; Rogers et al., 2021). When the EUV flux from young
host stars is high enough to allow for radiation recombination equilibrium in
planetary winds, the escaping flow is limited by radiative cooling via Lyman-
alpha lines (e.g., Murray et al., 2009; Chen & Rogers, 2016). Recent
hydrodynamical simulations suggested that $\eta$ can be expressed by an
analytical function of gravitational potential and/or incident EUV flux to
cover both energy- and radiation recombination-limited regimes of mass loss
(Salz et al., 2016; Caldiroli et al., 2022). In Appendix C, we show that our
model with $\eta=0.08$-0.75 does not significantly change the results because
the escape rate is limited by $f_{r}$ for the transonic flow driven by high
EUV flux (Hu et al., 2015). Therefore, we expect that applying the radiation
recombination-limited model to our grid of simulations would not dramatically
change the results of the H-D-He fraction. Future comparison studies based on
the aforementioned photoevaporation models are worthwhile to test our
expectation.

Kubyshkina et al. (2018a); Kubyshikina et al. (2018b) revised the
photoevaporation rate based on a grid of hydrodynamic models for hydrogen-only
atmospheres. The authors demonstrated that the mass-loss rate based on the
energy-limit model is underestimated by several orders of magnitude in the
regime of the low Jeans escape parameter $\Lambda$, which is a ratio of the
gravitation energy to the intrinsic thermal energy for a hydrogen atom in an
atmosphere. In this regime, the escape is primarily caused by low gravity and
high equilibrium temperature of planets and depends weakly on stellar EUV
flux. Hence, under the assumption that the revised model works properly for
standard atmospheres containing both hydrogen and helium, the mass-loss rate
of less massive planets should be drastically enhanced, resulting in different
H-D-He fractionations in the atmospheres. Nevertheless, the observed bimodal
size distribution and the radius valley for Kepler small planets can be
reproduced with the revised photoevaporation model as long as the planet
population shifts to a high-mass distribution (Mordirrousta-Galian et al.,
2021; Ketzer & Poppenhaeger, 2023). We then expect that the evolved D/H ratio
would be similar to what we derive and still be highest along the upper edge
of the radius valley if high-mass planets are considered. A study will be
conducted to examine the speculation using the latest version of MESA to
possibly cope with the issue of the equation of state in this fast-wind model.

Additionally, while our study focuses on solar-type stars with a fiducial
surface temperature of 6000 K and adopts the EUV flux described in Appendix B,
it is noteworthy that the time-integrated X-ray exposure of a planet at a
fixed incident bolometric flux decreases with stellar mass (McDonald et al.,
2019). To better understand the nature of the radius valley, this effect
depending on the stellar mass was considered to distinguish the
photoevaporation from the core-powered mass-loss model (Rogers et al., 2021;
Berger et al., 2023). The stronger mass loss driven by more intense X-rays may
be hostile to deuterium fractionation. Investigating the D/H ratio of planets
around various stellar types in relation to the radius valley in multiple-
parameter space is outside the scope of this work and will be studied in
future work.

###  4.3 Core-powered mass-loss model

In our grid of planet models with the same initial entropy specified by Malsky
& Rogers (2020), the initial temperature of the envelope base is a few
thousand kelvins. If we assume that the entire core temperature is the same as
that at the envelope base, our planet model would be consistent with the
molten rocky core scenario to drive a significant atmospheric escape, i.e.,
core-powered mass loss (Ginzburg et al., 2016, 2018). As described in the
Introduction, this alternative model also successfully produces the observed
radius valley (Gupta & Schlichting, 2019; Berger et al., 2023). When a low-
mass planet is young and thus large such that $\Lambda$ is small, the core-
powered mass loss of gas envelope can be much more substantial than the
classic energy-limited photoevaporation (Kubyshkina & Fossati, 2021). The
core-powered mass loss arises from the slow release of the core internal
energy through the optically thick envelope on $\sim$ Gyr timescales (e.g.,
Gupta & Schlichting, 2020; Rogers et al., 2021). Although heating due to
radioactive decay in the core has been implemented in our MESA module for
photoevaporative escape (Chen & Rogers, 2016), core-powered mass loss is not
included in our current model. It remains to be examined whether the escape
rate of the prolonged mass loss at later phases can decline close to the
critical diffusive escape rate of deuterium $\Phi_{crit,D}$ to enable
significant fractionation. In principle, both photoevaporation and core-
powered mass-loss mechanisms are expected to occur during the lifetime of a
sub-Neptune (Mordirrousta-Galian & Korenaga, 2023).

###  4.4 Deuterium homopause

In this work, $r_{0}$ in Equation(5) is simply taken to be $R_{h}$, which is
determined by the homopause of the hydrogen-helium mixture, i.e.,
$K_{zz}=\mathcal{D}_{H,He}$. In principle, the deuterium homopause does not
necessarily lie at the helium homopause in a hydrogen-rich atmosphere.
Specifically, since $X_{D}\ll X_{H}$ and $X_{He}$, we can consider D in the
binary mixture of H and He and calculate the mixing diffusivity of deuterium
in the mixture as follows (Tang et al., 2014):

${1\over\mathcal{D}_{D,H,He}}={X_{H}\over\mathcal{D}_{D,H}}+{X_{He}\over\mathcal{D}_{D,He}}.$
(14)

Hence, the location of the deuterium homopause is determined by the relation
$K_{zz}=\mathcal{D}_{D,H,He}$ rather than $K_{zz}=\mathcal{D}_{H,He}$. With
Equation(14), we can estimate the difference between $\mathcal{D}_{D,H,He}$
and $\mathcal{D}_{H,He}$ at a given altitude (i.e. given temperature and
pressure) to assess our simplified approach. To estimate $\mathcal{D}_{D,H}$
and $\mathcal{D}_{D,He}$ in Equation(14), the atomic diffusion volume of
deuterium $V_{D}$ is needed (see Equation(A5)), which is unknown. We expect
that the value of $V_{D}$ lies between those of $V_{H}$ and $V_{He}$. It then
follows that at a given altitude, $\mathcal{D}_{D,H,He}/\mathcal{D}_{H,He}$ is
about 1.1 and 0.98 for $V_{D}=V_{H}$ and $V_{He}$, respectively; namely, the
two binary diffusivities are almost identical. Furthermore, the H-He
fractionation hardly changes when the eddy diffusivity $K_{zz}$ alters from
$10^{7}$-$10^{11}$ cm2/s (Malsky & Rogers, 2020), and the same consequence
applies to the H-D fractionation in our model (see Appendix C). All of the
above reasons suggest that considering $R_{h}$ as the homopause for deuterium
is reasonable in our problem.

###  4.5 Observational prospects

Here we describe some potential observational avenues to test the suggested
D/H signatures. Assuming that the initial sub-Neptunes in the solar
neighborhood are close to the protosolar abundance, i.e., D/H$\approx 2\times
10^{-5}$ and CH3D/CH${}_{4}\approx 8\times 10^{-5}$, it would be possible for
future ground-based telescopes, such as the Extremely Large Telescope, using
the cross-correlation technique (Molliére & Snellen, 2019) to probe
vibrational emissions from deuterated methane CH3D at $\sim 4.7\mu$m toward
these planets with equilibrium temperatures below 600 K (i.e. roughly
corresponding to $d>0.2$ au in this study).

Previous work sampling a variety of Earth-like abundance isotopologues on
Venus-like atmospheres has shown that fractionation signals from species such
as HDO can be accomplished with JWST’s NIRSpec in as few as 10 transits
(Lincowski et al., 2019). Similarly, infrared spectroscopy on the atmospheres
of sub-Neptune-sized planets should aim to characterize HDO features at 3.7,
2.4, and 1.5 $\mu$m. Despite the lower degrees of fractionation in our
simulations (i.e., for a Venus-like atmosphere, D/H in water vapor is assumed
to be $\sim 100\times$ (D/H)ocean) compared to our case of $\gtrsim 2\times
1.65\times$ (D/H)protosolar), the much larger scale heights of extended sub-
Neptune planet envelopes might aid the detection of deuterium and other
species associated with photoevaporation-induced fractionation, particularly
for the strongly irradiated cases (e.g., $d=0.1$ au, Figure 2). An
investigation into the detectability of deuterated water through the
transmission spectra of nearby sub-Neptunes would be encouraged to constrain
our model. While significant deuterium fractionation could indicate a more
evolved atmosphere, even a non-detection of deuterium enhancement may be
attributed to planet age, other unknown assumptions (such as mass-loss
efficiency, initial disk D/H, and metallicity), or additional considerations
not included in this work (such as nonthermal escape, inclusion of interactive
chemistry, and the possibility of deep H2O reservoirs).

##  5 Summary

We further extend the MESA module based on the EUV photoevaporation model
adopted by Hu et al. (2015) and Malsky & Rogers (2020). For the first time, we
conduct compositional simulations coupled with the thermal interior structure
to study the evolution of D/H of photoevaporating sub-Neptunes. The key
results are summarized as follows.

  * •

We derive the critical diffusive escape rate for deuterium $\Phi_{crit,D}$ in
a hydrogen- and helium-rich atmosphere. The deuterium ceases to escape when
the atmospheric escape rate becomes smaller than $\Phi_{crit,D}$.

  * •

The planets with smaller $f_{env}$ and thus smaller $R_{p}$ can lead to
smaller $\dot{M}_{p}$, which facilitates fractionation. Specifically, in our
grid of simulations with low $f_{env}<0.005$, a low-mass sub-Neptune
(4-$5M_{\oplus}$) at $d\approx$ 0.25-0.4 au or a high-mass sub-Neptune
(10-$15M_{\oplus}$) at $d\approx$ 0.1-0.25 au can increase the D/H values by
20% or more over 7.5 Gyr.

  * •

Analogous to the helium-enhanced planets (Malsky et al., 2023), the planets
along the upper boundary of the radius valley are the best targets to detect
high D/H ratios in their thin atmospheres. The ratio can rise by a factor of
$\lesssim$ 1.65 within 7.5 Gyr in our grid of evolutionary calculations,
independent of the initial D/H value.

  * •

A few cases of small $f_{env}$ are removed in our grid of simulations due to
the numerical limitation. Therefore, the D/H ratio is expected to increase
more from the initial value. Assuming the initial D/H in the gas envelope of
sub-Neptunes is protosolar, Figure 5 illustrates the D/H value of the planets
along the radius valley from our study in comparison with those in Earth’s
ocean, as well as in the gas and icy planets in the Solar System.

Figure 5: The D/H ratio of giant and icy planets in the Solar System (J:
Jupiter, S: Saturn, U: Uranus, and N: Neptune), along with the simulated ratio
of the photoevaporating sub-Neptunes along the upper edge of the radius
valley. The D/H values of Earth’s ocean and the protosolar nebula are shown by
blue and gray horizontal lines, respectively. Analogous to the helium-enhanced
planets due to atmospheric escape (Malsky et al., 2023), the planets along the
upper boundary of the radius valley generally exhibit the largest D/H value in
their thin atmospheres (see Figure 2). The simulated D/H value along the upper
radius valley only provides the lower limit.

We thank Issac Malsky, Peter Bodenheimer, James Owen, and Teppei Okumura for
the informative discussions. We also thank the anonymous referee for helpful
comments that improved the quality of the manuscript. P.-G.G. acknowledges
support from the National Science and Technology Council in Taiwan through
grants NSTC 111-2112-M-001-037 and 112-2112-M-001-035.

##  Appendix A Fractionation due to H-rich atmospheric escape

We adopt the assumption that the lowermost subsonic regions of the wind
control differential escape (Zahnle et al., 1986, 1990). Therefore, in a great
depth of planetary winds where the flow is subsonic, its kinetic energy can be
ignored for our interest in species fractionation. Under this circumstance,
the equation of motion for an isothermal wind for the species $j$ in a steady
state, with H as the major constituent, may be approximated to (Zahnle et al.,
1990)

${1\over X_{j}}{dX_{j}\over dr}=-{GM_{p}(m_{j}-m_{H})\over
kTr^{2}}+\sum_{i}{r_{0}^{2}\over r^{2}b_{i,H}}\left({X_{i}\over
X_{H}}\phi_{H}-\phi_{i}\right)+\sum_{i}{r_{0}^{2}\over
r^{2}b{i,j}}\left(\phi_{i}-{X_{i}\over X_{j}}\phi_{j}\right),$ (A1)

where $X$ is the mixing ratio,777Note that $X$ in this paper is defined as the
mixing ratio to the entire gas mixture, while it is defined as the mixing
ratio relative to the H abundance in Zahnle et al. (1990). and $r$ is the
radial coordinate. In the above equation, the species $i$ exerts a drag force
on the flow of the species $j$, characterized by the binary diffusion
coefficient $b$ defined by $kT/(\mu_{ij}k_{ij})$, where $\mu_{ij}$ is the
reduced mass, and $k_{ij}$ is the collision rate. For a non-negligible flow
for the species $j$, its abundance would be much more uniform over one density
scale height of H, i.e., $d\ln X_{j}/dr\ll d\ln n_{H}/dr$. It follows from the
derivation of Equation(A1) that the term $d\ln X_{j}/dr$ in the above equation
may be ignored to allow for a simple analytical solution.

In the escaping flow consisting only of H and He from a photoevaporating
planet envelope, the number fluxes of H and He are governed by mass and
momentum conservation evaluated at the homopause $r_{0}=R_{h}$ (Hu et al.,
2015; Malsky & Rogers, 2020; Malsky et al., 2023):

$\displaystyle\Phi$ $\displaystyle\approx$
$\displaystyle\Phi_{H}+\Phi_{He}=4\pi
R_{h}^{2}(\phi_{H}m_{H}+\phi_{He}m_{He}),$ (A2) $\displaystyle{\phi_{He}\over
X_{He}}$ $\displaystyle=$ $\displaystyle{\phi_{H}\over X_{H}}-\phi_{DL,He},$
(A3)

Equation(A3) is derived from Equation(A1) by neglecting $d\ln X_{He}/dr$,
using $X_{H}+X_{He}\approx 1$, and employing the diffusion-limited escape flux
for He given by

$\phi_{DL,He}={GM_{p}(m_{He}-m_{H})b^{\prime}_{H,He}\over R_{h}^{2}kT_{h}},$
(A4)

where $T_{h}$ is the homopause temperature, $b^{\prime}_{H,He}$ is the binary
diffusion coefficient for H and He taking into account the H ionization.
Moreover, $R_{h}$ for H and He is determined by the location where the binary
diffusivity $\mathcal{D}_{H,He}$ equals the eddy diffusivity $K_{zz}$, for
which a nominal value $10^{9}$ cm2/s adopted by Malsky & Rogers (2020) is
used. $\mathcal{D}_{H,He}$ is determined by Fuller’s method (Fuller et al.,
1966; Tang et al., 2014; Malsky & Rogers, 2020):

$\mathcal{D}_{H,He}={10^{-3}T_{h}^{1.75}\left({m_{H}m_{He}\over
m_{H}+m_{He}}\right)^{-1/2}\over
P_{h}\left(V_{H}^{1/3}+V_{He}^{1/3}\right)^{2}},$ (A5)

where $P_{h}$ is the gas pressure in units of atm at $R_{h}$, and
$V_{H}(=1.98)$ and $V_{He}(=2.88)$ are the atomic diffusion volumes.

Once $P_{h}$ is identified, the homopause radius $R_{h}$ is obtained in a
hydrostatic atmosphere using the gas pressure, as well as constant values for
gravitational acceleration $g$, molecular weight $\mu$, and the pressure scale
height at MESA’s outermost zone where the optical depth is $2/3$ (i.e.,
photosphere). Moreover, we use the transit radius, defined by a pressure level
of 1 mbar, as the planetary radius $R_{p}$. Detailed descriptions of these
radii and atmospheric boundary conditions for an irradiated sub-Neptune can be
found in Malsky & Rogers (2020) and Malsky et al. (2023).

##  Appendix B Model for EUV photoevaporation

We adopt the mass-loss rate due to EUV photoevaporation given by Hu et al.
(2015), Malsky & Rogers (2020), and Malsky et al. (2023)

$\Phi=f_{r}\Phi_{EL},$ (B1)

where $\Phi_{EL}=L_{EUV}\eta a^{2}R_{h}^{3}/(4Kd^{2}GM_{p})$ is the energy-
limited mass-loss rate for a planet at an orbital distance $d$. Here $L_{EUV}$
is the host-star EUV luminosity evolving with the stellar age $\tau$ described
by the power law $\log_{10}(L_{EUV}/{\rm
J\,s^{-1}})=22.12-1.24\log_{10}(\tau/{\rm Gyr})$ for GK dwarfs at an age
$t\geq 10^{8}$ years (Sanz-Forcada et al., 2011). When $t<10^{8}$ years, the
stellar high-energy flux is almost saturated (Owen & Wu, 2017, references
therein). Hence, $L_{EUV}$ for $t<10^{8}$ years is set to be the same as the
power law of $L_{EUV}$ at a stellar age of $10^{8}$ years. Moreover in
Equation(B1), the heating efficiency $\eta=0.1$, the fraction of heating
radius $a=1$, and the Roche potential reduction factor $K$ are taken to be the
same values and expression as those in Hu et al. (2015) and Malsky & Rogers
(2020). The $f_{r}$ is introduced in Equation(B1) to recover the neglected
terms associated with the kinetic energy and $d\ln X_{j}/dr$ (thermal energy)
according to the DSMC, which corrects the overestimate of $\Phi$ due solely to
$\Phi_{EL}$ (Johnson et al., 2013; Hu et al., 2015).

##  Appendix C Sensitivity of D/H to free parameters

In our fiducial model, the free parameters $\eta=0.1$, $K_{zz}=10^{9}$ cm2/s,
and $T_{h}=10000$ K are adopted to study the D/H evolution of close-in sub-
Neptunes due to EUV photoevaporation. In this section, we investigate whether
the D/H values are sensitive to these free parameters. Besides the fiducial
values, we follow Malsky & Rogers (2020) and Malsky et al. (2023) and consider
the eddy diffusion coefficient $K_{zz}=10^{9}$ and $10^{11}$ cm2/s, and
homopause temperature $T_{h}=3000$ K. We also consider $\eta=0.08$ and $0.75$
for the low and high photoevaporation efficiency, respectively. The
investigation is carried out such that the value of one parameter is changed
at a time, while everything else remains the same as the fiducial values in
each simulation. Figure 6 shows the resulting D/H evolution (dashed and dotted
curves) compared with the typical evolution shown in Figure 1 (solid curves).

Overall, the D/H values are insensitive to these free parameters, except in
extreme cases for photoevaporation with $\eta=0.75$. For the influence of
$K_{zz}$ and $T_{h}$, it has been shown that the mass-loss rate increases with
$K_{zz}$ (Malsky & Rogers, 2020) and the level of H-He fractionation in the
escaping wind decreases with $T_{h}$ (Malsky et al., 2023). These effects
apply to H-D fractionation too, resulting in the outcome that the D/H values
slightly increase with $K_{zz}$ and decrease with $T_{h}$, as shown in Figure
6. However, the influences on the D/H values are insignificant.

On the other hand, the models with $\eta=0.75$ yield a much faster increase in
D/H than those in the fiducial model in late evolution, as shown in the lower
left panel of Figure 6 for the more massive sub-Neptune of 15$M_{\oplus}$ with
smaller initial envelope mass fractions, e.g., $f_{env}=0.005$. What happens
is that the mass-loss rates in most cases are approximately similar, even when
the mass loss is significant due to the suppression factor $f_{r}$ (Hu et al.,
2015). For a massive sub-Neptune with a small $f_{env}$, the mass-loss rate
starts to drop after $t\sim 1$ Gyr due to the weak EUV flux. Hence, $f_{r}$
becomes unity; thus, the magnitude of $\eta$ becomes influential during the
late stage. When $\eta=0.08$ and 0.1, deuterium can stop escaping during this
late phase due to $\Phi<\Phi_{D,crit}$, hence increasing the D/H value (see
the lower left panel of Figure 1). However, this increase is still less than
the rise of D/H ratio due to the faster diffusive escapes of both D and H when
$\eta=0.75$ even though $\Phi$ has not reached $\Phi_{D,crit}$. This effect
leads to the slightly large departure of the D/H values among different $\eta$
in the cases of $f_{env}=0.005$ and 0.1 at $t\gtrsim 2$-3 Gyr as shown in the
lower left panel of Figure 6. Despite this subtlety, we conduct the same grid
of simulations and find that the largest D/H value in the models with
$\eta=0.75$ and 0.08 is similar to that with $\eta=0.1$. This may infer that
our results are insensitive to $\eta=0.08$(low)-0.75(high photoevaporation
efficiency). With that being said, we caution that the parameter space of the
cases with the various values of $\eta$ is different due to numerical
limitations. For example, as shown in the upper left panel of Figure 6, the
case for the set of free parameters given by $M_{p}=5M_{\oplus}$,
$f_{env}=0.005$, and $\eta=0.75$ is disregarded due to the unphysical radius
anomaly.

Figure 6: Dependence of D/H evolution on the free parameters $\eta$, $K_{zz}$,
and $T_{h}$. The solid curves illustrate the typical evolution with the
fiducial parameters shown in Figure 1, while the dashed and dotted curves
present those with other values of free parameters. The curve for $\eta=0.75$
in the case of $f_{env}=0.005$ is not plotted in the upper left panel due to
the unphysical radius anomaly in the simulation. The evolution is insensitive
to the free parameters except in extreme cases for the high photoevaporation
efficiency $\eta=0.75$ shown in the lower left panel, where the D/H value of
the massive sub-Neptune with an initial thin atmosphere (e.g., $f_{env}=0.005$
and $M_{p}=10M_{\oplus}$) can be more enhanced at $t=7.5$ Gyr due to the
faster diffusive escapes of both H and D during the late stage, when
$t\gtrsim$1-2 Gyr.

## References

  * Atreya et al. (2020) Atreya, S. K., Hofstadter, M. H., Joong Hyun In, Mousis,, O., Reh K. & Wong, M. H. 2020, Space Science Reviews, 216, 18
  * Batalha et al. (2015) Batalha, N., Domagal-Goldman, S. D., Ramirez, R., & Kasting, J. F. 2015, Icarus, 258, 337
  * Berger et al. (2023) Berger, T. A., Schlieder, J. E., Huber, D., et al. 2023, arXiv:2302.00009
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
  * Ceccarelli et al. (2014) Ceccarelli, C., Caselli, P., Bockelée-Morvan, D., et al. 2014, Protostars and Planets VI, 859
  * Caldiroli et al. (2022) Caldiroli, A., Haardt, F., Gallo, E., et al. 2022, A&A, 663, A122
  * Chen & Rogers (2016) Chen, H., & Rogers, L. A. 2016, ApJ, 831, 180
  * Chen & Jacobson (2022) Chen, H., Jacobson, S. A. 2022, Earth and Planetary Science Letters, 594, p.117741
  * Cleeves et al. (2014) Cleeves, L. I., Bergin, E. A., Alexander, C. M. O., et al. 2014, Science, 345, 1590
  * Donahue & Pollack (1983) Donahue, T. M., & Pollack, J. B. 1983, in Venus, ed. D. M. Hunten et al. (Tucson, AZ: Univ. Arizona Press), 1003
  * Drake (2005) Drake, M. J. 2005, M&PS, 40, 519
  * Drouart et al. (1999) Drouart, A., Dubrule, B., Gauiter, D., & Robert, F. 1999, Icarus, 140, 129
  * Ehrenreich et al. (2015) Ehrenreich, D., Bourrier, V,, Wheatley, P. J,, Lecavelier des Etangs A., Hébrard, G, et al. 2015, Nature, 522, 459
  * Feuchtgruber et al. (2013) Feuchtgruber, H., Lellouch, E., Orton, G., et al. 2013, A&A, 551, 126
  * Fuller et al. (1966) Fuller, E. N., Schettler, P. D., & Giddings, J. C. 1966, Industrial & Engineering Chemistry, 58, 18
  * Fulton & Petigura (2018) Fulton, B. J., & Petigura, E. A. 2018, AJ, 156, 264
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A, Howard, A. W., Isaacson, H, Marcy, G. W., et al. 2017, AJ, 154, 109
  * Genda & Ikoma (2008) Genda, H., & Ikoma, M. 2008, Icarus, 194, 42
  * Geiss & Gloeckler (1998) Geiss, J. & Gloeckler, G. 1998, Space Sci. Rev., 84, 239
  * Geiss & Reeves (1981) Geiss, J., & Reeves, H. 1981, A&A, 93, 189
  * Ginzburg et al. (2016) Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29
  * Ginzburg et al. (2018) Ginzburg, S., Schlichting, H. E., & Sari, R. 2018, MNRAS, 476, 759
  * Gomes et al. (2005) Gomes, R., Levison, H.F., Tsiganis, K., Morbidelli, M., 2005. Nature 435, 466
  * Gupta & Schlichting (2019) Gupta, A., & Schlichting. H. E. 2019, MNRAS, 487, 24
  * Gupta & Schlichting (2020) Gupta, A., & Schiichting, H. E. 2020, MNRAS, 493, 792
  * Guillot & Gautier (2015) Guillot, T. & Gautier, D. 2015, Treatise on Geophysics, 529
  * Hersant et al. (2001) Hersant, F., Gautier, D., & Huré, J.-M. 2001, ApJ, 554, 391
  * Hu et al. (2015) Hu, R., Seager, S., Yung, Y. L. 2015, ApJ, 807, 8
  * Ikoma & Genda (2006) Ikoma, M., & Genda, H. 2006, ApJ, 648, 696
  * Jin et al. (2014) Jin, S., Mordasini, C., Parmentier, V., et al. 2014, ApJ, 795, 65
  * Johnson et al. (2013) Johnson, R. E., Volkov, A. N., & Erwin, J. T. 2013, ApJ, 768, L4
  * Kasting & Pollack (1983) Kasting, J. F., & Pollack, J. B., 1983, Icarus, 53, 479
  * Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108
  * Kasting et al. (2015) Kasting, J. F., Chen, H., & Kopparapu, R. K. 2015, ApJ, 813, L3
  * Ketzer & Poppenhaeger (2023) Ketzer, L., & Poppenhaeger, K. 2023, MNRAS, 518, 1683
  * Kofman & Villanueva (2019) Kofman, V., & Villanueva, G. L. 2019, EPSC-DPS Joint Meeting, 13, 685
  * Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R. & Kasting, J. F. 2013, ApJ, 765, 131
  * Krasnopolsky (2015) Krasnopolsky, V. A. 2015, Icarus, 257, 377
  * Kubyshkina et al. (2018a) Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018a, A&A, 619, A151
  * Kubyshikina et al. (2018b) Kubyshikina, D., Fossati, L., Erkaev, N. V., et al. 2018b, ApJ, 886, 18
  * Kubyshkina & Fossati (2021) Kubyshkina, D. I. & Fossati, L. 2021, Research Notes of the American Astronomical Society, 5, 74.
  * Kulikov et al. (2006) Kulikov, Y. N., Lammer, H., Lichtenegger, H. I. M., et al. 2006, Planet. Space Sci., 54, 1425
  * Lecluse et al. (1996) Lecluse, C., Robert, F., Gautier, D., & Guiraud, M. 1996, P&SS, 44, 1579
  * Lee & Chiang (2016) Lee, E. J., & Chiang, E. 2016, ApJ, 817, 90
  * Lellouch et al. (2001) Lellouch, E., Bézard, B., Fouchet, T., et al. 2001, A&A, 370, 610
  * Lincowski et al. (2019) Lincowski, A. P., Lustig-Yaeger, J., & Meadows, V. S. 2019, AJ, 158, 26
  * Linsky et al. (1995) Linsky, J. L., Diplas, A., Wood, B. E. et al. 1995, ApJ, 451, 335
  * Malsky & Rogers (2020) Malsky, I., & Rogers, L. A. 2020, ApJ, 896, 48
  * Malsky et al. (2023) Malsky, I., Rogers, L., Kempton, E. M.-R., et al. 2023, Nature Astronomy, 7, 57
  * Marty (2012) Marty, B. 2012, Earth and Planetary Science Letters, 313, 56
  * Matsui & Abe (1986) Matsui, T., & Abe, Y. 1986, Nature, 322, 526
  * Matsumoto et al. (2021) Matsumoto, Y., Kokubo, E., Gu, P.-G., et al. 2021, ApJ, 923, 81
  * McDonald et al. (2019) McDonald, G. D., Kreidberg, L., & Lopez, E. 2019, ApJ, 876, 22
  * Mordirrousta-Galian & Korenaga (2023) Mordirrousta-Galian, D., & Korenaga, J. 2023, ApJ, 934, 11
  * Mordirrousta-Galian et al. (2021) Mordirrousta-Galian, D., Locci, D., Micela, G. 2021, ApJ, 891, 158
  * Molliére & Snellen (2019) Molliére, P., & Snellen, I. A. G. 2019, A&A, 622, 139
  * Morley et al. (2019) Morley, C. V., Skemer, A. J., Miles, B. E. et al. 2019, ApJ, 882, 29
  * Murray et al. (2009) Murray-Clay R. A., Chiang E. I., Murray N., 2009, ApJ, 693, 23
  * Nomura et al. (2023) Nomura, H., Furuya, K., Cordiner, M. A., et al. 2023, Protostars and Planets VII, Univ. of Arizona Press, eds. Shu-ichiro Inutsuka, Yuri Aikawa, Takayuki Muto, Kengo Tomida, and Motohide Tamura
  * Owen & Wu (2013) Owen, J. E. & Wu, Y. 2013, ApJ, 775, 105
  * Owen & Wu (2017) Owen J. E., & Wu Y. 2017, ApJ, 847, 29
  * Owen (2019) Owen J. E., 2019, Annual Review of Earth and Planetary Sciences, 47, 67
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Piani et al. (2020) Piani, L., Marrocchi, Y., Rigaudier, T. et al. 2020, Science, 369, 1110
  * Pierel et al. (2017) Pierel, J. D. R., Nixon, C. A., Lellouch, E., et al. 2017, AJ, 154, 178
  * Planck Collabration (2016) Planck Collaboration 2016, A&A, 594, 13
  * Ramirez et al. (2014) Ramirez, R.M. et al., 2014, Nature Geoscience, 7, 59
  * Raymond & Izidoro (2017) Raymond, S. N., & Izidoro, A. 2017, Icarus, 297, 134
  * Rogers et al. (2021) Rogers, J. G., Gupta, A., Owen, J. E., & Schlichting, H. E. 2021, MNRAS, 508, 5886
  * Salz et al. (2016) Salz, M., Schneider, P. C., Czesla, S., & Schmitt, J. H. M. M. 2016, A&A, 585, L2
  * Sanz-Forcada et al. (2011) Sanz-Forcada, J., Micela, G., Ribas, I., et al. 2011, A&A, 532, A6
  * Schunk & Nagy (1980) Schunk, R. W., & Nagy, A. F. 1980, RvGeo, 18, 813
  * Spack et al. (2018) Spake, J. J., Sing, D. K., Evans, T. M., Oklopĉić, A., Bourrier, V, et al. 2018, Nature, 557, 68
  * Tytler et al. (1996) Tytler, D., Fan, X.-M., & Burles, S. 1996, Nature, 381, 207
  * Van Eylen et al. (2018) Van Eylen V., Agentoft C., Lundkvist M. S., Kjeldsen H., Owen J. E., Fulton B. J., Petigura E., Snellen I., 2018, MNRAS, 479, 4786
  * Venturini et al. (2020) Venturini, J. E., Guilera, O. M., Haldemann, J., Ronco, M. P., & Mordasini, C. 2020, Astronomy and astrophysics, 643, L1
  * Vidal-Madjar et al. (2003) Vidal-Madjar A, Lecavelier des Etangs A, Désert J. M, Ballester G. E, Ferlet R. et al. 2003, Nature, 422, 143
  * Villanueva et al. (2015) Villanueva, G. L., Mumma, M. J., Novak, R. E., et al. 2015, Science, 348, 218
  * Tang et al. (2014) Tang, M. J., Cox, R. A., & Kalberer, M. 2014, Atoms. Chem. Phys., 14, 9233
  * Watson (1974) Watson, W. D. 1974, ApJ, 188, 35
  * Weiss et al. (2023) Weiss, L. M., Millholland, S. C., Petigura, E. A., et al. 2023, Protostars and Planets VII, Univ. of Arizona Press, eds. Shu-ichiro Inutsuka, Yuri Aikawa, Takayuki Muto, Kengo Tomida, and Motohide Tamura
  * Wu et al. (2018) Wu, J., Desch, S. J., Schaefer, L., et al. 2018, Journal of Geophysical Research (Planets), 123, 2691
  * Zahnle et al. (1986) Zahnle, K., & Kasting, J. F. 1986, Icarus, 462, 480
  * Zahnle et al. (1990) Zahnle, K., Kasting, J. F., & Pollack, J. B. 1990, Icarus, 84, 502
  * Zeng et al. (2019) Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, Proceedings of the National Academy of Sciences, 116, 9723


# Novel Real-Time EMT-TS Modeling Architecture for Feeder Blackstart
Simulations

Victor Paduani, Bei Xu, David Lubkeman, Ning Lu North Carolina State
University, Raleigh, NC
<EMAIL_ADDRESS>This research is supported by the U.S. Department of Energy’s
Office of Energy Efficiency and Renewable Energy (EERE) under the Solar
EnergyTechnologies Office Award Number DE-EE0008770. Department of Electrical
and Computer Engineering

###### Abstract

This paper presents the development and benchmarking of a novel real-time
electromagnetic-transient and transient-stability (EMT-TS) modeling
architecture for distribution feeder restoration studies. The work presents
for the first time in the literature a real-time EMT-TS testbed in which the
grid-forming unit is simulated in EMT domain, operating as the slack bus of
the phasor domain while including unbalanced voltage conditions. To evaluate
the performance and limitations of the proposed model, an equivalent EMT
testbed is developed to serve as a benchmark. First, the co-simulation
framework and the domain coupling methodology are introduced. Next, the
steady-state operation of the EMT-TS and EMT models are matched. Then, tests
are conducted to analyze the transient performance of the proposed EMT-TS
model when compared to its EMT benchmark. Results reveal that when utilizing a
battery energy storage system (BESS) as the grid-forming unit, the EMT-TS
testbed can maintain high accuracy for typical load steps.

###### Index Terms:

Co-simulation, EMT-TS, hardware-in-the-loop, real-time simulation, transients.

##  I Introduction

As the participation of inverter-based resources (IBRs) in the grid increases,
new simulation methods that can integrate the fast response of power
electronic devices must be developed [1]. Conventional transient stability
(TS) software packages that run at the millisecond level are sufficient to
simulate synchronous generator dynamic behaviors in a larger network. However,
as synchronous generators are being displaced by energy systems powered by
IBRs with advanced grid-support functionalities [2], [3], electromagnetic
transient (EMT) simulations are required for capturing IBR dynamic responses
at the microsecond level.

EMT tools provide a greater depth of analysis for a wide frequency range,
whereas TS is less accurate if there are high frequency dynamics in the grid.
Nevertheless, the computational effort associated with EMT simulations make
them an unviable solution when simulating realistic distribution networks [4].
Therefore, by integrating EMT-modeled power electronic systems into the TS-
modeled power transmission and distribution networks, co-simulations can
become a reliable alternative for modeling grids with high penetration of
IBRs. As a consequence, the interest in EMT-TS co-simulations have been
growing in both industry and academia in recent years [5].

In addition, by developing co-simulation models in real-time simulators,
hardware-in-the-loop (HIL) systems can be used to test protection schemes in
large networks while including the dynamics of IBRs spread through the grid.
Because an HIL testbed can also represent realistic communication links
between distributed controllable devices and their centralized controllers,
energy management systems (EMS) coordinating the operation of IBRs can be
developed and tested in real-time or faster-than-real-time, day-long settings
[6].

In [7], Mongrain _et al._ introduced a real-time simulation test system with
established TCP/IP communication links on the eMEGASIM platform from OPAL-RT,
but no phasor domain is included. A EMT-TS model built also on eMEGASIM is
presented by Athaide _et al._ in [8], but it is not a real-time model.
Researchers from RTDS Technologies recently presented a model with a real-time
EMT-TS coupling via a Dynamic Phasor Line (DPL) [9]. However, the method
requires the DPL to have a propagation delay larger than the EMT domain
timestep with only a maximum of four coupling points allowed. Most
importantly, there is a need for quantitative studies that can highlight how
much information is lost by moving from an EMT model into an EMT-TS model,
especially for real-time operation. In [10], Song _et al._ compared the co-
simulation model with the network reduction model and analyzed the
computational costs of the two methods. Still, not enough is presented related
to the EMT and EMT-TS performance benchmarking.

To resolve the above issues, in this paper, we present a novel real-time EMT-
TS co-simulation modeling architecture, and conduct tests to validate the
performance of the proposed testbed against an equivalent model built in EMT
domain. The paper presents the following contributions to the literature.

  * •

Introduces for the first time in the literature a real-time EMT-TS testbed in
which the grid-forming unit is modeled in EMT domain, operating as the slack
bus of the phasor domain including unbalanced voltage conditions.

  * •

Analyzes the transient performance of the proposed model for different grid
operation conditions, and compares the results with an equivalent model
developed in EMT domain.

  * •

Presents a coupling method for multiple coupling points between the EMT and
phasor domains for moving devices across the feeder into the EMT domain.

##  II Methodology

Figure 1 displays the EMT-TS testbed, which consists of a microgrid connected
to a distribution feeder via a point of common coupling (PCC). The components
are split into two subsystems. The first subsystem includes distributed energy
resources, e.g., a grid-forming battery energy storage system (BESS), utility-
scale PV systems, diesel generators (DGs), and the grounding transformer. This
subsystem is simulated at the microsecond level in eMEGASIM. The second
subsystem includes distributed rooftop PVs, shunt capacitor banks, voltage
regulators, ZIP load models, and the unbalanced IEEE 123-bus network model
[11] used to emulate the opration of realistic distribution feeders. The
second subsystem is simulated in ePHASORSIM at the milisecond level. Note that
a subsystem can be assigned to more than one core, if needed, but each
subsystem requires at least one core for itself.

The EMS algorithms are running externally and can read measurements, set
devices status, and send power setpoint commands using Modbus registers. Note
that this is a function that a real-time simulation platform will normally
provide for modeling communication links between the real-time simulator and
an external computer.

Figure 1: Proposed EMT-TS co-simulation testbed framework.

###  II-A EMT-TS Co-Simulation Coupling

Since both subsystems are simulated in parallel and exchange information with
each other in real-time, there is a delay for events to propagate from one
susbystem to another. In this paper, we introduce a coupling method in which
the microgrid operates as the slack bus of the system, represented as a
voltage source in phasor domain. This allows the simulation of feeder
blackstart scenarios, in which a microgrid powered by BESS, utility-scale PVs,
and/or DGs can be used for picking up sections of the feeder during system
restoration.

####  II-A1 Time Interpolation Coupling

Time interpolation is a commonly used technique for coupling EMT and phasor
domains in co-simulations. The method consists of linearly interpolating both
magnitude and phases of the signals coming from phasor domain for each EMT
domain timestep [12]. Thus, in the proposed EMT-TS testbed, the distribution
feeder is represented as a current source, with current magnitude ($I_{x}$)
and angle ($\phi_{x}$) calculated from the phasor domain current ($I_{ph}$) as
follows.

$I_{x}(t+nh)=\lvert I_{\mathrm{ph}}(t)\rvert+\dfrac{n}{N}\Big{(}\lvert
I_{\mathrm{ph}}(t+H)\rvert-\lvert I_{\mathrm{ph}}(t)\rvert\Big{)}$ (1)

$\phi_{x}(t+nh)=\angle I_{\mathrm{ph}}(t)+\dfrac{n}{N}\Big{(}\angle
I_{\mathrm{ph}}(t+H)-\angle I_{\mathrm{ph}}(t)\Big{)}$ (2)

where the index $x$ represents phases ($a/b/c$), $h$ and $H$ correspond to the
EMT and phasor timesteps, respectively, $N=H/h$, and $n$ ranges from zero to
$N$. Thus, the current waveform for each phase is obtained by (3).

$i_{x}(t+nh)=\sqrt{2}I_{x}(t+nh)\cos{[\omega_{\mathrm{pcc}}(t+nh)+\phi_{x}(t+nh)]}$
(3)

where $\omega_{\mathrm{pcc}}$ is the angular frequency measured at the PCC
with a phase-locked loop.

As will be shown in the simulation results section, the time interpolation
coupling can smoothen the impact of fast current transients applied in the EMT
domain, causing a larger deviation for the response between the EMT and EMT-TS
testbeds. Therefore, instead of time interpolation, in this work, we directly
update the waveforms with signals coming from the phasor domain with (4). In
that case, to avoid numerical issues in the solver, a first-order low-pass
filter with time constant equal to the EMT timestep ($T_{\mathrm{s,emt}}$) is
applied on the recreated current signals. Moreover, due to the delay
introduced by the filter, a feedforward compensation ($T_{\mathrm{s,emt}}$) is
added to (4).

$i_{x}(t)=\sqrt{2}I_{x}\cos{\big{[}\omega_{\mathrm{pcc}}(t+T_{\mathrm{s,emt}})+\phi_{x}\big{]}}$
(4)

####  II-A2 Phasor Extraction

Phasor domain signals are obtained from the three-phase voltages measured at
the PCC at a rate given by the phasor domain timestep ($T_{\mathrm{s,ph}}$).
The voltages magnitudes are obtained with the true root mean square (RMS)
value of the corresponding waveform, whereas the angles are obtained with its
fundamental value, which can be obtained with Fourier analysis (7). Both the
RMS and phase calculations are performed over a running average window of one
cycle of the signal’s fundamental period ($T$).

$a=\frac{2}{T}\int_{(t-T)}^{t}v(t)\cos(\omega_{0}t)\,dt$ (5)

$b=\frac{2}{T}\int_{(t-T)}^{t}v(t)\sin(\omega_{0}t)\,dt$ (6)

$\phi_{\mathrm{ph}}=a\mathrm{tan}{(b/a)}$ (7)

###  II-B EMT Domain Benchmark Model

To benchmark the EMT-TS testbed performance, an equivalent EMT model with the
same network parameters and loading conditions is built using the combined
state-space nodal (SSN) modeling method introduced in [13]. SSN divides the
state-space matrix utilized by Simulink’s solver into separate matrices based
on nodes defined by the user. As shown in Fig. 2, the nodal voltage error
distribution between the EMT and EMT-TS testbeds is negligible in steady state
operation with the maximum error below 0.0006 p.u. and the interquartile range
within $\pm$0.0003 p.u.. In addition, the active and reactive power flow
measured at the PCC (EMT-domain side) between both testbeds has the highest
mismatch of 1.0 kW (on phase $b$) and 2.6 kVAR (on phase $a$) when supplying
peak load.

Figure 2: A comparison of the nodal voltage magnitude between the full EMT and
EMT-TS testbeds.

##  III Simulation Results

In this section, we present the following simulation results: (A) transient
performance evaluation during load steps, (B) the impact of the TS domain
timestep selection on propagation delays, (C) the impact of the time
interpolation coupling, (D) a capacitor bank switching event, and (E)
computational cost analysis.

###  III-A Transient Performance During Load Steps

The load steps applied during testing are presented in Table I. Two different
types of grid-forming units are analyzed: a 2-MVA BESS (parameters from [14])
and a 3.125-MVA diesel generator (parameters from [15]). The voltage and
frequency transient responses of the proposed EMT-TS and the full EMT model
are shown in Fig. 3. Note that a delay of 2 ms is added to the load steps in
the EMT model so that the disturbances are applied at the same point-on-wave
(POW) for both testbeds. The RMS error (RMSE) between EMT and EMT-TS curves is
included in the figure of each test, with the highest RMSE value corresponding
to the phase with the largest deviation.

Because the response of the DG is above the millisecond level, the modeling
error of the EMT-TS testbed is negligible. For the BESS unit, the transient
error is also negligible in the first load step-change (approximately 25% of
the BESS rated power). In the second load step-change (approximately 50% of
the BESS rated power), a maximum error of 0.006 p.u. is observed. During the
third load step change (close to the BESS rated power), a maximum RMS voltage
error of 0.03 p.u. is observed. Note that this is an extreme scenario and
thus, represents the upper bound of the modeling error. It is worth mentioning
that for the third load step-change in the diesel generator scenarios, (Fig.
3(c)), phase b (in red) presents a voltage rise at the beginning of the
transient. This is because the large, unbalanced load step causes a numerical
issue in the DG model. This can be fixed by changing the synchronous machine
solver from trapezoidal non-iterative to backward Euler robust (or trapezoidal
robust). However, trapezoidal non-iterative was the only solver compatible
with the SSN modeling (from the EMT benchmark). Therefore, the solver is kept
the same to ensure a consistent comparison between the two testbeds.

TABLE I: The Setup of Load Step-change Tests | P (kW) | Q (kVAr)
---|---|---
Step | S (kVA) | A | B | C | A | B | C
I | 448 | 196 | 39 | 166 | 97 | 20 | 81
II | 937 | 379 | 174 | 285 | 191 | 85 | 141
III | 1978 | 734 | 490 | 640 | 378 | 243 | 324
Figure 3: Transient response comparison between the EMT-TS and the EMT
testbeds. Plots on the left: 2-MVA BESS. Plots on the right: 3.125-MVA diesel
generator. (a) Load step change I; (b) Load step change II, and (c) Load step
change III. Refer to Table I for the test setup.

###  III-B Influence of Coupling Delay

In this test, the timestep of the phasor domain solver is increased from 1 to
10 ms. As discussed in [16], the propagation delay of events between the EMT
and TS domains can be as high as (2Ts,ph\+ T${}_{\mathrm{s,emt}})$. Hence, the
timestep change increases the maximum coupling delay between the EMT and co-
simulation testbeds from 2.1 to 20.1 ms. The delayed transient can be observed
in Fig. 4(a). Furthermore, since the transient is affected by the POW when the
load-step is applied, the voltage mismatch between the EMT and co-simulation
testbeds will increase accordingly.

Figure 4: Voltage and frequency transient responses during Step III from Table
I due to: (a) increasing the communication delay between EMT and phasor
domains; (b) utilizing the time interpolation coupling technique Figure 5:
Capacitor bank switching event: (a) peak-voltage energization (worst-case
scenario); (b) Zero-voltage-crossing energization.

###  III-C Influence of Interpolation

Figure 4(b) displays the transient response of the co-simulation testbed
during the third load step when including the interpolation method described
in Section II-C. Clearly, the interpolation considerably smoothed the
response, causing a significant reduction of both voltage and frequency
transients. This is an important finding that must be considered when
developing a co-simulation testbed in which the grid-forming unit (or slack
bus) is in EMT domain.

###  III-D Capacitor Switching Transients

A capacitor bank is turned on at bus 83 from the feeder. In EMT domain, the
bank is modeled as ideal Yg capacitances, whereas in the EMT-TS testbed it is
modeled inside ePHASORSIM as shunt loads. Figure 5 displays the voltages at
bus 83 when the capacitors are switches at (a) peak grid voltage (worst-case
scenario), and (b) zero-voltage crossing (ZVC).

Evidently, the co-simulation does not capture the peak voltage of 1.6 p.u.
observed in the EMT model in case (a). In fact, there is no significant
difference between (a) or (b) for the EMT-TS model, since the phasor domain
solver does not simulate the capacitors’ differential equations. Thus, as long
as the capacitor bank of the system in study contains protective measures such
as ZVC, the co-simulation can be utilized. Next, we present a method to
further improve the performance of the co-simulation for case (b).

By adding a second EMT-TS coupling (at bus 83), the capacitors can be modeled
in the EMT domain. First, the equivalent impedance between the feeder head and
bus 83 is found, and a series RLC circuit is built. Then, the domains are
coupled by (i) inserting a voltage source in EMT domain based on voltage
phasors measured at bus 83, and (ii) inserting a current source to bus 83 of
the phasor domain corresponding to the capacitors’ currents in EMT domain.
Figure 6 displays a comparison between the RMS currents measured for the EMT
and EMT-TS testbeds for case (b), whereas Fig. 7 compares the feeder bus
voltages during the switching event. Clearly, when simulating the capacitors
in EMT domain, their impact on the feeder bus voltages can be better captured.
This demonstrates the possibility of utilizing multi-EMT-TS couplings to move
devices distributed across the feeder into the EMT domain as needed.

Figure 6: Current injection of the capacitor bank: (a) EMT testbed; (b) EMT-TS
testbed with the shunt capacitors modeled in EMT domain. Figure 7: Feeder bus
voltages during capacitor bank switching: (a) EMT model; (b) EMT-TS model with
shunt capacitors modeled in phasor domain; (c) EMT-TS model with shunt
capacitors modeled in EMT domain.

###  III-E Computational Cost Analysis

The simulator used is OPAL-RT OP5700. The EMT testbed is simulated on one core
whereas the EMT-TS testbed is modeled on two cores (minimum required). For a
fair comparison, we set the EMT model to run four times slower than real-time
in one core, and the EMT-TS testbed to run two times slower than real-time in
two cores. The EMT model required an average of
$295.4\text{\,}\mathrm{\SIUnitSymbolMicro s}$ to process each
$100\text{\,}\mathrm{\SIUnitSymbolMicro s}$ step, consuming approximately 73%
of one core when running four times slower than real-time. The EMT-TS testbed
required an average of $21.2\text{\,}\mathrm{\SIUnitSymbolMicro s}$ per
$100\text{\,}\mathrm{\SIUnitSymbolMicro s}$ step, consuming around 11.4% of
two cores when running two times slower than real-time.

##  IV Conclusion

A novel real-time EMT-TS co-simulation testbed is presented, in which a grid-
forming unit modeled in EMT domain operates as a slack bus in the phasor
domain. The proposed testbed is validated against an equivalent EMT model
developed as a benchmark. Results demonstrate that the proposed model can
capture the transient response of BESS units for typical load steps.
Furthermore, we demonstrate that (i) the major impact of increasing the phasor
domain timestep is the POW time-drift due to the propagation delay; (ii) using
EMT-TS time interpolation coupling method will cause a smoothing effect when
modeling fast changing waveforms, and (iii) the framework can be used for
multi-EMT-TS coupling so that devices can be moved to the EMT domain such that
their transient responses can be captured in higher detail.

## References

  * [1] J. Matevosyan, J. MacDowell, N. Miller, B. Badrzadeh, D. Ramasubramanian, A. Isaacs, R. Quint, E. Quitmann, R. Pfeiffer, H. Urdal, T. Prevost, V. Vittal, D. Woodford, S. H. Huang, and J. O’Sullivan, “A future with inverter-based resources: Finding strength from traditional weakness,” _IEEE Power and Energy Magazine_ , vol. 19, no. 6, pp. 18–28, 2021.
  * [2] V. Paduani, L. Song, B. Xu, and N. Lu, “Maximum power reference tracking algorithm for power curtailment of photovoltaic systems,” in _2021 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 2021, pp. 1–5.
  * [3] V. Paduani and N. Lu, “Implementation of a two-stage PV system testbed with power reserves for grid-support applications,” in _2021 North American Power Symposium_ , 2021, pp. 1–6.
  * [4] H. Jain, B. A. Bhatti, T. Wu, B. Mather, and R. Broadwater, “Integrated transmission-and-distribution system modeling of power systems: State-of-the-art and future research directions,” _Energies_ , vol. 14, no. 1, p. 12, 2021.
  * [5] D. Shu, X. Xie, Q. Jiang, Q. Huang, and C. Zhang, “A novel interfacing technique for distributed hybrid simulations combining emt and transient stability models,” _IEEE Transactions on Power Delivery_ , vol. 33, no. 1, pp. 130–140, 2017.
  * [6] A. Shirsat, V. Muthukaruppan, R. Hu, N. Lu, M. Baran, D. Lubkeman, and W. Tang, “Hierarchical multi-timescale framework for operation of dynamic community microgrid,” _arXiv preprint arXiv:2011.10087_ , 2020.
  * [7] R. S. Mongrain, Z. Yu, and R. Ayyanar, “A real-time simulation testbed for hierarchical control of a renewable energy-based microgrid,” in _2019 IEEE Texas Power and Energy Conference (TPEC)_. IEEE, 2019, pp. 1–6.
  * [8] D. Athaide, J. Qin, and Y. Zou, “Matlab/simulink-based electromagnetic transient-transient stability hybrid simulation for electric power systems with converter interfaced generation,” in _2019 IEEE Texas Power and Energy Conference (TPEC)_. IEEE, 2019, pp. 1–6.
  * [9] H. Konara Mudiyanselage, “Interfacing a transient stability model to a real-time electromagnetic transient simulation using dynamic phasors,” 2020.
  * [10] X. Song, H. Cai, T. Jiang, T. Sennewald, J. Kircheis, S. Schlegel, L. N. Martinez, Y. Benzetta, and D. Westermann, “Research on performance of real-time simulation based on inverter-dominated power grid,” _IEEE Access_ , vol. 9, pp. 1137–1153, 2020.
  * [11] IEEE 123 Node Test Feeder. [Online]. Available: http://sites.ieee.org/pes-testfeeders/resources/
  * [12] F. Plumier, P. Aristidou, C. Geuzaine, and T. Van Cutsem, “Co-simulation of electromagnetic transients and phasor models: A relaxation approach,” _IEEE Transactions on Power Delivery_ , vol. 31, no. 5, pp. 2360–2369, 2016\.
  * [13] C. Dufour, J. Mahseredjian, and J. Bélanger, “A combined state-space nodal method for the simulation of power system transients,” _IEEE Transactions on Power Delivery_ , vol. 26, no. 2, pp. 928–935, 2010.
  * [14] B. Xu, V. Paduani, D. Lubkeman, and N. Lu, “A novel grid-forming voltage control strategy for supplying unbalanced microgrid loads using inverter-based resources,” _arXiv preprint arXiv:2111.09464_ , 2021.
  * [15] V. Paduani, H. Yu, B. Xu, and N. Lu, “A Unified Power-Setpoint Tracking Algorithm for Utility-Scale PV systems with Power Reserves and Fast Frequency Response Capabilities,” _IEEE Transactions on Sustainable Energy_ , pp. 1–12, 2021.
  * [16] F. Xie, C. McEntee, M. Zhang, and N. Lu, “An asynchronous real-time co-simulation platform for modeling interaction between microgrids and power distribution systems,” in _2019 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 2019, pp. 1–5.