of $\mathcal{N}$, and all pieces are atoroidal.

  4. (4)

Existence of a non-positively curved metric on $\mathcal{N}$: As the manifold
$\mathcal{N}$ is $Haken$, i,e it is irreducible, and contains a closed
incompressible surface, and it has one atoroidal piece, we can define a non-
positively curved Riemannian metric on $\mathcal{M}$ by [Lee95, Theorem 3.3].
This opens a potential new direction to explore as described below.

###  6.1. Future questions

By item (4) above, we can define a non-positively curved Riemannian metric on
$\mathcal{N}$. Then the universal cover $\widetilde{\mathcal{N}}$ with the
path metric induced by the non-positively curved Riemannian metric is
$CAT(0)$. For a $CAT(0)$-space we can define different types of ‘boundaries at
infinity’,for example the Tits boundary or the Morse boundary and we can
define topologies on the boundary at infinity,
$\partial_{\infty}\widetilde{\mathcal{N}}$. For example, one question to ask
is whether the weak-stable or weak-unstable foliations of $\Psi_{t}$ satisfy
the Continuous Extension Property:

######  Question 6.1.

Suppose $F\in\widetilde{\mathcal{F}}^{ws}$ or $\widetilde{\mathcal{F}}^{wu}$
of $\widetilde{\Psi}_{t}$. It is known that leaves of the weak-stable or weak-
stable foliations of Anosov flows are Gromov hyperbolic and we can define the
Gromov boundary $S^{1}(F)$ (it is a circle). Then does the inclusion
$i\mathrel{\mathop{\ordinarycolon}}F\rightarrow\widetilde{\mathcal{N}}$
extends continuously to a map $\hat{i}\mathrel{\mathop{\ordinarycolon}}F\cup
S^{1}(F)\rightarrow\widetilde{\mathcal{N}}\cup\partial_{\infty}\widetilde{\mathcal{N}}$?

Similar type of questions have been extensively studied for Anosov flows on
hyperbolic $3$-manifolds. A big difference in the case of Franks-Williams
manifolds is that the Morse lemma is not true in general, but in some
particular cases there is possibility that the quasigeodesic flow rays are
boundedly away from actual geodesic rays. Hence this question is relevant for
quasigeodesic Anosov flow on non-positively curved manifolds.

## References

  * [Ano69] D. V. Anosov, Geodesic flows on closed Riemannian manifolds with negative curvature, Proc. Steklov Inst. Math. 90 (1969).
  * [Bar96] Barbot, T., Flots d’Anosov sur les variétés graphées au sens de Waldhausen, Ann. Inst. Fourier (Grenoble) 46 (1996) 1451-1517.
  * [BBY17] François Béguin, Christian Bonatti, and Bin Yu. Building Anosov flow on $3$-manifolds, Geom. Topol. 21 (2017) 1837-1930.
  * [Cal06] D. Calegari. Universal circles for quasigeodesic flows. Geom. Topol., 10 (2006) 2271-2298.
  * [doC88] do Carmo, Manfredo, Riemannian geometry, Birkhäuser, 1988.
  * [Fen94] S. Fenley, Anosov flows in $3$-manifolds, Ann. of Math. 139 (1994) 79-115.
  * [Fen95] S. Fenley, Quasigeodesic Anosov flows and homotopic properties of flow lines, Jour. Diff. Geom. 41 (1995) 479-514.
  * [Fen16] S. Fenley, Quasigeodesic pseudo-Anosov flows in hyperbolic $3$-manifolds and connections with large scale geometry, Advances in Math. 303 (2016) 192-278.
  * [Fen22] S. Fenley, Non R-covered Anosov flows in hyperbolic 3-manifolds are quasigeodesic, math.arxiv.2210.09238.
  * [FM01] S. Fenley and L. Mosher, Quasigeodesic flows in hyperbolic $3$-manifolds, Topology 40 (2001) 503-537.
  * [FW80] John Franks and Bob Williams. Anomalous Anosov flows. Global theory of dynamical systems (Proc. Internat. Conf., Northwestern Univ., Evanston, Ill., 1979), volume 819 of Lecture Notes in Math., pages 158-174. Springer, Berlin, 1980.
  * [Ghy84] E. Ghys, Flots d’Anosov sur les 3-variétés fibrées en cercles, Ergodic Theory Dynam. Systems 4 (1984), no. 1, 67-80.
  * [Gro87] M. Gromov, Hyperbolic groups, in Essays in group theory, 75–263, Math. Sci. Res. Inst. Publ. 8 (1987) Springer, New York.
  * [Hem76] J. Hempel, 3-manifolds, Ann. of Math. Studies 86, Princeton Univ. Press, 1976.
  * [KH95] Anatole Katok and Boris Hasselblatt.Introduction to the modern theory of dynamical systems, volume 54, Encyclopedia of Mathematics and its Applications. Cambridge University Press, Cambridge, 1995. With a supplementary chapter by Katok and Leonardo Mendoza.
  * [KL98] M. Kapovich, B. Leeb, $3$-manifold groups and nonpositive curvature, Geom. Funct. Anal. 8 (1998) 841-852.
  * [Lee95] B. Leeb. $3$-manifolds with(out) metrics of nonpositive curvature. Invent. Math., 122 (1995) 277-289.
  * [Mos92a] L. Mosher, Examples of quasigeodesic flows on hyperbolic $3$-manifolds, Proceedings of the Ohio State Research semester on Low-Dimensional topology, W. de Gruyter, 1992.
  * [Mos92b] L. Mosher, Dynamical systems and the homology norm of a $3$-manifold II, Inven. Math. 107 (1992) 243-281.
  * [Ngu19] H. T. Nguyen, Distortion of surfaces in $3$-manifolds, J. of Topol. 12 (2019) 1115-1145.
  * [Sma67] S. Smale, Differential dynamical systems, Bull. A.M.S. 73 (1982) 357-381.
  * [Thu82] W. Thurston, The geometry and topology of 3-manifolds, Princeton University Lecture Notes, 1982.
  * [Wal68] Waldhausen, Friedhelm, On irreducible 3-manifolds which are sufficiently large. Ann. of Math. (2) 87 (1968), 56-88.
  * [Wil70] R. F. Williams. The DA maps of Smale and structural stability. In Global Analysis (Proc.Sympos. Pure Math., Vol. XIV, Berkeley, California 1968), pages 329-334. Amer. Math. Soc., Providence, R.I., 1970.
  * [Zeg93] A. Zeghib, Sur les feuilletages géodésiques de variétés hyperboliques, Inven. Math. 114 (1993) 193-206.


#  Stable Energy Distribution of Weakly Dissipative Gasses
under Collisional Energy Cascades

Keisuke Fujii<EMAIL_ADDRESS>Jun Imano  Arseniy Kuzmin  Taiichi
Shikama  Masahiro Hasuo  Department of Mechanical Engineering and Science,
Graduate School of Engineering, Kyoto University Kyoto 615-8540, Japan

###### Abstract

Collisional thermalization of a particle ensemble under the energy dissipation
can be seen in variety of systems, such as heated granular gasses and
particles in plasmas. Despite its universal existence, analytical descriptions
of the steady-state distribution have been missing. Here, we show that the
steady-state energy distribution of the wide class of collisional energy
cascades can be well approximated by the generalized Mittag-Leffler
distribution, which is one of stable distributions. This distribution has a
power-law tail, as similar to Levy’s stable distribution, the index of which
is related to the energy dissipation rate. We demonstrate its universality by
comparing Mont-Carlo simulations of dissipative gasses as well as the
spectroscopic observation of the atom velocity distribution in a low-
temperature plasma.

Nonthermal energy distributions have been observed in many nonequilibrium
systems. One of typical classes is the collisional energy cascade, where a
high-temperature particle is injected at a certain rate into an ensemble of
particles, and this energy is distributed to many other particles through
collisions and consumed by the dissipation process in the system. Because of
its simplicity, many systems are categolized into this class. One example is
the granular gasses, which have been modeled with their inelasticity of the
collision Aranson and Tsimring (2006); Villani (2006). With the energy input,
a nontrivial steady state is formed Rouyer and Menon ; van Zon and MacKintosh
(2004); Ben-Naim and Machta (2005); Kang  _et al._ (2010), where the input
energy is balanced with the energy dissipation by the inelastic collisions.
Another example is an ensemble of atoms in plasmas. Energetic atoms are
generated by several processes Corrigan (1965); Hey  _et al._ (2004); Scarlett
_et al._ (2017); McConkey  _et al._ (2008); Starikovskiy (2015) and this input
energy is balanced with collisional energy loss with other particles and walls
having lower temperature.

