AA 2022

# Atomic Hydrogen in the Milky Way: A Stepping Stone in the Evolution of
Galaxies

Naomi M. McClure-Griffiths  1 Snežana Stanimirović2 and Daniel R. Rybarczyk2
1Research School of Astronomy & Astrophysics, Australian National University,
Canberra ACT, Australia 2610 email<EMAIL_ADDRESS>2Department of Astronomy, University of Wisconsin–Madison, Madison WI, USA

###### Abstract

Atomic hydrogen (H i) is a critical stepping stone in the gas evolution cycle
of the interstellar medium (ISM) of the Milky Way. H i traces both the cold,
pre-molecular state before star-formation and the warm, diffuse ISM before and
after star-formation. This review describes new, sensitive H i absorption
surveys, which together with high angular and spectral resolution H i emission
data, have revealed the physical properties of H i, its structure and its
association with magnetic fields. We give an overview of the H i phases,
discuss how H i properties depend on environment and what its structure can
tell us about feedback in the ISM. Key findings include:

  * •

The mass fraction of the cold neutral medium is $\lesssim 40$% on average,
increasing with $A_{V}$ due to the increase of mean gas density.

  * •

The cold disk extends to at least $R\sim 25$ kpc.

  * •

Approximately 40% of the H i is warm with structural characteristics that
derive from feedback events.

  * •

Cold H i is highly filamentary, whereas warm H i is more smoothly distributed.

We summarize future observational and simulation opportunities that can be
used to unravel the 3-D structure of the atomic ISM and the effects of heating
and cooling on H i properties.

###### keywords:

atomic hydrogen (H i), interstellar medium, galaxy evolution, magnetic fields

††journal: Annu. Rev. Astron. Astrophys.

###### Contents

  1. 1 INTRODUCTION
  2. 2 FUNDAMENTALS
    1. 2.1 Overview of phases
    2. 2.2 Observable quantities
  3. 3 NUMERICAL SIMULATIONS OF H i
    1. 3.1 The H i mass fraction as a function of phase
    2. 3.2 Synthetic H i line profiles
  4. 4 OBSERVATIONS OF H i
    1. 4.1 H i in Emission
      1. 4.1.1 Summary of 21-cm Emission Surveys
    2. 4.2 H i Continuum Absorption
      1. 4.2.1 Summary of 21-cm Continuum Absorption Surveys
    3. 4.3 H i Self-Absorption
  5. 5 THE NATURE OF H i IN THE MILKY WAY
    1. 5.1 CNM temperature
      1. 5.1.1 CNM and H2 temperature
    2. 5.2 CNM fraction
    3. 5.3 Environmental Dependencies of Cold H i Properties
      1. 5.3.1 Does the CNM temperature vary across the Milky Way?
      2. 5.3.2 What physical conditions are essential for setting H i properties?
    4. 5.4 The UNM temperature and fraction
    5. 5.5 The WNM temperature and fraction
    6. 5.6 The “Dark” Neutral Medium
  6. 6 CHARACTERIZING THE STRUCTURE OF Hi ON ALL SCALES
    1. 6.1 Deterministic Structure: From Bubbles & Shells to Clouds, Filaments, & Chimneys
    2. 6.2 Linking magnetic fields to H i structure
    3. 6.3 Statistical description of the H i distribution
      1. 6.3.1 H i spatial power spectrum (SPS)
    4. 6.4 3-D Distribution
      1. 6.4.1 Radial Distribution
      2. 6.4.2 Scale Height
  7. 7 FUTURE SURVEYS OF ATOMIC HYDROGEN IN THE MILKY WAY
    1. 7.1 Future complementary data

##  1 INTRODUCTION

Neutral hydrogen (H i) gas is the fundamental building block of galaxies. In
the local Universe fully three-quarters of the neutral gas of galaxies is in
the atomic phase (Carilli & Walter, 2013). H i defines a galaxy’s morphology
on the evolutionary path between the ionized intergalactic medium and the
molecular clouds where stars are formed. H i is observed to have structure on
all spatial scales, from tiny, AU scales observed in the Milky Way (reviewed
in Stanimirović & Zweibel, 2018), to Galactic-scale coherent structures (e.g.
McClure-Griffiths et al., 2004). And, unlike many standard tracers of the
structure of a galaxy, H i is not biased towards tracing star formation - it
follows spiral arms with and without stars, probes gaseous disks extending
beyond the stellar disk (Saintonge & Catinella, 2022) and even traces the
interaction between a galaxy and its starless circumgalactic environment
(Cortese et al., 2021).

Around galaxies, H i exists in the form of highly diffuse, infalling, gas that
replenishes the disk H i reservoir and provides an important energy
contribution. Cosmological simulations predict diffuse, infalling H i with
column density $<10^{18}$ cm-2, and large concentrations originating from
galaxy interactions (Stevens et al., 2019). Around the Milky Way, large
concentrations of H i are associated with the Magellanic System and smaller,
dense concentrations are scattered throughout the halo and disk-halo interface
as High and Intermediate Velocity Clouds. In this review, we will not focus on
the Milky Way’s circumgalactic H i, which is well reviewed by D’Onghia & Fox
(2016) and Putman et al. (2012).

Instead, we will focus on the properties of H i in the disk of the Milky Way
where it is perpetually in the act of “flowing” between its dense, cold state
and its warm, diffuse state. H i is present over a range of temperatures
between $\sim 20~{}{\rm K}$ and $\sim 8000~{}{\rm K}$ (Heiles & Troland,
2003b), tracing most of the semi-quiescent conditions of a galaxy. While the
star formation rate of a galaxy is tightly linked to its molecular hydrogen
(H2) content, the star formation rate per unit H i mass varies widely among
and within galaxies (Leroy et al., 2008). The standard interpretation is that
the star formation in galaxies is separated into two processes: one in which
molecular clouds are assembled from H i and the second in which stars form
from H2 (Schruba et al., 2011). In fact, the key regulator of star formation
appears not to be the formation of H2 itself, but rather the formation of gas
that is sufficiently shielded to reach very low temperatures so that thermal
pressure cannot prevent collapse (Krumholz et al., 2011, Glover & Clark,
2012). H i can therefore be considered the fuel for star formation and the key
ingredient responsible for the shielding and survival of molecular gas. The
transition from warm, unshielded to cold, UV-shielded H i might act as a
throttle in the process of converting the interstellar medium (ISM) into
stars, making it essential to the evolution of galaxies.

In this review we will expand on the role of H i in a galaxy using the Milky
Way as the template. We particularly focus on the properties of the cold
neutral medium (CNM) to highlight the first step in the gas evolution cycle
from diffuse, warm to cold and molecular. Over the last 40 years, H i has been
primarily studied in emission and the focus has been on the large-scale
structure and morphology of the warm neutral medium (WNM). Recent improvements
in observational techniques for H i absorption, together with improved spatial
and spectral resolution revealing cold, small-scale filamentary and fibrous
structures in emission, have redirected Galactic H i studies towards the CNM
as a stepping stone towards molecular cloud and star formation. Whereas the
WNM is smooth and ubiquitous, the CNM has abundant small-scale structure and
is better observed on small scales that are easily accessed in the Milky Way.
Only in the Milky Way can we observe H i with sufficient resolution to probe
the physical scales that are relevant to the cooling, condensing, and
turbulent dissipation that affect how H i converts to H2, as well the feedback
scales on which H i is heated and turbulence driven. Understanding the CNM as
a phase of H i is essential for understanding H i as a path to molecular
clouds and star formation.

The topic of H i in the Milky Way disk is vast and has been periodically
reviewed (e.g., Kulkarni & Heiles, 1987, Burton, 1988, Dickey & Lockman, 1990,
Kalberla & Kerp, 2009). We strive to use this review to emphasize the
importance of H i in understanding how galaxies evolve and what H i can and
cannot tell us about the ISM. We will describe what has been revealed about
the fundamental properties of H i over the last decade, what we still don’t
know, and finally the prospects for the field over the coming decade. While
most of us learn about the hyper-fine transition of atomic hydrogen at
1420.405752 MHz in our undergraduate physics classes, few learn about how it
is measured and the radiative transfer that leads to its widespread usage as a
tracer of mass and dynamics within galaxies. We will start by revisiting some
of the radiative transfer that is needed to understand how physical properties
are derived from observations of H i in emission and absorption. Our vision is
that this article can simultaneously serve as an overview for the new
researcher on the Milky Way and external galaxies, demonstrate the constant
evolution of the field, and convey the excitement of what new telescopes and
simulations will reveal about the Milky Way atomic hydrogen in the coming
decade. For extragalactic researchers this review should complement the review
of the Cold Interstellar Medium of Galaxies in the Local Universe by Saintonge
& Catinella (2022).

##  2 FUNDAMENTALS

Neutral hydrogen, as traced by the hyperfine transition at 21-cm, is usually
described by its spin or excitation temperature, $T_{s}$ which is related to
the H number density of two states of the hyper-fine splitting, $n_{2}$ and
$n_{1}$, by the Boltzmann equation as

$\frac{n_{2}}{n_{1}}=\frac{g_{2}}{g_{1}}e^{-h\nu/kT_{s}},$ (1)

where $g_{2}=3$ and $g_{1}=1$ are the statistical weights for the H i line.
The spin temperature is not necessarily the same as the kinetic temperature,
$T_{k}$ of the gas. Instead, $T_{s}$ is determined by the Lyman $\alpha$
(Ly-$\alpha$) photon field, the pervasive radiation field around the rest
wavelength of the hyperfine transition (such as the CMB), and the kinetic
temperature and density, which set the collisional excitation rate (Field,
1958, Liszt, 2001). Connecting $T_{s}$ and $T_{k}$ requires detailed
understanding of all excitation processes for the 21-cm line. In the high-
density CNM, the 21-cm transition is thermalized by collisions with electrons,
ions, and other H i atoms, which drives $T_{s}$ towards $T_{k}$. In general,
observers assume that $T_{s}\approx T_{k}$ for most observations of
temperatures up to $\sim 1000$ K.

At the low densities typical for the WNM, collisions cannot thermalize the
21-cm transition and therefore $T_{s}<T_{k}$ (Field, 1958, Wouthuysen, 1952,
Deguchi & Watson, 1985, Kulkarni & Heiles, 1988, Liszt, 2001). The Ly-$\alpha$
radiation field can help to thermalize the transition, but a very large
optical depth and many scatterings of Ly-$\alpha$ photons are required to
bring the radiation field and the gas into local thermal equilibrium. While
the underlying atomic physics of the thermalization process is understood
(Wouthuysen, 1952, Field, 1958, Pritchard & Loeb, 2012), the details of
Ly-$\alpha$ radiative transfer are complicated and depend on the topology and
strength of the Ly-$\alpha$ radiation field, which are complex and poorly
constrained in the multi-phase ISM. It is still customary to assume a constant
Ly-$\alpha$ radiation field, with a density of $10^{-6}$ Ly-$\alpha$ photons
${\rm cm^{-3}}$(Liszt, 2001), in theoretical and numerical calculations.

###  2.1 Overview of phases

The balance of heating and cooling processes in the ISM results in a medium
with a range of temperatures and densities. Over a narrow range of pressures H
i can be approximated as a two-phase medium (Field et al., 1969, Wolfire et
al., 1995, Cox, 2005) in which the warm and cold phases of H i coexist in
pressure equilibrium. Of course, the ISM is not entirely atomic; McKee &
Ostriker (1977) proposed the three-phase model in which CNM clouds are
surrounded by the WNM and those neutral clouds are embedded within a supernova
heated hot, ionized medium (Cox & Smith, 1974) \- all in pressure equilibrium
set by the supernovae. Subsequent refinements to the theoretical models
noting, for example, the effects of turbulent pressure have led to a much more
complex structure of phases in the ISM and the understanding that the CNM plus
WNM can occupy as much volume as the hot, ionized medium (e.g. Wolfire et al.,
2003, Bialy & Sternberg, 2019, Audit & Hennebelle, 2005, Hill et al., 2018) .

CNMCold Neutral Medium; $T_{s}=25-250$ K, 28-40% of the observed H i mass
UNMUnstable Neutral Medium; $T_{s}=\sim 250-4000$ K, $\sim 20-28$% of the
observed H i mass WNMWarm Neutral Medium; $T_{s}=4000-8000$ K, $\sim 32-64$ %
of the observed H i mass

The thermal equilibrium range of pressures where the WNM and CNM can both
exist is set by the balance of heating, $\Gamma$, and cooling, $\Lambda$,
leading to the well-known phase diagram showing thermal pressure as a function
of density ($P=T\times\Gamma/\Lambda(T)$). As shown in Figure 1a (black line),
for solar metallicity, the range of stable pressures where both WNM and CNM
can exist is $1500<P_{k}<10^{4}~{}{\rm cm^{-3}~{}K}$, with the upper limit
typically closer to $\sim 3000~{}{\rm cm^{-3}~{}K}$ (Field et al., 1969,
Wolfire et al., 2003). Along the curve where $dP/dn>0$ gas is thermally
stable, meaning that if a gas parcel is perturbed from the equilibrium, it
will either cool or heat to return to the same equilibrium position. These
positions define the WNM at the low density end and the CNM at the high
density end. There is another equilibrium point on the phase diagram where
$dP/dn<0$ that corresponds to the intermediate temperature and density regime.
This region is thermally unstable, meaning that if a H i parcel is perturbed
out of equilibrium it will either cool or heat such that it moves to one of
the two stable solutions. This regime defines the unstable neutral medium
(UNM). Based on theoretical and numerical studies, the UNM temperature is
expected to be in the range $\sim 250-4000$ K, although this depends on the
details of the local heating and coolin (Field et al., 1969, Wolfire et al.,
1995, 2003).

The dominant heating mechanism for diffuse ISM is the photoelectric effect on
dust grains. While the heating rate per unit volume depends on the starlight
intensity and is a function of grain properties, including charging (which
depends on the gas ionization), it does not depend strongly on temperature. As
discussed by Bialy & Sternberg (2019), at very high densities or very low
metallicities H2 heating and cooling becomes important. On the other side of
the balance equation, the key cooling mechanism for cool H i (10–104 K) is the
[CII] 158 $\mu$m fine structure line emission, followed by the [OI] 63 $\mu$m
fine structure line emission. As these transitions are excited by collisions
with H atoms and electrons, the cooling rate depends on the fractional
ionisation. At warmer temperatures ($\gtrsim 10^{4}$ K) Ly-$\alpha$ emission
is the dominant coolant. As these are collisional processes, the cooling rate
per volume depends strongly on temperature, $\Lambda(T)$.

Figure 1: (a): Phase diagram for heating and cooling balance in the ISM. The
solid black line is for solar neighborhood metallicity and UV radiation field.
The orange and blue curves are for UV radiation fields that are 0.1 and 10
times the Solar neighborhood levels. For the multi-phase H i to exist, thermal
pressure is required to be within a well-defined range, $P_{\rm min}$ to
$P_{\rm max}$ denoted by the horizontal lines. Reproduced with permission from
Bialy & Sternberg (2019). (b): The expected CNM kinetic temperature based on
the heating-cooling equilibrium (Bialy & Sternberg, 2019, Bialy, 2020). At
higher optical extinction, H i is shielded from external radiation by dust,
which results in a lower radiation field and reduces the range of CNM
temperatures allowed. Figures provided by S. Bialy.

