The Diffusion Limit of Photoevaporation in Primordial Planetary Atmospheres

Diffusion-limited flux is the maximum rate at which a species can traverse through an atmosphere when advection and diffusion are considered concurrently. Despite its name, diffusion-limited mass loss does not imply a lack of advection. Instead, it describes how momentum diffusion between species acts against bulk flow in an advecting binary gas mixture. At the microscopic scale, the fast moving light species collide with slower moving heavy species, resulting in the light species decelerating and the heavy species accelerating. At the macroscopic scale, it appears as if the heavy species exert drag on the lighter species while the lighter species pull on the heavy species. The diffusion-limited mass loss model extends the Euler equations by incorporating interactions between different species (Chapman & Cowling 1970, p. 107). In the traditional approach of Hunten (1973), the diffusion-limited flux is estimated by evaluating the relative average velocity of a binary gas mixture in a gravitational field (Chapman & Cowling 1970),

$$\vec{v}_{1}-\vec{v}_{2}=-D\left[\frac{n^{2}}{n_{1}n_{2}}\nabla\left(\frac{n_{1}}{n}\right)+\frac{\mu_{2}-\mu_{1}}{\bar{\mu}}\nabla\left(\ln{P}\right)+\alpha_{T}\nabla\left(\ln{T}\right)-\frac{\mu_{1}\mu_{2}}{\bar{\mu}k_{\rm B}T}\left(\vec{a}_{1}-\vec{a}_{2}\right)\right],$$

(1)

where subscripts 1 and 2 are for the light and heavy constituents, respectively, $D$ is the diffusion coefficient, $n$ is the particle number density, $\bar{\mu}$ is the mean molecular mass, $P$ is the pressure, $\alpha_{T}$ is the thermal diffusion factor (not to be confused with the coefficient of thermal diffusivity; Leuenberger & Lang 2002), $T$ is the temperature, $k_{\rm B}$ is the Boltzmann constant, and $a$ is the acceleration acting on the particles from external forces. The first term is the concentration gradient, the second is the mass gradient, the third is the temperature gradient, and the fourth is the external force gradient. Equation 1 is evaluated for Earth at the homopause, which is very close to the mesopause, where there is a temperature inversion and $\nabla\left(\ln{T}\right){=}0$. The homopause is the level below which an atmosphere is well-mixed whereas the mesopause is the boundary between the mesosphere and the thermosphere; the homopause and the mesopause are both located at approximately $85~{}{\rm km}$, though the homopause is known to vary in altitude from $80{-}120~{}{\rm km}$ (Salinas et al. 2016; Liu 2021; Swenson et al. 2021).

The planetary magnetic field is small, so the external acceleration for both particles is set by gravity, which is $g{=}GM_{\rm p}/r^{2}$ ($G$ is the gravitational constant, $M_{\rm p}$ is the planetary mass, and $r$ is the radial distance) and therefore cancels out. In its current form, Equation 1 only includes molecular diffusion, which involves the random movement of individual molecules from areas of high concentration to low concentration. In contrast, eddy diffusion involves mixing by turbulence and acts on the compositional gradient term (Catling & Kasting 2017), thus it has to be included in its numerator. Introducing the mole fraction $\chi$ so that $n_{1}{=}\chi n$ and $n_{2}{=}\left(1{-}\chi\right)n$, and considering only the radial profile, we thus have

$$v_{1}-v_{2}=-D\left[\frac{1+\frac{K_{\rm zz}}{D}}{\chi\left(1-\chi\right)}\frac{{\rm d}\chi}{{\rm d}r}+\frac{\mu_{2}-\mu_{1}}{\bar{\mu}}\frac{1}{P}\frac{{\rm d}P}{{\rm d}r}\right].$$

(2)

According to Hunten (1973), the compositional gradient term is negligible. Moreover, gas is assumed to be in hydrostatic equilibrium at the homopause, so $1/P\left({\rm d}P/{\rm d}r\right)$ is equal to ${-}1/H$, where $H$ is the scale height defined as $H{=}k_{\rm B}T/\left(g\bar{\mu}\right)$,

$$v_{1}-v_{2}=\frac{Dg}{k_{\rm B}T}\left(\mu_{2}-\mu_{1}\right).$$

(3)

In the limit when $\mu_{1}{\ll}\mu_{2}$, Graham’s law of diffusion suggests that $v_{1}{\gg}v_{2}$, and $v_{1}{-}v_{2}{\approx}v_{1}$. Instead of employing Chapman-Enskog theory (Chapman & Cowling 1970) for calculating the molecular diffusion coefficient, Hunten (1973) substitutes $D$ with the binary diffusion coefficient, $b_{\rm ij}{=}nD$, which is determined experimentally (Marrero & Mason 1972) and takes the form $b_{\rm ij}{=}AT^{B}$, with $A$ and $B$ being constants. Multiplying both sides of Equation 3 by the number density of the light species, $n_{1}$, yields the diffusion-limited flux,

$$\phi_{1}=n_{1}v_{1}=\frac{\chi b_{\rm ij}g}{k_{\rm B}T}\left(\mu_{2}-\mu_{1}\right),$$

(4)