Despite the simplicity and universal existence of the collisional cascade,
only a few of its statistical properties have been revealed. Ben-Naim et al.
have pointed out that the velocity distribution of granular gasses under the
collisional energy cascade shows a power-law tail Ben-Naim and Machta (2005);
Kang  _et al._ (2010). Although they have reported the analytical
representation of the power-law index, the shape of the overall velocity
distribution including the low-energy part is not clarified. Particularly, the
convergence to the Maxwell distribution in the no-dissipation limit is not
obvious as the power-law index approaches to a finite value in their
representation.

The non-thermal velocity distribution of atoms in plasmas have been observed
experimentally for a long time Vrhovac  _et al._ (1991); Amorim  _et al._
(2000); Samm  _et al._ (1989); Hey  _et al._ (1999); Shikama  _et al._ (2004).
As the existence of the high-energy tail significantly changes the chemical
reaction rates Wakelam  _et al._ (2012, 2015); Flower (1997), the
understandings of the energy distribution are demanded Adamovich  _et al._ .
Many groups have empirically approximated such non-thermal energy
distributions by a sum of a few Maxwell distributions Amorim  _et al._ (2000);
Hey  _et al._ (1999). Some Monte-Carlo simulations have reproduces the
observed non-thermal velocity distribution Sommerer and Kushner (1991);
Starikovskiy (2015); Ponomarev and Aleksandrov (2017). However, the analytical
representation is still missing and thus it is difficult to extract physical
knowledge from the observation.

In this Letter, we point out that the generalized Mittag-Leffler (G-ML)
distribution well approximates the steady-state energy distribution in these
collisional energy cascades. Although the G-ML distribution has no analytic
representation except for a few special cases, its Laplace transform can be
simply written by
$\mathcal{L}_{f_{\operatorname{G-ML}}}(s)=\int_{0}^{\infty}f_{\operatorname{G-ML}}(E)\;e^{-sE}dE=[1+(\epsilon_{0}s)^{\alpha}]^{-\nu}$.
Here, $\epsilon_{0}>0$ is the energy scale, $\nu>0$ is related to the degree
of freedom of the system, and $0<\alpha\leq 1$ is the stability parameter
related to relative importance of the dissipation process. The G-ML
distribution has a power-law tail,
$f_{\operatorname{G-ML}}(E)\approx\nu\alpha\epsilon_{0}^{\alpha}E^{-\alpha-1}/\Gamma(1-\alpha)$,
where $\Gamma(x)=\int_{0}^{\infty}t^{x-1}e^{-t}dt$ is the Gamma function. The
existence of the power-law energy tail is consistent with the previous reports
Ben-Naim and Machta (2005); Kang  _et al._ (2010), but at the same time this
distribution naturally converges to the Maxwell distribution as
$\alpha\rightarrow 1$.

Figure 1: The steady-state energy distribution with the collisional energy
cascade computed by the Monte-Carlo simulation. (a) The energy distributions
of inelastic gas. The upper and lower curves show the results for the Maxwell
gas ($x=0$) and the classical gas ($x=1/2$), respectively. The restitution
coefficients used for the simulations are $e=$ 0.99, 0.9, 0.5, and 0.0. The
best fits by the G-ML distribution are shown by bold curves. The values of
$\alpha$ estimated from the best fit is shown in (c) by stars as a function of
$\overline{\delta}$, which is calculated from Eq. (7). (b) The energy
distributions of the elastic gas colliding with heavier species. The upper and
lower curves show the results for the Maxwell gas ($x=0$) and the Van-der-
Waals gas ($x=1/6$), respectively. Thin solid and dashed curves (they are
almost identical and hardly distinguishable) show the results with different
mass ratio $\mu=10^{-1}$ and $10^{-2}$, respectively. Bold curves are the best
fits by the G-ML distribution. The best-fit $\alpha$ is shown in (c) by
circles as a function of $\overline{\delta}$, which is calculated from Eq.
(8). The theoretical prediction of $\alpha$ (Eq. (5)) is show by dashed curves
in (c). All the estimated values of $\alpha$ are consistent with the theory.

Throughout this Letter, we focus on spatially homogeneous and isotropic
systems. We start from the kinetics of Maxwell-particles in $d$-dimensional
space with no energy dissipation. Let $f(E)$ be the number density of the
particle with the kinetic energy $E$. As we assume the Maxwell particle, the
collision rate among them has no $E$-dependence Maxwell (1867). Therefore, we
can consider the energies of two colliding particles as random samples from
$f(E)$. Let $E_{1}$ and $E_{2}$ be the energies of the colliding particles.
After the collision, the total kinetic energy $E_{1}+E_{2}$ are partitioned
and distributed to the two particles. With the partition ratio $u\in[0,1]$,
which is also a random variable, the post-collision energy of one particle
$E^{\prime}$ can be written as

$\displaystyle E^{\prime}=(E_{1}+E_{2})u.$ (1)

At the steady state, $E^{\prime}$ also follows $f(E)$, i.e., $f(E)$ is one
kind of stable distributions against the operation Eq. (1).

The distribution of $u$, $p(u)$, should be symmetric against the exchange of
$u$ and $(1-u)$. From the statistical weight of the $d$-dimensional space
$\propto E^{d/2-1}$, we expect the beta distribution
$B_{d/2,d/2}(u)=u^{d/2-1}(1-u)^{d/2-1}/B(d/2,d/2)$ as a reasonable choice for
the partition distribution $p(u)$. Here,
$B(a,b)=\int_{0}^{1}x^{a-1}(1-x)^{b-1}dx$ is the beta function. Note that in
this form of $p(u)$ we implicitly assume the no-memory limit, where the
energies of the two colliding particles are completely randomized by the
collision and thus $u$ does not depend on $E_{1}$ nor $E_{2}$. See
Supplemental Material for the correction of this assumption111For the
theoretical and experimental details, see the Supplemental Material, which
includes Kremer (2010); Fantz  _et al._ (2006); Wood (1920); McNeill and Kim
(1982); Baravian  _et al._ (1987); PetroviÄ  _et al._ ; CvetanoviÄ  _et al._
as references.

Because at this point we have no dissipation processes, $f(E)$ should be a
Maxwellian. This is easily seen by taking the Laplace transform of the
distribution, $\mathcal{L}_{f}(s)=\int_{0}^{\infty}f(E)e^{-sE}dE$; the sum of
two random variables (i.e., the convolution) is written as a product of
$\mathcal{L}_{f}$, and the multiplication of the random variables can be
written as a scale mixture. Equation (1) is equivalent with

$\displaystyle\mathcal{L}_{f}(s)=\int_{0}^{1}\bigl{[}\mathcal{L}_{f}(us)\bigr{]}^{2}\;p(u)\;du.$
(2)

Indeed, the Maxwell distribution (the Laplace transform of which is
$(1+\epsilon_{0}s)^{-d/2}$) satisfies Eq. (2).

Let us additionally consider an energy dissipation process. We introduce
another random variable $\delta\geq 0$ that represents the fractional energy
dissipation per one self-collision. The energy of our particle after one self-
collision and the dissipation process is $E^{\prime}=(E_{1}+E_{2})u^{\prime}$
with $u^{\prime}=ue^{-\delta}$. Because of the dissipation process, the
distribution $p(u^{\prime})$ is not symmetric anymore and is skewed toward
$u^{\prime}=0$. We find that for many dissipation processes this skewed
distribution can be approximated by the generalized beta distribution Note
(1),

$\displaystyle
p(u^{\prime})=\frac{\alpha}{B(\nu,\nu)}u^{\prime\alpha\nu-1}(1-u^{\prime\alpha})^{\nu-1},$
(3)

where $0<\alpha\leq 1$ represents the asymmetricity of the distribution.

The distribution stable against the operation
$E^{\prime}=(E_{1}+E_{2})u^{\prime}$ with Eq. (3) is the generalized Mittag-
Leffler (G-ML) distribution, the Laplace transform of which is

$\displaystyle\mathcal{L}_{f_{\operatorname{G-ML}}}(s)=\left[1+(\epsilon_{0}s)^{\alpha}\right]^{-\nu}.$
(4)