The details of the classical steady-state heating and cooling processes and
the resulting H i phases were revised in recent modern treatments by Wolfire
et al. (2003) and Bialy & Sternberg (2019). Bialy & Sternberg (2019) explored
a broad parameter space of metallicity, interstellar radiation field, and the
cosmic-rate ionization rate and showed how these effect the thermal structure
of H i. As local conditions of the ISM change, so do the H i phases. For
example, at low metallicity (0.1 Solar) and the case of dust and metal
abundance being reduced by the same amount, the photoelectric heating is less
effective resulting in the expectation that the CNM should be colder. As shown
in Figure 1a (orange curve), when the UV radiation field increases by a factor
of 10, the equilibrium curve shifts to the right and the CNM and WNM co-exist
at a higher pressure.

Other local conditions, such as column density, can also modify the H i
thermal structure. Figure 1(b) shows how varying the column density or optical
extinction $A_{V}$ has a strong effect on the range of predicted CNM
temperatures, $T_{\rm min}$ to $T_{\rm max}$. At low $A_{V}$ the CNM
temperature has a broad distribution ($\sim 70-350$ K), while at $A_{V}>0.2$
the CNM temperature range becomes very narrow, settling around 50 K for
$A_{V}\sim 1$. This is due to dust shielding effectively reducing the ambient
radiation field.

Most of the detailed analytical prescriptions for heating and cooling
processes consider steady-state thermal equilibrium, e.g. McKee & Ostriker
(1977), Wolfire et al. (2003), Bialy & Sternberg (2019). However, supernovae,
stellar winds and many other dynamical processes constantly disrupt the
heating-cooling balance and broaden the pressure distribution function.
Furthermore, turbulent processes constantly drive gas into the unstable gas
regime. To observe how various dynamical sources set the CNM/UNM/WNM
distributions we require numerical simulations.

###  2.2 Observable quantities

The H i line is observed in either emission or absorption, depending on the
optical depth $\tau(v)$ and, by extension, the temperature of the gas,
$T_{s}$. Based on the physics of spectral line profiles (and assuming a
Gaussian line profile), the peak optical depth of the H i line is related to
the gas spin temperature, the (FWHM) linewidth, $\delta v$, and the column
density, $N_{HI}$ as:

$\tau_{0}=\frac{N_{HI}~{}[{\rm cm^{-2}}]}{1.94\times 10^{18}~{}[{\rm
cm^{-2}~{}K^{-1}~{}km^{-1}~{}s}\,]}\,\frac{1}{T_{s}\,[{\rm K}]\,\delta
v\,[{\rm km~{}s^{-1}}]}.$ (2)

For a given column density of H i, $\tau_{0}$ is inversely proportional
$T_{s}$. In fact, because $\delta v$ for a Boltzmann thermal distribution goes
as $\delta v=0.213\sqrt{T_{k}}~{}{\rm km~{}s^{-1}}$, when $T_{s}\approx T_{k}$
then $\tau_{0}\propto T_{s}^{-3/2}$ . If we assume a WNM temperature of
$T_{s}\sim 8000~{}{\rm K}$, which has a thermal velocity linewidth of $\delta
v\sim 19~{}{\rm km~{}s^{-1}}$, we can easily see that for any column density
observed in the Milky Way ($N_{HI}<10^{22}~{}{\rm cm^{-2}}$) the WNM is
optically thin ($\tau_{0}<10^{-2}$) and would require sensitive observations
to be detected in absorption.

A primary observational property of H i in emission is the brightness
temperature, $T_{b}(v)$, which is measured as a function of Doppler velocity
and related to the optical depth, $\tau(v)$, and the spin temperature,
$T_{s}$, as:

$T_{b}(v)=T_{s}(1-e^{-\tau(v)}).$ (3)

Again considering WNM temperatures of 8000 K and if $\tau\ll 1$ the brightness
temperatures are approximated as $T_{b}\approx T_{s}\tau$ and can range from
very small to $\sim 100$ K. Similar arguments about temperature and density
show that H i with non-negligible optical depths, observed via absorption, is
mostly CNM. For most places in the Galaxy, gas is only atomic if it has
$T_{s}>25$ K, implying a thermal linewidth of $\delta v\sim 1~{}{\rm
km~{}s^{-1}}$ and so $\tau>1$ for column densities as low as $1\times
10^{21}~{}{\rm cm^{-2}}$.

It can be tempting to assume that H i absorption mostly traces CNM and that H
i emission traces WNM, but the definitions of phases in emission is not nearly
so clear. If $T_{s}\sim 50$ K, optical depths of $\tau\approx 0.5$ still
produce H i brightness temperatures that are significantly larger than the
Kelvin or sub-Kelvin observational limits. Conversely, observations with very
high optical depth sensitivity can detect the Warm, or at least the Unstable,
Neutral Medium in absorption (Roy et al., 2013a, Murray et al., 2018). The
separation between the WNM and CNM in both emission and absorption therefore
relies on high sensitivity observations and high spectral resolution for
resolving narrow linewidths.

To calculate the column density of a single homogeneous H i feature we need
its excitation temperature and the optical depth:

$N_{HI}=C\,T_{s}\int\tau(v)\,dv,$ (4)

where $C=1.823\times 10^{18}~{}[{\rm cm^{-2}(km~{}s^{-1}~{}K)^{-1}}]$. If
there are several H i features along the line of sight, the column density for
each feature must be treated individually with its own $T_{s}$ (e.g. Heiles &
Troland, 2003a). In the absence of full information about the distribution of
spin temperatures with respect to velocity, it is often necessary and prudent
to assume that the gas is isothermal, such that only one temperature is
represented at each velocity channel, where

$T_{s}(v)\approx\frac{T_{b}(v)}{1-e^{-\tau(v)}}.$ (5)

The column density of the entire sight-line can then be written

$N_{HI,\rm{iso}}=C\int\frac{\tau(v)T_{b}(v)}{(1-e^{-\tau(v)})}dv,$ (6)

as proposed by Dickey & Benson (1982) and developed by Chengalur et al.
(2013). Clearly, for small optical depths Equation 6 reduces to the well known
simplification

$N_{HI,\rm{thin}}=C\,\int T_{b}(v)dv.$ (7)

As described above, H i in emission can trace both the CNM and the WNM. H i
emission on its own provides an estimate of total column density if we assume
optically thin emission ($\tau\ll 1$, Equation 7), and an upper limit on the
kinetic temperature of the gas, $T_{k}$ through its FWHM linewidth, $\delta
v$:

$T_{k}[K]\leq 22\,(\delta v[{\rm km~{}s^{-1}}])^{2}.$ (8)

The CNM temperatures of $\sim 20$ – $200$ K produce thermal linewidths of only
$1$ – $2~{}{\rm km~{}s^{-1}}$, while WNM temperatures of $\sim 3000$ – $8000$
K produce thermal linewidths of $11$-$20~{}{\rm km~{}s^{-1}}$.

H i observers directly measure optical depth $\tau(v)$ through absorption, and
brightness temperature, $T_{b}(v)$ through emission, plus the linewidth,
$\delta v$, of individual features in the spectrum. However, H i emission and
absorption can be combined to solve for the spin temperature, $T_{s}$ of an H
i feature, as described in Equation 3. From here, and by using the linewidth
from H i emission, the contributions from turbulent motions of the gas can be
estimated (Heiles & Troland, 2003b).

##  3 NUMERICAL SIMULATIONS OF H i

Numerical simulations of the ISM are invaluable for understanding how the
many, varied physical processes of the ISM shape H i over a wide range of
environments. Huge advances in computational facilities and techniques in
recent years have moved from 2-D idealized boxes with basic feedback processes
driven artificially at a fixed rate (Audit & Hennebelle, 2005) to 3-D, self-
consistent, parsec-scale, multi-phase simulations that track the effects of
different feedback sources, e.g. SNe, stellar winds, cosmic rays, and
radiation (Kim et al., 2013, Kim & Ostriker, 2017, Rathjen et al., 2021). Even
still, numerical simulations are not able to provide the dynamic range needed
to self-consistently track the H i from large-scale infalling, diffuse gas all
the way to molecular clouds within a larger galactic context (outflows, disk-
halo interface regions, galactic rotation). However, many simulations have
demonstrated the importance of including the large-scale fundamental processes
to shape H i in realistic ways. Hopkins et al. (2012, 2018, 2020), show that
stellar feedback is essential to reproduce realistic H i morphology observed
in the Milky Way and other galaxies, without it simulated gas density
distributions are too high and lead to runway star formation and unrealistic
phase distributions. Similarly, the distribution of spatial power, represented
by the spatial power spectrum, requires stellar feedback to match observations
(Iffrig & Hennebelle, 2017, Grisdale et al., 2018).

Over the last two decades many simulations have focused on high resolution
models of the highly dynamic and turbulent character of the ISM trying to
reproduce the observed thermal and morphological structure of H i. For
example, Audit & Hennebelle (2005) showed that a collision of incoming
turbulent flows can initiate fast condensation from WNM to CNM. Koyama &
Inutsuka (2002) and Mac Low et al. (2005) explored how shocks driven into
warm, magnetized, and turbulent gas by supernova explosions create dense, cold
clouds. The latter study showed a continuum of gas temperatures, with the
fraction of the thermally-unstable H i constrained by the star formation rate.
Ntormousi et al. (2011) modeled the formation of cold clouds in a realistic
environment of colliding superbubbles, and showed that cold, dense and
filamentary structures form naturally in the collision zone. As H i is the
seed for the formation of H2 many simulations start with initial conditions
for the diffuse H i and follow development of the CNM first, and then H2
(Glover et al., 2010, Clark et al., 2012, Valdivia et al., 2016, Vázquez-
Semadeni et al., 2010, Seifried et al., 2022). Simulations such as Vázquez-
Semadeni et al. (2010) and Goldbaum et al. (2011) have suggested that the CNM
plays a central role in the formation and evolution of molecular clouds via
accretion flows. However, exactly how the H i cycles through various phases in
the process of building molecular clouds is still not clear. Numerical
simulations by Dobbs et al. (2012) suggested that it is the UNM that is direct
fuel for molecular clouds.

Figure 2: Fraction of total column density at a given spin temperature as
predicted from the 3-D simulations TIGRESS-NCR (Kim et al., 2022) and SILCC
(Rathjen et al., 2021) and the 2-D simulation by Audit & Hennebelle (2007,
labeled A+H 07). Data were kindly provided to the authors by C.-G. Kim, T.-E.
Rathjen & P. Hennebelle.

A key power of numerical simulations is in providing data for a direct
comparison with observations. We summarize two topics where such comparisons
are important: H i mass fraction and synthetic line profiles.

###  3.1 The H i mass fraction as a function of phase

The mass fraction of different H i phases is a useful quantity to predict
because it can be constrained observationally, which helps to clarify the
processes of importance. The existence of the bi-modal H i distributions has
become a stable feature in recent simulations. While the inclusion of
different feedback sources (supernovae, stellar winds, UV radiation from HII
regions) initially resulted in a broad density distribution without two well-
separated H i phases (e.g., Mac Low et al., 2005), more recent magneto-
hydrodynamic simulations easily produce the bi-stable H i structure with the
CNM and WNM as distinct density peaks (Kim et al., 2014, Kim & Ostriker, 2017,
Hill et al., 2018). Two suites of simulations with the highest resolution that
focus on the thermal state of H i are TIGRESS (Kim & Ostriker, 2017) and SILCC
(Walch et al., 2015, Rathjen et al., 2021). The TIGRESS simulations show that
although turbulence, temporal fluctuations of the heating rate, and expanding
superbubbles continuously populate the UNM phase, thermal instability and
rapid cooling reduce the amount of UNM and result in the majority of the H i
following the thermal equilibrium and being in the two-phase state (Figure 2).
Even when varying the heating rate significantly, Hill et al. (2018) showed
that the two-phase H i persisted and the mass fraction of the thermally
unstable H i was $<20$%.

While recent simulations agree about the existence of the bi-stable H i,
significant differences exist in the mass fraction across different phases.
Figure 2 shows predictions for the H i fraction of total column density at a
given $T_{s}$ from several simulations. SILCC predicts a significantly higher
CNM mass fraction ($\sim 50$%) relative to TIGRESS ($<20$%). The UNM fractions
seem to agree at $\sim 20-30$%. Surprisingly, different feedback channels do
not have a large effect on the H i phase fractions in TIGRESS (Kim, private
communication), while slightly larger differences are seen in SILCC (Rathjen
et al., 2021).

###  3.2 Synthetic H i line profiles

Synthetic H i spectra in emission and absorption are essential for comparing
observations with simulations. Through direct comparison it is possible to
test the reliability of simulations and investigate biases in observational
processing techniques (Haworth et al., 2018). As an example, using TIGRESS Kim
et al. (2014) provided thousands of synthetic spectra that Murray et al.
(2018) analyzed exactly as observed H i spectra. By comparing identified
spectral features in the synthetic spectra with 3D density structures, it was
possible to estimate the completeness and accuracy of the radiative transfer
approach used to observationally estimate $T_{s}$. Similarly, considering
correspondence between the true positions and observed radial velocities of
molecular clouds, Beaumont et al. (2013) showed that the superposition of
clouds along the line of sight introduces significant uncertainty to
observational estimates of cloud mass, size and velocity dispersion. While
this study specifically focused on CO and denser environments, masses and
sizes of H i structures are similarly affected by the line-of-sight
complexities. Simulated spectra will become even more important in the future
to assess the accuracy of automated data handling routines that will be used
by large surveys.

The comparison of simulated and observed spectra is important to fully
constrain the excitation processes for the WNM. A large uncertainty exists in
excitation temperature of H i ($T_{s}$) especially for the WNM where the
collisional excitation is insignificant. Excitation by Ly-$\alpha$ resonant
scattering can be a dominant excitation mechanism (Liszt, 2001, Kim et al.,
2014). Recently, direct Ly-$\alpha$ radiation transfer in a realistic ISM
reveals that the Ly-$\alpha$ excitation in the solar vicinity is efficient
enough to make $T_{s}$ as high as gas kinetic temperature $T_{k}\sim
7000-8000$ K (Seon & Kim, 2020), in agreement with the observational result of
Murray et al. (2014).

##  4 OBSERVATIONS OF H i

Observations of H i have overwhelmingly focused on H i emission, which is used
to study the structure of H i, measure column density with assumptions about
optical depth, and measure velocity linewidths. Recently, H i absorption
surveys have increased in both coverage and sensitivity. However, absorption
observations are targeted towards a limited set of directions but provide
measurements of optical depth as well as linewidth. The power of H i is really
unlocked when absorption and emission are brought together.

###  4.1 H i in Emission

