#  Cosmic Hydrogen and Ice Loss Lines
人能常清靜，天地悉皆歸

Li Zeng Department of Earth and Planetary Sciences, Harvard University, 20
Oxford Street, Cambridge, MA 02138, USA  Stein B. Jacobsen Department of Earth
and Planetary Sciences, Harvard University, 20 Oxford Street, Cambridge, MA
02138, USA

###### Abstract

We explain the overall equilibrium-temperature-dependent trend in the
exoplanet mass-radius diagram, using the escape mechanisms of hydrogen and
relevant volatiles, and the chemical equilibrium calculation of molecular
hydrogen (H2) break-up into atomic hydrogen (H). We identify two Cosmic
Hydrogen and Ice Loss Lines (CHILLs) in the mass-radius diagram. Gas disks are
well known to disperse in ten million years. However, gas-rich planets may
lose some or almost all gas on a much longer timescale. We thus hypothesize
that most planets that are born out of a hydrogen-gas-dominated nebular disk
begin by possessing a primordial H2-envelope. This envelope is gradually lost
due to escape processes caused by host-stellar radiation.

Extrasolar Planet (251)— Interdisciplinary astronomy(804)

##  1 Introduction

Planetary systems are thought to form in a residual disk after the formation
of a star. This residual disk is thought to be of the same composition as the
star: primarily hydrogen and helium, but a small component (1%) can form ices
(H2O, NH3 and CH4) and an even smaller component (0.5%) can form silicate rock
and metals. One possible outcome would be that planets of roughly stellar
composition (gas giants) formed at various distances in the residual
protoplanetary disk. However, we know from exoplanet observations that this is
not the case, in fact gas giants makes up a relatively small portion
($\lesssim$10%) of the observed exoplanet inventory. With recent observations
of exoplanets around small M-dwarf host stars (Luque & Pallé, 2022; Cherubim
et al., 2022; Piaulet et al., 2022; Diamond-Lowe et al., 2022), as well as
atmospheric characterization of puffy exoplanets by JWST (O’Callaghan, 2022;
Alderson et al., 2023; Rustamkulov et al., 2023; Feinstein et al., 2023; Ahrer
et al., 2023), the water world hypothesis becomes viable. Likely, the most
common exoplanet ($\gtrsim$60%) appears be what we call water worlds,
consisting dominantly of H2O, NH3 and CH4 ice components (Zeng et al., 2019),
which is also consistent with the statistics on the planet radius versus flux
(Rp-f) diagram. Earth-like rock-metal planets are the second most common
($\sim$30%) of observed exoplanets. This raises the question of how almost all
of the H2-He component was lost from most planets and how in addition the ice
components were lost from Earth-like rock-metal planets. We thus infer that
there must be processes resulting in the Cosmic Hydrogen and Ice Loss Lines
that we call CHILLs. In this paper we make use of the exoplanet systematics in
the mass-radius diagram to infer, discuss, and model these processes.

We thus hypothesize that the formation of most planets greater than about one
Earth mass (1 M⊕) would start by acquiring a thick primordial hydrogen-helium
envelope while being embedded in the residual gaseous nebular disk. Then, this
primordial hydrogen-helium envelope is eroded away over time by various
mechanisms primarily caused by the host star, at various escape rates
depending on the host-stellar irradiation level and other factors. Therefore,
when we look at the distribution of exoplanets in the mass-radius diagram, we
see two Cosmic Hydrogen and Ice Loss Lines (CHILLs) of hydrogen and related
volatile escape.

The CHILL#1 may correspond to the complete removal of the primordial hydrogen-
helium envelopes from smaller planets including super-Earths and mini-
Neptunes. This removal helps to reveal the nature of their core materials
which are made of elements heavier than H/He. This can explain the observed
bi-modal radius distribution of small exoplanets separated by the _small
planet radius gap_ as core materials dominated by the elements of oxygen,
carbon, nitrogen, and their chemical compounds, versus core materials
dominated by magnesium, silicon, iron, and their chemical compounds.

On the other hand, the CHILL#2 may correspond to the complete loss of the
hydrogen-helium envelope from the more massive gas giant exoplanets, which are
dominated by hydrogen and helium in their composition. Their masses typically
range from a few tens to a few hundred Earth masses (M⊕). The current study of
their atmospheric composition may add in more constraints and inputs in their
models.

##  2 Escape Mechanism

###  2.1 Geological Evidence

Nebular gas disks typically remove on the 5-10 million years timescale
(Armitage, 2010; Hartmann, 2008). On the other hand, the noble gas abundances
sourced from the Earth’s deep mantle suggest an early existence of a massive
H2-He proto-atmosphere/envelope of the order of $\gtrsim 10^{3}$bar (0.1 GPa,
or equivalently, $10^{3}$ times the mass of the current Earth’s N2-O2
atmosphere) in equilibrium with the Earth’s early molten magma ocean surface
(Harper & Jacobsen, 1996). This massive H2-He proto-atmosphere/envelope
presumably has lasted long enough to allow the in-gassing of its helium (3He)
and neon (20Ne, 22Ne) content into the Earth’s early molten mantle (magma
ocean), which later freezes (solidifies) to become the present-day Earth’s
mantle. Later on, this massive proto-atmosphere is removed. Subsequently, a
meteoritic (planetary) atmospheric component is accreted, and is reflected by
the abundances of heavier noble gases including argon (36Ar) and xenon
(130Xe). Therefore, the geological evidence points to the existence of two
chemical reservoirs available during planet formation: the ” _solar_ ” and the
” _planetary_ ” (=” _meteoritic_ ”), corresponding to the free gas versus the
gas trapped in solid phases in the nebula, respectively. The Earth’s mantle
plume source 20Ne-22Ne ratio and 3He-20Ne ratio preserve the ” _solar_ ”-like
relative abundance ratio, while the Earth’s mantle plume source 36Ar/130Xe
ratio is close to the ” _planetary_ ” ($=$” _meteoritic_ ”) abundance ratio
(Jacobsen, 2005; Williams & Mukhopadhyay, 2019).

In order for significant in-gassing to occur over a limited time duration of
planet formation, it requires enough pressure built up at the interface
between the molten magma ocean and the overlaying gaseous layer. To the first
order approximation, the in-gassing rate is proportional to the number density
of gaseous species present, which is in turn proportional to the pressure
directly overlying the molten surface. Small bodies do not have enough gravity
to build up that pressure required for in-gassing for the limited lifetime of
the gas disk.

Thus, the overall noble gas abundances of the light versus the heavy noble gas
species paint a general two-stage picture of Earth formation: First, a co-
accretion of a massive H2-He proto-atmosphere/envelope with the proto-Earth,
and Second, the subsequent removal and replacement of such proto-
atmosphere/envelope. Likewise, other planets more massive than one Earth mass
(1 M⊕) would likely acquire a H2-He proto-atmosphere/envelope during their
formation while being embedded in a H2-He gas-rich disk (Hayashi et al., 1979;
Mizuno, 1980; Sekiya et al., 1980; Hartmann, 2008). Now, the question is the
timescale that this primordial atmosphere/envelope can last on exoplanets
under various host-stellar irradiation conditions, resulting from different
physical and initial conditions compared to our own solar system planets.

Hydrogen (H and H2) is thereby inferred as the dominant species in the upper
atmosphere of many of the close-in exoplanets, due to (1) the
availability/abundance of hydrogen and (2) its lightest weight and thus its
largest scale height compared to the other atmospheric species. In any
atmosphere above a certain critical level (the von Karman line), the gas
density becomes low enough, and the mean free path of molecular collision
becomes large enough, so that each gas species separates out according to its
own scale height determined by its own molecular weight at a given
temperature. For example, this is the case for Earth’s atmosphere, which
orbits the Sun at 1AU. This physics should also be true for exoplanets which
are closer in and simultaneously more H2-rich or H2O-rich. Therefore, we
expect the escape and dissociation of hydrogen and hydrogen-bearing species as
one of the main physical mechanisms sculpting the overall distribution of
exoplanets on the mass-radius diagram. The escape process is ultimately
powered by the intense host-stellar electromagnetic irradiation at all
wavelengths (Owen, 2019).

###  2.2 Atmospheric Layering

Under the influence of planet gravity field, a planetary atmosphere typically
exhibits layering structure. This atmospheric layering induced by
gravitational potential as well as material property of gas will have strong
influence when we consider the physics of the escape processes occurring on
these planets. For simplicity, let us first consider the thermal escape
process (see Appendix C and Appendix D) under the relatively low host-stellar
photon flux level.

The thermal escape process is typically modelled by two limiting scenarios,
(i) Jeans escape (see Appendix C), and (ii) Hydrodynamic escape (see Appendix
D) (Jeans, 1925; Chapman & Cowling, 1939; Chapman, 2000; Cercignani, 2000;
Parker et al., 1958; Hunten & Strobel, 1974; Hunten, 1971, 1973, 1993, 2002;
Chamberlain & Hunten, 1989; Walker, 1977; Goody & Walker, 1972; Goody, 1995;
Anderson & Donahue, 1975; Pudovkin & Golovchanskaya, 1983; Zahnle & Kasting,
1986; Zahnle & Catling, 2017; Catling & Kasting, 2017; Lehmer et al., 2017;
Owen, 2019; Gronoff et al., 2020; Mierkiewicz, 2021; Kang et al., 2021). For
Scenario (i) Jeans escape, the atom or molecule escapes mostly by its own
(random) thermal velocity ($v_{\text{random}}$) from the most tenuous part of
the upper atmosphere (virtually collisionless), with typically a near-
Maxwellian velocity distribution assumed for the purpose of calculation. For
Scenario (ii) Hydrodynamic escape, the atoms or molecules escape collectively
as a bulk outflow characterized by an altitude-dependent bulk outflow (mean)
velocity ($\overline{v}$) sourced from a certain region in the atmosphere,
with the flow itself driven by an _effective pressure gradient_ , powered by
certain physical mechanism such as the host-stellar irradiation.

In reality, both Scenario (i) escape due to the random thermal velocity
($v_{\text{random}}$) and Scenario (ii) escape due to the average group flow
velocity ($\overline{v}$), should occur simultaneously. Thus, they are merely
two limiting cases of the same thermally-driven escape process. Collectively
speaking, they are both considered thermal escape processes. Therefore, it is
important to identify the _source region_ and the _heating mechanism_ of this
thermal escape process.

As mentioned before, due to the planet surface gravity and the compressibility
of the atmospheric gas species, the planet atmosphere typically exhibits a
multi-layered structure with different critical number densities at different
altitude levels in the atmosphere, which correspond to different physics
operating at each level. This number density can vary across many orders of
magnitude across the thickness of the planetary atmosphere. It is also
important to view the escape process in the context of this multi-layered
atmosphere structure. According to the literature and previous works by other
researchers, we tentatively identify four important layers in the planet
atmosphere as follows: (1) the von Karman line, (2) the exobase, (3) the
transonic front, and (4) the H-H dissociation front.

(1) The first critical level in the atmosphere is the von Karman line (Von
Karman & Edson, 1967; McDowell, 2018), or the physical edge of the planetary
atmosphere. It is also a critical level below which the atmosphere is well-
mixed (Wallace & Hobbs, 2006). At this altitude level the molecular and eddy
diffusion coefficients are approximately equal. Below this level the mean free
path of molecular collisions becomes small enough ($\lambda\lesssim 1$ cm) to
allow effective mixing of different gas species so that they all have the same
scale height, which is calculated by the average molecular mass of the mixture
at a given temperature. Above this level, each gas species then assumes its
own scale height according to its own molecular weight. For example, this
level is about $\sim$80 km height above sea-level in the Earth’s atmosphere.
It is also closely related with the _homopause_ , _turbopause_ , or
_mesopause_ (boundary between mesosphere and thermosphere). It is considered
as the physical edge of the atmosphere and defines its boundary with the outer
space.

Below the von Karman line, $\lambda$ is short enough so that each atmospheric
particle may experience many collisions in its neighborhood before
transporting to other levels. Thus, Local Thermal Equilibrium (LTE) can be
safely assumed (Goody, 1995). Here $n$ is the number density (cm-3), $\sigma$
is the molecular collisional cross-section (cm2), and $\lambda$ is the mean
free path for molecular collision (cm),

$1\text{~{}cm}\sim\lambda_{\text{von Karman line}}\sim\frac{1}{n_{\text{von
Karman line}}\cdot\sigma}$ (1)

The collisional cross-section $\sigma$ is typically the atomic or molecular
dimension ($\sim 1~{}\AA$) squared:

$\sigma\sim(1\AA)^{2}\sim 10^{-16}\text{~{}cm},\text{~{}so
that~{}}n_{\text{von Karman line}}\sim 10^{16}\text{~{}cm}^{-3}.$ (2)

Generally speaking, there are two types of collisional cross-sections to
consider: low-energy elastic collisions versus energetic inelastic collisions.
So that the collisional cross-section is energy-dependent. Thus, the cross-
section for energetic in-elastic collision required for molecular dissociation
may be one order or even two orders of magnitude smaller than the low-energy
elastic collisional cross-section (Landau & Lifshitz, 1981).

In Appendix A, we explain that the H2 molecule are ”squeezy”, in fact, any
molecules are ”squeezy”. They would appear to have different collisional
cross-sections under different energy transfer situations (Landau, 1932a, b).

The typical van der waal radius for the low-energy ”elastic” collisional
cross-sections (two molecules just gently touch each other and bounce away) is
usually a few times larger, giving rise to the so-called _kinetic cross-
section_ , being about 10-100 times larger, than the _reaction cross-section_
for inelastic collisions when the energy exchange is so great (on the order of
a few electron-volts in order to rip a chemical bond).

For our purpose, we need high-energy ”inelastic collisions” that can deposit
the kinetic energy into the vibrational mode, that is the stretching of the
H-H (hydrogen) molecule or the H-O-H (water) molecule and eventually break the
X-H covalent bond. We need ”head-on” collisions, not the ”sideway” collisions.
High-energy ”head-on” ”inelastic” collisional cross-section is smaller
compared to that calculated from the gentler low-energy ”sideway” collisional
cross-section from the van der waal radius. On the other hand, sideway
collisions are capable of maintaining the rotational mode of (H-H) molecule in
equilibrium with its translational mode down to much lower number density.

(2) The second critical layer in the upper atmosphere is the exobase
(p${}_{\text{exob}}$) where the mean free path (MFP or $\lambda$) of
atom/molecule becomes so large as comparable to the local atmospheric scale
height ($\mathcal{H}\approx 10^{2}\sim 10^{3}$km) of the upper atmosphere,
which is of course temperature-dependent and species-dependent. There the
atmosphere becomes extremely tenuous and virtually collisionless (Goody,
1995),

$\mathcal{H}\sim\lambda_{\text{exobase}}\sim\frac{1}{n_{\text{exobase}}\cdot\sigma}$
(3)

Recall the definition of the molecular or atomic collisional cross-section
$\sigma$:

$\sigma\sim(1\AA)^{2}\sim 10^{-16}\text{~{}cm},\text{~{}so
that~{}}n_{\text{exobase}}\sim 10^{8}-10^{9}\text{~{}cm}^{-3}.$ (4)

(3) The third critical layer in the upper atmosphere is the sonic level (rc)
or transonic front where the outward escaping flow transitions from _sub-
sonic_ to _super-sonic_ , which is a characteristic of the flow enforced by
the continuity equation. Let’s define
$u_{0}\equiv\sqrt{k_{B}T/\mu_{\text{H}}}$ as the local isothermal sound speed
for the gas species (neutral atomic hydrogen) considered (Lehmer et al., 2017;
Walker, 1977), and let $T_{3}\equiv T/10^{3}$K, we have,

$r_{c}=\frac{GM_{p}}{2\cdot
u_{0}^{2}}=\frac{GM_{p}\cdot\mu_{\text{H}}}{2k_{B}T}=\frac{R_{p}^{2}}{2\cdot\mathcal{H}}=\bigg{(}\frac{R_{p}^{2}}{2\cdot(k_{B}T)/(\mu_{\text{H}}\cdot
g)}\bigg{)}=3.75R_{\oplus}\cdot\bigg{(}\frac{M_{p}/M_{\oplus}}{T/10^{3}\text{K}}\bigg{)}\propto\frac{M_{p}}{T_{3}}$
(5)

Due to its inverse-proportionality to Temperature T, the sonic level (rc) is
lowest and closest to the planet surface where the temperature is highest,
which occurs typically near the sub-stellar region of the planet—hot spot.

(4) The fourth critical layer to consider is the H-H bond _dissociation front_
(see the following discussion as well as Appendix E and Appendix F). It is
considered as the source region of the neutral atomic hydrogen outflow (Figure
4). Because a hydrogen atom is the lightest species among all atmospheric
species, which is only half the molecular mass of an H2-molecule, and thus, a
hydrogen atom is expected to escape into outer space more easily than an
H2-molecule under similar conditions such as temperature and planetary
gravitational field (see Appendix C and Appendix D).

Under higher host-stellar photon flux level, the atmosphere will become much
more dynamic as compared to the static case. Layering is still expected to
exist due to planet gravity, however, it geometry may change, especially
within the hot spot region where the interaction between the host-stellar
photon and planetary matter is the most intense. Some of the layers may become
less important or even irrelevant, while some of the layers, in particular,
the _sonic front_ and the _dissociation front_ may even merge to form a _shock
front_ , so that the escape process becomes much more energetic and is no-
longer diffusion-limited. We will discuss this photon-driven expansive shock
escape mechanism in Section 2.7.