Although the G-ML distribution has no analytic representation, the precise
numerical approximation is available Haubold  _et al._ (2011); Barabesi  _et
al._ (2016); Korolev  _et al._ (2020); Note (1). The G-ML distribution has a
power-law tail at large $E$,
$f_{\operatorname{G-ML}}(E)\approx\nu\alpha\epsilon_{0}^{\alpha}E^{-\alpha-1}/\Gamma(1-\alpha)$,
which is consistent with the observation by Ben-Neim et al. Ben-Naim and
Machta (2005); Kang  _et al._ (2010). With the no-dissipation limit
($\alpha\rightarrow 1$ and $\nu\rightarrow d/2$), the G-ML distribution
naturally converges to the Maxwell distribution.

The values of $\alpha$ and $\nu$ may depend on dissipation processes. However,
we find that, in many dissipative systems with $\delta\ll 1$, they are
represented by the following forms Note (1)

$\displaystyle\alpha=\frac{1}{1+\overline{\delta}/\log 2},$ (5)

and $\nu=\frac{d/2}{\alpha^{2}+(1-\alpha)^{2}d/2}$, where $\overline{\delta}$
is the mean of the fractional energy dissipation.

These discussions can be easily extended to more general gasses, where the
collision rate has the energy dependence of $\propto E^{x}$. For example, the
collision rate of the classical balls is proportional to $E^{1/2}$; the
particles interacting through the Van-der-Waals potential collide with
$E^{1/6}$ Massey (1934); Flannery (2006). In order to deal with such systems,
we consider the weighted distribution, $\hat{f}(E)=E^{x}f(E)/Z$, with $Z$ the
normalization constant. Based on an approximation
$(E_{1}/\epsilon_{0}+E_{2}/\epsilon_{0})^{x}\approx(E_{1}E_{2}/\epsilon_{0}^{2})^{x}$,
which is valid if $|x|\ll d/2$, this weighting approximately represents the
energy dependence of the collision rate. Although this weighting changes the
statistical weight of the $d$-dimensional space to $\propto E^{d/2+x-1}$, the
rest of the discussion is the same and we find

$\displaystyle\nu=\frac{d/2+x}{\alpha^{2}+(1-\alpha)^{2}(d/2+x)}.$ (6)

As a numerical demonstration of the above discussion, we carried out several
Monte-Carlo simulations, one of which is similar to that presented by Ben-Neim
et al Ben-Naim and Machta (2005). We prepared $10^{4}$ inelastic particles
with the restitution coefficient $e\in[0,1]$. With a certain rate we heat one
of the particles to the temperature $1$, while this input energy is balanced
with the energy dissipation by the inelastic collisions. A pair of colliding
particles are chosen based on their velocities. Their post-collision
velocities are determined by the differential cross section and the
restitution coefficient Villani (2006); Note (1). Note that in this
simulation, we do not assume the non-memory limit of the collision, but we
consider the actual collision process taking the momentum conservation and the
collision geometry into account. We ran two simulations, one is for the
Maxwell gas ($x=0$) and the other is for the classical gas ($x=1/2$).

In Fig. 1 (a), we show the steady-state energy distributions of these
inelastic particles. The distribution with $e=0.99$, which is close to the
elastic limit, has a similar profile to the Maxwellian in the low-energy
region but has a power-law tail in the high-energy side. With more
inelasticity (smaller $e$ value), the power-law tail becomes larger and the
low-energy distribution deviates from the Maxwellian more significantly. The
bold curves in the figure show the best fits by the G-ML distributions with
taking Eq. (6) into account but with $\alpha$ being left adjustable. All the
results in both the low- to high-energy sides are well represented by the G-ML
distribution.

For the hard-sphere inelastic gasses in 3-dimensional space, we find Note (1)

$\displaystyle\overline{\delta}\approx\frac{1-e^{2}}{8(1+e/3)}.$ (7)

Star markers in Fig. 1 (c) shows the value of $\alpha$ that gives the best fit
to the simulated distribution, as a function of $\overline{\delta}$ which is
computed from the value of $e$ used in the simulation. The values of $\alpha$
at the best fit are consistent with Eq. (5) (the dashed curve).

Another simulation was carried out to study the atom kinetics in plasmas,
where the atoms experience elastic collisions not only among the same species
but also with heavier particles [we can also consider collision to a material
surface instead of the heavy particles]. We assume that the temperature of the
heavier particles is zero and unchanged, i.e., the heavy-particle collision
acts as the energy dissipation. With this setting, the mean of the fractional
energy dissipation $\delta$ can be written as Note (1)

$\displaystyle\overline{\delta}=2^{1-x}\mu\frac{N\sigma^{N}_{\mathrm{mt}}}{n\sigma^{n}_{\mathrm{vi}}},$
(8)

where $\mu=mM/(m+M)^{2}\ll 1$ is the reduced mass ratio ($m$ and $M$ are the
masses of our atoms and heavier particles, respectively), $n$ is the density
of our atoms, $N$ is the density of the heavy particles, and
$\sigma^{n}_{\mathrm{vi}}$ and $\sigma^{N}_{\mathrm{mt}}$ are the viscosity
cross section of the self-collision and the momentum transfer cross section of
the heavy-particle collision, respectively.

Similar to the previous simulation, we prepare $10^{4}$ atoms, generate a high
energy atom with the temperature $1$ at a certain rate, and simulate their
collisional energy cascade. We assume the isotropic-type collision Note (1)
and use the identical cross section for the self-collision and the heavy-
particle collision. Two simulations with the Maxwell gas ($x=0$) and the Van-
der-Waals gas ($x=1/6$) were carried out.

Figure 1 (b) shows the steady-state energy distributions of these atoms in
plasmas. The results with several values of $\mu N/n$ (0.003, 0.03 and 0.3)
are shown in the figure. Also, the simulated distributions under the two $\mu$
values ($10^{-1}$ and $10^{-2}$) but with the same $\mu N/n$ are plotted by
black and gray curves. They are almost identical and hardly distinguished, as
predicted by Eq. (8). As similar to the inelastic-gas simulation, all the
distributions have a power-law tail. With more dissipation, the tail becomes
more significant.

The bold curves in the figure show the best fit by the G-ML distributions with
taking Eq. (6) into account. The G-ML distribution well reproduces the
results. The circle markers in Fig. 1 (c) show the values of $\alpha$ at the
best fit as a function of $\overline{\delta}$, which is computed from Eq. (8)
and the value of $\mu N/n$ used in the simulation. Results with $\mu=10^{-1}$
and $10^{-2}$ are shown by filled and open circles, respectively. All of them
are close to the theoretical prediction, Eq. (5).

Figure 2: Hydrogen Balmer-$\alpha$ spectra observed from the low-pressure glow
discharge tube at different positions. The inset of (a) shows the discharge
tube and the $z$-axis along the capillary. (b) The spectra in the log-log plot
as a function of the wavelength shift $\Delta\lambda$, where $\Delta\lambda$
is defined as the wavelength shift from the nearest fine structure (see the
axis in (a)). The solid curves are the Gaussian-convoluted G-ML distribution
taking the fine structure into account. In the inset of (b), the
$n_{\mathrm{H}}$-dependence of the best-fit $\alpha$ values are shown, where
$n_{\mathrm{H}}$ is estimated from the intensity ratio of the atomic and
molecular lines of hydrogen Note (1).