where $\phi_{1}$ is the escaping particle flux. Mass conservation ensures equal mass flow through all atmospheric spherical shells, thus the global mass loss rate is set by mass flow at the homopause (subscript H),

$$\dot{M}_{1}=\chi_{\rm H}\frac{4\pi GM_{\rm p}\mu_{1}b_{\rm ij,H}}{k_{\rm B}T_{\rm H}}\left(\mu_{2}-\mu_{1}\right).$$

(5)

Equation 5 assumes a binary gas mixture, with a light major component (atomic hydrogen) and a heavy minor component. Atmospheres are, however, composed of various species, and it is therefore necessary to define the mean molecular mass of the average effective heavy component,

$$\mu_{2}=\frac{\bar{\mu}-\chi\mu_{1}}{1-\chi},$$

(6)

where $\bar{\mu}$ is the local mean molecular mass. Collectively, one gets,

$$\dot{M}_{1}=\frac{\chi_{\rm H}}{1{-}\chi_{\rm H}}\frac{4\pi GM_{\rm p}\mu_{1}b_{\rm ij,H}}{k_{\rm B}T_{\rm H}}\left(\bar{\mu}-\mu_{1}\right).$$

(7)

Equations 5 and 7 apply to planets that have achieved a diffusion-limited steady state (section 3.1), such as Earth (Joshi et al. 2019). These equations are independent of X-ray and ultraviolet irradiation because they describe the maximum rate at which hydrogen can be transported within the atmosphere, thus setting the limit for how much hydrogen can be lost at the top of the atmosphere. As mass loss is evaluated at the homopause, and as the homopause is neutral, the effects of photochemistry are unimportant.

The crossover mass is defined as the threshold mass above which particles are too heavy to be dragged along with other escaping species. Fractionation occurs when the crossover mass is low enough for a gas mixture to undergo differential separation. Hunten et al. (1987) finds the crossover mass by solving the diffusion-limited flux for $\mu_{\rm cr}$ (${=}\mu_{2}$ in Equations 4 and 5),

$$\mu_{\rm cr}=\mu_{1}+\frac{k_{\rm B}T\phi_{1}}{\chi b_{\rm ij}g},$$

(8)

or, equivalently,

$$\mu_{\rm cr}=\mu_{1}+\frac{k_{\rm B}T\dot{M}_{1}}{4\pi\chi b_{\rm ij}GM_{\rm p}\mu_{1}}.$$

(9)

Equations 8 and 9 derive from Equation 1, which assumes that both gas species experience the same acceleration as the bulk fluid, that is, ${\rm d}v_{1}/{\rm d}t{=}{\rm d}v_{2}/{\rm d}t{=}{\rm d}v/{\rm d}t$ (Chapman & Cowling 1970, p. 107). However, this assumption does not apply when evaluating fractionation because, under such circumstances, ${\rm d}v_{2}/{\rm d}t{=}0$. To address this, we use Equation 6.62, 9 of Chapman & Cowling (1970) instead, with ${\rm d}v_{2}/{\rm d}t$ in place of ${\rm d}v/{\rm d}t$:

$$\rho_{2}\frac{{\rm d}v_{2}}{{\rm d}t}=-\rho_{2}g-\frac{{\rm d}P_{2}}{{\rm d}r}+\frac{\rho_{1}\rho_{2}}{\rho\tau_{1,2}}\left(v_{1}-v_{2}\right),$$

(10)

where $\tau_{1,2}$ is the time between particle collisions. From Equations 6.62, 6 and 6.62, 7 of Chapman & Cowling (1970),

$$\tau_{1,2}=\frac{\rho_{1}\rho_{2}P}{P_{1}P_{2}\rho}D,$$

(11)

which upon inserting into Equation 10 and simplifying with the ideal gas equation,

$$\frac{{\rm d}v_{2}}{{\rm d}t}=-g-\frac{c_{2}^{2}}{P_{2}}\frac{{\rm d}P_{2}}{{\rm d}r}+\frac{\chi c_{2}^{2}}{D}\left(v_{1}-v_{2}\right),$$

(12)

where $c_{2}{=}\sqrt{k_{\rm B}T/\mu_{2}}$ is the isothermal sound speed of the heavy species. Applying chain rule and dividing through by $c_{2}^{2}$,

$$\frac{1}{c_{2}^{2}}\frac{{\rm d}v_{2}}{{\rm d}t}=-\frac{GM_{\rm p}\mu_{2}}{k_{\rm B}Tr^{2}}+\frac{1}{1-\chi}\frac{{\rm d}\chi}{{\rm d}r}-\frac{1}{P}\frac{{\rm d}P}{{\rm d}r}+\frac{\chi}{D}\left(v_{1}-v_{2}\right).$$

(13)

