#  Experimental Validation of Collision-Radiation Dataset for Molecular
Hydrogen in Plasmas

Keisuke Fujii<EMAIL_ADDRESS>Fusion Energy Division, Oak Ridge National
Laboratory, Oak Ridge, TN 37831-6305, United States of America  Keiji Sawada
Faculty of Engineering, Shinshu University, 4-17-1 Wakasato, Nagano 380-8553,
Japan  Kuzmin Arseniy  Department of Mechanical Engineering and Science,
Graduate School of Engineering, Kyoto University Kyoto 615-8540, Japan
Motoshi Goto  Masahiro Kobayashi  National Institute for Fusion Science, Toki,
Gifu, 5909-5292, Japan  Liam H. Scarlett  Dmitry V. Fursa  Igor Bray
Department of Physics and Astronomy, Curtin University, Perth, Western
Australia 6102, Australia  Mark C. Zammit  Theoretical Division, Los Alamos
National Laboratory, Los Alamos, NM 87545, United States of America  Theodore
M. Biewer  Fusion Energy Division, Oak Ridge National Laboratory, Oak Ridge,
TN 37831-6305, United States of America

###### Abstract

Quantitative spectroscopy of molecular hydrogen has generated substantial
demand, leading to the accumulation of diverse elementary-process data
encompassing radiative transitions, electron-impact transitions,
predissociations, and quenching. However, their rates currently available are
still sparse and there are inconsistencies among those proposed by different
authors. In this study, we demonstrate an experimental validation of such
molecular dataset by composing a collisional-radiative model (CRM) for
molecular hydrogen and comparing experimentally-obtained vibronic populations
across multiple levels.

From the population kinetics of molecular hydrogen, the importance of each
elementary process in various parameter space is studied. In low-density
plasmas (electron density $n_{\mathrm{e}}\lesssim 10^{17}\;\mathrm{m^{-3}}$)
the excitation rates from the ground states and radiative decay rates, both of
which have been reported previously, determines the excited state population.
The inconsistency in the excitation rates affects the population distribution
the most significantly in this parameter space. On the other hand, in higher
density plasmas ($n_{\mathrm{e}}\gtrsim 10^{18}\;\mathrm{m^{-3}}$), the
excitation rates from excited states become important, which have never been
reported in the literature, and may need to be approximated in some way.

In order to validate these molecular datasets and approximated rates, we
carried out experimental observations for two different hydrogen plasmas; a
low-density radio-frequency (RF) heated plasma ($n_{\mathrm{e}}\approx
10^{16}\;\mathrm{m^{-3}}$) and the Large Helical Device (LHD) divertor plasma
($n_{\mathrm{e}}\gtrsim 10^{18}\;\mathrm{m^{-3}}$). The visible emission lines
from $\mathrm{EF^{1}\Sigma_{g}^{+}}$, $\mathrm{H\bar{H}^{1}\Sigma_{g}^{+}}$,
$\mathrm{D^{1}\Pi_{u}^{\pm}}$, $\mathrm{GK^{1}\Sigma_{g}^{+}}$,
$\mathrm{I^{1}\Pi_{g}^{\pm}}$, $\mathrm{J^{1}\Delta_{g}^{\pm}}$,
$\mathrm{h^{3}\Sigma_{g}^{+}}$, $\mathrm{e^{3}\Sigma_{u}^{+}}$,
$\mathrm{d^{3}\Pi_{u}^{\pm}}$, $\mathrm{g^{3}\Sigma_{g}^{+}}$,
$\mathrm{i^{3}\Pi_{g}^{\pm}}$, and $\mathrm{j^{3}\Delta_{g}^{\pm}}$ states
were observed simultaneously and their population distributions were obtained
from their intensities. We compared the observed population distributions with
the CRM prediction, in particular the CRM with the rates compiled by Janev et
al., Miles et al., and and those calculated with the molecular convergent
close-coupling (MCCC) method. The MCCC prediction gives the best agreement
with the experiment, particularly for the emission from the low-density
plasma. On the other hand, the population distribution in the LHD divertor
shows a worse agreement with the CRM than those from low-density plasma,
indicating the necessity of the precise excitation rates from excited states.
We also found that the rates for the electron-attachment is inconsistent with
experimental results. This requires further investigation.

##  I Introduction

Molecular hydrogen emission appears in various low-temperature plasmas, such
as interstellar media Rosenthal  _et al._ (2000), process plasmas Lieberman
and Lichtenberg (2005), and cold-temperature regions of fusion plasmas Janev
_et al._ (2003); Fantz  _et al._ (2001); Ishihara  _et al._ . In particular,
in fusion devices, hydrogen (or its isotopologues) molecules play an important
role in the particle balance chain Stangeby and Others (2000), as well as one
of efficient heat exhaust by molecular assisted processes, such as molecular-
assisted recombination Janev  _et al._ (2003); Ohno  _et al._ (1998); Verhaegh
_et al._ (2022). Despite its importance, its quantitative diagnostics have
been rather unestablished.

A key technique to bridge the experimentally observable emission lines and
particle (ions, atoms, and molecules) behavior in plasmas is collisional-
radiative model (CRM), which solves the excited-state population based on the
rate equations. Since CRMs require a complete dataset of elementary process
rates, such as the electron-impact transition and spontaneous decay, a huge
amount of effort has been paid to measure, calculate, compile, and validate
these rates. Thanks to the efforts, the recommended rates for atomic species
have been compiled in a database, e.g., ADAS Summers , and have been in active
use for diagnostic purposes.

However, the situation for molecules is behind that for atomic species. The
essential difficulty of molecular spectroscopy is in the complex quantum
structure of molecules; the interplay between electron and nuclei motion makes
the energy structure more complex. This complexity leads to a higher hurdle in
composing the elementary-process rates, both experimentally and theoretically.
Hydrogen molecule, which is the simplest neutral molecule, has been studied
for a long time, and many rates have been reported, but there is still
inconsistency between the results among authors (see Sec. II.3 later). CRMs
for hydrogen molecules have also been developed by several authors based on
these rates Sawada and Fujimoto (1995); Greenland (2002); Lavrov  _et al._
(2006); Guzmán  _et al._ (2013); Shakhatov  _et al._ (2016); Sawada and Goto
(2016); Wünderlich  _et al._ (2021); yac , however, their validity remains
unclear and an experimental evidence to solve the inconsistencies are still
lacking.

The purpose of this paper is to compose a CRM for molecular hydrogen based on
the currently-available datesets and to validate these dataset by comparing
our CRM predictions with experimental observations.

In the next section, we will present the basic principle of the CRM and the
important elementary processes for hydrogen molecule. Inconsistencies among
the datasets reported by different authors, such as those proposed by Miles
Miles  _et al._ (1972), Janev Janev  _et al._ (2003), and those theoretically
calculated with the recently-developed molecular convergent close-coupling
(MCCC) method mcc ; Scarlett  _et al._ (2017, 2021a, 2021b, 2021c, 2021d,
2021d, 2022, 2022, 2023) will be pointed out. Based on our CRM, in Sec. III,
we will discuss the population kinetics and its dependence on the plasma
parameters. In particular, we will show that the population kinetics is
essentially different in low-density ($n_{\mathrm{e}}\lesssim
10^{17}\;\mathrm{m^{-3}}$, coronal phase) and high-density regions
($n_{\mathrm{e}}\gtrsim 10^{18}\;\mathrm{m^{-3}}$, saturation phase). The
important elementary processes are also different in these phases. Finally, in
Sec. IV we demonstrate the experimental validation of our CRM and point out
the current limitation in molecular-hydrogen spectroscopy.

##  II CRM for Molecular Hydrogen and the Elementary Processes

A CRM is essentially a set of rate equations for the excited-state population
with various excitation / deexcitation processes taken into account. For
molecular hydrogen, the important processes are radiative decay (Sec. II.2),
electron-impact excitation, deexcitation, ionization (Sec. II.3),
predissociation (Sec. II.4), quenching (Sec. II.5), and electron attachment
(Sec. II.6). These processes are schematically illustrated in Fig. 1.

Figure 1: A schematic illustration of the elementary processes of molecules in
plasmas.

By including these processes, the rate equation for the population in the
excited state $p$ can be written as follows:

$\displaystyle\frac{d}{dt}n_{p}$ $\displaystyle=\sum_{q}[C_{p\leftarrow
q}n_{e}+A_{p\leftarrow q}]n_{q}$ (1) $\displaystyle-\sum_{q}[C_{q\leftarrow
p}n_{e}+A_{q\leftarrow p}]n_{p}$ (2) $\displaystyle-\sum_{k}[C_{k\leftarrow
p}n_{e}+C_{k\leftarrow p}^{\mathrm{attach}}n_{e}+$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;A_{k\leftarrow p}+P_{k\leftarrow
p}+Q_{k\leftarrow p}n_{\mathrm{H_{2}}}]n_{p},$ (3)

where the first term on the right hand side (Eq. (1)) is the population influx
from $q$ state to $p$ state, while the second (Eq. (2)) and third terms (Eq.
(3)) are the population outflux from $p$ state to $q$ and $k$ state,
respectively. Here, $q$ indicates stable states of molecular hydrogen, and $k$
indicates other states, such as dissociative unstable states and bound states
of molecular ions. For the sake of simplicity, we assume the ionizing plasma,
so that the recombination is not important. This means that there is no influx
from the dissociative and ionized states. This assumption may be valid since
the volume association and recombination are often negligible in the density
range considered here ($n_{\mathrm{e}}$ $\lesssim 10^{20}\mathrm{\;m^{-3}}$),
since the dominant generation process of molecules is the surface-assisted
association.

$A_{p\leftarrow q}$ indicates the radiative decay rate, while $P_{k\leftarrow
p}$ is the predissociation rate, both of which are spontaneous processes.
$C_{p\leftarrow q}$ indicates the rate coefficient by electron impact
(excitation, deexcitation, and ionization). $C_{k\leftarrow
q}^{\mathrm{attach}}$ represents the rate coefficient of the electron-
attachment process

$\displaystyle\mathrm{H}_{2}^{\star}+e^{-}\rightarrow\mathrm{H}^{-}+\mathrm{H}$
(4)

in which the electron will attach to an excited molecule eventually leading to
the dissociation.

$Q_{k\leftarrow p}$ is the rate coefficient for the quenching process,

$\displaystyle\mathrm{H}_{2}^{\star}+\mathrm{H}_{2}\rightarrow
2\mathrm{H}+\mathrm{H}_{2},$ (5)

which also leads a dissociation eventually.

By assuming the quasi steady state for all the excited-state population
$dn_{p}/dt=0$, we obtain the population density

$\displaystyle n_{p}=R_{1}^{p}(n_{e},T_{e})n_{e}n_{1},$ (6)

where $n_{1}$ is the ground state density, and $R_{1}^{p}(n_{e},T_{e})$ is
called collisional-radiative population coefficient.

###  II.1 Quantum structure of hydrogen molecule

Figure 2: (a) A simplified level diagram of molecular hydrogen. (b), (c)
Potential energy curves for the hydrogen molecules. The vibrational levels are
indicated by horizontal bars for $\mathrm{X^{1}\Sigma_{g}^{+}}$ states.

