Deuterium Escape on Photoevaporating Sub-Neptunes

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 Collaboration 2016). The D/H ratio in the present local interstellar medium was inferred through Lyα 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., Ceccarelli et al. 2014; Cleeves 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 of ≈2 × 10⁻⁵ (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 ∼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 (≈1.56 × 10⁻⁴) by about a factor of 100. Previous work has shown that the rapid evaporation of the putative early ocean via a 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, 2015; Kopparapu et al. 2013). 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; Drake 2005; Gomes et al. 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). ⁴

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; Spake et al. 2018). The transit depths indicated by the broad profiles of Lyα 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 ∼1.8 R_⊕ within the orbital period of ∼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). ⁵ 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 during planetary evolution (Zeng et al. 2019; Venturini et al. 2020). Given the composition degeneracy, Mordirrousta-Galian et al. (2020) 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 postevolution 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 (aka 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 CH₄ 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 MESA ⁶ 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.

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 the 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 ≲ 10 M_⊕) 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 also can 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) become ⁷

$$\begin{eqnarray}&&{\rm{\Phi }}\approx {{\rm{\Phi }}}_{{\rm{H}}}+{{\rm{\Phi }}}_{{\rm{D}}}+{{\rm{\Phi }}}_{\mathrm{He}}=4\pi {R}_{h}^{2}({\phi }_{{\rm{H}}}{m}_{{\rm{H}}}+{\phi }_{{\rm{D}}}{m}_{{\rm{D}}}+{\phi }_{\mathrm{He}}{m}_{\mathrm{He}}),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\phi }_{\mathrm{He}}}{{X}_{\mathrm{He}}}\approx \displaystyle \frac{{\phi }_{{\rm{H}}}}{{X}_{{\rm{H}}}}-{\phi }_{\mathrm{DL},\mathrm{He}}+\left(\displaystyle \frac{1}{{b}_{{\rm{D}},\mathrm{He}}}-\displaystyle \frac{1}{{b}_{{\rm{D}},{\rm{H}}}}\right){b}_{{\rm{H}},\mathrm{He}}^{{\prime} }{\phi }_{{\rm{D}}},\end{eqnarray} \tag{ 2 }$$

where the escape rate Φ is given by Equation (B1) in Appendix B, ϕ 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}_{{\rm{H}},\mathrm{He}}^{{\prime} }$ is b_H,He taking into account the H ionization, and ϕ_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 }}_{{\rm{H}},\mathrm{He}}$. Detailed descriptions of these quantities can be found in Appendix B. In deriving Equation (2), we have assumed X_D ≪ (b_D,He/b_He,H)X_H, X_D ≪ (b_D,H/b_He,H)X_H, and X_H + X_He ≈ 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 ϕ_DL,He. Furthermore, in Equation (2), the term due to ϕ_D is likely much smaller than ϕ_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

$$\begin{eqnarray}&&\displaystyle \frac{{\phi }_{{\rm{D}}}}{{X}_{{\rm{D}}}}\approx \displaystyle \frac{{\phi }_{{\rm{H}}}-\tfrac{{{GM}}_{p}({m}_{{\rm{D}}}-{m}_{{\rm{H}}})}{{{kTr}}_{0}^{2}}{b}_{{\rm{D}},{\rm{H}}}+\tfrac{{b}_{{\rm{D}},{\rm{H}}}}{{b}_{\mathrm{He},{\rm{H}}}}\left(\tfrac{{X}_{\mathrm{He}}}{{X}_{{\rm{H}}}}-\tfrac{{\phi }_{\mathrm{He}}}{{\phi }_{{\rm{H}}}}\right){\phi }_{{\rm{H}}}+\tfrac{{b}_{{\rm{D}},{\rm{H}}}}{{b}_{{\rm{D}},\mathrm{He}}}{\phi }_{\mathrm{He}}}{{X}_{{\rm{H}}}+\tfrac{{b}_{{\rm{D}},{\rm{H}}}}{{b}_{{\rm{D}},\mathrm{He}}}{X}_{\mathrm{He}}},\end{eqnarray} \tag{ 3 }$$