Hunten (1973) evaluates the diffusion-limited flux at the homopause, but it is unclear if this is a suitable location for evaluating the crossover mass, which depends on local system properties. In fact, the properties of the upper atmosphere differ substantially from those in the deeper regions, and, if we are interested in knowing the compositional evolution of the bulk atmosphere, the most relevant place would be the upper edge of the atmosphere where atmospheric particles escape to space. It is thus preferable to evaluate the crossover mass at the exobase (subscript x), which is defined as the location above which gas becomes rarefied and collisions no longer dominate particle dynamics. The exobase height varies from $500{-}1000~{}{\rm km}$ depending on solar activity (Emmert 2015); beyond the exobase, a lower crossover mass is less significant because the reduced particle density gives rise to greater statistical noise. This statistical noise invalidates the continuity assumption because individual particle motions are too significant for gas to be treated as a coherent flow (Oran et al. 1998; Shematovich et al. 2015). Evaluating Equation 13 at the exobase where hydrostatic equilibrium applies (Modirrousta-Galian & Korenaga submitted) reduces the $(1/P){\rm d}P/{\rm d}r$ term to the negative inverse of the scale height of the bulk fluid,

$$\frac{1}{c_{2}^{2}}\frac{{\rm d}v_{2}}{{\rm d}t}=\frac{GM_{\rm p}\left(\bar{\mu}-\mu_{2}\right)}{k_{\rm B}T_{\rm x}R_{\rm x}^{2}}+\frac{1}{1-\chi_{\rm x}}\left(\frac{{\rm d}\chi}{{\rm d}r}\right)_{\rm x}+\frac{\chi_{\rm x}}{D_{\rm x}}\left(v_{\rm 1,x}-v_{\rm 2,x}\right).$$

(14)

Fractionation occurs for species heavier than the crossover mass ($\mu_{2}{=}\mu_{\rm cr}$) when $v_{2}{=}{\rm d}v_{2}/{\rm d}t{=}0$ so that

$$\frac{GM_{\rm p}\left(\bar{\mu}-\mu_{\rm cr}\right)}{k_{\rm B}T_{\rm x}R_{\rm x}^{2}}+\frac{1}{1-\chi_{\rm x}}\left(\frac{{\rm d}\chi}{{\rm d}r}\right)_{\rm x}+\frac{\chi_{\rm x}v_{\rm 1,x}}{D_{\rm x}}=0.$$

(15)

Multiplying the numerator and denominator of the rightmost term by $n$ yields,

$$\frac{GM_{\rm p}\left(\bar{\mu}-\mu_{\rm cr}\right)}{k_{\rm B}T_{\rm x}R_{\rm x}^{2}}+\frac{1}{1-\chi_{\rm x}}\left(\frac{{\rm d}\chi}{{\rm d}r}\right)_{\rm x}+\frac{\phi_{\rm x}}{b_{\rm ij,x}}=0,$$

(16)

which can also be expressed in terms of the mass flow rate,

$$\frac{GM_{\rm p}\left(\bar{\mu}-\mu_{\rm cr}\right)}{k_{\rm B}T_{\rm x}R_{\rm x}^{2}}+\frac{1}{1-\chi_{\rm x}}\left(\frac{{\rm d}\chi}{{\rm d}r}\right)_{\rm x}+\frac{\dot{M}_{1}}{4\pi R_{\rm x}^{2}\mu_{1}b_{\rm ij,x}}=0.$$

(17)

The compositional term is small and can be discarded (Zahnle & Kasting 1986). Solving for $\mu_{\rm cr}$,

$$\mu_{\rm cr}=\bar{\mu}+\frac{k_{\rm B}T_{\rm x}\dot{M}_{1}}{4\pi b_{\rm ij,x}GM_{\rm p}\mu_{1}},$$

(18)

where we see that Equations 9 and 18 share a similar functional form but with different parameter values, and $\chi$ not being in the denominator. When a diffusion-limited steady state applies (section 3.1), the hydrogen mass flow rate, $\dot{M}_{1}$, is set by Equation 7 and cannot be replaced by alternative models such as the energy-limited (Watson et al. 1981) or hydro-based (Kubyshkina et al. 2018) approximations (discussed in section 4). As mentioned previously, diffusion-limited mass loss combines advection and momentum diffusion. This interaction creates a negative feedback mechanism in multicomponent gas mixtures where extreme mass loss becomes unsustainable because heavy species slow down escaping hydrogen, reducing the crossover mass, and leading to increased fractionation and further decreased hydrogen outflow. It would therefore be misleading to use an extreme mass loss model that assumes no fractionation because this is the very reason why they exhibit extreme mass loss in the first place. Such models inherently preclude fractionation because they do not incorporate the drag effect imposed by heavy species on escaping hydrogen. To self-consistently quantify fractionation in a multicomponent system, we therefore use Equation 7 and combine it with $b_{\rm ij}{\mathchoice{\mathrel{\raise 1.07639pt\hbox{\hbox to0.0pt{\hbox{$\displaystyle\propto$}\hss}\lower 4.03563pt\hbox{$\displaystyle\sim$}}}}{\mathrel{\raise 1.07639pt\hbox{\hbox to0.0pt{\hbox{$\textstyle\propto$}\hss}\lower 4.03563pt\hbox{$\textstyle\sim$}}}}{\mathrel{\raise 0.75346pt\hbox{\hbox to0.0pt{\hbox{$\scriptstyle\propto$}\hss}\lower 2.82494pt\hbox{$\scriptstyle\sim$}}}}{\mathrel{\raise 0.5382pt\hbox{\hbox to0.0pt{\hbox{$\scriptscriptstyle\propto$}\hss}\lower 2.0178pt\hbox{$\scriptscriptstyle\sim$}}}}}T^{3/4}$ (Mason & Marrero 1970; Marrero & Mason 1972),