###  2.3 Covalent bond

One H-H molecule can be visualized as a dumbbell consisting of two weights
(protons) held together by a spring (a covalent bond): Fig. 1. The elastic
modulus and tensile strength of this spring are two physical parameters which
characterize this molecule in study.

Figure 1: Schematic diagram of a hydrogen molecule, showing two hydrogen atoms
connected by a spring (covalent bond).

On one hand, the elastic modulus or spring constant ($k$) of this spring
determines the fundamental vibrational frequency of this molecule. For small
displacement $r$ from equilibrium, the potential energy $V$ can be
approximated as (Tsien, 1965; Lewis, 2004; Draine, 2011): $V(r)\approx
V(0)+\frac{1}{2}\times k\cdot r^{2}$, with fundamental frequency
$\nu_{0}=\frac{1}{2\pi}\cdot\sqrt{k/m_{r}}$, with $m_{r}$ being the _reduced
mass_ of this molecule $m_{r}=m_{1}\cdot m_{2}/(m_{1}+m_{2})$. Quantum-
mechanically, this molecule would have vibrational energy mode spacing as
follows:

$E_{n}=\left(n+\frac{1}{2}\right)\cdot h\cdot\nu_{0},\text{~{}where
$n=0,1,2,3,...$}$ (6)

On the other hand, its tensile strength determines the _maximum strength_
($\sim\frac{1}{2}k\cdot{r_{\text{max}}}^{2}$) of this spring beyond which it
is broken, where $r_{\text{max}}$ being the maximum displacement. Then, the
molecule become dissociated into its constituent atoms and fly apart.
Interestingly enough, this bond strength of a H-H bond is very similar to that
of an O-H bond, or even a N-H bond or a C-H bond (Kim & Le, 2020;
Karapetyants, 1978; Tsien, 1965), all being at about 4.5 eV.

Both $k$ and $r_{\text{max}}$ stay roughly the same for any X-H single
covalent bond. The main difference comes about from the _reduced mass_
$m_{r}$. Hydrides (X-H) have the smallest reduced masses, with H2 molecule
(H-H) being the extreme. Thus, H2 has an un-usually high $\nu_{0}\approx
1.5\times 10^{14}$ Hz corresponding to a wavelength $\lambda\approx 2\mu$m, or
an energy $h\nu_{0}\approx 0.6$ eV, which is equivalent to the energy of a
near-infrared photon (Lewis, 2004; Draine, 2011).

Thus, for H2 molecule, there exist about $10-15$ vibrational energy levels
(counting from ground state up and taking into account _anharmonicity_ which
diminishes the spacing of higher-lying energy levels), before it is fully
dissociated into hydrogen atoms. Therefore, H2 is the most quantum-mechanical
and least classical of all molecules. For other molecules, the spacing between
energy levels becomes tighter due to its larger _reduced mass_ $m_{r}$, and
thus, smaller $\nu_{0}$, so that more vibrational energy levels can fit in
before molecular dissociation.

Just as in a _tug-of-war_ game in which two teams pull at opposite ends of a
rope, if the rope breaks all of a sudden, then the two teams would usually not
stay still but would fall apart with extra kinetic energy.

Likewise, if a bond breaks or fractures, then its constituents would typically
not stand still but would fly apart with extra kinetic energy corresponding to
the same order of magnitude as the bond strength (Fig. 2).

Figure 2: Schematic diagram of a hydrogen molecule dissociation, showing bond
breakage and constituents flying apart.

If it is the same physical mechanism which breaks the covalent bond that also
launches the H-atom into space, then, one can estimate the fly-apart velocity
of escaping atom with the following _back-of-envelope calculation_ using
energy conservation:

$m_{\text{H-atom}}\sim m_{\text{proton}}\sim 1.67\times 10^{-27}\text{~{}g}$
(7)

$2\cdot\left(\frac{1}{2}\cdot m_{\text{H-atom}}\cdot v^{2}\right)\sim
4.5\text{~{}eV}$ (8)

$v\sim 20\text{~{}km/s}$ (9)

This value ($\sim$20 km/s) falls right within the ballpark of the observed
range of ($\sim$10-30 km/s) of the high-speed neutral H-atom outflow observed
on several exoplanets (Zhang et al., 2022, 2023). This is a first piece of
supportive evidence towards our hypothesis of escape due to _molecular
dissociation_.

The discussion here also applies to hydrogen atom sourcing from the
dissociation of other hydrides X-H such as H2O molecule. However, in those
scenarios, due to the heavy mass of X compared to H and _conservation of
momentum_ , the single H-atom after dissociation would carry away most of the
kinetic energy.

###  2.4 Collisional Dissociation

So far we have discussed the dissociation of a H2 molecule from the point-of-
view of a single molecule. Now, we will place the molecule into a
thermodynamic setting as surrounded by numerous other neighboring molecules
constantly interacting and exchanging energy with each other, and discuss the
physical and chemical process from the point-of-view of an assembly, where
both temperature and pressure can be well-defined.

The dissociation of a hydrogen molecule into two hydrogen atoms is favored at
increasing temperature or decreasing pressure within the assembly. The
relative proportion between the neutral atomic hydrogen (H-atom) and the
molecular hydrogen (H2) will change at different altitude levels in the
atmosphere. The dissociation of H-H bond involves a reversible chemical
reaction under the consideration of collisional equilibrium (see Appendix E):

$\text{H}_{2}\rightleftharpoons\text{H}+\text{H}$ (10)

Assuming Local Thermodynamic Equilibrium (LTE) and Local Chemical Equilibrium
(LCE), the relative proportion between the number of hydrogen atoms and the
number of hydrogen molecules can be estimated from the partial pressure
equilibrium constant ($K_{p}$) according to the _law of mass action_ as:

$K_{p}(T)=\frac{(p_{\text{H}})^{2}}{p_{\text{H}_{2}}}$ (11)

$p$’s are partial pressures of each species are in the unit of atm
($\sim$bar=$10^{5}$Pa).

Equivalently, $K_{p}$ can be expressed in terms of the degree of dissociation
$\alpha$, and the total pressure $P$ as:

$K_{p}(T)=\frac{((2\alpha)\cdot nRT)^{2}}{(1-\alpha)\cdot
nRT}=\frac{(2\alpha)^{2}}{(1-\alpha)}\cdot\bigg{(}nRT\bigg{)}=\frac{(2\alpha)^{2}}{(1-\alpha)}\cdot\bigg{(}\frac{P}{1+\alpha}\bigg{)}=\frac{4\alpha^{2}}{(1-\alpha^{2})}\cdot
P$ (12)

where $n$ is the number density of total H2 _originally_ present, $P$ is the
total pressure of the mixture under _current_ equilibrium condition:
$P=p_{\text{H}}+p_{\text{H}_{2}}$, and $\alpha$ is the degree of dissociation
defined as:

$\alpha=\frac{p_{\text{H}_{2}\text{(decomposed)}}}{p_{\text{H}_{2}\text{(total
H${}_{2}$ originally present)}}}$ (13)

For the temperature range from 1000 K to 4000 K, the value of $K_{p}$ is
approximately (Lewis, 2004):

$\log_{10}K_{p}(T)=6.16-\frac{23,500}{T(\text{K})}=6.16-\frac{23.5}{T_{3}}$
(14)

The first term (6.16) reflects the entropy difference between the products and
the reactants measured in the appropriate unit, whereas the numerator above
the temperature in the second term basically reflects the bond strength of an
H-H bond measured in a temperature scale, which is about 4.5 eV, or 2.25 eV
per H-atom liberated from this H-H chemical bond. This energy (enthalpy) is
much higher compared to the typical thermal energy of an atom or molecule in
the temperature range ($\sim 10^{3}$ K) that we are interested. Therefore,
only the most energetic molecules in the tail-end of the Maxwellian velocity
distribution can participate in this dissociation reaction.

One one hand, a higher total pressure ($P=p_{\text{H}}+p_{\text{H}_{2}}$)
would favor a shift towards more H2-molecules or more associated hydrogen-
bearing molecules (hydrides) X-H within the assembly.

On the other hand, a higher temperature would favor a shift towards more
dissociated molecules, fractals, and single hydrogen atoms within the
assembly. The _effective_ collisional cross-sections ($\sigma$’s) of all the
different species participating in the dissociation reaction network will
definitely affect the _rate_ of the dissociation reaction in either the
forward or the backward direction.

###  2.5 Thermal Dissociation Line of H2 Molecule

Below $\sim$2000 K, hydrogen is present almost exclusively as molecular H2
(Lewis, 2004). The _transition_ from hydrogen-molecule-dominated regime to
hydrogen-atom-dominated regime occurs along a narrow zone in the P-T diagram.
This narrow zone can be approximated as a narrow curve in the lgP-lgT diagram.
Let’s now find such a curve.

$\log_{10}\left[\frac{4\alpha^{2}}{(1-\alpha^{2})}\cdot
P_{\text{atm}}\right]=6.16-\frac{23.5}{T_{3}}$ (15)

For small $\alpha$, we have:

$\log_{10}(4\alpha^{2})+\log_{10}P_{\text{atm}}\approx
6.16-\frac{23.5}{T_{3}}$ (16)

so,

$\log_{10}(\alpha^{2})\approx
6.16-\log_{10}{4}-\log_{10}P_{\text{atm}}-\frac{23.5}{T_{3}}$ (17)

or,

$\log_{10}\alpha\approx
2.78-\frac{1}{2}\cdot\log_{10}P_{\text{atm}}-\frac{11.75}{T_{3}}$ (18)

Consider the temperature-dependent critical pressure
($P_{\text{diss1\%H${}_{2}$}}$(T)) at which dissociation reaches one-percent
level ($\alpha=1$%). It traces out a curve in P-T space (Fig. 3). Along this
curve approximately 1% hydrogen molecules are dissociated into hydrogen atoms
under equilibrium condition. Equivalently, $P_{\text{diss1\%H${}_{2}$}}$ can
be expressed as:

$\log_{10}P_{\text{diss1\%H${}_{2}$}}(\text{atm})\approx
9.56-\frac{23.5}{T_{3}}$ (19)

Alternative pathways of forming atomic hydrogen through photo-dissociation:
the direct dissociation of hydrogen-bearing molecule species by the incoming
energetic UV photons from the host-star are further explained in Appendix F.
If those hydrogen atoms formed from photo-dissociate make frequent enough
collisions with other atmospheric species before their escape, as they have
cascaded their energy to other species we can consider them being
_thermalized_ and include them in our discussion of thermal escape here.

In the lower atmosphere, collisions are so frequent as to ensure the local
thermodynamic equilibrium (LTE), which then allows the exchange of energy
across different energy modes of a molecule—translational, vibrational,
rotational, through frequent collisions. This is indeed needed to convert the
translational kinetic energy and deposit it through collisions into the
molecular vibrational energy in order to break up the H-H bond. Due to the
symmetry of a diatomic H2-molecule, it has zero dipole-moment, and thus, is
almost transparent to the electromagnetic radiation field by itself in the
visible or infrared wavelengths of either stellar or planetary origin (see
Appendix G). However, the presence of trace amount of other polar molecules
such as H2O can have strong absorptions and transport the energy to H2 through
collision cascades.

In the upper atmosphere, collisions become rarer so that a liberated hydrogen
atom can deviate from its collisional/chemical equilibrium with H2 and other
hydrogen-bearing species, and easily make escape. Due to the reduced
collisional frequency, the atmospheric temperature above the von Karman line
also starts to deviate from the planet equilibrium temperature, and becomes
highly sensitive and variable as influenced by the host-stellar electro-
magnetic activities and also space environments. Thus, this is another reason
for us to evaluate the thermal escape process right at the von Karman line,
which is the physical boundary between the atmospheric body and outer space,
where we can still rely upon the temperature calculated from the planet
equilibrium temperature calculation (see Appendix B).

Under the more energetic and dynamic settings with higher incident stellar
photon flux, the planetary atmospheric layering geometry considered under the
more static scenario in Section 2.2 may be _warped_ significantly within the
hot spot region. The actual escape geometry and mechanism may then involve a
local _shock front_ formed from the merging of two or more critical layers in
the atmosphere driven by the strong interaction between the incoming stellar
photons and the planetary surface matter, which we will discuss in detail in
the following Section 2.7.

However, the discussions in Section 2.2 and the calculations in Appendix D and
Appendix H are still valuable and meaningful because those arguments boil down
to the general principle of energy conservation (_Bernoulli’s Principle_ for
flow and the exploitation of the degree-of-freedom of molecular dissociation),
such that one can evaluate the temperature-dependence of the escape rate near
the source.

###  2.6 Thermal Dissociation Line of Water Molecule

Likewise, the dissociation of a water molecule into hydrogen and oxygen (or OH
radical) proceeds at increasing temperature or decreasing pressure. This
chemical equilibrium has been experimentally investigated by _Walther Nernst_
and his colleagues (Jellinek, 1986; Karapetyants, 1978). For not-too-high-
temperatures where further dissociation into H atoms and OH radicals does not
play important role, we have:

$2~{}\text{H}_{2}\text{O}(\text{g})\rightleftharpoons
2~{}\text{H}_{2}(\text{g})+\text{O}_{2}(\text{g})+57~{}\text{kcal/(mol of
H${}_{2}$O)}$ (20)

This enthalpy of 57 kcal/(mol of H2O) translates to 114 kcal/(mol of
H2O-pairs), which is $\sim$4.9 eV per pair of H2O from the stoichiometry of
the formula above. The partial pressure (in unit of atm) equilibrium constant
$K_{p}$ is:

$K_{p}(T)=\frac{(p_{\text{H${}_{2}$}})^{2}\cdot(p_{\text{O${}_{2}$}})}{(p_{\text{H${}_{2}$O}})^{2}}$
(21)

Equivalently, $K_{p}$ can be expressed in terms of the degree of dissociation
$\alpha$, and the total pressure $P$ as:

$K_{p}(T)=\frac{(\alpha\cdot nRT)^{2}\cdot(\frac{1}{2}\alpha\cdot
nRT)}{((1-\alpha)\cdot
nRT)^{2}}=\frac{\alpha^{3}}{2\cdot(1-\alpha)^{2}}\cdot\bigg{(}nRT\bigg{)}=\frac{\alpha^{3}}{2\cdot(1-\alpha)^{2}}\cdot\bigg{(}\frac{P}{1+\frac{1}{2}\alpha}\bigg{)}=\frac{\alpha^{3}\cdot
P}{(2+\alpha)\cdot(1-\alpha)^{2}}$ (22)

where $n$ is the number density of total H2O _originally_ present, $P$ is the
total pressure of the mixture under _current_ equilibrium condition:
$P=p_{\text{H}_{2}\text{O}}+p_{\text{H}_{2}}+p_{\text{O}_{2}}$, and $\alpha$
is the degree of dissociation defined as:

$\alpha=\frac{p_{\text{H}_{2}\text{O}\text{(decomposed)}}}{p_{\text{H}_{2}\text{O}\text{(total
H${}_{2}$O originally present)}}}$ (23)

For the temperature range from 1000 K to 3000 K, the value of $K_{p}$ is
approximately (which agree with the experimental data within $\sim 1$% of
$\log_{10}K_{p}(T)$) (Karapetyants, 1978):

$\log_{10}K_{p}(T)=\bigg{(}-0.24+1.75\cdot\log_{10}\left[T(\text{K})\right]\bigg{)}-\frac{25,266}{T(\text{K})}$
(24)

Let’s again introduce $T_{3}\equiv T/10^{3}$K to simplify the above
expression:

$\log_{10}K_{p}(T_{3})=\bigg{(}5.01+\underbrace{1.75\cdot\log_{10}T_{3}}_{\text{due
to the entropy difference between products and reactants:
}\int_{0}^{T}\frac{\Delta C_{p}\cdot dT}{R\cdot
T}}\bigg{)}-\frac{25.266}{T_{3}}$ (25)

For 1000K-3000K, $K_{p}$ can be further simplified to:

$\log_{10}K_{p}(T_{3})\approx 5.5-\frac{25.0}{T_{3}}$ (26)

again,

$\log_{10}\left[\frac{\alpha^{3}}{(2+\alpha)(1-\alpha)^{2}}\cdot
P_{\text{atm}}\right]=5.5-\frac{25.0}{T_{3}}$ (27)

For small $\alpha$, we have:

$\log_{10}\left(\frac{\alpha^{3}}{2}\right)+\log_{10}P_{\text{atm}}\approx
5.5-\frac{25.0}{T_{3}}$ (28)

so,

$\log_{10}\left(\alpha^{3}\right)\approx
5.5+\log_{10}{2}-\log_{10}P_{\text{atm}}-\frac{25.0}{T_{3}}$ (29)

or,

$\log_{10}\alpha\approx
1.93-\frac{1}{3}\cdot\log_{10}P_{\text{atm}}-\frac{8.3}{T_{3}}$ (30)

Therefore, on the surface of some water worlds, the strong interaction with
host-stellar photons would lead to gradual dissociation of water molecules and
the escaping of liberated hydrogen over time.

Figure 3: 1%-dissociation front of H2 (blue) and H2O (orange) in the P-T
diagram under equilibrium consideration.

The critical pressure $P_{\text{diss1\%H${}_{2}$O}}$ at which dissociation
reaches one-percent level ($\alpha=1$%) in P-T space (Fig. 3) is:

$\log_{10}P_{\text{diss1\%H${}_{2}$O}}(\text{atm})\approx
11.8-\frac{25.0}{T_{3}}$ (31)

1-bar or 1-atm level is approximately the photosphere of these planets we
consider, where collisionally-induced opacity becomes significant to absorb
most incoming stellar photons.

Notice that although the H2O dissociation front penetrates deeper to reach
higher pressure-levels under the same temperature condition, it is more
difficult to _completely_ dissociate a neutral hydrogen atom out of a water
molecule, because the dissociation of water generally proceeds in two steps:
first one would have to liberate the hydrogen from the strong chemical
affinity of oxygen, then one would have to break the H-H covalent bond.

###  2.7 Phase Transition: Thermal (entropy-limited) versus Shock-driven
(energy-limited) Escape

Figure 4: Schematic diagrams illustrate the escape process of neutral H-atoms
due to the interplay of critical layerings in the planetary atmosphere, in
particular, the sonic line and the dissociation line of X-H bond, and the
emergence of a shock front and critical phenomenon. In (a), the lenticular
body (lens-shaped region colored yellow) is presumably the source region of
hot neutral hydrogen atoms prepared ready for thermal escape as powered by the
incident host-stellar photons. Subplot (a) is applicable to the low-flux
_entropy-limited_ escape scenario. Subplot (b)(c)(d) show the transition into
the _energy-limited_ regime when the flow structure undergoes a _phase
transition_ and is accelerated by an expansive shock front, at either near the
_tropopause_ -level ($\sim 0.1$-bar) for small exoplanets or near the
_photosphere_ -level ($\gtrsim 1$-bar) for giant exoplanets. Green arrows show
the escaping H-atoms (see Appendix H). Subplot (b) shows that the _sonic line_
bends inward under higher photon flux to catch up with the _dissociation
front_ , when they touch they form a shock front. Subplot (c) shows that when
this shock front touches the _tropopause_ (dashed line), it then has gained
direct access to the bulk circulating material of planetary body. Subplot (d)
shows that the escape process is now coupled to the global planet circulation
and is no longer limited by any diffusive process. Matter flowing into this
_shock front_ is dissociated and accelerated across this thin shock front to
high velocity of $\sim$10-30 km/s. The outflow likely becomes beamed also.
Effects due to planet rotation and orbital revolution have been ignored in
this simplified picture. If included, however, observations show that the hot
spot or the source region of escaping flux can either lead ahead or lag behind
the direct sub-stellar point as a result of planetary-scale rotation and fluid
motion.

Now, let’s visualize a planetary surface rich in H2 or H2O subject to host-
stellar photons’ bombardment. When we tune up the photon flux, both the
_dissociation line_ and the _sonic line_ would push deeper into the planet
fluid body, but at somewhat different paces. Under low-flux, the _sonic line_
always lies higher above. However, under higher flux, _sonic line_ would
quickly catch up with the the _dissociation line_. When they touch, it would
lead to the formation of a _shock front_. Thus, this leads to the emergence of
a critical phenomenon at a critical photon flux level (Fig. 4). Essentially,
this creates a planetary-scale rocket engine (Tsiolkovsky et al., 2004).

Under low-flux (Fig. 4(a)), the _lens-shaped region_ is like a combustion
chamber which sources the hot expanding neutral atomic hydrogen gas, and is
powered by incoming host-stellar photons. The nozzle of this rocket engine is
like the sonic level where the escaping outflow transitions from being sub-
sonic to super-sonic due to the conservation of mass flux (continuity
equation), and achieves a velocity of the same order-of-magnitude as the
planet escape velocity ($v_{\text{esc}}\equiv\sqrt{2GM_{p}/r_{c}}$) or even
higher. Additional heating of the escaping flux along its course of outflow
may occur due to other physical mechanisms which would increase its degree of
ionization.

Under high-flux (Fig. 4(b,c,d)), the _lens-shaped region_ in (a) would
transform into a shock front.

This escaping flux of material provides an alternative pathway to cool the
planet surface in addition to the escaping flux of thermal infrared photons
from the planetary surface/photo-sphere. A detailed energy balance can be
shown that the escape flux of material is highly sensitive to the planet
equilibrium temperature, which can be calculated from the bolometric flux-
level of its host star upon the planet (see Appendix H).

The stellar photons which the planet intercepts were in thermal equilibrium
with the stellar surface at a much higher temperature (a few thousand Kelvin),
while they are far far away from equilibrium with the planetary surface, in
terms of their spatial distribution as well as energy distribution when they
reach the planetary surface.

Therefore, it is exactly this dis-equilibrium which drives flows and all kinds
of phenomena occurring on the planetary surface. Thus, the entire planet
surface with its fluid can be viewed as a heat engine which absorbs the
stellar photons from a high temperature heat source and do all kinds of useful
work such as moving fluids around and dissociating molecules, while eventually
cooling off to space with infrared photons as the low temperature heat sink.
One example of such work is to _liberate_ certain atoms or certain molecules
from the bondage of chemical bond as well as the gravitational bond into
space. In particular, we consider the liberation of single neutral hydrogen
atom (H-atom), the lightest atom in our universe, from the bondage of hydrogen
molecules (H-H) and other hydrides (X-H), and from the planetary surface into
space.

Because of that, the escape process is fundamentally driven by an _asymmetric
driving force_ , and we expect a regime-shift to occur as flux level is
increased above a _critical point_ , which is analogous to a _phase
transition_ potentially describable by a theory of _symmetry-breaking_. The
geometry of layering under static atmosphere would morph locally in the hot
spot region, to change its curvature and even merge to produce discontinuity
under more dynamic setting with higher incoming stellar photon flux.

One can visualize that as the incoming photon flux level increases, the _sonic
line_ would bend inward and retreat backwards towards the planetary surface
and eventually merge with the _dissociation front_. When these two fronts
merge, then a _shock front_ forms. Matter which flows into this shock front
are quickly (partially)-dissociated within this very thin shock front and then
accelerated to high escaping velocity ($\sim$10-30 km/s) which is much higher
than typical sonic speed ($\sim$1-3 km/s) expected on a planetary surface (see
Appendix H.1).

Then, the depth that this _shock front_ eats into the planet body depends on
the incoming photon flux level. When the _shock front_ reaches the
_tropopause_ , then the escaping neutral H-atom flux is no longer limited by
diffusion but is directly coupled to the winds and circulation within the
troposphere of the planetary body, and greatly enhanced. This transition may
be paralleled with a _phase transition_ from an _entropy-limited_ regime to an
_energy-limited_ regime.

At even higher incoming photon flux level, this _shock front_ would touch the
_photosphere_ directly so that photon energy absorbed there is directed
converted into the dissociation energy as well as kinetic energy of escaping
atoms.

Ordinarily, such an _”expansive” shock front_ is impossible in nature, meaning
that the inflowing material into a shock front is _subsonic_ while the
outflowing material away from a shock front is _supersonic_ , because it
violates the _Second Law of Thermodynamics_ that there is a net entropy
decrease across the shock front if there are no other physical mechanisms
involved (Landau & Lifshitz, 2013). Precisely, here the dissociation of
molecules and the breakage of H-H or X-H covalent bond provides this physical
mechanism that it is the extra _degree-of-freedom_ in this system which
ultimately generates more entropy by the absorption and conversion of stellar
photons into liberating neutral H-atom.

Yet another alternative way to visualize the physics is that as the incoming
photon flux is increased, the planet’s atmosphere is trying to catch up with
higher and higher wind speed. However, there is a physical limit (sonic speed)
in the material property of any fluid which can mobilize and carry energy
around so much fast without forming a shock front. So, at a _critical point_
of incoming photon flux level, the matter flow can no longer catch up with its
limited sonic speed and a _shock front_ emerges in the hot spot region. This
_shock front_ is again interpreted as a _phase transition_ which accelerates
the exiting outflow from the ordinary thermal escape scenario into strong
hydrodynamic escape scenario. Within this shock front, which is very thin
compared to the physical dimension of the planet, the hydrogen-molecules and
other hydrogen bearing molecules such as water-molecules dissociate and
produce hydrogen atoms to escape at high velocity of $\sim$10-30 km/s.

Now this (1) _phase transition_ in the escape geometry directly correlates
with the (2) _phase transition_ of exoplanet temperature-dependent
distribution in the mass-radius diagram, which is ultimately tied to the (3)
_phase transition_ from hydrogen molecule H-H and hydride molecules X-H into
free hydrogen atom. Marvelously, the three _phase transitions_ tie together in
one physical problem, a trinity! Also, this shock-front driven escape
mechanism involves three fronts: transonic front, collisional front (von
Karman line), and dissociation front merge into one single shock front
(another trinity!) to drive a beamed outflow. Thus, there is double trinity in
this physical problem.

This _phase transition_ may help explain the peculiar _negative slopes_ of the
two observed cut-offs in the mass-radius diagram to be discussed in the next
section, as compared to positive slopes expected in the ordinary thermal
escape scenario with contours more or less parallel to the gravitational
potential contour M/R.

The winds at the level of _troposphere_ and _photosphere_ would be first
driven by pressure gradient which is induced by the thermal gradient and
converge towards the hot spot from all directions (Fig. 4(c,d)). Along the
way, the matter in the winds gradually pick up more and more thermal energy
and also speed by continuously absorbing host-stellar photons. Of course, the
winds also radiate back into space but there is a net increase in their energy
as they _sail_ towards the hot spot. Eventually, when all these winds converge
and collide in the hot spot they dump their thermal energy (characterized by
small-scale random motion) as well as kinetic energy (characterized by large-
scale uniform motion) entering the _shock front_. Then, this converging energy
dissociates molecules and drive escape. Furthermore, the compensating
(returning) material flux is radiatively-cooled in a compensating flow away
from hot spot region at a somewhat higher level in the atmosphere (Fig. 4(d)).
This is a planetary-scale heat engine which uses photosphere and its working
material—fluid on planet surface as conveyor-belt to concentrate stellar
energy absorbed into a small region, liberate neutral H-atoms, and produce
collimated outflow.

Further investigations of the escape flux and its dependence on (1) the
planetary system architecture (2) and the host stellar spectral type are
discussed in Appendix I and Appendix J.

##  3 Mass-Radius Diagram

An updated version of the _Mathematica_ tool for synthesizing and analyzinng
the multi-dimensional exoplanets’ data (Zeng et al.,
2021),(https://community.wolfram.com/groups/-/m/t/2445247) is utilized to
study the well-studied transiting exoplanets in the Transiting Exoplanet
Property Catalogue (TEPCat) (Southworth, 2011),
(www.astro.keele.ac.uk/jkt/tepcat/), to interpret the Exoplanet Systematics in
the Mass-Radius Diagram.

Here, we consider two sets of mass-radius diagrams: one set for smaller
exoplanets which include the super-Earths/Mini-Neptunes with typical mass
range of 0-20 Earth masses (M⊕) and radius range of 0-5 Earth radii (R⊕), and
the other set for gas giant exoplanets with typical mass range of 50-2000
Earth masses (M⊕) and radius range of 8-20 Earth radii (R⊕). Since we consider
hydrogen as one of the main components of the atmosphere and envelope of some
exoplanets, we need to exploit the hydrogen Equation-of-States (EOS) in
detail, which is done in Appendix A.

The mass-radius distribution of exoplanets is highly sensitive to the planet
equilibrium temperature, or equivalently, the flux-level that planet receives
from its host star. The actual temperature distribution on planet surface may
not be uniform due to geometric factors as well as the efficiency of energy
transport on planet surface, and planetary rotation. The details are discussed
in Appendix B. The average planet equilibrium temperature (T${}_{\text{eq}}$)
have a simple correlation with the temperature of sub-stellar hot spot—the
highest temperature achievable on such a planet surface, and also the source
region of the material escape flux considered in our model.

In the process of analyzing the mass-radius diagram, we identify two cut-offs
as straight lines with negative slope in the lgM-lgR plot, one cut-off for the
small exoplanets which we call the Cosmic Hydrogen and Ice Loss Line #1 (CHILL
#1), and the other cut-off for the giant exoplanets which we call the Cosmic
Hydrogen and Ice Loss Line #2 (CHILL #2).

###  3.1 CHILL#1 for small exoplanets

First, let us discuss the CHILL#1 for small exoplanets, as shown in Figure 5
and Figure 6.

These small exoplanets possess a certain degree of their primordial hydrogen
envelope or thick atmosphere, which place them at their corresponding
positions on the mass-radius diagram. Their envelope or thick atmosphere is
expected to undergo escape process as discusses earlier. They are also
expected to possess a core made of either rock-metals or hydrogen-bearing
ices.

Due to the selective escape of lighter species, their planetary atmosphere is
expected to be gradually enriched in heavier species over the course of time,
which increases the mean molecular weight of the atmosphere and decreases the
scale height. Thus, eventually, their planetary radii will become close to
that of bare cores, either rock-metal or water-ice, as calculated from theory.

Therefore, we interpret the observed bi-modal distribution of small exoplanet
radii and the small exoplanet radius gap at around 2 R⊕ as caused by the
density difference between a rock-metal core and a water-ice-rich core (Zeng
et al., 2017, 2018, 2019, 2021).

Other possibilities of explaining the small exoplanet radius gap can be found
in the following literature (Berger et al., 2018, 2020; Fulton et al., 2017;
Fulton & Petigura, 2018; Petigura et al., 2018, 2022; Owen & Wu, 2013, 2017;
Wu, 2019; Van Eylen et al., 2018; Ginzburg et al., 2018; Rogers et al., 2021,
2023; Izidoro et al., 2022; Lehmer & Catling, 2017).

Furthermore, we consider the CHILL#1 as a physical boundary differentiating
two group of small exoplanets: one group which still partially retains their
primordial hydrogen envelope of near _nebular_ composition, and the other
group which has completely lost their primordial hydrogen envelope due to
higher planet equilibrium temperature, but may possess a highly-evolved
secondary envelope/atmosphere of near _meteoritic_ /_planetary_ composition.

Thus, there is a significant temperature contrast for small exoplanets on
either side of this physical boundary. The critical temperature of this
physical boundary CHILL#1 is estimated at T${}_{\text{eq}}\approx$1000K, or
equivalently, a flux-level of two hundred times the Earth’s flux-level
($f/f_{\oplus}\approx 200$).

This flux cut-off is further confirmed and elucidated in additional plots in
Appendix I and Appendix J.

This physical boundary CHILL#1 in the mass-radius space can be approximated by
the following analytical relation:

$(M_{p}/10M_{\oplus})^{2}\times(R_{p}/2R_{\oplus})^{3}\approx 1$ (32)

where the product of the planet mass squared and the planet radius cubed
equals a constant approximately along this physical boundary.

Another equivalent expression is:

$(M_{p}/10M_{\oplus})^{2/3}\times(R_{p}/2R_{\oplus})\approx 1$ (33)

where the product of the planet mass to the two-thirds power and the planet
radius equals a constant along this physical boundary.

This physical boundary CHILL#1 is represented as the thick dashed curve in
Figure 5 and Figure 6. The most interesting feature of it is its negative
slope, as opposed to the positive slope of a typical mass-radius curve of a
uniform composition in the regime that we consider. This negative slope is
most likely due to the intrinsic physics and chemistry of the escape process.

Figure 5: CHILL#1 (thick dashed line) for Small Exoplanets: Super-Earths/Mini-
Neptunes/Cores with H2-Envelopes. Linear-Linear Plot in Planet Masses
($M_{p}$) and Radii ($R_{p}$). Planets are color-coded logarithmically
according to their average global equilibrium temperatures
(T${}_{\text{eqavg}}$) assuming bond-albedo = 0 and efficient heat transport
and re-distribution. The purple mass-radius curve for H2 is for cold dense
molecular fluid which is approximately one-sixteenth the density of cold H2O
fluid at similar pressures or planet masses. The four thin dashed mass-radius
curves in between that of cold H2 fluid and cold H2O fluid are for the 10%,
20%, 30%, and 50% mass-mixtures of the two end members. The subplot at the
top-right corner shows the mass-radius relation of the host stars color-coded
according to the stellar luminosity. These host-stars are Main-Sequence (MS)
stars hotter, bigger, and more massive than M-dwarfs.  Figure 6: CHILL#1
(thick dashed line) for Small Exoplanets: Super-Earths/Mini-Neptunes/Cores
with H2-Envelopes. Linear-Linear Plot in Planet Masses ($M_{p}$) and Radii
($R_{p}$), Zoom-in. Planets are color-coded logarithmically according to their
incident bolometric host-stellar flux ($f/f_{\oplus}$) where $f_{\oplus}=1.36$
kW/m2 is the standard flux received at Earth-distance from the Sun (solar
constant). According to the planet radius ($R_{p}$) histogram attached to the
right of the main mass-radius diagram, it seems there is a possibility to
further sub-divide the planets which are in between 2-3 $R_{\oplus}$ into two
sub-populations. Furthermore, the clustering of exoplanets near the
intersection of CHILL#1 and pure-H2O mass-radius curve is interesting.

###  3.2 CHILL#2 for giant exoplanets

Then, let us consider the CHILL#2 of gas giant exoplanets, as shown in Figure
7 and Figure 8.

These gas giant exoplanets are dominated by hydrogen and helium in their
composition. Therefore, the CHILL#2 represents very significant mass loss or
even complete destruction of these exoplanets. A simultaneous physical factor
to consider is the significant inflation of the planet radius due to its host-
stellar insolation, which is also the physical mechanism which drives the
escape process. Thus, a strong correlation is expected between the measured
planet (inflated) radius, planet mass, and its flux-level.

Another important physical factor which affects the measured planet radius is
the metallicity of the planet hydrogen envelope. Due to internal convective
motion and mixing, a significant portion of metals, elements heavier than
hydrogen and helium, is likely dispersed throughout the envelope, with a
possible gradient concentrating towards the planet core. It has been
demonstrated by calculation that what matters is the total metallicity or
total mass of metals contained within the planetary body, regardless of its
distribution in the envelope or core (Stevenson, 2017). Thus, to the leading
order, we can express the relationship between planet radius, planet mass,
flux-level, and planet metallicity in the following functional form:

$R_{p}=\mathcal{F}(M_{p},\text{flux},\text{metallicity})$ (34)

where the terms on the right-hand side inside the function $\mathcal{F}$ are:
(1) the planet mass (Mp); (2) the incident bolometric host-stellar flux; (3)
the metallicity.

Here we choose the host-stellar [Fe/H] to represent the planet metallicity
because it is most readily measurable. However, we need to keep in mind that
there may be significant difference between the host-stellar metallicity and
the planet metallicity as caused by the ongoing escape process. In other
words, the planets are expected to be gradually enriched in metals due to the
escape process.

Also, an implicit assumption of the correlation between the iron abundance
with the abundance of other important planet-building elements such as
characterized by [C/O], [Mg/Si], [Fe/Si], [O/H], etc.) are made.