The energy structure of molecules is dominantly represented by three types of
motions: electron motion, nuclear vibration, and nuclear rotation. Since the
time scale of the electron motion ($\approx 10^{-15}$ s) is faster than the
nuclear vibration ($\approx 10^{-13}$ s) and nuclear rotation ($\approx
10^{-12}$ s), the nuclei rarely move while an electron rotates around them. In
the Born-Oppenheimer approximation, we simply solve the Schöldinger equation
only for electrons with fixed nuclear positions, based on the time-scale
difference of their motion.

Figure 2 (a) shows a level diagram of the electronic structure of molecular
hydrogen. On one hand, since two electrons are orbiting around the nuclei in a
hydrogen molecule, the electron energy structure has some similarity with the
energy structure of atomic helium; there are two energy manifolds depending on
the relative directions of two electron spins, i.e., singlet ($S=0$) and
triplet ($S=1$) states. Here, $S$ indicates the spin quantum number of the
electrons. On the other hand, the axially-symmetric potential of a diatomic
molecule violates the conservation of the angular momentum of the electrons,
but instead its amplitude of the projection to the molecular axis is
conserved. A quantum number $\Lambda=|M_{L}|$ is assigned to this projection
and the states with different values of $\Lambda$ have different energies,
where $L$ is the orbital quantum number for electron, and $M_{L}$ is the
projection to the molecular axis. Capital greek letters
$\Sigma,\Pi,\Delta,\ldots$ are assigned to $\Lambda=0,1,2,\ldots$ states,
respectively (see the horizontal axis in Fig. 2 (a)). For the states with
$\Lambda>0$, there are twofold degeneracies, corresponding to $M_{L}=-\Lambda$
and $M_{L}=+\Lambda$ states. This degeneracy is related to the symmetry
property of the electron wavefunction. The electron wavefunction must be
either symmetric or anti-symmetric against the reflection at any plane passing
through both nuclei. If the sign of the wavefunction remains unchanged by this
reflection, a superscript $+$ is assigned to this state, while for the other
case a superscript $-$ is assigned.

For the diatomic molecules with the same-charge nuclei (e.g., $\mathrm{H_{2}}$
and HD), the electron wavefunction has another symmetry property, where the
wavefunction is either symmetric or anti-symmetric against the reflection
around the center of the two nuclei. If the sign of the wavefunction remains
unchanged by this reflection $g$ is assigned to this state, for while for the
other case $u$ is assigned (from the German gerade and ungerade).

Conventionally, a capital letter X is assigned to the electronic ground state,
and B, C, $\dots$ are used for the 1st, 2nd, $\dots$ electronic excited states
in the same multiplet to the X state, while lower case letters a, b, $\dots$
are used for the other multiplet states [ for historical reasons, there are
some irregularities in the labelling of hydrogen molecule states in these
conventional names, such as the EF state and B’ state, as shown in Fig. 2
(b)]. Each electronic state has been specified by the notation
$\mathrm{[assigned\;letter]}\;^{(2S+1)}\Lambda_{\mathrm{[g\;or\;u]}}^{+/-}$.
To avoid the confusion, in this paper, we use full notation for each
electronic state, e.g., $\mathrm{EF^{1}\Sigma_{g}^{+}}$ and
$\mathrm{d^{3}\Pi_{u}^{-}}$.

The Born-Oppenheimer approximation gives the excited energy of electrons with
a fixed nuclear position. Its dependence on the inter-nuclear distance can be
interpreted as the potential surfaces for nuclear motion. Figure 2 (b) and (c)
show the potential energy curves as a function of the internuclear distance
$R$ for the singlet and triplet states, respectively Nakashima and Nakatsuji
(2018); Wolniewicz (1993). $\mathrm{b^{3}\Sigma_{u}^{+}}$ state in the triplet
does not have a local minimum, meaning a purely dissociative state. The
vibrational and rotational motions of nuclei are obtained by solving the
Schrödinger equation for nuclear motions with these potential energy curves.
In Fig. 2 (b), the vibrational energy levels for
$\mathrm{X^{1}\Sigma_{g}^{+}}$ are indicated by horizontal lines. The
vibrational and rotational energy intervals for hydrogen molecules are
$\approx 0.5$ eV and $\approx$ 0.05 eV, respectively. In our CRM, only the
electronic and vibrational states will be considered, and the rotational
levels are ignored. On the other hand, the rotational states are needed to be
taken into account in the analysis of an observed emission spectrum.

###  II.2 Radiative transition

As similar to an atomic system, an excited molecule decays to a lower state by
emitting a photon. The selection rules for the electronic states is
$|\Lambda^{\prime}-\Lambda^{\prime\prime}|\leq 1$,
$S^{\prime}-S^{\prime\prime}=0$, and $\mathrm{g\leftrightarrow u}$. Figure 3
(a) shows the allowed transitions among each state by markers. Note that
according to the convention, a single prime ′ and double prime ′′ are used to
indicate the initial and final states of each transition, respectively.

Figure 3: Transitions available in each dataset. (a) Radiative decay rates
($A$) by Fantz et al. Fantz and Wünderlich (2006). Filled markers indicate the
available data, while the open markers indicate the optically-allowed
transitions but the data is not available. (b) the Franck-Condon factors (FCF)
by Fantz et al. Fantz and Wünderlich (2006). Note that the Franck-Condon
factors are symmetric, i.e.,
$q_{\alpha^{\prime}v^{\prime}-\alpha^{\prime\prime}v^{\prime\prime}}=q_{\alpha^{\prime\prime}v^{\prime\prime}-\alpha^{\prime}v^{\prime}}$.
(c) Availability of the electron-impact excitation cross-sections. The circle
markers indicate the available data in MCCC database mcc , while plus markers
and dots indicate those available in Janev dataset Janev  _et al._ (2003) and
Miles Miles  _et al._ (1972), respectively.

The radiative decay rate among electronic and vibronic states,
$A^{\alpha^{\prime}v^{\prime}}_{\alpha^{\prime\prime}v^{\prime\prime}}$, have
been computed by Funtz et al Fantz and Wünderlich (2006) based on the
electronic dipole transition moment available in literature, not only for
$\mathrm{H_{2}}$, but also all the possible isotopologues. The transitions
with the data available are shown by filled markers in Fig. 3 (a).

Among allowed transitions among the listed levels, the one for
$\mathrm{c^{3}\Pi_{u}^{\pm}}$ $\leftarrow$$\mathrm{j^{3}\Delta_{g}^{\pm}}$ is
missing in their data set, because of the unavailability of the electronic
dipole transition moment. We approximate the transition rate based on the
experimental value by Astashkevich et al Astashkevich  _et al._ (1996), and
the method based on Franck-Condon factor and Hönl-London factor.

###  II.3 Electron-impact transition

Figure 4: A comparison of the excitation cross sections from
$\mathrm{X^{1}\Sigma_{g}^{+}}$ $(v=0)$ state to the $n=2$ and $n=3$ excited
states. The vibrational quantum number for final state are not distinguished
(the cross-sections are summed up). The cross-sections by MCCC mcc , Janev
Janev  _et al._ (2003), and Miles Miles  _et al._ (1972) are shown by
different lines.

Electron-impact transition is the most important process that generates
excited-state molecules in many plasmas. Several groups have proposed
different sets of cross-sections. Miles et al. have assumed an analytical form
for the vibrationally-resolved electron-impact cross sections based on the
generalized oscillator strength method, and adjust a few parameters so that
this form fits experimental data available at that time Miles  _et al._
(1972). Later, Janev et al. have also compiled various electron-impact cross
sections for molecular hydrogen Janev  _et al._ (2003). However, there has
been inconsistency among them. As seen in Fig. 4, a significant difference can
be seen in some of the cross sections, such as $\mathrm{X^{1}\Sigma_{g}^{+}}$
$\rightarrow$ $\mathrm{H\bar{H}^{1}\Sigma_{g}^{+}}$.

Also recently, a systematic theoretical method based on MCCC method has been
developed which is able to calculate the vibrational-state-resolved electron-
impact transition cross-sections mcc ; Scarlett  _et al._ (2017, 2021a, 2021b,
2021c, 2021d, 2021d, 2022, 2022, 2023). The MCCC rates are also shown in Fig.
4. These three still show inconsistencies, although the experimental
comparison for the $\mathrm{d^{3}\Pi_{u}^{-}}$ state population by Wünderlich
et al. suggests a better performance of the MCCC cross sections Wünderlich
_et al._ (2021).

MCCC also provides some of the cross sections among excited states. The
available rates are shown in Fig. 3 (c). However, many cross sections from
excited states are still missing, which are important for high-density plasmas
as described below. The following two approximation methods for such missing
rates can be considered.

####  II.3.1 Helium approximation

In the first CR model for hydrogen molecule Sawada and Fujimoto (1995), the
authors have used the cross-section for the transition among helium excited
states. For example, $\mathrm{EF^{1}\Sigma_{g}^{+}}$ has $2s$ electron orbit,
while $\mathrm{D^{1}\Pi_{u}^{-}}$ has $3p$ orbit. Thus, this rate coefficient
may be approximated as follows,

$\displaystyle
C_{\mathrm{D^{-}}(v^{\prime\prime})\leftarrow\mathrm{EF}(v^{\prime})}\approx\frac{1}{3}C^{\mathrm{He}}_{3^{1}\mathrm{P}\leftarrow
2^{1}\mathrm{S}}\;q_{\mathrm{EF}(v^{\prime})-\mathrm{D}(v^{\prime\prime})},$
(7)

where $C^{\mathrm{He}}_{3^{1}\mathrm{P}\leftarrow 2^{1}\mathrm{S}}$ is the
excitation rate of helium atom for its $3^{1}\mathrm{P}\leftarrow
2^{1}\mathrm{S}$ transition. The factor $1/3$ accounts for the three
$3^{1}\mathrm{P}$ states in hydrogen molecule (i.e., $\mathrm{B^{\prime
1}\Sigma_{u}^{+}}$, $\mathrm{D^{1}\Pi_{u}^{+}}$, $\mathrm{D^{1}\Pi_{u}^{-}}$).
$q_{\mathrm{EF}(v^{\prime})-\mathrm{D}(v^{\prime\prime})}$ is the Franck-
Condon factor between $\mathrm{EF^{1}\Sigma_{g}^{+}}$ and
$\mathrm{D^{1}\Pi_{u}^{\pm}}$ states.

In order to see the accuracy of this approximation, we show in Fig. 5 the rate
coefficients from each $n=2$ states ($\mathrm{EF^{1}\Sigma_{g}^{+}}$,
$\mathrm{B^{1}\Sigma_{u}^{+}}$, $\mathrm{C^{1}\Pi_{u}^{\pm}}$,
$\mathrm{a^{3}\Sigma_{g}^{+}}$, $\mathrm{c^{3}\Pi_{u}^{\pm}}$) to any of $n=3$
states. The bold blue curves show the calculation result by MCCC, while the
orange chain curves show those for helium atoms. These two sets of cross
sections agree within a factor of $\approx 10$.