As an experimental demonstration, we carried out spectroscopic observations of
a low-temperature hydrogen plasma generated inside a commercially available
discharge lamp (Edmund optics, #60-906, shown in the inset of Fig. 2 (a)) with
a DC high-voltage power supply (Spellman SP16122A) Note (1). Hydrogen atoms in
the discharge plasma are mainly generated through electron-impact
dissociation, where the generated atoms typically have 3 eV kinetic energy
Corrigan (1965); Hey  _et al._ (2004); Scarlett  _et al._ (2017). This kinetic
energy is distributed to other atoms by self-collisions and dissipated by the
collision with the walls and other species. We observed the hydrogen atom
Balmer-$\alpha$ line (central wavelength $656.280$ nm) with a homemade high-
resolution spectrometer with the wavelength resolution of 8 pm at $\lambda\sim
650$ nm. The observation was carried at different positions $z$ of the lamp
capillary (see Fig. 2 (a) inset). The small inner diameter of the capillary (1
mm) compared with the mean free path of atoms makes the diffusion coefficient
small. Therefore, the dissociation ratio is higher at deeper position of the
capilarly. This tendency was confirmed by the intensity ratio against the
molecular Fulcher-$\alpha$ line Note (1).

Figure 2 (a) shows the observed Balmer-$\alpha$ spectrum at $z=4$ and 13 mm.
In Fig. 2 (b), we show the same spectrum in the log-log plot as a function of
the wavelength shift. Note that since the Balmer-$\alpha$ mainly consists of
two components, we define the wavelength shift $\Delta\lambda$ from the
nearest line center as shown in Fig. 2 (a). The wing profile is far from the
Gaussian but rather having power-law tails. A steeper tail (i.e., closer to
the Gaussian) is found in the spectrum measured at $z$ = 13 mm than that at
$z$ = 4 mm, while the dissociation ratio at $z=$ 13 mm ($\approx 10^{0}$) is
larger than that in $z=4$ mm ($\approx 10^{-1}$). This spectral profile mostly
reflects the velocity distribution of hydrogen atoms through the Doppler
effect. We note that although the Lorentz component due to the Stark effect
and the diffraction effect inside the spectrometer Fujii  _et al._ (2014),
both of which scale $|\Delta\lambda|^{-2}$ in the tail, are small, we
estimated them from the further wing region of the observed spectra and
subtracted.

The interaction between two hydrogen atoms can be approximated by the Van-der-
Waals potential, thus $x=1/6$ int (1999); Chapman (1991); Note (1). A hydrogen
atom also collides with the inner wall of the capillary and hydrogen
molecules, which should act as the energy dissipation. In our experiment they
should have a finite temperature, though we have assumed the zero temperature
for the heavier species in the above theory. The instrumental broadening
(which is close to the Gaussian profile) also affects the observed profiles.
Therefore, we consider the convolution of the Gaussian distribution and the
velocity distribution derived from the G-ML energy distribution Note (1).

The best fit with the Gaussian-convoluted G-ML distribution is shown by solid
curves in Fig. 2 (b). The observed profiles are well represented by this
distribution up to $|\Delta\lambda/\lambda_{0}|\lesssim 7\times 10^{-4}$,
which corresponds to 9 eV kinetic energy. The values of $\alpha$ at $z$ = 13
mm and 4 mm are 0.83 and 0.74, respectively.

In the inset of Fig. 2 (b), we show the estimated values of $\alpha$ as a
function of the atom density $n_{\mathrm{H}}$, which is estimated from the
dissociation ratio Note (1). Here, we assume that the $\overline{\delta}$ is
inversely proportional to $n_{\mathrm{H}}$. Furthermore, since neither the
absolute atom density nor the restitution coefficient of the wall are clear,
we scaled the horizontal values with a single constant to match the
theoretical prediction Eq. (5). The positive dependence of $1-\alpha$ on
$\overline{\delta}$ is consistent with the theoretical prediction.

Figure 3: Emission spectra of neutral carbon atoms measured for carbon-
nitrogen-argon plasma at the National Solar Observatory NSO . $\lambda_{0}=$
909.483, 1068.308, 1068.535, 1069.125, 1070.734, 1072.953 and 1454.249 nm
lines are plotted with the area normalized to unity. The thin and bold curves
show the best fit by the Gaussian and the Gaussian-convoluted G-ML
distribution, respectively. Eight neutral argon emission spectra are also
shown by gray points, which indicate the negligible effect of the instrumental
profile on the carbon profiles.

This G-ML velocity distribution is also expected for other atoms. Figure 3 (b)
shows the neutral carbon atom emission profiles measured at the Kitt Peak
Solar Observatory NSO  from a microwave discharge generated with the mixture
of carbon, nitrogen, and argon to study cyanide spectra. The spectrum in
$\lambda=$ 830–2710 nm was observed with a high-resolution Fourier transform
spectrometer at the observatory. In the figure, we show seven different
emission lines of neutral carbon measured simultaneously as a function of
$|\Delta\lambda/\lambda_{0}|$. All the profiles show a similar distribution,
indicating the dominant contribution of the Doppler broadening over other
broadening effects, such as Stark broadening or pressure broadening. From the
neutral argon emission lines observed at the same time (gray points in Fig.
3), we find that the instrumental broadening has a negligible effect on the
carbon line profiles. As seen in the figure, the carbon profiles show
significant tails compared with a Gaussian function (thin curve in Fig. 3).
From the fit by the Gaussian-convoluted G-ML distribution (bold curve in the
figure), we obtain $\alpha=0.77\pm 0.02$, suggesting significant self-
collisions in this plasma.

In this Letter, we show that the G-ML distribution universally represents the
steady-state energy distribution of the collisional energy cascade. It is
validated by the numerical simulation for the dissipative gasses, as well as
the spectroscopic observation of the velocity distribution of atoms in
plasmas. Because the stationary parameter $\alpha$ depends on the energy
dissipation ratio, i.e., how close the system is to the thermal equilibrium,
this provides a new diagnostic tool for various systems.

Here, we discussed a fundamental property of spatially homogeneous and
isotropic systems, and thus further investigations are necessary to apply our
findings to non-uniform systems Note (1). However, our assumptions may be
valid if the mean free path of the particles is comparable to or larger than
the system size where the spatial temperature gradient is not significant,
such as the divertor plasmas of nuclear fusion reactors. Although typical
divertor simulation codes assume the Maxwellian for the local velocity
distributions of atoms to compute their transport Kukushkin  _et al._ (2005),
the G-ML distribution can be another reasonable choice for such systems.

###### Acknowledgements.

This work was supported by JSPS KAKENHI Grant Number 19K14680 and 19KK0073. K.
F. thanks Dr. Keiichiro Urabe for the fruitful discussions. K. F. also thanks
an anonimous person with the username vitamin d, who gave us an essential
suggestion in mathoverflow https://mathoverflow.net/questions/401835/inverse-
laplace-transform-of-frac1sa-1-with-0-a-leq-1

## References

  * Aranson and Tsimring (2006) I. S. Aranson and L. S. Tsimring, Reviews of Modern Physics 78, 641 (2006), publisher: American Physical Society.
  * Villani (2006) C. Villani, Journal of Statistical Physics 124, 781 (2006).
  * (3) F. Rouyer and N. Menon,  85, 3676.
  * van Zon and MacKintosh (2004) J. S. van Zon and F. C. MacKintosh, Physical Review Letters 93, 038001 (2004).
  * Ben-Naim and Machta (2005) E. Ben-Naim and J. Machta, Physical Review Letters 94, 138001 (2005), publisher: American Physical Society.
  * Kang  _et al._ (2010) W. Kang, J. Machta, and E. Ben-Naim, EPL (Europhysics Letters) 91, 34002 (2010), arXiv:1002.0995v2 .
  * Corrigan (1965) S. J. B. Corrigan, The Journal of Chemical Physics 43, 4381 (1965).
  * Hey  _et al._ (2004) J. D. Hey, C. C. Chu, P. Mertens, S. Brezinsek, and B. Unterberg, Journal of Physics B: Atomic, Molecular and Optical Physics 37, 2543 (2004).
  * 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).
  * McConkey  _et al._ (2008) J. McConkey, C. Malone, P. Johnson, C. Winstead, V. McKoy, and I. Kanik, Physics Reports 466, 1 (2008).
  * Starikovskiy (2015) A. Y. Starikovskiy, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 373, 20140343 (2015).
  * Vrhovac  _et al._ (1991) S. Vrhovac, S. Radovanov, S. Bzenić, Z. Petrović, and B. Jelenković, Chemical Physics 153, 233 (1991).
  * Amorim  _et al._ (2000) J. Amorim, G. Baravian, and J. Jolly, Journal of Physics D: Applied Physics 33, R51 (2000).
  * Samm  _et al._ (1989) U. Samm, P. Bogen, H. Hartwig, E. Hintz, K. Höthker, Y. Lie, A. Pospieszczyk, D. Rusbüldt, B. Schweer, and Y. Yu, Journal of Nuclear Materials 162-164, 24 (1989).
  * Hey  _et al._ (1999) J. D. Hey, C. C. Chu, and E. Hintz, Journal of Physics B: Atomic, Molecular and Optical Physics 32, 3555 (1999).
  * Shikama  _et al._ (2004) T. Shikama, S. Kado, H. Zushi, A. Iwamae, and S. Tanaka, Physics of Plasmas 11, 4701 (2004).
  * Wakelam  _et al._ (2012) V. Wakelam, E. Herbst, J.-C. Loison, I. W. M. Smith, V. Chandrasekaran, B. Pavone, N. G. Adams, M.-C. Bacchus-Montabonel, A. Bergeat, K. BÃ©roff, V. M. Bierbaum, M. Chabot, A. Dalgarno, E. F. van Dishoeck, A. Faure, W. D. Geppert, D. Gerlich, D. Galli, E. HÃ©brard, F. Hersant, K. M. Hickson, P. Honvault, S. J. Klippenstein, S. Le Picard, G. Nyman, P. Pernot, S. Schlemmer, F. Selsis, I. R. Sims, D. Talbi, J. Tennyson, J. Troe, R. Wester, and L. Wiesenfeld, The Astrophysical Journal Supplement Series 199, 21 (2012).
  * Wakelam  _et al._ (2015) V. Wakelam, J.-C. Loison, E. Herbst, B. Pavone, A. Bergeat, K. BÃ©roff, M. Chabot, A. Faure, D. Galli, W. D. Geppert, D. Gerlich, P. Gratier, N. Harada, K. M. Hickson, P. Honvault, S. J. Klippenstein, S. D. L. Picard, G. Nyman, M. Ruaud, S. Schlemmer, I. R. Sims, D. Talbi, J. Tennyson, and R. Wester, The Astrophysical Journal Supplement Series 217, 20 (2015).
  * Flower (1997) D. R. Flower, Monthly Notices of the Royal Astronomical Society 288, 627 (1997).
  * (20) I. Adamovich, S. D. Baalrud, A. Bogaerts, P. J. Bruggeman, M. Cappelli, V. Colombo, U. Czarnetzki, U. Ebert, J. G. Eden, P. Favia, D. B. Graves, S. Hamaguchi, G. Hieftje, M. Hori, I. D. Kaganovich, U. Kortshagen, M. J. Kushner, N. J. Mason, S. Mazouffre, S. M. Thagard, H.-R. Metelmann, A. Mizuno, E. Moreau, A. B. Murphy, B. A. Niemira, G. S. Oehrlein, Z. L. Petrovic, L. C. Pitchford, Y.-K. Pu, S. Rauf, O. Sakai, S. Samukawa, S. Starikovskaia, J. Tennyson, K. Terashima, M. M. Turner, M. C. M. v. d. Sanden, and A. Vardelle,  50, 323001, publisher: IOP Publishing.
  * Sommerer and Kushner (1991) T. J. Sommerer and M. J. Kushner, Journal of Applied Physics 70, 1240 (1991).
  * Ponomarev and Aleksandrov (2017) A. A. Ponomarev and N. L. Aleksandrov, Plasma Sources Science and Technology 26, 044003 (2017).
  * Maxwell (1867) J. C. Maxwell,  157, 49 (1867), publisher: Royal Society.
  * Note (1) For the theoretical and experimental details, see the Supplemental Material, which includes Kremer (2010); Fantz  _et al._ (2006); Wood (1920); McNeill and Kim (1982); Baravian  _et al._ (1987); PetroviÄ  _et al._ ; CvetanoviÄ  _et al._ as references.
  * Haubold  _et al._ (2011) H. J. Haubold, A. M. Mathai, and R. K. Saxena, Journal of Applied Mathematics , Art. ID 298628, 51 (2011).
  * Barabesi  _et al._ (2016) L. Barabesi, A. Cerasa, A. Cerioli, and D. Perrotta, Electronic Journal of Statistics 10, 3871 (2016), publisher: Institute of Mathematical Statistics and Bernoulli Society.
  * Korolev  _et al._ (2020) V. Korolev, A. Gorshenin, and A. Zeifman, Journal of Mathematical Sciences 246, 503 (2020).
  * Massey (1934) H. S. W. Massey, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 144, 188 (1934).
  * Flannery (2006) M. Flannery, in  _Springer Handbook of Atomic, Molecular, and Optical Physics_ (Springer New York, New York, NY, 2006) pp. 659–691.
  * Fujii  _et al._ (2014) K. Fujii, S. Atsumi, S. Watanabe, T. Shikama, M. Goto, S. Morita, and M. Hasuo, Review of Scientific Instruments 85, 023502 (2014).
  * int (1999) _Atomic and PlasmaâMaterial Interaction Data for Fusion_, Atomic and PlasmaâMaterial Interaction Data for Fusion No. 8 (INTERNATIONAL ATOMIC ENERGY AGENCY, Vienna, 1999).
  * Chapman (1991) S. Chapman,  _The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases_ (Cambridge University Press, Cambridge ; New York, 1991).
  * (33) “NSO’s historical archive,” https://nso.edu/data/historical-archive/, [Online; accessed 01-Feb-2021].
  * Kukushkin  _et al._ (2005) A. S. Kukushkin, H. D. Pacher, V. Kotov, D. Reiter, D. Coster, and G. W. Pacher,  45, 608 (2005), publisher: IOP Publishing.
  * Kremer (2010) G. M. Kremer, in  _An Introduction to the Boltzmann Equation and Transport Processes in Gases_, Interaction of Mechanics and Mathematics, edited by G. M. Kremer (Springer, Berlin, Heidelberg, 2010) pp. 81–107.
  * Fantz  _et al._ (2006) U. Fantz, H. Falter, P. Franzen, D. WÃ¼nderlich, M. Berger, A. Lorenz, W. Kraus, P. McNeely, R. Riedl, and E. Speth, Nuclear Fusion 46, S297 (2006).
  * Wood (1920) R. W. Wood, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 97, 455 (1920), publisher: Royal Society.
  * McNeill and Kim (1982) D. H. McNeill and J. Kim, Physical Review A 25, 2152 (1982).
  * Baravian  _et al._ (1987) G. Baravian, Y. Chouan, A. Ricard, and G. Sultan, Journal of Applied Physics 61, 5249 (1987).
  * (40) Z. L. PetroviÄ, B. M. JelenkoviÄ, and A. V. Phelps,  68, 325, publisher: American Physical Society.
  * (41) N. CvetanoviÄ, B. M. ObradoviÄ, and M. M. Kuraica,  105, 043306.
  * (42) B. P. Lavrov, A. V. Pipa, and J. RÃ¶pcke,  15, 135, publisher: IOP Publishing.
  * Amorim  _et al._ (1996) J. Amorim, G. Baravian, and G. Sultan, Applied Physics Letters 68, 1915 (1996).