$$\mu_{\rm cr}=\bar{\mu}_{\rm x}+\frac{\chi_{\rm H}}{1{-}\chi_{\rm H}}\left(\frac{T_{\rm x}}{T_{\rm H}}\right)^{\frac{1}{4}}\left(\bar{\mu}_{\rm H}-\mu_{1}\right).$$

(19)

For well-mixed gas, such as exoplanets with primordial atmospheres (i.e., $K_{\rm zz}{\gg}D$; Modirrousta-Galian & Korenaga 2023), $\bar{\mu}_{\rm x}{\approx}\bar{\mu}_{\rm H}{\approx}\bar{\mu}$ and,

$$\mu_{\rm cr}=\bar{\mu}\left[1+\frac{\chi_{\rm H}}{1{-}\chi_{\rm H}}\left(\frac{T_{\rm x}}{T_{\rm H}}\right)^{\frac{1}{4}}\left(1-\frac{\mu_{1}}{\bar{\mu}}\right)\right].$$

(20)

Equations 19 and 20 apply to all planets with compositionally stratified and well-mixed ideal gas atmospheres, respectively, that have achieved a diffusion-limited steady state (section 3.1). They suggest that in the diffusion-limited mass loss framework, the crossover mass is independent of the hydrogen outflow rate but depends on the hydrogen mole fraction, the temperature ratio to the one-fourth power, and the atmospheric composition. The photosphere (where optical depth $\tau{=}2/3$) is a suitable alternative to the homopause for exoplanets with hydrogen-rich atmospheres because it is also a temperature minimum (i.e., ${\rm d}T/{\rm d}r{=}0$; Modirrousta-Galian & Korenaga 2023). In the following section, we compare our approach to that of Hunten et al. (1987).

Our model for the crossover mass (Equations 19 and 20) and that of Hunten et al. (1987) (Equations 8 and 9) are based on different assumptions. We compare their predictions with the observations of Earth’s atmosphere, using the NRLMSIS 2.0 code (Emmert et al. 2021), which provides the average observed behavior of an atmospheric column at a given latitude, longitude, and elevation. Athens, Greece (37.9838^∘ N, 23.7275^∘ E) is chosen as the reference location because of its proximity to multiple continents and its mean annual temperature being similar to Earth’s average temperature (${\sim}288~{}{\rm K}$). From the output data, we extract the mole fraction of hydrogen (H), helium (He), oxygen (O and O₂), nitrogen (N and N₂), and argon (Ar), along with the mean molecular mass, pressure, and temperature for elevations of $0{-}1000~{}{\rm km}$. Figure 1 shows the parameter and compositional profiles adopted in this study.

Equations 8 and 19 are applied to the exobase ($500{-}1000~{}{\rm km}$), with the former using the binary diffusion coefficients of Hunten (1973), Hunten et al. (1987), and Catling & Kasting (2017). We obtain the required particle flux, $\phi_{1}$, using the parameterization $\phi_{1}{=}1.2{\times}10^{12}\left(r/6837~{}{\rm km}\right)^{-2}~{}{\rm m^{-2}~{}s^{-1}}$, which derives from observations of Earth’s thermosphere (Joshi et al. 2019). Figure 2 shows that the crossover mass of Equation 8 remains relatively constant at around $1~{}{\rm amu}$, with little variation across different binary diffusion coefficient values. In contrast, our crossover mass formulation (Equation 19) yields a value of $9{\pm}6~{}{\rm amu}$, and varies with altitude because of the changing mean molecular mass across the heterosphere, suggesting that fractionation occurs primarily at higher altitudes.

The crossover mass value we obtain when using the approach of Hunten et al. (1987), $\mu_{\rm cr}{\approx}1~{}{\rm amu}$, differs from the value given in their study, $\mu_{\rm cr}{=}2.5~{}{\rm amu}$, and in other iterations of their work (e.g., $\mu_{\rm cr}{=}2.25~{}{\rm amu}$; Catling & Kasting 2017, Table 5.3). This discrepancy can be understood through differences in the assumptions and parameter values adopted when evaluating the crossover mass. First, these studies derive the crossover mass from Equation 1, which assumes that all gas species have equal acceleration (Chapman & Cowling 1970, p. 107). This is physically inconsistent because fractionation can occur only when different species accelerate at different rates and, thus, Equation 7 should be used instead. Second, the energy-limited approximation (Watson et al. 1981) is used to estimate hydrogen flux whereas we use present-day observations of leaking hydrogen (Joshi et al. 2019). The use of the energy-limited approximation is inadequate because it assumes free advection, making it inherently incapable of evaluating fractionation. Fractionation involves differential separation resulting from momentum diffusion exchange between different species, which is not incorporated in the energy-limited approximation. Third, they assume a hydrogen mole fraction close to one ($\chi{\approx}1$), which is valid only for hydrogen-rich atmospheres and not for the present-day Earth (Figure 1).