Based on the relations in Section 2.2, the brightness temperature of H i
emission is proportional to column density, resulting in the H i emission
being a powerful tracer of the distribution of H i. Furthermore, it is
possible for modern observational surveys to measure H i emission at all
spatial positions, producing a fully-sampled atlas of the distribution of H i.
Historically, this emission has been used to describe the 3-D distribution of
the Milky Way (Section 6.4), as well as the overall morphology of the atomic
ISM (Section 6). Recently, the spatial structure of H i emission has been
extended to trace the disk-halo interaction (Section 6) and to show its
connection with magnetic fields (Section 6.2). Relying on its complete
sampling of all spatial scales, statistical studies of H i emission have been
used to elucidate the turbulent properties of the medium (see Section 6.3.1).

For many years H i observers have been working to decompose the WNM and CNM
components of the emission on the basis of velocity dispersion (e.g. Kulkarni
& Fich, 1985, Verschuur, 1995). The decomposition is challenging because the
WNM is bright, spatially pervasive and has broad-linewidths, whereas the CNM
has narrow-linewidths and is weak unless $\tau\approx 1$, meaning that its
contribution to $T_{b}(v)$ is easily swamped by the WNM. Recently,
sophisticated techniques for semi-autonomous decomposition have been used on
large-area H i emission surveys to derive estimates of the gas phase fractions
of CNM, UNM and WNM and to consider differences in the structure of the WNM
and CNM (Haud & Kalberla, 2007, Kalberla & Haud, 2018, Marchal et al., 2019).
For example, Kalberla & Haud (2018) used Gaussian decomposition on the HI4PI
survey (HI4PI Collaboration et al., 2016) to show that the CNM is largely
filamentary. Marchal & Miville-Deschênes (2021) used a different decomposition
technique called ROHSA (Marchal et al., 2019) that makes use of spatial
continuity between adjacent spectra to derive a regularized Gaussian fit over
a high-latitude H i field. Murray et al. (2020) took a different approach and
used machine learning with a Convolutional Neutral Network, informed by H i
absorption spectra, to separate the CNM from WNM in emission. These varied
techniques provide estimates of CNM fraction and show structural differences
between the phases.

####  4.1.1 Summary of 21-cm Emission Surveys

H i surveys of emission in the Milky Way have traditionally been carried out
with single dish radio telescopes. With each decade since the discovery of the
H i line in 1951 by Ewen & Purcell (1951) the angular resolution and sky
coverage of H i surveys has increased by roughly a factor of two. Last century
ended with the famous Leiden-Argentine-Bonn (LAB) compilation (Kalberla et
al., 2005) providing a uniform survey at 30′ angular resolution, $1~{}{\rm
km~{}s^{-1}}$ velocity resolution and $\sim 100$ mK brightness sensitivity
over the whole sky. As reviewed in Kalberla & Kerp (2009), the LAB survey
formed the basis for much of what we know about the large-scale H i structure
of the Milky Way.

At the end of the 20th century two technical advances, multi-beam receivers
and interferometric mosaicing, came to radio telescopes and pushed H i surveys
into a new era. The first of these, the multi-beam receiver, decreased the
survey time required to reach high sensitivity over extensive sky areas with
large ($60$ \- $300$ m) telescopes. The multi-beam receivers on the Parkes (13
beams), Effelsberg (7 beams), and Arecibo (7 beams) telescopes all produced a
new generation of surveys covering large areas at resolutions less than 16′.
The Parkes Galactic All-Sky Survey (GASS; McClure-Griffiths et al., 2009,
Kalberla et al., 2010, Kalberla & Haud, 2015) and Effelsberg-Bonn H i survey
(EBHIS; Winkel et al., 2010) were combined together to produce HI4PI (HI4PI
Collaboration et al., 2016) at 16′. The all-sky composite HI4PI shows
structures to be traced across the whole sky and enables statistically
complete studies of phase distributions to be conducted in different regions
with the same data. The most recent, and highest resolution, of the single
dish surveys is GALFA-HI conducted with Arecibo (Peek et al., 2018). GALFA-HI
has pushed single dish H i emission surveys into a new realm of sensitivity
coupled with angular and spectral resolution. Unlike many of the previous
surveys of H i emission, GALFA-HI has a velocity resolution of $\sim 0.2$ km
s-1, revealing very fine velocity gradients across the H i emission sky.

Simultaneously with multi-beam receivers, interferometric mosaicing came into
common practice. Interferometric mosaicing, which combines observations of
many overlapping fields-of-view, enabled observations of a large area of sky
and the recovery of larger angular scales than accessible through a single
pointing. Mosaicing was originally proposed by Ekers & Rots (1979), but was
not widely adopted for diffuse imaging until the 1990’s when Sault (1994)
perfected the technique for the Australia Telescope Compact Array (ATCA)
starting with the Small Magellanic Cloud (Staveley-Smith et al., 1997). The
next development, combining interferometric mosaics with image cubes from
single dish telescopes (Stanimirović et al., 1999, Stanimirović, 2002)
resulted H i cubes that recovered emission on all angular scales from many
degrees to the resolution limit of the array, typically 1 - 3 arcminutes. The
technique was adopted by the Dominion Radio Astrophysical Observatory’s
Synthesis Telescope with the Canadian Galactic Plane Survey (CGPS; Taylor et
al., 2003) and the DRAO H i Intermediate Galactic Latitude Survey (Blagrave et
al., 2017), the ATCA’s Southern Galactic Plane Survey (SGPS; McClure-Griffiths
et al., 2005) and, finally, the Very Large Array (VLA)’s Galactic Plane Survey
(VGPS; Stil et al., 2006).

The time required to reach a given brightness sensitivity goes as the inverse
square of the longest baseline of the interferometer forcing a constant trade-
off between resolution and brightness sensitivity. Because the time needed to
survey large areas at high resolution and with high surface brightness
sensitivity with interferometers is prohibitive, the high angular resolution
interferometric surveys were limited to much smaller aerial regions than
single-dish surveys, focusing primarily on the Galactic Plane and a single
intermediate latitude patch. Recently, the Very Large Array conducted a new
Galactic Plane survey, THOR (Beuther et al., 2016, Wang et al., 2020b), aimed
at improving the angular resolution of the VGPS to 20′′ in angular resolution.
Together, the H i surveys of the last twenty years have demonstrated the
inherent wealth of structure in the Galactic H i over all size scales, which
we will discuss in Section 6 below.

###  4.2 H i Continuum Absorption

The most powerful way to constrain the CNM, with its moderately high optical
depths, is via absorption against continuum background sources. To see how
this proceeds let us extend the radiative transfer discussion from Section
2.2. In general, the simple radiative transfer equation for a single H i
feature given in Equation 3 will include a contribution of the diffuse radio
continuum emission temperature of the sky, $T_{c,sky}$ that incorporates the
CMB and the Galactic synchrotron emission, varying with position. Given a
bright background radio source with a continuum spectrum, $T_{c,bg}$, (usually
a quasar or a radio galaxy) absorption measurements can be made on-source. The
radiative transfer equation in the direction of the background source then
becomes:

$T_{b}(v)=T_{s}(1-e^{-\tau(v)})+(T_{c,sky}+T_{c,bg})e^{-\tau(v)}.$ (9)

Comparing H i spectra obtained in the direction of the background source and
spectra close to the background source (“off-source”, where $T_{c,bg}=0$ K),
under the assumption of uniform foreground H i emission, we can solve Equation
9 for both $T_{s}$ and $\tau(v)$. Equation 4 then provides an estimate of the
column density, $N_{H\textsc{i}}$ of H i structures along the line of sight.

For the uniform H i emission assumption to work, high angular resolution is
required so that the off-source spectrum samples a very similar sight-line to
the on-source spectrum. In this respect, interferometers are highly
advantageous and they resolve out the large-scale structure to help simplify
Equation 9 with $T_{c,sky}\approx 0$ K. Because interferometric mapping with
the inclusion of all spatial scales is time consuming, single-dish telescopes
are often used for the off-source measurements of $T_{b}(v)$. However, as the
emission fluctuations within the single-dish beam can contaminate absorption
measurements, more complex strategies for estimating spatial variations of H i
emission are required (e.g. Heiles & Troland, 2003a). A complicating factor
when using single-dish telescopes to provide off-source spectra is the
mismatch of beam sizes of emission and absorption spectra when solving the
radiative transfer equations for $T_{s}$ and $\tau(v)$. The H i absorption
spectrum samples a very small solid angle occupied by the background (often
point) source, while the H i emission is typically measured on arcminutes
scales (Section 4.1). This can result in absorption features not having clear
corresponding emission components.

The biggest complication in deriving both $T_{s}$ and $\tau(v)$ is that both
$T_{b}(v)$ and $\tau(v)$ spectra usually contain multiple velocity components,
resulting in the need for more sophisticated radiative transfer calculations.
Several different approaches have been used to handle this, e.g. fitting
individual components with Gaussian functions (e.g., Heiles & Troland, 2003a),
the slope method (e.g. Mebold et al., 1997, Dickey et al., 2003) or simply
working with integrated spectral quantities (e.g. Equations 5 and 6). This
problem is exacerbated at low Galactic latitudes where lines of sight have
many velocity components. A related concern is whether selected velocity
structures correspond to real physical structures in the ISM or are possibly
seen as superpositions of many unresolved structures. Such superpositions
introduce biases in observational estimates of $T_{s}$ and $\tau(v)$ and have
been investigated by numerical simulations (e.g., Hennebelle & Audit, 2007,
Kim et al., 2014). Murray et al. (2017) used synthetic spectra from Kim et al.
(2014) and showed that the completeness of recovering $T_{s}$ with Gaussian
fits depends strongly on the complexity of H i spectra. For simulated high-
latitude lines-of-sight ($|b|>50$ deg) 99% of identified structures in radial
velocity had corresponding density features. However, the recovery
completeness dropped to 67% for 20 $<|b|<50$ deg and 53% for 0 $<|b|<20$ deg.
As line-of-sight complexity increases, completeness decreases, reflecting the
difficulty in unambiguously associating spectral features in emission and
absorption in the presence of line blending and turbulence.

####  4.2.1 Summary of 21-cm Continuum Absorption Surveys

Figure 3: Compilation of H i absorption spectra from the bighicat catalog
overlaid on a map of H i emission integrated from -75 km s-1 to +75 km s-1
(HI4PI Collaboration et al., 2016). Non-detections are shown as crosses;
detections are sized according to the number of absorption components.
Coordinates are Galactic.

Due to the high optical depth of the CNM, 21-cm absorption signatures have
been easy to detect, even with low sensitivity observations (e.g. Hughes et
al., 1971, Knapp & Verschuur, 1972, Radhakrishnan & Goss, 1972, Lazareff,
1975, Dickey et al., 1977, 1978, Payne et al., 1978, Crovisier et al., 1978,
Dickey & Benson, 1982, Belfort & Crovisier, 1984, Braun & Walterbos, 1992,
Heiles & Troland, 2003a). Perhaps the most influential single-dish absorption
survey was the Millennium Arecibo 21 Centimeter Absorption-Line Survey,
comprised of 79 H i absorption and emission spectral pairs spread over the
full Arecibo Observatory sky (Heiles & Troland, 2003a).

The optical depth of the UNM and WNM is very low ($\tau\leq 10^{-2.5}$) and
highly sensitive observations are required to detect these phases in
absorption. A few studies have targeted individual detections of the WNM in
absorption to directly measure $T_{s}$ (e.g., Carilli et al., 1998,
Dwarakanath et al., 2002, Murray et al., 2014), while many other surveys
estimated WNM $T_{s}$ from upper limits (i.e. $T_{k}$) or as line-of-sight
averages in the presence of strongly absorbing CNM gas (Mebold et al., 1982,
Heiles & Troland, 2003a, Kanekar et al., 2003, 2011, Roy et al., 2013b). To
clearly detect broad and weak absorption lines from the WNM requires excellent
spectral baselines and is best done using interferometers. The upgrade of the
Karl G. Jansky Very Large Array (JVLA), resulted in the bandpass stability
high enough to detect shallow ($\tau_{0}\sim 10^{-3}$), wide ($\delta v\sim
7-8\rm\,km\,s^{-1}$) absorption lines (Begum et al., 2010). In combination
with H i emission from the Arecibo Observatory, Begum et al. identified
individual absorption lines in the UNM regime with $T_{s}=400-900\rm\,K$. Roy
et al. (2013a,b) used the Westerbork Synthesis Radio Telescope (WSRT), Giant
Metrewave Radio Telescope (GMRT) and the ATCA for another deep H i absorption
survey detecting the UNM. These studies emphasized the need for larger samples
of interferometric detections of H i absorption at high sensitivity to
constrain the fractions of gas in all H i phases.

The 21-SPONGE project (Murray et al., 2014, 2015, 2018) was the one of the
most sensitive absorption line surveys. It targeted 58 bright background
sources reaching sensitivity of $\sigma_{\tau}<1\times 10^{-3}$ to detect in
absorption H i in all stable and thermally unstable phases. The H i absorption
spectra were complemented with H i emission from Arecibo, and a streamlined,
reproducible fitting and radiative transfer approach used to constrain
physical properties of H i.

While targeted H i absorption surveys are time consuming, large
interferometric surveys of the Milky Way plane (CGPS, SGPS, VGPS, THOR) have
had high enough spatial resolution to allow extraction of both H i absorption
and emission spectra in the direction of background sources (Strasser et al.,
2007, Dickey et al., 2003, 2009). For example, this approach enabled a large
sample of $\sim 300$ absorption-emission pairs that was used in Dickey et al.
(2009). While providing important constraints about the spatial distribution
of the CNM, these interferometric surveys mainly used integrated properties as
H i spectra within $\pm$ few degrees from the plane are very complex.

Recent 21-cm H i absorption and emission surveys are summarized in Table 1 of
the Supplementary Material. From a subset of these we have compiled a catalog
containing key observed properties of H i absorption in the Milky Way. The
compilation, which we call bighicat and describe in the Supplementary
material, combines publicly available spectral Gaussian decompositions of
several surveys. In total, bighicat comprises 372 unique lines of sight and
1223 Gaussian absorption components giving $T_{s}$, peak $\tau$ and other
properties. Figure 3 shows the distribution of these lines of sight.

###  4.3 H i Self-Absorption

Another manifestation of H i absorption is H i self-absorption (HISA), which
occurs when cold H i lies in front of warmer background H i emission at the
same velocity. HISA is not, strictly speaking, self-absorption - the two H i
features are not co-spatial, but they are co-spectral. H i in the inner
Galaxy, where velocities correspond to two distances, provides an ideal
background for detecting foreground cold H i as HISA. Broadly, HISA is evident
in two situations: high optical depth environments with a moderate background,
such as molecular clouds (often referred to as HINSA, e.g. Li & Goldsmith
2003, 2012) or moderate optical depth CNM with a high brightness backgrounds,
such as the Galactic Plane (Gibson et al., 2005, e.g.). HISA has historically
been used to study the structure and temperature of cold clouds (e.g.
Heeschen, 1955, Knapp, 1974, Baker & Burton, 1979), but its widespread
presence became particularly clear with the high-resolution afforded by the
interferometric Galactic Plane surveys (Gibson et al., 2000, 2005, Wang et
al., 2020a, c, Beuther et al., 2020).