Stable Energy Distribution of Weakly Dissipative Gasses
under Collisional Energy Cascades

##  I Theoretical Details

###  I.1 Elastic Gas with Heavy-Particle Collision

In this subsection, we consider the elastic gas undergoing self-collisions and
heavy-particle collisions. The densities of our particle and the heavy
particle are $n$ and $N$, respectively. In the following, we consider the
kinetic energy loss and randomization based on the collision theory. Cross
sections for the self-collision and the heavy-particle collision are
distinguished by the superscripts $n$ and $N$, respectively.

Figure S1:  Schematic illustration of the two-body collision in the center-of-
mass frame. (a) Elastic collision and (b)inelastic collision.

####  I.1.1 Energy loss by heavy-particle collision

Let us consider an atom with mass $m$ moving toward the $z$ axis with velocity
$\mathbf{u}=(0,0,u_{z})$ and elastically colliding with a heavier particle at
rest with mass $M\gg m$. The center of mass (CM) velocity is written as
$\mathbf{V}_{\mathrm{CM}}=\frac{m}{m+M}\mathbf{u}$. With the scattering angle
$\Omega$ in the CM frame, the atom velocity after the elastic collision can be
written as follows:

$\displaystyle v_{x}$ $\displaystyle=(u-V_{\mathrm{CM}})\sin\Omega$ (S1)
$\displaystyle v_{z}$
$\displaystyle=V_{\mathrm{CM}}+(u-V_{\mathrm{CM}})\cos\Omega.$ (S2)

Here, we assume that the scattering takes place in the $x$-$z$ plane without
any loss of generality. The amount of the kinetic energy of our particle lost
by the collision is