In the context of exoplanets with primordial atmospheres, Equation 20 can be used because the hydrogen mole fraction is high across the entire atmosphere ($\chi{\approx}0.9$) and mixing from eddy diffusion is efficient throughout (Parmentier et al. 2013; Charnay et al. 2015). If compositional stratification occurs, Equation 19 should be employed at the top of the atmosphere (i.e., the exobase) instead. Equation 19 is independent of the chemical gradient in the deeper sections of the atmosphere because the crossover mass at the exobase sets the fractionation bottleneck for the entire atmosphere. Figure 3 shows that the crossover mass increases with $\chi$, suggesting that the preferential loss of hydrogen leads to a progressively more stable (i.e., more fractionating) and hydrostatic atmosphere composed of heavier species.

In the previous section we discussed how planets with hydrogen-rich upper-atmospheres are prone to compositional fractionation because the low crossover mass prevents heavy species from escaping, leading to local hydrogen depletion. However, information about this depletion is not felt instantaneously across the atmosphere because only the region immediately below this depleted layer can feel the effect of the compositional gradient. It is for this reason that the diffusion limit applies only when hydrogen flow is uniform across all atmospheric spherical shells, and only through achieving such a steady state does mass flow become limited by the maximum rate at which it can be transported within the atmosphere.

In this section, we explore and apply this concept to exoplanets with primordial envelopes experiencing photoevaporation. For mass loss to continue, there must be a preferential transport of hydrogen from deeper layers; otherwise the mole fraction of hydrogen, $\chi$, approaches zero, and mass loss ceases. This transport can occur either through advection or diffusion. If the mechanism is advection, the upper atmosphere will become exceedingly more enriched as more hydrogen is lost and more heavier species are left behind. In other words, the combination of preferential hydrogen removal at the top of the atmosphere and advection with the average composition would lead to further enrichment, preventing a steady state. Therefore, if advection were the restocking mechanism, the concentration gradient would keep growing until diffusive transport restores equilibrium. Consequently, a steady state can be reached only by balancing the preferential loss of hydrogen with the preferential restocking of hydrogen, which can occur only through diffusion. The balance between diffusion and advection determines whether chemical inhomogeneities form. If diffusion is efficient, these inhomogeneities will disperse rapidly. To evaluate whether this scenario is possible, we consider a hydrogen-rich hydrodynamic atmosphere that has just experienced preferential hydrogen loss. We follow the framework of Modirrousta-Galian & Korenaga (submitted), where it was shown that the exobase and quasi-sonic point coincide for hydrodynamic atmospheres, and we therefore refer to the highest point in the atmosphere as the quasi-sonic point (subscript QS) because we assume the atmosphere is initially hydrodynamic. This point represents the velocity maximum in the atmosphere, which, in the limit of $v{=}c_{\rm s}$, corresponds to the sonic point. Figure 4 is a schematic diagram of the configuration being examined.

After experiencing mass loss, diffusion will act on the locally formed compositional gradient within a timescale of,

$$t_{\rm df}=\frac{\delta_{\rm DL}^{2}}{2\left(K_{\rm zz}{+}D\right)},$$

(21)

where $\delta_{\rm DL}$ is the depleted layer depth. Assuming the atmosphere begins in a highly hydrodynamic state, eddy diffusion greatly exceeds molecular diffusion (i.e., $K_{\rm zz}{\gg}D$). The eddy diffusion coefficient is defined as $K_{\rm zz}{=}\langle vL_{\rm m}\rangle$ where $v{=}v_{\rm QS}$ is the bulk vertical wind velocity at the quasi-sonic point and $L_{\rm m}$ is the mixing length. The mixing length is not well understood because it is subject to highly nonlinear dynamical processes, such as gravitational wave breaking (e.g., Lindzen 1971, 1981). Experiments suggest that $L_{\rm m}{\sim}C_{\rm vK}L_{\rm c}$, where $C_{\rm vK}{\approx}0.4$ is the von Kármán constant (Pope 2000) and $L_{\rm c}$ is the characteristic length, which is sensitive to local system properties. In this context, the characteristic length is the depleted layer size, obtaining $K_{\rm zz}{\sim}\langle C_{\rm vK}v_{\rm QS}\delta_{\rm DL}\rangle$. Evaluating the above

$$t_{\rm df}>\frac{5\delta_{\rm DL}}{4v_{\rm QS}}=\frac{5}{4}t_{\rm a},$$

(22)

where $t_{\rm a}$ is the advection timescale, demonstrating that diffusion is less efficient than advection. Under such circumstances, a locally depleted layer must form because advection removes hydrogen faster than it is restocked. This will cause the hydrogen mole fraction at the quasi-sonic point to decrease, creating a depleted layer through which hydrogen diffuses. As mass loss continues, the hydrogen content in the region immediately below the depleted layer diffuses through it and is subsequently lost at the quasi-sonic point. The size of the depleted layer will therefore increase with further mass loss (Figure 4).

This gives rise to three equations for mass loss: one for advection at the quasi-sonic point,

$$\phi_{\rm QS}=\chi_{\rm QS}n_{\rm QS}v_{\rm QS},$$