Figure 4: HISA in the Riegel-Crutcher Cloud (McClure-Griffiths et al., 2006).
The spectra on the right are from the positions marked. The top spectrum is at
the position $\otimes$ and the bottom spectrum, containing HISA, from the
position marked $\oplus$.

An advantage of HISA over other observations of H i emission is that it
clearly separates CNM from WNM. Whereas H i emission is usually dominated by
the WNM, and the physical properties of the CNM are primarily studied in
discrete lines-of-sight towards continuum sources (an HI-absorption “grid”),
HISA can give continuous spatial sampling of the CNM. In practice, however,
the efficacy of spatial recovery of the cloud properties ($N_{H}$, $\tau$, and
$T_{s}$) is not trivial because variations in background H i emission can
obfuscate the properties of the absorbing cold cloud.

Solving the radiative transfer equation (Equation 9) for HISA is even more
complicated than for continuum absorption. We have to consider not just the
continuum background, $T_{c}$ but also the spin temperature and optical depth
of all of the emitting and absorbing H i along the line of sight, including
the foreground ($T_{s,{\rm f}}$, $\tau_{\rm f}$) and background ($T_{s,{\rm
b}}$, $\tau_{\rm b}$), and the HISA feature itself ($T_{s,{\rm HISA}}$,
$\tau_{\rm HISA}$). The expected brightness temperature, $T_{b}$, observed in
the direction of a cold HISA structure is therefore:

$\displaystyle T_{b}\,=$ $\displaystyle T_{s,{\rm f}}\,(1-e^{-\tau_{\rm f}})$
(10) $\displaystyle+$ $\displaystyle\,T_{s,{\rm HISA}}\,(1-e^{-\tau_{\rm
HISA}})\,e^{-\tau_{\rm f}}$ $\displaystyle+$ $\displaystyle\,T_{\rm s,{\rm
b}}\,(1-e^{-\tau_{\rm b}})\,e^{-(\tau_{\rm f}+\tau_{\rm HISA})}$
$\displaystyle+$ $\displaystyle\,T_{c}\,e^{-(\tau_{\rm f}+\tau_{\rm
HISA}+\tau_{\rm b})}.$

All quantities except $T_{c}$ are functions of radial velocity. Unlike
continuum absorption where the brightness temperature without any absorption
can be easily measured adjacent to the source, for spatially extended HISA
$T_{\rm off}$ is usually estimated by interpolating over the absorption
feature in velocity space. In all HISA observations, $T_{s,{\rm HISA}}$ and
$\tau_{\rm HISA}$ are degenerate. Various approaches have been proposed to
solve Equation 10 for $T_{s,{\rm HISA}}$ and $\tau_{\rm HISA}$ (Knapp, 1974,
Gibson et al., 2000), but remains a challenge.

Despite the challenges of robustly solving for $T_{s}$ with HISA, H i self-
absorption provides a valuable probe of the transition zone between atomic and
molecular gas in molecular clouds. Syed et al. (2020) and Wang et al. (2020c),
for example, use the combination of traditional H i emission, HISA-traced gas,
and CO-traced H2 to measure the full hydrogen column density probability
density function for Giant Molecular Filaments. The addition of HISA-traced
gas allows for an optical depth correction to the H i emission, giving a more
realistic total column density (subject to uncertainties in assumed values for
$T_{s}$ and the fraction of foreground to background $T_{b}$). However, recent
synthetic observations of SILCC-Zoom simulations have raised some concerns
about the reliability of using HISA to measure $N_{H}$ (Seifried et al.,
2022). They found that observations of HISA tend to under-estimate the column
density by factors of as much as 3 - 10 because HISA clouds rarely have a
single temperature. Further studies of simulations will help guide observers
in how best to incorporate HISA in $N_{H}$ measurements.

Perhaps the most important, unique information that HISA gives is the
spatially resolved kinematics of the CNM. To spatially resolve velocity
gradients in the CNM from emission observers have to separate the narrow-line
CNM from the pervasive broad-linewidth WNM. Spectral decomposition of complex
H i emission profiles is generally not unique (see Marchal et al., 2019). CNM
traced by HISA, on the other hand, is relatively easy to separate spectrally
from the confusing WNM, making it useful for deriving velocity fields of CNM
structures. Beuther et al. (2020) used this technique to show that the HISA-
traced CNM around the infrared dark cloud G28.3 is kinematically decoupled
from the denser gas as traced by 13CO and [CI]. It seems certain that with
increased angular resolution and sensitivity of future surveys HISA will be
used more extensively.

##  5 THE NATURE OF H i IN THE MILKY WAY

While H i emission surveys have been invaluable in revealing detailed
properties of the Galaxy’s gaseous structure (see reviews by Kulkarni &
Heiles, 1988, Dickey & Lockman, 1990, Kalberla & Kerp, 2009), perhaps the
greatest advances in understanding the physics of Galactic H i over the last
decade have come from studying the physical properties, overall distribution,
and structure of the absorbing H i. The absorbing H i, as we discuss in this
Section, represents mainly the CNM, with a very small fraction of absorption
corresponding to the UNM and WNM. Because of observational limitations, our
understanding of the nature of the CNM of galaxies almost exclusively comes
from the Milky Way. And yet, the CNM is a critical step towards the formation
of molecular clouds and therefore star formation. While many sensitive
observations of H i absorption have been undertaken over the last 20 years, we
are still scratching the surface on understanding CNM properties. Small
samples still limit studies of the thermally unstable H i and a significant
debate persists about what fraction of H i is in the UNM. As we move forward,
the Milky Way will remain the key place for testing details of the multi-phase
physics and we anticipate great advances in this area with upcoming large
surveys.

Following the first astronomical observations of absorption and emission via
the 21 cm transition of H i (Ewen & Purcell, 1951, Muller & Oort, 1951, Hagen
et al., 1955), clear differences in the observed velocity structure between
emission and absorption spectra were attributed to significant variations in
the temperature and density of the gas along the line of sight (Clark, 1965,
Dickey et al., 1978). The distinction between H i structures of different
phase (density and temperature) is still complex. As discussed in Section 4.2,
this complexity comes from the often large overlap of different H i structures
in the radial velocity space, as well as the limitations posed by the
observational sensitivity.

Below we summarize what is currently known about the temperatures and mass
fractions of the different H i phases. We make use of bighicat to examine
environmental dependencies for the CNM.

Figure 5: Fraction of total column density detected as a function of spin
temperature in the bighicat H i absorption compilation. The overall fraction
(black) is divided on the basis on visual extinction, $A_{V}$. The median
temperature is $T_{s}\sim 50$ K. For high extinction environments
($A_{V}>1.0$, orange) the temperature range observed in the CNM is narrow and
peaked at $T_{s}\approx 50$ K, whereas the lower extinction environments
($A_{V}<1.0$, blue) show a broader range of CNM temperatures peaked closer to
$T_{s}\approx 100$ K. The $T_{s}$ cumulative distribution function (CDF) for
each sample is shown on the top. Uncertainties, calculated by bootstrapping
over 1000 trials, are shown as shaded regions around the different histograms.
The sample is incomplete in the range $500~{}{\rm K}\lesssim T_{s}\lesssim
4000$ K where observational sensitivity limits the number of measurements.

###  5.1 CNM temperature

The CNM appears ubiquitous in the Milky Way. As seen in Figure 3, it is easily
detected at low Galactic latitudes, where absorption spectra often have
multiple components, and high Galactic latitudes, where spectra are simple,
often with a single, low optical depth CNM feature. In the bighicat 306/372
directions show significant absorption in at least one component, and many H i
absorption surveys, including 21-SPONGE, have shown a similarly high detection
rate of H i absorption, over 80%.

In Figure 5 we use entries from the bighicat database where individual
velocity components were fitted with Gaussian functions and $T_{s}$ was
derived to show the fraction of the total H i column density as a function of
$T_{s}$ (black line). Murray et al. (2018) used synthetic H i spectra to
demonstrate that the approach of fitting Gaussian functions to individual
velocity components, in combination with radiative transfer calculations, is
successful and complete at recovering the overall fraction of the H i mass
that is in the CNM. The bottom panel of Figure 5 can be compared with the
numerical expectations shown in Figure 2. The observed H i gas as a fraction
of total column density (or mass) peaks around 50-100 K, has a broad shoulder
up to 500 K, and an extended tail up to $T_{s}=3000$ K. Based on theoretical
and numerical studies, the peak of the H i distribution corresponds to the CNM
and has been observationally straightforward to detect (many references listed
in Section 3.2, also Dickey & Lockman, 1990). Clearly, the (absorption)
observed H i gas fraction is missing the WNM portion of the expected
distribution. At this stage the decline in the gas fraction about $T_{s}=500$
K is an observational limitation; even with long integrations it has been very
hard to detect the UNM and WNM in absorption (see Sections 5.4 & 5.5).

The observed peak of CNM temperature distribution of $\sim 50-200$ K is
largely in agreement with theoretical expectations for Solar metallicity and
radiation field (Figure 2), at low optical extinction. The bighicat catalog is
biased towards higher latitude observations (due to simplicity of H i spectral
lines). As a result, about $41\%$ of the bighicat sample probes environments
with $A_{V}<0.2$. In Figure 5 we also show gas fraction histograms by
splitting the sample into $A_{V}<1$ and $A_{V}>1$ (shown in blue and orange
lines, respectively). The $A_{V}<1$ sub-sample peaks at 80-100 K, while the
$A_{V}>1$ sub-sample peaks at a lower temperature ($\sim 50$ K). This figure
suggests that the broad CNM distribution comes from a combination of denser
and more diffuse environments, with different physical environments resulting
in the different range of temperatures.

Overall, we conclude that the observed CNM $T_{s}$ distribution is in line
with theoretical expectations. It is broad, but the width appears in line with
expectations for the environmental dependence of $T_{s}$, as well as how local
turbulent perturbations affect H i temperature. Building larger CNM samples to
probe even more diverse ISM environments is highly important for the future.

####  5.1.1 CNM and H2 temperature

Because it is expected that molecular hydrogen (H2) forms largely out of the
CNM (e.g., Krumholz et al., 2009) an important question is: how do the CNM and
H2 temperatures compare? For example, Heiles & Troland (2003b) compared the
CNM temperature with the H2 temperature measured by FUSE for four directions
where radio and UV sources were spatially close. Their conclusion was that H2
and CNM temperatures did not agree. If we perform a similar analysis using the
sightlines in the bighicat within $\sim 1^{\circ}$ of FUSE sightlines in Shull
et al. (2021), we find that the H2 temperature, as measured through the lowest
three rotational states (J = 0, 1, 2) of H2, always lies between the minimum
and maximum CNM temperature along the line of sight. In 3/4 cases, the H2
temperature is consistent with the temperature of the CNM component with the
highest H i column density. Yet, as Heiles & Troland (2003b) noted, the
temperatures derived for H2 and H i are not directly comparable. The H2
temperatures are calculated as a weighted mean over all velocity components
because of the strong line saturation.

The mean excitation temperatures in Shull et al. (2021) are $\langle
T_{01}\rangle=88\pm 20$ K and $\langle T_{02}\rangle=77\pm 18$ K, which is
very close to the CNM peak seen in Figure 5. For sight lines with $E(B-V)>0.5$
(corresponds to $A_{V}>1.5$) and $N_{H}>10^{20.7}$ cm -2, the H2 temperatures
decreased to 50–70 K. This is qualitatively in agreement with Figure 5 for
$A_{V}>1$ where we see that $T_{s}$ shifts to lower temperature, supporting
the idea that H2 and the CNM are thermally coupled. Detailed comparisons of
larger samples of H i velocity components with H2 measurements remain as an
important future task. Homogeneous samples of other molecular species found in
the diffuse ISM are starting to emerge thanks to ALMA and NOEMA. For example,
Rybarczyk et al. (2022) found that the molecular gas traced by HCO+, C2H, HCN,
and HNC is associated only with H i structures that have an HI optical depth
$>0.1$, a spin temperature $<80$ K, and a turbulent Mach number $>2$. This
result hints that not all CNM is useful for forming H2, with only colder and
denser CNM being thermally coupled to H2.

Figure 6: Left: line-of-sight mean spin temperature, $\langle T_{s}\rangle$,
reproduced from data in Dickey et al. (2009). The mean spin temperature is
remarkably constant with Galactic radius in the outer Galaxy. Right: the
implied CNM fraction $f_{\rm CNM}=T_{s}/\langle T_{s}\rangle$ as a function of
galactic radius. The CNM fraction shown here is dependent on the assumed cold
gas temperature. The dark points assume $T_{s}$ varies with radius based on
Wolfire et al. (2003) (described in the text) and the lighter points assume
$T_{s}=50$ K. The VGPS data are affected by the Galactic warp, which takes the
cold disk out of the survey area.

###  5.2 CNM fraction

As the CNM is the key building block for H2, constraining how the CNM fraction
varies across the Milky Way is essential for understanding star formation
efficiencies. The Millennium Arecibo survey showed that the majority of 79
observed directions had a CNM column density (or mass) fraction $<0.3$, with
only a few lines of sight being dominated by the CNM (Heiles & Troland,
2003b). In a few directions, this study found essentially no CNM and suggested
that the CNM was destroyed in these directions by superbubbles. They also
noticed that the CNM fraction increased with the total H i column density up
to $\sim 1.2\times 10^{20}$ cm-2, and then leveled off. By comparing the
WNM/CNM fraction at $|b|<10$ degrees, where lines of sight are mainly in the
Galactic disk, with higher latitudes, Heiles & Troland (2003b) tested the
hypothesis that the WNM fraction should be lower in the plane where thermal
pressure is the highest but found no difference.

With larger samples and more sensitive data, observers are starting to find
regions with CNM fractions that clearly deviate from typical ISM values.
Murray et al. (2021) showed that the Galactic H i in the foreground of Complex
C is particularly under-abundant in the CNM. The whole region has a relatively
low column density ($<10^{20}$ cm-2 on average) and in spite of narrow H i
linewidths shows a very low CNM fraction. Although the line-of-sight CNM
fraction in the region is different from typical ISM fields, individual
velocity components have similar $T_{s}$ and $\tau_{0}$, suggesting that this
region is particularly quiescent. Murray et al. (2021) concluded that the
region may be sampling an area that has not been recently disturbed by
supernova shocks, leading its H i properties to be dominated by thermal
motions rather than nonthermal, turbulent motions.

At the other extreme, observations near several GMCs have shown on average
higher CNM fractions. For example, Stanimirović et al. (2014) investigated 26
H i absorption and emission pairs obtained in the direction of radio sources
all in the vicinity of Perseus molecular cloud. Their CNM fraction is in the
range 0.2-0.5, with the median CNM fraction being $0.33$ (in comparison to
$0.22$ from Heiles & Troland 2003b). Similarly, Nguyen et al. (2019) used 77 H
i emission-absorption pairs in the vicinity of Taurus, California, Rosette,
Mon OB1, and NGC 2264 and found the CNM fraction in the range $\sim 0.2-0.75$
with the median value of $0.35$. These results suggest a scenario in which a
high CNM fraction is required for molecule formation, and GMCs are built up
stage by stage—from WNM-rich gas to CNM-rich gas to molecular clouds. Another
reason for the high CNM fraction around GMCs could be effective CNM accretion.
An important future task will be to provide even denser grids of H i
absorption sources and map spatially and kinematically the distribution of the
CNM in the vicinity of GMCs.