Nonetheless, adding in the host-stellar [Fe/H] as an extra dimension in the
mass-radius plot gives us insight into the physics which shapes the planet
distribution, as demonstrated in Figure 9.

Detailed examinations of Figure 7, Figure 8, and Figure 9 show that
significant inflations of giant exoplanet radii begin to occur in the planet
equilibrium temperature range of 1000-1500K (T${}_{\text{eq}}$).

Moreover, complete destruction of giant exoplanets begins to occur on a grand
scale in the planet equilibrium temperature range of 1500K-2000K
(T${}_{\text{eq}}$). In Figure 8, we find a sharp truncation or physical
boundary in the mass-radius diagram which we identify as the CHILL#2.

This physical boundary or the CHILL#2 has a well-defined negative slope of
-2/3 in the log-log plot of mass-radius space, again. We can also express this
relation analytically as:

$(M_{p}/300M_{\oplus})^{2}\times(R_{p}/10R_{\oplus})^{3}\approx 1$ (35)

or equivalently, in unit of Jupiter mass (M${}_{\text{J}}$) and Jupiter radius
(R${}_{\text{J}}$),

$(M_{p}/M_{\text{J}})^{2}\times(R_{p}/0.9R_{\text{J}})^{3}\approx 1$ (36)

In particular, Figure 9 and additional plots in Appendix K show that the
CHILL#2 is located at the planet equilibrium temperature of 1500K
(T${}_{\text{eq}}$) assuming zero bond albedo, or equivalently, eight hundred
times the Earth’s flux-level ($f/f_{\oplus}\approx 800$). This suggests that
the giant exoplanets which lie along the edge of this physical boundary in
Figure 8 receives about eight hundred times the Earth’s flux. When converted
into standard units, this quantity is $10^{6}$W/m2, or one million Watts per
square meter on a surface perpendicular to the host-stellar irradiation summed
over the entire electro-magnetic spectrum. Thus, we believe this strong
insolation serves as the major power source which drives the planet envelope
and atmosphere escape.

Figure 7: CHILL#2 for Large Exoplanets: Gas Giants/Hot Jupiters. 1000K-1500K
(200$f_{\oplus}$-800$f_{\oplus}$). Significant inflation of planet radii occur
within this equilibrium temperature range. Notice the approximate triangular
shape of planet distribution region in lgM-lgR, which tapers off towards
higher masses. This overall distribution is a consequence of three factors
combined: (1) planet mass (Mp); (2) insolation received ($f/f_{\oplus}$) or
equilibrium temperature (T${}_{\text{eq}}$); (3) metallicity [Fe/H].  Figure
8: CHILL#2 for Large Exoplanets: Gas Giants/Hot Jupiters. 1500K-2000K
(800$f_{\oplus}$-3000$f_{\oplus}$). Significant evaporative loss and
destruction of gas giants occur within this equilibrium temperature range. The
cut-off with (-2/3) slope in lgM-lgR is evident. Notice the change of the
temperature scale.  Figure 9: CHILL#2 for Large Exoplanets: Gas Giants/Hot
Jupiters. 1500K-2000K. This is a re-rendering of Figure 8, with a stricter
selection in planet errorbars, and a new subplot showing the insolation flux
received ($f/f_{\oplus}$) versus planet radii ($R_{p}/R_{\oplus}$). The flux
is measured in the unit of what Earth receives from the Sun, i.e., the solar
constant ($f_{\oplus}=1.36$ kW/m2). The cut-off in the lgM-lgR diagram shows
up as a cut-off in flux of ($f_{\text{crit}}/f_{\oplus}\approx 800$). This is
expected because the lower limit of the T${}_{\text{eq}}$-range is set at
1500K which exactly corresponds to this $f_{\text{crit}}$. The gas giants
which can survive at or above this threshold flux level generally are more
massive than Saturn ($\gtrsim$ 100 M⊕). Furthermore, the planets in the main
plot as well as the subplot are now color-coded according to their host-
stellar metallicity [Fe/H]. Their inflated radii correlate with their
metallicity [Fe/H] with the general trend of lower-metallicity ones having
larger radii (more-inflated) and higher-metallicity ones having smaller radii
(less-inflated), as expected.

###  3.3 Discussion

The key feature that we want to solve and explain by theory is the negative
slope of the CHILL#1 and the CHILL#2.

Contrary to intuition, the negative slope suggests along the physical boundary
of escape, the less massive and more puffier exoplanets which have a shallower
gravitational potential can somehow last for the same timescale as compared to
the more massive and less puffier exoplanets which have a deeper gravitational
potential. Moreover, the puffier exoplanets have proportionally larger surface
area to intercept the incoming host-stellar radiation which should result in
more escape. Therefore, we need to consider and calculate the efficiency of
the aforementioned planet-scale rocket engine under different parameters, that
the rocket engine on puffier exoplanet may be of poorer efficiency as compared
to the rocket engine on a less puffier and more compact exoplanet, due to
self-shielding effect and metallicity effect.

We speculate there is a _shock front_ emerging in the critical scenario which
accelerates the flow to super-sonic speed (please see discussion in Appendix
H.1). Thus, there exists a _phase transition_ in the escape scenario from
ordinary thermal escape to _shock-wave_ accelerated scenario. This might
explain the change of the positive-slope boundary delineating the boundary of
hot Saturns to the negative-slope boundary cutting off the hot Jupiter
population as well as the small exoplanet population in the mass-radius
diagram.

##  4 Conclusion