$\displaystyle\epsilon$
$\displaystyle=\frac{m}{2}(v_{x}^{2}+v_{z}^{2})=\frac{m}{2}u^{2}\left(1-\frac{2mM}{(m+M)^{2}}(1-\cos\Omega)\right)=E\left(1-2\mu(1-\cos\Omega)\right),$
(S3)

with $\mu=\frac{mM}{(m+M)^{2}}$. By integrating over the scattering angle
$\Omega$, we obtain the average value of the energy loss. In the 3-dimensional
space, this is

$\displaystyle
1-\frac{\overline{\epsilon}}{E}=2\mu\int_{0}^{\pi}(1-\cos\Omega)\frac{\partial\sigma^{N}(\Omega)}{\partial\Omega}2\pi\sin\Omega
d\Omega=2\mu\sigma^{N}_{\mathrm{mt}},$ (S4)

where, $\partial\sigma^{N}(\Omega)/\partial\Omega$ is the differential cross
section for this heavy-particle collision. Here we see that the mean energy
loss can be written with the momentum transfer cross section
$\sigma^{N}_{\mathrm{mt}}$. The list of $\sigma_{\mathrm{mt}}$ for several
types of the collisions can be found in Table S1.

Figure S2:  Cross section the elastic collision between two hydrogen atoms.
(a) Differential cross-section and (b) the energy dependence of the viscosity
cross section.

####  I.1.2 The energy randomization

In Eq. (1), we assume that in one self-collision the kinetic energies of two
collided particles are completely randomized with keeping the total kinetic
energy $E_{1}+E_{2}$. However, the randomization of the kinetic energies is
not generally perfect. For example, the collision with $\Omega\approx\pi$ does
not change their velocity significantly and with $\Omega\approx 0$ only swaps
their energies. As an example of the collision of quantum particles, we show
the differential cross section of elastic collision between two hydrogen atoms
in Fig. S2 (a). The significance of the forward ($\Omega\approx\pi$) and
backward ($\Omega\approx 0$) parts of the cross section is aparent.

In order to quantify the randomization efficiency for each collision, the
viscosity cross section $\sigma_{\mathrm{vi}}$ has been used Kremer (2010)

$\displaystyle\sigma^{n}_{\mathrm{vi}}=\int_{0}^{\pi}\sin^{2}\Omega\;\;\frac{\partial\sigma^{n}(\Omega)}{\partial\Omega}2\pi\sin\Omega
d\Omega,$ (S5)

where $\frac{\partial\sigma^{n}(\Omega)}{\partial\Omega}$ is the differential
cross section for the self-collision. Table S1 shows the value of
$\sigma_{\mathrm{vi}}$ for several types of the collisions. Figure S2 (b)
shows the viscosity cros ssection of hydrogen atom-atom elastic collision. As
the Van-der-Waals interaction is dominant at $E\lesssim 3$ eV, the cross
section behaves as $\propto E^{-1/3}$. Note that in general, the cross section
of the elastic collision of two particles interacting with $r^{-a}$-potential
has the energy dependence of $\propto E^{-a/2}$ Chapman (1991), where $r$ is
the internuclear distance.

####  I.1.3 Average of the fractional energy dissipation $\overline{\delta}$

The mean of the fractional energy dissipation $\overline{\delta}$ for the
elastic gas with heavy particle collision is

$\displaystyle\overline{\delta}=\left(1-\frac{\overline{\epsilon}}{E}\right)\frac{N}{n\sigma^{n}_{\mathrm{vi}}}=2\mu\frac{N\sigma^{N}_{\mathrm{mt}}}{2^{x}n\sigma^{n}_{\mathrm{vi}}}.$
(S6)

Here, we assume that $\sigma^{N}_{\mathrm{mt}}$ and $\sigma^{n}_{\mathrm{vi}}$
have the same energy dependence, $\propto E^{x-1/2}$. The term $2^{x}$ in the
denominator comes from the difference in the relative velocity of the self-
collision ($\propto(2E)^{1/2}$) and the heavy-particle collision ($\propto
E^{1/2}$).

Table S1: A list of crosssections for basic elastic collisions in 3-dimensional space. type |  differential crosssection $\frac{\partial\sigma}{\partial\Omega}$ |  momentum transfer crosssection $\sigma_{\mathrm{mt}}$ |  viscosity crosssection $\sigma_{\mathrm{vi}}$
---|---|---|---
hard sphere | $\frac{\sigma_{\mathrm{el}}}{4\pi}$ | $\sigma_{\mathrm{el}}$ | $\frac{2}{3}\sigma_{\mathrm{el}}$
isotropic | $\frac{\sigma_{\mathrm{el}}}{\pi\sin\Omega}$ | $\sigma_{\mathrm{el}}$ | $\frac{1}{2}\sigma_{\mathrm{el}}$

###  I.2 Inelastic gas

When two particles undergo an inelastic collision, the scattering angle
$\Omega$ depends on the restitution coefficient $0\leq e\leq 1$ (Fig. S1 (b)).
Because of the inelasticity, the momentum normal to the collision direction
(vector $\mathbf{n}$ in the figure) changes Villani (2006),

$\displaystyle v\cos\theta=eu\cos\theta,$ (S7)

with the post-collision velocity $u$ in the CM frame. The momentum
perpendicular to $\mathbf{n}$ is conserved. $u$ and the scattering angle
$\Omega$ can be written as

$\displaystyle u$ $\displaystyle=v\sqrt{\sin^{2}\theta+e^{2}\cos^{2}\theta},$
(S8) $\displaystyle\Omega$
$\displaystyle=\cos^{-1}\left(-\frac{e\cos\theta}{\sqrt{\sin^{2}\theta+e^{2}\cos^{2}\theta}}\right)-\theta.$
(S9)

In this work, we use the hard-sphere cross section for the inelastic gas. By
averaging Eq. (S8), we obtain $\langle u\rangle=v\frac{1-e^{2}}{2}$. Since in
an isotropic system, the kinetic energy should be shared equally by the
kinetic energy in the CM frame and that of the center of mass, i.e., $\langle
v^{2}\rangle=\langle V_{\mathrm{CM}}^{2}\rangle$, the average of the
fractional energy loss per one collision is $(1-e^{2})/4$.

As similar to the discussion for the elastic gas, we use the viscosity cross
section to find the effective collision rate. With the elastic limit
$e\rightarrow 1$, we have $\sigma_{\mathrm{vi}}=2\sigma_{\mathrm{el}}/3$ as
shown in Table S1. Here $\sigma_{\mathrm{el}}$ is the total crosssection. In
the complete inelastic limit
$e=0$,$\sigma_{\mathrm{vi}}=\sigma_{\mathrm{el}}/2$. Although no analytic
forms are available for an arbitrary value of $e$, we adopt the linear
approximation $\sigma_{\mathrm{vi}}\approx 1/2+e/6$, which gives Eq. (7).

###  I.3 Derivation of the generalized beta distribution

Figure S3:  Comparison between Eq. (S11) and Eq. (3). (a) The distributions of
Eq. (S11) and Eq. (3), with $\overline{\delta}=0.03,0.1$, and $0.3$ and $d=3$.
The dashed curves represent Eq. (S11) while the solid curves represent Eq. (3)
with $\alpha$ and $\nu$ determined from Eq. (5) and Eq. (6), respectively. The
gray curve is the distribution for $\overline{\delta}=0$, i.e., the standard
beta distribution, which is symmetric against $u^{\prime}=1/2$. In the inset
of (a), the distribution of $\delta$ (Eq. (S7)) is plotted. (b) Comparison of
their moments. The values of Eq. (S13) are plotted for
$\overline{\delta}=0.03,0.1$ and 0.3. The dashed curves are from Eq. (S11)
while the solid curves are from Eq. (3).

Let us consider the Maxwell atoms experiencing heavy-particle collisions with
$\mu\ll 1$. The time interval $\tau$ between the two successive self-collision
follows $p(\tau)=e^{-\tau/(n\sigma_{\mathrm{vi}})}/n\sigma_{\mathrm{vi}}$.
After one heavy-particle collision, the kinetic energy of atoms changes as a
factor of $\approx e^{-2\mu\sigma_{\mathrm{mt}}^{N}}$. Therefore, the
distribution of the fractional energy loss $\delta\in[0,1]$ between two
successive self-collisions can be written as

$\displaystyle
p(\delta)=\frac{1}{\overline{\delta}}(1-\delta)^{1/\overline{\delta}-1},$
(S10)