Figure 5: The total rate coefficients from $n=2$ states
($\mathrm{EF^{1}\Sigma_{g}^{+}}$, $\mathrm{B^{1}\Sigma_{u}^{+}}$,
$\mathrm{C^{1}\Pi_{u}^{\pm}}$, $\mathrm{a^{3}\Sigma_{g}^{+}}$,
$\mathrm{c^{3}\Pi_{u}^{\pm}}$) to all of the $n=3$ states. The bold curves are
from the MCCC database mcc , while the helium approximation and Mewe’s
approximation are shown by chain curves and thin curves, respectively. See the
text for details.

####  II.3.2 Mewe’s approximation

Another possible approximation would be the one proposed by Mewe Mewe , based
on the line-strength for dipole transitions. With this approximation, the rate
coefficient for the electron-impact excitation from $q$ state to $p$ state can
be written as

$\displaystyle C_{p\leftarrow
q}(T_{e})\approx\frac{1}{g_{q}}\frac{\beta}{\sqrt{kT_{e}}}S_{qp}\exp\left[-\frac{\omega_{pq}}{kT_{e}}\right],$
(8)

with

$\displaystyle\beta=\frac{2^{7/2}}{3\sqrt{3}}\pi^{3/2}\alpha
a_{0}^{2}c\sqrt{E_{\mathrm{H}}}\frac{1}{e^{2}a_{0}^{2}}\xi$ (9)

where $g_{p}$ is the statistical weight of the initial state, $k$ is the
Boltzmann constant, $T_{e}$ is the electron temperature, and
$\omega_{pq}=E_{q}-E_{p}$ is the energy difference between states $p$ and $q$,
$E_{\mathrm{H}}\approx$ 27.2 eV is the Hartree energy, $\alpha$ is the fine
structure constant, $e$ is the elementary charge, $a_{0}$ is the Bohr radius,
and $c$ is the light speed. $S_{qp}$ is the line strength of the transition,
which is evaluated from the A coefficients for the transition,

$A_{p\leftarrow q}=\gamma\frac{1}{g_{q}}\omega_{pq}^{3}S_{pq}.$ (10)

with

$\gamma=\frac{4}{3}\alpha^{4}\frac{c}{a_{0}}\frac{1}{E_{\mathrm{H}}^{3}}\frac{1}{e^{2}a_{0}^{2}}$
(11)

$\xi$ is the gaunt factor, which is represented by the following equation

$\displaystyle\xi=a+b\exp\left(\frac{\omega_{pq}}{kT_{e}}\right)\exp_{1}\left(\frac{\omega_{pq}}{kT_{e}}\right),$
(12)

where $a\approx 0$ and $b\approx 1$ are parameters that depend on the
particular system. We assume $a=-0.03$ and $b=0.28$ so that Eq. (8) best
represents the MCCC rate coefficients. Orange curves in Fig. 5 show the rate
coefficients by Eq. (8). Equation (8) (orange thin curve) represents the MCCC
rate coefficients better than the corresponding helium rates (green chain
curve). This is partly because of the degrees of freedom in $\xi$ (we
basically fit by adjusting $a$ and $b$) but also because it takes into account
the energy difference of the actual system.

###  II.4 Predissociation

Predissociation is a spontaneous process where an excited molecule undergoes
an internal conversion to another state, typically leading to a dissociation.
One of the important dissociation paths is

$\displaystyle\mathrm{H}_{2}(c^{3}\Pi^{+})\rightarrow\mathrm{H}_{2}(b^{3}\Sigma^{+})\rightarrow
2\mathrm{H},$ (13)

the decay rate of which is $1.6\times 10^{8}\;\mathrm{s}^{-1}$ and $7.0\times
10^{8}\;\mathrm{s}^{-1}$ for $v=0,N=1$ and $v=1,N=1$ states, respectively
Comtet and De Bruijn (1985). Although $\mathrm{c^{3}\Pi_{u}^{+}}$ is a
radiatively metastable state, this predissociation process reduces its
lifetime significantly.

The lifetime of an isolated excited molecules $\tau$ can be written as

$\displaystyle\frac{1}{\tau_{p}}=\sum_{q}A_{q\leftarrow
p}+\sum_{q}P_{q\leftarrow p}.$ (14)

The values of $\tau_{p}$ have been compiled by Astashkevich et al Astashkevich
and Lavrov (2015). We evaluate $\sum_{q}P_{q\leftarrow p}$ for each excited
states by the table in their paper. Note that although the predissociation
rate weakly depends on the rotational quantum states, we use the value for the
smallest $N$ state for our CR model.

###  II.5 Quenching rates

Quenching is a similar process to the predissociation, but induced by a heavy-
particle collision (Eq. (5)). Wedding et al has reported the rate coefficient
$(1.88\pm 0.10)\times 10^{-15}\mathrm{\;m^{3}s^{-1}}$ for 300-K
$\mathrm{H_{2}}$ collisions Wedding and Phelps (1988). As this number is
$\approx 10^{2}$ times smaller than the electron-impact excitation rate (Fig.
5), this process becomes important for low-ionization-degree plasmas,
$n_{\mathrm{e}}/n_{\mathrm{H_{2}}}\lesssim 10^{-3}$. In our experimental
conditions described below, this process is not significant.

###  II.6 Dissociative attachment

In Ref. Datskos  _et al._ (1997), it has been suggested that the dissociative
attachment (Eq. (4)) is very efficient for high-Rydberg state of molecualr
hydrogen, with the rate coefficient of $C^{\mathrm{attach}}\approx 6\times
10^{-11}\mathrm{\;m^{3}/s}$. However, the original paper Datskos  _et al._
(1997) has not specified which states are regarded as a high-Rydberg state. In
Refs. Wünderlich  _et al._ (2020); yac , the authors have assumed this value
for $n\geq 3$ states. In our paper, we also examine the validity of this
assumption.

##  III CRM predictions and population dynamics

Our CRM described in this paper considers the population of molecule in each
electronic- and vibrational-states up to $n=4$ states, while the rotational
states are not resolved. Also, we assume that the population in the
electronic-ground state can be represented by the Boltzmann distribution,

$\displaystyle
n_{\mathrm{X}}(v)=n_{\mathrm{H_{2}}}\frac{1}{Z}\exp\left(-\frac{E_{v}}{kT_{v}}\right),$
(15)

where $\mathrm{H_{2}}$ is the molecular density in
$\mathrm{X^{1}\Sigma_{g}^{+}}$ state, $E_{v}$ is the energy of vibrational
state $v$ in $\mathrm{X^{1}\Sigma_{g}^{+}}$, and $T_{v}$ is the vibrational
temperature. $k$ is the Boltzmann constant. Here, $Z$ is the normalization
constant,

$\displaystyle Z=\sum_{v}\exp\left(-\frac{E_{v}}{kT_{v}}\right).$ (16)

The solution of the CRM can be represented by a similar manner to Eq. (6), but

$\displaystyle n_{p}=R_{1}^{p}(n_{e},T_{e},T_{v})n_{e}n_{\mathrm{H_{2}}}.$
(17)

Figure 6: (Top panel) Population coefficients $R_{1}^{p}$ for some excited
states of molecular hydrogen calculated with our CRM, and its
$n_{\mathrm{e}}$-dependence. $T_{\mathrm{e}}=10$ eV and $T_{v}=0.1$ eV are
assumed. $R_{1}$ for $\mathrm{EF^{1}\Sigma_{g}^{+}}$ has a significant
$n_{\mathrm{e}}$-dependence, while that for $\mathrm{D^{1}\Pi_{u}^{\pm}}$ is
small. (2nd, 3rd, and 4th rows): $n_{\mathrm{e}}$-dependence for $n=2$ and
$n=3$ states. The population for the vibrational states are summed up. The
thick solid curves show the population coefficients calculated with MCCC
cross-sections and Mewe’s approximations for the missing transitions, while
thin solid curves show the population coefficients calculated with helium
cross-sections for missing transitions. These two shows a significant
difference in the high-$n_{\mathrm{e}}$ space, which is the systematic
uncertainty in the current CRM. Thin dotted curve shows the result with the
dissociative attachment included. Since the rate coefficient for the
dissociative attachment recommended in Ref. yac  is large compared with the
electron-impact rate, the inclusion of this process significantly changes the
population. In particular, the critical density between the coronal phase
(where $R_{1}^{p}\approx n_{\mathrm{e}}^{0}$) and the saturation phase (where
$R_{1}^{p}\approx n_{\mathrm{e}}^{-1}$) moves to the lower density side.

Figure 6 (a) shows examples of $R_{1}^{p}(n_{e},T_{e},T_{v})$ with
$T_{\mathrm{e}}=10\mathrm{\;eV}$ and $T_{\mathrm{v}}=0.1\mathrm{\;eV}$, for
several different values of $n_{\mathrm{e}}$. Horizontal position shows the
excited states of hydrogen molecule; each column separated by dotted lines
show an electronic state, while in each column the vibrational quantum number
dependence is shown. Note that for this calculation, we ignore the quenching
and electron-attachment processes. MCCC excitation rates are used and Mewe’s
approximation for the dipole transition where MCCC is unavailable.

Different $n_{\mathrm{e}}$-dependence can be seen for different electronic
states. For example, $R_{1}^{p}$ for $\mathrm{D^{1}\Pi_{u}^{\pm}}$ states
shows little $n_{\mathrm{e}}$-dependence, while that for
$\mathrm{EF^{1}\Sigma_{g}^{+}}$ state shows a significant $n_{\mathrm{e}}$
dependence. Such population dynamics can be understood from the dominant
elementary process for each state.

Figure 7 shows the elementary-process composition of influx and outflux. In
the low-$n_{\mathrm{e}}$ limit, the influx to each state is dominated by the
collisional excitation from the ground state (denoted by $C_{n=1}$), while the
dominant outflux is the radiative decay to the lower state (denoted by $A$ in
Fig. 7). Thus, the population can be written as

$\displaystyle n_{p}\approx\frac{\sum_{v}C_{p\leftarrow
v}n_{\mathrm{X}}(v)}{\sum_{q}A_{q\leftarrow p}}n_{\mathrm{e}},$ (18)

where $R_{1}^{p}\propto n_{\mathrm{e}}^{0}$. This phase is called coronal
phase, and the population coefficient has no $n_{\mathrm{e}}$-dependence.

As $n_{\mathrm{e}}$ increases, the dominant outflux changes to the collisional
transition to other states ($C_{n=2}$ and $C_{n=3}$ in Fig. 7). In this
region, the population may be written as follows;

$\displaystyle n_{p}\approx\frac{\sum_{v}C_{p\leftarrow
v}n_{\mathrm{X}}(v)n_{\mathrm{e}}}{\sum_{q}C_{q\leftarrow p}n_{\mathrm{e}}},$
(19)

where the dominant influx (numerator) and outflux (denominator) are both the
electron-impact excitation. Thus, the $n_{\mathrm{e}}$-dependences on the
numerator and denominator are canceled out and $R_{1}^{p}\propto
n_{\mathrm{e}}^{-1}$. This phase is called saturation phase.