(23)

one for hydrogen diffusion through the depleted layer,

$$\phi_{\rm DL}=-\left(\bar{K}_{\rm zz}+\bar{D}\right)\frac{\chi_{\rm QS}n_{\rm QS}-\chi_{\rm CB}n_{\rm CB}}{\delta_{\rm DL}},$$

(24)

and the other for the growth rate of the depleted layer,

$$\phi_{\rm CB}=\chi_{\rm CB}n_{\rm CB}\frac{{\rm d}\delta_{\rm DL}}{{\rm d}t}.$$

(25)

The parameters $\chi_{\rm CB}$ and $n_{\rm CB}$ are the mole fraction of hydrogen and the volumetric total particle number density at the bottom of the depleted layer (labeled the compositional boundary), and $\bar{K}_{\rm zz}$ and $\bar{D}$ are the average eddy and molecular diffusion coefficients,

$$\bar{K}_{\rm zz}{+}\bar{D}=\frac{1}{R_{\rm QS}-R_{\rm CB}\left(t\right)}\int^{R_{\rm QS}}_{R_{\rm CB}\left(t\right)}K_{\rm zz}\left(r\right){+}D\left(r\right){\rm d}r.$$

(26)

The above equations can be solved analytically if a plane-parallel atmosphere and constant average diffusion coefficients are assumed,

$$\bar{K}_{\rm zz}{+}\bar{D}\approx\frac{1}{R_{\rm QS}-R_{0}}\int^{R_{\rm QS}}_{R_{0}}K_{\rm zz}\left(r\right){+}D\left(r\right){\rm d}r,$$

(27)

with the lower limit set to the planetary surface ($R_{0}$). Balancing Equations 24 and 25,

$$\delta_{\rm DL}\frac{{\rm d}\delta_{\rm DL}}{{\rm d}t}=\left(\bar{K}_{\rm zz}+\bar{D}\right)\left(1-\frac{\chi_{\rm QS}n_{\rm QS}}{\chi_{\rm CB}n_{\rm CB}}\right).$$

(28)

The mole fraction and number density ratio is found by balancing Equations 23 and 24,

$$\frac{\chi_{\rm QS}n_{\rm QS}}{\chi_{\rm CB}n_{\rm CB}}=\frac{\bar{K}_{\rm zz}+\bar{D}}{v_{\rm QS}\delta_{\rm DL}+\bar{K}_{\rm zz}+\bar{D}},$$

(29)

which can be inserted into Equation 28 and integrated,

$$\int^{R_{\rm QS}{-}R_{0}}_{0}\left(\frac{\delta_{\rm DL}}{\bar{K}_{\rm zz}+\bar{D}}+\frac{1}{v_{\rm QS}}\right){\rm d}\delta_{\rm DL}=\int^{t}_{0}{\rm d}t,$$

(30)

yielding,

$$\begin{split}t&=\frac{\left(R_{\rm QS}-R_{0}\right)^{2}}{2\left(\bar{K}_{\rm zz}+\bar{D}\right)}+\frac{R_{\rm QS}-R_{0}}{v_{\rm QS}}\\ &\approx\frac{\left(R_{\rm QS}-R_{0}\right)^{2}}{2{K}_{\rm zz}}\end{split}$$

(31)

Evaluating Equation 31 for typical values relevant to Earth’s atmosphere, with $\bar{K}_{\rm zz}{\sim}10^{2}~{}{\rm m^{2}~{}s^{-1}}$ (Liu 2021), suggests that the time required for steady state to apply is geologically negligible (${\lesssim}40~{}{\rm years}$). In fact, general circulation atmospheric models suggest that exoplanets may possess significantly higher eddy diffusion coefficients because they reside in more extreme environments (e.g., Parmentier et al. 2013; Charnay et al. 2015), indicating that steady state may be established in even shorter timescales. The amount of atmospheric depletion that occurs during the onset of steady state depends on the hydrogen outflow rate at the quasi-sonic point. Even with inviscid XUV-induced photoevaporation, the overall impact is likely negligible because total mass loss requires $10^{3}{-}10^{5}$ years (Figure 5). Higher mass loss rates, however, may result in complete atmospheric loss during the steady-state onset, and we discuss this possibility in section 4.1.

Having assessed the rapid onset of steady state in planets with hydrogen-rich atmospheres, we now present several representative simulations for exoplanets experiencing mass loss. Because hydrogen loss at the quasi-sonic point has to be supplied from below, diffusion-limited mass loss at any point is constrained by the region just below it. The global diffusion-limited mass loss rate is therefore equal at all points in the atmosphere. Unlike Hunten (1973) and Hunten et al. (1987), we evaluate mass loss at the photosphere (subscript p; where optical depth $\tau{=}2/3$) because it is the temperature minimum of the atmosphere and because the compositional gradient is small ($K_{\rm zz}{\gg}D$; Modirrousta-Galian & Korenaga 2023), allowing us to use Equation 7. We replace the experimentally-derived binary diffusion coefficient, $b_{\rm ij}$, with the more general molecular diffusion coefficient and photospheric number density, $n_{\rm p}D$. To estimate the diffusion coefficient we apply the Chapman-Enskog hard sphere approximation (Chapman & Cowling 1970) with the average particle kinetic radius taken as that of an atomic hydrogen-helium mixture ($\sigma_{\rm 12}{=}2{\times}10^{-10}~{}{\rm m}$),