where G is the gravitational constant; M_p is the planetary mass; T is the temperature; k is the Boltzmann constant; r₀ is the homopause radius for evaluating ϕ_H, ϕ_D, and ϕ_He; and X_D ≪ X_H has been applied to simplify the equation. If ϕ_He ≠ 0, the above equation along with Equation (2) can be rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{{\phi }_{{\rm{D}}}}{{X}_{{\rm{D}}}}\approx \displaystyle \frac{{\phi }_{{\rm{H}}}-{\phi }_{\mathrm{DL},{\rm{D}}}+{\alpha }_{2}{\phi }_{\mathrm{DL},\mathrm{He}}{X}_{\mathrm{He}}+{\alpha }_{3}{\phi }_{\mathrm{He}}}{{X}_{{\rm{H}}}+{\alpha }_{3}{X}_{\mathrm{He}}},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{\phi }_{\mathrm{DL},{\rm{D}}}\equiv \displaystyle \frac{{{GM}}_{p}({m}_{{\rm{D}}}-{m}_{{\rm{H}}}){b}_{{\rm{D}},{\rm{H}}}^{{\prime} }}{{{kTr}}_{0}^{2}},\ {\alpha }_{2}\equiv \displaystyle \frac{{b}_{{\rm{D}},{\rm{H}}}^{{\prime} }}{{b}_{\mathrm{He},{\rm{H}}}^{{\prime} }},\ {\alpha }_{3}\equiv \displaystyle \frac{{b}_{{\rm{D}},{\rm{H}}}^{{\prime} }}{{b}_{{\rm{D}},\mathrm{He}}}.\end{eqnarray} \tag{ 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 ϕ_D/X_D → ϕ_H/X_H as ${b}_{{\rm{D}},{\rm{H}}}^{{\prime} }\to 0$ or as both ϕ_DL,D and ϕ_DL,He → 0 (i.e., no H–D fractionation when D and H are strongly coupled). In addition, when X_He = 0 and thus X_H ≈ 1, Equation (4) becomes the same as Equation (A3) with He replaced by D, as expected. For simplicity, we assume r₀ 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)–(3) (or Equation (4) if Φ_He is nonzero) can be solved for ϕ_H, ϕ_He, and ϕ_D evolving with time as the envelope escapes at the rate given by Φ. Since deuterium is a trace element with X_D ∼ 10⁻⁵–10⁻⁴, it is expected that Φ_H and Φ_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 ${\rm{\Phi }}\leqslant {{\rm{\Phi }}}_{\mathrm{crit},\mathrm{He}}\equiv {\phi }_{\mathrm{DL},\mathrm{He}}{X}_{{\rm{H}}}{m}_{{\rm{H}}}4\pi {R}_{h}^{2}$,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{H}}}\approx {\rm{\Phi }},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{He}}\approx 0,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{D}}}\approx {X}_{{\rm{D}}}{m}_{{\rm{D}}}\displaystyle \frac{(1+{\alpha }_{2}{X}_{\mathrm{He}}/{X}_{{\rm{H}}}){{\rm{\Phi }}}_{{\rm{H}}}/{m}_{{\rm{H}}}-{\phi }_{\mathrm{DL},{\rm{D}}}4\pi {R}_{h}^{2}}{{X}_{{\rm{H}}}+{\alpha }_{3}{X}_{\mathrm{He}}};\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\rm{h}}{\rm{o}}{\rm{w}}{\rm{e}}{\rm{v}}{\rm{e}}{\rm{r}},\,{\rm{i}}{\rm{f}}\,{\rm{\Phi }}\leqslant {{\rm{\Phi }}}_{{\rm{c}}{\rm{r}}{\rm{i}}{\rm{t}},{\rm{D}}}\equiv \displaystyle \frac{1}{1+{\alpha }_{2}{X}_{{\rm{H}}{\rm{e}}}/{X}_{{\rm{H}}}}\\ \,\times \,{\phi }_{{\rm{D}}{\rm{L}},{\rm{D}}}4\pi {R}_{h}^{2}{m}_{{\rm{H}}},{\,{\rm{\Phi }}}_{{\rm{D}}}=0.\end{array}\end{eqnarray} \tag{ 9 }$$