Since $\mathrm{D^{1}\Pi_{u}^{\pm}}$ has shorter radiative lifetime ($\approx
3$ ns), the critical value of $n_{\mathrm{e}}$ between the coronal to
saturation phases is higher. At $n_{\mathrm{e}}$ $\approx
10^{19}\;\mathrm{m^{-3}}$, this state is still in coronal phase and thus the
population coefficient has no $n_{\mathrm{e}}$-dependence (see Fig. 6). On the
other hand, since $\mathrm{EF^{1}\Sigma_{g}^{+}}$ has longer radiative
lifetime ($\approx 100$ ns), this transition happens in a lower
$n_{\mathrm{e}}$, and above this critical density, $R_{1}^{p}\propto
n_{\mathrm{e}}^{-1}$.

Figure 7: Compositions of the population influx and outflux for some levels,
as a function of $n_{\mathrm{e}}$. $C_{X}$ indicates the influx from the
ground state, while $C_{n=2}$ and $C_{n=3}$ indicates flux from / to $n=2$ and
$n=3$ states, respectively. In the lower density region, the outflux from all
of these states is the radiative decay (indicated by $A$) to lower states,
while in higher density region collisional transition dominates the outflux.
The critical value of $n_{\mathrm{e}}$ is different for these states depending
on the relative importance of the radiative decay and the collisional
transition. All the values were calculated with our CRM under $T_{\mathrm{e}}$
= 10 eV and $T_{v}=0.1$ eV.

Lower panels in Fig. 6 shows the $n_{\mathrm{e}}$-dependence of the population
coefficient for several electronic states (vibrational populations are summed
up). In the low-density region $R_{1}^{p}\approx n_{\mathrm{e}}^{0}$ while in
high-density regions $R_{1}^{p}\approx n_{\mathrm{e}}^{-1}$, which is
consistent with the above discussion.

$R_{1}^{p}$ values for $\mathrm{g^{3}\Sigma_{g}^{+}}$ and
$\mathrm{i^{3}\Pi_{g}^{\pm}}$ levels show a positive
$n_{\mathrm{e}}$-dependence in a certain density range in contrast to those
for the other levels. This can be understood from the contribution of the two-
step excitation ($C_{n=2}$-influx in Fig. 7). Since $(n=2)$-state population
($n_{2}$) is proportional to $n_{\mathrm{e}}^{1}$ in the low-density region,
the two-step excitation flux is proportional to $n_{2}n_{\mathrm{e}}\propto
n_{\mathrm{e}}^{2}$. Because of the contribution of this process, $R_{1}^{p}$
values can have a positive $n_{\mathrm{e}}$-dependence. The population
converges to $n_{p}\propto n_{\mathrm{e}}^{0}$ at the high-density limit as
similar to the other states.

###  III.1 Dependence on different datasets

####  III.1.1 Difference due to the rate-approximation methods

Current database, including MCCC, do not provide the cross-sections from the
excited states with $n\geq 3$, and there are uncertainty which approximations
are more appropriate, i.e., Mewe’s approximation (Sec. II.3.2) or helium
approximation (Sec.II.3.1). Two solid curves in each lower panel of Fig. 6
show the population coefficients predicted by our CRM, but with (bold) Mewe’s
approximation and (thin) helium-crosssections for unavailable transitions in
MCCC database. As in the higher-$n_{\mathrm{e}}$ limit the influx and outflux
are dominated by the electron-impact transitions among excited states (see
Fig. 7 too), the difference in the values of $R_{1}^{p}$ becomes bigger. On
the other hand, in the lower-$n_{\mathrm{e}}$ limit the influx is the
excitation from the ground state (the cross sections of which are available in
the MCCC database mcc ) while the dominant outflux is the radiative decay.
Thus the uncertainty due to the cross-sections for excited states is smaller.
We will examine this effect later.

####  III.1.2 The effect of the dissociative attachment

Thin dotted curves in Fig. 6 shows the same calculation but with the
dissociative attachment process included. Since the proposed rate coefficients
for the dissociative attachment for $n=3$ states in Ref. yac  is much higher
than the electron-impact excitation rate (i.e., $\sum_{q}C_{q\rightarrow
p}\gg\sum_{k}C_{k\rightarrow p}^{\mathrm{attach}}$), all the $(n=3)$-states,
including $\mathrm{D^{1}\Pi_{u}^{\pm}}$, fall in the saturation phase at
$n_{\mathrm{e}}>10^{17}\mathrm{\;m^{-3}}$ if this process is included.

This effect may be examined by a $n_{\mathrm{e}}$-dependence of the population
ratio, e.g., that between $\mathrm{EF^{1}\Sigma_{g}^{+}}$ and
$\mathrm{D^{1}\Pi_{u}^{\pm}}$. If the proposed rate is correct and both the
states are in the saturation phase, the population ratio should have only a
weak $n_{\mathrm{e}}$-dependence. On the other hand, if the proposed rate is
overestimated, $\mathrm{D^{1}\Pi_{u}^{\pm}}$ is in coronal phase and thus the
population ratio should have a strong $n_{\mathrm{e}}$-dependence. This will
be discussed later in Sec. IV.2 and Fig. 16.

####  III.1.3 The difference between the excitation-rate datasets

In Fig. 8, we show the results with different sets of the excitation rates,
i.e., those by Janev (orange), Miles (green), and MCCC (blue), without
including the dissociative attachment process. Because of the difference in
the rate coefficients for electron-impact transitions, the population
coefficients are different. In higher-temperature plasma (Fig. 8 (a)), the
population by Janev and MCCC are similar, while that by Miles shows a
significant difference from the other two. On the other hand, in lower-
temperature plasma, the population in $\mathrm{EF^{1}\Sigma_{g}^{+}}$ show a
significant difference among these three datasets.

Figure 8: Population coefficients $R_{1}^{p}$ for some excited states of
molecular hydrogen calculated with our CRM, with different datasets for the
electron-impact transitions. Top and bottom panels show the prediction under
different values of $n_{\mathrm{e}}$ and $T_{\mathrm{e}}$. By different
colors, the results with different datasets (MCCC mcc , Janev Janev  _et al._
(2003), and Miles Miles  _et al._ (1972)) are shown.

##  IV Experimental comparisons

In order to compare with the CRM predictions described above, we carry out
experimental observations of emission spectra from molecular hydrogen. To
cover a wide range of the parameter space, we measure the emission from a low-
density RF plasma in Shinshu university (Sec. IV.1), and a divertor region of
Large Helical Device (LHD, Sec. IV.2). As described below, the RF plasma has
$n_{\mathrm{e}}\approx 10^{16}\mathrm{\;m^{-3}}$ and $T_{\mathrm{e}}\approx
10^{0}\mathrm{\;eV}$, while the LHD divertor has much higher temperature and
density range of $n_{\mathrm{e}}\approx 10^{19}\mathrm{\;m^{-3}}$ and
$T_{\mathrm{e}}\approx 10^{1}\mathrm{\;eV}$. Furthermore, the parameter in the
LHD divertor can be varied by changing the heating and fueling conditions.
This even widens the parameter space we consider here.

###  IV.1 Low-density plasma experiment in Shinshu university

We measured an emission spectrum from a low-temperature RF plasma at Shunshu
University, Japan. A hydrogen plasma was generated in a vacuum chamber made of
Pyrex glass, with 50-mm inner diameter and 1100-mm length. This chamber is
located inside two solenoids, which generates 0.012 T magnetic field at the
plasma center. A $\approx$100-W RF power is applied to generate the plasma.
This experimental setup is similar to that shown in Ref. Sawada  _et al._
(2010), however instead of the pure-helium plasma in the previous work, we
generated a plasma with helium-hydrogen gas mixture in this work. The partial
gas pressure of hydrogen molecule and helium atom are 0.07 and 0.02 torr,
respectively.

An echelle spectrometer with a crossed disperser configuration (EMP-200 Bunko-
Keiki, 202-mm focal length) and a charge-coupled-device (CCD) detector
(DV420A-OE, Andor inc., $255\times 1024$ pixels, $26\times 26\mathrm{\;\mu
m^{2}}$ pixel size) were used to measure the emission spectrum. This
spectroscopic system can simultaneously measure the spectrum in the wavelength
range of 400 – 800 nm with 0.05-nm resolution.

The values of $n_{\mathrm{e}}$ and $T_{\mathrm{e}}$ in the plasma were
estimated as $T_{\mathrm{e}}\approx 3\;\mathrm{eV}$ and $n_{\mathrm{e}}\approx
10^{16}\;\mathrm{m^{-3}}$, based on the helium-line-ratio method described in
Ref. Sawada  _et al._ (2010). In this method, the absolute emission
intensities from multiple excited states of helium atoms (
$3s\,^{1}\mathrm{S}$, $3p\,^{1}\mathrm{P}$, $3d\,^{1}\mathrm{D}$,
$4s\,^{1}\mathrm{S}$, $4p\,^{1}\mathrm{P}$, $4d\,^{1}\mathrm{D}$,
$5s\,^{1}\mathrm{S}$, $5d\,^{1}\mathrm{D}$, $6d\,^{1}\mathrm{D}$,
$3s\,^{3}\mathrm{S}$, $3p\,^{3}\mathrm{P}$, $3d\,^{3}\mathrm{D}$,
$4s\,^{3}\mathrm{S}$, $4d\,^{3}\mathrm{D}$, $5s\,^{3}\mathrm{S}$, and
$5d\,^{3}\mathrm{D}$) are fitted by a prediction of collisional-radiative
model for helium atoms. By the fit, several adjustable parameters are
determined, which are $n_{\mathrm{e}}$, $T_{\mathrm{e}}$, density of
metastable states ($2s^{1}\mathrm{S}$ and $2s^{3}\mathrm{S}$), and the
radiative reabsorption effect by $3p^{1}\mathrm{P}$ and $4p^{1}\mathrm{P}$
states. Note that although many adjustable parameters are used to fit multiple
emission intensities, the analysis is robust enough. The value of
$T_{\mathrm{e}}$ is basically determined by the population ratio among
$3s\,^{1}\mathrm{S}$, $4s\,^{1}\mathrm{S}$, $3s\,^{3}\mathrm{S}$ and
$3s\,^{3}\mathrm{P}$ states, where the other effects (excitation from
metastable and radiation reabsorption) are negligible in this condition. The
value of $n_{\mathrm{e}}$ is in turn determined by comparing the absolute
intensities with the helium atom density and the excitation rates, which are
the function of $T_{\mathrm{e}}$.

Figure 9: (top) Emission spectrum observed from the low-density RF plasma with
helium-hydrogen mixture. To clearly present the emission lines from molecular
hydrogen, strong atomic lines and helium lines are omitted. (bottom) The
composition of the hydrogen-molecule spectrum for each upper-electronic
states, which is obtained from the spectral fit.

Figure 9 shows the observed emission spectrum. To clearly present the
molecular emission, strong atomic helium and hydrogen emission lines as well
as continuum emission were excluded from the figure. Three clusters of lines
can be found in the visible range; 450–500 nm, 570–630 nm, 720–800 nm.

####  IV.1.1 Analysis method

The observed emission corresponds to transitions with the change in
electronic, vibrational and rotational states, while the CRM we have developed
in Sec. II resolves only the electronic and vibrational states. In order to
obtain the population from the experiment in each electronic and vibrational
state but integrated over the rotational states, we assume the Boltzmann
distribution for the rotational population,