$$\begin{split}\dot{M}_{1}&=1.1{\times}10^{5}\frac{\chi_{\rm p}}{1{-}\chi_{\rm p}}\left(\frac{\bar{\mu}_{\rm p}}{\mu_{1}}{-}1\right)\left(\frac{\mu_{1}}{\mu_{2}}+1\right)^{\frac{1}{2}}\left(\frac{M_{\rm p}}{M_{\oplus}}\right)\left(\frac{T_{\rm p}}{1000~{}{\rm K}}\right)^{-\frac{1}{2}}~{}{\rm kg~{}s^{-1}}\\ &\approx 1.1{\times}10^{5}\frac{\chi_{\rm p}}{1{-}\chi_{\rm p}}\left(\frac{\bar{\mu}_{\rm p}}{\mu_{1}}{-}1\right)\left(\frac{M_{\rm p}}{M_{\oplus}}\right)\left(\frac{T_{\rm p}}{1000~{}{\rm K}}\right)^{-\frac{1}{2}}~{}{\rm kg~{}s^{-1}},\end{split}$$

(32)

where $\left(\mu_{1}/\mu_{2}{+}1\right)^{1/2}{\approx}1$. Simulations suggest that eddy diffusion is consistently greater than molecular diffusion in hydrogen-rich atmospheres (Modirrousta-Galian & Korenaga 2023), so $\chi_{\rm p}$ can be approximated with $M_{\rm H}/\left(M_{\rm H}{+}M_{\rm z}\right)$, where $M_{\rm H}$ and $M_{\rm z}$ are the total atmospheric hydrogen and heavy components. For non-gas giant planets, the heavy component originates mainly from outgassed carbon dioxide during magma ocean solidification (Elkins-Tanton 2008; Lebrun et al. 2013; Salvador et al. 2017; Bower et al. 2019), with water being released more gradually because of its high solubility in magma (Lichtenberg et al. 2021; Miyazaki & Korenaga 2022). The initial total non-hydrous volatile budget of silicate planets is thought to be of the order ${\sim}0.01\%$ by weight (Hirschmann & Dasgupta 2009; Marty 2012) so that $M_{\rm z}{\sim}10^{-4}M_{\rm p}$. The initial accreted primordial atmosphere is less well known, with simulations providing a wide range of estimates (Lee & Chiang 2015; Ginzburg et al. 2016; Mordasini 2020), so we set it as a free parameter. We therefore approximate the mean molecular mass as,

$$\bar{\mu}=\max{\left[2.2~{}{\rm amu},\chi_{\rm p}\mu_{\rm H}{+}(1{-}\chi_{\rm p})\mu_{\rm CO_{2}}\right]},$$

(33)

with $\mu_{\rm H}$ and $\mu_{\rm CO_{2}}$ being the molecular masses of atomic hydrogen and carbon dioxide. Our test planet is GJ 357 b, which is a super-Earth with a mass, radius, and equilibrium temperature of $1.84{\pm}0.31~{}M_{\oplus}$ (Luque et al. 2019), $1.20{\pm}0.06~{}R_{\oplus}$ (Oddo et al. 2023), and $525~{}{\rm K}$ (Luque et al. 2019) respectively. The X-ray luminosity and age of the host star were determined through XMM-Newton observations to be $L_{\rm x}{=}10^{18.73}~{}{\rm W}$ and $5~{}{\rm Gyr}$ (Modirrousta-Galian et al. 2020). The X-ray and ultraviolet luminosity evolution of the star are estimated through empirical relations (Penz & Micela 2008; Sanz-Forcada et al. 2011, see also Chadney et al. 2015; Johnstone et al. 2021), which we extend back to right after the T Tauri stage (age${\sim}10~{}{\rm Myr}$). For the atmospheric size we employ the model of Lopez & Fortney (2014), which provides the atmospheric radius as a function of the planetary mass, the atmospheric mass fraction, the equilibrium temperature, and the planetary age. To provide a basis for comparison with the diffusion-limited mass loss model, we also use the energy-limited approximation of Watson et al. (1981) and the hydro-based model of Kubyshkina et al. (2018). Both models assume inviscid free advection, but the latter considers thermal heating and chemistry, usually resulting in higher mass loss rates.

Figure 5 shows the mass loss evolution of GJ 357 b for initial atmospheric hydrogen budgets of $M_{\rm H}{=}0.1M_{\rm p}$ and $0.01M_{\rm p}$. The diffusion-limited model predicts lower mass loss rates than the hydro-based and energy-limited models, enabling the retention of atmospheric hydrogen over significantly longer timescales. Assuming the stellar age of $5~{}{\rm Gyr}$ is coeval with that of the planet (Modirrousta-Galian et al. 2020), GJ 357 b probably accreted an atmosphere less than ${\sim}0.16~{}M_{\rm p}$; otherwise, it would have survived photoevaporation and still existed. This model, therefore, provides tighter constraints on the mass loss evolution of GJ 357 b than the hydro-based model, which predicts an initial hydrogen reservoir less than ${\sim}21~{}{\rm M_{p}}$ (Modirrousta-Galian et al. 2020).