If Φ > Φ_crit,He,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{H}}}\approx \displaystyle \frac{{\rm{\Phi }}{m}_{{\rm{H}}}{X}_{{\rm{H}}}+{\phi }_{\mathrm{DL},\mathrm{He}}{m}_{{\rm{H}}}{m}_{\mathrm{He}}{X}_{{\rm{H}}}{X}_{\mathrm{He}}4\pi {R}_{h}^{2}}{{m}_{{\rm{H}}}{X}_{{\rm{H}}}+{m}_{\mathrm{He}}{X}_{\mathrm{He}}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{He}}\approx \displaystyle \frac{{\rm{\Phi }}{m}_{\mathrm{He}}{X}_{\mathrm{He}}-{\phi }_{\mathrm{DL},\mathrm{He}}{m}_{{\rm{H}}}{m}_{\mathrm{He}}{X}_{{\rm{H}}}{X}_{\mathrm{He}}4\pi {R}_{h}^{2}}{{m}_{{\rm{H}}}{X}_{{\rm{H}}}+{m}_{\mathrm{He}}{X}_{\mathrm{He}}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{D}}}\approx {X}_{{\rm{D}}}{m}_{{\rm{D}}}\displaystyle \frac{{{\rm{\Phi }}}_{{\rm{H}}}/{m}_{{\rm{H}}}-{\phi }_{\mathrm{DL},{\rm{D}}}4\pi {R}_{h}^{2}+{\alpha }_{2}{\phi }_{\mathrm{DL},\mathrm{He}}{X}_{\mathrm{He}}4\pi {R}_{h}^{2}+{\alpha }_{3}{{\rm{\Phi }}}_{\mathrm{He}}/{m}_{\mathrm{He}}}{{X}_{{\rm{H}}}+{\alpha }_{3}{X}_{\mathrm{He}}}.\end{eqnarray} \tag{ 12 }$$

Note that if we ignore the He species (i.e., X_He = Φ_He = 0), Equations (8) and (12) are identical and have the same form as Equation (11) with He replaced by D. In this case, without He, Φ_D = 0 when ${\rm{\Phi }}\leqslant {\phi }_{\mathrm{DL},{\rm{D}}}4\pi {R}_{h}^{2}{m}_{{\rm{H}}}$, analogous to the criterion for the He escape with He replaced by D and X_H ≈ 1.