$\displaystyle
n_{\alpha^{\prime}v^{\prime}N^{\prime}}\approx\frac{n_{\alpha^{\prime}v^{\prime}}}{Z_{\alpha^{\prime}v^{\prime}}}(2S^{\prime}+1)(2N^{\prime}+1)g^{\mathrm{as}}_{N^{\prime}}\exp\left(-\frac{E_{\alpha^{\prime}v^{\prime}}(N^{\prime})}{kT_{\alpha^{\prime}v^{\prime}}}\right),$
(20)

where $\alpha^{\prime}$, $v^{\prime}$, $S^{\prime}$, and $N^{\prime}$ indicate
the electronic state, vibrational quantum number, electron spin quantum
number, and rotational quantum number of the upper level, respectively, and
$E_{\alpha^{\prime}v^{\prime}}(N^{\prime})$ is the excitation energy.
$g^{\mathrm{as}}_{N^{\prime}}$ is the nuvlear spin statistical weight.
$Z_{\alpha^{\prime},v^{\prime}}$ is the partition function. We assume an
independent rotational temperature $T_{\alpha^{\prime},v^{\prime},N^{\prime}}$
for each electronic and vibrational state. The population and rotational
temperature in each electronic and vibronic state,
$n_{\alpha^{\prime},v^{\prime}}$, were optimized so that the predicted
emission intensity
$I^{\alpha^{\prime}v^{\prime}N^{\prime}}_{\alpha^{\prime\prime}v^{\prime\prime}N^{\prime\prime}}=A^{\alpha^{\prime}v^{\prime}N^{\prime}}_{\alpha^{\prime\prime}v^{\prime\prime}N^{\prime\prime}}n_{\alpha^{\prime}v^{\prime}N^{\prime}}$
represents the observed spectrum the best. We use the literature values for
the A coefficient if available, otherwise we use the Hönl-Londom approximation
assuming Hund’s case b with the vibrationally-resolved A coefficient values
Fantz and Wünderlich (2006).

We fit the entire spectrum based on the above assumption. The decomposition of
the spectrum is shown in the bottom part of Fig. 9. It is clear that the three
clusters have different upper-electronic states; 450 – 500 nm lines are mainly
originated from $\mathrm{GK^{1}\Sigma_{g}^{+}}$ and
$\mathrm{I^{1}\Pi_{g}^{\pm}}$ states, 570 – 630 nm lines are from
$\mathrm{d^{3}\Pi_{u}^{\pm}}$, $\mathrm{i^{3}\Pi_{g}^{\pm}}$, and
$\mathrm{j^{3}\Delta_{g}^{\pm}}$ states, and 720 – 800 nm lines are from
$\mathrm{EF^{1}\Sigma_{g}^{+}}$, $\mathrm{H\bar{H}^{1}\Sigma_{g}^{+}}$,
$\mathrm{D^{1}\Pi_{u}^{\pm}}$, $\mathrm{I^{1}\Pi_{g}^{\pm}}$, and
$\mathrm{J^{1}\Delta_{g}^{\pm}}$ states.

Figure 10: Distribution of the line-integrated population of molecular
hydrogen observed for the low-density RF plasma (markers). Kinked lines for
each electronic state are the prediction by the CRM with three different
excitation cross section (blue: MCCC, green: Miles, orange: Janev). The values
of $n_{\mathrm{e}}$ and $T_{\mathrm{e}}$ to calculate the population were
estimated based on the helium line-ratio method Sawada  _et al._ (2010), which
are $T_{\mathrm{e}}\approx 3\;\mathrm{eV}$ and $n_{\mathrm{e}}\approx
10^{16}\;\mathrm{m^{-3}}$. Equation (21) is used to compute $T_{v}\approx
2400$ K. The only one adjustable parameter is the line-integrated molecular
density, $n_{\mathrm{H_{2}}}L$, is adjusted so that the population in
$\mathrm{EF^{1}\Sigma_{g}^{+}}$ $(v=6)$ is identical between the observation
and the prediction, which are $n_{\mathrm{H_{2}}}L\approx
10^{20}\;\mathrm{m^{-2}}$.

Markers in Fig. 10 shows the population obtained from the analysis. We also
plot the prediction by our CRM with the values of $n_{\mathrm{e}}$ and
$T_{\mathrm{e}}$ obtained from the helium-line-ratio method. The value of
$T_{v}$, which is the vibrational temperature of the ground state molecules,
are obtained from the empirical relation Brezinsek  _et al._ (2002),

$\displaystyle\frac{T_{v}}{\mathrm{K}}\approx 2400+2.6\times
10^{-16}\left[\frac{n_{\mathrm{e}}}{\mathrm{m^{-3}}}\right].$ (21)

For this low-density plasmas, this means essentially $T_{v}\approx
2400\mathrm{\;K}$. Only one adjustable parameter in the comparison is the
value of $n_{\mathrm{H_{2}}}L$ (where $L$ is the effective diameter of the
plasma). We scale the prediction by adjusting this value so that the
population in $\mathrm{EF^{1}\Sigma_{g}^{+}}$ $v=6$ states becomes identical
between the experiment and the simulation (the square marker in Fig. 10). From
this scaling, $n_{\mathrm{H_{2}}}L\approx 10^{20}\;\mathrm{m^{-2}}$ is
obtained, which is consistent with the gas density ($\approx
10^{22}\;\mathrm{m^{-3}}$) and the plasma column diameter $L\approx
10^{-2}\;\mathrm{m}$.

The three types of lines in Fig. 10 show the CRM predictions with different
cross-section datasets, i.e., those by MCCC (blue), by Miles (green), and
Janev (orange). It is clearly seen that the prediction by the MCCC cross-
sections represents the experimental observation the best, in particular the
populations in $\mathrm{EF^{1}\Sigma_{g}^{+}}$,
$\mathrm{H\bar{H}^{1}\Sigma_{g}^{+}}$, $\mathrm{J^{1}\Delta_{g}^{\pm}}$,
$\mathrm{h^{3}\Sigma_{g}^{+}}$, $\mathrm{e^{3}\Sigma_{u}^{+}}$,
$\mathrm{d^{3}\Pi_{u}^{\pm}}$, and $\mathrm{j^{3}\Delta_{g}^{\pm}}$. Janev and
Miles cross-sections particularly overestimate the triplet population.

However, even the prediction by MCCC shows a significant underestimation in
the population in $\mathrm{D^{1}\Pi_{u}^{\pm}}$ states. The overpopulation
might be caused by the radiation reabsorption, since
$\mathrm{D^{1}\Pi_{u}^{\pm}}$ states are optically allowed levels from
$\mathrm{X^{1}\Sigma_{g}^{+}}$ states. From the ground state density ($\approx
10^{22}\mathrm{\;m^{-3}}$), the mean free path of the photon from
$\mathrm{D^{1}\Pi_{u}^{\pm}}$ states is in the order of $\lesssim 10^{-3}$ m.
This is much shorter than the plasma size of this experiment ($\approx
10^{-2}$ m), and thus this radiation reabsorption effectively decreases the
radiative decay rate from $\mathrm{D^{1}\Pi_{u}^{\pm}}$ state and increases
its population. Furthermore, the $v$-dependence in
$\mathrm{D^{1}\Pi_{u}^{\pm}}$ is also consistent with the trend of the Franck-
Condon factor with $\mathrm{X^{1}\Sigma_{g}^{+}}$, i.e.,
$q_{\mathrm{X^{1}\Sigma_{g}^{+}}(v=0)-\mathrm{D^{1}\Pi_{u}^{-}}(v=0,1,2)}\approx[0.10,0.17,0.18]$;
since the emission from $v=2$ state of $\mathrm{D^{1}\Pi_{u}^{\pm}}$ has
bigger absorption rate than that from $v=0$ state, the overpopulation in $v=2$
state should be more significant. These are all consistent with the
experiment, although this effect is not included in our CRM analysis.

###  IV.2 LHD experiment

In order to examine the dataset in other parameter spaces, we carried out the
emission observation from LHD divertor region, with the same setting described
in Ref. Ishihara  _et al._ . The light from the divertor region is collected
by an optical lens, introduced into an optical fiber, and guided to the
entrance slit of another echelle spectrometer (home-made, 300-mm focal length
Tanaka  _et al._ (2018)). This spectrometer can simultaneously measure an
emission spectrum in the entire visible range (430 – 800 nm) with the high-
wavelength resolution (0.06 nm).

Figure 11 shows the summary of a typical LHD experiment. The temporal
evolutions of $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$ on the plasma axis of the
LHD are shown in Fig. 11 (a), while Fig. 11 (b) and (c) show the radial
distributions of $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$, respectively, for two
measurement timings. The values of $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$ are
high in the confined region (inside the last closed flux surface, LCFS), while
they become much smaller in the divertor. In Ref. Kobayashi  _et al._ (2010),
it has been found that the values of $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$ at
the divertor, $T_{\mathrm{e}}^{\mathrm{div}}$ and
$n_{\mathrm{e}}^{\mathrm{div}}$, respectively, approximately have the
following relationship with the values on the LCFS,

$\displaystyle\frac{T_{\mathrm{e}}^{\mathrm{div}}}{\mathrm{eV}}$
$\displaystyle\approx
30\times\left[\frac{n_{\mathrm{e}}^{\mathrm{LCFS}}}{10^{19}\mathrm{\;m^{-3}}}\right]^{-0.5}$
(22)
$\displaystyle\frac{n_{\mathrm{e}}^{\mathrm{div}}}{10^{19}\mathrm{\;m^{-3}}}$
$\displaystyle\approx
0.08\times\left[\frac{n_{\mathrm{e}}^{\mathrm{LCFS}}}{10^{19}\mathrm{\;m^{-3}}}\right]^{1.5}$
(23)

It is important to note that, besides the above empirical relation, the exact
$n_{\mathrm{e}}$ and $T_{\mathrm{e}}$ values at the emission locations are not
available. Thus, in the following we use this empirical relation with the
uncertainty of 30%.

Figure 11: (a) A temporal evolution of $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$
at the plasma center of the LHD. The radial distributions of (b)
$T_{\mathrm{e}}$ and (c) $n_{\mathrm{e}}$ for two different timings $t=3.25$ s
and $t=5.25$ s. The position of the LCFS are indicated by the dotted line.

Figure 12 shows the two emission spectra observed for LHD experiment #143973
at $t=3.25$ s and $t=5.25$ s. The shapes of the spectrum are different between
the two timings; the two strong lines in 720–800 nm in the $t=3.25$ s spectrum
disappear in $t=5.25$ s, while new lines appear in 570–600 nm.