Two Cosmic Hydrogen and Ice Loss Lines (CHILLs #1 & #2) have been identified
for small exoplanets and gas giant exoplanets in the mass-radius diagram. The
CHILL#1 corresponds to a critical boundary of planetary equilibrium
temperature (T${}_{\text{eq}}$) of 1000 K. The CHILL#2 corresponds to a
critical boundary of planetary equilibrium temperature (T${}_{\text{eq}}$) of
1500 K.

This negative slope requires further investigation, and calls for the
attention of the study of the _planetary-scale rocket engine_. In conclusion,
we hypothesize that many planets ($\gtrsim 1$M⊕) have initially acquired
various degrees of thick primordial H2-He envelopes from being embedded in
gas-rich disks. This primordial H2-He envelope on any close-in planet is then
gradually lost over time which is powered by the host-stellar irradiation. A
viable escape mechanism is the _shock-front-driven_ hydrodynamic outflow of
neutral hydrogen atoms which are dissociated from hydrogen molecules (H-H) and
other hydrogen-bearing molecules (X-H) such as the gaseous form of cosmic ices
including water, ammonia, and methane, due to the absorption of host-stellar
photons and subsequent heating, molecular dissociation, and acceleration.

##  5 Acknowledgement

The author Li Zeng would like to thank Robin Wordsworth, Junjie Dong, Collin
Cherubim, Kaitlyn Loftus, Michail I. Petaev, Joe Gonzeles, and James G.
Anderson of Harvard University Department of Earth and Planetary Sciences for
fruitful discussions.

##  Appendix A The Hydrogen EOS

The H2 EOS is needed to understand the hydrogen-dominated envelopes and
atmospheres on planets.

Figure 10: Schematic Representation of Hydrogen EOS in three regimes: (1)
Gaseous, (2) Dense (super-critical) Molecular Fluid, and (3) Metallic Fluid.
The hydrogen EOS is commonly encountered in exoplanet interior modelling,
because hydrogen is the most abundant element available in the universe. The
(1) gaseous EOS can be first approximated by an ideal gas law, and then with a
modification term taking into account inter-molecular collisions and molecular
dimensions. The (3) metallic EOS can be first approximately by a power law
with power-index of two which takes into account the energy of de-localized
electrons. The (2) molecular fluid EOS bears most of the complexity (grey-
shaded region) because molecules strongly interact with neighbors that their
average distance is shorter than the sum of their van der Waal radii and the
pressure effect and temperature effect are of equal importance.

The EOS is best-described by considering three different regimes from the
structure of matter perspective (see Figure 10): (1) gaseous, (2) high-density
molecular (super-critical) fluid, (3) high-density metallic fluid.

(1) Gaseous EOS. This is applicable when the distance of neighboring molecules
is large compared to the linear dimension of molecule ($\sim 10^{-8}$cm). The
ideal gas law is the simplest EOS to describe such a state ($\rho\ll
0.07$g/cc), while modification of which gives better approximation to the EOS
at relatively high density taking into account the interactions of neighboring
molecules. One convinient form of such a real gas EOS is given by (Tsien,
1965), by defining a characteristic temperature
$\Theta_{1}\equiv\epsilon/k_{B}$, where $\epsilon$ is the characteristic
energy of inter-molecular potential, $v\equiv\mu/\rho$ is the average total
volume assigned to each molecule, and $v^{*}\equiv D^{3}$ is the volume
occupied by each molecule’s physical dimension (effectively the reach of their
outermost occupied electronic orbitals), that is, the cubic power of the
equivalent molecular diameter $D$,

$\underbrace{\frac{P\cdot v}{k_{B}T}=1}_{\text{ideal gas
EOS}}+\underbrace{\left(\frac{1}{0.278\cdot\left(T/\Theta_{1}\right)^{1/6}\cdot\left(v/v^{*}\right)-0.177}\right)}_{\text{non-
ideal correction: up to $\sim$few to ten, but reduces to 0 at high T}}$ (A1)

Notice that the first part of such an EOS is simply the ideal gas EOS, while
the second part gives the non-ideal correction primarily due to occupied
volume of molecule. This correction becomes smaller at higher temperatures
because the effective size of a molecule in collision becomes smaller due to
higher kinetic energy of the molecules involved in a collision. The effective
collisional cross-section $\sigma$ for molecules decreases with increasing
temperature or increasing collisional kinetic energy, because molecules are
not rigid spheres but are _squeezy_. In other words, the effective size of
each molecule which makes its way into the correction term decreases.

(2) High-density molecular (super-critical) fluid hydrogen EOS. This EOS is
generally complicated because molecules strongly interact with each other with
partially-overlapping outer electronic orbitals with neighbors in this
physical regime. The general feature is that this EOS quickly stiffens when a
critical density is reached. This critcal density for hydrogen fluid is about
1/16-th that of water at ambient conditions. Then, this fluid EOS follows the
general compression pattern of other materials, with similar slopes and curve
shapes on the mass-radius diagram. The EOS eventually softens when the matter
metallizes by further compression for a factor of about ten-fold in volume
(from $\sim 0.07$g/cc compresses to $\sim 0.7$g/cc), or equivalently, a
reduction of a factor of two in each linear dimension. At this point, the
electronic orbitals of neighboring molecules overlap so much that electrons
become de-localized.

(3) High-density metallic fluid hydrogen EOS. Ultimately the compression leads
to the metallization where the EOS again softens and becomes simple due to the
presence of de-localized electrons. For a significant range of density and
pressure this metallic fluid EOS can be approximated by a power-law with
power-index of two (Stevenson, 2017),

$\displaystyle\underbrace{P\approx K\cdot\rho^{2}\approx
K^{\prime}\cdot\frac{1}{v^{2}}}_{\text{characteristic of any metal under
extreme compression}},~{}$ where $K$ and $K^{\prime}$ are proportionality
constants. (A2) $\displaystyle K\approx 200\text{~{}GPa}\approx 2\cdot
10^{11}\text{~{}Pa},~{}$ where $\rho$ is measured in g/cc (A3)

This value of proportionality constant $K$ can change if heavier elements or
components with higher mean molecular weights are mixed into this envelope.

Assume _volume additivity_ , and define the mass fraction of total heavier
elements other than hydrogen-helium in the mixture as $Z$ (Stevenson, 2017),
we have,

$\frac{1}{\rho(P)}=\left(\frac{1-Z}{\sqrt{P/K}}+\frac{Z}{\rho_{\text{heavy}}}\right)\approx\frac{1-Z}{\sqrt{P/K}}$
(A4)

This simply suggests that we can define an effective $K_{\text{eff}}$ so that,

$P=K_{\text{eff}}\cdot\rho^{2},\text{~{}where $K_{\text{eff}}=(1-Z)^{2}\cdot
K$}$ (A5)

The radius of a planet in this regime does not vary much with the change of
mass and can be approximated as:

$R_{p}\approx\sqrt{\frac{\pi K_{\text{eff}}}{2G}}\approx\sqrt{\frac{\pi
K}{2G}}\cdot(1-Z)$ (A6)

where $G$ is the gravitational constant.

Thus, the planet radius or envelope thickness is reduced by the factor
$(1-Z)$. (Stevenson, 2017) also shows that $R_{p}$ remains essentially
_unchanged_ if this much $Z$ heavier elements all concentrate into a core. So
$R_{p}$ is only related to the total amount ($Z$) but not sensitive to the
exact distribution of heavier elements within the planet interior profile. We
will then use this relation in a later discussion of the negative slope of
CHILL in the mass-radius plot.

##  Appendix B On Planet Temperature

This average equilibrium temperature ($T_{\text{eq}}$) for a planet on an
ecliptic orbit can be calculated as follows.

First, let’s consider the total amount of stellar insolation the planet
receives averaged over one full orbital period:

$\bar{f}=\frac{1}{P}\cdot\int_{0}^{P}\frac{L_{\star}}{4\pi r^{2}}\,dt$ (B1)

Then, we can substitute the $r^{2}$ in the denominator with the angular
velocity $\dot{\theta}$ using the Kepler’s Second Law.

$\text{const}=\frac{1}{2}r^{2}\dot{\theta}=\frac{1}{2}r^{2}\frac{d\theta}{dt}=\frac{\pi
ab}{P}=\frac{\pi a^{2}\sqrt{1-e^{2}}}{P}$ (B2)

Then, recall the Kepler’s Third Law we can re-write $P$ as:

$P=\frac{2\pi}{\sqrt{G(M_{\star}+M_{p})}}\cdot
a^{3/2}\approx\frac{2\pi}{\sqrt{GM_{\star}}}\cdot a^{3/2}$ (B3)

Then, we can figure out the ”const” as:

$\text{const}=\frac{\pi
a^{2}\sqrt{1-e^{2}}}{P}\approx\frac{1}{2}\sqrt{GM_{\star}\cdot a(1-e^{2})}$
(B4)

Thus, the inverse of orbital distance ($r$) squared is directly proportional
to the angular velocity ($\dot{\theta}$):

$\frac{1}{r^{2}}=\frac{\dot{\theta}}{\sqrt{GM_{\star}\cdot a(1-e^{2})}}$ (B5)

Then, the integral can be easily carried out because a planet simply revolves
about $2\pi$-angle in each orbital period:

$\bar{f}=\frac{1}{P}\cdot\int_{0}^{P}\frac{L_{\star}}{4\pi}\cdot\frac{1}{r^{2}}\,dt=\frac{1}{P}\cdot\int_{0}^{P}\frac{L_{\star}}{4\pi}\cdot\frac{\dot{\theta}\,dt}{\sqrt{GM_{\star}\cdot
a(1-e^{2})}}=\frac{1}{P}\cdot\frac{L_{\star}}{4\pi}\cdot\frac{2\pi}{\sqrt{GM_{\star}\cdot
a(1-e^{2})}}=\frac{L_{\star}}{4\pi a^{2}\sqrt{1-e^{2}}}$ (B6)

Notice the extra factor of $1/\sqrt{1-e^{2}}\approx
1+\frac{1}{2}e^{2}+\frac{3}{8}e^{4}+...$, is the fundamental argument behind
_Milankovitch cycle_. The orbit eccentricity change affects the average
insolation. In the long-run, it builds up and affects the overall climate of
the planet. On the other hand, the change of planet _obliquity_ only affects
the re-distribution of stellar insolation over different areas of the planet
surface without changing the total amount of energy received, assuming albedo
does not change.

Thus, then the equilibrium temperature of the planet can be calculated as
follows:

First, the _average_ equilibrium temperature $T_{\text{eq}}$ of the planet
surface can be considered as the temperature if the total stellar insolation
received over $\pi R_{p}^{2}$ is transported and uniformly distributed over
its entire surface area of $4\pi R_{p}^{2}$:

$T_{\text{eq}}=\left(\frac{\bar{f}}{4\sigma}\right)^{1/4}=\left(\frac{L_{\star}}{4\sigma\cdot
4\pi a^{2}\sqrt{1-e^{2}}}\right)^{1/4}=\left(\frac{4\pi R_{\star}^{2}\sigma
T_{\star}^{4}}{4\sigma\cdot 4\pi a^{2}\sqrt{1-e^{2}}}\right)^{1/4}$ (B7)

where $\sigma$ is the _Stefan-Boltzmann constant_ , $a$ is orbit semi-major
axis, and the effect of orbit eccentricity $e$ is included, assuming bond-
albedo = 0 (meaning complete absorption of incident host-stellar bolometric
flux) and efficient heat transport and re-distribution over the planet
surface. Therefore,

$T_{\text{eq}}=T_{\text{star}}/(2a/R_{\text{star}})^{1/2}/(1-e^{2})^{1/8}=\frac{T_{\star}}{\sqrt{2a/R_{\star}}}\cdot\left(1+\frac{1}{8}e^{2}+\frac{9}{128}e^{4}+...\right)$
(B8)

In this paper, when we refer to the equilibrium temperature we always refer to
$T_{\text{eq}}$, i.e., the temperature averaged over the entire planet sphere
and over one full orbital period, assuming zero-albedo. Locally, certain areas
on planet surface can become hotter or colder due to inefficient or incomplete
or different mode of heat re-distribution.

Now, if this close-in planet is tidally locked to its host star, and if its
eternal day-side never communicates heat with its eternal night-side but
averages within its hemisphere, and then re-radiate its received energy as
infrared photons into space, then, we can estimate its day-side average
hemispherical temperature as:

$T_{\text{dayside-hemisphere}}=2^{1/4}\cdot T_{\text{eq}}$ (B9)

Furthermore, if the sub-stellar point of this planet, which is the point of
the planet under direct vertical insolation from its host star, never
communicates heat with other areas, and then re-radiate its received energy as
infrared photons back into space, then, we can estimate its sub-stellar point
temperature as (the highest possible temperature achievable on such a planet):

$T_{\text{sub-stellar}}=2^{1/2}\cdot T_{\text{eq}}$ (B10)

Likewise, in such a scenario of no heat communication, for a spot which
subtends angle $\theta$ from the sub-stellar point on the day-side its
temperature is estimated to be:

$T_{\text{dayside}}\left(\theta\right)=T_{\text{sub-
stellar}}\cdot\left(\cos{\theta}\right)^{1/4}=2^{1/2}\cdot
T_{\text{eq}}\cdot\left(\cos{\theta}\right)^{1/4}$ (B11)

Figure 11: This $\left(\cos{\theta}\right)^{1/4}$-dependence. Rapid change
occurs near the rim $\frac{5}{12}\pi$-$\frac{1}{2}\pi$.

Its night-side may have some arbitrarily low temperature.

Thus, for more detailed quantitative calculation, one needs to consider which
hemisphere or portion of the planet is heated and is responsible for the main
escape flux. This also depends on the geometry of heat re-distribution on such
a planet.

In terms of orbital parameters alone, if not including planet albedo effect,
the equilibrium temperature of a planet primarily depends on its orbit semi-
major axis $a$, and secondarily on the effect of its orbit eccentricity $e$,

$T_{\text{eq}}=T_{\text{star}}/(2a/R_{\text{star}})^{1/2}/(1-e^{2})^{1/8}$
(B12)

Observationally speaking, the $T_{\text{eq}}$-characteristic of Cosmic
Hydrogen-Ice Loss Line (CHILL#1) and Hydrogen Loss Line (CHILL#2) are:
T${}_{\text{eq}}\sim$1000K (or approximately 200 times the Earth’s insolation
($f_{\oplus}$)) for Super-Earths/Mini-Neptunes—planets of masses on the order
of a few up to a few tens of Earth masses and T${}_{\text{eq}}\sim$1500K (or
approximately 800 times the Earth’s insolation ($f_{\oplus}$)) for massive gas
giants—planets of masses on the order of a few hundred Earth masses or one
Jupiter mass. This strong T${}_{\text{eq}}$-dependence of the escape process
can be thought of as an On-Off switch which determines whether a planet with
significant volatile envelope can survive over the planetary system age at a
specified host-stellar irradiation level.

The difference in this critical value between the small and large exoplanets
(200 $f_{\oplus}$ versus 800 $f_{\oplus}$) may be explained by different
escaping geometry required for substantial depletion of the hydrogen envelope
for a small planet versus for a massive gas giant. Assuming inefficient heat
circulation (or zero heat re-distribution), the sub-stellar hot spot of a
highly-irradiated exoplanet can get as high as $\sqrt{2}$ of its globally-
averaged equilibrium temperature T${}_{\text{eq}}$,

$\text{T}_{\text{hotspot}}\approx\sqrt{2}\cdot\text{T}_{\text{eq}}\sim\begin{cases}1500\text{~{}K},&\text{$f/f_{\oplus}\approx
200$ for small planet}.\\\ 2000\text{~{}K},&\text{$f/f_{\oplus}\approx 800$
for large planet}.\end{cases}$ (B13)

The depth of penetration of dissociation fronts is likely different for small
planet versus large planet, and thus, results in different loss rates.

##  Appendix C Jeans Escape

We can express the pressure at the _exobase_ (p${}_{\text{exob}}$) as a
function of collisional cross-section ($\sigma$), the mean-molecular weight of
the escaping species ($\mu$), and the surface gravity of the planet at the
level of exobase ($g_{\text{exob}}$):

$\bigg{(}p_{\text{exob}}=n_{\text{exob}}\cdot k_{B}\cdot
T\bigg{)}\sim\frac{k_{B}\cdot T}{\mathcal{H}\cdot\sigma}\sim\frac{\mu\cdot
g_{\text{exob}}}{\sigma}=\bigg{(}\frac{\mu\cdot
GM_{p}/r_{\text{exob}}^{2}}{\sigma}\bigg{)}$ (C1)

Collisional cross-sections $\sigma$ for atoms (e.g. H, O, …) or small
molecules (e.g. H2, OH, O2, etc.) are similar and reflect their characteristic
dimensions of $\sim 10^{-8}$cm and are commonly near $\sim 10^{-16}$cm2. This
is the same cross-section for photon interaction (scattering and absorption)
with the electron cloud in atoms and small molecules (Lewis, 2004, p.79)
(Tsien, 1965, Chap. 1).

Therefore, we can plug in this number of $\sigma\approx 10^{-16}$cm2 and
obtain:

$p_{\text{exob}}\sim\frac{(1\text{g}/N_{A})\cdot
10^{3}\text{cm}/\text{s}^{2}}{10^{-16}\text{cm}^{2}}\cdot\frac{M_{p}/M_{\oplus}}{(r_{\text{exob}}/R_{\oplus})^{2}}\approx
10^{-6}\text{Pa}\cdot\frac{M_{p}/M_{\oplus}}{(r_{\text{exob}}/R_{\oplus})^{2}}\approx
10^{-11}\text{atm}\cdot\frac{M_{p}/M_{\oplus}}{(r_{\text{exob}}/R_{\oplus})^{2}}$
(C2)

Therefore, the pressure (p${}_{\text{exob}}$) which defines the _exobase_ is a
low number ($10^{-11}$ atm), modulated by a factor proportional to surface
gravity.

The escape flux ($\Phi_{\text{Jeans}}$ in cm-2 s-1): number of atoms/molecules
escaping into outer space per unit surface area at the _exobase_ -level per
unit time) can be found by integrating over the Maxwellian velocity
distribution of the escaping species at the _exobase_ (Walker, 1977; Goody,
1995; Catling & Kasting, 2017):

$\Phi_{\text{Jeans}}=\frac{n_{\text{exob}}\cdot u}{2\cdot\pi^{1/2}}\cdot
e^{-v_{\text{esc}}^{2}/u^{2}}\cdot\left(\frac{v_{\text{esc}}^{2}}{u^{2}}+1\right)$
(C3)

Here $n_{\text{exob}}$ is the number density of escaping species evaluated at
the escape level, that is, the exobase-level.

$v_{\text{esc}}\equiv\sqrt{2GM_{p}/r_{\text{exob}}}=\sqrt{2}\cdot\sqrt{GM_{p}/r_{\text{exob}}}$
is the escape velocity and $u\equiv\sqrt{2k_{B}T/\mu}=\sqrt{2}\cdot u_{0}$ is
the most probable velocity of the Maxwellian velocity distribution under
thermal equilibrium111where $u_{0}\equiv\sqrt{k_{B}T/\mu}$ is the isothermal
sound speed, where $\mu$ is the molecular weight for the gas species
considered., both evaluated at the escape level.

It is convenient to define a dimensionless _Jeans parameter_ $\lambda$ as
(Goody, 1995; Catling & Kasting, 2017; Owen, 2019):

$\lambda\equiv\frac{v_{\text{esc}}^{2}}{u^{2}}=\frac{GM_{p}}{r_{\text{exob}}\cdot
u_{0}^{2}}=\frac{GM_{p}/r_{\text{exob}}}{k_{B}T/\mu}=\frac{g_{\text{exob}}\cdot
r_{\text{exob}}}{k_{B}T/\mu}=\left(\frac{r_{\text{exob}}}{\mathcal{H}_{\text{exob}}}\right)\propto\mu$
(C4)

_Jeans parameter_ $\lambda$ measures how difficult it is for a gas species to
escape, the higher the value of $\lambda$, the more difficult it is to escape.
Since $\lambda\propto\mu$, it is harder for heavier species to escape compared
to lighter species, under the same temperature T condition.

For example, because $\mu_{\text{H}_{2}}=2\cdot\mu_{\text{H}}$, the _Jeans
parameter_ $\lambda$ doubles for a H2-molecule when it is compared to that of
an H-atom: $\lambda_{\text{H}_{2}}=2\cdot\lambda_{\text{H}}$. Thus,
theoretically it is much easier for an H-atom to thermally escape than a
H2-molecule.

The last expression of
$\lambda=\left(r_{\text{exob}}/\mathcal{H}_{\text{exob}}\right)=\left(\mathcal{H}_{\text{exob}}/r_{\text{exob}}\right)^{-1}$
suggests that $\lambda$ can be viewed over the planetary scale, as the ratio
of planet radius at the exobase-level over the scale height of the escaping
species evaluated at the exobase-level.

In other words, $\lambda$ is the reciprocal of the scale height $\mathcal{H}$
measured in the unit of planet radius $\left(R_{p}\sim
r_{\text{exob}}\right)$. The smaller the value of $\lambda$, the higher the
ratio of $\left(\mathcal{H}_{\text{exob}}/r_{\text{exob}}\right)$.

The escape flux ($\Phi_{\text{Jeans}}$) in terms of $\lambda$ is:

$\Phi_{\text{Jeans}}=\bigg{(}\frac{n_{\text{exob}}\cdot
u}{2\cdot\pi^{1/2}}\bigg{)}\cdot\bigg{(}\left(\lambda+1\right)\cdot
e^{-\lambda}\bigg{)}$ (C5)

Ignoring all the factors of the order of unity including the geometric factor,
the escape flux $\Phi_{\text{Jeans}}\propto\exp{(-\lambda})$.

$\Phi_{\text{Jeans}}\sim n_{\text{exob}}\cdot u\cdot e^{-\lambda}$ (C6)

The smaller the value of $\lambda$, the more intense the escape. In fact, this
formula is not sensitive to the exact level (as long as it is below the sonic
level and above the von Karman line, in the upper atmosphere) at which both
$n$ and $\lambda$ are evaluated because the product $n\cdot e^{-\lambda}$ is
largely conserved due to the _barometric law_ (Catling & Kasting, 2017). Later
we will show the advantage of evaluating it right at the von Karman line.

On the other hand, the ratio
$\left(\Phi_{\text{Jeans}}/n_{\text{exob}}\right)$ is often referred to as the
effusion velocity (Goody, 1995):

$u\cdot e^{-\lambda}\sim u_{0}\cdot e^{-\lambda}$ (C7)

Therefore, the effusion velocity is of the same order of magnitude as the
sound speed $u_{0}\equiv\sqrt{k_{B}T/\mu}$ attenuated by the exponential
factor $\exp{(-\lambda})$.

With the definition of $\lambda$, the critical number density of H-atom at the
exobase can be expressed as:

$n_{\text{H,exob}}=\frac{p_{\text{H,exob}}}{k_{B}T}\sim\frac{p_{\text{exob}}}{k_{B}T}=\frac{1}{k_{B}T}\cdot\frac{\mu_{\text{H}}}{\sigma}\cdot\frac{GM_{p}}{r_{\text{exob}}^{2}}=\frac{1}{u_{0}^{2}}\cdot\frac{1}{\sigma}\cdot\frac{GM_{p}}{r_{\text{exob}}^{2}}=\frac{\lambda_{\text{H}}}{\sigma\cdot
r_{\text{exob}}}$ (C8)

In all cases, it is much easier for an H-atom to escape compared to a
H2-molecule by a factor of $\sim
e^{\left(-\lambda_{\text{H}}\right)}/e^{\left(-\lambda_{\text{H}_{2}}\right)}\sim
e^{\left(\lambda_{\text{H}_{2}}-\lambda_{\text{H}}\right)}\sim
e^{\lambda_{\text{H}}}$ evaluated at the _exobase_ , which amounts to a factor
of at least a few hundred $\sim$ a few thousand times and sometimes much more.
Therefore, the escaping flux of H2 is often negligible compared to that of the
H-atom. Meanwhile, H2-molecule converts into H-atom by molecular dissociation
occurring at high Temperature and low Pressure regions in the upper
atmosphere.

##  Appendix D Hydrodynamic Escape

Applying _Bernoulli’s Principle_ (Conservation of Energy) to a compressible
flow, we have:

$\Delta\left[\underbrace{\frac{1}{2}\cdot\overline{v}^{2}}_{\text{kinetic
energy}}+\underbrace{\int\frac{dP}{\rho(P)}}_{\text{internal
energy}}+\underbrace{\Psi}_{\text{gravitational
potential}}\right]=0\text{,~{}"along a streamline"}$ (D1)

There is exchange in between the three types of energy, and therefore, the
total energy is always conserved in the flow. Here $\overline{v}$ refers to
the group (average) velocity of the flow, not of its individual constituent
atom or molecule. The random (thermal) velocity distribution is included in
the internal energy part.

The escape flux ($\Phi_{\text{hydrodynamic}}$) due to hydrodynamic escape is
of similar functional form as the Jeans flux attenuated by the same
exponential factor $\exp{(-\lambda})$, to within a factor of unity including
geometric factor:

$\Phi_{\text{hydrodynamic}}\sim n_{\text{source}}\cdot u_{0}\cdot
e^{-\lambda}$ (D2)

This is easy to understand because both the hydrodynamic escape and the Jeans
escape are end-member cases of thermal escape, of which the _Boltzmann_ factor
$\exp{(-\lambda})$ is characteristic.

The only important difference is that the number density in
$\Phi_{\text{hydrodynamic}}$ is no longer the number density of the escaping
species evaluated at the exobase, but substituted by the number density
evaluated at the origin of the hydrodynamic escape flow: $n_{\text{source}}$,
is the number density (cm-3) of the escaping species in the source region
where the hydrodynamic outflow originates. In our consideration, we equate
$n_{\text{source}}$ as the the number density of H-atom at the von Karman
line. So,

$\Phi_{\text{hydrodynamic}}\sim n_{\text{H,von Karman line}}\cdot u_{0}\cdot
e^{-\lambda}$ (D3)

Then, this $\lambda$ is evaluated at the von Karman line as well. This
$n_{\text{H,von Karman line}}$ can be viewed as sourcing an effective pressure
which drives the outflow into space.

Due to the Continuity Equation or (Conservation of Mass Flux), this outflow
must be initially sub-sonic when the number density is high, but eventually
transforms into super-sonic when the number density is low. Thus, a transonic
front is expected to occur somewhere along the streamline.

##  Appendix E H2 Collisional Dissociation

Here we consider the direct dissociation of H2 in a hydrogen-dominated
atmosphere:

$\text{H}_{2}\rightleftharpoons\text{H}+\text{H}$ (E1)

The pressure equilibrium constant $K_{p}$:

$\lg
K_{p}=\lg\frac{p_{\text{H}}^{2}\text{(atm.)}}{p_{\text{H}_{2}}\text{(atm.)}}=\left[6.16-\frac{23,500}{T(\text{K})}\right]$
(E2)

where

$\displaystyle p_{\text{H}}$ $\displaystyle=n_{\text{H}}\cdot k_{B}T$ (E3)
$\displaystyle p_{\text{H}_{2}}$ $\displaystyle=n_{\text{H}_{2}}\cdot k_{B}T$
(E4)

and

$\displaystyle n_{\text{H}}+2\cdot n_{\text{H}_{2}}$
$\displaystyle=\text{total number of H-atoms in both species, conserved}$ (E5)
$\displaystyle n_{\text{H}_{2}}+n_{\text{H}_{2}}$
$\displaystyle=n\text{,~{}total number of particles, not conserved}$ (E6)

Let

$\displaystyle n_{\text{H}}$ $\displaystyle\equiv\mathcal{X}\cdot n$ (E7)
$\displaystyle n_{\text{H}_{2}}$ $\displaystyle\equiv(1-\mathcal{X})\cdot n$
(E8)

So that,

$\lg
K_{p}=\lg\left(\frac{n_{\text{H}}^{2}\cdot(k_{B}T)^{2}}{n_{\text{H}_{2}}\cdot(k_{B}T)}\right)=\lg\left(\frac{n_{\text{H}}^{2}}{n_{\text{H}_{2}}}\cdot
k_{B}T\right)=\lg\left(\frac{\mathcal{X}^{2}}{(1-\mathcal{X})}\cdot\frac{n\cdot
k_{B}T}{\text{1 atm}}\right)$ (E9)

Thus,

$n\approx\left(\frac{(1-\mathcal{X})}{\mathcal{X}^{2}}\right)\cdot\left(\frac{\text{1
atm}}{k_{B}T}\right)\cdot
10~{}\widehat{}\left[6.16-\frac{23,500}{T(\text{K})}\right]$ (E10)

and,

$\displaystyle n_{\text{H}}\equiv\mathcal{X}\cdot n$
$\displaystyle\approx\left(\frac{1-\mathcal{X}}{\mathcal{X}}\right)\cdot\left(\frac{\text{1
atm}}{k_{B}T}\right)\cdot
10~{}\widehat{}\left[6.16-\frac{23,500}{T(\text{K})}\right]$ (E11)
$\displaystyle n_{\text{H}_{2}}\equiv(1-\mathcal{X})\cdot n$
$\displaystyle\approx\left(\frac{1-\mathcal{X}}{\mathcal{X}}\right)^{2}\cdot\left(\frac{\text{1
atm}}{k_{B}T}\right)\cdot
10~{}\widehat{}\left[6.16-\frac{23,500}{T(\text{K})}\right]$ (E12)

At the H2 dissociation front, $\left(\frac{1-\mathcal{X}}{\mathcal{X}}\right)$
is of the order of unity, and thus,

$n_{\text{H,dissociation-front}}\sim\left(\frac{\text{1
atm}}{k_{B}T}\right)\cdot
10~{}\widehat{}\left[6.16-\frac{23,500}{T(\text{K})}\right]$ (E13)

For T$\sim 10^{3}$K,

$\left(\frac{\text{1
atm}}{k_{B}T}\right)\sim\left(\frac{10^{5}\text{Pa}}{k_{B}\cdot
10^{3}\text{K}}\right)\approx 7.25\cdot 10^{18}\text{cm}^{-3}\approx
10^{18.86}\text{cm}^{-3}$ (E14)

Thus, to within a factor of unity, for T$\sim 10^{3}$K-$4\times 10^{3}$K,

$n_{\text{H,dissociation-
front}}\sim\left(\frac{1}{T/(10^{3}\text{K})}\right)\cdot
10~{}\widehat{}\left[25.0-\frac{23,500}{T(\text{K})}\right]\text{cm}^{-3}$
(E15)

Notice the strong temperature-dependence of $n_{\text{H}}$, assuming chemical
equilibrium, we have the following table:

Table 1: $n_{\text{H}}$ dependence on temperature T T (K) |  $n_{\text{H}}$(cm-3)
---|---
1000 K | $3.\times 10^{1}$
1200 K | $2.\times 10^{5}$
1500 K | $1.5\times 10^{9}$
2000 K | $1.\times 10^{13}$
3000 K | $5.\times 10^{16}$

Therefore, starting around $\sim 1500$K, the H2 dissociation front penetrates
into the exobase ($n\sim 10^{(8\sim 9)}$cm-3). Starting around $\sim 2000$K,
the H2 dissociation front penetrates deeper into the atmosphere and approaches
the level of the von Karman line ($n\sim 10^{16}$cm-3) below which different
gas species are well-mixed with a common scale-height. For comparison, the
mass density of super-critical fluid H2 is approximately
$\rho_{\text{H}_{2}\text{~{}fluid}}\sim 0.07$g cm-3 or $\sim\frac{1}{16}$g
cm-3. The corresponding number density is then
$n_{\text{H}_{2}\text{~{}fluid}}\sim 2\times 10^{22}$cm-3. Approaching this
number density, EOS needs to take into account the strong inter-molecular
repulsion.

However, a major flaw in this calculation is whether the local chemical
equilibrium can be achieved when the number densities of H2 and thus the
collisional rates fall below a certain value in the upper atmosphere. One
needs to carefully consider the reaction kinetics of a uni-molecular
decomposition of H2 through molecular collisions (Jordan, 1979):

$\begin{split}\text{H}_{2}+\text{H}_{2}&\underset{k_{2}}{\stackrel{{\scriptstyle
k_{1}}}{{\rightleftharpoons}}}\text{H}_{2}^{*}+\text{H}_{2}\\\
\text{H}_{2}^{*}&\stackrel{{\scriptstyle
k_{d}}}{{\rightarrow}}\text{H}+\text{H}\\\
\hline\cr\text{\emph{net}~{}~{}~{}}\text{H}_{2}&\rightarrow\text{H}+\text{H}\end{split}$
(E16)

A H2 molecule is activated by collision with another H2 molecule. In the
collisional process the relative kinetic energy of the two colliding molecules
is transformed, in part, to (vibrational) internal energy of one of the
molecules, and thus, activating it. Such inelastic collisions provide the
basis for understanding the activation process. The excited molecule
H${}_{2}^{*}$ either de-excites through another collision, or, decomposes into
the products—atomic hydrogen.

Above the von Karman line ($n\sim 10^{16}$cm-3) the gas becomes highly
rarefied so that vibrational energy of a molecule decouples from the thermal
equilibrium, because a large number of collisions is necessary to excite the
vibrational degree of freedom. On the other hand, rotational degree of freedom
remains in equilibrium to much lower density because it is more easily
perturbed by collisions. Therefore, it is more reasonable to evaluate the
number density of H-atom ($n_{\text{H}}$) at the von Karman line where the
H2-molecule ($n_{\text{H}_{2}}$) is still the dominating species and the
collisions among them are frequent enough. Recall that,

$\lg
K_{p}=\lg\left(\frac{n_{\text{H}}^{2}\cdot(k_{B}T)^{2}}{n_{\text{H}_{2}}\cdot(k_{B}T)}\right)=\lg\left(\frac{n_{\text{H}}^{2}}{n_{\text{H}_{2}}}\cdot\frac{k_{B}T}{\text{1
atm}}\right)=\left[6.16-\frac{23,500}{T(\text{K})}\right]$ (E17)

Thus,

$n_{\text{H}}^{2}\sim\left(n_{\text{H}_{2}}\cdot\frac{\text{1
atm}}{k_{B}T}\right)\cdot
10~{}\widehat{}\left[6.16-\frac{23,500}{T(\text{K})}\right]$ (E18)

Thus, for T$\sim 10^{3}$K-$4\times 10^{3}$K,

$n_{\text{H}}^{2}\sim\left(\frac{n_{\text{H}_{2}}}{T/(10^{3}\text{K})}\right)\cdot
10~{}\widehat{}\left[25.0-\frac{23,500}{T(\text{K})}\right]$ (E19)

Now, for convenience of calculation, let:

$T_{3}\equiv T/(10^{3}\text{K})$ (E20)

Then,

$n_{\text{H}}^{2}\sim\left(\frac{n_{\text{H}_{2}}}{T_{3}}\right)\cdot
10~{}\widehat{}\left[25.0-\frac{23.5}{T_{3}}\right]$ (E21)

So,

$n_{\text{H}}\sim\sqrt{\left(\frac{n_{\text{H}_{2}}}{T_{3}}\right)\cdot
10~{}\widehat{}\left[25.0-\frac{23.5}{T_{3}}\right]}\text{cm}^{-3}$ (E22)

Equivalently,

$n_{\text{H}}\sim\sqrt{\frac{n_{\text{H}_{2}}}{T_{3}}}\cdot
10~{}\widehat{}\left[12.5-\frac{11.75}{T_{3}}\right]\text{cm}^{-3}$ (E23)

Finally,

$n_{\text{H,von Karman line}}\sim\sqrt{\frac{n_{\text{H}_{2\text{,von Karman
line}}}}{T_{3}}}\cdot
10~{}\widehat{}\left[12.5-\frac{11.75}{T_{3}}\right]\text{cm}^{-3}$ (E24)

Plug in the number density estimated at the von Karman line
$n_{\text{H}_{2\text{,von Karman line}}}\sim 10^{16}\text{~{}cm}^{-3}$,

$n_{\text{H,von Karman line}}\sim\frac{1}{\sqrt{T_{3}}}\cdot
10~{}\widehat{}\left[20.5-\frac{11.75}{T_{3}}\right]\text{cm}^{-3}$ (E25)

Note that this formula will inevitably fail at some high temperature when the
H2 dissociation front penetrates deep enough so that the number density of the
escape flux is no longer limited by the temperature itself but instead limited
by the energy input. Because now the escape flux produces a non-negligible
effect of _cooling_ , this _cooling_ has to be compensated by a significant
portion of the energy input of host-stellar incident radiation energy. In
other words, breaking H-H bond provides effective cooling. At this point, the
escape flow saturates and transitions from an _entropy-limited_ to an _energy-
limited_ regime. We will re-visit this point when we discuss the energetics,
and show that this is a _phase transition_ from ordinary thermal escape
scenario to _shock-front_ -driven hydrodynamic high-speed outflow.

##  Appendix F H2 Photo-Dissociation

The other possibility is the ability of energetic host-stellar photons with
wavelengths less than 277 nm (corresponding to H-H bond strength of 4.5 eV) to
rupture the strong H-H bond in molecular hydrogen. The sequence of events is:

$\begin{split}\text{H}_{2}+\gamma&\rightarrow\text{H}_{2}^{*}\\\
\text{H}_{2}^{*}&\rightarrow\text{H}+\text{H}\\\
\hline\cr\text{\emph{net}~{}~{}~{}}\text{H}_{2}+\gamma&\rightarrow\text{H}+\text{H}\end{split}$
(F1)

$\gamma$ is a UV photon with wavelength shorter than 277 nm. The decay of the
excited state H${}_{2}^{*}$ is rapid, and the reaction that will be observed
is the _net_ reaction. Thus, the rate equation becomes,

$\frac{d[\gamma]}{dt}=\frac{d[\text{H}_{2}]}{dt}=-\frac{1}{2}\cdot\frac{d[\text{H}]}{dt}=-J_{\text{H}_{2}}\cdot[\text{H}_{2}]$
(F2)

where the square brackets indicate concentrations. $J_{\text{H}_{2}}$ is the
_photochemical loss rate_ of molecular hydrogen per hydrogen molecule or the
_inverse lifetime_ of a hydrogen molecule for photo-dissociation for our
purpose.

The rate of change of photon density is the rate of change of radiant energy
density divided by $(h\cdot\nu)$. Thus,

$\displaystyle\frac{d[\gamma]}{dt}$
$\displaystyle=\int_{\nu_{0}}^{\infty}\frac{d\nu}{h\nu}\,\Big{(}\mathbf{\nabla}\cdot\mathbf{F}_{\nu}\Big{)}$
($\mathbf{F}_{\nu}$ is the monochromatic flux)
$\displaystyle=\int_{\nu_{0}}^{\infty}\frac{d\nu}{h\nu}\,\bigg{(}\int_{4\pi}d{\omega}_{\vec{l}}\cdot[\text{H}_{2}]\cdot\sigma_{\nu}\cdot(-I_{\nu})\bigg{)}$
(integrate over all solid angle of direction $\vec{l}$) (F3)
$\displaystyle=-[\text{H}_{2}]\cdot\underbrace{\int_{\nu_{0}}^{\infty}\frac{d\nu}{h\nu}\,\bigg{(}\int_{4\pi}d{\omega}_{\vec{l}}\cdot\sigma_{\nu}\cdot
I_{\nu}\bigg{)}}_{\text{photochemical loss rate constant $J_{\text{H}_{2}}$}}$
($[\text{H}_{2}]$ is the number density of H2) (F4)
$\displaystyle=-[\text{H}_{2}]\cdot\underbrace{\int_{\nu_{0}}^{\infty}\frac{d\nu}{h\nu}\,\bigg{(}\sigma_{\nu}\cdot
f_{\nu,0}\cdot\exp{\bigg{(}-\int_{z}^{\infty}\sigma_{\nu}(z^{\prime})\cdot
n(z^{\prime})\cdot\frac{dz^{\prime}}{\xi_{\star}}\bigg{)}}\bigg{)}}_{\text{photochemical
loss rate constant $J_{\text{H}_{2},d}\Big{(}z\Big{)}$ evaluated at height
$z$}}$ ($f_{\nu,0}$ is host-stellar irradiance before attenuation) (F5)

$\nu_{0}$ is the minimum frequency that can give rise to a dissociation,
corresponding to 277 nm for molecular hydrogen. $\sigma_{\nu}$ is the
monochromatic _absorption coefficient per molecule_ which by definition is
again a collisional cross-section, but this time for a photon of negligible
size to interact with a molecule with an effective frequency-dependent
collisional cross-section. Thus, $\sigma_{\nu}$ is the cross-section for
light-matter interaction compared to $\sigma$ for molecular collision.
Usually,

$\sigma_{\nu}<\sigma$ (F6)

$I_{\nu}$ is the monochromatic radiance which originates from the direct host-
stellar beam within the small solid angle ($\Delta\omega_{\star}$) of host-
stellar disc but has already been attenuated by the planetary atmosphere from
$\infty$ until height $z$. Thus, before any planet atmospheric attenuation, we
have,

$\displaystyle f_{\nu,0}$ $\displaystyle=I_{\nu,0}\cdot\Delta\omega_{\star},$
where $\Delta\omega_{\star}$ is the small solid angle of host-stellar disc
viewed from planet (F7) $\displaystyle\xi_{\star}$
$\displaystyle=\cos{\theta_{\star}},$ where $\theta_{\star}$ is the zenith
angle subtended by the host-star as viewed from a location on planet. (F8)

The frequency-dependent absorption cross-section $\sigma_{\nu}$ for H2 and He
from 0-300 nm wavelength is needed to model the photo-dissociation process—the
absorption bands/spectrum of H2 and He in the UV. Fortunately, H2 ceases to
absorb at 100 nm, and He, 50 nm (Lewis, 2004; Draine, 2011). If the product
from photo-dissociation collides with many other atmospheric species and
shares its energy widely, then its energy enters the thermal reservoir and is
considered thermalized.

##  Appendix G H2 Opacity

In order for the dissociation to occur, the energy must be deposited into the
vibrational mode (degree-of-freedom) of H2-molecule, through collisions, but
NOT through absorption of visible or (thermal) infrared photons. A symmetric
diatomic molecule such as H2 has zero electric-dipole moment. Its dimension
($\sim 1\AA$) is much smaller compared to the wavelength of an infrared photon
($\sim 10^{4}\AA$), so the effect of the changing E-field is seen as the same
by both H-atoms in a H2-molecule. So it cannot absorb an infrared photon even
if the photon is precisely at the natural vibration frequency of the molecule.
Thus, the direct interaction of H2-molecule with an infrared photon is
_forbidden_.

On the other hand, electronic transitions of H2 lie in the ultra-violet part
of the spectrum (Draine, 2011). Thus, H2 is mostly transparent to the
planetary radiation field which emits infrared photons, as well as the stellar
radiation field except at discrete lines due to the electronic transitions in
the ultra-violet, which can only provide very limited opacity. Thus, the
stellar bolometric radiation can penetrate all the way down to great depth to
reach $\gtrsim 1$ bar level ($n\sim 10^{19}$cm-3) where the pressure-induced
(i.e. collisionally-induced) absorption cross-section and the Rayleigh
scattering cross-section for H2 become significant (Hayashi et al., 1979;
Goody, 1995; Lewis, 2004). There H2-fluid becomes dense enough to absorb all
incoming photons as a continuum. This level would define the photo-sphere of
such a planet.

Now, if there is trace amount of polar molecules with non-zero dipole moments
such as H2O and CO mixed in with the H2-gas, then, H2O and CO molecules can
absorb the infrared photons through their vibrational mode, and collisionally
cascade this energy to the vibrational and rotational mode of H2-molecule.
Depending on the mixing ratio, the presence of these polar molecules can
contribute significant opacity and bring the photo-sphere to a shallower depth
and lower pressure level. However, all this can only occur in the lower
atmosphere (homosphere) where the mean free path ($\lambda$) is short enough
($\lambda\lesssim 1$cm) to allow effective mixing of gas species of different
molecular weights. This $\lambda\sim 1$cm roughly corresponds to the pressure
level of $\gtrsim 1$ Pascal to the order-of-magnitude. This opacity
($\kappa_{\text{polar-molecules}}$) contributed by H2O and CO molecules is
important for the temperature range of 1500K-3000K, and can be crudely given
by:

$\kappa_{\text{polar-molecules}}\approx 0.1\cdot
Z,~{}~{}~{}[\text{cm}^{2}\text{g}^{-1}]$ (G1)

The strongest absorber, the one with the greatest dipole moment and also the
most abundant among the available polar molecules, is always H2O gas (vapor).
On the other hand, the opacity ($\kappa_{\text{H}^{-}}$) contributed by the
presence of negative hydrogen ion, H-, only becomes important at higher
temperatures of 4000K-8000K, such as in stellar atmospheres of Main-Sequence
(MS) stars (Schwarzschild, 1958):

$\kappa_{\text{H}^{-}}\approx 1.1\times 10^{-25}\cdot
Z^{0.5}\cdot\rho^{0.5}\cdot T^{7.7},~{}~{}~{}[\text{cm}^{2}\text{g}^{-1}]$
(G2)

where $X$, $Y$, and $Z$ denote the hydrogen, helium, and heavy element
fractions by mass, respectively.

##  Appendix H Energetic View

Generally, for thermal escape, it is particle which lies in the energetic tail
of the _Boltzmann_ distribution that can overcome the energy barriers and
escape. For our purpose, the energy barriers refer to (1) the barrier to break
the H-H bond, and (2) the barrier to escape the gravitational potential of the
planet.

These particles which lie in the energetic tail of the _Boltzmann_
distribution have much higher energy compared to the average thermal energy
($k_{B}T$) of the ensemble. Thus, when they escape, they also carry away their
high energy and result in _evaporative cooling_ of the planet surface and
planet atmosphere. Ultimately, this energy loss due to escape is supplied by a
fraction of incoming stellar radiation. This process can be compared to the
kinetics and energetics of evaporation (Lu et al., 2019). Therefore, it is
worthwhile to look at the energetics of the escape process. A simple back-of-
envelope calculation would illuminate this point. First, let’s consider the
bolometric solar luminosity:

$L_{\odot}=3.828\times 10^{26}\text{J/s}$ (H1)

Earth only intercepts a small fraction of this energy flow:

$L_{\oplus}=L_{\odot}\cdot\frac{\pi R_{\oplus}^{2}}{4\pi
a^{2}}=\frac{L_{\odot}}{4\cdot(a/R_{\oplus})^{2}}=1.75\times
10^{17}\text{J/s}$ (H2)

If all this energy goes into dissociate H2 into H-atoms, we can estimate how
long it takes to dissociate 1 Earth mass ($M_{\oplus}$) equivalent of H2:

$\tau=\frac{\Delta H\cdot
M_{\oplus}/\mu_{\text{H}_{2}}}{L_{\oplus}}=\frac{440\text{kJ/mol}\times 6\cdot
10^{24}\text{kg}/(2\text{g/mol})}{1.75\times
10^{17}\text{J/s}}=\frac{1.32\times 10^{33}\text{J}}{1.75\times
10^{17}\text{J/s}}=0.75\times 10^{16}\text{s}\approx 0.2\times
10^{9}\text{years}$ (H3)

Therefore, at Earth’s distance, it would take $\sim$0.2 billion years (our
estimate) to dissociate one Earth-mass equivalent of H2 into H-atoms, if all
the energy that Earth receives from the Sun over this duration goes into this
process. Marvelously, this timescale is shorter but not too different from the
age of our solar system!

In reality, only a fraction of the incoming stellar radiation can contribute
to the H-H bond dissociation. Thus, it would be extremely important to
determine this fraction. Now, many exoplanets we consider are significantly
H2-rich compared to our own Earth, that some of them do have $\sim
1M_{\oplus}$ or more equivalent of H2 mass stored in their envelopes. On the
other hand, they are also significantly closer to their host stars, so they
typically receive a hundred-fold ($f/f_{\oplus}\sim$100) or more of incident
stellar radiation compared to our own Sun. Thus, we consider this fraction to
be of the order of percent-level, in order for these planets to last over
billion year timescale.

Therefore, we need to consider the following energy balance equation, for
instance, for the region near the sub-stellar hot spot where the escape is
expected to strongest, and ignore any albedo effect for now:

$\underbrace{f_{0}}_{\text{absorbed stellar
flux}}=\underbrace{\left(\frac{\Delta
H}{\mu_{\text{H}_{2}}}+\frac{GM_{p}}{R_{p}}\right)\cdot\mu_{\text{H}}\cdot\Phi_{\text{H,escape}}}_{\text{(A):
evaporative cooling by escaping H}}+\underbrace{\sigma\cdot T^{4}}_{\text{(B):
radiative cooling by escaping photons}}$ (H4)

Therefore, the _cooling_ in order to balance the incoming stellar radiation
can be accomplished in two ways: (A) through escaping particles (H-atoms
considered here) which carry away energy, (B) through escaping photons
(thermal infrared radiation considered here) which carry away energy. The key
difference between (A) and (B) is that for (B) photons can readily escape from
planet surface, while for (A) particles need to overcome significant energy
barriers so only a very small fraction of them with high enough energy can
escape.

(A) in turn contains two parts: (A1) energy to break up H-H bond, (A2) energy
to overcome gravitational potential.

There is a very intricate and subtle relation between term (A) and term (B),
because the escaping flux $\Phi_{\text{H,escape}}$ in (A) depends (strongly)
on T (thermally driven in our assumption) in (B).

Thus, there must exist a feedback mechanism between (A) and (B). This feedback
mechanism can be understood as follows: if temperature T increases a little
bit, then $\Phi_{\text{H,escape}}$ would increase rapidly (because it is
strongly temperature-dependent). This rapid increase in
$\Phi_{\text{H,escape}}$ would result in even more cooling which would help
bring down T, and vice versa. We consider this feedback mechanism acts as a
_peak-holding valve_ which controls the escape rate.

In other words, we need to figure out the functional dependence of escape flux
$\Phi_{\text{H,escape}}$ on temperature T:

$\Phi_{\text{H,escape}}=\Phi_{\text{H,escape}}(T)$ (H5)

Let’s turn back to the equation of thermal (hydrodynamic) escape. Recall the
earlier formula:

$n_{\text{H,dissociation-
front}}\sim\left(\frac{10^{25}\text{cm}^{-3}}{T/(10^{3}\text{K})}\right)\cdot
10~{}\widehat{}\left[-\frac{23,500}{T(\text{K})}\right]$ (H6)

equivalently, converting to natural-based logarithm by multiplying a factor of
$\ln{10}\approx 2.303$ inside the power-index,

$n_{\text{H,dissociation-
front}}\sim\left(\frac{10^{25}\text{cm}^{-3}}{T/(10^{3}\text{K})}\right)\cdot\exp{\left[2.303\cdot\bigg{(}-\frac{23,500}{T(\text{K})}\bigg{)}\right]}$
(H7)

$n_{\text{H,dissociation-
front}}\sim\left(\frac{10^{25}\text{cm}^{-3}}{T/(10^{3}\text{K})}\right)\cdot\exp{\left[-\frac{5.4\times
10^{4}}{T(\text{K})}\right]}$ (H8)

On the other hand,

$n_{\text{H,von Karman
line}}\sim\left(\frac{10^{20.5}\text{cm}^{-3}}{\sqrt{T/(10^{3}\text{K})}}\right)\cdot
10~{}\widehat{}\left[-\frac{11,750}{T(\text{K})}\right]$ (H9)

equivalently, converting to natural-based logarithm by multiplying a factor of
$\ln{10}\approx 2.303$ inside the power-index,

$n_{\text{H,von Karman
line}}\sim\left(\frac{10^{20.5}\text{cm}^{-3}}{\sqrt{T/(10^{3}\text{K})}}\right)\cdot\exp{\left[2.303\cdot\bigg{(}-\frac{11,750}{T(\text{K})}\bigg{)}\right]}$
(H10)

$n_{\text{H,von Karman
line}}\sim\left(\frac{10^{20.5}\text{cm}^{-3}}{\sqrt{T/(10^{3}\text{K})}}\right)\cdot\exp{\left[-\frac{2.7\times
10^{4}}{T(\text{K})}\right]}$ (H11)

The reason for this conversion is show that the numerator inside the power-
index is just a measure of the bond strength of H-H covalent bond, that is,
the enthalpy of H-H covalent bond dissociation $\Delta H\approx 440$kJ/mol
$\approx 100$kcal/mol divided by the gas constant
R$=8.314\text{J}\text{mol}^{-1}\text{K}^{-1}$:

$\frac{\Delta H}{\text{R}}=\frac{440\times
10^{3}\text{J}\text{mol}^{-1}}{8.314\text{J}\text{mol}^{-1}\text{K}^{-1}}\approx
5.4\times 10^{4}\text{K}$ (H12)

(1) Thus, per bond $\Delta E_{\text{H-H bond}}\approx 4.5$eV:

$\Delta E_{\text{H-H bond}}\text{(eV)}=\frac{\Delta H}{N_{A}\times
e}=\frac{440\times 10^{3}\text{J}\text{mol}^{-1}}{6.022\cdot 10^{23}\times
1.602\times 10^{-19}\text{coulombs}}\approx 4.5\text{eV}$ (H13)

This is typical value for a covalent bond (Speight & Lange, 2016). This value
is almost exactly one-third (1/3) of the ionization energy of an H-atom:
$\Delta E_{\text{H-ionization}}\approx 13.5$eV.

(2) The energy required for one H-atom to overcome the gravitational potential
of the planet is:

$\Delta E_{\text{grav}}\text{(eV)}=\frac{GM_{p}\cdot\mu_{\text{H}}}{R_{p}\cdot
e}=\frac{M_{p}/M_{\oplus}}{R_{p}/R_{\oplus}}\cdot\frac{(11.2\cdot
10^{3}\text{m/s})^{2}/2\times 1.6735\cdot 10^{-27}\text{kg}}{1.602\times
10^{-19}\text{coulombs}}\approx\frac{M_{p}/M_{\oplus}}{R_{p}/R_{\oplus}}\cdot
0.65\text{eV}\sim 1.\text{eV}$ (H14)

(3) In contrast, the thermal energy $k_{B}T=k_{B}\cdot(10^{3}\text{K})\approx
0.086$eV$\sim 0.1$eV,

$\Delta
E_{\text{thermal}}\text{(eV)}=\frac{k_{B}T}{e}=\frac{T}{10^{3}\text{K}}\cdot\frac{1.38\cdot
10^{-23}\text{J}\text{K}^{-1}\times 10^{3}\text{K}}{1.602\times
10^{-19}\text{coulombs}}\approx\frac{T}{10^{3}\text{K}}\cdot
0.086\text{eV}\sim 0.1\text{eV}$ (H15)

which is about two orders-of-magnitude smaller compared to the energy needed
to break the H-H bond, and about one order-of-magnitude smaller compared to
the energy needed to escape the gravitational potential well of the planet.

Thus, we perceive the _hierarchy_ of energy here, and re-confirms our earlier
point that only the energetic tail of the thermal distribution can escape.
Recall that,

$\Phi_{\text{H-atom, hydrodynamic escape}}\sim n_{\text{H,von Karman
line}}\cdot u_{0}\cdot e^{-\lambda}$ (H16)

Therefore, by collecting the terms for the two energy barriers together, we
can write the escape flux generally as:

$\Phi_{\text{H,escape}}\sim\left(\frac{10^{20.5}\text{cm}^{-3}}{\sqrt{T/(10^{3}\text{K})}}\right)\cdot
u_{0}\cdot\exp{\left(-\frac{\Delta E_{\text{H-H bond}}/2+\Delta
E_{\text{grav}}}{\Delta E_{\text{thermal}}}\right)}$ (H17)

plug-in the isothermal sound speed expression
$u_{0}\equiv\sqrt{k_{B}T/\mu_{\text{H}}}\approx\sqrt{T/(10^{3}\text{K})}\cdot
3$km/s, we have

$\Phi_{\text{H,escape}}\sim\left(\frac{10^{20.5}\text{cm}^{-3}}{\sqrt{T/(10^{3}\text{K})}}\right)\cdot\sqrt{T/(10^{3}\text{K})}\cdot\left(3\text{km/s}\right)\cdot\exp{\left(-\frac{\Delta
E_{\text{H-H bond}}/2+\Delta E_{\text{grav}}}{k_{B}T}\right)}$ (H18)

or, equivalently,

$\Phi_{\text{H,escape}}\sim\bigg{(}1.0\times
10^{26}\text{cm}^{-2}\text{s}^{-1}\bigg{)}\cdot\exp{\left(-\frac{\Delta
E_{\text{H-H bond}}/2+\Delta E_{\text{grav}}}{k_{B}T}\right)}$ (H19)

To simplify the expression, we can write H-atom number escape flux
$\Phi_{\text{H,escape}}$ as:

$\Phi_{\text{H,escape}}\sim
C\cdot\exp{\left(-\frac{a/2+b\cdot(m/r)}{T}\right)}$ (H20)

where,

$\displaystyle a$ $\displaystyle\equiv\frac{E_{\text{H-H bond}}}{k_{B}}\approx
5.4\times 10^{4}\text{K}$ (H21) $\displaystyle b$
$\displaystyle\equiv\frac{E_{\text{grav}}}{k_{B}}\approx 7.6\times
10^{3}\text{K}$ (H22) $\displaystyle m$ $\displaystyle\equiv M_{p}/M_{\oplus}$
(H23) $\displaystyle r$ $\displaystyle\equiv R_{p}/R_{\oplus}$ (H24)
$\displaystyle C$ $\displaystyle\sim\bigg{(}1.0\times
10^{26}\text{~{}cm}^{-2}\text{s}^{-1}\bigg{)}$ (H25)

Now, we can convert this expression to the mass flux of H-atom escape as:

$\mu_{\text{H}}\cdot\Phi_{\text{H,escape}}\sim
C^{\prime}\cdot\exp{\left(-\frac{a/2+b\cdot(m/r)}{T}\right)}$ (H26)

where $\mu_{\text{H}}=1$ g/mol, so that,

$C^{\prime}\sim\bigg{(}1.66\times
10^{2}\text{~{}g}\cdot\text{cm}^{-2}\text{s}^{-1}\bigg{)}$ (H27)

Now, the mass loss rate can be estimated by multiplying the mass flux with an
effective surface area of escape ($\eta\cdot R_{p}^{2}$) where $\eta$ is a
geometric factor of the order of unity depending on the geometry of the escape
flow:

$\dot{M_{p}}=\mu_{\text{H}}\cdot\Phi_{\text{H,escape}}\cdot\eta\cdot
R_{p}^{2}$ (H28)

We can non-dimensionalize both sides of the equation so that,

$\dot{m}=\mu_{\text{H}}\cdot\Phi_{\text{H,escape}}\cdot\frac{R_{\oplus}^{2}}{M_{\oplus}}\cdot\eta\cdot
r^{2}$ (H29)

Now,

Recall the energy balance equation of the absorbed stellar flux ($f_{0}$). In
term (A), we have:

$\displaystyle\frac{\Delta H}{\mu_{\text{H}_{2}}}$ $\displaystyle=\frac{a\cdot
k_{B}\cdot
N_{A}}{\mu_{\text{H}_{2}}}=\frac{a\cdot\text{R}}{\mu_{\text{H}_{2}}}=\frac{a\cdot\text{R}}{2\cdot\mu_{\text{H}}}=\frac{a}{2}\cdot\frac{\text{R}}{\mu_{\text{H}}}=\frac{a}{2}\cdot\left(8.314\text{~{}J}\text{~{}g}^{-1}\text{K}^{-1}\right)$
(H30) $\displaystyle\frac{GM_{p}}{R_{p}}$ $\displaystyle=\frac{b\cdot
k_{B}\cdot
N_{A}}{\mu_{\text{H}}}\cdot\frac{m}{r}=\frac{b\cdot\text{R}}{\mu_{\text{H}}}\cdot\frac{m}{r}=b\cdot\frac{m}{r}\cdot\frac{\text{R}}{\mu_{\text{H}}}=b\cdot\frac{m}{r}\cdot\left(8.314\text{~{}J}\text{~{}g}^{-1}\text{K}^{-1}\right)$
(H31)

So that,

$\underbrace{f_{0}}_{\text{absorbed stellar
flux}}=\underbrace{\left(8.314\text{~{}J}\text{~{}g}^{-1}\text{K}^{-1}\right)\cdot\left(\frac{a}{2}+b\cdot\frac{m}{r}\right)\cdot\mu_{\text{H}}\cdot\Phi_{\text{H,escape}}}_{\text{(A):
evaporative cooling by escaping H}}+\underbrace{\sigma\cdot T^{4}}_{\text{(B):
radiative cooling by escaping photons}}$ (H32)

$\underbrace{f_{0}}_{\text{absorbed stellar
flux}}=\underbrace{\left(8.314\text{~{}J}\text{~{}g}^{-1}\text{K}^{-1}\right)\cdot\left(\frac{a}{2}+b\cdot\frac{m}{r}\right)\cdot
C^{\prime}\cdot\exp{\left(-\frac{a/2+b\cdot(m/r)}{T}\right)}}_{\text{(A):
evaporative cooling by escaping H}}+\underbrace{\sigma\cdot T^{4}}_{\text{(B):
radiative cooling by escaping photons}}$ (H33)

$\underbrace{f_{0}}_{\text{absorbed stellar
flux}}=\underbrace{\bigg{(}1.38\times
10^{3}\text{~{}J}\text{~{}cm}^{-2}\text{K}^{-1}\text{s}^{-1}\bigg{)}\cdot\left(\frac{a}{2}+b\cdot\frac{m}{r}\right)\cdot\exp{\left(-\frac{a/2+b\cdot(m/r)}{T}\right)}}_{\text{(A):
evaporative cooling by escaping H}}+\underbrace{\sigma\cdot T^{4}}_{\text{(B):
radiative cooling by escaping photons}}$ (H34)

$\underbrace{f_{0}}_{\text{absorbed stellar
flux}}=\underbrace{\bigg{(}1.38\times
10^{7}\text{~{}W}\text{~{}m}^{-2}\text{K}^{-1}\bigg{)}\cdot\left(\frac{a}{2}+b\cdot\frac{m}{r}\right)\cdot\exp{\left(-\frac{a/2+b\cdot(m/r)}{T}\right)}}_{\text{(A):
evaporative cooling by escaping H}}+\underbrace{\sigma\cdot T^{4}}_{\text{(B):
radiative cooling by escaping photons}}$ (H35)

$\underbrace{f_{0}}_{\text{absorbed stellar
flux}}=\underbrace{\bigg{(}1.38\times
10^{7}\text{~{}W}\text{~{}m}^{-2}\text{K}^{-1}\bigg{)}\cdot\left(\frac{a}{2}+b\cdot\frac{m}{r}\right)\cdot\exp{\left(-\frac{a/2+b\cdot(m/r)}{T}\right)}}_{\text{(A):
evaporative cooling by escaping H}}+\underbrace{\left(5.67\times
10^{-8}\text{~{}W}\text{~{}m}^{-2}\text{K}^{-4}\right)\cdot T^{4}}_{\text{(B):
radiative cooling by escaping photons}}$ (H36)

Now, for convenience of calculation, let’s define:

$T_{3}\equiv T/(10^{3}\text{K})$ (H37)

and,

$\displaystyle a_{3}$ $\displaystyle\equiv a/(10^{3}\text{K})=54.$ (H38)
$\displaystyle b_{3}$ $\displaystyle\equiv b/(10^{3}\text{K})=7.6$ (H39)

so that,

$\underbrace{f_{0}}_{\text{absorbed stellar
flux}}=\underbrace{\left(1.38\times
10^{10}\text{~{}W}\text{~{}m}^{-2}\right)\cdot\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3}}\right)}}_{\text{(A):
evaporative cooling by escaping H-atoms}}+\underbrace{\left(5.67\times
10^{4}\text{~{}W}\text{~{}m}^{-2}\right)\cdot T_{3}^{4}}_{\text{(B): radiative
cooling by escaping photons}}$ (H40)