We adopt the following values for the diffusion coefficients: ${b}_{{\rm{H}},\mathrm{He}}=1.04\times {10}^{18}{T}_{h}^{0.732}$ cm⁻¹ s⁻¹ (Hu et al. 2015) and ${b}_{{\rm{D}},{\rm{H}}}=7.183\times {10}^{17}{({T}_{h})}^{0.728}$ cm⁻¹ s⁻¹ (Genda & Ikoma 2008). Here 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 ≈ $\sqrt{\left(\tfrac{{m}_{{\rm{D}}}+{m}_{\mathrm{He}}}{{m}_{{\rm{D}}}{m}_{\mathrm{He}}}\right)/\left(\tfrac{{m}_{{\rm{D}}}+{m}_{{\rm{H}}}}{{m}_{{\rm{D}}}{m}_{{\rm{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}_{{\rm{D}},{\rm{H}}}^{{\prime} }$, is given by the relation ${kT}/{b}_{{\rm{D}},{\rm{H}}}^{{\prime} }=(1-x){kT}/{b}_{{\rm{D}},{\rm{H}}}+{{xm}}_{{{\rm{H}}}^{+}}{\nu }_{{\rm{D}},{{\rm{H}}}^{+}}/{n}_{{\rm{D}}}$, where x is the ionization fraction of H, and the momentum transfer collisional frequency ${\nu }_{i,{{\rm{H}}}^{+}}={n}_{i}\langle \sigma v{\rangle }_{i,{{\rm{H}}}^{+}}{m}_{i}/({m}_{i}+{m}_{{{\rm{H}}}^{+}})$ for the neutral species denoted by i. We use ${\nu }_{\mathrm{He},{{\rm{H}}}^{+}}/{n}_{\mathrm{He}}=10.6\,\times {10}^{-10}$ cm³ s^–1 (Schunk & Nagy 1980; Hu et al. 2015) to estimate ${\nu }_{{\rm{D}},{{\rm{H}}}^{+}}\approx \sqrt{\tfrac{{m}_{{{\rm{H}}}^{+}}{m}_{{\rm{D}}}/({m}_{{{\rm{H}}}^{+}}+{m}_{{\rm{D}}})}{{m}_{{{\rm{H}}}^{+}}{m}_{\mathrm{He}}/({m}_{{{\rm{H}}}^{+}}+{m}_{\mathrm{He}})}}{\nu }_{\mathrm{He},{{\rm{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:

$$\begin{eqnarray}&&{X}_{{\rm{D}},n}=\displaystyle \frac{{M}_{\mathrm{env},n-1}{X}_{{\rm{D}},n-1}-({{\rm{\Phi }}}_{{\rm{D}}}{dt})}{{M}_{\mathrm{env},n-1}-({{\rm{\Phi }}}_{\mathrm{He}}+{{\rm{\Phi }}}_{{\rm{H}}}+{{\rm{\Phi }}}_{{\rm{D}}}){dt}},\end{eqnarray} \tag{ 13 }$$

where n is the index for the time step. We consider the same elemental species as those in Malsky & Rogers (2020) ⁸ and add deuterium to the reaction network of MESA. The ³He 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.

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 (η = 0.1) and homopause conditions (K_zz = 10⁹ cm² s⁻¹ and T_h = 10,000 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

The inclusion of deuterium escape into the compositional evolution of super-Earths and sub-Neptunes allow us to simultaneously track the variation of hydrogen, helium, and deuterium with time. Example evolution tracks of a 5 and 10 M_⊕ 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 Φ decreases with time according to Equation (B1).

Zoom In Zoom Out Reset image size

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 a 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 an age of ∼a few gigayears, the escape of deuterium and helium halts with the escape of H in the atmosphere of the planet with 10 M_⊕. This is because the total mass-loss rate Φ starts to become smaller than both the critical values for the helium loss rate Φ_crit,He and the deuterium loss rate Φ_crit,D at this point, as described by Equations (6)–(9). In comparison, Φ has not reached its critical value for deuterium in the case of 5 M_⊕ 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 ≤ 0.01 of a planet with 5 (10) M_⊕, and deuterium can be lost by more (less) than 10% for f_evn ≤ 0.01 of a planet with 5 (10) M_⊕. Consequently, for a planet of 5 (10) M_⊕ with f_evn = 0.005–0.01, the helium abundance Y increases to ≈0.26–0.29 (0.25–0.26) from the initial value of 0.24, and the D/H ratio increases to ≈2.12–2.3 (2.08–2.17) × 10⁻⁵ from the initial ratio of 2 × 10⁻⁵ within 7.5 Gyr, as illustrated in the middle and right panels of Figure 1. Although all of our models start from 2 × 10⁻⁵ 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 Φ_D ≪ Φ_H + Φ_He and Φ_D ∝ X_D,n−1. Therefore, for a planet of 5 (10) M_⊕ with f_env = 0.005–0.01, the D/H ratios increase by a factor of ≈6%–15% (4%–8.5%) within 7.5 Gyr, regardless of the initial deuterium abundance.

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_⊕, 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 ≈ 1.8–2 R_⊕, 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 ≈ 1.8–2 R_⊕ 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 ≈ 1.8–2 R_⊕.

Zoom In Zoom Out Reset image size

Our calculated D/H ratios evolve in a similar manner to Y. This can be seen by comparing the top and bottom panels of Figure 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 evolution shown in Figure 1. Specifically, the D/H ratio of the planets evolving to R_p ≈ 1.8–2 R_⊕ can rise to ≈2.5, 3.0, and 3.3 × 10⁻⁵ or increase by a factor of 1.25, 1.5, and 1.65 at ages of 2.5, 5, and 7.5 Gyr, 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 to be higher than those shown in the bottom panels of Figure 2, as long as these planets do not evolve to bare rocky cores.

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 ≥2.4 × 10⁻⁵ to be deuterium-enriched, which corresponds to an increase of 20% or 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–5 M_⊕) at ≈0.25–0.4 au or a high-mass sub-Neptune (10–15 M_⊕) at ≈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.

Zoom In Zoom Out Reset image size

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 of 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 a 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 ≤ 2.5 R_⊕ 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 ≈ 1.8 Gyr in the simulations, Figure 4 presents the results in 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.6 R_⊕ ≲ R_p ≲ 2 R_⊕) 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). ⁹

Zoom In Zoom Out Reset image size

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–10 M_⊕), R_p (1.5–1.8 R_⊕), 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⁸ yr (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 a 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.

4.2. Photoevaporation Models

4.3. Core-powered Mass-loss Model

4.4. Deuterium Homopause

In this work, r₀ 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 }}_{{\rm{H}},\mathrm{He}}$. In principle, the deuterium homopause does not necessarily lie at the helium homopause in a hydrogen-rich atmosphere. Specifically, since X_D ≪ 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):

$$\begin{eqnarray}&&\displaystyle \frac{1}{{{ \mathcal D }}_{{\rm{D}},{\rm{H}},\mathrm{He}}}=\displaystyle \frac{{X}_{{\rm{H}}}}{{{ \mathcal D }}_{{\rm{D}},{\rm{H}}}}+\displaystyle \frac{{X}_{\mathrm{He}}}{{{ \mathcal D }}_{{\rm{D}},\mathrm{He}}}.\end{eqnarray} \tag{ 14 }$$

Hence, the location of the deuterium homopause is determined by the relation ${K}_{{zz}}={{ \mathcal D }}_{{\rm{D}},{\rm{H}},\mathrm{He}}$ rather than ${K}_{{zz}}={{ \mathcal D }}_{{\rm{H}},\mathrm{He}}$. With Equation (14), we can estimate the difference between ${{ \mathcal D }}_{{\rm{D}},{\rm{H}},\mathrm{He}}$ and ${{ \mathcal D }}_{{\rm{H}},\mathrm{He}}$ at a given altitude (i.e., a given temperature and pressure) to assess our simplified approach. To estimate ${{ \mathcal D }}_{{\rm{D}},{\rm{H}}}$ and ${{ \mathcal D }}_{{\rm{D}},\mathrm{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 }}_{{\rm{D}},{\rm{H}},\mathrm{He}}/{{ \mathcal D }}_{{\rm{H}},\mathrm{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⁷ to 10¹¹ cm² 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

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 the D/H ratio of photoevaporating sub-Neptunes. The key results are summarized as follows.

• 1.  
  We derive the critical diffusive escape rate for deuterium Φ_crit,D in a hydrogen- and helium-rich atmosphere. The deuterium ceases to escape when the atmospheric escape rate becomes smaller than Φ_crit,D.
• 2.  
  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–5 M_⊕) at d ≈ 0.25–0.4 au or a high-mass sub-Neptune (10–15 M_⊕) at d ≈ 0.1–0.25 au can increase the D/H values by 20% or more over 7.5 Gyr.
• 3.  
  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 ≲1.65 within 7.5 Gyr in our grid of evolutionary calculations, independent of the initial D/H value.
• 4.  
  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.

Zoom In Zoom Out Reset image size

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.

We adopt the assumption that the lowermost subsonic regions of the wind control differential escape (Zahnle & Kasting 1986; Zahnle et al. 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)

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{{X}_{j}}\displaystyle \frac{{{dX}}_{j}}{{dr}}=-\displaystyle \frac{{{GM}}_{p}({m}_{j}-{m}_{{\rm{H}}})}{{{kTr}}^{2}}+\displaystyle \sum _{i}\displaystyle \frac{{r}_{0}^{2}}{{r}^{2}{b}_{i,H}}\\ \times \,\left(\displaystyle \frac{{X}_{i}}{{X}_{{\rm{H}}}}{\phi }_{{\rm{H}}}-{\phi }_{i}\right)+\displaystyle \sum _{i}\displaystyle \frac{{r}_{0}^{2}}{{r}^{2}{bi},j}\left({\phi }_{i}-\displaystyle \frac{{X}_{i}}{{X}_{j}}{\phi }_{j}\right),\end{array}\end{eqnarray} \tag{ A1 }$$

where X is the mixing ratio, ¹⁰ 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/(μ_ij k_ij ), where μ_ij is the reduced mass, and k_ij is the collision rate. For a nonnegligible flow for the species j, its abundance would be much more uniform over one density scale height of H, i.e., $d\mathrm{ln}{X}_{j}/{dr}\ll d\mathrm{ln}{n}_{{\rm{H}}}/{dr}$. It follows from the derivation of Equation (A1) that the term $d\mathrm{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₀ = R_h (Hu et al. 2015; Malsky & Rogers 2020; Malsky et al. 2023):

$$\begin{eqnarray}&&{\rm{\Phi }}\approx {{\rm{\Phi }}}_{{\rm{H}}}+{{\rm{\Phi }}}_{\mathrm{He}}=4\pi {R}_{h}^{2}({\phi }_{{\rm{H}}}{m}_{{\rm{H}}}+{\phi }_{\mathrm{He}}{m}_{\mathrm{He}}),\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\phi }_{\mathrm{He}}}{{X}_{\mathrm{He}}}=\displaystyle \frac{{\phi }_{{\rm{H}}}}{{X}_{{\rm{H}}}}-{\phi }_{\mathrm{DL},\mathrm{He}}.\end{eqnarray} \tag{ A3 }$$

Equation (A3) is derived from Equation (A1) by neglecting $d\mathrm{ln}{X}_{\mathrm{He}}/{dr}$, using X_H + X_He ≈ 1, and employing the diffusion-limited escape flux for He given by

$$\begin{eqnarray}&&{\phi }_{\mathrm{DL},\mathrm{He}}=\displaystyle \frac{{{GM}}_{p}({m}_{\mathrm{He}}-{m}_{{\rm{H}}}){b}_{{\rm{H}},\mathrm{He}}^{{\prime} }}{{R}_{h}^{2}{{kT}}_{h}},\end{eqnarray} \tag{ A4 }$$

where T_h is the homopause temperature, and ${b}_{{\rm{H}},\mathrm{He}}^{{\prime} }$ is the binary diffusion coefficient for H and He taking into account the H ionization. Moreover, the R_h for H and He is determined by the location where the binary diffusivity ${{ \mathcal D }}_{{\rm{H}},\mathrm{He}}$ equals the eddy diffusivity K_zz , for which a nominal value of 10⁹ cm² s⁻¹ adopted by Malsky & Rogers (2020) is used. The ${{ \mathcal D }}_{{\rm{H}},\mathrm{He}}$ is determined by Fuller's method (Fuller et al. 1966; Tang et al. 2014; Malsky & Rogers 2020),

$$\begin{eqnarray}&&{{ \mathcal D }}_{{\rm{H}},\mathrm{He}}=\displaystyle \frac{{10}^{-3}{T}_{h}^{1.75}{\left(\tfrac{{m}_{{\rm{H}}}{m}_{\mathrm{He}}}{{m}_{{\rm{H}}}+{m}_{\mathrm{He}}}\right)}^{-1/2}}{{P}_{h}{\left({V}_{{\rm{H}}}^{1/3}+{V}_{\mathrm{He}}^{1/3}\right)}^{2}},\end{eqnarray} \tag{ A5 }$$

where P_h is the gas pressure in units of atmospheres 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 μ, 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).

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

$$\begin{eqnarray}&&{\rm{\Phi }}={f}_{r}{{\rm{\Phi }}}_{\mathrm{EL}},\end{eqnarray} \tag{ B1 }$$

where ${{\rm{\Phi }}}_{\mathrm{EL}}={L}_{\mathrm{EUV}}\eta {a}^{2}{R}_{h}^{3}/(4{{Kd}}^{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 τ described by the power law ${\mathrm{log}}_{10}({L}_{\mathrm{EUV}}/{\rm{J}}\,{{\rm{s}}}^{-1})\,=22.12\mbox{--}1.24{\mathrm{log}}_{10}(\tau /\mathrm{Gyr})$ for GK dwarfs at an age of t ≥ 10⁸ yr (Sanz-Forcada et al. 2011). When t < 10⁸ yr, the stellar high-energy flux is almost saturated (Owen & Wu 2017, and references therein). Hence, the L_EUV for t < 10⁸ yr is set to be the same as the power law of L_EUV at a stellar age of 10⁸ yr. Moreover, in Equation (B1), the heating efficiency η = 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\mathrm{ln}{X}_{j}/{dr}$ (thermal energy) according to the DSMC, which corrects the overestimate of Φ due solely to Φ_EL (Johnson et al. 2013; Hu et al. 2015).

In our fiducial model, the free parameters η = 0.1, K_zz = 10⁹ cm² s⁻¹, and T_h = 10,000 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⁹ and 10¹¹ cm² s⁻¹ and homopause temperature T_h = 3000 K. We also consider η = 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).

Zoom In Zoom Out Reset image size

Overall, the D/H values are insensitive to these free parameters, except in extreme cases for photoevaporation with η = 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 η = 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_⊕ 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 ∼ 1 Gyr due to the weak EUV flux. Hence, f_r becomes unity; thus, the magnitude of η becomes influential during the late stage. When η = 0.08 and 0.1, deuterium can stop escaping during this late phase due to Φ < Φ_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 the D/H ratio due to the faster diffusive escapes of both D and H when η = 0.75, even though Φ has not reached Φ_D,crit. This effect leads to the slightly large departure of the D/H values among different η in the cases of f_env = 0.005 and 0.1 at t ≳ 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 η = 0.75 and 0.08 is similar to that with η = 0.1. This may infer that our results are insensitive to η = 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 η 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 = 5 M_⊕, f_env = 0.005, and η = 0.75 is disregarded due to the unphysical radius anomaly.