The line-of-sight mean spin temperature $\langle T_{s}\rangle$ can also be
used as an indicator of the CNM fraction where $f_{\rm CNM}=T_{s}/\langle
T_{s}\rangle$, assuming $T_{s}$ as an average cloud temperature. Dickey et al.
(2009) compiled H i absorption-emission measurements in the Galactic plane
from the three Galactic Plane Surveys (CGPS, SGPS and VGPS) to measure
$\langle T_{s}\rangle$ as a function of Galactic radius, $R$. We have
reproduced those data in Figure 6. We also show in Figure 6 the implied CNM
fraction $f_{\rm CNM}$ with two representative assumptions of $T_{s}$: average
cold temperatures from Wolfire et al. (2003) heating and cooling models
($T_{s}=85$ K at $R=8.5$ kpc to $T_{s}=44.1$ K at $R=18$ kpc) or a constant
temperature of $T_{s}=50~{}{\rm K}$. If the typical cold cloud temperature
decreases with Galactic radius then $f_{\rm CNM}$ also decreases from $\sim
0.3$ near the Solar circle to less than $0.05$ at $R\sim 25$ kpc. On the other
hand, if we assume a flat $T_{s}$, as indicated by Strasser et al. (2007), the
CNM fraction is remarkably flat to large radii. Clearly the key observational
element required to understand how the CNM fraction varies will be the typical
cold cloud temperature. For completeness we note that although these curves
were calculated under the assumption of smooth circular rotation, which breaks
down for individual sight-lines, the trends largely hold for the extended
azimuthal ($>70\mbox{${}^{\circ}$}$) averages of the data shown here.

H i emission surveys have also been used, together with spectral
decomposition, to determine the fraction of gas in the CNM. For example,
Kalberla & Haud (2018) found that for Galactic latitudes
$|b|>20\mbox{${}^{\circ}$}$ the column density fraction of CNM was 25$\pm 5$%
at local velocities, whereas Marchal & Miville-Deschênes (2021) estimated that
in their high latitude field the average CNM mass fraction was $8\pm 6$%,
showing large excursions (up to $\sim 30$%) along filamentary structures.
Using a convolutional neutral network applied to the large-area survey GALFA-
HI Murray et al. (2020) found that the CNM fraction varies from $<10$% to
$\sim 40$%, depending on sky location.

In the full bighicat compilation the mean CNM fraction (of $T_{s}<250$ K)
along the line-of-sight is $0.35$ and the median is $0.34$. However, the
bighicat compilation is not uniform and contains numerous fields selected to
be near the Galactic plane or molecular clouds. Limiting the sample to
$|b|>20\mbox{${}^{\circ}$}$ gives a mean of $0.31$ and a median of $0.27$ and
might be more representative of the general ISM. The bighicat CNM mass
fraction is 40% if we assume that 50% of the total H i mass is detected only
in emission (Murray et al., 2018).

###  5.3 Environmental Dependencies of Cold H i Properties

Theoretical models and numerical simulations show that the thermal pressure at
which cold and warm H i co-exist depend on key physical properties such as the
ambient interstellar radiation field, metallicity, dust properties
(composition, charge, efficiency etc), and interstellar turbulence. As these
physical quantities vary across galaxies, we expect to see variations of H i
gas properties. For example, photoelectric heating by dust grains may be
enhanced in particularly dust-rich environments. Large-scale gradual
metallicity variations are common across spiral galaxies (e.g. Shaver et al.,
1983). Even small-scale metallicity variations could be common as suggested by
De Cia et al. (2021) who showed that many regions in the Solar vicinity have
low metallicity, down to about 17 per cent Solar metallicity and possibly
below. Similarly, as shown in Figure 1, the CNM temperature should be affected
by dust shielding in higher $A_{V}$ regions. However, besides several sporadic
studies (some listed in the previous sections), the observational evidence of
the regional diversity of CNM properties across the Milky Way is still
lacking. This is largely due to the limited number of H i absorption spectra,
as well as large measurement uncertainties in deriving key physical properties
($T_{s}$, $f_{\rm CNM}$). With the systematic observational and data analyses
approaches, and larger data samples in recent years, some interesting trends
are starting to emerge.

####  5.3.1 Does the CNM temperature vary across the Milky Way?

Over the years, sporadic observations have shown occasional directions with
$T_{s}\sim 20-30$ K (Knapp & Verschuur, 1972, Meyer et al., 2006, 2012). These
temperatures are very low and in the realm of what is usually found for
molecular gas. Explaining theoretically such low temperatures requires the
absence of photoelectric heating (e.g., Spitzer Jr, 1978, Wolfire et al.,
2003), and/or significant shielding (Glover et al., 2010, Gong et al., 2017).
On the other hand, many studies have noticed relatively uniform $T_{s}$ for
the CNM in the inner and outer Galaxy (Strasser et al., 2007), in and around
the Perseus molecular cloud and several additional GMCs (Stanimirović et al.,
2014, Nguyen et al., 2019). Using line-of-sight integrated properties, Dénes
et al. (2018) hinted at the existence of a slightly warmer $T_{s}$ close to
the Galactic center, as did (Bihr et al., 2015) near the W43 “mini-starburst”
star-forming region.

Theories of heating and cooling balance (e.g. Figure 1) show that regions with
high dust column density (high $A_{V}$) have a low CNM temperature, whereas in
low $A_{V}$ regions the CNM temperature takes on a wide range of values (Bialy
& Sternberg, 2019, Bialy, 2020). We use all available H i absorption spectra
compiled in the bighicat dataset to test whether $T_{s}$ depends on the
optical extinction. As shown in Figure 5, the H i components in the directions
with $A_{V}<1$ occupy the higher-$T_{s}$ portion of the distribution relative
to H i components in the directions with $A_{V}>1$ which predominantly sample
the lower-$T_{s}$ portion of the $T_{s}$ distribution. This demonstrates that
higher $A_{V}$ regions provide more intense shielding, which result in lower
$T_{s}$ (see also Figure 10 in Kim et al. 2022). While this result is in line
with theoretical expectations the observed distributions of the cold H i are
broader than theoretical expectations. This could be due to the presence of
turbulence which spreads out the temperature distribution.

Clearly, these observations are just the beginning and future H i surveys with
larger samples of absorption lines are needed to confirm these trends. These
results demonstrates the power of large H i absorption samples to, for the
first time, depict differences of H i temperature distributions caused by the
underlying physical conditions.

Figure 7: Left: The optical-depth weighted spin temperature from the bighicat
datasets as a function of optical extinction from Planck. Points are coloured
by the dust temperature and lines are from the PDR model by Gong et al.
(2017). Right: The CNM fraction from the bighicat dataset. The points are
colored by mean volume density $\langle n_{\rm{H}}\rangle$, extracted from the
3D Lallement et al. (2019) dust map. The H i absorption data probe arcsecond
angular scales, while $A_{V}$ from Planck probes angular scales of $\sim
5^{\prime}$. The high-$A_{V}$ limit in this figure is not well defined as H i
spectra can be saturated and complex for interpretation.

####  5.3.2 What physical conditions are essential for setting H i properties?

The $T_{s}$ distribution in Figure 5 suggests that the H i temperature depends
on optical extinction. To investigate a physical origin of the dependence of
temperature on extinction in Figure 7 (left) we show the optical depth-
weighted average spin temperature, $\langle
T_{s}\rangle=\int\tau(v)T_{s}(v)dv/\int\tau(v)dv$, as a function of $A_{V}$.
Both $\langle T_{s}\rangle$ and $A_{V}$ are line-of-sight quantities. We see
that cold H i exists over a very broad range of conditions: from very diffuse
environments $A_{V}\sim 0.03$ all the way to $A_{V}\sim 10$, which is often
used to characterize dense molecular clouds (Snow & McCall, 2006). However,
the observed distribution is smooth, suggesting a continuous range of H i
structures, instead of well-segregated groups (e.g. diffuse atomic, diffuse
molecular, translucent, dense molecular). We overplot in this figure two color
lines showing predictions from the PDR model from Gong et al. (2017) for
particular input values $G\sim 1$ (i.e. solar UV radiation field) and $n=1-10$
cm-3. These lines nicely bracket the observed $\langle T_{s}\rangle$ values
showing that the local density plays a crucial role in constraining H i
properties.

As discussed above, $\langle T_{s}\rangle$ is related to the CNM fraction
along the line-of-sight, $f_{\rm CNM}$111The values of $f_{\rm CNM}$ shown in
Figure 7 are calculated using Equation 18 of Kim et al. (2014). Uncertainties
are calculated by varying the warm and cold cloud temperatures, as in Murray
et al. (2021). We find the same qualitative results if we take $f_{\rm{CNM}}$
as $N|_{T_{S}<250~{}\rm{K}}/N_{\rm{LOS}}$. as $\langle T_{s}\rangle\sim
T_{s}/f_{\rm CNM}$ so the trend of $\langle T_{s}\rangle$ shown in Figure 7
(left) can reflect variations in the cold cloud temperature, $T_{s}$, or the
CNM fraction, $f_{\rm CNM}$. However, for $A_{V}>0.2$ (see Figure 1) the CNM
temperature has a very narrow distribution ($\sim 50$ K). Assuming a well
constrained CNM temperature for $A_{V}>0.2$, Figure 7 (right) shows that
$f_{\rm CNM}$ increases with $A_{V}$ from $\sim 10$ to $\sim 80$%. The color
of the points in the right panel of Figure 7 shows the mean volume density
along the line of sight, showing that the sightlines with higher $A_{V}$—and
higher $f_{\rm{CNM}}$—are those with the highest mean densities. We conclude
that local density plays a crucial role in constraining the CNM fraction.
Higher local density implies more shielding, as well as cooling, and therefore
more CNM.

While some particular conditions could result in the H i being entirely CNM,
this is not the case in the Milky Way. The mean points in the right panel of
Figure 7 (shown as black squares), show that the CNM fraction flattens at
about 60-80%. Some of the highest-$A_{V}$ directions certainly probe
environments around dense molecular clouds so the flattening of the CNM
fraction at $A_{V}\sim 1$ could signify the transition of H i into H2.
However, as Stanimirović et al. (2014) pointed out, even lines of sight that
probe deep inside the GMCs have a contribution from the WNM. The geometry and
the level of mixing of the CNM and WNM are still not understood; it is not
clear if the WNM is located primarily in outer regions of the H i envelope as
theorized by McKee & Ostriker (1977) or if WNM is well-mixed by turbulence. A
priori one might assume that WNM mixed into high-pressure molecular clouds
would rapidly cool, but Hennebelle & Inutsuka (2006) showed that the
dissipation of magnetic waves can provide substantial heating to maintain the
WNM inside molecular clouds. It therefore seems unsurprising that the CNM and
WNM are well mixed in all environments and $f_{\rm CNM}<1$.

The low-$A_{V}$ end of Figure 7 is also interesting, but poorly constrained by
the observational data. Kanekar et al. (2011) noticed that sightlines with
column density $>2\times 10^{20}$ cm-2 had low optical-depth averaged spin
temperature, which they interpreted as a evidence of a threshold for CNM
formation. Kim et al. (2014) suggested that instead of a column density
threshold, $A_{V}\sim 0.1$ may simply correspond to a characteristic length
scale for CNM structures. In spite of the larger sample in bighicat, the jury
is still out on whether column density threshold exists for CNM formation.

###  5.4 The UNM temperature and fraction

It has been difficult to constrain the observational properties of H i at
temperatures $T_{s}>250~{}{\rm K}$. Until recently most estimates of the UNM
temperature were made as upper limits from the line-width-based kinetic
temperatures of either emission or absorption. UNM mass fraction estimates
made from H i emission give fractions between $28\pm 11\%$ (Marchal & Miville-
Deschênes, 2021) and $41$% (Kalberla & Haud, 2018). Using upper limits to
kinetic temperatures of H i absorption, Kanekar et al. (2003) found UNM
fractions of: 77% and 72% by mass for two LOS. By contrast, Begum et al.
(2010) estimated $27\%$ across five lines-of-sight and Roy et al. (2013b) set
a lower limit of $\sim 28$% mass fraction with values as high as $44\%$. The
lower limit is based on the conservative assumptions that all detected CNM has
$T_{s}=200\rm\,K$ and all non-detected WNM has $T_{s}=5000\rm\,K$, and not on
direct measurements of the absorbing and emitting properties of the unstable
gas. Limits on the fraction of the medium in the unstable state based on spin
temperature should be more robust, but have been difficult because of the low
optical depths inherent to the UNM and WNM.

Only a handful of observations have achieved the optical depth sensitivity
required to measure $T_{s}$ in the UNM, e.g., Heiles & Troland (2003b)
detected 13 UNM components giving an estimated column density fraction of 29%
for spin temperatures measured out of the plane. The best estimates for UNM
fraction have come from the 21-SPONGE survey (Murray et al., 2018), which
conducted very deep H i absorption observations and was sensitive to H i with
spin temperature up to $\sim 4000$ K. Not only this was the first
statistically significant survey designed to detect low optical depth
($\tau\sim\mathrm{few}\times 10^{-3}$) broad lines, but this survey also took
exceptional care of observational biases introduced in the analysis method
used to constrain H i spin temperature and recover the H i mass distribution.
21-SPONGE tripled the number of UNM components and found $20\%$ of the total H
i mass, and $41\pm 10$% of the H i mass detected in absorption, to correspond
to the UNM (Murray et al., 2018). The result agreed with the Heiles & Troland
(2003b) sample and the Begum et al. (2010) and Roy et al. (2013b) studies. The
high sensitivity Nguyen et al. (2019) survey subsequently found a similar
fraction.

A seemingly trivial, but important, distinction in the various observational
and theoretical studies is the chosen boundary between the CNM and UNM. By
combining several H i absorption surveys in bighicat, we have the largest
database of high latitude ($b>15\mbox{${}^{\circ}$}$) spin temperatures to
search for a possible observational boundary between CNM and UNM. We find no
clear observational boundary, but suggest that observers and theorists adopt
the value of $250$ K as an nominal boundary between CNM and UNM simply to ease
cross-comparison.

Driven largely by 21-SPONGE, the consensus seems to be pointing towards about
20 - 30% of the H i by mass being in the unstable phase. The result is similar
to the H i emission derived UNM fractions (28% and 41%) by Marchal & Miville-
Deschênes (2021) and Kalberla & Haud (2018) given the errors on those
estimates. The observed UNM fraction agrees well with the mass fraction at UNM
temperatures from the TIGRESS and SILCC numerical simulations shown in Figure
2 (Kim & Ostriker, 2017, Rathjen et al., 2021).

