# The sub-critical illusion: synthetic Zeeman effect observations from
galactic zoom-in simulations

Zipeng Hu,${}^{\href https://orcid.org/0000-0002-3758-552X\,1}$ Benjamin D.
Wibking,${}^{\href https://orcid.org/0000-0003-3175-2291\,1,2}$ and Mark R.
Krumholz${}^{\href https://orcid.org/0000-0003-3893-854X\,1,2}$
1Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia
2ARC Centre of Excellence for Astronomy in Three Dimensions (ASTRO-3D),
Canberra, ACT 2611, Australia
<EMAIL_ADDRESS>(ZPH)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Mass-to-flux ratios measured via the Zeeman effect suggest the existence of a
transition from a magnetically sub-critical state in H i clouds to a super-
critical state in molecular clouds. However, due to projection, chemical, and
excitation effects, Zeeman measurements are subject to a number of biases, and
may not reflect the true relations between gravitational and magnetic
energies. In this paper, we carry out simulations of the formation of
magnetised molecular clouds, zooming in from an entire galaxy to sub-pc
scales, which we post-process to produce synthetic H i and OH Zeeman
measurements. The mass-to-flux ratios we recover from the simulated
observations show a transition in magnetic criticality that closely matches
observations, but we find that the gravitational-magnetic energy ratios on
corresponding scales are mostly super-critical, even in the H i regime. We
conclude that H i clouds in the process of assembling to form molecular clouds
are already super-critical even before H2 forms, and that the apparent
transition from sub- to super-criticality between H i and H2 is primarily an
illusion created by chemical and excitation biases affecting the Zeeman
measurements.

###### keywords:

ISM: clouds – ISM: magnetic fields

††pubyear: 2023††pagerange: The sub-critical illusion: synthetic Zeeman effect
observations from galactic zoom-in simulations–The sub-critical illusion:
synthetic Zeeman effect observations from galactic zoom-in simulations

##  1 Introduction

The dynamical importance of magnetic fields (B-field) relative to gravity
during the star formation process remains an open question. Quantitatively,
for a uniform spherically symmetric cloud of mass $M$, volume $V$, and radius
$R$, threaded by a uniform magnetic field of intensity $B$, the ratio of
magnetic to gravitational energy is

$\frac{E_{\rm mag}}{E_{\rm
grav}}=\frac{B^{2}V}{8\pi}\times\frac{3R}{5GM^{2}}\propto\frac{B^{2}R^{4}}{M^{2}}\propto\left(\frac{\Phi}{M}\right)^{2},$
(1)

where $\Phi=\pi R^{2}B$ is the magnetic flux, which is conserved under ideal
magnetohydrodynamics (MHD). Mouschovias & Spitzer (1976) show that clouds can
collapse only if their mass-to-flux ratio $M/\Phi$ exceeds a critical value

$\left(\frac{M}{\Phi}\right)_{\rm crit}=\frac{1}{3\pi}\sqrt{\frac{5}{G}}.$ (2)

A cloud with $\mu_{B}\equiv(M/\Phi)/(M/\Phi)_{\mathrm{crit}}>1$ is called
super-critical, implying a weak magnetic field compared to gravity, while a
cloud with $\mu_{B}<1$ has strong magnetic field support against gravity and
is called sub-critical. In a sub-critical cloud, the magnetic field strength
will rapidly grow during a uniform gravitational collapse, preventing further
contraction long before stars form (Mouschovias & Spitzer, 1976). Thus if the
bulk of the interstellar medium (ISM) is magnetically sub-critical, non-ideal
MHD effects are required to allow stars to form. One candidate effect is the
ambipolar diffusion (Shu, 1977; Shu et al., 1987), whereby in the very poorly
ionised regions found at high densities, neutral gas becomes decoupled from
ions, allowing it to cross the field lines and further collapse. Alternately,
turbulent motion of the gas might lead to magnetic reconnection and change the
$B$-field topology, releasing magnetic flux and enabling collapse (Vishniac &
Lazarian, 1999). Understanding whether such non-ideal processes are needed,
and if so which ones are dominant, requires observational work to precisely
measure the magnetic field structure both in and around molecular clouds.

The most direct method of measuring magnetic fields, and by far the least
uncertain when it comes to measuring field strengths, is via the Zeeman
effect, which reveals the strength of magnetic fields along the line of sight
($B_{\rm LOS}$) via the frequency splitting these fields induce between the
components of emission or absorption lines with different polarisations. By
using a range of Zeeman-sensitive tracers (most commonly H i, OH, and CN), the
method can cover a large range of total hydrogen column density ($N_{\rm H}$).
Crutcher (2012) summarizes the results of 5 Zeeman effect surveys (Crutcher et
al., 1999; Bourke et al., 2001; Heiles & Troland, 2004; Falgarone et al.,
2008; Troland & Crutcher, 2008), which show a clear transition at $N_{\rm
H}=10^{21-22}\>\rm cm^{-2}$: diffuse H i regions with lower column densities
are mostly sub-critical, while OH and CN measurements that trace denser
regions indicate that they are mostly super-critical. Troland et al. (2016)
make much more precise Zeeman effect measurements with both H i and OH tracers
in the Orion Veil, and Thompson et al. (2019) conduct an OH absorption line
Zeeman effect survey in the molecular cloud envelopes. Both these studies find
results consistent with the Crutcher (2012) review, showing a sub-critical to
super-critical transition in the same column density range. The mechanism
behind such a transition, however, is still unclear, partially due to the lack
of observations in the atomic-molecular transition regime. In an early attempt
to explore this regime, Ching et al. (2022) perform H i narrow self-absorption
(HINSA) Zeeman effect observation on L1544, and their result suggests that the
cloud is already super-critical when H2 molecules start to form.

Partly in an attempt to probe the sub- to super-critical transition, quite a
few theoretical works have also investigated magnetic fields in the cloud
formation process using MHD simulations. Bertram et al. (2012) determine the
mass-to-flux ratios in cores ($\sim$ 0.1 pc) and clumps ($\sim$ 1 pc) from
their simulations starting with a range of initial magnetic field strengths.
They find that the mass-to-flux ratio of the structures they form is only
loosely related to the initial criticality of the large scale environment. The
recent work by Ibáñez-Mejía et al. (2021) starts from a galactic-scale
simulation with a box size of 1 kpc, then zooms into 3 dense clouds, reaching
a maximum resolution of 0.1 pc. They report the criticality transition happens
when the hydrogen number density $n_{\rm H}\approx 100\;\rm cm^{-3}$. Although
these simulations produce results similar to observations, it is still unclear
how the sub- to super-critical transition occurs, particularly since their
simulations use ideal MHD. Priestley et al. (2022) produce synthetic molecular
line observations from the prestellar cores in their non-ideal MHD
simulations, and argue that magnetically sub-critical cores can actually
appear super-critical in observations.

One significant limit of these works, however, is that they begin from an
arbitrarily prescribed uniform magnetic field morphology, rather than a
realistic galactic magnetic field structure. Given that observations indicate
that magnetic fields are continuously connected from galactic to cloud scales
(e.g., Li et al., 2014), this is a potentially significant issue. Moreover,
none of the works published to date have solved the radiative transfer problem
when producing H i 21 cm line synthetic observations, or fully solved the
problem of non-local thermodynamic equilibrium (non-LTE) excitation for the
Zeeman-sensitive OH lines. Since the regions of interest are precisely those
where H i self-absorption is expected to be strongest (as illustrated by the
Ching et al. 2022 HINSA observations), and where the density is well below the
OH critical density, this is another potentially significant limitation.

In this paper we study the transition between the sub-critical and super-
critical ISM using ideal MHD simulations that start from the scale of the
whole galaxy and then zoom in to reach maximum resolutions in dense regions up
to $\approx 0.4$ pc. These simulations allow us to for the first time combine
realistic galactic magnetic field structures with detailed dynamics on the
cloud scale in a single self-consistent calculation. We post-process our
simulations to derive realistic chemical compositions, and we then solve the
radiative transfer equation including non-LTE excitation effects along
different lines of sight to produce synthetic H i and OH observations. Our
goal is to understand the extent to which observed mass-to-flux ratio
accurately reflect the ratio of magnetic to gravitational energy in the 3D
volume they probe, and, to the extent that they do so, to understand how the
sub- to super-critical transition that appears in the observations arises in
the true 3D structure of the gas and magnetic field as they assemble to form
molecular clouds.

This paper is structured as follows: Section 2 summarises the numerical method
and simulation data. Section 3 describes our post-processing analysis method
to derive chemical composition, calculate radiative transfer, and produce the
synthetic observations. Section 4 presents the comparison between synthetic
observations and true 3D energy relations. Section 5 discusses the physical
meaning behind the results, and possible future work in this area. Section 6
concludes the work done in the paper.

##  2 Simulation

Our work begins from a simulation of an entire galaxy. To produce synthetic OH
Zeeman effect observations from this, we need to locate molecular clouds with
strong OH emission. Then we can produce synthetic H i observations in the
atomic regions surrounding the selected molecular clouds. To achieve this
goal, the scheme of the simulations in this work has 3 main steps: (1) carry
out MHD simulations of a full galaxy; (2) identify regions with strong OH
emission; (3) re-simulate the formation history of the selected regions with
higher resolution. We discuss both the numerical method and the detailed
scheme of simulations in this section.

###  2.1 Numerical method

Our full-galaxy simulations are based on the galactic magnetic field structure
simulation by Wibking & Krumholz (2022, hereafter WK22). Here we simply
summarise the physics and numerical methods, and highlight areas where our
treatment differs from that of WK22. For details on the remainder of the
methods, we refer readers to WK22. We solve the equations of ideal MHD with
the GIZMO code (Hopkins, 2015; Hopkins & Raives, 2016; Hopkins, 2016), which
uses a mesh-free, Lagrangian finite-mass Godunov method. The simulations
include a time-dependent chemistry network that tracks the abundances of H i,
H ii, He i, He ii, He iii, and free electrons, and uses these abundances to
calculate the atomic cooling rates due to hydrogen and helium. We also include
metal line cooling, using cooling rates interpolated from a prescribed table
as a function of density and temperature under the assumption of ionization
equilibrium and solar metallicity. The table is computed with CLOUDY (Ferland
et al., 1998), using the method of Smith et al. (2008) as applied in the
GRACKLE chemistry and cooling library (Smith et al., 2017).

For gas particles with density $\rho_{\rm g}$ larger that a threshold
$\rho_{\rm crit}$, we calculate a local star formation rate density

$\rho_{\rm SFR}=\epsilon_{\rm ff}\frac{\rho_{\rm g}}{t_{\rm ff}},$ (3)

where $\epsilon_{\rm ff}$ is the star formation efficiency, $\rho_{\rm g}$ is
the local gas density, and $t_{\rm ff}$ is the local gas free-fall time. As
summarized in the review by Krumholz et al. (2019), $\epsilon_{\rm ff}$ has a
mean value $\approx 0.01$ across a wide density range, which is also the value
we adopt. We set the value of $\rho_{\mathrm{crit}}$ to be approximately equal
to the Jeans density at our mass resolution; in the simulations of WK22, the
mass resolution $\Delta M=859.3$ M⊙, and we therefore set
$\rho_{\mathrm{crit}}=100$ H cm-3, but we increase the mass resolution and
thus $\rho_{\mathrm{crit}}$ in the zoom-in stages as described below. If
$\rho_{\rm g}>100\rho_{\rm crit}$, we increase $\epsilon_{\rm ff}$ to 1 in
order to limit the computational cost, since high densities require
correspondingly small time steps. This effectively sets an upper limit of
$100\rho_{\mathrm{crit}}$ on the gas density resolved by the simulation, since
gas particles with higher density will instantly be converted into stellar
particles. Gas particles for which $\rho_{\mathrm{SFR}}>0$ have probability

$P=1-\text{exp}(-\epsilon_{\rm ff}\Delta t/t_{\rm ff})$ (4)

of being converted to a star particle per time step of size $\Delta t$.
Stellar particles represent clusters that only interact with the galactic
environment via gravity and feedback.

While our treatment of star formation and cooling is identical to that of WK22
(other than our use of a higher $\rho_{\mathrm{crit}}$ at higher resolution),
our treatment of stellar feedback is necessarily slightly different because
our resolution will be high enough that stellar particles no longer represent
clusters massive enough that we expect them to sample the full initial mass
function (IMF). This matters because it means that we cannot compute stellar
feedback simply by integrating over the IMF. Instead, for each stellar
particle, we stochastically draw its stellar population from a Chabrier IMF
(Chabrier, 2005) through the stellar population synthesis code SLUG (da Silva
et al., 2012; Krumholz et al., 2015), following the approach to using SLUG as
a sub-grid model described by Fujimoto et al. (2018) and Armillotta et al.
(2019). We compute stellar evolution using the Padova stellar tracks (Bressan
et al., 2012), and derive stars’ ionising luminosities from their masses and
ages using the “starburst99” spectral synthesis method (Leitherer et al.,
1999). We determine which stars end their lives as type II supernovae, and the
stellar ages at which those supernovae occur, using the tabulated results of
Sukhbold et al. (2016). Our numerical treatment of ionisation and supernova
feedback is identical to that of WK22; our method in the present paper differs
only in how we determine the ionising luminosity, number, and timing of
supernovae produced by a given stellar particle.

###  2.2 Zoom-in method

The original WK22 simulation follows a Milky Way-like galaxy at a mass
resolution $\Delta M=859.3$ M⊙, insufficient to resolve the small, dense
regions from which the OH lines used in Zeeman effect measurements arise. In
order to capture these regions, we carry out two additional simulation stages
to zoom into molecular clouds from the galactic scale. For the first stage, we
use the $t=0.6$ Gyr snapshot in WK22’s MHD simulation as our initial
condition. We choose this snapshot because both the mass-weighted mean
magnetic field strength ($6\;\rm\mu G$) and total SFR (2 $\rm
M_{\odot}yr^{-1}$) are stable and close to the values observed for the
present-day Milky Way. WK22 explore the properties of the simulated galaxy at
this time in detail, so we will not discuss them further here.

As a first step toward increasing the resolution, at this point we switch the
feedback prescription over to the SLUG-based method described in Section 2.1
and continue the simulation at the same resolution for another 0.1 Gyr. During
this time the SFR undergoes a transient increase, since no supernovae occur
for $\approx 4$ Myr after the change in feedback prescription, but thereafter
the SFR stabilises at around $3\;\rm M_{\odot}$ yr-1 for the last 0.05 Gyr.
Once the simulation reaches 0.7 Gyr of evolution, our next zoom-in step is
that we replace each gas particle with 10 smaller particles, which are
randomly distributed around the center of mass of their parent particle with a
distance of 0.2 times the parental smoothing length. All other physical
parameters remain the same during the split, except that the particle mass is
divided by 10, and the smoothing length is divided by $10^{1/3}$. In this way
we increase the mean gaseous mass resolution to $\approx 90\;\rm M_{\odot}$.
During this stage we also increase $\rho_{\rm crit}$ to $10^{4}\;\rm
m_{H}/cm^{3}$, which is approximately the Jeans density at the mass resolution
we will reach during the next zoom-in stage; we use this value rather than the
Jeans density during this simulation stage to ensure that there are no
transients after we zoom in again. We continue this mid-resolution run for
another 21 Myr, to $t=0.721$ Gyr. We show the face-on gas surface density
$\Sigma_{\rm gas}$ at this time in Figure 1.

In order to select regions for further zoom-in, we use the chemistry post
analysis code DESPOTIC to derive the OH 1665 and 1667 MHz emission line
luminosities of each gas particle; we provide details of our method for doing
so in Section 3.1. For each gas particle we first determine its total
luminosity $L_{\rm OH}$ as the sum of these 2 lines, then randomly select 7
particles with $L_{\rm OH}>(2/3)\max\left(L_{\rm OH}\right)$, where the
maximum is computed over all gas particles. We then draw a sphere with a
radius of 0.3 kpc around each selected particle; we show these regions as
white circles in Figure 1. We choose this radius to be close to the gaseous
disc scale height 0.34 kpc. The free-fall times of the 7 spheres range from
1.2 Myr to 21 Myr, with a median value of 3.8 Myr. Therefore, simulating the
last 10 Myr of these regions will be sufficient to reconstruct their formation
history.

Figure 1: Face-on gas surface density map of the last snapshot ($t=0.721$ Gyr)
of the mid-resolution run. The 7 white circles, each of radius 0.3 kpc, are
the high OH luminosity regions we select for the highest resolution simulation
stage.

We record the particle IDs of all gas particles within the 7 spheres, then
find them in an earlier simulation snapshot 10 Myr before the one shown in
Figure 1 ($t=0.711$ Gyr). Using the method described above, we split these gas
particles into 10 again, reaching a mass resolution of 8.9 $\rm M_{\odot}$;
the corresponding Jeans length is $\approx 0.4$ pc, which we can use as an
estimate of the spatial resolution in the dense regions. We then re-simulate
the last 10 Myr ($t=0.711-0.721$ Gyr), so that the regions we will use in our
final analysis have all been simulated for $\approx 1$ free-fall time at our
highest resolution. We refer to this final stage as the high-resolution run.

To check if the simulation result is consistent after splitting selected gas
particles, we compare the face-on surface density maps of the mid- and high-
resolution runs at t = 0.721 Gyr for one example region in Figure 2. We see
that, as expected, the large-scale density structure is the same, but that the
higher resolution case shows more structure at high $\Sigma_{\rm gas}$.

To derive some basic physical parameters on the molecular cloud scale ($\sim$
100 pc), we draw (100 pc)3 cubes around the densest particle for each group of
split particles in the last high-resolution snapshot. This size is close to
the typical size of the largest molecular cloud complexes (Krumholz, 2014b).
In Table 1, we tabulate the total mass in the form of H i, $M_{\rm HI}$, total
mass in H2, $M_{\rm H_{2}}$, volume-weighted mean hydrogen number density,
$\bar{n}_{\rm H}$, and volume-weighted and mass-weighted mean magnetic field
strengths, $\bar{B}_{V}$ and $\bar{B}_{M}$, in each of these 7 cubes at
$t=0.721$ Gyr. We explain how we derive $M_{\rm HI}$ and $M_{\rm H_{2}}$ in
Section 3.1.

Figure 2: Face-on gas surface density maps of cloud 4 in the mid-resolution (top panel) and high-resolution (bottom panel) simulations at time $t=0.721$ Gyr. Cloud | $M_{\rm HI}$ | $M_{\rm H_{2}}$ | $\bar{n}_{\rm H}$ | $\bar{B}_{V}$ | $\bar{B}_{M}$
---|---|---|---|---|---
| ($10^{5}\;\rm M_{\odot}$) | ($10^{5}\;\rm M_{\odot}$) | $(\rm cm^{-3})$ | $(\rm\mu G)$ | $(\rm\mu G)$
1 | 0.35 | 0.23 | 5.0 | 2.7 | 5.1
2 | 0.26 | 0.54 | 5.8 | 1.4 | 4.9
3 | 0.35 | 0.72 | 6.6 | 1.2 | 9.9
4 | 0.55 | 0.75 | 8.3 | 1.2 | 7.9
5 | 1.1 | 1.4 | 16 | 15 | 22
6 | 1.0 | 3.0 | 20 | 7.2 | 20
7 | 1.6 | 10 | 53 | 17 | 48
Table 1: Physical parameters derived from (100 pc)3 cubes drawn around the
centres of mass of the high resolution regions at simulation time $t=0.721$
Gyr. $M_{\rm HI}$ and $M_{\rm H_{2}}$ are the total mass in the form of H i
and H2, respectively. $\bar{n}_{\rm H}$ is volume weighted mean hydrogen
number density. $\bar{B}_{V}$ and $\bar{B}_{M}$ refer to volume-weighted and
mass-weighted mean magnetic field strength, respectively.

##  3 Synthetic Observations

To study the mechanism behind the observed mass-to-flux ratio transition, we
need to produce synthetic observations using methods that match real
observational methods as closely as possible. This process consists of 3 main
steps: (1) determine the chemical state of each gas particle, and its
absorption and emission coefficients in the H i and Zeeman-sensitive OH lines,
respectively; (2) solve the equation of radiative transfer to obtain synthetic
spectra; (3) produce synthetic emission maps and B-field maps from the
spectra, using the same techniques used to analyse observed spectra. We now
proceed to describe these steps in detail.

###  3.1 Chemistry and excitation calculation

Since our simulation does not include molecular chemistry, and even
simulations that do include chemistry (e.g., Glover et al., 2010) do not self-
consistently compute non-LTE excitation (which is critical for the OH
emission), we post-process our simulations using the DESPOTIC astrochemistry
and radiative transfer code (Krumholz, 2014a). DESPOTIC derives the
equilibrium excitation, thermal, and chemical states of clouds from their
elemental abundances (assumed to be Solar for the present work), radiation
field environment, and physical properties. The thermal equilibrium
calculation includes cosmic ray heating, molecule line cooling, the grain
photoelectric effect, and collisional energy exchange between gas and dust,
while the chemical equilibrium calculation adopts the H-C-O chemical network
by Gong et al. (2017). The excitation calculation, which also determines
optical depths for line cooling, uses the large velocity gradient
approximation. We refer the readers to Krumholz (2014a) for more details on
DESPOTIC.

We perform the chemistry post-analysis following the methodology of Armillotta
et al. (2020). The first step is to generate tabulations of the H i and H2
mass fractions, $X_{\rm HI}$ and $X_{\rm H_{2}}$, the OH transition line
luminosities at 1665 and 1667 MHz, $L_{\rm OH,1665}$ and $L_{\rm OH,1667}$,
and the gas temperature,111Note that the temperature returned by the DESPOTIC
calculation does not match the temperature found by GIZMO because the GIZMO
calculation does not include molecular cooling, while DESPOTIC does; thus
DESPOTIC generally yields lower temperatures in dense, molecular regions. This
difference is unimportant for the dynamics, because in such cold regions
thermal pressure is generally an unimportant force, but it matters a great
deal for the excitation and radiative transfer calculations. $T$, as a
function of local velocity gradient, $dv/dr$, hydrogen number density, $n_{\rm
H}$, and hydrogen column density, $N_{\rm H}$; note that $n_{\mathrm{H}}$ and
$N_{\mathrm{H}}$ are the number of hydrogen nuclei per unit volume and per
unit area, so these quantities are invariant under different chemical
compositions. Since all 7 clouds are on spiral arms away from the galactic
center, we assume the interstellar radiation field (ISRF) strength to be
$\chi=1\;G_{0}$, where 1 $G_{0}$ corresponds to the solar neighbourhood ISRF
strength (Draine, 1978). We set the primary cosmic ray ionization rate to a
Solar neighbourhood-like value $\zeta=10^{-16}$ s-1. We run DESPOTIC on a grid
of models covering the parameter range $n_{\rm H}=10^{-2}-10^{6}\;\rm
cm^{-3}$, $N_{\rm H}=10^{19}-10^{25}\;\rm cm^{-2}$, and
$dv/dr=10^{-3}-10^{2}\;\rm km/s/pc$. The final grid is $21\times 16\times 15$
points in size in $(n_{\rm H},N_{\rm H},dv/dr)$ space, and is uniformly spaced
in logarithm. We have verified by visual inspection that quantities of
interest vary smoothly over our grid, so interpolation errors are small.

The second step is to assign values of $X_{\mathrm{HI}}$,
$X_{\mathrm{H_{2}}}$, $L_{\mathrm{OH,1665}}$, $L_{\rm OH,1667}$, and $T$ to
each gas particle in the simulations by performing trilinear interpolation in
the variables $dv/dr$, $n_{\rm H}$, and $N_{\rm H}$. We compute the hydrogen
number density for each particle as

$n_{\rm H}=\frac{X_{\rm H}\rho_{\rm g}}{m_{\rm H}},$ (5)

where $X_{\rm H}=0.71$ is the mass fraction of hydrogen assuming Solar
metallicity, $\rho_{\rm g}$ is the local gas density, and $m_{\rm
H}=1.67\times 10^{-24}$ g is the mass of the hydrogen nucleus. For $N_{\rm
H}$, we adopt the definition in Safranek-Shrader et al. (2017) as

$N_{\rm H}=n_{\rm H}L_{\rm shield},$ (6)

where $L_{\rm shield}$ is the shielding length. Safranek-Shrader et al. (2017)
show that the local approximation to the shielding length that most closely
reproduces the results of detailed ray-tracing radiative transfer calculations
is to set $L_{\mathrm{shield}}$ to the local Jeans length $L_{\rm J}$, with
the Jeans length computed using the smaller of the local gas temperature and
$T=40$ K. Finally, we calculate the velocity gradient as

$\frac{dv}{dr}=\|\nabla\times\textit{v}\|,$ (7)

where v is the local gas velocity. Our grid covers the full range of gas
particle properties in the high-resolution simulations, so extrapolation off
the grid is not necessary.

To illustrate the outcome of our interpolation procedure, in Figure 3 we plot
the distribution of gas properties along a 100 pc ray through the minimum of
gas temperature in cloud 4. This ray is along the x-axis as shown in Figure 2,
and we set its coordinates in the transverse direction. In Figure 3 we sample
the ray at 200 positions spaced at intervals of $\Delta L=0.5$ pc. We compute
the properties shown at each sampling point as a sum over nearby gas
particles, weighted by their smoothing kernels following the GIZMO kernel
weighting scheme described by Hopkins (2015). The upper panel shows the H i
number density $n_{\rm HI}$ (red), twice the H2 number density $2n_{\rm
H_{2}}$ (green), and the gas temperature $T$ (blue), while the bottom panel
shows the angle-averaged, frequency-integrated emissivity in the OH 18 cm
lines $j_{\rm OH,1665}$ (orange) and $j_{\rm OH,1667}$ (purple), together with
line-of-sight velocity $V_{\rm LOS}$ (black), where we take the LOS direction
to be along the ray. In these plots there is a clear transition from H i to H2
in the central cold region, and we see that OH emission is confined to this
region.

Figure 3: Gas properties along a 100 pc ray through the point of the lowest
gas temperature for cloud 4; the ray is parallel to the $x$-axis shown in
Figure 2. Top panel: number density of H nuclei $n_{\rm H}$ in the form of H i
(red) and H2 (green), together with local gas temperature $T$ (blue). Bottom
panel: frequency- and angle-integrated emissivity $j_{\rm OH,1665}$ (orange
line) and $j_{\rm OH,1667}$ (purple line) in the OH 18 cm lines, and line-of-
sight velocity $V_{\rm LOS}$ (black line).

###  3.2 Radiative transfer

Our next step is to compute observable spectra. For OH this is
straightforward, since Zeeman measurements are made using OH lines seen in
emission, and observations indicate that the OH lines used are invariably
optically thin. Thus we can compute the OH intensity as a function of velocity
simply by integrating the emissivities returned by the DESPOTIC calculation
along the line of sight. In practice, we do this by interpolating the particle
data onto a 2003 cube with a $\Delta L=0.5$ pc voxel size centred on the
densest particle; at each voxel we compute the OH emissivities in the 1665 and
1667 MHz lines by interpolating the particle data as described in Section 3.1
and used to construct Figure 3. We now consider lines of sight (LOSs) through
each pixel on the face of the cube in the $x$ and $y$ directions, i.e., the
two directions that lie within the galactic plane; thus we have a total of
$2\times 200^{2}$ pixels for each of our 7 clouds. For each pixel, we
determine its mass weighted mean line-of-sight velocity $\bar{V}_{\rm LOS}$,
and consider a grid of velocities over a range $\bar{V}_{\mathrm{LOS}}$ $\pm$
10 km/s with 0.1 km/s resolution. We then integrate the intensity in each
pixel and at each velocity to produce our final OH spectra.

The situation for H i is more complex, because H i Zeeman measurements are
generally made using absorption against background radio continuum sources
(e.g., Heiles & Troland, 2004), and comparing the on-source absorption profile
against an off-source emission profile. We therefore require a full radiative
transfer calculation for H i. Our first step in this calculation is to
interpolate the particle data onto a regular cube exactly as for OH; at each
pixel we have the H i and H2 number densities $n_{\rm HI}$ and $n_{\rm
H_{2}}$, along with the gas temperature $T$, gas velocity $\mathbf{v}$, and
the local sound speed $c_{s}=\sqrt{k_{B}T/m_{\rm H}}$. Then we compute the H i
emissivity $j_{\nu}$ and attenuation coefficient $\kappa_{\nu}$ as

$\displaystyle j_{\nu}$ $\displaystyle=$
$\displaystyle\frac{3}{16\pi}A_{ul}E_{ul}n_{\rm HI}\phi_{\nu}$ (8)
$\displaystyle\kappa_{\nu}$ $\displaystyle=$
$\displaystyle\frac{3A_{ul}hc\lambda_{ul}n_{\rm HI}}{32\pi k_{B}T}\phi_{\nu},$
(9)

where $\lambda_{ul}=21.12$ cm is the transition rest wavelength,
$E_{ul}=5.87\times 10^{-6}$ eV is the transition energy, $A_{ul}=2.88\times
10^{-15}\;\rm s^{-1}$ is the Einstein coefficient, and $\phi_{\nu}$ is the
line shape function, given by

$\phi_{\nu}=\frac{\lambda_{ul}}{\sqrt{2\pi}c_{s}}e^{-(V_{c}-V_{\rm
LOS})^{2}/c_{s}^{2}},$ (10)

where $V_{c}$ is the channel velocity and $V_{\mathrm{LOS}}$ is the component
of the gas velocity $\mathbf{v}$ parallel to the LOS.

To produce absorption spectra against background radio sources, we assume that
H i emission is negligible, and add up the optical depth along each LOS and
for each velocity channel as

$\tau_{\rm\nu}=\sum_{i}\kappa_{\nu,i}\Delta L,$ (11)

where $i$ denotes the pixel number, increasing from $i=1$ for the farthest
pixel to $i=200$ for the nearest one. To produce the corresponding H i
emission spectra for an off-source observation, we solve the radiation
transfer equation along the same lines of sight by iteratively evaluating

$I_{{\rm HI},\nu,i+1}=e^{-\kappa_{\nu,i}\Delta L}I_{{\rm
HI},\nu,i}+j_{\nu,\mathrm{i}}\Delta L,$ (12)

for each velocity channel, where $i$ again runs for 1 to 200. Our final
emission spectrum is the intensity $I_{\rm HI}$ that emerges from the last
pixel.

In Figure 4 we show an example result of applying this procedure to the LOS
shown in Figure 3; in the figure we plot the brightness temperature $T_{\rm
B}=\lambda_{ul}^{2}I_{\rm HI}/(2k_{B})$ (top) and absorption fraction
$e^{-\tau}$ (bottom) as functions of the channel relative velocity $\Delta
V_{\rm LOS}\equiv V_{c}-\bar{V}_{\rm LOS}$.

Figure 4: Top panel: the synthetic H i emission spectrum (off-source)
corresponding to the LOS illustrated in Figure 3. $T_{\rm B}$ is the
brightness temperature corresponding to the observed intensity, and $\Delta
V_{\rm LOS}$ is the difference between the velocity channel and the mass-
weighted mean $\bar{V}_{\rm LOS}$ of the whole LOS. The velocity resolution is
0.1 km/s. Bottom panel: absorption depth from the same LOS assuming the
presence of a background radio continuum emitter (on-source). The blue line is
the observed absorption fraction $e^{-\tau_{\rm obs}}$, while the red dash
lines show a decomposition of the absorption spectrum into Gaussian
components.

###  3.3 Emission and B-field maps

After solving equation 12 for all LOSs, we can integrate the derived emission
spectra over velocity to produce 2 maps of total H i 21 cm line intensity,
$I_{\mathrm{H\textsc{i}}}$, normal to $x$\- and $y$-axis, respectively. We
similarly produce maps of OH velocity-integrated intensity, $I_{\mathrm{OH}}$,
averaging the intensities of the 1665 and 1667 MHz lines. We compare the
resulting H i and OH total intensity maps to the true column density map for
cloud 4 (as shown in Figure 2) in Figure 5; the true column density
$N_{\mathrm{H}}$ we show is integrated over the same 100-pc lines of sight
used for the simulated emission maps. The three maps do not perfectly mimic
each other. The $N_{\rm H}$ and $I_{\rm HI}$ maps differ because low
$N_{\mathrm{H}}$ regions can be ionised and high $N_{\rm H}$ regions can be
optically thick and can be dominated by H2, all of which decrease the observed
intensity. The OH map strongly picks out only the densest region, and due to
varying excitation the OH intensity distribution does not perfectly mirror the
total density distribution even in regions where there is strong OH
emission.222Careful readers may note that the OH intensity $I_{\mathrm{OH}}$
shown in the bottom panel of Figure 5 can be slightly negative. This is a real
effect, which occurs because differential selection rules for the upper and
lower hyperfine states of the OH 1665 and 1667 MHz lines can drive the
excitation temperatures below the CMB temperature, so the line is seen in
absorption against the CMB. As a result the total intensity of the line, after
subtracting off the CMB, can be slightly negative along some lines of sight.

Figure 5: Maps from the x-axis projection of the (100 pc)3 cube extracted
around the densest gas particle in cloud 4. Top panel: total column density of
hydrogen nuclei $N_{\rm H}$. Middle panel: H i emission line velocity-
integrated intensity $I_{\rm HI}$. Bottom panel: total intensity $I_{\rm OH}$
in the OH 1665 and 1667 MHz lines; the quantity we plot is the average of the
two lines.

Our final task is to derive magnetic field measurements from the synthetic
spectra. We do so using three of the Zeeman techniques described in Section 1:
OH emission, H i absorption against a background radio source, and H i narrow
self-absorption (HINSA).

When one uses OH emission lines to measure the magnetic field strength through
the Zeeman effect, the emission regions that contribute more photons to the
observer will contribute more to the measured $B$-field strength. Since the
lines are optically thin and single-peaked, the contribution is simply linear
in the emission. Therefore, we approximate the magnetic field that would be
inferred from the OH lines as simply the $L_{\rm OH}$-weighted mean line of
sight magnetic field for both OH 18 cm lines. We then take the average of
these two values as the final result $B_{\rm OH}$. In Figure 6, we compare the
true, volume-weighed mean LOS magnetic field for cloud 4 projected along the
x-axis with maps generated from synthetic OH observations. For $B_{\rm OH}$ we
show only pixels where $I_{\rm OH}>I_{\rm OH,max}/100$ as a rough proxy for
where the OH lines would be bright enough to be detectable.

We next consider H i absorption against a background radio continuum source,
which is somewhat more complex since each H i absorption spectrum typically
contains several cold neutral medium (CNM) components with different $V_{\rm
LOS}$. To make separate measurements for each component, following the
approach of Heiles & Troland (2004), we perform least-squares fits of at most
5 Gaussians to the absorption spectrum. The optical depth of a single
component $\tau_{\rm comp,\nu}$ is

$\tau_{\rm comp,\nu}=\frac{\tau_{\rm
norm}}{\sqrt{2\pi\sigma_{v}^{2}}}\;e^{-(V_{c}-V_{\rm
comp})^{2}/(2\sigma_{v}^{2})},$ (13)

where $\tau_{\rm norm}$, $\sigma_{v}$ and $V_{\rm comp}$ are parameters to be
fit. We show the example of such a multi-component Gaussians as red dash lines
in the bottom panel of Figure 4. Due to computational resource limitations, we
do not perform this decomposition for all of our $5.6\times 10^{5}$ spectra;
instead, for each projection map we decompose 1000 randomly-selected LOSs,
chosen out of the set with $N_{\rm HI}\in(10^{19},\;10^{21.5})$ cm-2, the same
column density range observed in Heiles & Troland (2004).

Once we have decomposed the H i spectra, we next compute the synthetic Zeeman-
inferred magnetic field for each. This requires some subtlety, because
components in velocity space do not necessarily map onto real structures in
physical space, as is implicitly assumed by the decomposition procedure used
by Heiles & Troland (2004) and others. In such situations, the marginal
contribution of the magnetic field in a particular fluid element in space to
the magnetic field inferred for a particular Gaussian component will depend to
some extent on the details of the decomposition and fitting procedure.
However, generically we expect that the contribution of a particular voxel to
the magnetic field inferred for a particular component to scale with both the
total amount of absorption it contributes and with its marginal contribution
to that particular component. We there define a weight $w_{i}$ for each voxel
$i$’s contribution to a particular component of central velocity
$V_{\mathrm{comp}}$ and dispersion $\sigma_{v}$ as

$w_{i}=e^{-(V_{\rm LOS,i}-V_{\rm comp})^{2}/(2\sigma_{v}^{2})}\int
k_{\nu,i}d\nu,$ (14)

and the resulting synthetic $B$-field measurement for one Gaussian component
as

$B_{\rm HI,comp}=\frac{\sum_{i}w_{i}B_{\rm LOS,i}}{\sum_{i}w_{i}}.$ (15)

We only perform the calculations above for Gaussian components with
$\sigma_{v}<2$ km/s, to make sure that most of the component lies within the
spectral range we have extracted. From all selected LOSs we find 30853
components and derive $B_{\rm HI,comp}$ values for them.

Our final task is to make synthetic HINSA magnetic field measurements. Since
HINSA only occurs where there is cold H i ($T\lesssim 100$ K), we first
project only cold H i into column density maps, then find the pixel with the
highest column density, and determine the $B$-field in the cold H i $B_{\rm
HI,cold}$ from the 100-pc long LOS projected onto this pixel. From Figure 3 we
can see that the cold H i region has nearly constant $V_{\rm LOS}$, thus we
can assume the Doppler profile shape is constant for the whole cold H i
region, and determine $B_{\rm HI,cold}$ as

$B_{\rm HI,cold}=\frac{\sum_{\rm i,cold}B_{\rm LOS,i}\int
k_{\nu,i}d\nu}{\sum_{\rm i,cold}\int k_{\nu,i}d\nu},$ (16)

where $\sum_{\rm i,cold}$ means summation over pixels with $T\leq 100$ K.
After 2 projections per cube, we have 14 $B_{\rm HI,cold}$ results. Note that
we make only one measurement toward each cloud in order to mimic the approach
of the published HINSA measurement (Ching et al., 2022), which is typically
conducted along the line of sight through the HINSA maximum.

It is clear that the B-field strength and morphology highlighted by OH
concretely differ from the mean B-field over the 100-pc-long LOS, indicating a
possible inconsistency between the observed mass-to-flux ratio and the real
overall magnetic-gravity energy relation. Since for each LOS we make several
$B_{\rm HI}$ measurements from the fitted Gaussian components, individually,
and the selected LOSs do not cover the whole (100 pc)2 area, we can not show a
$B_{\rm HI}$ map similar to Figure 6.

Figure 6: Line-of-sight magnetic field maps from the $x$-axis projection of
the (100 pc)3 cube extracted around the densest gas particle in cloud 4. Top
panel: true, volume weighted mean line-of-sight magnetic field $B_{\rm vol}$
map. Bottom panel: magnetic field inferred from synthetic OH Zeeman
observations, $B_{\mathrm{OH}}$. Note that we only show $B_{\rm OH}$ for
pixels with total OH intensity $I_{\rm OH}>I_{\rm OH,max}/100$. Also note that
the colorbar range differs in the 2 panels.

##  4 Results

In this section, we first reproduce the observed $B$-$N$ relation with
synthetic observations in Section 4.1, then explore how well this relation
captures the true ratio of gravitational to magnetic energy from the full 3D
simulation data in Section 4.2.

###  4.1 Magnetic field versus column density

After making the synthetic observations, we can now compare magnetic field to
column density following Crutcher (2012), and look for the magnetic
criticality transition. The pixel size is 0.5 pc in both Figure 5 and Figure
6, which is close to the $3^{\prime}$ beam width of the Arecibo telescope when
applied to a cloud at a distance of $\sim$ 300 pc. Therefore, we simply treat
each pixel in the OH emission map as a synthetic beam, and directly take OH-
inferred magnetic field values $B_{\rm OH}$ from the pixel values derived in
Section 3.1. We similarly derive the corresponding column density from the
emission intensity maps. Assuming the gas is optically thin (as is assumed in
observations), we convert OH intensity into hydrogen column density as

$N_{\rm H,OH}=\frac{a\lambda^{3}I_{\rm OH}}{2k_{B}},$ (17)

where $\lambda$ is the emission line wavelength, $I_{\rm OH}$ is the velocity-
integrated intensity, and $a$ is an empirical constant. For the two OH 18 cm
lines, we have $a=2.12\times 10^{22}$ cm-2 / (K km s-1) for the 1665 MHz line
and $a=1.18\times 10^{22}$ cm-2 / (K km s-1) for the 1667 MHz line,
respectively (Troland & Crutcher, 2008).

We again follow observational procedure in deriving the column density related
to the H i absorption in each Gaussian component. For an isothermal cloud, the
H i column density is related to the velocity-integrated optical depth by
(e.g., Draine, 2011)

$N_{\rm HI}=\left(1.82\times 10^{18}\frac{\mbox{cm}^{-2}}{\mbox{K km
s}^{-1}}\right)T_{\rm s,comp}\int\tau_{\rm comp,\nu}\;dV,$ (18)

where $T_{\rm s,comp}$ is the spin temperature, and the integral is over the
velocity range of the absorption spectrum. With the combination of the off-
source emission spectrum $T_{\rm B}$ and the on-source absorption spectrum
$\tau$, we can determine $T_{\rm s,comp}$ as

$T_{\rm s,comp}=\frac{T_{\rm B,peak}}{1-\rm exp(-\tau_{\nu,\rm peak})},$ (19)

where the subscript peak indicates values measured at the center of the
Gaussian component where the line-of-sight velocity $V_{\rm LOS}=V_{\rm
comp}$. This gives column densities for the H i absorption against background
radio source measurements.

Finally, since we did not produce synthetic spectra for the HINSA
observations, we simply take the hydrogen column density in the form of H2 as
$N_{\rm H,cold}$; the approach matches that used in Ching et al. (2022),
though in their case they derive the H2 column from an independent estimate.
In our case we take the H2 column directly from the simulation data.

For consistent comparison with Crutcher (2012), we select H i synthetic beams
with a converted column density $N_{\rm HI}\in(10^{19},\;10^{21})\;\rm
cm^{-2}$, and OH synthetic beams with a converted column density $N_{\rm
H,OH}>10^{21.5}\;\rm cm^{-2}$. In Figure 7 we show the median and 16th to 84th
percentile range of measured magnetic field strength from synthetic H i
absorption and OH measurements in bins of column density 0.25 dex wide; the
blue band shows the former, the yellow band shows the latter. We also plot the
14 HINSA measurements toward the densest point in each synthetic map as
individual red points. The dashed green line shows the sub- / super-critical
transition, defined by

$\left|B_{\mathrm{LOS,crit}}\right|=3\pi\sqrt{\frac{G}{5}}\mu
m_{\mathrm{H}}N_{\rm H},$ (20)

where $\mu=1.4$ is the mean mass per H nucleon for gas with the usual
composition of $\approx 75\%$ H and $\approx 25\%$ He by mass. Finally, we
show the observations compiled by Crutcher (2012) as grey dots with error bars
for comparison.

The primary point to take from Figure 7 is that our synthetic observations do
an excellent job of reproducing the real observations. The median LOS field
strength is $\sim$ 1 $\mu$G for H i observations, and rises to $\sim$ 10
$\mu$G for OH observations. We can also see a clear transition in magnetic
criticality at column densities $N_{\rm H}\in(10^{21},\;10^{22})\;\rm
cm^{-2}$, such that the majority of the H i absorption measurements indicate
sub-criticality, while the OH measurements are almost exclusively in the
super-critical region. Again consistent with the result of Ching et al.
(2022), most HINSA synthetic observations also lie in the super-critical
region.

Figure 7: Magnetic field strengths inferred from synthetic Zeeman effect
observations using H i absorption against background radio sources (blue
band), OH emission (yellow band) and H i narrow self-absorption (red points),
as a function of inferred column density. For the yellow and blue bands, lines
with points show 50th percentile values of the $|B_{\rm LOS}|$ in bins of
column density 0.25 dex wide, while the band edges illustrate the 16th to 84th
percentile range. Grey points show real observations, from the compilation of
Crutcher (2012). The green dashed line represents $\mu_{B}=1$ (equation 20).
Measurements below this line correspond to super-critical conditions, while
those above are sub-critical.

###  4.2 Gravitational versus magnetic energy

Although the synthetic Zeeman effect measurements agree with previous
observational works, these may not reveal the true relation between
gravitational energy $E_{g}$ and magnetic energy $E_{B}$. To investigate
these, we compute the total gravitational potential and magnetic energies
within some volume of interest $V$ as

$\displaystyle E_{g}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{p}\sum_{q,q\neq p}\frac{GM_{p}M_{q}}{r_{pq}},$
(21) $\displaystyle E_{B}$ $\displaystyle=$
$\displaystyle\int_{V}\frac{B^{2}}{8\pi}\,dV.$ (22)

In the equations above, $G$ is the gravitational constant, $p$ and $q$ are
indices that run over all the GIZMO particles within $V$, $M_{p}$ and $M_{q}$
are the particle masses, $r_{pq}$ is the distance between particle $p$ and
particle $q$, and $B$ is the local total magnetic field strength. We evaluate
these energies as a function of size scale by identifying the densest gas
particle in each of our seven sample clouds, then drawing four spheres with
radii $r=10,25,50,$ and 100 pc. We choose these size scales to span the
atomic-to-molecular transition that is, based on the observational
diagnostics, closely related to the sub- / super-critical transition: for all
seven samples, the H2 mass fraction is $\approx 0.9$ for $r<10$ pc, but drops
to around $0.3$ for $90<r<100$ pc, indicating that our sample regions span
from the molecular core to the cloud edge. The 50 pc sphere is also close in
size to the 100 pc region onto which we interpolate onto a grid when computing
the synthetic Zeeman effect, though we will see below that the synthetic
observations are actually sensitive to considerably smaller scales than this.

In Figure 8 we plot the $E_{g}/E_{B}$ versus sphere radius for all seven
clouds. It is clear that $E_{g}/E_{B}>1$ for all clouds on most scales. The
only exception is cloud 1, which has $E_{g}/E_{B}\approx 0.6$ on $r=25$ and 50
pc scales, which may be related to the fact that it has the smallest total
mass in our sample; however, even these sub-critical regions are embedded
within a larger 100 pc-scale super-critical one. Therefore, from the 3D
information, we conclude that all seven clouds are already super-critical on
scales where the H i synthetic observations seem to indicate that they are
sub-critical. Using the full 3D information, we find no transition from a sub-
to a super-critical state across the H i \- H2 transition. This suggests that
the observed transition in magnetic criticality across the column density is
most likely a result of inconsistency between 2D projected observations and
real 3D structures. We explore the nature of this mismatch in the next
section.

Figure 8: Ratios of gravitational energy $E_{g}$ to magnetic energy $E_{B}$
measured in spheres with different radii $r$, centred on the densest point in
each sample cloud.

##  5 Discussion

We have seen that, although both actual Zeeman measurements and our synthetic
versions of them indicate the existence of a transition between sub- and
super-critical states as material transitions from H i to H2, no such
transition is actually found in our simulations; instead, even on 100 pc
scales around giant molecular clouds where the mass is predominantly atomic,
clouds are super-critical. Our goal in this discussion is to understand why
the H i Zeeman observations seem to indicate the presence of a sub-critical
state where none actually exists. We find that this sub-critical illusion has
two contributing factors. First, in previous observational work, the column
density $N$ associated with a particular H i Zeeman measurement is directly
converted from the observed absorption line depth, which measures only H i
material, and misses other chemical states of gas. Second, Zeeman effect
measurements only provide information on the magnetic field in regions of
strong absorption (or emission for OH), and such regions may not be
representative of the true magnetic field strength, leading to a biased
estimate of $E_{g}/E_{B}$.

Starting with the first of these factors, to study how accurate the $N_{\rm
HI}$ synthetic measurement is, we compare it with results directly derived
from the simulation data. For each Gaussian component for which we make a
synthetic Zeeman measurement, we first find the related line of sight, then
compute the density $n_{\rm HI}$ at positions along that LOS where the LOS gas
velocity $V_{\rm LOS}$ lies in the range $V_{\rm comp}\;\pm\sigma_{v,\rm
comp}$, i.e., we find all the fluid elements whose velocities lie within one
standard deviation of the central velocity of each component. We then compute
the column density $N_{\rm HI,1\sigma_{v}}$ of these fluid elements. We then
define the intrinsic H i column density related to the Gaussian component as
$N_{\rm HI,sim}=N_{\rm HI,1\sigma_{v}}/0.683$, where 0.683 is the fraction of
the mass that lies within $\pm 1\sigma$ of the centre for a Gaussian
distribution. We also use an identical procedure to derive the total hydrogen
column density $N_{\rm H,sim}$ by projecting the total density of H nuclei
$n_{\rm H}$; $N_{\rm H,sim}$ differs from $N_{\rm HI,sim}$ in that the former
includes all hydrogen regardless of its chemical state, while the latter
includes only hydrogen that is chemically in the form of H i.

We can compare these true column densities to the columns $N_{\mathrm{HI}}$
derived from the 21 cm spectra using observational methods, as described in
Section 4.1, equation 15 and equation 19. The ratios $(N_{\rm HI}/N_{\rm
HI,sim})$ and $(N_{\rm HI}/N_{\rm H,sim})$ then tell us how far the
observational measurements lie from the H i mass and the total mass. In Figure
9 we plot the 16th, 50th and 84th percentiles of $(N_{\rm HI}/N_{\rm HI,sim})$
and $(N_{\rm HI}/N_{\rm H,sim})$, binning the data by observationally-
estimated column density $N_{\mathrm{HI}}$ in the same pattern as in Figure 7.
From the figure we can see that the observational procedure recovers about
$\sim 40\%$ of H i mass, but only $\sim 10\%$ of the total mass. The reason is
that most mass along the line of sight is in the form of H2 and ionised
hydrogen. In retrospect this result is not surprising: observational probes of
H i around molecular clouds do not represent a randomly-selected sample of
lines of sight through the ISM. Instead, they are probing regions where it is
reasonable to believe that dark gas composed mostly of H2 might exist. This
result shows that observationally-inferred mass-to-flux ratios likely suffer
from underestimates of $N_{\rm HI}$ due to chemical composition variations.

Figure 9: The blue band shows the ratio of inferred to true H i column
density, $(N_{\rm HI}/N_{\rm HI,sim})$, while the grey band shows the ratio of
inferred H i column density to true column density of H nuclei in all chemical
states, $(N_{\rm HI}/N_{\rm H,sim})$. For each band, lines with dots show the
50th percentile values in bins of observationally-inferred H i column density
$N_{\mathrm{HI}}$, while the band edges illustrate the 16th and 84th
percentiles.

The second factor contributing to the inconsistency between the observed and
true levels of magnetic criticality we identify in Section 4 is the limited
spatial scale of the magnetic field probed by Zeeman effect measurements. From
equation 15 we learn that $B_{\rm HI}$ measurements are dominated by regions
with large absorption coefficients, which may only reveal only a small part of
the overall magnetic field structure. To quantify this effect, we compute an
effective length scale probed by H i absorption measurements as follows. For
each LOS, we assign the $i$th sample point on that LOS a position
$x_{i}=i\Delta L$. We then define the length scale of H i the B-field
measurement $l_{\rm HI}$ to be the weighted standard deviation of $x_{i}$,
where the weight $w_{i}$ is defined by that sample point’s contribution to the
synthetic Zeeman measurement, as defined in equation 14. Quantitatively, we
compute $l_{\rm HI}$ as

$\displaystyle x_{\rm mean}$ $\displaystyle=$
$\displaystyle\frac{\sum_{i}w_{i}x_{i}}{\sum_{i}w_{i}},$ (23) $\displaystyle
l_{\rm HI}$ $\displaystyle=$
$\displaystyle\left[\frac{\sum_{i}w_{i}(x_{i}-x_{\rm
mean})^{2}}{\sum_{i}w_{i}}\right]^{1/2},$ (24)

where the sum runs over all the sample points $i$ along a given line of sight.
Intuitively, $l_{\rm HI}$ calculated in this way quantifies the length along
the line of sight of the region that contributes to a Zeeman measurement; for
uniform $w_{i}$, given our 100 pc sampling region, $l_{\rm HI}$ equals 57.7
pc, while the smallest $l_{\mathrm{HI}}$ that we can resolve given our
sampling, corresponding to a case where the region that produces the
absorption is captured by exactly two sample points that both contribute
equally, is 0.25 pc. Performing this calculation over all LOS’s, we find that
the 16th, 50th, and 84th percentiles of $l_{\rm HI}$ are 4.11 pc, 11.4 pc, and
29.3 pc, respectively. Therefore, our synthetic H i observations generally
measure magnetic fields in regions of strong H i absorption whose
characteristic size scale is $\sim$ 10 pc. We have perform a similar
calculation on OH emission as in equation 24, but we change the weight
function to $L_{\rm OH}$, and derive the characteristic size scale of OH
Zeeman effect measurement as $\sim$ 3 pc.

To see why selecting H i-dominated regions of this size scale might produce a
bias, for each (100 pc)3 cube we identify all sample points where at least 75%
of hydrogen atoms are in the form of H i; then we randomly select 20 such
points that are at least 10 pc from each other, and draw $r=10$ pc spheres
around them. We determine the $E_{g}/E_{B}$ ratio for these spheres, yielding
a total of 140 samples – 20 from each of the seven sample clouds. For these
samples we find a median $E_{g}/E_{B}=0.25$, and the ratio is less than unity
for 63% of them. This result indicates that $\sim 10$-pc scale regions of H
i-dominated gas are likely to be locally sub-critical, which is consistent
with the observations of Heiles & Troland (2004) and others. However, this is
an illusion created by sampling unrepresentative regions: when considering the
whole cloud from which these H i regions are drawn on 100 pc scales, or when
selecting regions centred on the density maximum rather than maxima of H i
absorption, we find that the gas is super-critical, as shown in Figure 8.

This analysis leads us to propose a scenario for the relation between gravity
and magnetic fields in molecular clouds and their neutral hydrogen
surroundings illustrated in Figure 10. For a cloud $\sim$ 100 pc in size, its
dense molecular core is super-critical, and is revealed as such when it is
measured using OH emission lines. Small H i sub-regions in the cloud outskirts
are generally locally sub-critical, in the sense that their magnetic energy is
larger than their gravitational self-energy, consistent with the H i
absorption observations. However, even for regions that are genuinely sub-
critical on $\sim 10$ pc scales, if we consider the large scale cloud
structure as a whole, the molecular mass ignored in H i observations can
greatly increase the total gravitational energy and make the whole cloud
super-critical. Magnetic fields may be able to locally support H i regions
against gravity, but cannot resist gravitational collapse on the large scale.
One can make a useful analogy here to the competition between thermal and
gravitational energy in the envelope of the Sun – if one considers only the
outer layers of the Sun, say the outer half in radius, its thermal energy is
far greater than its gravitational binding energy. However, this is a
misleading statistic, because the dominant source of gravitational binding for
the outer 20% of the Sun is not its self-energy, but instead the gravitational
attraction provided by the remaining 80% of the mass. A measurement that is
only sensitive to the ratio of thermal to gravitational energy in the envelope
will completely miss this effect.

In this view, the apparent transition between sub-critical and super-critical
gas seen in Zeeman measurements, and described in Crutcher (2012), is a
confusion caused by mixing local and cloud-scale energy relations. Correct
recovery of the level of magnetic criticality requires the ability to measure
both the magnetic field and the true column across a whole molecular cloud.
Zeeman measurements do not provide this ability.

Figure 10: Schematic view of regions detected by Zeeman effect observation in
the H i cloud. The dashed circle labeled as H i A.C. is the H i sub-region
measured by an absorption component (A.C.) in observation. Note that the
length scale of different regions is not strictly proportional.

##  6 Conclusions

This work aims at understanding the apparent transition in magnetic
criticality between H i\- and H2-dominated regions found in Zeeman effect
observations. Such study requires a combination of high-resolution simulations
of the molecular cloud formation process and synthetic observations that
capture the complexities of Zeeman measurements in different chemical tracers,
following real observational procedures as closely as possible. To this end,
we perform the first MHD simulation of molecular cloud formation out of
realistic full-galactic magnetic field structure. Starting from the low-
resolution galactic simulation, we zoom into 7 molecular clouds, then trace
and re-simulate their formation history at high resolution ($\sim 0.4$ pc). We
then post-process the simulations to derive realistic chemical compositions,
gas temperatures, and non-LTE level populations, yielding physically accurate
calculations of emissivities and absorption coefficients in the most commonly-
used Zeeman-sensitive lines. We use these to solve the radiative transfer
equation along 100-parsec-long lines of sight to produce synthetic H i 21 cm
line absorption and emission spectra, which mimic the absorption / emission
strategy used in real H i Zeeman measurements, as well as synthetic OH
emission and H i narrow self-absorption (HINSA) observations. We process these
results using the same procedures used in observations to produce realistic
synthetic measurements of magnetic field strengths and column densities in
different tracers.

We find that our synthetic $B$ vs. $N$ relations are quantitatively consistent
with the observed trends summarised by Crutcher (2012): molecular clouds
appear to be sub-critical in H i observations, and become super-critical in OH
observations; the transition of criticality happens in the hydrogen column
density range of $N_{\rm H}\in[10^{21},10^{22}]$ cm-2. HINSA measurements
probe this column density regime, and generally indicate that clouds are
super-critical.

Such results, however, may not reveal the real ratio between gravitational
energy and magnetic energy $E_{g}/E_{B}$. We compare out synthetic Zeeman
measures to ratios of $E_{g}/E_{B}$ computed directly from the 3D simulation
data, evaluated over a range of size scales and for a sample of different
molecular clouds. On scales from 10 pc to 100 pc we find that the real ratio
$E_{g}/E_{B}>1$ for most clouds, which indicates that gravity dominates over
magnetic fields. This applies even when averaging over size scales large
enough that the material we are probing is largely atomic, indicating that the
H i Zeeman observations do not accurately reflect the true dynamical state of
the gas.

We identify two main contributors to this inconsistency. First, column
densities used in H i Zeeman observations only measures the hydrogen mass in
the neutral atomic form, which underestimates the true total column density
along lines of sight where the chemical composition is not uniformly atomic
hydrogen. Second, the regions probed by H i absorption Zeeman measurements
capture density and magnetic field structure in $\sim$ 10 pc-scale regions on
the outskirts of molecular clouds. Although these local H i sub-regions may be
sub-critical as observed, they are nonetheless part of larger structures that
are still super-critical once we include the gravitational influence of the
central dense molecular region. We conclude that accurate observational
determination of the relative importance of magnetic fields and gravity in the
molecular cloud formation process requires careful attention to ensure that
the observations accurately capture all the self-gravitating mass, and that
the regions that the observations capture are representative of entire clouds,
rather than small sub-regions within them.

## Acknowledgements

The authors acknowledge Prof. Di Li, Dr. Tao-chung Ching, and Dr. Hiep Nguyen
for suggestions on synthetic observations. We further acknowledge funding from
the Australian Research Council through its Future Fellowship and Laureate
Fellowship funding schemes, awards FT180100375 and FL220100020, and high-
performance computing resources provided by the Australian National
Computational Infrastructure (grants jh2 and ek9) through the National and ANU
Computational Merit Allocation Schemes, and by the Leibniz Rechenzentrum and
the Gauss Centre for Supercomputing (grant pr32lo).

## Data Availability

The data underlying this article will be shared upon reasonable request to the
corresponding author.

## References

  * Armillotta et al. (2019) Armillotta L., Krumholz M. R., Di Teodoro E. M., McClure-Griffiths N. M., 2019, MNRAS, 490, 4401
  * Armillotta et al. (2020) Armillotta L., Krumholz M. R., Di Teodoro E. M., 2020, MNRAS, 493, 5273
  * Bertram et al. (2012) Bertram E., Federrath C., Banerjee R., Klessen R. S., 2012, MNRAS, 420, 3163
  * Bourke et al. (2001) Bourke T. L., Myers P. C., Robinson G., Hyland A. R., 2001, ApJ, 554, 916
  * Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  * Chabrier (2005) Chabrier G., 2005, in Corbelli E., Palla F., Zinnecker H., eds, Astrophysics and Space Science Library Vol. 327, The Initial Mass Function 50 Years Later. p. 41 (arXiv:astro-ph/0409465), doi:10.1007/978-1-4020-3407-7_5
  * Ching et al. (2022) Ching T. C., Li D., Heiles C., Li Z. Y., Qian L., Yue Y. L., Tang J., Jiao S. H., 2022, Nature, 601, 49
  * Crutcher (2012) Crutcher R. M., 2012, ARA&A, 50, 29
  * Crutcher et al. (1999) Crutcher R. M., Troland T. H., Lazareff B., Paubert G., Kazès I., 1999, ApJ, 514, L121
  * Draine (1978) Draine B. T., 1978, ApJS, 36, 595
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium. Princeton University Press: Princeton, NJ
  * Falgarone et al. (2008) Falgarone E., Troland T. H., Crutcher R. M., Paubert G., 2008, A&A, 487, 247
  * Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761
  * Fujimoto et al. (2018) Fujimoto Y., Krumholz M. R., Tachibana S., 2018, MNRAS, 480, 4025
  * Glover et al. (2010) Glover S. C. O., Federrath C., Mac Low M., Klessen R. S., 2010, MNRAS, 404, 2
  * Gong et al. (2017) Gong M., Ostriker E. C., Wolfire M. G., 2017, ApJ, 843, 38
  * Heiles & Troland (2004) Heiles C., Troland T. H., 2004, ApJS, 151, 271
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53–110
  * Hopkins (2016) Hopkins P. F., 2016, MNRAS, 462, 576
  * Hopkins & Raives (2016) Hopkins P. F., Raives M. J., 2016, MNRAS, 455, 51
  * Ibáñez-Mejía et al. (2021) Ibáñez-Mejía J. C., Mac Low M.-M., Klessen R. S., 2021, arXiv e-prints, p. arXiv:2108.04967
  * Krumholz (2014a) Krumholz M. R., 2014a, MNRAS, 437, 1662
  * Krumholz (2014b) Krumholz M. R., 2014b, Physics Reports, 539, 49–134
  * Krumholz et al. (2015) Krumholz M. R., Fumagalli M., da Silva R. L., Rendahl T., Parra J., 2015, MNRAS, 452, 1447
  * Krumholz et al. (2019) Krumholz M. R., McKee C. F., Bland-Hawthorn J., 2019, ARA&A, 57, 227
  * Leitherer et al. (1999) Leitherer C., et al., 1999, ApJS, 123, 3
  * Li et al. (2014) Li H. B., Goodman A., Sridharan T. K., Houde M., Li Z. Y., Novak G., Tang K. S., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 101 (arXiv:1404.2024), doi:10.2458/azu_uapress_9780816531240-ch005
  * Mouschovias & Spitzer (1976) Mouschovias T. C., Spitzer L. J., 1976, ApJ, 210, 326
  * Priestley et al. (2022) Priestley F. D., Yin C., Wurster J., 2022, MNRAS, 515, 5689
  * Safranek-Shrader et al. (2017) Safranek-Shrader C., Krumholz M. R., Kim C.-G., Ostriker E. C., Klein R. I., Li S., McKee C. F., Stone J. M., 2017, MNRAS, 465, 885
  * Shu (1977) Shu F. H., 1977, ApJ, 214, 488
  * Shu et al. (1987) Shu F. H., Adams F. C., Lizano S., 1987, ARA&A, 25, 23
  * Smith et al. (2008) Smith B., Sigurdsson S., Abel T., 2008, MNRAS, 385, 1443
  * Smith et al. (2017) Smith B. D., et al., 2017, MNRAS, 466, 2217
  * Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H. T., 2016, ApJ, 821, 38
  * Thompson et al. (2019) Thompson K. L., Troland T. H., Heiles C., 2019, ApJ, 884, 49
  * Troland & Crutcher (2008) Troland T. H., Crutcher R. M., 2008, ApJ, 680, 457
  * Troland et al. (2016) Troland T. H., Goss W. M., Brogan C. L., Crutcher R. M., Roberts D. A., 2016, ApJ, 825, 2
  * Vishniac & Lazarian (1999) Vishniac E. T., Lazarian A., 1999, ApJ, 511, 193
  * Wibking & Krumholz (2022) Wibking B. D., Krumholz M. R., 2022, MNRAS, stac2648
  * da Silva et al. (2012) da Silva R. L., Fumagalli M., Krumholz M., 2012, ApJ, 745, 145