This is a _transcendental_ equation if one wants to exactly solve for $T_{3}$
for a given $f_{0}$. It can be solved by recursive method, similar to the way
to solve the _Kepler equation_ : $E-e\sin{E}=M$.

However, one can quickly estimate at what $T_{3\text{,crit}}$ term (A) and
term (B) become comparable in magnitude:

$\left(1.38\times
10^{10}\text{~{}W}\text{~{}m}^{-2}\right)\cdot\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit}}}\right)\sim\left(5.67\times
10^{4}\text{~{}W}\text{~{}m}^{-2}\right)\cdot T^{4}_{3\text{,crit}}}$ (H41)

so,

$\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit}}}\right)\sim\left(4.1\times
10^{-6}\right)\cdot T^{4}_{3\text{,crit}}}$ (H42)

or, equivalently,

$\left(\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit}}}\right)\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit}}}\right)\sim\left(4.1\times
10^{-6}\right)\cdot T^{3}_{3\text{,crit}}}$ (H43)

if we take natural logarithm on both sides:

$\ln{\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)}-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit}}}\sim-12.4+4\cdot\ln{T_{3\text{,crit}}}$
(H44)

re-arranging the terms, we have:

$\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit}}}\sim\ln{\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)}+12.4-4\cdot\ln{T_{3\text{,crit}}}$
(H45)