The spatial distribution of the UNM is currently uncertain. For example,
Murray et al. (2018) noticed that highest Galactic latitudes are dominated by
WNM, while the CNM and UNM dominate low Galactic latitudes. However, the
Gaussian decomposition is least reliable at low latitudes. On the other hand,
Kalberla & Haud (2018) who studied UNM via H i emission lines, concluded that
CNM structures are surrounded by the UNM which often has filamentary
morphology. As shown in Figure 5, from the bighicat catalog we see that
regions with $A_{V}>1$ have slightly colder peak $T_{s}$, due to increased
shielding, and also exhibit less UNM (at $T_{s}\sim 200-600$ K) relative to
the more diffuse H i environments with $A_{V}<1$. This is likely caused by
more pronounced dynamical processes at low-$A_{V}$, and more shielding at
high-$A_{V}$, in agreement with suggestions by Wolfire (2015). Future
observational constraints of the UNM spatial distribution will be very
important to help determine the processes that drive H i out of equilibrium.

###  5.5 The WNM temperature and fraction

As the WNM usually dominates H i emission spectra, its properties are
sometimes derived purely from H i emission data. Kalberla & Haud (2015, 2018)
decomposed the HI4PI survey spectra into Gaussian components and selected
three phases (CNM/UNM/WNM) based on the width of Gaussian functions. They
found that the WNM was smoothly distributed with a mass fraction of $\sim
32$%. Similarly, Marchal & Miville-Deschênes (2021) decomposed the high
resolution $12\times 12$ degree GHIGLS field into Gaussian components and
found that 64% of the H i mass is in the WNM. However, the standard deviation
of that value is relatively high (35% of the median value). By comparing H i
with extinction data, they estimated the WNM temperature $T_{k}=(6.0\pm
1.3)\times 10^{3}$ K.

To robustly estimate the excitation temperature of the WNM, without a
corruption by turbulent motions, absorption-emission pairs are required. High
sensitivity interferometric H i absorption observations have been used to
measure the WNM temperature. For example, Carilli et al. (1998) measured
$T_{s}=5500-8700$ K in the direction of Cygnus A. Dwarakanath et al. (2002)
and Kanekar et al. (2003) found lower values of 3000-3600 K.

The very sensitive H i absorption survey, 21-SPONGE, detected a handful of
components with $T_{s}>1000$ K (e.g., $<10$% by number in Murray et al.,
2015). To improve sensitivity to shallow, broad absorption features further,
Murray et al. (2014) employed a spectral stacking analysis on 1/3 of the
survey data and detected a pervasive population of WNM gas with
$T_{s}=7200_{-1200}^{+1800}$ K. This result was confirmed using data from the
entire 21-SPONGE survey in Murray et al. (2018), where residual absorption was
stacked and binned based on residual emission, demonstrating the existence of
a significant absorption feature with a harmonic mean $T_{s}\sim 10^{4}$ K.

The range of predicted kinetic temperatures for the WNM from the most detailed
ISM heating and cooling considerations, is $T_{k}\sim 5000–8800$ K (Wolfire et
al., 2003). If the H i is only collisionally excited, as shown by Kim et al.
(2014), this range implies $T_{s}\sim 1000–4000$ K (their Figure 2). When,
however, the common prescription is used to account for the Ly$\alpha$ flux,
it is expected that $T_{s}\sim 3500–5000$ K. The Murray et al. (2014) and
Carilli et al. (1998) results are significantly higher than predictions from
standard ISM models based on collisional H i excitation. As can be seen in
Figure 5, even all combined H i absorption measurements together are not
enough to detect the WNM peak of the H i gas distribution. This suggests that
the WNM spin temperature is likely higher than standard analytical and
numerical predictions. Another problem could be the implementation of the
$T_{k}$ to $T_{s}$ conversion and the use of a constant Ly$\alpha$ flux. In a
recent study by Seon & Kim (2020), where a multi-phase model was used and
photons originating from Hii regions were tracked to produce the Ly$\alpha$
radiation field, the Ly$\alpha$ flux was strong enough to bring the 21 cm spin
temperature of the WNM close to the kinetic temperature. This encouraging
result demonstrates that a careful treatment of Ly$\alpha$ photons needs to be
incorporated in numerical simulations, instead of commonly used uniform
Ly$\alpha$ flux.

In summary, the WNM remains difficult to detect in absorption and requires
long integration times or stacking analyses. The detection difficulty likely
stems from the possibility that it has a higher excitation temperature than
what was expected. Future surveys of H i absorption will be essential to probe
spatial variations of the WNM temperature via stacking. Constraining methods
that use Gaussian decomposition to separate H i phases using simulated data is
essential as studies of external galaxies rely on this method exclusively.

###  5.6 The “Dark” Neutral Medium

The transition from H i to H2 is a critical step in the evolution of galaxies.
The efficiency with which this transition occurs appears to influence star
formation efficiency globally. While the transition from H i to H2 is
extremely important, it is also not easy to study. H2 is impossible to observe
directly in typical star-forming conditions (Carilli & Walter, 2013). Carbon
monoxide (CO) is traditionally used as a proxy for H2, but because it cannot
self-shield it is easily dissociated. In the dense conditions of star-forming
molecular clouds H2 and dust shield CO so that it is readily detectable, but
in the early stages of the atomic-to-molecular transition CO is not observed.
At the same time, the CNM that goes on to make H2 is studied via absorption
measurements that are largely limited to sight-lines in the direction of
background sources. Excluding the dominant tracers of CO and H i, how do we
find gas in its transition from H i to H2?

The Fermi $\gamma$-ray telescope and the Planck telescope independently
discovered a component of so-called “dark” gas, not detectable by H i or CO
but inferred by its excess $\gamma$-ray emission and dust opacity (Grenier et
al., 2005, Planck Collaboration et al., 2011). Astonishingly, the “dark” gas
mass seems to be between 20 and 50% of the total H i and CO detected gas in
the solar neighborhood. It was thought that this might be the vast reservoir
of transition gas, caught between its purely atomic and self-shielding
molecular states. The discovery prompted a flurry of research into whether the
“dark” gas was primarily optically thick atomic gas or primarily CO-dark
molecular gas (e.g. Fukui et al., 2014, 2015, Lee et al., 2015, Murray et al.,
2018, Sofue, 2018, Nguyen et al., 2019).

To resolve the nature of dark gas, indirect tracers are often used. For
example, Fukui et al. (2014, 2015) assumed that, if dust and neutral gas are
well mixed and that the specific dust opacity is constant, the H i optical
depth and spin temperature can be approximated from the observed Planck dust
properties. Under these assumptions they estimated that the H i optical depth
must exceed one along most local sight-lines but that low resolution data
missed significant unresolved dense gas. But conversely, Lee et al. (2015)
used direct measurements of H i absorption to show that in the Perseus
molecular cloud less than 20% of sight-lines show $\tau>1$. Murray et al.
(2018) revisited this problem using some of the techniques of Fukui et al.
(2015) applied to higher angular resolution GALFA-HI data and found that the H
i is largely diffuse at high latitudes down to 4′ scales and the lower angular
resolution observations are not missing optically thick H i. Nguyen et al.
(2018) showed that the H i opacity effects only become important above $N_{\rm
HI}>5\times 10^{20}~{}{\rm cm^{-2}}$.

Meanwhile, several studies have investigated alternative molecular gas tracers
(HCO+ & OH) instead of CO. Liszt et al. (2018) observed HCO+—known to trace H2
even in diffuse environments—in addition to CO towards 13 background sources
in the direction of the Chamaeleon cloud. Whereas CO emission was detected in
only 1 of these directions, HCO+ absorption was detected in 12/13 directions,
and the H2 column densities inferred from HCO+ explained the discrepancy
between the column densities derived from $\gamma$ rays and dust emission with
those estimated from H i and CO. While they found that optically thick H i
could contribute to the dark gas in several directions, their results showed
that the dark gas could primarily be accounted for by dark molecular gas
rather than dark atomic gas (although Hayashi et al. 2019 presented a
contradictory view that Chamaeleon has more optically think H i rather than
CO-dark H2). Independently, Nguyen et al. (2018) showed that the OH could be
used as a proxy for H2 over a broad range of column densities and that
optically thick H i is not a major contributor. In the end, it seems that a
combination of factors including a small amount of optically thick H i ($\sim
10-30\%$ increase in column density), evolution of dust properties and CO-dark
molecular gas can resolve the dust opacity excess.

##  6 CHARACTERIZING THE STRUCTURE OF Hi ON ALL SCALES

The morphology of the H i in the Milky Way is primarily studied through H i in
emission and therefore our knowledge of structure has been dominated by the
WNM. Observational advances with spatial and spectral resolution over the past
twenty years have started to reveal some hints about the CNM structure.
Extending our understanding of the statistical thermal properties of the two
stable phases to an understanding of the differences in morphology of the
Galaxy’s structure in the H i phases is still underway.

The structure of H i in emission can be characterized in two ways:
deterministic and scale-free. Through the ever-improving observational data,
supported by numerical simulations (e.g. Kim & Ostriker, 2018, Hill et al.,
2018), a picture has developed of the H i gas cycle in which the diffuse,
scale-free H i is repeatedly ordered by large-scale compression into dense,
often filamentary and magnetized (McClure-Griffiths et al., 2006),
deterministic structures that catalyze cooling (Dawson et al., 2011) and may
drive turbulence. Deterministic structures are those that likely originate
from a specific event or series of events; of those the most common are the
shells, bubbles and chimneys that dominate the H i sky on angular scales
between one and hundreds of degrees (Heiles, 1979, Hu, 1981, Heiles, 1984,
McClure-Griffiths et al., 2002, Ehlerová & Palouš, 2005).

###  6.1 Deterministic Structure: From Bubbles & Shells to Clouds, Filaments,
& Chimneys

While some of the H i shells observed in the Milky Way are clearly associated
with relics of massive stars (e.g Normandeau et al., 1996, McClure-Griffiths
et al., 2001), many of the largest H i shells have not been clearly connected
to their parent stars. In fact, H i shells are rarely observed around
currently active star-forming regions or even young supernova remnants. It
seems that the H i shells are more often associated with larger, older
objects, although there are some notable exceptions, e.g. GLIMPSE bubble N107
(Sidorin et al., 2014).

Identifying H i shells in the Milky Way is not easy. Until recently the
majority of shells had been cataloged manually (Heiles, 1979, 1984, McClure-
Griffiths et al., 2002). Some attempts had been made to catalog shells semi-
autonomously (e.g. Ehlerová & Palouš, 2005) with varying success. Issues
affecting the success of automated searches relate to the velocity ambiguity
of the inner Galaxy, the prevalence of arcs or shell-like features even in
turbulent data (a quick look at terrestrial clouds can confirm!) and complex
overlapping structures. Most recently, Ehlerová & Palouš (2013) refined
automated shell searches with requirements for completeness of the shell,
velocity widths, etc. The resultant catalog from the LAB survey has 333
entries, distributed across the sky. Considering only shells outside the solar
circle, the fraction of the disk’s surface area covered by shells drops from
$\sim 80$% at the solar circle to $\sim 10$% at 15 kpc. Apart from a handful
of newly identified shells, work on H i shells has largely stagnated over the
last decade. Progress has been partially limited by the lack of distances to H
i structure and the difficulty in defining three-dimensional structures. This
is set to change in the coming years with three-dimensional dust maps helping
to overcome problems with velocity ambiguities. Despite the difficulties in
robustly identifying H i shells, they do have important roles in the evolution
of H i structure and density.

From the various observational studies published on H i supershells throughout
the late 1990’s and 2000’s it became clear that the role of supershells is to
order the diffuse H i into large-scale, filamentary structures visible
throughout the Galactic plane. Does the H i also cool within these structures?
Simulations have shown that supershells sweep up the warm, ambient medium into
dense, structured filamentary walls where the gas can cool to form CNM and
molecular clumps (Koyama & ichiro Inutsuka, 2000, Mac Low et al., 2005,
Hopkins et al., 2012). From an observational perspective this role is less
well documented although there certainly have been examples of supershells
with associated molecular gas (e.g. Fukui et al., 1999, Dawson et al., 2008).
McClure-Griffiths et al. (2003) found clumps of narrow line-width H i in the
walls of a large-supershell, suggesting that the structures were formed by
Rayleigh-Taylor instabilities. Knee & Brunt (2001) and Moss et al. (2013) each
used HISA to show that the walls of two supershell structures were dominated
by cold gas. The most comprehensive analysis of supershells as collectors and
coolers of gas was produced by Dawson et al. (2011), who showed that within
two supershells the walls were dominated by cold gas with estimated
temperatures and densities of $T_{k}\sim 100~{}{\rm K}$ and $n_{H}\sim
10~{}{\rm cm^{-3}}$. Comparing the walls to the surrounding medium, they found
statistically significant indicators of narrow linewidths consistent with
enhanced CNM fraction as well as molecular fraction.

Filamentary structures. Feedback events, like those observed in supershells,
indicate that star-formation helps establish the large-scale structure of the
H i that may ultimately be devolved to smaller scales. Much of that large-
scale structure after the star formation tracers are gone is filamentary. For
example, Soler et al. (2020) used the THOR survey at 40′′ angular resolution
to show that within the Galactic plane the H i has a filamentary nature with
the filaments predominantly aligned with the Galactic plane. In some places
where star formation appears to have been active, the filamentary structure
abruptly switches direction to out of the plane. Extending this work across
the Galactic disk, Soler et al. (2022) showed that in the outer Galaxy the
filamentary H i structures are predominantly parallel to the Galactic plane
and inside the Solar circle they are mostly perpendicular to the plane or
without a preferred orientation. Their interpretation is that the filaments
show the imprint of supernova feedback in the inner Galaxy and Galactic
rotation and shear in the outer Milky Way.

Observationally, distinguishing the characteristic morphology of the CNM from
that of the WNM has been difficult to disentangle because the large linewidth
and diffuse character of the WNM blurs structure within position-position-
velocity cubes of H i emission. The morphology of the CNM, with its generally
smaller brightness temperatures and narrow linewidths, is easily obscured
within spectral channels dominated by warm gas. But spectacular examples of
discrete, cold features that appear cloud-like or filamentary abound. For
example, the very local Leo cold cloud’s exceptionally narrow velocity width
indicates a kinetic temperature of between 11 and 20 K (Peek et al. 2011;
supported by absorption measurements (Murray et al. 2018)). Recent work making
use of the combination of spectral and angular resolution, and sensitivity
that have been afforded by the GALFA-HI survey has revealed small-scale,
narrow linewidth, fibrous structure (Clark et al., 2019). The inference that
the fibrous structures are cold has been confirmed through various techniques
including correlation with NaI absorption (Peek & Clark, 2019), convolutional
neural network CNM separation (Murray et al., 2020) and dust-to-gas comparison
(Kalberla et al., 2020). All three papers hint that the morphology of the CNM,
as observed in emission, is distinct from the WNM. Dense grids of H i
absorption spectra, as well as new techniques that perform spectral
decomposition of H i emission for separating the WNM and CNM by using H i
absorption data as learning templates (e.g. Murray et al., 2020), will be
important for further characterizing the CNM morphology.