with $\overline{\delta}$ is defined by Eq. (S6). The distribution of
$u^{\prime}=u(1-\delta)$ can be derived from the product distribution,

$\displaystyle p(u^{\prime})=\frac{\frac{1}{\overline{\delta}}u^{\prime
1/\overline{\delta}-1}}{B\left(\frac{d}{2},\frac{d}{2}\right)}\left[B\left(\frac{d}{2}-\frac{1}{\overline{\delta}},\frac{d}{2}\right)-B\left(u^{\prime};\frac{d}{2}-\frac{1}{\overline{\delta}},\frac{d}{2}\right)\right],$
(S11)

where $B(x;a,b)=\int_{x}^{1}t^{a-1}(1-t)^{b-1}dt$ is the incomplete beta
function. Dashed curves in Fig. S3 (a) show the distribution by Eq. (S11) for
several values of $\overline{\delta}$. With a larger value of
$\overline{\delta}$, the distribution is skewed more significantly.

The solid curves in Fig. S3 (a) show the generalized beta distribution Eq.
(3). Both distributions are very close to each other. In order to find the
best values of $\alpha$ and $\nu$ that give the closest profile to Eq. (S11)
with given $\overline{\delta}$, we compare the $t$-th moments around
$u^{\prime}=0$ for both the distributions,

$\displaystyle\mathcal{M}_{p(u^{\prime})}(t)=\frac{1}{\overline{\delta}t+1}B\left(\frac{d}{2}+t,\frac{d}{2}\right)\approx
B\left(\frac{\nu+t}{\alpha},\frac{\nu}{\alpha}\right).$ (S12)

We take the logarithm of Eq. (S12) and differentiate by $t$,

$\displaystyle\frac{d}{dt}\left[-\log\mathcal{M}_{p(u^{\prime})}(t)\right]=\frac{1}{t+1/\overline{\delta}}+\psi\left(d+t\right)-\psi\left(\frac{d}{2}+t\right)\approx\frac{1}{\alpha}\psi\left(2\nu+\frac{t}{\alpha}\right)-\frac{1}{\alpha}\psi\left(\nu+\frac{t}{\alpha}\right),$
(S13)

where $\psi(x)=\frac{d}{dx}[\log\Gamma(x)]$ is the digamma function. By
matching the zero-th and first derivatives of Eq. (S13) at $t=0$ with
substituting the approximtion $\psi(x)\approx\log(x)-\frac{1}{2x}$, we obtain
Eq. (5) and Eq. (6). The values of the both sides of Eq. (S13) are shown in
Fig. S3 (b) by dashed and solid curves, respectively. The generalized beta
distribution very well approximates Eq. (S11).

###  I.4 Validity Conditions

In the main text, we adopted the following assumptions to develop our theory,

  1. 1.

The systems are in the steady state.

  2. 2.

The systems are spatially homogeneous and isotropic.

  3. 3.

The collision rate is proportional to $E^{x}$.

  4. 4.

The particle source is located at the high-energy side.

  5. 5.

(for the heavy-particle collision system) The heavy-particle temperature is
zero.

  6. 6.

The fractional energy dissipation is small, i.e., $\delta\ll 1$.

Because these assumptions are relatively strong, the direct application of our
theory to general nonthermal systems is not obvious and further investigations
are necessary.

The spatial homogeneousity (2) may be often difficult to be satisfied, since
in many nonthermal systems, the heating and dissipating sources are spatially
separated from each other. For example, if the heating source is an injection
of the energetic particle to a system, the entrance point of the energetic
particles may be localized at the boundary of the system. In the low-pressure
discharges, such as the one we described in the main text, the dominant
dissipation process is the collision with walls, which are spatially located
at the system boundary. However, even in such inhomogeneous systems if the
mean free path of the particles are comparable with or longer than the system
size, the system can be approximated as homogeneous. This is the case for our
low temperature discharge.

The cross sections of many collision processes have the power dependence on
collision energy (3), but often in the limited energy range. For example, the
hydrogen atom-atom collision behaves as $\propto E^{-1/3}$ at $E\lesssim 3$
eV, as can be seen in Fig. S2 (b). On the other hand, the cross section at
$E\gg 3$ eV has a different energy dependence. Therefore, the G-ML
distribution should be only observed at $E\lesssim 3$ eV for hydrogen atoms,
although the generation of higher energy atoms has been discussed and observed
in a variety of plasmas PetroviÄ  _et al._ ; CvetanoviÄ  _et al._ .

If the heavy-particle energy is much smaller than the typical energy
$\epsilon_{0}$, the zero-temperature assumption (5) may be satisfied. Even if
the heavy-particle temperature is finite, this effect may be corrected by
convoluting the Maxwell distribution to the G-ML distribution, as we did in
the main text.

##  II Numerical evaluation of G-ML distribution

Although closed analytical forms are not available for the G-ML distribution,
a convenient mixture representation has been reported Haubold  _et al._
(2011); Barabesi  _et al._ (2016); Korolev  _et al._ (2020),

$\displaystyle
f_{\operatorname{G-ML}}(E)=\frac{1}{\pi}\frac{\nu^{1/\alpha}}{\epsilon_{0}}\int_{0}^{\infty}\frac{\exp\left(-y{\nu}^{1/\alpha}\frac{E}{\epsilon_{0}}\right)\sin\left(\pi\alpha\nu
F_{\alpha}(y)\right)}{\left(y^{2\alpha}+2y^{\alpha}\cos(\pi\alpha)+1\right)^{2\nu}}dy,$
(S14)

with

$\displaystyle
F_{\alpha}(y)=1-\frac{1}{\pi\alpha}\cot^{-1}\left(\cot(\pi\alpha)+\frac{y^{\alpha}}{\sin(\pi\gamma)}\right).$
(S15)

We evaluate values of the G-ML distribution by numerically integrating Eq.
(S14).

The velocity distribution corresponding to G-ML distribution can be evaluated
by substituting $E=\frac{1}{2m}(v_{x}^{2}+v_{y}^{2}+v_{z}^{2})$ and integrate
it over $v_{x}$ and $v_{y}$ by taking the statistical weight of the space into
account. Since only the term depending on $E$ in Eq. (S14) is $\exp(-cE)$ with
$c=y\nu^{1/\alpha}/\epsilon_{0}$, we only need the integration of this term.

Let us consider the energy distribution $g(E|c)=c\exp(-cE)$ in 3-dimensional
space. Since the statistical weight of the space is $\sqrt{2mE}$, the
distribution of $v_{z}$ is

$\displaystyle g(v_{z}|c)$
$\displaystyle=c\int_{-\infty}^{\infty}\exp(-cE)\frac{1}{\sqrt{2mE}}dv_{x}dv_{y}$
(S16)
$\displaystyle=\frac{1}{2}\sqrt{\frac{c}{2m}}\;\Gamma\left(\frac{1}{2},c\frac{v_{z}^{2}}{2m}\right),$
(S17)

where $\Gamma(s,x)=\int_{x}^{\infty}t^{s-1}e^{-t}dt$ is the lower incomplete
gamma function. By substituting Eq. (S17) into Eq. (S14), we obtain the
velocity distribution of particles with the kinetic energy following the G-ML
distribution.

For $x\neq 0$ case, the energy distribution can be obtained by simply
multiplying $E^{-x}$ to Eq. (S14).

$\displaystyle
f(E)\propto\frac{1}{\pi}\frac{\nu^{1/\alpha}}{\epsilon_{0}}\int_{0}^{\infty}\frac{E^{-x}\exp\left(-y{\nu}^{1/\alpha}\frac{E}{\epsilon_{0}}\right)\sin\left(\pi\alpha\nu
F_{\alpha}(y)\right)}{\left(y^{2\alpha}+2y^{\alpha}\cos(\pi\alpha)+1\right)^{2\nu}}dy.$
(S18)

The corresponding velocity distribution can be obtained by replacing
$g_{x}(E|c)=cE^{-x}\exp(-cE)$ by

$\displaystyle g_{x}(v_{z}|c)$
$\displaystyle=\frac{c^{x}}{2}\sqrt{\frac{c}{2m}}\;\Gamma\left(\frac{1}{2}-x,c\frac{v_{z}^{2}}{2m}\right).$
(S19)

##  III Experimental Details

###  III.1 Spectroscopic Observation of a Commercial Hydrogen Gas Discharge
Tube