Having evaluated diffusion-limited mass loss (section 2.1), the crossover mass (section 2.2), and the transition to the diffusion-limited regime (section 3.1), we now discuss the self-consistency of mass loss and fractionation.

In Modirrousta-Galian & Korenaga (2023), atmospheric evaporation is categorized into three regimes. In the first regime, mass loss arises from internal energy generated during accretion; in the second regime, mass loss occurs through a combination of internal energy and incoming thermal radiation; and in the third regime, mass loss is driven by X-ray and ultraviolet irradiation. Selecting the appropriate mass loss regime depends on the thermodynamic properties of the system, which we explored in our previous paper. To keep our explanation concise and avoid reiterating what has already been discussed in Modirrousta-Galian & Korenaga (2023), we focus only on regimes one and three, which represent endmember cases when mass loss is driven by either internal energy or incoming X-ray and ultraviolet irradiation.

Regime one (also known as core-powered mass loss) is very efficient and can remove a primordial atmosphere in a few days (Modirrousta-Galian & Korenaga 2023), which is too fast for steady state to activate (${\lesssim}40~{}{\rm years}$) and diffusion-limited mass loss to apply. This suggests that fractionation is unlikely to occur as a result of geological events that bring a planet into regime one, such as giant impacts (e.g., Canup 2008; Lock et al. 2018) or late accretion (e.g., Marchi et al. 2018, 2023). In contrast, regime three (also known as XUV photoevaporation) proceeds more slowly, requiring $10^{3}{-}10^{5}$ years to remove an atmosphere from super-Earth and sub-Neptune exoplanets (Figure 5). These longer timescales allow for steady state to activate and diffusion-limited mass loss to apply. Consequently, whereas a planet in regime three may initially experience mass loss according to traditional inviscid free advection models (e.g., the energy-limited and the hydro-based models), it will rapidly shift to diffusion-limited mass loss after the onset of steady state.

Of the three regimes of atmospheric evaporation, the third is the most widely discussed because of its prolonged duration and observability through Lyman $\alpha$ spectroscopy (Linsky et al. 2010; Lecavelier des Etangs et al. 2012; Vidal-Madjar et al. 2013; Rockcliffe et al. 2023). The prevailing model employed for assessing mass loss in the third regime is the energy-limited model (Watson et al. 1981). Despite its apparent simplicity and intuitive nature, this model is beset by several well-documented limitations (e.g., Kubyshkina et al. 2018; Krenn et al. 2021; Modirrousta-Galian & Korenaga 2023). Most notably, it neglects the influence of incoming thermal radiation and erroneously assumes that incoming XUV energy is deposited at the XUV-photosphere ($\tau_{\rm XUV}{=}2/3$). These assumptions are hard to justify because exoplanets are often highly irradiated, and in the case of primordial atmospheres, the XUV photosphere is usually above the sonic point (Sekiya et al. 1980, 1981). Therefore, the sonic point (or more correctly, the quasi-sonic point; Modirrousta-Galian & Korenaga submitted) is a better location for assessing XUV energy deposition.

To overcome these limitations, Kubyshkina et al. (2018) proposed the Hydro-based model, which incorporates a more realistic chemical analysis and addresses the limitations of the energy-limited approximation. However, this model also makes several assumptions that are likely to be invalid. First, they assume that gas accelerates indefinitely, which is not always accurate (Modirrousta-Galian & Korenaga submitted). Second, it assumes that mass loss occurs through inviscid free advection for prolonged durations, and it is therefore efficient in stripping hydrogen from planetary atmospheres. Consequently, the question still remains as to why some exoplanets maintain their hydrogen-rich atmospheres while others do not. If, indeed, fractionation is significant, this could provide a plausible explanation for the observed diversity in exoplanet atmospheres because extreme hydrogen loss is unsustainable if fractionation takes place and mass loss becomes diffusion-limited. In other words, the energy-limited and hydro-based models are both limited in their ability to describe mass loss, and they should be applied judiciously. In contrast, the diffusion-limited mass loss model incorporates the effects of chemical fractionation, and it therefore provides a more realistic description of mass loss.

Our model suggests that planetary atmospheres become chemically fractionated because of XUV-induced photoevaporation, transitioning from a primordial hydrogen-rich atmosphere to a secondary heavier one. The specific details of this evolution depend on the planet’s bulk properties and its environment, requiring a case-by-case evaluation. This evolution may be detectable using atmospheric spectroscopy (Burrows 2014; Madhusudhan 2019). However, observations are limited to cloudless planets because clouds obstruct visibility into deeper atmospheric layers (Bétrémieux & Swain 2017, 2018) and create pronounced chemical stratification, likely masking the effects of diffusion-limited chemical fractionation. The recent launch of the James Webb Space Telescope and upcoming missions like Ariel and Twinkle hold promise for advancing our understanding of the link between chemical fractionation and planetary evolution.