Disk-halo clouds. The familiar paradigm that supershells can drive hot gas and
metals is clearly illustrated by many examples of H i supershells that
completely disrupt the H i disk (Normandeau et al., 1996, Heiles, 1998,
McClure-Griffiths et al., 2002, Pidopryhora et al., 2007, Dawson et al.,
2008). They also seem to contribute to defining of the structure the disk-halo
interface region into “clouds” or “cloudlets”. In their 1990 review of H i in
the Milky Way Dickey & Lockman (1990) asserted that the concept of H i clouds
was not a useful one. Indeed, a quick examination of H i emission images, such
as those shown in the movie in the Supplementary material will confirm that
assertion - how could any of this structure be defined as a “cloud”?
Nevertheless, this description of H i has persisted in the ISM lexicon. In
spite of its general appearance there are a few examples of warm H i in the
disk-halo interface region that can reasonably be described as “cloud”-like.
Lockman (2002) identified an unexpected population of H i structures that are
clearly discrete and bounded and therefore cloud-like. These parsec-scale
clouds are near the tangent point of the inner Galaxy, where small offsets
from circular velocity place them at nearly forbidden velocities. Subsequent
studies suggest that the disk-halo cloud population is not completely
ubiquitous throughout the Milky Way (Stil et al., 2006, Stanimirović et al.,
2006, Ford et al., 2008, 2010, Pidopryhora et al., 2015) but instead they may
be a result of ISM structure driven in some indirect way from star formation.
Ford et al. (2010) connected the clouds with broken remnants of shells and
supershells driven away from the disk.

Figure 8: Top: H i image at $v=-7~{}{\rm km~{}s^{-1}}$ from GALFA-HI (Peek et
al., 2018) showing fibrous structure overlaid with starlight polarization
vectors (Heiles, 2000). Bottom: Rolling Hough Transform (Clark et al., 2014)
backlit image of filamentary H i structure from the same H i data, overlaid
with the same polarization vectors.

###  6.2 Linking magnetic fields to H i structure

Although the link between atomic gas and magnetism is not necessarily
intuitive, early work on the correlation of starlight polarization with many
large-scale local H i structures, such as Loop I and the Orion-Eridanus bubble
(Heiles, 1997), indicated that magnetic fields are associated with the
morphological structure of H i. Using observations of Zeeman splitting of the
H i line in absorption, Heiles (1989) showed that for a sample of lines-of-
sight through local supershells the magnetic field strengths were larger than
the typical ISM values. The implication from the Zeeman work and the
correlation of starlight polarization angle with H i loops is that the
magnetic field and H i were co-structured through a common formation
mechanism, such as a supershells (Heiles, 1989, 1998).

H i Zeeman work culminated with the Millennium survey of H i absorption,
giving 69 measurements on arbitrary lines-of-sight with magnetic field
strengths of $|B_{los}|\leq 10~{}{\rm\mu G}$ (Heiles & Troland, 2004, 2005).
Through comparison of $|B_{los}|$ with column density they showed that the
field strengths were not consistent with flux freezing, as is typically found
in molecular clouds (Crutcher, 1999). The technical challenges of H i Zeeman
observations have limited its extensive application and large-scale
conclusions about the connection of magnetic field strengths with the
structure of the neutralmedium.

Over the last decade there has also been a growing body of work indicating
that H i structure is linked to magnetic field structure, as traced by
starlight and dust polarization (e.g. McClure-Griffiths et al., 2006, Clark et
al., 2014, Clark & Hensley, 2019). Figure 8 shows an example H i velocity
channel image with the fibrous structure emphasized by the Rolling-Hough
Transform overlaid withstarlight polarization. So far, most of the
observational evidence linking H i structure to magnetic fields has focused on
the CNM. Clearly the H i Zeeman measurements, by nature of tracing absorption,
exclusively measure CNM. This connection with the CNM was borne out in Clark
et al. (2019) and Kalberla et al. (2020), who showed that the magnetically
correlated fibrous structure of the H i is dominated by narrow linewidths.
Similarly, the structure of the cold Riegel-Crutcher HISA cloud (discussed in
Section 4.3) was shown to have outstanding agreement with starlight
polarization (McClure-Griffiths et al., 2006). Extending to the whole sky,
Kalberla et al. (2021) found remarkably good agreement between far-infrared
polarization angle and H i filament orientation angle for velocity widths of
$1~{}{\rm km~{}s^{-1}}$, but that for larger velocity widths the agreement
disappeared. The conclusion seems to be that H i filamentary structure is cold
and magnetically dominated with the field orientation parallel to the
filaments.

Turning the technique around, the fibrous H i structure has been used as a
proxy for direct observations of magnetic fields, adding a velocity dimension
to the observational technique of dust polarization. Clark & Hensley (2019)
performed a full-sky decomposition of the fibrous H i that led to a velocity
resolved projection of the plane-of-sky magnetic field. An alternative method
proposed by González-Casanova & Lazarian (2017) and Yuen & Lazarian (2017)
utilizes the connection of turbulent motion to the magnetic field by using
gradients of the velocity amplitude, which are expected to be perpendicular to
the magnetic field direction, to trace the magnetic field orientation. In a
series of papers, Hu et al. (2021) employed another technique focusing on the
structure function of simulated H i to reveal both the magnetic field
direction and its strength. The underlying correlation of H i structure and
magnetic fields is an area of great activity and will likely develop rapidly
over the coming years.

The fact that the two ISM elements - magnetic fields and CNM - correlate so
well is suggestive that the magnetic field is either dynamically important in
setting the structure of H i in many environments or that the same dynamical
events (e.g. supershells) that help produce CNM also enhance magnetic field
strengths. While most observations linking magnetic fields and H i have
focused on starlight or dust polarization as the tracer of magnetic field,
there have been several studies in the past few years that have used Faraday
rotation of radio continuum as the magnetic field tracer (Van Eck et al.,
2019, Bracco et al., 2020, Turić et al., 2021). Radio continuum studies may
help to break the bias of magnetic fields traced in cold gas by dust, instead
showing magnetic fields in the warm, ionized medium probed by Faraday
rotation. However, recent low radio frequency polarimetric observations of
Faraday rotation have generally found morphological association with cold H i
structures. For the first time, Campbell et al. (2022) used mid-frequency
radio polarization Faraday rotation to show that the H i appears clearly
connected to the warm ionized medium via a common magnetic field. While they
consider whether a compressive event triggered the connection of the H i and
magnetic field, they cannot definitively conclude that it did. Future studies
of Faraday rotation and H i will be important to reveal whether magnetism is
as crucial to the structure of the WNM as the CNM.

###  6.3 Statistical description of the H i distribution

Statistical studies provide an important and complementary way to describe and
characterize complex, scale-free structure seen in high-resolution H i images.
As interstellar turbulence is known to be ubiquitous throughout the ISM and
greatly responsible for the H i distribution in galaxies, e.g. Elmegreen &
Scalo (2004), statistical descriptors of H i are important for constraining
the drivers and the nature of interstellar turbulence. In addition, the H i
turbulence constraints are essential for understanding the process of
molecular cloud formation and evolution as studies of molecular clouds (e.g.,
Hennebelle & Falgarone, 2012) show that most of the turbulent energy is
injected from outside.

Turbulence can be driven by many energy sources that operate on very different
spatial scales. The key candidate is stellar feedback – a term used to express
the combined effect of proto-stellar jets, stellar winds, photo-heating, and
supernovae (Offner & Arce, 2014, Grisdale et al., 2017). Turbulent drivers
operating mainly on larger spacial scales include gravitational instabilities
(Elmegreen et al., 2003, Krumholz & Burkhart, 2016), thermal instabilities
(Kritsuk et al., 2007, Kim & Kim, 2013, Iwasaki & Inutsuka, 2014), and
magneto-rotational instabilities (Piontek & Ostriker, 2004, 2005, 2007).
However, clear observational signatures of different drivers of turbulence are
still lacking, partially due to the need for extremely high spatial dynamic
range (from few pc all the way to many kpc).

Many different statistical approaches have been used on H i observations to
characterize complex structure. Some of these approaches were summarized in
Elmegreen & Scalo (2004) and Scalo & Elmegreen (2004). Here we focus on the
spatial power spectrum of H i as one of the common tools used to analyze H i
emission, with a goal of compiling published results and pulling together some
common trends.

####  6.3.1 H i spatial power spectrum (SPS)

Figure 9: The SPS slope obtained for different H i regions. All data points
were derived from the H i column density images integrated over a certain
range of velocities, mostly focusing on the low-velocity HI. Black data points
are from Mittal et al. (in prep). Green: Pingel et al. (2013), Pingel et al.
(2018), Choudhuri & Roy (2019), Marchal et al. (2021), Miville-Deschênes et
al. (2003), Blagrave et al. (2017), Martin et al. (2015), Khalil et al.
(2006), Dickey et al. (2001), Kalberla & Haud (2019), Marchal & Miville-
Deschênes (2021), Green (1993), Crovisier & Dickey (1983), Chepurnov et al.
(2010).

In the presence of turbulence the 2D spatial power spectrum (SPS) of H i
intensity or column density fluctuations is expected to have a single or
double power-law functional form (Scalo, 1987). The SPS of a 2D image,
$I(\mathbf{x})$, is defined as

$P(k)=\int\int\langle
I(\mathbf{x})I(\mathbf{x}^{\prime})\rangle\exp^{i\mathbf{L}\cdot
k}d\mathbf{L},$ (11)

where $k$ is the spatial frequency and
$\mathbf{L}=\mathbf{x}-\mathbf{x}^{\prime}$ is the distance between two points
in the image. The power-law slope of the SPS, $P(k)\propto k^{\gamma}$,
derived from 2D H i images reflects contributions from the underlying 3D
velocity and/or density fluctuations (Lazarian & Pogosyan, 2000).
Theoretically, the SPS can illuminate the spatial scales at which energy is
injected or dissipated, as well as the mechanisms by which that energy
cascades from large to small scales or vice versa.

Practically, however, the interpretation of SPS results is complex due to many
observational and theoretical effects. The SPS slope depends on the line-of-
sight complexity (e.g. multiple H i structures in the same direction), the
multi-phase nature of H i, and optical-depth effects. Furthermore, the disk
geometry and the H i scale height can result in the transition between 2D and
3D turbulence, producing breaks in the SPS. Recent studies have shown that the
interpretation of SPS breaks can be greatly affected by telescope instrumental
effects (beam and noise), and the multi-phase nature of the medium (Koch et
al., 2019, Körtgen et al., 2021). Similarly, it has been shown that separating
density and velocity fluctuations from observed intensity fluctuations
requires a careful consideration. The approach for doing this by varying the
width of velocity channels (Lazarian & Pogosyan, 2000) has been recently
updated with a more complex treatment by Yuen et al. (2021). However, some
debate still exists. Clark et al. (2019) and Kalberla et al. (2016) suggested
that H i intensity fluctuations seen at different velocities are phase
dependent, with narrow velocity channels being dominated by the CNM and
causing the H i SPS slope to be more shallow than what is found for the
integrated intensity images. Finally, it is currently not clear whether the
SPS of H i emission is sensitive in distinguishing between different turbulent
modes within a multi-phase medium.

From the observational point of view, although there are many measurements of
the H i SPS, most studies employed different analysis techniques, resulting in
heterogeneous samples that are hard to compare, making it difficult to search
for spatial variations of turbulent properties with interstellar environments.
Early H i SPS studies by Crovisier & Dickey (1983) and Green (1993) fit a
single power-law with indices of $-3.0$ and $-2.8$, respectively. Subsequent
large-area (many degrees across) H i maps have also been used to calculate a
single SPS. As many studies of H i SPS showed no evidence for a turnover on
large spatial scales, the interpretation has been that H i turbulence in the
Milky Way is being driven on scales larger than the size of individual H i
structures, (e.g. Pingel et al., 2013, 2018). This is somewhat puzzling
considering the existence of many star-forming regions and supernovae in the
disk, which we expect would modify turbulent properties of H i.

In Figure 9 we attempt to compare values of the H i SPS slope obtained from
the literature. We note that all these measurements were obtained from H i
images integrated over a range of velocity channels that roughly corresponds
to the local H i ($\pm 20-40$ km s-1) and likely represent estimates of
density fluctuations (based on Lazarian & Pogosyan, 2000, Yuen et al., 2021).
We also show the application of the rolling SPS analysis on the entire region
observed by the GALFA-HI survey from Mittal et al. (in prep). The rolling SPS
method, introduced by Szotkowski et al. (2019), calculates SPS locally by
rolling a spatial kernel across an image. This allows for a consistent
handling of instrumental effects.

The GALFA-HI survey covers about 13,000 square degrees and includes a region
in the inner Galaxy, a region in the outer Galaxy, as well as some high-
latitude H i. The SPS slope, calculated over the angular range of scales from
1.5 degrees to 16 arcmin, varies between $-2.6$ and $-3.5$. Close to the
Galactic plane, $|b|<10$ degrees, the SPS slope is the steepest, while at
higher Galactic latitudes the SPS slope gradually becomes more shallow
reaching to $\sim-2.7$ at $|b|>70$ degrees (Figure 9). While the gradual
change of SPS slope with Galactic latitude cannot be explained with the
optical-depth effects, it becomes less significant when the H i column density
is corrected for the difference in the line-of-sight length caused by the
plane parallel geometry of the HI disk.

It is interesting to see that the H i SPS slope appears relatively uniform
across large spatial areas, and is only consistently slightly more shallow at
high Galactic latitudes. The GALFA-HI survey samples the range of H i column
densities from $\sim 10^{20}$ to $\sim 10^{22}$ cm-2, yet the SPS slope
remains largely around $-2.8$ for $|b|<50$ degrees. This could be due to the
dominance of large-scale turbulent driving, on Galaxy- or kpc-wide spatial
scales, possibly caused by a combination of Galactic rotation, gravitational
instabilities, infall onto the Galaxy, and even energy injection from large
supershells (Wada et al., 2002, Bournaud et al., 2010, Krumholz & Burkhart,
2016). Numerical simulations by Yoo & Cho (2014) showed that the turbulent
energy spectrum is very sensitive to large-scale driving. Unless the energy
injection rates on small scales are much higher than the energy injection rate
on large scales, the large-scale driving will always dominate. Therefore, even
in the presence of significant stellar feedback in the Galactic plane, large-
scale turbulent driving might be the dominant mechanism. The finding of a
relatively uniform SPS slope for the Milky Way resembles the situation with
the H i in the Small Magellanic Cloud where Szotkowski et al. (2019) found a
very uniform SPS slope across the entire galaxy.

Hennebelle & Falgarone (2012) showed that the kinetic energy transfer rate
measured for the population of molecular clouds traced by 12CO showed no
variation with spatial scale ranging from $\sim 0.01$ to $\sim 100$ pc. They
concluded that this result suggests that molecular clouds are part of the same
turbulent cascade as H i. From this point of view, the relatively uniform H i
SPS slope may not be so surprising – if molecular clouds of varying levels of
star formation and stellar feedback sources do not show significant changes of
their internal turbulent properties, the H i may be also not be strongly
affected by small-scale turbulent driving.