The same analysis to that described in the previous section is performed for
all the exposure frames in this experiment. Figure 13 shows the population
obtained from the analysis for the two timings. From $t=3.25$ s to 5.25 s, the
population in $\mathrm{EF^{1}\Sigma_{g}^{+}}$ decreases while those in
$\mathrm{g^{3}\Sigma_{g}^{+}}$, $\mathrm{i^{3}\Pi_{g}^{\pm}}$, and
$\mathrm{j^{3}\Delta_{g}^{\pm}}$ increase. Note that although the values of
$T_{\mathrm{e}}$ and $n_{\mathrm{e}}$ changes during LHD experiments, as the
timescale of the parameter change in the plasma is much longer ($\gtrsim
10^{-3}$ s) than that of the population change ($\lesssim 10^{-7}$ s), the
quasi steady-state approximation for the excited state population is still
valid. The populations in $\mathrm{D^{1}\Pi_{u}^{\pm}}$ and
$\mathrm{d^{3}\Pi_{u}^{\pm}}$ stay almost the same in these two timings. This
is consistent with the parameter dependence on $R_{1}$ as shown in Fig. 6.

Figure 12: A similar figure to Fig. 9, but for the LHD divertor plasma. (top)
Emission spectrum observed at (gray) $t=3.25$ s and (magenta) $t=5.25$ s. To
clearly present the emission lines from molecular hydrogen, strong atomic
lines are omitted. (bottom) The composition of the hydrogen-molecule spectrum
observed at $t=3.25$ s for each upper-electronic states, which is obtained
from the spectral fit.  Figure 13: Population distributions observed for the
LHD divertor at (grey) $t=3.25$ s and (magenta) $t=5.25$ s. Different behavior
in each state can be seen, e.g., the population in
$\mathrm{EF^{1}\Sigma_{g}^{+}}$ state decreases while those in the triplet
states increase.  Figure 14: A similar figure to Fig. 10 but for the LHD
divertor plasmas. Experimentally observed population distribution is shown by
markers, while kinked lines are the prediction by the CRM with three different
excitation cross section (blue: MCCC, green: Miles, orange: Janev). The values
of $n_{\mathrm{e}}$, $T_{\mathrm{e}}$, and $T_{v}$ to calculate the population
were estimated based on Eq. (11) and Eq. (21). The only one adjustable
parameter is the line-integrated molecular density, $n_{\mathrm{H_{2}}}L$, is
adjusted so that the population in $\mathrm{EF^{1}\Sigma_{g}^{+}}$ $(v=6)$ is
identical between the observation and the prediction.

Figure 14 shows the comparison with our CRM predictions. As similar to the
discussion for Fig. 10, we used fixed $n_{\mathrm{e}}$, $T_{\mathrm{e}}$, and
$T_{v}$ values based on Eq. (22) and Eq. (21). To account for the uncertainty
in these parameters, we assume independent 10% error bands for these three
parameters. The value of $n_{\mathrm{H_{2}}}L$ is estimated so that the
$\mathrm{EF^{1}\Sigma_{g}^{+}}$ ($v=6$) population by the experiment and
predictions match.

The predictions with the cross sections by MCCC and Janev equally well
represent the experimental observations. Miles’ prediction show the
significant underestimation on the $\mathrm{J^{1}\Delta_{g}^{\pm}}$,
$\mathrm{h^{3}\Sigma_{g}^{+}}$, and $\mathrm{d^{3}\Pi_{u}^{\pm}}$ populations.

Even with the CRM using the MCCC dataset, the fit is worse than that in the
low-density RF plasma Fig. 10. The populations in
$\mathrm{GK^{1}\Sigma_{g}^{+}}$, $\mathrm{i^{3}\Pi_{g}^{\pm}}$, and
$\mathrm{j^{3}\Delta_{g}^{\pm}}$ are not well reproduced by the CRM, in
particular in the higher density plasma (the lower panel of Fig. 14). The
reason might be the uncertain cross sections from the excited state. We will
discuss this issue in the next section.

Figure 15: (a) A temporal evolution of $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$
on the LCFS for LHD experiment #143973. (b) The temporal evpolutions of the
populations in $\mathrm{EF^{1}\Sigma_{g}^{+}}$ $(v=6)$,
$\mathrm{D^{1}\Pi_{u}^{-}}$ $(v=0)$, and $\mathrm{d^{3}\Pi_{u}^{-}}$ $(v=0)$.
Not only the population but also the population ratio changes in time,
according to the change in $n_{\mathrm{e}}$ and $T_{\mathrm{e}}$.

We also conduct the above emission measurement and analysis for the spectra
observed in other timings. Figure 15 shows the temporal evolution of the
population of several excited states, for this particular LHD experiment. The
population and the population ratios change in time. In the first frame where
the plasma has the lowest density and highest temperature during the
discharge, the population in $\mathrm{EF^{1}\Sigma_{g}^{+}}$ is highest while
that in $\mathrm{D^{1}\Pi_{u}^{-}}$ is lowest. While $n_{\mathrm{e}}$
increases as time goes by, the $\mathrm{EF^{1}\Sigma_{g}^{+}}$ population
decreases while $\mathrm{D^{1}\Pi_{u}^{-}}$ population increases. At the final
frame, where the temperature dropped significantly, the
$\mathrm{D^{1}\Pi_{u}^{-}}$ population decreases significantly while the
$\mathrm{d^{3}\Pi_{u}^{-}}$ population stays almost the same. This is
consistent with the energy dependence of the cross sections; the cross section
to the triplet state has higher efficiency at the lower temperature (see Fig.
4 and Fig. 5).

Figure 16: A distribution of population ratios observed for the LHD divertor
and the low-density RF plasma. (a) A correlation between the two population
ratios, $\mathrm{EF^{1}\Sigma_{g}^{+}}$ $(v=6)$ / $\mathrm{D^{1}\Pi_{u}^{-}}$
$(v=0)$ and $\mathrm{d^{3}\Pi_{u}^{-}}$ $(v=0)$ / $\mathrm{D^{1}\Pi_{u}^{-}}$
$(v=0)$. (b, c) Histograms of the observed population ratios and the predicted
population ratios by our CRM. On the left of each panel, the predictions by
the three different cross-section datasets (blue: MCCC, green: Miles, orange:
Janev) are shown. On the right, we compare the predictions by the CRM (blue)
without and (black) with the recommended electron-attachment rates. We assume
$T_{v}=0.2$–$1.0\;\mathrm{eV}$ and (b) $T_{\mathrm{e}}=5$–$20\;\mathrm{eV}$
and (c) $n_{\mathrm{e}}=(0.5$–$2)\times 10^{19}\;\mathrm{m^{-3}}$ for the CRM.

The same analysis are repeated for multiple LHD experiments spanning wide
range in the parameter space. Figure 16 shows the distribution of some
population ratios. The star markers in the figure shows the same population
ratios from the low-temperature RF plasma.

In Fig. 16 (b) and (c) we compare the observed population ratios and the CRM
prediction. On the left side of each panel, the CRM prediction with the three
different excitation cross-sections datasets are shown in different colors. On
the right half of the panels, the CRM predictions with and without the
dissociative electron-attachment are plotted. Note that we assume
$T_{v}=0.2$–$1.0\;\mathrm{eV}$ and (b) $T_{\mathrm{e}}=5$–$20\;\mathrm{eV}$
and (c) $n_{\mathrm{e}}=(0.5$–$2)\times 10^{19}\;\mathrm{m^{-3}}$ for the CRM.

The population ratios between $\mathrm{EF^{1}\Sigma_{g}^{+}}$ and
$\mathrm{D^{1}\Pi_{u}^{-}}$ shows a sensitivity on $n_{\mathrm{e}}$, as
expected from the discussion in Sec.III. The CRM with the Miles’ dataset or
Janev’s dataset fail to reproduce the range of the observed population ratios.
This $n_{\mathrm{e}}$-sensitivity disappears if we include the electron-
attachment process, and fails to reproduce the observed population ratios.

The population ratios between $\mathrm{d^{3}\Pi_{u}^{-}}$ and
$\mathrm{D^{1}\Pi_{u}^{-}}$ shows a $T_{\mathrm{e}}$-sensitivity. Although
this is not perfectly independent from $n_{\mathrm{e}}$, this suggests that
this ratio might be useful to diagnose the value of $T_{\mathrm{e}}$ from the
molecular emission spectrum.

##  V Summary and conclusions

We developed a CRM for hydrogen molecules, and discussed the effects of
different datasets. The CRM prediction was compared with two different
experiments, which is the lower-density RF plasmas and higher-density LHD
divertor plasmas. From the population comparisons the following is found;

  1. 1.

Compared with Miles’ and Janev’s cross sections, MCCC cross sections show a
better agreement with the experiment.

  2. 2.

CRM prediction is more accurate in the lower density plasma, where the
excitation from excited states is negligible. In the higher density plasmas
this effect is more significant. Cross-sections from excited states might be
necessary to improve for the high density plasma diagnostics.

  3. 3.

Dissociative electron-attachment rate proposed in Ref. yac  might be too high.

  4. 4.

Some pairs of the emission lines of hydrogen molecule show $n_{\mathrm{e}}$
and $T_{\mathrm{e}}$ dependencies, e.g.,
$\mathrm{EF^{1}\Sigma_{g}^{+}}$-$\mathrm{D^{1}\Pi_{u}^{\pm}}$ has a
$n_{\mathrm{e}}$ dependence, and
$\mathrm{D^{1}\Pi_{u}^{\pm}}$-$\mathrm{d^{3}\Pi_{u}^{\pm}}$ has a
$T_{\mathrm{e}}$ dependence. This may be useful for the plasma diagnostics for
low-temperature hydrogen plasmas.

Since most of the important data is already available for other isotopologues,
such as $\mathrm{D}_{2}$, the same analysis can be done also. This is left for
the future study.

###### Acknowledgements.

This work was partly supported by the U.S. D.O.E contract DE-AC05-00OR22725,
the Australian Government through the Australian Research Council’s Discovery
Projects funding scheme (project DP240101184), and by the United States Air
Force Office of Scientific Research. K.F. would like to specifically
acknowledge Oak Ridge National Laboratory’s Laboratory Directed Research and
Development program Project No. 11367. L.H.S is the recipient of an Australian
Research Council Discovery Early Career Researcher Award (project number
DE240100176) funded by the Australian Government. HPC resources were provided
by the Pawsey Supercomputing Research Centre and the National Computational
Infrastructure, with funding from the Australian Government and the Government
of Western Australia, and the Texas Advanced Computing Center (TACC) at the
University of Texas at Austin. M.C.Z would like to specifically acknowledge
the support of the Los Alamos National Laboratory (LANL) ASC PEM Atomic
Physics Project. LANL is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of the U.S. Department of Energy
under Contract No. 89233218NCA000001. M.C.Z. would like to specifically
acknowledge Los Alamos National Laboratory’s Laboratory Directed Research and
Development program Project No. 20240391ER.