finally,

$T_{3\text{,crit}}\sim\frac{a_{3}/2+b_{3}\cdot(m/r)}{\ln{\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)}+12.4-4\cdot\ln{T_{3\text{,crit}}}}$
(H46)

by simple inspection, one can see that $T_{3\text{,crit}}\sim 2$ to 3 because
the dominating term in the numerator is $a_{3}/2\approx 27.$, and the
dominating term in the denominator is $12.4$. Also, out of coincidence, the
term $\ln{\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)}$ and the term
$\left(4.\cdot\ln{T_{3\text{,crit}}}\right)$ roughly cancel each other out in
magnitude. Therefore, a crude approximation gives:

$T_{3\text{,crit}}\sim\frac{a_{3}/2+b_{3}\cdot(m/r)}{12.4}$ (H47)

plug in the numbers: $a_{3}=54$ and $b_{3}=7.6$, we have

$T_{3\text{,crit}}\sim\frac{27+7.6\cdot(m/r)}{12.4}$ (H48)

or,

$T_{3\text{,crit}}\sim 2.18+\underbrace{0.61\cdot\frac{m}{r}}_{\text{due to
gravity potential}}$ (H49)

Now, the physical intuition is that the second term due to gravity potential
can be dropped when it approaches the critical scenario, because as the escape
process occurs it may _absorb_ more energy in various modes along the way of
escape, and the escape process may be accelerated by a _shock wave front_ ,
that is, there is a _phase transition_ in the flow. We will re-visit this
point when we discuss the _tropopause_. Thus,