In summary, while the interpretation of the H i SPS remains complicated, some
trends are starting to emerge by analyzing large surveys with consistent
methodologies. In moving forward, it is essential to develop a consistent and
uniform approach for handling data, as well as more sophisticated methods for
disentangling geometry, phase, and turbulent mode effects.

###  6.4 3-D Distribution

####  6.4.1 Radial Distribution

By nature of the large-scales involved, our overall understanding of the 3-D
distribution of the Milky Way in H i has not changed much since the excellent
review by Kalberla & Kerp (2009). They carefully explored the radial
distribution, the warp and the scale height of the warm disk. Key
characteristics included an asymmetric warp beyond $R\approx 9$ kpc and a
rapidly flaring disk beyond about $R>18$ kpc.

Like most H i rich galaxies, the Milky Way has an extended H i disk. What is
the true extent of the H i disk and what form does it take? It is difficult to
associate a distance for the outer Milky Way because observational methods
are, in general, limited to kinematic distances. With the assumption of a
mostly flat rotation curve Kalberla & Dedes (2008) found that the H i disk
extends beyond $R\sim 30$ kpc, but the H i surface density there is low:
$N_{H}\sim 1.5\times 10^{19}~{}{\rm cm^{-2}}$ or $0.1~{}{\rm
M_{\odot}~{}pc^{-2}}$. At even larger galactocentric radius, the disk likely
extends as far as 60 kpc but with a very high velocity dispersion
($\sigma_{v}=74~{}{\rm km~{}s^{-1}}$), which cannot be thermal and therefore
must suggest a clumpy distribution with a high cloud-to-cloud velocity
dispersion. A large radial extent to the H i disk has also been observed in
other galaxies, for example NGC 2997 (Pisano, 2014) shows clumps of H i up to
60 kpc.

The tail end of the H i radial distribution, however, is at odds with previous
theoretical estimates of the column density cutoff due to ionization from the
CMB (e.g. Dove & Shull, 1994, Maloney, 1993), which suggest that the H i
column density should drop precipitously after $R\sim 35~{}{\rm kpc}$. Bland-
Hawthorn et al. (2017) revisited the extended H i disk models including the
effects of non-uniform gas distribution or reducing the scale height. They
showed that if the H i distribution was compressed into a narrower scale-
height, increasing the mean volume density by a factor of 100 to $n_{H}\sim
10~{}{\rm cm^{-3}}$, the disk could survive as H i with $N_{H}>10^{18}~{}{\rm
cm^{-2}}$ to $R>55$ kpc. But the radial extent of the disk is strongly
dependent on the assumed mean volume density; reducing to $n_{H}\sim 1~{}{\rm
cm^{-3}}$ decreases the disk radius at the same sensitivity limit to $R\sim
45$ kpc.

Explaining the Milky Way’s apparently large, low column, radial extent relies
on both increasing the mean volume density and also explaining the high
velocity dispersion. The latter can be possibly attributed to the discrete
extended H i line wings like those observed by Kang & Koo (2007) and
tentatively associated with old supernova remnants. The evidence for the large
volume density predicted by Bland-Hawthorn et al. (2017) is not as clear.
Certainly, we know from Strasser et al. (2007) and Dickey et al. (2009) that
dense material exists at radii $15<R<25$ kpc with estimated number densities
of $n_{H}\approx 0.7-2.4~{}{\rm cm^{-3}}$.

It is worth noting that there is some direct evidence for the H i extending at
least as far as $R=23$ kpc. This comes from the interaction between the
Magellanic Leading Arm and the Milky Way disk (McClure-Griffiths et al., 2004,
Price-Whelan et al., 2019) at $R=23$ kpc, which has resulted in an orphan star
cluster associated with the Magellanic Leading Arm and formed from the
interaction with the disk. Nidever et al. (2019) used the full three-
dimensional velocity vector and star cluster’s orbit to predict the gas
density of the disk through which the cluster passed. This interesting example
provides an unusual method to estimate the density of the H i disk at the
midplane at $R=20$ kpc: $n_{H}=6.0\times 10^{-2}~{}{\rm cm^{-3}}$.

The CNM and WNM have different spatial extents. Unlike the warm gas, the scale
height and extent of the CNM disk has not been well characterized. For the
most part this has been an observational limitation; H i emission is dominated
by warm gas and the only places where the profiles are simple enough to allow
CNM-WNM decomposition is at high latitudes and dominated by the local gas.
Nonetheless, the sheer number of H i absorption measurements enabled by the
Galactic plane surveys (CGPS, VGPS, CGPS, THOR; Taylor et al., 2003, Stil et
al., 2006, McClure-Griffiths et al., 2005, Wang et al., 2020c), has given some
indications of the extent of the cold disk.

The cold H i disk extends to at least 25 kpc and shows a distinct pattern in
longitude-velocity space that nominally traces the spiral arm distribution
(Strasser et al., 2007). As we already discussed in Section 5.2, Dickey et al.
(2009) showed that the CNM fraction remained mostly flat with Galactic radius
out to at least 20 kpc, with the assumption that the average temperature of
the cold component does not change. Lately, Dickey et al. (2022) has revisited
that result by stacking 175 H i absorption spectra towards extragalactic
sources from the new Galactic ASKAP (GASKAP) Pilot Survey and found that the
2009 result stands with a constant CNM fraction out to $R\sim 40$ kpc. These
results are somewhat surprising; the outer Galaxy environment has a much
metallicty (possibly $<1/10$ the Solar circle gas phase O/H abundance (Wenger
et al., 2019)) and the balance of heating and cooling is likely quite
different.

####  6.4.2 Scale Height

The scale height of the H i emission layer and its variations with
Galactocentric radius exterior to the solar circle have been well measured by
Kalberla & Dedes (2008) and summarized in the Kalberla & Kerp (2009) review.
Kalberla & Kerp (2009) described the radial dependence of the full-width at
half-max, FWHM, of the emission layer between $R=5~{}{\rm kpc}$ and
$R=35~{}{\rm kpc}$ as $f=f_{0}e^{(R-R_{0})/R_{0}}$, where $f_{0}=300~{}{\rm
pc}$ and $R_{0}=9.8$ kpc. Interior to the solar circle, from $4\lesssim
R\lesssim 8$ kpc, the density, $n(z)$ of the H i disk in the inner Galaxy has
long been characterized by a dual Gaussian with FWHM of 212 pc and 530 pc,
with an additional low density exponential tail (Lockman et al., 1986, Dickey
& Lockman, 1990, Savage & Wakker, 2009). The full distribution has a FWHM of
430 pc but Lockman (1984) also noticed that inside about $R=2$ kpc, the scale
height appeared to decrease. This was confirmed in detail by Lockman &
McClure-Griffiths (2016), who showed that from $R\approx 2.1$ kpc on the
negative longitude side of the Galactic center to $R\approx 2.4$ kpc on the
positive longitude side the large scale height layer is missing and the disk
is fit very well by a single Gaussian with FWHM $\sim 125$ pc. This diminished
disk is matched with a dramatic reduction in the H i surface density in the
inner 2.5 kpc of the Galaxy. At $R=3.5-3.75$ kpc the surface density is
$\Sigma_{HI}=3.4~{}{\rm M_{\odot}~{}pc^{-2}}$, dropping to
$\Sigma_{HI}=0.4~{}{\rm M_{\odot}~{}pc^{-2}}$ in the annulus from $1.5$ kpc to
$1.75$ kpc. The cause of the significant drop in H i surface density is not
clear, with possible causes related to the conversion of H i to molecular
material, evacuation due to Galactic centre activity or simply an
observational effect of non-circular orbits (Lockman & McClure-Griffiths,
2016).

The scale height of the H i emission disk may also vary in the outer Galaxy.
Working from a flat rotation curve Levine et al. (2006) transformed the
brightness temperature distribution of H i emission, $T(l,b,v)$ to an
estimated number density distribution $n(R,z,\phi)$ and found variations in
the scale height of 15-20% correlating with variations in the H i gas surface
density $\Sigma_{HI}$ of 30-50%. Using a similar technique, but fitting for
the rotation curve as well, Kalberla et al. (2007) showed the same general
trends. The conclusion from these works is that the gas disk is thinner where
the surface density is higher, and that these variations correlate with spiral
structure. It should be noted, however, that these conclusions are based on
nearly flat, azimuthally smooth rotation curves. We know that dense spiral
arms can drive up to 10% deviations from a flat rotation curve (Burton, 1971,
McClure-Griffiths & Dickey, 2016, Peek et al., 2022), which can drive
systematic variations in distance dependent properties. The conclusions of
varying H i scale height correlating spiral arms should therefore be viewed
with some caution.

While the distribution of the H i emission layer has been well characterized
across the Galaxy, the distribution of the cold H i disk has been much more
difficult to constrain. Focusing on the solar neighborhood, Crovisier (1978)
pioneered a technique to ascertain the scale height of the cold H i disk at
the solar circle from observations of H i in absorption. By comparing the
radial velocities in a catalog of absorbing H i clouds at $|b|>10^{\circ}$
(isolating only the local clouds) to those expected for clouds in a
differentially rotating disk, Crovisier (1978), and later Belfort & Crovisier
(1984), estimated the mean distance of local H i clouds from the plane,
$\langle|z|\rangle\approx 100~{}\rm{pc}$. Both Crovisier (1978) and Belfort &
Crovisier (1984) found that $\langle|z|\rangle$ increased at higher latitudes
($|b|\gtrsim 20^{\circ}$), interpreted as a result of the low abundance of H i
clouds very near to the Sun. For a Gaussian distribution
$\langle|z|\rangle=100~{}{\rm pc}$ corresponds to a FWHM of almost 300 pc.

Figure 10: Normalized distributions of average fractional absorption,
$\left<1-e^{-\tau}\right>$, and average emission, $\left<T_{b}\right>$, versus
Galactic latitude between $l=339^{\circ}$ and $l=342\mbox{${}^{\circ}$}$ for
velocities close to the inner Galaxy tangent point at $R\sim 2.8$ kpc (left)
and in the outer Galaxy at $R\sim 12.5$ kpc (right). The upper axis in $z$ is
for an assumed distance of $d=7.3$ kpc and $d=20$ kpc for the tangent point
and assuming a flat rotation curve with $R_{\odot}=8.31$ kpc and
$\Theta_{0}=240~{}{\rm km~{}s^{-1}}$. Absorption data were obtained from
Dickey et al. (2022) and emission data from GASS (McClure-Griffiths et al.,
2009, Kalberla et al., 2010).

The irregular aerial sampling of H i absorption through the low Galactic
latitudes, which are most useful for probing the scale height of the cold
disk, has been a limitation for understanding the CNM disk structure away from
the solar neighborhood. Dickey et al. (2009) used the Galactic Plane surveys
to produce a projection of H i opacity as a function $R$ and $\phi$ outside
the Solar circle. Recently, the Dickey et al. (2022) measurements of H i
absorption towards 175 extragalactic continuum sources in the GASKAP Pilot
region have provided the first indications of the scale height of the CNM disk
in several regions away from the solar circle. Dickey et al. (2022) stacked H
i absorption spectra from the region $338^{\circ}\leq l\leq 342^{\circ}$,
focusing on regions exterior to the solar circle and in the inner galaxy at
the sub-central point where the line-of-sight is tangent to a circular orbit
around the Galactic center and the H i reaches its terminal velocity (e.g.
McClure-Griffiths & Dickey, 2016). These measurements, together with Solar
neighborhood measurements (Crovisier, 1981, Belfort & Crovisier, 1984), are
starting to give a clear picture of variations in the height cold H i disk
with Galactocentric radius.

The $z$ structure of the emission and absorption disks in the inner and outer
Galaxy are compared in Figure 10. In this Figure we have compiled the Dickey
et al. (2022) data, showing the comparison of the averaged H i absorption
$\left<1-\exp(-\tau)\right>\Delta v$ to the averaged H i emission
$\left<T_{b}\right>$ as a function of latitude in a region in the inner
($R\sim 2.7~{}{\rm kpc}$) and outer ($R\sim 12.8~{}{\rm kpc}$) Galaxy. The
inner Galaxy region shows only the absorption and emission at the terminal
velocity $-145~{}{\rm km~{}s^{-1}}<v_{LSR}<-125~{}{\rm km~{}s^{-1}}$, where
the distance is geometrically defined at $d\sim 7.8~{}{\rm kpc}$ for
$R_{\odot}=8.3~{}{\rm kpc}$. The outer Galaxy region shows absorption and
emission between $15~{}{\rm km~{}s^{-1}}<v_{LSR}<35~{}{\rm km~{}s^{-1}}$,
where for a flat rotation curve $d\approx 20$ kpc and $R\sim 12.5$ kpc.
Fitting to the normalized $z$ distributions, we find that the FWHM of the
absorption disk is 185 pc at $R\sim 2.8~{}{\rm kpc}$ and 490 pc at $R\sim
12.5~{}{\rm kpc}$, while the inner and outer Galaxy emission FWHM’s are 245 pc
and 710 pc, respectively. The fits are consistent with Dickey et al. (2022),
yet we also confirm the classic view that for the inner Galaxy the emission is
better fit with two Gaussians of FWHM 110 and 340 pc, while the low resolution
of HI4PI limits the outer Galaxy fit to a single Gaussian component. The key
conclusions from the comparison are: the absorption disk is narrower than the
emission disk at all Galactic radii and the CNM and emission disks are
narrower in the inner Galaxy than in the outer Galaxy. When combined with the
Crovisier (1981) a picture develops of the cold H i disk increasing in
thickness from a FWHM of $\sim 185$ pc at $R\sim 2.8~{}{\rm kpc}$, through 300
pc around the Sun, to $\sim 490$ pc at $R\sim 12.5$ kpc (see Figure 11). These
conclusions are not necessarily surprising given the higher pressure for the
inner Galaxy, but suggest that the cold H i gas disk more closely matches the
$H_{2}$ disk than the warm H i.

Figure 11: Scale height comparison of cold and warm H i with CO disk. Blue
points are the FWHM of the absorption disk, orange points are FWHM of the
emission disk and purple points are the FWHM of the CO from Heyer & Dame
(2015) and the dashed line is the Kalberla & Kerp (2009) fit for $n(z)$ based
on the assumption of hydrostatic equilibrium. Absorption data are from
Crovisier (1981) and Dickey et al. (2022), emission FWHM were calculated from
Gaussian fits to HI4PI data (HI4PI Collaboration et al., 2016).

Ideally, analyses like Dickey et al (2009, 2022) will be extended across the
Galactic disk to help us realize an azimuthally resolved view of the disk
scale height. If non-circular motions could be accounted for with future
stellar derived mass models for the Milky Way we might also hope to be able to
incorporate those and overcome problems related to the rotation curve.

##  7 FUTURE SURVEYS OF ATOMIC HYDROGEN IN THE MILKY WAY