## References

  * Rosenthal  _et al._ (2000) D. Rosenthal, F. Bertoldi, and S. Drapatz, Astronomy and astrophysics 356, 705 (2000).
  * Lieberman and Lichtenberg (2005) M. A. Lieberman and A. J. Lichtenberg,  _Principles of Plasma Discharges and Materials Processing_ (Wiley, 2005).
  * Janev  _et al._ (2003) R. K. Janev, U. Samm, and D. Reiter,  _Collision processes in low-temperature hydrogen plasmas_, Tech. Rep. PreJuSER-38224 (Institut für Plasmaphysik, 2003).
  * Fantz  _et al._ (2001) U. Fantz, D. Reiter, B. Heger, and D. Coster, Journal of Nuclear Materials 14th Int. Conf. on Plasma-Surface Interactions in Controlled Fusion D evices, 290-293, 367 (2001).
  * (5) H. Ishihara, A. Kuzmin, M. Kobayashi, T. Shikama, K. Sawada, S. Saito, H. Nakamura, K. Fujii, and M. Hasuo,  267, 107592.
  * Stangeby and Others (2000) P. C. Stangeby and Others,  _The plasma boundary of magnetic fusion devices_, Vol. 224 (Institute of Physics Pub. Philadelphia, Pennsylvania, 2000).
  * Ohno  _et al._ (1998) N. Ohno, N. Ezumi, S. Takamura, S. I. Krasheninnikov, and A. Y. Pigarov, Physical review letters 81, 818 (1998).
  * Verhaegh  _et al._ (2022) K. Verhaegh, B. Lipschultz, J. Harrison, N. Osborne, A. Williams, P. Ryan, J. Clark, F. Federici, B. Kool, T. Wijkamp, A. Fil, D. Moulton, O. Myatra, A. Thornton, T. Bosman, G. Cunningham, B. Duval, S. Henderson, R. Scannell, and the MAST Upgrade team,  (2022), 10.1088/1741-4326/aca10a, arXiv:2204.02118 [physics.plasm-ph] .
  * (9) H. P. Summers, “The adas user manual, version 2.6,” http://www.adas.ac.uk.
  * Sawada and Fujimoto (1995) K. Sawada and T. Fujimoto, Journal of applied physics 78, 2913 (1995).
  * Greenland (2002) P. T. Greenland, Contributions to Plasma Physics 42, 608 (2002).
  * Lavrov  _et al._ (2006) B. P. Lavrov, A. V. Pipa, and J. Röpcke, Plasma Sources Science and Technology 15, 135 (2006).
  * Guzmán  _et al._ (2013) F. Guzmán, M. O’Mullane, and H. P. Summers, Journal of Nuclear Materials 438, S585 (2013).
  * Shakhatov  _et al._ (2016) V. A. Shakhatov, Y. A. Lebedev, A. Lacoste, and S. Bechu, High Temperature 54, 124 (2016).
  * Sawada and Goto (2016) K. Sawada and M. Goto, Atoms for Peace, an International Journal 4, 29 (2016).
  * Wünderlich  _et al._ (2021) D. Wünderlich, L. H. Scarlett, S. Briefi, U. Fantz, M. C. Zammit, D. V. Fursa, and I. Bray, Journal of physics D: Applied physics 54, 115201 (2021).
  * (17) “Yacora on the Web,” https://yacora.de/, [Accessed: 3-Oct-2023].
  * Miles  _et al._ (1972) W. T. Miles, R. Thompson, and A. E. S. Green, Journal of applied physics 43, 678 (1972).
  * (19) “”mccc database”,” https://mccc-db.org.
  * Scarlett  _et al._ (2017) L. H. Scarlett, J. K. Tapley, D. V. Fursa, M. C. Zammit, J. S. Savage, and I. Bray, PHYSICAL REVIEW A 96, 62708 (2017).
  * Scarlett  _et al._ (2021a) L. H. Scarlett, J. S. Savage, D. V. Fursa, I. Bray, M. C. Zammit, and B. I. Schneider, Physical review. A 103, 032802 (2021a).
  * Scarlett  _et al._ (2021b) L. H. Scarlett, I. Bray, and D. V. Fursa, Physical review letters 127, 223401 (2021b).
  * Scarlett  _et al._ (2021c) L. H. Scarlett, D. V. Fursa, M. C. Zammit, I. Bray, and Y. Ralchenko, Atomic Data and Nuclear Data Tables 139, 101403 (2021c).
  * Scarlett  _et al._ (2021d) L. H. Scarlett, D. V. Fursa, M. C. Zammit, I. Bray, Y. Ralchenko, and K. D. Davie, Atomic Data and Nuclear Data Tables 137, 101361 (2021d).
  * Scarlett  _et al._ (2022) L. H. Scarlett, D. K. Boyle, M. C. Zammit, Y. Ralchenko, I. Bray, and D. V. Fursa, Atomic Data and Nuclear Data Tables 148, 101534 (2022).
  * Scarlett  _et al._ (2023) L. H. Scarlett, E. Jong, S. Odelia, M. C. Zammit, Y. Ralchenko, B. I. Schneider, I. Bray, and D. V. Fursa, Atomic Data and Nuclear Data Tables , 101573 (2023).
  * Nakashima and Nakatsuji (2018) H. Nakashima and H. Nakatsuji, The Journal of chemical physics 149, 244116 (2018).
  * Wolniewicz (1993) L. Wolniewicz, The Journal of chemical physics 99, 1851 (1993).
  * Fantz and Wünderlich (2006) U. Fantz and D. Wünderlich, Atomic Data and Nuclear Data Tables 92, 853 (2006).
  * Astashkevich  _et al._ (1996) S. A. Astashkevich, M. Käning, E. Käning, N. V. Kokina, B. P. Lavrov, A. Ohl, and J. Röpcke, Journal of quantitative spectroscopy & radiative transfer 56, 725 (1996).
  * (31) R. Mewe,  .
  * Comtet and De Bruijn (1985) G. Comtet and D. P. De Bruijn, Chemical physics 94, 365 (1985).
  * Astashkevich and Lavrov (2015) S. A. Astashkevich and B. P. Lavrov, Journal of Physical and Chemical Reference Data 44, 023105 (2015).
  * Wedding and Phelps (1988) A. B. Wedding and A. V. Phelps, The Journal of chemical physics 89, 2965 (1988).
  * Datskos  _et al._ (1997) P. G. Datskos, L. A. Pinnaduwage, and J. F. Kielkopf, Physical review. A 55, 4131 (1997).
  * Wünderlich  _et al._ (2020) D. Wünderlich, M. Giacomin, R. Ritz, and U. Fantz, Journal of quantitative spectroscopy & radiative transfer 240, 106695 (2020).
  * Sawada  _et al._ (2010) K. Sawada, Y. Yamada, T. Miyachika, N. Ezumi, A. Iwamae, and M. Goto, Plasma and Fusion Research 5, 001 (2010).
  * Brezinsek  _et al._ (2002) S. Brezinsek, P. Mertens, A. Pospieszczyk, G. Sergienko, and P. T. Greenland, Contributions to Plasma Physics 42, 668 (2002).
  * Tanaka  _et al._ (2018) M. Tanaka, D. Kato, G. Gaigalas, P. Rynkun, L. RadžiÅ«tÄ, S. Wanajo, Y. Sekiguchi, N. Nakamura, H. Tanuma, I. Murakami, and H. A. Sakaue, The Astrophysical Journal 852, 109 (2018), arXiv:1708.09101 .
  * Kobayashi  _et al._ (2010) M. Kobayashi, Y. Feng, S. Morita, S. Masuzaki, N. Ezumi, T. Kobayashi, M. B. Chowdhuri, H. Yamada, T. Morisaki, N. Ohyabu, M. Goto, I. Yamada, K. Narihara, A. Komori, and O. Motojima, Fusion Science and Technology 58, 220 (2010).


# Ordinals and recursive algorithms

Lior Shiboli, Gabriel Nivasch Ariel University, Ariel, Israel.
<EMAIL_ADDRESS>University, Ariel, Israel<EMAIL_ADDRESS>
###### Abstract

We prove the following result, which generalizes results of Erickson et al.
(2022) and Bufetov et al. (2024): Call a function $f:\mathbb{R}\to\mathbb{R}$
_ordinal decreasing_ if for every infinite decreasing sequence
$x_{0}>x_{1}>x_{2}>\cdots$ there exist $i<j$ such that $f(x_{j})\geq
f(x_{i})$. Given ordinal decreasing functions $f,g_{0},\ldots,g_{k},s$ that
are larger than $0$, define the recursive algorithm “$M(x)$: if $x<0$ return
$f(x)$, else return $g_{0}(-M(x-g_{1}(-M(x-\cdots-
g_{k}(-M(x-s(x)))\cdots))))$”. Then $M(x)$ halts and is ordinal decreasing for
all $x\in\mathbb{R}$.

##  1 Introduction

Erickson, Nivasch and Xu [3, 4, 5] studied the following recursive algorithm
$M$:

$M(x)=\begin{cases}-x,&\text{if $x<0$;}\\\ \frac{M(x-M(x-1))}{2},&\text{if
$x\geq 0$.}\end{cases}$ (1)

For example it can be checked that $M(1)=\frac{M(1-M(0))}{2}=\frac{1}{8}$,
$M(\frac{3}{2})=\frac{1}{32}$, $M(2)=2^{-10}$, $M(\frac{5}{2})=2^{-31}$ and
$M(3)=2^{-1,541,023,937}$ (see Figure 1).

Erickson, Nivasch and Xu [3, 4, 5] proved that $M$ terminates on all real
inputs, although Peano Arithmetic cannot prove that $M$ terminates on all
natural inputs. PA-independence was shown by proving that $\frac{1}{M(n)}$
grows as fast as $F_{\epsilon_{0}}(n-7)$ for integers $n\geq 8$.

The motivation for algorithm $M$ lies in the set of _fusible numbers_. As
Erickson et al. [3] showed, $M(x)$ returns the distance between $x$ and the
smallest “tame fusible number” larger than $x$. However, algorithm $M$ is
worth studying on its own right, since it is a simple algorithm for which its
termination is not so easy to prove.

Figure 1: The graph of $M(x)$ up to $x=2.1$.

In a follow-up paper, Bufetov et al. [1] studied a generalization of fusible
numbers to _$n$ -fusible numbers_ and a corresponding generalization of
algorithm $M$ to the following algorithm $M_{n}$:

$M_{n}(x)=\begin{cases}-x,&\text{if $x<0$;}\\\ \frac{M_{n}(x-M_{n}(x-\cdots-
M_{n}(x-1)\cdots))}{n},&\text{if $x\geq 0$;}\end{cases}$ (2)

where $M_{n}(x-\cdots)$ repeats $n$ times. They showed that for every $n\geq
1$, $M_{n}$ terminates on all real inputs. In this paper we study the
following question: Which modifications can be done to algorithm $M$ such that
it will still halt on all real inputs?

We identify a large class of algorithms that halt on all real inputs. The
description of the algorithms, as well as the proof that they halt on all real
inputs, involve a property of functions, which we call _ordinal decreasing_.

###  1.1 Ordinal decreasing functions

######  Definition 1.

Given a function $f:D\to\mathbb{R}$ for some $D\subseteq\mathbb{R}$, we call
$f$ _ordinal decreasing_ if for every infinite decreasing sequence
$x_{0}>x_{1}>x_{2}>\cdots$ in $D$ there exist $i<j$ (so that $x_{j}<x_{i}$)
such that $f(x_{j})\geq f(x_{i})$. A function is called ordinal decreasing up
to $a$ if its ordinal decreasing in the interval $(-\infty,a]\subseteq D$.

Some examples of ordinal decreasing functions are:

  1. 1.

Every nonincreasing function.

  2. 2.

A function $f$ in the domain $[0,1)$ that is divided into a sequence of
intervals $\\{I_{n}\\}$ where $I_{n}=[1-2^{-n},1-2^{-n+1})$ and where $f$ is
decreasing in each interval (see Figure 2, left).

  3. 3.