$T_{3\text{,crit}}\sim 2.18$ (H50)

Therefore, we can assume that the critical scenario occurs right around the
transition from the _entropy-limited_ regime into the _energy-limited_ regime.
This coincides with the observed cut-off temperature for large exoplanets at
$T_{3\text{,crit}}\sim 2$. On the other hand, the second term
$\left(\sim\frac{m}{r}\right)$ may be invoked to explain the ordinary
positive-slope boundary bordering _hot Saturns_ in lg-lg mass-radius plot
under normal non-limiting non-critical thermal escape scenario. Thus, in this
respect of the critical _phase transition_ of flow status, these two terms can
be separated.

###  H.1 Tropopause

Tropopause is the level which separates the troposphere, where convective
upwelling dominates the energy transport mode, and the stratosphere, where
radiative cooling dominates the energy transport mode. Because the visible
opacity of H2 gas goes like $\sim P^{2}$.

The tropopause is typically located at ten percent the pressure level of
photosphere. Say, if photosphere is at 1 bar level, then tropopause is
typically at 100 mb level. The thermal flux involved in the stratosphere and
all the atmospheric mass above the tropopause is typically on the order of one
percent level ($0.1^{2}=0.01$) of the total thermal flux received from the
host star.

So, one should also estimate at what $T_{3\text{,crit2}}$ term (A) equals one
percent (1%) of term (B) in magnitude:

$\left(1.38\times
10^{10}\text{~{}W}\text{~{}m}^{-2}\right)\cdot\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit2}}}\right)\sim
0.01\times\left(5.67\times 10^{4}\text{~{}W}\text{~{}m}^{-2}\right)\cdot
T^{4}_{3\text{,crit2}}}$ (H51)

so,

$\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit2}}}\right)\sim\left(4.1\times
10^{-8}\right)\cdot T^{4}_{3\text{,crit2}}}$ (H52)

or, equivalently,

$\left(\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit2}}}\right)\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit2}}}\right)\sim\left(4.1\times
10^{-8}\right)\cdot T^{3}_{3\text{,crit2}}}$ (H53)

if we take natural logarithm on both sides:

$\ln{\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)}-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit2}}}\sim-17.0+4\cdot\ln{T_{3\text{,crit2}}}$
(H54)

re-arranging the terms, we have:

$\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3\text{,crit2}}}\sim\ln{\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)}+17.0-4\cdot\ln{T_{3\text{,crit2}}}$
(H55)

finally,

$T_{3\text{,crit2}}\sim\frac{a_{3}/2+b_{3}\cdot(m/r)}{\ln{\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)}+17.0-4\cdot\ln{T_{3\text{,crit2}}}}$
(H56)

again, by simple inspection, one can see that $T_{3\text{,crit2}}\sim 2$
because the dominating term in the numerator is $a_{3}/2\approx 27.$, and the
dominating term in the denominator is $17.0$.

Also, out of coincidence, the term
$\ln{\left(\frac{a_{3}}{2}+b_{3}\cdot\frac{m}{r}\right)}$ and the term
$\left(4.\cdot\ln{T_{3\text{,crit2}}}\right)$ roughly cancel each other out in
magnitude. Therefore, a crude approximation gives:

$T_{3\text{,crit2}}\sim\frac{a_{3}/2+b_{3}\cdot(m/r)}{17.0}$ (H57)

plug in the numbers: $a_{3}=54.0$ and $b_{3}=7.6$, we have

$T_{3\text{,crit2}}\sim\frac{27+7.6\cdot(m/r)}{17.0}$ (H58)

or,

$T_{3\text{,crit2}}\sim 1.58+0.44\cdot\frac{m}{r}$ (H59)

Similar argument as in the previous discussion can be invoked to drop the
second term so that,

$T_{3\text{,crit2}}\sim 1.58$ (H60)

Thus, this $T_{3\text{,crit2}}$ is close to the observed actual cut-off
flux/temperature of small exoplanets at $T_{3\text{,crit2}}\sim 1.5$, which
means that their limiting scenario is that their escaping matter flux is about
one-percent ($\sim 1\%$) level of their total energy flux.

So the rate of rise of temperature in stratosphere if without cooling is:

$\Delta T/\Delta t=(\Delta q/\Delta t)/(mC_{p})=0.01\cdot
f/(C_{p}\cdot\rho_{\text{tropopause}}\cdot\mathcal{H})=0.01\cdot f/(C_{p}\cdot
P_{\text{tropopause}}/g)$ (H61)

Taking a stratospheric temperature of $\sim 1000$K, $C_{p}=3.5$ R,
$P=0.1\text{bar}=10^{5}\text{dyn cm}^{-2}$, $f=100f_{\oplus}=1.37\cdot
10^{8}\text{erg/s cm}^{-2}$, $g=g_{\oplus}\approx 10^{3}\text{cm s}^{-2}$,
R$=8.314\cdot 10^{7}\text{erg/K mol}^{-1}$, so, $\Delta T/\Delta t\approx
0.5\times 10^{-4}$K/s.

Given these rates of absorption and emission of heat above the tropopause, the
time required for the temperature to change by one-degree Kelvin would be
$1\text{K}/(\Delta T/\Delta t)\sim\Delta t\sim 10^{4}$ seconds. Thus, we say
the _radiative time constant_ in this region, $\tau_{\text{rad}}$, is $10^{4}$
s (3 hours), for the 100-mbar level in the stratosphere assumed.

Therefore, in order for the temperature to be controlled radiatively in this
layer, the time constant for turbulent mixing and convective overturn,
$\tau_{\text{dyn}}$, the _dynamic time constant_ , must be longer than
$\tau_{\text{rad}}$.

At the tropopause, a typical sample of a cloudtop material may travel a
significant portion of the circumference of the planet $\sim 10^{4}$km before
encountering a downdraft, which requires a horizon mean wind velocity on order
of $10^{4}\text{km}/10^{4}\text{s}\sim 1$ km/s, which is comparable to the
sound speed $c_{s}$ in a solar-composition gas in $500\sim 2000$ Kelvin
temperature range.

$c_{s}=\sqrt{\left(\frac{\partial
P}{\partial\ln{\rho}}\right)_{\text{adiabatic}}}=\sqrt{\frac{\gamma
RT}{\mu}}=\sqrt{\frac{C_{p}RT}{\mu(C_{p}-R)}}\sim 1.5-3\text{~{}km/s}$ (H62)

Therefore, we expect to encounter sub-sonic (Mach$<1$), trans-sonic (Mach$\sim
1$), and then super-sonic (Mach$>1$), and then even hyper-sonic (Mach$\gg 1$)
flow regimes with a _shock-wave front_ present in certain part of the
exoplanetary atmosphere, which are generated and powered by the strong host-
stellar irradiation absorbed in the exoplanetary atmosphere.

This scenario should be explored with detailed hydrodynamic calculation in
future. Indeed this is a planetary-scale jet engine and then a rocket engine!
We expect that as the irradiation power is increased, the atmosphere is trying
to catch up with higher speed wind to carry away absorbed energy from sub-
stellar host spot. However, at a critical point, the flow can no longer catch
up and a _shock front_ emerges, and this _shock front_ is understood as a
_phase transition_ of the flow structure which accelerates the exiting outflow
from the ordinary thermal escape scenario into strong hydrodynamic jet with
very high-speed shooting towards zenith.

Our physical intuition is that for small exoplanets, the efficiency of this
_planetary-scale rocket engine_ is on the order of one-percent ($1\%$) in
terms of converting the absorbed host-stellar flux into the jet outflow. So
this shock front is expected to occur right around _tropopause_.

On the other hand, for large exoplanets the efficiency of this _planetary-
scale rocket engine_ is approaching the order of unity in terms of converting
the absorbed host-stellar flux into super-sonic outflow accelerated by the
_shock wave front_. This shock front is expected to occur then near the level
of _photosphere_.

###  H.2 Mass Loss Rate Estimate under non-limiting scenario

Turning back to the mass loss rate equation:

$\dot{m}=\bigg{(}1.66\times
10^{2}\text{~{}g}\cdot\text{cm}^{-2}\text{s}^{-1}\bigg{)}\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3}}\right)}\cdot\frac{R_{\oplus}^{2}}{M_{\oplus}}\cdot\eta\cdot
r^{2}$ (H63)

plug in the number for $M_{\oplus}$ and $R_{\oplus}$, we have,

$\dot{m}=\bigg{(}1.66\times
10^{2}\text{~{}g}\cdot\text{cm}^{-2}\text{s}^{-1}\bigg{)}\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3}}\right)}\cdot\frac{(6.371\times
10^{8}\text{cm})^{2}}{6\times 10^{27}\text{g}}\cdot\eta\cdot r^{2}$ (H64)

so that,

$\dot{m}=\bigg{(}1.1\times
10^{-8}\text{~{}s}^{-1}\bigg{)}\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3}}\right)}\cdot\eta\cdot
r^{2}$ (H65)

$\dot{m}=\left(\frac{1}{3.0\text{~{}year}}\right)\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3}}\right)}\cdot\eta\cdot
r^{2}$ (H66)

$\frac{\dot{m}}{m}\approx\left(\frac{1}{3.0\text{~{}year}}\right)\cdot\exp{\left(-\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3}}\right)}\cdot\frac{\eta\cdot
r^{2}}{m}$ (H67)

now the characteristic lifetime $\tau$ of such a planet can be estimated as
$m/\dot{m}$:

$\tau\sim\frac{m}{\dot{m}}\sim
3\text{~{}yr}\cdot\frac{1}{\eta}\cdot\exp{\left(\frac{a_{3}/2+b_{3}\cdot(m/r)}{T_{3}}\right)}\cdot\frac{m}{r^{2}}$
(H68)

$\tau\sim\frac{m}{\dot{m}}\sim 3\text{~{}yr}\cdot\frac{1}{\eta}\cdot
10~{}\widehat{}\left[\frac{a_{3}/2+b_{3}\cdot(m/r)}{2.303\cdot
T_{3}}\right]\cdot\frac{m}{r^{2}}$ (H69)

$\tau\sim\frac{m}{\dot{m}}\sim 3\text{~{}yr}\cdot\frac{1}{\eta}\cdot
10~{}\widehat{}\left[\frac{27.+7.6\cdot(m/r)}{2.303\cdot
T_{3}}\right]\cdot\frac{m}{r^{2}}$ (H70)

Therefore, the escape process and the lifetime of such a planet is highly
sensitive to the planet temperature.

##  Appendix I On Puffy and TTV Planets

Figure 12 shows an additional subplot of planet radius versus stellar
insolation flux. A critical cut-off flux of about 200 times the equivalent of
Earth’s insolation flux ($f_{\text{crit}}/f_{\oplus}\approx 200$) is clearly
identified in this subplot. Beyond this cut-off flux the planets bear no
H2-envelope. To make a comparison, this cut-off flux level is equivalent to
moving current Earth 15 times closer to the Sun so there is $15^{2}=225$ times
increase in flux. The overall picture of the flux-dependent mass-radius
diagram is consistent with a flux-driven H2-envelope escaping scenario.

One multi-planet system (HIP41378) of relatively low insolation flux is shown
for comparison. Many planets there bear significant H2-envelope and are
consistent with our hypothesis (Santerne et al., 2019) (Vanderburg et al.,
2016).

Figure 12: CHILL#1 for Small Exoplanets: Super-Earths/Mini-Neptunes/Cores with
H2-Envelopes. Log-Log Plot in Planet Masses ($M_{p}$) and Radii ($R_{p}$),
Zoom-out. Here we focus on the higher-flux portion of the planet radius versus
flux subplot. A cut-off flux of ($f_{\text{crit}}/f_{\oplus}\approx 200$) is
identified in the $(f/f_{\oplus})-R_{p}$ subplot, as measured in the unit of
solar constant ($f_{\oplus}=1.36$ kW/m2). Below this cut-off flux
($f_{\text{crit}}$), some planets can retain a low-mean-molecular-weight
H2/He-Envelope, which are consistent with their locations in the planet mass-
radius main diagram. Above this cut-off flux level ($f_{\text{crit}}$), no
planets bigger than $\sim 3R_{\oplus}$ exist. This suggests a likely loss of
such low-mean-molecular-weight H2/He-Envelope. An example of a multi-planet
system: HIP41378 Planetary System is shown here for comparison (Santerne et
al., 2019). Several planets in this system have received very low flux
compared to many other exoplanets shown in this figure, and thus, are expected
to retain their envelopes. The equilibrium temperatures of Planet d, e, and f
actually lie below the lower limit of the temperature-scale, and thus, they
show up as purple altogether. Planet b, c, and d seem to lie along the CHILL#1
critical for envelope loss.