Figure S4:  A schematic illustration of the experimenta setup for the hydrogen
plasma measurement. (a) A commercially-available discharge lamp is mounted on
an automatic translational stage. The emission from the lamp is collected by a
pair of achromatic lenses. (b) A schematic illustration of the high-resolution
spectrometer.

Figure S4 shows a schematic illustration of our experimental system to observe
hydrogen velocity distribution. A commercially available hydrogen discharge
lamp (Edmund optics, #60-906) is connected to a DC high-voltage power supply
(Spellman SP16122A). The discharge lamp is mounted on an automated
translational stage (SUS Corporation, XA-28L-50E). The light from the
discharge lamp is collected by a pair of achromatic lenses (focal lengths of
80 mm and 100 mm) and focused on the entrance slit of a homemade spectrometer.

In this spectrometer, the light entering through the slit is collimated by an
achromatic lens (1525 mm focal length, Edmund optics #54-568) and dispersed by
a diffraction grating (2400 grooves/mm, Jovin Ybon). The diffracted light is
focused by another achromatic lens (1525 mm focal length, Edmund optics
#54-568) on the image sensor of a cooled CCD camera (Andor DV435). The
diffraction grating is mounted on an automated rotational stage (OSMS-60YAW,
Sigma Koki) so that the observed wavelength can be chosen remotely. Also, the
pair of achromatic lenses are mounted on the automated translational stage
(HPS60-20x, Sigma Koki) to finely adjust the position and compensate the
residual wavelength dependence of their focal lengths. The resolution was
estimated by an iron emission line as 8 pm at the wavelength $\lambda\sim 650$
nm.

We observed emission from the discharge lamp with 0.8, 1.5, and 3 mA discharge
current and from several different positions. In Fig. S5 (a), we show the
observed intensities of hydrogen atomic Balmer-$\alpha$ line and molecular
Fulcher-$\alpha$ $Q1$ line of the $v^{\prime}-v^{\prime\prime}=2-2$
transition. The edge of the capilarly (1 mm diameter) is at $z=0$ mm, while in
$z<0$ mm the inner diameter of the discharge tube is $\approx 10$ mm. Both the
intensities are almost constant in $z\lesssim 0$ mm, while in $z\gtrsim 0$ mm
the Balmer-$\alpha$ becomes stronger while the Fulcher-$\alpha$ line intensity
remains roughly the same.

Several groups have reported a convenient method to estimate the dissociation
ratio of hydrogen from the emission line intensities Lavrov  _et al._ ; Fantz
_et al._ (2006). We use the result by Lavrov et al. Lavrov  _et al._ , where
the dissociation ratio is estimated from the line intensities of
Balmer-$\alpha$ ($I_{\mathrm{H_{\alpha}}}$) and Fulcher-$\alpha$ $Q1$ line of
the $v^{\prime}-v^{\prime\prime}=2-2$ transition ($I_{\mathrm{ful(2-2)Q1}}$).
Although the value of the electron temperature $T_{\mathrm{e}}$ is necessary
to accurately estimate the dissociation ratio, which is not available in our
discharge, we assume $T_{\mathrm{e}}=10^{0}$ – $10^{1}$ eV, which may be
reasonable for typical glow discharges. Based on these values, the density
ratio can be approximated by $\approx\eta
I_{\mathrm{H_{\alpha}}}/I_{\mathrm{ful(2-2)Q1}}$ with the proportional
coefficient $0.5\times 10^{-2}\lesssim\eta\lesssim 2\times 10^{-2}$.

In Fig. S5 (b), we show the dissociation ratio estimated from the line
intensity. The error bar represented in the figure reflects the uncertainty in
the value of $\eta$. The dissociation ratio starts to increase at the
capillary edge. This tendency can be explained by the small diffusion
coefficient of atoms and molecules in the capillary. Since on the glass
surface the association rate of atoms to molecules is smaller than that on the
metal surface Wood (1920), most of the molecules may be generated on the
electrodes located outside the capillary. On the other hand, the atoms are
mainly generated through electron-impact dissociation inside the capillary.
Therefore, the dissociation ratio is higher in the deeper inside the
capillary. From the pressure balance, the atom density in the discuarge tube
$n_{\mathrm{H}}$ is estimated by $\propto\eta I_{\mathrm{H_{\alpha}}}/(\eta
I_{\mathrm{H_{\alpha}}}+I_{\mathrm{ful(2-2)Q1}})$.

In Fig. S5 (c), we show the value of $\alpha$ that gives the best fit to the
observed Balmer-$\alpha$ profile. This value is almost constant outside the
capillary and decreases deeper inside the capillary. This is consistent with
the behavior of the dissociation ratio, i.e., because of the higher
dissociation ratio in the capillary the rate of the self-collisions of atoms
become more significant.

Figure S5:  Observation result of the hydrogen low-pressure discharge tube.
(a) The observed intensity of hydrogen atomic Balmer-$\alpha$ line (solid
lines) and molecular Fulcher-$\alpha$ $Q1$ line of the
$v^{\prime}-v^{\prime\prime}=2-2$ transition (dotted lines). (b) The estimated
density ratio based on Ref. Lavrov  _et al._

###  III.2 Velocity Distribution of Ground-State Hydrogen Atoms in Plasmas

In the above analysis, we observed the velocity distribution of the excited
hydrogen atoms with assuming that it reflects that of the ground state atoms,
although the direct generation of excited atoms via dissociative excitation
McNeill and Kim (1982); Baravian  _et al._ (1987) may make these two
distributions different. In contrast to our experiment, Amorim et al. have
directly measured the velocity distribution of ground-state hydrogen atoms in
a discharge tube with hydrogen–nitrogen mixture by the two-photon laser
fluorescence method Amorim  _et al._ (2000). We extracted the data from their
electronic manuscript, where the data are embedded as an xml format. Their
spectrum is shown in Fig. S6 on a double-logarithmic scale. The Doppler
profile of the ground-state hydrogen atoms also exhibits a power-law
dependence in its tails. From the fit with by the Gaussian-convoluted G-mL
distribution (the bold curve in Fig. S6), we find $\alpha=0.42\pm 0.02$. From
Eq. (5), a significant dissipation in this plasma compared with the self-
collision is suggested.

Figure S6: The laser-induced fluorescence spectrum of ground-state hydrogen
atoms reported by Amorim et al. Amorim  _et al._ (1996, 2000). The inset shows
the same spectrum in the linear scale. The thin and bold curves show the best
fit by the Gaussian and Gaussian-convoluted G-ML distribution, respectively.

###  III.3 Details of the NSO spectrum

In the main text, we analyzed a high-resolution spectrum by National Solar
Observatory NSO . This was originally measured to study cyanide spectra in
1977 from a microwave discharge with the mixture of carbon, nitrogen, and
argon. This measurement was carried out with a high-resolution Fourier
transform spectrometer with a 1-m optical path difference in the wavelength
range $\lambda=$ 830–2710 nm. The wavelength resolution is 3.5 pm (as the full
width half maximum) for this measurement. Figure S7 (a) shows the original
data and Fig. S7 (b) and (c) show expanded views for argon lines at $\lambda=$
1350, 1372, and 1337 nm and carbon lines at $\lambda=$ 1069.4, 1068.6, and
1454 nm, respectively.

The observed line width of the argon lines is 5.2 pm. This is consistent with
the convolution of the argon Doppler width at the room temperature (3.7 pm)
and the instrumental width. The oscillation seen in the argon line wings,
i.e., the side robe, originates from the instrumental function of this
spectrometer. Since this spectrometer is based on the Fourier transform, the
spectrum is affected by the window function, such as the sinc function
$\frac{\sin\left(\pi\Delta\lambda/\delta_{\lambda}\right)}{\pi\Delta\lambda/\delta_{\lambda}}$
for a rectangular window. Although other line broadenings, such as the Doppler
broadening, averages out this oscillation in the instrumental side robes, this
is still apparent in the argon lines because of their similar Doppler widths
to the instrumental width. However, this effect is negligible for the carbon
lines owing to their much larger Doppler widths.

Figure S7: (a) Spectrum observed in NSO NSO . Markers show the wavelengths of
the argon and carbon lines analyzed in the main text. (b) Expanded spectra of
three argon lines. The light gray markers show the negative points (multiplied
by -1), representing the oscillation of the instrumental side robes. (c)
Expanded spectra of three carbon lines.