A function $f$ in the domain $(0,1]$ that is divided into a sequence of
intervals $\\{I_{n}\\}$ where $I_{n}=(2^{-(n+1)},2^{-n}]$, such that $f$ is
decreasing in each $I_{n}$ but such that for each $n$ there exist only
finitely many $n<m$ such that $2^{-m}>\lim_{x\to(2^{-n})^{+}}(f(x))$ (see
Figure 2, right).

Figure 2: Examples of ordinal decreasing functions.

###  1.2 Real induction

In this paper we will use the following result, which is called _real
induction_ (see Clark [2] for a survey).

######  Lemma 1.

Let $S\subset\mathbb{R}$ be a set that satisfies:

(R1) There exists $a\in\mathbb{R}$ such that $(-\infty,a)\subset S$.

(R2) For all $x\in\mathbb{R}$, if $(-\infty,x)\subset S$, then $x\in S$.

(R3) For all $x\in S$, there exists $y>x$ such that $(x,y)\subset S$.

Then $S=\mathbb{R}$ .

###### Proof.

Suppose $S\neq\mathbb{R}$. Let $a=\inf(\mathbb{R}\setminus S)$. By (R1)
$a\neq-\infty$. Therefore by (R2), $a\in S$. Therefore (R3) yields a
contradiction. ∎

It is also worth noting for our purposes that since Peano Arithmetic is built
upon the natural numbers, we cannot use real induction within Peano
Arithmetic, but must rely on Second Order Arithmetic.

###  1.3 Our result

Now we can state the main result of our paper.

######  Theorem 1.

Consider the recursive algorithm:

$M(x)=\begin{cases}f(x)&\text{if $x<0$;}\\\ g_{0}(-M(x-g_{1}(-M(x-\cdots-
g_{k}(-M(x-s(x)))\cdots))))&\text{if $x\geq 0$;}\end{cases}$ (3)

where the functions $s(x)$, $f(x)$ and $g_{i}(x)$ for all $i$ are all ordinal
decreasing and larger than 0 for every $x$ in the appropriate ranges:
$(-\infty,0)$ for $f,g_{i},h$, and $[0,\infty)$ for $s$.

Then $M(x)$ halts and is ordinal decreasing for all $x\in\mathbb{R}$.

Note that Theorem 1 covers the cases mentioned above, by taking $k=n$,
$g_{0}(x)=-\frac{x}{n}$, $g_{i}(x)=-x$ for $1\leq i\leq n$ and $s(x)=1$.

##  2 Proof of Theorem 1

We start by proving some properties of ordinal decreasing functions.

######  Claim 1.

Suppose $f(x)$ is ordinal decreasing up to $y$. Then for every infinite
decreasing sequence $\\{x_{n}\\}$ with $x_{0}<y$ there is an infinite
subsequence $\\{x^{\prime}_{n}\\}$ for which $\\{f(x^{\prime}_{n})\\}$ is
nondecreasing.

###### Proof.

By the infinite Ramsey’s theorem. Define an infinite complete graph in which
there is a vertex for each $x_{i}$ and color each edge $\\{x_{i},x_{j}\\},i<j$
red if $f(x_{i})>f(x_{j})$ and green otherwise.

Since $f(x)$ is ordinal decreasing up to $y$ our graph cannot contain a
monochromatic red infinite complete subgraph. Therefore there exists a
monochromatic green infinite subgraph, and thus the original sequence contains
an infinite nondecreasing subsequence (comprised of all the vertices in the
subgraph). ∎

######  Claim 2.

Suppose $g(x)$ is ordinal decreasing up to $y$ and $f(x)$ is ordinal
decreasing up to $\max_{x<y}(-g(x))$. Then $f(-g(x))$ is ordinal decreasing up
to $y$.

###### Proof.

Consider an infinite decreasing sequence $\\{x_{n}\\}$ with $x_{0}<y$. By
Claim 1, there exists a nondecreasing subsequence of $\\{g(x_{n})\\}$. If this
subsequence is not strictly increasing, there exists $i<j$ such that
$g(x_{i})=g(x_{j})$, and so is $f(-g(x_{i}))=f(-g(x_{j}))$. Otherwise we have
a strictly decreasing sequence of $-g(x)$, and since $f(x)$ is ordinal
decreasing up to $\max_{x<y}(-g(x))$, we can find $f(-g(x_{j}))\geq
f(-g(x_{i}))$, and we are done. ∎

######  Claim 3.

Suppose $f(x)$ and $g(x)$ are ordinal decreasing up to $y$. Then $f(x)+g(x)$
is ordinal decreasing up to $y$.

###### Proof.

By Claim 1, for every strictly decreasing sequence $\\{x_{n}\\}$ with
$x^{\prime}_{0}<y$ we can find a subsequence $\\{x^{\prime}_{n}\\}$ such that
$\\{f(x^{\prime}_{n})\\}$ is nondecreasing. For that subsequence we can find
by definition $i<j$ such that $g(x_{j})\geq g(x_{i})$. Hence,
$f(x^{\prime}_{j})+g(x^{\prime}_{j})\geq f(x^{\prime}_{i})+g(x^{\prime}_{i})$
as desired. ∎

######  Claim 4.

Let $\beta>0$, and suppose $f(x)$ is ordinal decreasing up to $y+\beta$ and
larger than $0$. Then there exists an $0<\epsilon<\beta$ such that $x-f(x)<y$
for every $y<x<y+\epsilon$.

###### Proof.

Suppose for a contradiction that for every $0<\epsilon<\beta$ we have a
counterexample $x-f(x)\geq y$ with $y<x<y+\epsilon$. Then, we have an infinite
sequence $\\{x_{n}\\}$ of such counterexamples with $\lim(x_{n})=y$, but
because $0<f(x_{n})\leq x_{n}-y$ there exists an infinite subsequence
$\\{x^{\prime}_{n}\\}$ for which $f(x^{\prime}_{n})$ is decreasing, in
contradiction to $f(x)$ being ordinal decreasing up to $y+\beta$. ∎

We are ready to prove our main result.

###### Proof of Theorem 1.

Consider the recursive algorithm

$M(x)=\begin{cases}f(x),&\text{if $x<0$}\\\ g_{0}(-M(x-g_{1}(-M(x-\cdots-
g_{k}(-M(x-s(x)))\cdots))))&\text{if $x\geq 0$}\end{cases}$ (4)

where $s(x)$, $f(x)$ and $g_{i}(x)$ for all $i$ are all ordinal decreasing and
larger than 0 for every $x$ in the appropriate ranges: $(-\infty,0)$ for
$f,g_{i},h$, and $[0,\infty)$ for $s$.

We claim that $M(x)$ halts and is ordinal decreasing for every $x$. We will
prove this by real induction (Lemma 1).

Assuming otherwise, let

$S=\\{x\mid\text{$M$ is defined and ordinal decreasing up to $x$}\\}.$

Since for $x<0$ $M(x)$ is defined by $M(x)=f(x)$, we have $(-\infty,0)\in S$.
Hence, S satisfies property (R1).

Next, suppose that $(-\infty,y)\in S$. Then note that $M(y)$ is defined, since
for every $i$ $M(y-g_{i}(-M(y-\cdots-g_{k}(-M(y-s(y)))\cdots)))$ is defined by
induction, since functions $g_{i}$ and $s$ have output larger than $0$. Hence,
$M(x)$ is ordinal decreasing up to $y$ itself. Hence $y\in S$ as well, so $S$
satisfies property (R2).

Finally suppose $y\in S$. We will show that $(y,y+\epsilon)\subseteq S$ for
some $\epsilon>0$, meaning $S$ satisfies property (R3).

We will do so by induction on $i=k,\ldots,0$ that $M_{i}=g_{i}(M(x-\cdots-
g_{k}(M(x-s(x)))\cdots))$ is defined and ordinal decreasing up to
$y+\epsilon_{i}$ for some $\epsilon_{i}>0$. Let us start with the base case
$i=k$. In this case $M_{k}=g_{k}(-M(x-s(x)))$. By Claim 3, $-x+s(x)$ is an
ordinal decreasing function in $[0,\infty)$, hence by Claim 4 on $s(x)$, there
is an $\epsilon_{k}$ such that $x-s(x)<y$ for every $y<x<y+\epsilon_{k}$.
Hence, by assumption and Claim 2, $M(-(-x+s(x)))=M(x-s(x))$ is ordinal
decreasing and defined up to $y+\epsilon_{k}$ and so is $g_{k}(-M(x-s(x)))$,
as desired. For the induction step, suppose $M_{i}(x)$ is ordinal decreasing
and defined up to $y+\epsilon_{i}$. By Claim 3, $-x+M_{i}(x)$ is ordinal
decreasing and defined up to $y+\epsilon_{i}$. Furthermore, by Claim 4 there
exists some $\epsilon_{i-1}$ such that $x-M_{i}(x)<y$ for every
$y<x<y+\epsilon_{i-1}$. Hence, by assumption and Claim 2 $M(x-M_{i}(x))$ is
defined and ordinal decreasing up to $y+\epsilon_{i-1}$. Hence, so is
$M_{i}=g_{i}(M(x-M_{i}(x)))$. Hence by Lemma 1 we have $S=\mathbb{R}$. ∎

##  3 Discussion

In this paper we have shown that the algorithm $M$ of Theorem 1 halts for
every real input, but the structure of the function it returns is still left
to be studied.

######  Definition 2.

Given an ordinal decreasing function $f$, define $O_{f}$ to be the function
from the reals to the ordinals given by
$O_{f}(x)=\lim_{x^{\prime}<x,f(x^{\prime})<f(x)}(O_{f}(x^{\prime})+1)$. Let
$O_{f}(\infty)=\lim_{x\in\mathbb{R}}O_{f}(x)$ be the ordinal of $f$.

It is easy to show that $f:\mathbb{R}\to\mathbb{R}$ is ordinal decreasing up
to $y$ if and only if $O_{f}(x)$ is defined for every $x<y$.

It was shown by [4] that $O_{M_{2}}(\infty)=\epsilon_{0}$, and by [1] that
more generally $O_{M_{n}}(\infty)=\varphi_{n-1}(0)$. It remains to be studied
if given the ordinals of $s,g_{i}$ and $f$ we can determine an upper bound for
the ordinal of $M$.

## References

  * [1] Alexander I. Bufetov, Gabriel Nivasch, and Fedor Pakhomov.  Generalized fusible numbers and their ordinals.  Annals of Pure and Applied Logic, 175(1, Part A), 2024.
  * [2] Pete L. Clark.  The instructor’s guide to real induction, 2012.  arXiv e-prints, math.HO, 1208.0973.
  * [3] Jeff Erickson.  Fusible numbers.  https://www.mathpuzzle.com/fusible.pdf.
  * [4] Jeff Erickson, Gabriel Nivasch, and Junyan Xu.  Fusible numbers and Peano Arithmetic.  Logical Methods in Computer Science, 18(3), 2022.
  * [5] Junyan Xu.  Survey on fusible numbers, 2012.  arXiv e-prints, math.CO, 1202.5614.