# Turbulent Diffuse Molecular Media with Non-ideal Magnetohydrodynamics and
Consistent Thermochemistry: Numerical Simulations and Dynamic Characteristics

Nannan Yue The Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing 100871, China  Lile Wang The Kavli Institute for Astronomy
and Astrophysics, Peking University, Beijing 100871, China  Department of
Astronomy, School of Physics, Peking University, Beijing 100871, China  Thomas
Bisbas Zhejiang Laboratory, Hangzhou 311121, China  Donghui Quan Zhejiang
Laboratory, Hangzhou 311121, China  Di Li National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100101, China  Zhejiang Laboratory,
Hangzhou 311121, China Lile Wang<EMAIL_ADDRESS>
###### Abstract

Turbulent diffuse molecular clouds can exhibit complicated morphologies caused
by the interactions among radiation, chemistry, fluids, and fields. We
performed full 3D simulations for turbulent diffuse molecular interstellar
media, featuring time-dependent non-equilibrium thermochemistry co-evolved
with magnetohydrodynamics (MHD). Simulation results exhibit the relative
abundances of key chemical species (e.g., $\mathrm{C}$, $\mathrm{CO}$,
$\mathrm{OH}$) vary by more than one order of magnitude for the “premature”
epoch of chemical evolution ($t\lesssim 2\times 10^{5}~{}\mathrm{yr}$).
Various simulations are also conducted to study the impacts of physical
parameters. Non-ideal MHD effects are essential in shaping the behavior of
gases, and strong magnetic fields ($\sim 10~{}\mu{\rm G}$) tend to inhibit
vigorous compressions and thus reduce the fraction of warm gases ($T\gtrsim
10^{2}~{}\mathrm{K}$). Thermodynamical and chemical conditions of the gas are
sensitive to modulation by dynamic conditions, especially the energy injection
by turbulence. Chemical features, including ionization (cosmic ray and diffuse
interstellar radiation), would not directly affect the turbulence power
spectra. Nonetheless, their effects are prominent in the distribution profiles
of temperatures and gas densities. Comprehensive observations are necessary
and useful to eliminate the degeneracies of physical parameters and constrain
the properties of diffuse molecular clouds with confidence.

Interstellar medium(847), Molecular clouds(1072), Astrochemistry(75),
Hydrodynamical simulations(767), Magnetohydrodynamical simulations(1966)

##  1 Introduction

Molecular clouds are an essential phase of the interstellar medium (ISM),
whose diffuse phase consists mainly of relatively dense ($\gtrsim
10^{2}~{}\mathrm{cm}^{-3}$), cool ($\lesssim 10^{2}~{}\mathrm{K}$), and
predominantly molecular gases (Draine, 2011). Prominent morphological
features, including filaments, bundles, cloud cores, and occasionally star-
induced bubbles, are shaped by the turbulent fluid motions interacting with
magnetic fields and radiation in these clouds (McKee & Ostriker, 2007;
Hennebelle & Falgarone, 2012). Understanding the characteristics and evolution
tracks of the diffuse molecular clouds is crucial for complete and
comprehensive physical pictures of the evolution of the ISM in general. The
diffuse regions with densities $\rho\lesssim
10^{3}~{}m_{p}~{}\mathrm{cm}^{-3}$, which seem less relevant to subsequent
star formation processes, also draw researchers’ attention by exibiting
relatively complicated thermochemical details, including exitation of ro-
vibrational molecular species in these relatively cold (typically $T\lesssim
10^{2}~{}\mathrm{K}$) regions indicated by existing and incoming observations
(e.g. Nehmé et al., 2008; Goldsmith et al., 2010; Langer et al., 2010; Ingalls
et al., 2011; Rice et al., 2018).

Multiple numerical studies have been conducted regarding interstellar
molecular gas’s morphological and dynamical characteristics. Pioneered by,
e.g., Mac Low et al. (1998); Mac Low (1999), a series of studies on molecular
clouds using 3D hydrodynamic or magnetohydrodynamic (MHD) simulations have
focused on the morphologies and statistical characteristics of decaying or
quasi-steady turbulence structures. Plenty of simulation works have assumed
isothermal equations of states, even some of them have adopted ray-tracing or
diffusive radiative transfer (Raskutti et al., 2016, 2017; Kim et al., 2018,
2019). Nonetheless, consistent cooling and heating processes are considered
essential to the reaction rates and excitation states of molecules, as well as
the dissipation conditions for the turbulence scaling laws (Neufeld et al.,
1995; Jappsen et al., 2007; Glover & Jappsen, 2007; Tilley & Balsara, 2011;
Qian et al., 2018; Whitworth & Jaffa, 2018). Simulating the turbulent media
with such global prescriptions, rather than real-time local calculations,
risks obscuring the interactions between thermochemistry and dynamics, leading
to potentially unreliable consequences in dynamics and observational
predictions.

Another critical factor is the non-ideal MHD effects. The weak ionization in
molecular clouds leads to finite tensorial conductivities, exhibiting
deviations from the flux-freezing ideal MHD. The tensor field of these
magnetic diffusivities depends on the spatial distributions of plasma
temperatures, free electrons, positive ions, and charged dust particles.
Because of the lack of detailed thermochemical information, most simulation
works tend to prescribe, rather than calculate in real-time, a set of
conductivity or diffusivity parameters (Downes & O’Sullivan, 2009; Tilley &
Balsara, 2011; Jones & Downes, 2011, 2012). Even in the star formation
simulations, magnetic diffusivities play crucial roles, Wurster et al. (2016)
adopted significantly simplified equilibrium chemistry for charge carriers and
ignored thermodynamics. However, the complications in charge-related reactions
(especially those involving molecules and dust grains) and radiation cooling
still necessitate real-time non-equilibrium thermochemistry despite the high
computational costs, as have been identified in multiple studies on
protoplanetary disks and molecular clouds (Xu & Bai, 2016; Wang et al., 2019;
Clark et al., 2019).

In order to address these critical issues for better intrinsic consistency and
predictability for observations, we adopt the GPU-accelerated co-evolving non-
ideal MHD-thermochemistry system that has been used in simulating
protoplanetary disks (e.g. Wang et al., 2019) and predicting exoplanetary
atmosphere observations (Wang & Dai, 2021). There are several recent
simulation studies that are conceptually similar. Works by e.g. Wurster et al.
(2019); Wurster & Rowan (2024) have delved into the use of magnetic
diffusivity profiles derived from chemical equilibrium to investigate the
formation of stellar cores and protostellar disks in relatively dense
molecular environments. Their findings underscore the critical role of
magnetic diffusivities in addressing the magnetic braking issue. Other
studies, such as Franeck et al. (2018); Clark et al. (2019); Lupi et al.
(2021), have examined various facets of molecular gases through the lens of
non-equilibrium chemistry, focusing on the evolution of observables and tracer
species. However, these studies have not taken into account the non-ideal MHD
effects. Our approach aims to bridge this gap by conducting a comprehensive
series of calculations within a dynamically quasi-steady turbulent framework.
We will employ consistent magnetic diffusivities that are grounded in more
extensive, non-equilibrium thermochemical networks in realtime. This endeavor
aims at providing a more thorough exploration covering the lifecycle of
typical diffuse molecular ISM, which spans approximately $\gtrsim
10^{6}~{}\mathrm{yr}$. Our goal is to reflect the long-term evolution of
thermochemical processes, thereby offering a more nuanced and accurate
representation of the ISM’s behavior over extended periods. We will portray
diffuse molecular clouds from diverse physical perspectives and systematically
explore the space spanned by key physical parameters, especially the
interstellar radiation field (ISRF), the cosmic rate (CR) ionization rate, the
mean-field magnetic intensity, the specific dust cross-section, metallicity,
and the turbulence energy injection rate. With detailed real-time
thermochemical networks evolved and recorded, extended explorations on the
observables and characterizations can also stem from these simulations in
subsequent future works.

This paper is structured as follows. Section 2 briefly introduces the methods
used in the simulations and describes the setups of several representative
simulations. Setups and results of the fiducial simulation are elaborated in
Section 3, and various models over the physical parameter space are explored
in Section 4. We discuss the possible future explorations based on current
results and conclude the paper in Section 5.

##  2 Methods and Setups

The computational schemes summarized in this section are primarily similar to
Wang et al. (2019), with modifications particular to molecular clouds. We will
also discuss the setups of numerical models in this paper.

###  2.1 Non-ideal Magnetohydrodynamics

The ionization of interstellar media, including the diffuse molecular gas
explored in this paper, is always significantly non-zero. This fact makes
their morphologies, evolution tracks, and thermochemistry profiles susceptible
to different interactions between fluids and fields. We use the grid-based
higher-order Godunov MHD code `Athena++` (Stone et al., 2020) on uniform
Cartesian grids, with HLLD approximate Riemann solver and piecewise linear
method (PLM) reconstruction. The MHD equations are solved in the conservative
form:

$\begin{split}&\partial_{t}\rho+\nabla\cdot(\rho\mathbf{v})=0\ ;\\\
&\partial_{t}(\rho\mathbf{v})+\nabla\cdot\left(\rho\mathbf{v}\mathbf{v}-\dfrac{\mathbf{B}\mathbf{B}}{4\pi}+P_{\mathrm{tot}}\mathbf{I}\right)=-\nabla\Phi\
;\\\
&\partial_{t}\mathbf{B}=\nabla\times(\mathbf{v}\times\mathbf{B}-c\mathbf{E}^{\prime})\
;\\\
&\partial_{t}\epsilon+\nabla\cdot\left[\left(\epsilon+P_{\mathrm{tot}}\right)\mathbf{v}-\dfrac{(\mathbf{B}\cdot\mathbf{v})\mathbf{B}}{4\pi}\right]=S\
,\end{split}$ (1)

where $\rho$, $\mathbf{v}$ and $p$ are the gas density, velocity and gas
thermal pressure respectively, $\mathbf{B}$ is the magnetic field,
$P_{\mathrm{tot}}\equiv p+B^{2}/(8\pi)$ is the total pressure, $\Phi$ is the
gravitational potential, $\mathbf{I}$ is the identity tensor, and

$\epsilon\equiv\dfrac{p}{\gamma-1}+\rho\left(\dfrac{v^{2}}{2}+\Phi\right)+\dfrac{B^{2}}{8\pi}$
(2)

is the total energy density. The extra source term $S$ in the energy equation
reflects non-adiabatic processes that affect gas energy, which is calculated
separately from the MHD solver in an operator-splitting manner (see §2.2). The
adiabatic index $\gamma$ is evaluated based on the real-time chemical
components in each simulation zone. Non-ideal MHD effects are also considered
by including an electric field in the local rest frame of fluids,

$\mathbf{E}^{\prime}=\dfrac{4\pi}{c^{2}}\left(\eta_{\mathrm{O}}\mathbf{J}+\eta_{\mathrm{H}}\mathbf{J}\times\mathbf{b}+\eta_{\mathrm{A}}\mathbf{J}_{\bot}\right)\
,$ (3)

where $\mathbf{b}\equiv\mathbf{B}/B$ is a unit vector along the local
direction of $\mathbf{B}$, $\mathbf{J}=c\nabla\times\mathbf{B}/4\pi$ is the
current density, and
$\mathbf{J}_{\bot}\equiv\mathbf{b}\times(\mathbf{J}\times\mathbf{b})$ is the
component of $\mathbf{J}$ perpendicular to the local magnetic field. Different
diffusivity components, including Ohmic ($\eta_{\mathrm{O}}$), ambipolar
($\eta_{\mathrm{A}}$), and Hall ($\eta_{\mathrm{H}}$), are calculated for each
simulated zone with the thermochemical conditions also elaborated in §2.2.

###  2.2 Non-equilibrium Thermochemistry with Consistency

The thermochemical processes significantly modulate the dynamic properties in
turbulent molecular zones, which are also important in observational
measurements and diagnostics. Technically, a “flow reaction” paradigm is
required, in which the MHD is in charge of advecting chemical species and
manipulating internal energy around the simulated space. Thermochemical
reaction networks are co-evolved with the MHD in every simulation zone.

The advection of chemicals with MHD is accomplished by specific methods so
that consistency of chemical elements and species fluxes in every timestep of
MHD evolution is ensured. The advected chemical abundances and internal energy
are then lead to the thermochemical solver for the current step, solving the
following coupled set of ordinary differential equations (ODEs) (in which the
Einstein summation convention applies):

$\dfrac{\mathrm{d}n^{i}}{\mathrm{d}t}=\mathcal{A}^{i}_{\;jk}n^{j}n^{k}+\mathcal{B}^{i}_{\;j}n^{j}\
;\quad\dfrac{\mathrm{d}\epsilon}{\mathrm{d}t}=\Gamma-\Lambda\ ;$ (4)

in which the terms $\\{\mathcal{A}^{i}_{\;jk}\\}$ describe two-body reactions,
those in $\\{\mathcal{B}^{i}_{\;j}\\}$ represent photoionization,
photodissociation, and spontaneous decays, and $\Gamma$ and $\Lambda$
represent all non-adiabatic heating, cooling, and heat transfer rates per unit
volume, respectively. We use a multi-step implicit method with error control
to solve since the ODEs (4) are usually stiff. The costs of co-evolving a
stiff network in each zone throughout the domain at every timestep could be
excessively high. For instance, a 3D simulation domain with $128^{3}$ zones
has to co-evolve $\sim 2\times 10^{6}$ stiff systems in each step. In order to
finish the simulations within reasonable “wall-clock time”, we use the
thermochemistry solver implemented on graphic processing units (GPUs) to
accelerate the simulations considerably for reasonable computing costs (see
also Wang & Goodman 2017; Wang et al. 2019).

The set of chemical species, reactions, and thermodynamic processes largely
follows Wang et al. (2019), with updates specific to diffuse molecular zones.
The 23 species are: $\mathrm{e}^{-}$ (free electrons), $\mathrm{H^{+}}$, H,
$\mathrm{H_{2}}$, $\mathrm{H_{2}^{+}}$, He, $\mathrm{He^{+}}$, O,
$\mathrm{O^{+}}$, OH, $\mathrm{H_{2}O}$, $\mathrm{H_{2}O^{+}}$,
$\mathrm{H_{3}O^{+}}$, C, $\mathrm{C^{+}}$, CO, $\mathrm{HCO^{+}}$,
$\mathrm{CH}$, $\mathrm{CH_{2}^{+}}$, Gr, $\mathrm{Gr^{+}}$,
$\mathrm{Gr^{-}}$. Here Gr and Gr± denote neutral and singly-charged dust
grains, respectively. This set covers the most concerned chemical species in
astronomical measurements. Internal energy density is also added as the last
“species” to the set, and the temperature dependences of
$\\{\mathcal{A}^{i}_{\;jk}\\}$ and $\\{\mathcal{B}^{i}_{\;j}\\}$ coefficients
are converted to internal energy dependences (via proper heat capacities) for
consistent thermochemical co-evolution. These “species” interact with each
other through a network with 146 chemical reactions and thermodynamic
processes (cooling and heating), most of which are included in the `UMIST`
database (McElroy et al. 2013), and other necessary processes are also
involved. To name a few mechanisms that are potentially important to the
diffuse molecular cloud evolution, we have:

  * •

Photoionizations by the interstellar radiation field (ISRF) and chemical
species, including the photoionization of $\mathrm{C}$, and photoelectric
effects of dust grains.

  * •

The ISRF photodissociation of multiple molecular species, e.g., $\mathrm{OH}$,
$\mathrm{H_{2}O}$, $\mathrm{CO}$, and $\mathrm{H_{2}}$.

  * •

Cosmic ray (CR) ionization of multiple neutral species, especially
$\mathrm{H_{2}}$, $\mathrm{H}$, $\mathrm{He}$, and $\mathrm{O}$.

  * •

Formation of important tracer molecules such as $\mathrm{CO}$ and
$\mathrm{OH}$, as subsequent products of ISRF photoionized $\mathrm{C^{+}}$,
and CR ionized $\mathrm{O^{+}}$ (see also Draine 2011).

  * •

Heating processes, including the energy injected via ISRF photoionization and
dissociation, CR ionization, and the $\mathrm{H_{2}}$ formation on the
surfaces of dust grains ($\sim 1~{}\mathrm{eV}$ per molecule formed).

  * •

Radiative cooling processes, including ro-vibrational transitions of
$\mathrm{CO}$, $\mathrm{OH}$, and $\mathrm{H_{2}}$ (Neufeld & Kaufman, 1993;
Omukai et al., 2010), and the fine structure transition of $\mathrm{C^{+}}$.

Note that the self-shielding effects for relevant photoionization and
photodissociation processes are included with the recipes in Heays et al.
(2017).

Both the radiation-related reactions and radiative cooling processes involve
the penetration of photons, and column densities of different chemical species
for photons to impinge and escape are the most important physical parameter.
For the local structure simulations using periodic boundary conditions, the
penetration columns for photons cannot be obtained consistently. Therefore, we
estimate the column densities of species X involved in shielding is estimated
by $N(\mathrm{X})\sim n(\mathrm{X})\bar{l}_{v}$, and the effective length
$l_{v}$ is approximated by

$\bar{l}_{v}\sim\mathrm{min}\left\\{l_{v},\ l_{\rm cap}\right\\}\ ;\quad
l_{v}\equiv\dfrac{|\mathbf{v}|}{|\nabla\cdot\mathbf{v}|}\ ,$ (5)

where $l_{\rm cap}$ approximates the distance to the nearest source of
radiation, describing case that the gas is largely uniform and static. There
are also alternatives to this Sobolev extimation of effective length, e.g.,
$l_{\rho}\equiv\rho/|\nabla\rho|$ and
$\bar{l}_{\rho}\sim\mathrm{min}\\{l_{\rho},\ l_{\rm cap}\\}$. However, the
estimates based on velocity best emulates the physical mechanisms involved for
the penetration of cooling and ionizing photons. For the diffuse region in the
Taurus molecular cloud (TMC hereafter), we assume that $l_{\rm cap}\sim
1~{}\mathrm{pc}$ represents the typical distance to one of the related
radiation sources (e.g. young stars) in the nearby regions. In addition, the
typical sizes of clumps inside the TMC diffuse regions are also considerably
smaller than $\sim 1~{}\mathrm{pc}$ (see also Goldsmith et al., 2010), making
$l_{\rm cap}\sim 1~{}\mathrm{pc}$ a resonably capping for the photon
penetration length. We will revisit and confirm this assumption when analyzing
the simulation results.

Similar to Wang et al. (2019), non-ideal MHD diffusivities are determined by
the gas temperature and the charge carriers’ abundances computed via the real-
time non-equilibrium thermochemical network. The general expressions for the
three diffusivities read (e.g. Bai, 2011),

$\begin{split}&\eta_{\mathrm{O}}=\dfrac{c^{2}}{4\pi}\left(\dfrac{1}{\sigma_{\mathrm{O}}}\right)\
,\quad\eta_{\mathrm{H}}=\dfrac{c^{2}}{4\pi}\left(\dfrac{\sigma_{\mathrm{H}}}{\sigma_{\mathrm{H}}^{2}+\sigma_{\mathrm{P}}^{2}}\right)\,\quad\\\
&\eta_{\mathrm{A}}=\dfrac{c^{2}}{4\pi}\left(\dfrac{\sigma_{\mathrm{P}}}{\sigma_{\mathrm{H}}^{2}+\sigma_{\mathrm{P}}^{2}}\right)-\eta_{\mathrm{O}}\
,\end{split}$ (6)

where $\sigma_{\mathrm{O}}$, $\sigma_{\mathrm{H}}$ and $\sigma_{\mathrm{P}}$
are Ohmic, Hall, and Pederson components of conductivity under magnetic
fields, and

$\begin{split}&\sigma_{\mathrm{O}}=\dfrac{ec}{B}\sum_{j}n_{j}Z_{j}\beta_{j}\
,\quad\sigma_{\mathrm{H}}=\dfrac{ec}{B}\sum_{j}\dfrac{n_{j}Z_{j}}{1+\beta_{j}^{2}}\
,\\\
&\sigma_{\mathrm{P}}=\dfrac{ec}{B}\sum_{j}\dfrac{n_{j}Z_{j}\beta_{j}}{1+\beta_{j}^{2}},\end{split}$
(7)

in which $Z_{j}e$ is the charge and $n_{j}$ the number density of the $j^{\rm
th}$ charged species, and the Hall parameter $\beta_{j}$ is the ratio of the
gyrofrequency to the collision rate with neutrals,

$\beta_{j}=\dfrac{Z_{j}eB}{m_{j}c}\dfrac{1}{\gamma_{j}\rho}\
;\quad\gamma_{j}\equiv\dfrac{\langle\sigma v\rangle_{j}}{\langle
m\rangle_{n}+m_{j}}\ ,$ (8)

where $m_{j}$ is the charged species’ molecular mass, $\langle m\rangle_{n}$
is the mean molecular mass of the neutrals, and $\langle\sigma v\rangle_{j}$
is the rate of collisional momentum transfer between the $j$th species and the
neutrals.

###  2.3 Turbulence Driving

From larger scales, turbulence energy cascades down to the simulated systems–a
localized, periodic, cubic domain–in this work. Various methods have been
adopted to emulate the injection of turbulence energy, and we adopt the mature
methods adopted by Mac Low (1999). The perturbations are inserted as a
globally uniform acceleration or inertial force. At each time step, a unity
vector with random direction $\hat{a}$ is assigned throughout the domain for
the direction of perturbations. Assuming the amplitude of perturbation is $A$,
the turbulence energy injection rate per by unit mass $\dot{\epsilon}_{\rm
turb}$ can be evaluated as,

$\dot{\epsilon}_{\rm turb}=A\left\\{[\langle\rho\rangle L_{\rm
box}^{3}]^{-1}\Delta t\ \sum_{\mathbf{i}}\delta
V_{\mathbf{i}}\;\rho_{\mathbf{i}}\mathbf{v}_{\mathbf{i}}\cdot\hat{a}\right\\}\
,$ (9)

where $L_{\rm box}$ is the box length, $\langle\rho\rangle$ is the mean mass
density throughout the domain, $\Delta t$ is the current timestep,
$\mathbf{i}$ is the cell index, $\delta V_{\mathbf{i}}$, $\rho_{\mathbf{i}}$
and $\mathbf{v}_{\mathbf{i}}$ are the cell volume, mass density, and velocity
vector for the $\mathbf{i}$-th cell, and the summation goes through all cells.
In practice, we specify a fixed $\dot{\epsilon}_{\rm turb}$ as one of the
model parameters throughout the evolution period. In each time step, we use
the summation in the braces of eq. (9) to determine the normalization factor
$A$, calculate the perturbation acceleration vector $A\hat{a}$, inject kinetic
energy into the system by
$\mathbf{v}_{\mathbf{i}}^{\prime}=\mathbf{v}_{\mathbf{i}}+A\hat{a}\Delta t$,
and adjust the total energy density accordingly.

###  2.4 Simulation Setups

Table 1: Properties of the Fiducial Model (§2.4)

Item | Value
---|---
Geometries |
Simulation domain | $(0.04~{}{\rm pc})^{3}$
Resolution | $N^{3}=128^{3}$
Irradiation |
Interstellar radiation field (ISRF) † | $0.3\;G_{0}$
Cosmic ray ionization ∗ | $\zeta=10^{-17}~{}\mathrm{s}^{-1}$
Initial MHD Parameters |
Temperature | $30~{}\mathrm{K}$
Density | $250~{}m_{p}~{}\mathrm{cm}^{-3}$
Magnetic Field | $3~{}\mu{\rm G}$
Initial abundances [$n_{\mathrm{X}}/n_{\mathrm{H}}$] |
$\mathrm{H_{2}}$ | 0.5
$\mathrm{He}$ | 0.1
$\mathrm{H}$ | $10^{-3}$
$e^{-}$ | $10^{-6}$
$\mathrm{H^{+}}$ | $10^{-6}$
$\mathrm{O}$ | $3.2\times 10^{-4}$
$\mathrm{C}$ | $1.4\times 10^{-4}$
Dust/PAH properties |
$a_{\mathrm{Gr}}$ |  $5$ Å
$\sigma_{\mathrm{Gr}}/\mathrm{H}$ | $\sim 10^{-23}~{}\mathrm{cm}^{2}/\mathrm{H}$

Note. — $\dagger$: $G_{0}$ denotes the Habing radiation flux (see also Draine
2011).

$*$: Total ionization rate per hydrogen atom by cosmic ray; the branching
ratio of ionization is proportional to the elemental abundances (see e.g. Bai
& Goodman (2009)).

In brief, we carry out 3D MHD non-equilibrium thermochemical simulations in a
(0.04 pc)3 box with periodic boundary conditions. Non-equilibrium
thermochemistry, whose network has 23 chemical species (plus internal energy
density) and 146 reactions, is co-evolved in every zone of the MHD simulation.
Non-ideal MHD effects with all components of tensorial conductivities are
included.

For the diffuse region outside the boundary in the TMC, the ISRF UV
photodissociation and photoionization parameter is $0.3~{}G_{0}$ (Pineda et
al., 2010). At $t=0$, we set the initial uniform mass density to be
$\rho_{t=0}=250~{}m_{p}~{}{\rm cm}^{-3}$ with $m_{p}$ the proton mass, and the
initial temperature to be uniformly $T_{t=0}=10~{}{\rm K}$ throughout the
entire simulation domain. We also put a $\mathbf{B}_{t=0}=B_{0}\mathbf{e}_{x}$
initial magnetic field, in which $B_{0}=3~{}\mu{\rm G}$ for this fiducial
model. We let the system evolve without additional energy input for
$10^{6}~{}{\rm yr}$ and examine the decay of the initial velocity
perturbations and the gas cooling to select a time that best represents the
thermochemical state of the TMC. These properties define the fiducial model,
summarized in Table 1. Using these properties on hydrostatic models, we obtain
two test cases that verify the applicability of our simulation system on
molecular clouds. These tests, which also serve as the reference point of our
subsequent analyses of thermochemistry for models with turbulences, are
elaborated in Appendix A.

##  3 Fiducial Model Results

Figure 1: Time evolution of key physical parameters of the fiducial model.
Upper panel (a): Domain-averaged specific energies (energy per unit mass;
$\epsilon_{k}$ for the kinetic energy, $\epsilon_{g}$ for the gas internal
energy, $\epsilon_{B}$ for magnetic energy, and $\epsilon_{\rm tot}$ for the
sum of all three components). Middle panel (b): Relative abundances of a few
key chemical species throughout the simulation domain (note that the
abundances of $\mathrm{OH}$ and $\mathrm{H_{2}O}$ are multiplied by $10^{4}$).
Lower panel (c): Average energy budgets per unit mass (solid lines for
cooling, dash and dotted lines for heating or injections).

The evolution history of the fiducial model (Model 0) is illustrated in Figure
1. The dynamical evolution, characterized by different components (internal
energy $\epsilon_{g}$, kinetic energy $\epsilon_{k}$, and magnetic energy
$\epsilon_{B}$) of the specific energy, has mostly reached its quasi-steady
state since $t\sim 10^{5}~{}\mathrm{yr}$. After that, the energy balance only
increases at a shallow rate ($\sim 10\%$ over $\sim 10^{6}~{}\mathrm{yr}$).
The thermochemistry evolution, nonetheless, has longer timescales. As the most
important coolants for material in the diffuse molecular phase, variations in
the abundances of a few key chemical species (e.g., $\mathrm{CO}$,
$\mathrm{C}$, $\mathrm{OH}$) fundamentally modulate the cooling channels in
the thermodynamic budgets of the simulated gas. Overall, the comparison in
timescales emphasizes the necessity of using non-equilibrium thermochemistry.
In Figure 2 we confirm that the effective penetration length estimated by
velocity divergences and density gradients yield results that are
qualitatively consistent, which fall well below the capping length $l_{\rm
cap}$. The visual extinction is mostly negligible ($A_{V}<0.5$ or even
$\lesssim 0.1$) for $\gtrsim 99.5\%$ of the simlated zones, which is
consistent with the studies for the diffuse regions TMC (e.g. Goldsmith et
al., 2010). In addition, the Jeans length is significantly greater than the
box scales even for compressed regions in the fiducial model. In addition, the
virial parameter

$\alpha_{\rm vir}\sim\dfrac{v^{2}+c_{s}^{2}}{G\rho l^{2}}\ ,$ (10)

are also large within the simulation domain. The median of $\alpha_{\rm vir}$
reads $\langle\alpha_{\rm vir}\rangle_{\rm median}\sim 3$ with $l\sim l_{J}$
over the Jeans length, and $\langle\alpha_{\rm vir}\rangle_{\rm median}\sim
10^{6}$ over the typical clump scales estimated with $l\sim l_{\rho}$. These
values indicate that the simulated diffuse regions are gravitationally
unbound, and neglecting the self-gravity is eligible for the analyses on
dynamics in this paper.

Figure 2: Upper panel: Distribution of volume fractions about the logarithm of
scale lengths ($\mathrm{d}f/\mathrm{d}\log_{10}l$), including extinction scale
lengths estimated by velocity divergence ($l_{v}$), density gradient
($l_{\rho}$), and the Jeans length ($l_{J}$). Lower panel: Distribution of
volume fractions about the logarithm of visual extinction, estimated by
$A_{V}\sim N_{\mathrm{H}}/(1.8\times 10^{21}~{}\mathrm{cm}^{-2})$ using
$l_{v}$ and $l_{\rho}$ for the effective column densities, respectively.

###  3.1 Fluid and Field Properties

Figure 3: Turbulence kinetic energy power spectra for fiducial model. Upper
row: Evolution of statistics over time; Lower row: Time-averaged statistics in
different evolution epochs; shaded bands illustrate the ranges of $2\sigma$
variation. The left column shows the power spectra (whose lower panel compares
them to the standard Kolmogorov $k^{-5/3}$ power law, and the $k^{-2}$ power
law), the middle column shows the distributions of spatial volume fraction
$f_{V}$ in the $\lg T$ space, and the right column shows the distributions in
the $\lg\rho$ space. Dotted curves on the lower-middle and lower-right panels
present the typical shapes of log-normal distributions for references.  Figure
4: Slices of physical parameters at the $z=0$ plane, illustrating the
$t=10^{6}~{}\mathrm{yr}$ snapshot for the fiducial model. The upper row shows
the key hydrodynamical parameters (mass density $\rho$, temperature $T$, and
${\rm Ma}_{\mathrm{A}}-1$ (${\rm
Ma}_{\mathrm{A}}\equiv|v|/[B/(4\pi\rho)^{1/2}]$ is the Alfvénic Mach number),
and the lower row shows relevant field parameters (magnetic field strength
$|B|$, plasma $\beta$ parameter, and magnetic Reynolds number ${\rm Re}_{\rm
p}$).

Figure 5: Slices at the $z=0$ plane for the fiducial for different evolution
time ($t=5\times 10^{4}~{}\mathrm{yr}$ for the upper row, and $5\times
10^{5}~{}\mathrm{yr}$ for the lower row), showing the relative abundances of
$\mathrm{C}$ (left column), $\mathrm{CO}$ (middle-left column), $\mathrm{CH}$
(middle-right column), and $\mathrm{OH}$ (right column), normalized to the
total density of hydrogen nuclei $n_{\mathrm{H}}$.

As the initial fields was set uniformly $B_{x0}=3~{}\mu{\rm G}$, the
corresponding initial plasma $\beta$ was $\beta_{0}=1.23$. Turbulent evolution
injects kinetic energy directly into the box, which subsequently branches into
thermal energy by shocks, viscous dissipations, and magnetic energy via flux
freezing. Similar to Figure 1, the evolution of kinetic energy spectrum and
distribution functions (for $\lg T$ and $\lg\rho$) illustrated in Figure 3
reaches quasi-steady states after $t\gtrsim 10^{5}~{}\mathrm{yr}$.

The kinetic energy spectrum does not show significant variation since
$t\gtrsim 10^{5}~{}\mathrm{yr}$, confirming that the quasi-steady state has
already been reached. Magnetization for the fiducial model is intermediate
($\beta\sim 1$), with most of the gas dynamics being sub-Alfvénic (Figure 4)
but supersonic (not shown in the figure). We observe a power law function
damped at $k\sim 1.3\times 10^{3}~{}\mathrm{pc}^{-1}$ (or equivalently $\sim
0.005~{}\mathrm{pc}$ spatial scale). At lower wavenumbers down to $k\sim
0.4\times 10^{3}~{}\mathrm{pc}^{-1}$ (lager scales at $\sim 1.5\times
10^{-2}~{}\mathrm{pc}$) the spetrum fits the Kolmogorov $E(k)\propto k^{-5/3}$
law. On a wider range of wavenumbers spanning roughly one decade (up to a
damping scale $k\sim 1.5~{}\times 10^{3}~{}\mathrm{pc}^{-1}$), the energy
power spectrum is comparable to the $E(k)\propto k^{-2}$ Burgers power law
(for weakly magnetized compressible plasma with shocks, e.g., Kowal et al.
2007; Boldyrev & Perez 2009). The fitting over the
$k/(10^{3}~{}\mathrm{pc})\in[0.5,1.5]$ range of wavenumbers prefers the
$k^{-2}$ power index ($\chi^{2}/N\sim 0.7$), and varying the wavenumber range
included yields power index fittings roughly ranging from $-1.9$ through
$-2.1$. Compared to the fluid morphologies in the configuration space (Figure
4), the damping scale is close to the sizes of a typical eddy ($\sim
0.003~{}\mathrm{pc}$). These sizes are significantly greater than the
numerical dissipation scale (roughly a few times the cell size, $\Delta
x\simeq 3\times 10^{-4}~{}\mathrm{pc}$). Therefore, the dissipation should be
physical rather than numerical or artificial. §4.1 re-examines such
dissipation scales from the aspect of fields.

The distribution functions of $\lg T$ and $\lg\rho$ exhibit more significant
variations than the power spectrum. The upper branch of the $\lg T$
distribution and the lower branch of the $\lg\rho$ distribution mostly fit the
log-normal distributions commonly observed in turbulences. In contrast, the
other sides (low temperatures and high densities) consist of the peaks above
the log-normal components that result directly from cooling processes. We
notice $\gtrsim 10\%$ of the simulated domain (volume fraction) whose
temperature is warm, $T>10^{2}~{}\mathrm{K}$. From the slices of the fiducial
model in Figure 4, we notice that these relatively warm regions usually have
lower (but not necessarily lowest) mass densities. Such temperatures will
likely populate the ro-vibrational energy levels of a few vital molecular
species. They may lead to important observational consequences, which will be
discussed in our following paper.

Compared to the initial conditions of $B_{0}=3~{}\mu\mathrm{G}$ and
$\beta_{0}=1.23$, the evolved magnetic fields depicted in Figure 4 demonstrate
a moderate enhancement due to turbulent dynamics, reaching magnitudes of
approximately $\sim 5-10~{}\mu\mathrm{G}$ predominantly in regions of
compression characterized by relatively high density ($\rho\gtrsim
10^{3}~{}m_{p}~{}\mathrm{cm}^{-3}$) and low temperature ($T\lesssim
20~{}\mathrm{K}$). Areas where the plasma beta ($\beta$) is significantly less
than unity typically exhibit magnetic field strengths close to the initial
value of $|B|\sim B_{0}$. These regions are predominantly composed of
kinematically decelerated gas with densities lower than the initial density,
($\rho\lesssim 10^{-1}\rho_{0}$), and temperatures well above average ($T\sim
10^{2}~{}\mathrm{K}$). The magnetization within these regions achieves a
quasi-steady state, representing a balance between the amplification effects
of turbulent motion and the counteracting influence of magnetic diffusion. The
lower-right panel presents ${\rm Re}_{\rm p}$, the magnetic Reynolds number
for the ambipolar (“Pederson”) component of magnetic diffusivities, defined
as,

${\rm Re}_{\rm p}\equiv\dfrac{|\mathbf{v}|l_{\rm dyn}}{\eta_{\rm
p}}\sim\dfrac{|\mathbf{v}|^{2}}{\eta_{\rm p}|\nabla\cdot\mathbf{v}|}\ ,$ (11)

in which we have used the uncapped
$l_{v}\equiv|\mathbf{v}|/|\nabla\cdot\mathbf{v}|$ to estimate the dynamical
length scale $l_{\rm dyn}$. The intensity of magnetic diffusion relative to
the inertia becomes more important with smaller ${\rm Re}_{\rm p}$. From the
figure, we observe a significant fraction of space where ${\rm Re}_{\rm
p}\lesssim 10$ or even less than unity, indicating that the ambipolar
diffusion becomes a major channel of dissipating the growths of fields through
turbulences. These ${\rm Re}_{\rm p}$ values emphasize the necessity of
including non-ideal MHD effects, which can only be accurately modeled by the
inclusion of real-time thermochemistry (see §2.1, §2.2).

###  3.2 Thermochemistry and Its Evolution

Profiles of chemical species vary in the simulation, from their abundances to
their forming and destruction pathways. Molecular hydrogen $\mathrm{H_{2}}$
are often traced by other molecules with relatively strong emission powers,
e.g., $\mathrm{CO}$, $\mathrm{OH}$, and $\mathrm{C^{+}}$. In Figure 5, we
present the abundances relative to hydrogen nuclei $n_{\mathrm{H}}$ of a few
key chemical species in the $z=0$ plane for the fiducial model at $t=5\times
10^{4}~{}\mathrm{yr}$, $2\times 10^{5}~{}\mathrm{yr}$, and
$10^{6}~{}\mathrm{yr}$, respectively.

####  3.2.1 Abundances of Molecules

The relative abundance of $X_{\mathrm{CO}}$ for CO, which is the most
frequently adopted $\mathrm{H_{2}}$ tracer, varies over time along the
evolution track. Formation of $\mathrm{CO}$ is lead by the chain from
$\mathrm{C^{+}}$ to $\mathrm{CH}$ (e.g. Draine, 2011), reaching its maximum
around $t\sim 2\times 10^{5}~{}\mathrm{yr}$ and starting to decline afterward
(which can also be observed in Figure 5). Such decline comes from the
increased abundance of atomic $\mathrm{H}$ competes with $\mathrm{O}$ in
consuming $\mathrm{CH}$, which is qualitatively similar to the tests in
Appendix A (see also the evolving $\mathrm{CH}$ abundances in Figure 5). The
timescale of overall $X_{\mathrm{CO}}$ variation is roughly
$\tau_{X_{\mathrm{CO}}}\gtrsim 2\times 10^{5}~{}\mathrm{yr}$. It is a
reasonable assumption that diffuse molecular gas coexists with the overall
molecular region, which may subsequently disperse due to intense
photodissociation following the formation of nearby young stellar objects. The
timescale over which $X_{\mathrm{CO}}$ varies is on par with the lifespan of
molecular regions in their entirety; therefore, this timescale merits close
consideration when one intends to investigate the long-term chemical evolution
of diffuse molecular gases.

We notice that the $\mathrm{CO}$ abundance is one order of magnitude above the
hydrostatic tests (Appendix A and Figure 13). Because the formation of
$\mathrm{CH_{2}^{+}}$, the rate-determining step in the reaction chain leading
to $\mathrm{CO}$, is primarily the two-body reaction
$\mathrm{C^{+}}+\mathrm{H_{2}}\rightarrow\mathrm{CH_{2}^{+}}$, the compression
caused by shocks and subsonic motions in turbulences largely enhances the
whole rate of formation. This enhancement also leads to the spatial variations
in $X_{\mathrm{CO}}$, as well as the stochasticity of $\langle
X_{\mathrm{CO}}\rangle$ even though the spatial average is taken across the
domain. In Figure 6, one can observe that the most of $\mathrm{CO}$ molecules
exist on the average temperature ($\sim 30~{}\mathrm{K}$) and relatively high
density, exhibiting distribution contours similar to $\mathrm{H_{2}}$ on the
$\lg n-\lg T$ plane. When coming to the $\lg n-\lg n(\mathrm{H_{2}})$ plane,
$X_{\mathrm{CO}}$ varies by about one decade at each $n_{\mathrm{H_{2}}}$, and
the median of $n_{\mathrm{CO}}$ increases faster than linear with greater
$n_{\mathrm{H_{2}}}$, which is another consequence of the two-body reaction
nature of the rate-determining step. Therefore, $X_{\mathrm{CO}}$ values
should be carefully adopted with their temporal and spatial variations taken
into account adequately when one wants to characterize diffuse molecular gases
with reasonable accuracy. Similar conclusions about $X_{\mathrm{CO}}$ have
emerged from various other simulation studies (e.g. Clark et al., 2019; Lupi
et al., 2021) . These studies emphasize the integration of co-evolved
chemistry with magnetohydrodynamics (MHD) when analyzing chemical tracers
within the molecular ISM. There are also extra complications at temperatures
below the threshold of $T\lesssim 20~{}\mathrm{K}$, CO molecules are prone to
further condensation onto the surfaces of dust grains (Cazaux et al., 2017).
Given that our study predominantly concentrates on the gas-phase
characteristics of diffuse molecular regions, we have opted to exclude these
additional complexities in the present paper. Instead, we intend to reserve
more exhaustive explorations for future research endeavors that delve into
multi-phase mechanisms, thereby providing a more comprehensive understanding
of the interplay between gas and dust in these environments.

The $\mathrm{OH}$ molecules show contrasts to $\mathrm{CO}$. Because of
turbulence compressions, they are $\sim 10^{4}$ more abundant than the
hydrostatic tests in Appendix A. Nonetheless, because $\mathrm{e^{-}}$ and
atomic $\mathrm{O}$ compete for the crucial progenitors ($\mathrm{H_{2}^{+}}$
and $\mathrm{H_{3}^{+}}$), regions with high densities could not always gain
more $\mathrm{OH}$ from increased two-body reaction rates. At higher
temperatures, the reaction between $\mathrm{e^{-}}$ and $\mathrm{H_{x}^{+}}$
($x=2,3$) is slightly suppressed (rate coefficient $\propto T^{-0.5}$; see
McElroy et al. 2013) while the channel with $\mathrm{O}$ leading to
$\mathrm{H_{2}O^{+}}$ and then $\mathrm{OH}$ is not affected. In addition, the
charge exchange channel
($\mathrm{H^{+}}+\mathrm{O}\rightarrow\mathrm{O^{+}}+\mathrm{H}$) that has
secondary importance in the formation of $\mathrm{OH}$ prefers higher
temperatures and higher $\mathrm{H^{+}}$ abundances, which mostly occur on the
edges of high-density regions where gas compression is likely to happen.
Therefore, the distribution of $\mathrm{OH}$ on the $\lg n$-$\lg T$ plane
prefers the compressed regions with intermediately low temperature and
densities, showing a “ridge”, roughly starts from the lower-left part (low
$n$, low $T$) to the upper-right corner (Figure 6). On the $\lg n-\lg
n(\mathrm{H_{2}})$ plane, the correlation between $\mathrm{OH}$ and
$\mathrm{H_{2}}$ is loose and should not be used alone to deduce the
properties of $\mathrm{H_{2}}$.

####  3.2.2 Energy Dissipation

Dissipation of energy is crucial to the quasi-steady turbulent structures in
the system. Short of cooling mechanisms, the fluids will eventually stop
behaving turbulently after their total internal energy ($\epsilon_{g}$)
exceeds the amount of energy injection from the largest scales. Enumerating
all prospectively important cooling processes in the simulation, we identify a
few major cooling mechanisms in the fiducial model and present their
corresponding heating and cooling timescale maps in the top and middle rows of
Figure 7.

It is observed that cooling via the fine structure transitions of
$\mathrm{C^{+}}$ mainly appears at boundaries separating dense and diffuse
regions, while the $\mathrm{CO}$ ro-vibrational transitions are responsible
for cooling the dense gas. The rotational transitions of $\mathrm{CO}$ are
characterized by lower energy separations. For instance, the transition
$J=1\rightarrow 0$ corresponds to a wavelength of $\lambda=2.6~{}{\rm mm}$.
This is in contrast to the fine structure transition of ionized carbon
($\mathrm{C^{+}}$), particularly the [C ii] line at $\lambda=158\mathrm{\mu
m}$ via which most of the $\mathrm{C^{+}}$ cooling takes place. The
distinction between these two cooling channels is more apparent at relatively
low temperatures $T\lesssim 30~{}\mathrm{K}$. This natural divergence at
cooler regimes has been previously discussed by e.g. Richings et al. (2014a,
b), who utilized thermochemical calculations in the absence of hydrodynamics
or MHD. In the context of co-evolved non-ideal MHD, this phenomenon is further
intertwined with the turbulence dynamics. The compression of gas through
turbulence facilitates more efficient production of $\mathrm{CO}$ from
$\mathrm{C^{+}}$ through a sequence of chemical reactions. Concurrently, the
increase in collisional partners augments the difference in cooling efficiency
between $\mathrm{CO}$ and $\mathrm{C^{+}}$ at low temperatures.

Due to relatively low temperature, the cooling mechanism on the third
panel–the ro-vibrational transitions of $\mathrm{CO}$–are $\sim 10^{-2}$ less
efficient than $\mathrm{C^{+}}$ and $\mathrm{H_{2}}$. The detectability of
these transitions could lead to significant observational effects, which will
be discussed in more detail in our following paper. In total, the evolution of
multiple cooling mechanisms is compared to the energy injection mechanisms:
(1) turbulent energy injection, and (2) the formation of $\mathrm{H_{2}}$ on
grain surfaces (each $\mathrm{H_{2}}$ formed will introduce $\sim
1~{}\mathrm{eV}$ into the gas thermal energy; e.g. Pantaleone et al. 2021).
The overall energy balance is close to equilibrium, yet with a small surplus
on the heating side (around $\sim
10^{-4}~{}\mathrm{erg}~{}\mathrm{g}~{}\mathrm{s}^{-1}$). Multiplying this
surplus by $\sim 10^{6}~{}\mathrm{yr}$ explains the gradual increase in the
total specific energy (Figure 1 a).

Figure 6: Distribution of spatial volume fraction $f_{V}$ in the spaces
spanned by $\\{\lg T\\}\times\\{\lg n({\rm X})\\}$ (upper row) and $\\{\lg
n(\mathrm{H_{2}})\\}\times\\{\lg n({\rm X})\\}$ (lower row), where X stands
for different chemical species indicated on each panel.  Figure 7: Important
thermochemical rates, quantified as the reciprocals of timescales
($\tau^{-1}$), illustrating the $z=0$ plane for the fiducial model at the
$t=2\times 10^{6}~{}\mathrm{yr}$ snapshot. The top and middle rows show the
heating and cooling processes. The upper-left panel illustrates the heating
via the ISRF dissociation and ionization, the upper-middel panel presents the
heating by cosmic rays, and the rest panels are for cooling: the fine
structures transitions of $\mathrm{O}$ (upper-right panel) and
$\mathrm{C^{+}}$ (middle-left panel), and the ro-vibrational transitions of
$\mathrm{CO}$ (middle panel) and $\mathrm{H_{2}}$ (middle-right panel). The
cooling rates by $\mathrm{OH}$ are relatively tiny and thus omitted. The
bottom row shows the rates related to $\mathrm{H_{2}}$ number densities: the
destruction rates via the reactions with $\mathrm{C^{+}}$ (left panel) and
photodissociation (middle panel), and the formation rates via the catalyzed
reaction by dust grains (right panel).

##  4 Parameter Studies: Various Models

Table 2: Physical parameters of various models (§4)

|  | $\langle\epsilon\rangle/(10^{9}~{}\mathrm{erg}~{}\mathrm{g}^{-1})$ | $\log_{10}[n({\rm X})/n_{\mathrm{H}}]$ | $\log_{10}[\langle\dot{\epsilon}\rangle_{\rm cool}/(\mathrm{erg}~{}\mathrm{g}^{-1}~{}\mathrm{s}^{-1})]$
---|---|---|---|---
Model | Description | Thermal | Kinetic | Magnetic | $e^{-}$ | $\mathrm{C^{+}}$ | $\mathrm{CO}$ | $\mathrm{OH}$ | $\mathrm{C^{+}}$ | $\mathrm{CO}$ | $\mathrm{H_{2}}$ | $\mathrm{O}$
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11) | (12) | (13)
0 | Fiducial model | $2.8$ | $3.0$ | $1.5$ | $-3.99$ | $-4.01$ | $-6.58$ | $-7.87$ | $-2.55$ | $-4.73$ | $-3.57$ | $-4.05$
I | Ideal MHD | $2.7$ | $3.2$ | $1.9$ | $-3.99$ | $-4.01$ | $-6.83$ | $-7.84$ | $-2.53$ | $-4.90$ | $-3.68$ | $-4.29$
B0 | No magnetic fields | $2.7$ | $3.4$ | $0.00$ | $-3.99$ | $-4.01$ | $-6.61$ | $-7.73$ | $-2.55$ | $-4.67$ | $-3.57$ | $-4.13$
B1 | $B_{x0}=1\ \mu{\rm G}$ | $2.7$ | $3.4$ | $0.5$ | $-4.00$ | $-4.01$ | $-6.33$ | $-7.84$ | $-2.57$ | $-4.46$ | $-3.52$ | $-4.04$
B10 | $B_{x0}=10\ \mu{\rm G}$ | $2.8$ | $1.7$ | $9.9$ | $-3.99$ | $-4.01$ | $-6.52$ | $-8.06$ | $-2.54$ | $-4.81$ | $-3.77$ | $-4.22$
T4 | $\dot{\epsilon}_{\rm turb}=3\times 10^{-4}{\rm\ erg\ s}^{-1}{\rm g}^{-1}$ | $1.9$ | $0.6$ | $1.1$ | $-3.98$ | $-4.01$ | $-7.05$ | $-9.35$ | $-3.35$ | $-5.53$ | $-4.91$ | $-6.43$
T2 | $\dot{\epsilon}_{\rm turb}=3\times 10^{-2}{\rm\ erg\ s}^{-1}{\rm g}^{-1}$ | $6.5$ | $13.1$ | $2.6$ | $-4.02$ | $-4.03$ | $-6.45$ | $-7.01$ | $-1.77$ | $-4.03$ | $-1.99$ | $-2.59$
N4 | $\rho_{0}=2.5\times 10^{4}\ m_{p}{\rm\ cm}^{-3}$ | $2.8$ | $3.3$ | $0.1$ | $-6.21$ | $-6.53$ | $-4.36$ | $-8.21$ | $-3.93$ | $-2.60$ | $-3.35$ | $-4.36$
N6 | $\rho_{0}=2.5\times 10^{6}\ m_{p}{\rm\ cm}^{-3}$ | $2.4$ | $3.2$ | $0.01$ | $-8.10$ | $-9.49$ | $-4.40$ | $-10.45$ | $-6.53$ | $-2.56$ | $-3.59$ | $-5.62$
G0 | $G=0$ | $3.9$ | $3.2$ | $0.9$ | $-5.36$ | $-5.82$ | $-4.56$ | $-6.32$ | $-3.29$ | $-2.97$ | $-3.02$ | $-3.33$
G3 | $G=3\ G_{0}$ | $3.0$ | $2.9$ | $1.5$ | $-3.98$ | $-4.00$ | $-8.79$ | $-8.69$ | $-2.45$ | $-6.67$ | $-3.52$ | $-3.91$
CR0 | $\xi_{\rm CR}=0$ | $2.7$ | $2.9$ | $1.4$ | $-4.02$ | $-4.02$ | $-5.94$ | $-10.04$ | $-2.57$ | $-4.32$ | $-3.59$ | $-4.05$
CR15 | $\xi_{\rm CR}=10^{-15}\ {\rm s}^{-1}$ | $3.8$ | $2.9$ | $1.7$ | $-3.53$ | $-4.00$ | $-7.48$ | $-6.74$ | $-2.11$ | $-5.32$ | $-3.27$ | $-3.73$
M01 |  $0.1\times$ Fiducial metal | $3.7$ | $2.9$ | $1.1$ | $-4.65$ | $-5.00$ | $-8.46$ | $-7.61$ | $-2.69$ | $-6.30$ | $-2.99$ | $-3.86$
M10 |  $10\times$ Fiducial metal | $2.8$ | $3.0$ | $1.6$ | $-3.44$ | $-3.44$ | $-5.93$ | $-8.06$ | $-2.55$ | $-4.09$ | $-3.62$ | $-4.27$

Note. — (1)–(2): Model codes and descriptions; each model is different from
the fiducial model by only one physical parameter.

(3)–(5) Average specific energy components (thermal
$\langle\epsilon_{g}\rangle$, kinetic $\langle\epsilon_{k}\rangle$, and
magnetic $\langle\epsilon_{b}\rangle$). Note that for a predominantly
molecular gas, the specific energy components can be converted to temperature
as $\langle T\rangle\simeq
11.4~{}\mathrm{K}~{}[\langle\epsilon_{g}\rangle/(10^{9}~{}\mathrm{erg}~{}\mathrm{g}^{-1})]$,
and root-mean-squared velocity $\langle v^{2}\rangle^{1/2}\simeq
0.45~{}\mathrm{km}~{}\mathrm{s}^{-1}[\langle\epsilon_{k}\rangle/(10^{9}~{}\mathrm{erg}~{}\mathrm{g}^{-1})]^{1/2}$.
The average plasma $\beta$ can also be estimated by
$\langle\beta\rangle\simeq\langle\epsilon_{g}\rangle/\langle\epsilon_{b}\rangle$.

(6)–(9): Logarithms of relative abundances of key chemical species.

(10)–(13): Specific cooling rates (energy loss rate per unit mass) via
different coolants.

As §3 discovered, various physical parameters are related to the behaviors in
the simulated diffuse cloud. In order to examine their prospective influences,
we conduct extra simulations examining the outcomes of modifying important
physical parameters. Table 2 lists relevant models marked by different
symbols; each varies only one parameter at a time.

###  4.1 Fields and Fluids

Figure 8: Statistics of the simulation series exploring the effect of non-
ideal MHD effects by comparing Models 0 (fiducial) and I (ideal MHD), showing
the kinetic energy spectra (top panel; an $E(k)\propto k^{-2}$ line is also
presented for reference), $\lg T$ distributions (middle panel), and $\lg\rho$
distributions (lower panel).  Figure 9: Similar to Figure 8 but for the series
with different initial magnetic field strengths (Models 0, B0, B1, and B10).

Turbulences can amplify the fields by flux freezing. This process is dispersed
and dissipated by the non-ideal MHD effects, and we confirm this point by
comparing Model 0 and Model I (ideal MHD). Figure 8 illustrates the difference
between these two models. Although the power spectra are identical at the
spatial scales of energy injection (viz. lowest wavenumbers), the inertial
range of Model I has greater kinetic power, $\sim 3$ times Model 0. The
damping wavenumber of the fiducial Model 0 is also lower ($k\simeq 1.3\times
10^{3}~{}\mathrm{pc}^{-1}$) compared to Model I ($k\simeq 2.2\times
10^{3}~{}\mathrm{pc}^{-1}$), exhibiting the effects of magnetic diffusion that
dissipates the kinetic energy at larger scales with greater efficiencies. The
distribution function in the $\lg T$ space for Model I shows that a
considerably smaller fraction of gas has temperatures above $\sim
10^{2}~{}\mathrm{K}$ compared to Model 0, which is attributed to the absence
of non-ideal MHD heating, as well as the tighter fluid-field coupling that
reduces the shocks and collisions. These differences confirm the importance of
including proper and consistent non-ideal MHD effects when simulating
magnetized molecular clouds.

Coupling between fluids and fields has another side effect: the energy
transfer to smaller scales is inhibited. Turbulent motions on smaller scales
tend to amplify the fields more efficiently, causing higher branching ratios
into magnetic energy. In the simulation series including Models B0 (a pure
hydrodynamic simulation without any fields), fiducial, B1, and B10, the power
spectra are very similar in the range $k\lesssim 1\times
10^{3}~{}\mathrm{pc}^{-1}$. At the same time, the amplitudes at higher
wavenumbers are anti-correlated to the initial field strengths. The situation
of Model B10 is more severe than others, whose strong magnetic fields suppress
the kinetic power by almost one decade and narrow the distribution functions
in the $\lg T$ and $\lg\rho$ spaces. This phenomenon can constrain the upper
limits of field intensities in observed molecular gases. For example, the
average field strength is unlikely to be greater than $\sim
10~{}\mu\mathrm{G}$ given the discovery of considerable fractions of $T\gtrsim
1.5\times 10^{2}~{}\mathrm{K}$ warm regions. With the rotational transition of
$\mathrm{H_{2}}$ observed identifying gas temperature $T\gtrsim
200~{}\mathrm{K}$ in the diffuse molecular gas in the TMC region (e.g.
Goldsmith et al., 2010), this result prospectively provide a constraint on the
strength of magnetization therein.

Recent simulations have delved into the intricacies of magnetized molecular
clouds, with notable studies such as that by Grete et al. (2021), which
explored the suppression of kinetic energy cascade due to magnetic tension.
This finding aligns well with the observed trends across a range of
simulations. However, it is crucial to acknowledge that the numerical
diffusivities associated with both fluids and fields, which vary with spatial
resolution in these simulations, may be interlinked with the identification of
damping effects attributed to magnetic tension. For a thorough and
quantitative examination of how small-scale structure growth is curtailed by
magnetic fields, it is essential to ensure a consistent coupling between
fields and fluids. This should be governed by appropriate magnetic
diffusivities and should include a consistent approach to the ultimate
dissipation of energy through thermodynamic processes. In the current work,
the suppression of energy cascade by field tensions is observed to influence
the power spectra at higher magnetic field strengths. The threshold for this
effect is distinctly identifiable between magnetic field strengths of
$B_{0}=3~{}\mu\mathrm{G}$ and $10~{}\mu\mathrm{G}$. Below this threshold, the
damping of fluid fluctuations predominantly occurs through non-ideal MHD
dissipation mechanisms. This can be inferred by comparing the damping
wavenumbers across Models B3, B10, 0 (Fiducial), and I (ideal MHD).

Higher average mass density yields large plasma $\beta$, lowering the relative
significance of magnetic pressure and stresses. As a result, Models N4 and N6
show very similar power spectra, which have considerably greater amplitudes
compared to the fiducial Model 0 (see Figure 10)–these spectra are also very
similar to the field-free Model B0. As the cooling rates are roughly
proportional to the square of mass density at a fixed temperature, while the
turbulence energy injection rate per unit volume is proportional to the mass
density itself (given fixed $\dot{\epsilon}_{\rm turb}$ per unit mass), the
temperature distributions of Models N4 and N6 appear on the lower end compared
to Model 0.

Another important process affected by the average density is the recombination
of ions. As the recombination rate coefficients are mostly the functions of
temperatures, assuming that the temperature is roughly fixed and that there is
only one ionizable species ($\mathrm{X}\leftrightarrow e^{-}+\mathrm{X^{+}}$)
thus $n_{e}=n(\mathrm{X^{+}})$, the ionization–recombination balance roughly
reads,

$\alpha x_{e}^{2}n_{\mathrm{H}}^{2}\sim\zeta n_{\mathrm{H}}\Rightarrow
x_{e}\sim n_{\mathrm{H}}^{-1/2},\ n_{e}\equiv x_{e}n_{\mathrm{H}}\sim
n_{\mathrm{H}}^{1/2}.$ (12)

In the series consisting of Models 0, N4, and N6, the fraction $x_{e}$ drops
faster than $\rho^{-1/2}$. The main charge carrier in Model N4 before $t\sim
2\times 10^{5}~{}\mathrm{yr}$ is still $\mathrm{C+}$, and changes to
$\mathrm{H^{+}}$ afterward due to efficient formation of $\mathrm{CO}$ that
uses up almost all carbon. Model N6 also has $\mathrm{H^{+}}$ as the main
charge carrier, but from a different origin: its self-shielding and cross-
shielding for the ISRF ionization of neutral $\mathrm{C}$ is significant that
the dominating ionization process is the cosmic ray ionization of
$\mathrm{H_{2}}$ in the first place. Without efficient $\mathrm{CO}$
formation, the specific energy keeps growing until $\sim
10^{5}~{}\mathrm{yr}$, after which $\langle X_{\mathrm{CO}}\rangle$ exceeds
$\sim 10^{-5}$, and the subsequent cooling causes a continuous drop in gas
temperature over the following $10^{6}~{}\mathrm{yr}$.

Figure 10: Similar to Figure 8 but for the series with different average mass
densities (Models 0, N4 and N6).

###  4.2 Turbulence Energy Injection

Turbulent energy injection (see §2.3) is the ultimate energy source in all
simulations throughout this paper. By turning the knob that controls the
injection rate $\dot{\epsilon}_{\rm turb}$ in Models T2 and T4 (injection
rates $3\times 10^{-2}$ and $3\times 10^{-4}$ in units
$\mathrm{erg}~{}\mathrm{g}^{-1}~{}\mathrm{s}^{-1}$, respectively), the
turbulence cascade and almost all phenomena attached to it will experience
considerable changes constrained by their energy budgets.

Figure 11: Similar to Figure 8 but for the series with different turbulence
energy injection rates (Models 0, T2, and T4). Note that the power spectra for
T2 and T4 have adjusted by scaling factors (indicated in the legend), so that
their spectra on the lowest wavenumbers overlap for better comparison.

The most prominent consequence of $\epsilon_{\rm turb}$ is the distribution of
gas dynamics and thermodynamics. Figure 11 compares Models T2 and T4 with the
fiducial Model 0 for their statistical functions. The shapes of the kinetic
energy spectra are almost identical at wavenumbers below the damping
wavenumber $k\simeq 1.5\times 10^{3}~{}\mathrm{pc}^{-1}$, while the amplitude
grows sublinearly with increasing intensity of turbulence energy cascades.
Above that wave number, the cases with stronger energy injection are
relatively less susceptible to the inhibition of the kinetic power by the
damping mechanisms including magnetic fields and diffusivities on smaller
spatial scales, and vice versa.

The density and temperature distributions in this series exhibit a clear trend
that higher $\dot{\epsilon}_{\rm turb}$ leads to “broader” distributions in
both $\lg T$ and $\lg\rho$ space by allowing for more vigorous shocks,
compressions, and adiabatic expansions. This fact suggests a promising way to
quantify the turbulence energy crossing rates via the measurements of
temperature and density distributions of concerned molecular gases. For
instance, the temperature distributions are reasonably sensitive to
$\dot{\epsilon}_{\rm turb}$, as Model T4 exhibits little gas warmer than $\sim
60~{}\mathrm{K}$, while Model T2 has a non-negligible fraction with $T\gtrsim
10^{3}~{}\mathrm{K}$. Admittedly, multiple other factors affect gas
temperature distributions, causing parametric degeneracies. Therefore,
quantitative measurement of the turbulence power relies on other constraints,
e.g., the abundance of $\mathrm{OH}$ is positively correlated to
$\dot{\epsilon}_{\rm turb}$ (Table2), due to the enhanced reaction rates at
higher density and temperatures. Synthetic mock observations are necessary in
the following works to confront the simulation results directly with
observations.

###  4.3 Thermochemical Conditions

Figure 12: Similar to the middle panel of Figure 8, with models varying the
cosmic ray intensities (left column; Models 0, CR0, and CR15), and the ISRF
(right column; Models 0, G0, and G3). The power spectra and $\lg\rho$
distributions are omitted as they are very similar to the fiducial model.

There are two major ionization channels in the simulation: (a) ISRF ionization
of carbon and (b) cosmic ray ionization (mainly on $\mathrm{H_{2}}$,
$\mathrm{H}$ and $\mathrm{He}$). These processes not only provide necessary
coupling between gas and fields but also effectively initiate the formation of
relevant chemical species, including $\mathrm{C^{+}}$, $\mathrm{CO}$,
$\mathrm{OH}$, which are the most important ones in terms of the gas
thermodynamics.

The statistics of simulations have revealed that the kinetic energy spectra
and the $\lg\rho$ distribution functions of these models are affected only
very slightly. Therefore, Figure 12 only emphasizes the comparisons in the
$\lg T$ space. Model G3 has a relatively similar $\lg T$ distribution function
compared to Model 0, with the lower end slightly warmer due to the increased
ISRF destruction of molecular coolants ($\mathrm{CO}$ and $\mathrm{OH}$).
Model G0, with the lack of production in the most critical coolants
$\mathrm{C^{+}}$ and $\mathrm{CO}$, exhibits a $\lg T$ profile that mostly
fits the standard log-normal distribution. With excessive cosmic ray
ionization, Model CR15 has a $\sim 2.5\times$ higher cutoff on the low-
temperature end due to extra cosmic ray heating, while its high end of $\lg T$
distribution is qualitatively similar to Model I due to higher abundances of
charge carriers and tighter fluid-field couplings. Since the most important
coolants are often metal-related, altering the gas’s metallicity impacts the
gas temperatures. Models 0, M01, and M10 have almost the same lower limit in
the temperature distribution functions (approximately $T\sim
10~{}\mathrm{K}$), while their high ends are anti-correlated with the
metallicity. In brief, the thermodynamic conditions can provide specific
constraints on the properties of the irradiation and cooling processes in
molecular gases.

##  5 Discussions and Summary

In this work, we have studied the characteristics and evolution of diffuse
($\rho\sim 250~{}m_{p}~{}\mathrm{cm}^{-3}$), magnetized ($B_{0}=3~{}\mu G$)
molecular clouds, using MHD simulations co-evolved with consistent non-
equilibrium thermochemistry. Although mostly numerical and theoretical, these
simulations have revealed prospective clues that can play crucial roles in
constraining diffuse molecular gas properties via observational measurements.

###  5.1 Fields, Fluids, and their Interactions

Co-evolved with thermochemistry, we observe from the simulation models that
the behaviors of gas dynamics are predominantly controlled by the turbulence
and MHD properties themselves. The $\lg T$ distributions (especially their
high-temperature cutoff) are hopefully the most sensitive indicator of the
turbulence energy injection rate. With consistent thermochemistry considered,
the thermal energy budgets at quasi-steady states are mostly balanced, and the
temperature profiles will reveal the dynamics relatively robustly. The $\lg T$
profiles are susceptible to other factors, including total densities, ionizing
radiations (and cosmic rays), field intensities, and metallicities.
Nevertheless, as discussed in §4, these factors only affect $\lg T$ with
secondary importance, especially regarding the high-temperature cutoffs in
these diffuse molecular gases. In observations aiming at constraining the gas
dynamic properties, various nebula diagnostics can determine the gas
temperature distribution in interested regions and constrain the turbulence
energy injections with relatively high confidence.

The field strengths we have studied are marginally crucial in curbing the
formation of turbulence structures, except for Model B10, whose average plasma
$\beta\lesssim 10^{-1}$ causes significant inhibition of turbulent structure
growths. The simulated models are mostly weakly ionized, and non-ideal MHD is
desired to model the kinematics and dynamics. Due to magnetic diffusions,
hydrodynamic profiles are not very sensitive probes in the fields. However,
the upper limit of fields can still be constrained by the temperature assisted
with other measurements.

###  5.2 The Necessity of Non-ideal MHD Studies with Consistent
Thermochemistry

####  5.2.1 Magnetic Fields

There have been many recent numerical simulations aiming at the turbulence
dynamics within the interstellar medium (ISM), including the behavior within
molecular clouds. A prime example is the work by Federrath et al. (2021),
which achieved remarkable resolution through hydrodynamic simulations, albeit
without incorporating magnetic fields. One of the major conclusions of this
study is the crucial role that magnetic fields play in sculpting the
turbulence power spectra, especially the damping wavenumber. The magnetic
pressure and tension are found to impede the development of structures at
smaller scales. In scenarios devoid of magnetic fields, the energy cascade is
observed to extend to the highest attainable wavenumbers, constrained only by
numerical dissipations at the resolution limits. However, this overlooks the
actual spatial scales at which damping occurs, given that magnetization is a
pervasive characteristic of the ISM under consideration. This oversight
underscores the necessity for simulations to account for the influence of
magnetic fields to more accurately reflect the ISM’s behavior.

####  5.2.2 Magnetic Diffusivities

Even the simulations incorporating MHD, such as Ostriker et al. (1999); Porter
et al. (2015); Bian & Aluie (2019); Grete et al. (2021), often overlook a
critical aspect: the influence of non-ideal MHD effects. As depicted in Figure
4 and corroborated by numerous other simulations, regions that are susceptible
to non-ideal MHD effects (indicated by ${\rm Re}_{\rm p}$ at the order of
unity) constitute a significant portion of the simulated volume, particularly
at the peripheries of compressed gas clumps. This phenomenon specifically
impacts the morphology at smaller scales, as illustrated by the comparison
between Model 0 and Model I in Figure 8. The damping scale, influenced by
magnetic diffusivities, is approximately twice as extensive as in scenarios
modeled with ideal MHD, where the effective magnetic diffusivities are
numerical artifacts contingent on resolution.

Also elaborated by Grete et al. (2023), the dissipation scales in MHD play a
pivotal role in determining how magnetic fields sculpt the dynamics of the
gas. The accurate incorporation of magnetic diffusivities in simulations that
resolve the spatial scales of diffusion is essential to circumvent such
artifacts. Moreover, magnetic diffusivities are sensitively dependent on the
actual thermochemical conditions of the gas, including the concentration and
types of charge carriers, the densities of neutral gas, and the prevailing
temperatures. These physical variables are contingent on the evolutionary
history of the gases and should not be reduced to mere status variables.
Consequently, non-ideal MHD simulations of molecular gases that employ
prescribed parametrization of diffusivities or highly simplified chemical
networks, as seen in Downes & O’Sullivan (2009); Tilley & Balsara (2011);
Jones & Downes (2011, 2012); Wurster et al. (2016), inevitably sacrifice
consistency. This is due to the fact that many reactions related to
ionization, recombination, and charge transfer involve complex molecular and
ionic species, and sometimes even surface chemical processes. This realization
has already been acknowledged in studies of protoplanetary disks (e.g., Bai &
Goodman, 2009; Xu & Bai, 2016; Wang et al., 2019) and has been identified as
being equally significant for the diffuse molecular gases discussed in this
paper, as exemplified in Figure 12.

####  5.2.3 Equation of State and Thermodynamics

The majority of prior research works, either with hydrodynamics or MHD, have
employed an isothermal equation of state (EoS), neglecting the essential
dissipation mechanisms carried out by comprehensive thermochemical approaches.
It is acknowledged that simulations relying solely on an adiabatic EoS in the
absence of cooling will eventually convert all injected energy into thermal
energy of the gas, overwhelming the kinetic and magnetic energy and finally
smearing out the turbulent structures. However, the use of an isothermal EoS
as a simplification also significantly compromises physical fidelity.

Consider the scenario of a gas clump that undergoes rapid compression prior to
cooling. The gas pressure, and consequently its internal energy density,
scales with the length scale $l$ of the clump as $p\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern
2.0pt\cr\sim\cr\kern-2.0pt\cr}}}l^{-3\gamma}$, while the density scales as
$\rho\mathrel{\vbox{ \offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern
2.0pt\cr\sim\cr\kern-2.0pt\cr}}}l^{-3}$. The cooling rate per unit volume at a
constant temperature is proportional to $\Lambda\rho^{2}\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern
2.0pt\cr\sim\cr\kern-2.0pt\cr}}}l^{-6}$, a steeper power-law than that of the
gas pressure. This indicates that the compressed gas patch will undergo
runaway cooling until it reaches a significantly lower temperature, at which
$\Lambda$ almost diminishes compared to its ambient value (for molecular gas,
this temperature is typically around $\sim 10~{}\mathrm{K}$ ), while other
regions of the gas remain at higher temperatures. This phenomenon results in a
broad temperature distribution that spans several orders of magnitude, as
observable in e.g., Figures 4 and Section 4.

The presence of a broad temperature distribution within turbulent materials is
essential for interpreting the diagnostics of observed nebulae. For instance,
the work by Goldsmith et al. (2010) has already identified the emission of
rotational transitions of $\mathrm{H_{2}}$ from diffuse molecular regions with
a median temperature $\lesssim 50~{}\mathrm{K}$. This finding underscores the
importance of regions heated to temperatures above $\gtrsim
10^{2}~{}\mathrm{K}$ through gas dynamics, thereby highlighting the need for
consistent thermodynamics in ISM simulations. In contrast, the simplified
dichotomy of “supersonic versus subsonic scales” as emphasized in isothermal
simulations, including those by Federrath et al. (2021), does not accurately
reflect the complexities of astrophysical phenomena taking place in the
diffuse molecular gas.

When aligning simulation results with ISM observations, the acquisition of
consistent temperature profiles is imperative for the precise interpretation
of spectral diagnostics pertaining to the properties of the gas. These
considerations underscore the necessity for the integration of robust
thermochemical computations within simulation models, moving beyond the
limitations of the overly simplified isothermal equation of state. This
approach is vital for capturing the nuanced interplay of physical processes
that govern the behavior of the interstellar medium.

####  5.2.4 Dimensionless versus Dimensional

Most of the simulations on turbulent ISM follow the dimensionless approach.
When the final dissipation is based on realistic microphysical mechanisms
(including magnetic diffusion, radiative cooling, etc.) rather than
mathematical or numerical dissipations, one can easily confirm with basic
dimensional analyses that the microphysical processes cannot be non-
dimensionalized simultaneously with the hydrodynamics or MHD. This fact
inevitably add to the complications and computational costs, yet this is the
inevitable costs that one has to pay to understand turbulent gases that are
coupled with physical mechanisms that are crucial to the ultimate dissipation
along the turbulence energy cascade, including magnetic fields and multi-scale
physics.

###  5.3 Thermochemical Probes of Diffuse Molecular Gases

Given relatively strong turbulence driving processes, constraining
thermochemistry via dynamic profiles are indirect and may suffer from
degeneracies of parameters—probably except for the influence by metallicities
on the higher limits of temperatures, which can be constrained in conjunction
with the specific turbulence powers via the upper limits in $T$. Instead,
thermochemistry itself is the probe of multiple dynamical processes. This
probe can only be forged with turbulent ISM simulations that can directly
yield analyses consistently with real-time thermochemistry. Adopting
isothermal EoS may mislead the interactions connecting the gas dynamics to
chemistry and subsequent observables, by providing inconsistent key parameters
including compression ratios and temperatures.

One prospectively interesting indicator that deserves extra effort is the
$\mathrm{CO}$ relative abundances, $X_{\mathrm{CO}}$. Presumably, this is
mostly a constant used to infer total gas densities in many studies, yet our
simulations suggest that $X_{\mathrm{CO}}$ varies from place to place, from
time to time. The conversion from $\mathrm{C^{+}}$ to $\mathrm{CO}$ takes
$\sim 2\times 10^{5}~{}\mathrm{yr}$ to equilibrate the photodissociation.
Before that, the reduced $\mathrm{CO}$ relative abundance could be a helpful
indicator limiting the evolution history of the gas. This constraint should be
confirmed with the measurement of $\mathrm{C^{+}}$ and $\mathrm{CH}$, which
are essential materials and intermediate products in the reaction chains for
$\mathrm{CO}$ formation.

Another promising probe is $\mathrm{OH}$, which is a sharp indication of
cosmic ionization of $\mathrm{H_{2}}$ (and atomic $\mathrm{O}$ as well, with
roughly half of the importance in the chemical network). Due to similar
ionization energy values for hydrogen and oxygen, efficient ionization can
only come from high-penetration sources, including the cosmic ray or the
external X-ray irradiation from, e.g., young protostars (not included in this
paper). Since the relative abundance of $\mathrm{OH}$ keeps climbing
throughout the $\sim 2\times 10^{5}~{}\mathrm{yr}$ period at the beginning of
evolution, the measurements could also help constrain the age of the system.
The reaction $\mathrm{H_{2}O^{+}}$ \+ $\mathrm{H_{2}}$ $\rightarrow$
$\mathrm{H_{3}O^{+}}$ \+ $\mathrm{H}$, one of the slowest steps in the
$\mathrm{OH}$ formation chain in such weakly ionized media, is insensitive to
temperature. In the meantime, the photodissociation of $\mathrm{OH}$ is
another limiting factor, as the process relies on the strong electric dipole
and is not susceptible to self or cross-shielding. Both factors make the final
product $\mathrm{OH}$ a possible indication of a few physical parameters,
including local density and the ISRF intensity. Consistent interpretations of
observations will require simulations for reasonable fittings on the
measurement results.

###  5.4 Future Works

Our current work focuses on numerical simulations, their intrinsic MHD
processes, and the properties based on the interactions between MHD, radiation
(ISRF and cosmic ray), and thermochemistry. In future works, we will elaborate
on possible predictions by synthesizing mock observations. One possibility is
the line intensities, either spatially resolved or integrated. Combined with
consistent MHD and thermodynamics, the spatial distribution features of
physical diffuse molecular cloud models can be converted into expectations on
actual astronomical measurements. Another possibility is the distribution
profiles in the position-position and position-position-velocity (PPV) spaces.
By confronting observables with synthetic profiles generated from an adequate
number of simulations, plausible interpretations of measurements will be
achieved via searching the parameter space, supported by proper microphysical
details that interact consistently with the dynamics and their observational
tracers.

N. Yue, L. Wang, T. Bisbas, and D. Quan are funded by Zhejiang Laboratory
Scientific Project # K2022PE0AB01. We also acknowledge the computing resources
provided by the Kavli Institute for Astronomy and Astrophysics in Peking
University. We thank our collegues for helpful discussions: Paul Goldsmith,
Xin Wang, and Haifeng Yang.

## References

  * Bai (2011) Bai, X.-N. 2011, ApJ, 739, 51, doi: 10.1088/0004-637X/739/1/51
  * Bai & Goodman (2009) Bai, X.-N., & Goodman, J. 2009, ApJ, 701, 737
  * Bian & Aluie (2019) Bian, X., & Aluie, H. 2019, Phys. Rev. Lett., 122, 135101, doi: 10.1103/PhysRevLett.122.135101
  * Boldyrev & Perez (2009) Boldyrev, S., & Perez, J. C. 2009, Phys. Rev. Lett., 103, 225001, doi: 10.1103/PhysRevLett.103.225001
  * Cazaux et al. (2017) Cazaux, S., Martín-Doménech, R., Chen, Y. J., Muñoz Caro, G. M., & González Díaz, C. 2017, ApJ, 849, 80, doi: 10.3847/1538-4357/aa8b0c
  * Clark et al. (2019) Clark, P. C., Glover, S. C. O., Ragan, S. E., & Duarte-Cabral, A. 2019, MNRAS, 486, 4622, doi: 10.1093/mnras/stz1119
  * Downes & O’Sullivan (2009) Downes, T. P., & O’Sullivan, S. 2009, ApJ, 701, 1258, doi: 10.1088/0004-637X/701/2/1258
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press)
  * Draine & Bertoldi (1996) Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269, doi: 10.1086/177689
  * Federrath et al. (2021) Federrath, C., Klessen, R. S., Iapichino, L., & Beattie, J. R. 2021, Nature Astronomy, 5, 365, doi: 10.1038/s41550-020-01282-z
  * Franeck et al. (2018) Franeck, A., Walch, S., Seifried, D., et al. 2018, MNRAS, 481, 4277, doi: 10.1093/mnras/sty2507
  * Glover & Jappsen (2007) Glover, S. C. O., & Jappsen, A. K. 2007, ApJ, 666, 1, doi: 10.1086/519445
  * Goldsmith et al. (2010) Goldsmith, P. F., Velusamy, T., Li, D., & Langer, W. D. 2010, ApJ, 715, 1370, doi: 10.1088/0004-637X/715/2/1370
  * Grete et al. (2021) Grete, P., O’Shea, B. W., & Beckwith, K. 2021, ApJ, 909, 148, doi: 10.3847/1538-4357/abdd22
  * Grete et al. (2023) —. 2023, ApJ, 942, L34, doi: 10.3847/2041-8213/acaea7
  * Heays et al. (2017) Heays, A. N., Bosman, A. D., & van Dishoeck, E. F. 2017, A&A, 602, A105, doi: 10.1051/0004-6361/201628742
  * Hennebelle & Falgarone (2012) Hennebelle, P., & Falgarone, E. 2012, A&A Rev., 20, 55, doi: 10.1007/s00159-012-0055-y
  * Ingalls et al. (2011) Ingalls, J. G., Bania, T. M., Boulanger, F., et al. 2011, ApJ, 743, 174, doi: 10.1088/0004-637X/743/2/174
  * Jappsen et al. (2007) Jappsen, A. K., Glover, S. C. O., Klessen, R. S., & Mac Low, M. M. 2007, ApJ, 660, 1332, doi: 10.1086/513085
  * Jones & Downes (2011) Jones, A. C., & Downes, T. P. 2011, MNRAS, 418, 390, doi: 10.1111/j.1365-2966.2011.19491.x
  * Jones & Downes (2012) —. 2012, MNRAS, 420, 817, doi: 10.1111/j.1365-2966.2011.20095.x
  * Kim et al. (2018) Kim, J.-G., Kim, W.-T., & Ostriker, E. C. 2018, ApJ, 859, 68, doi: 10.3847/1538-4357/aabe27
  * Kim et al. (2019) —. 2019, ApJ, 883, 102, doi: 10.3847/1538-4357/ab3d3d
  * Kowal et al. (2007) Kowal, G., Lazarian, A., & Beresnyak, A. 2007, ApJ, 658, 423, doi: 10.1086/511515
  * Langer et al. (2010) Langer, W. D., Velusamy, T., Pineda, J. L., et al. 2010, A&A, 521, L17, doi: 10.1051/0004-6361/201015088
  * Lupi et al. (2021) Lupi, A., Bovino, S., & Grassi, T. 2021, A&A, 654, L6, doi: 10.1051/0004-6361/202142145
  * Mac Low (1999) Mac Low, M.-M. 1999, ApJ, 524, 169, doi: 10.1086/307784
  * Mac Low et al. (1998) Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Phys. Rev. Lett., 80, 2754, doi: 10.1103/PhysRevLett.80.2754
  * McElroy et al. (2013) McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565, doi: 10.1146/annurev.astro.45.051806.110602
  * Nehmé et al. (2008) Nehmé, C., Gry, C., Boulanger, F., et al. 2008, A&A, 483, 471, doi: 10.1051/0004-6361:20078373
  * Neufeld & Kaufman (1993) Neufeld, D. A., & Kaufman, M. J. 1993, ApJ, 418, 263, doi: 10.1086/173388
  * Neufeld et al. (1995) Neufeld, D. A., Lepp, S., & Melnick, G. J. 1995, ApJS, 100, 132, doi: 10.1086/192211
  * Omukai et al. (2010) Omukai, K., Hosokawa, T., & Yoshida, N. 2010, ApJ, 722, 1793, doi: 10.1088/0004-637X/722/2/1793
  * Ostriker et al. (1999) Ostriker, E. C., Gammie, C. F., & Stone, J. M. 1999, ApJ, 513, 259, doi: 10.1086/306842
  * Pantaleone et al. (2021) Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2021, ApJ, 917, 49, doi: 10.3847/1538-4357/ac0142
  * Pineda et al. (2010) Pineda, J. L., Velusamy, T., Langer, W. D., et al. 2010, A&A, 521, L19, doi: 10.1051/0004-6361/201015089
  * Porter et al. (2015) Porter, D. H., Jones, T. W., & Ryu, D. 2015, ApJ, 810, 93, doi: 10.1088/0004-637X/810/2/93
  * Qian et al. (2018) Qian, L., Li, D., Gao, Y., Xu, H., & Pan, Z. 2018, ApJ, 864, 116, doi: 10.3847/1538-4357/aad780
  * Raskutti et al. (2016) Raskutti, S., Ostriker, E. C., & Skinner, M. A. 2016, ApJ, 829, 130, doi: 10.3847/0004-637X/829/2/130
  * Raskutti et al. (2017) —. 2017, ApJ, 850, 112, doi: 10.3847/1538-4357/aa965e
  * Rice et al. (2018) Rice, J. S., Federman, S. R., Flagey, N., et al. 2018, ApJ, 858, 111, doi: 10.3847/1538-4357/aabae7
  * Richings et al. (2014a) Richings, A. J., Schaye, J., & Oppenheimer, B. D. 2014a, MNRAS, 440, 3349, doi: 10.1093/mnras/stu525
  * Richings et al. (2014b) —. 2014b, MNRAS, 442, 2780, doi: 10.1093/mnras/stu1046
  * Stone et al. (2020) Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4, doi: 10.3847/1538-4365/ab929b
  * Tilley & Balsara (2011) Tilley, D. A., & Balsara, D. S. 2011, MNRAS, 415, 3681, doi: 10.1111/j.1365-2966.2011.18982.x
  * Wang et al. (2019) Wang, L., Bai, X.-N., & Goodman, J. 2019, ApJ, 874, 90, doi: 10.3847/1538-4357/ab06fd
  * Wang & Dai (2021) Wang, L., & Dai, F. 2021, ApJ, 914, 98, doi: 10.3847/1538-4357/abf1ee
  * Wang & Goodman (2017) Wang, L., & Goodman, J. 2017, ApJ, 847, 11.  http://stacks.iop.org/0004-637X/847/i=1/a=11
  * Whitworth & Jaffa (2018) Whitworth, A. P., & Jaffa, S. E. 2018, A&A, 611, A20, doi: 10.1051/0004-6361/201731871
  * Wurster et al. (2019) Wurster, J., Bate, M. R., & Price, D. J. 2019, MNRAS, 489, 1719, doi: 10.1093/mnras/stz2215
  * Wurster et al. (2016) Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037, doi: 10.1093/mnras/stw013
  * Wurster & Rowan (2024) Wurster, J., & Rowan, C. 2024, MNRAS, 528, 2257, doi: 10.1093/mnras/stae090
  * Xu & Bai (2016) Xu, R., & Bai, X.-N. 2016, ApJ, 819, 68, doi: 10.3847/0004-637X/819/1/68

##  Appendix A Hydrostatic Models

Figure 13 presents the two tests’ most abundant and essential chemical
species. These test runs are identical to the fiducial model except that they
are static (i.e., have no turbulence injection), and one of them has zero
cosmic ray intensity as well (the lower panel). In order to verify that our
simulation system—already verified in various non-ideal MHD simulations for
PPDs and planetary atmospheres—is also functioning correctly for the diffuse
molecular clouds, this appendix elaborates on the consistency of a few key
species by comparing the simulated key values of chemical kinetics with
analytic results.

We first examine the species related to hydrogen. Both models exhibit the
dissociation of hydrogen molecules over the $10^{6}~{}\mathrm{yr}$ evolution,
at the end of which the abundance of atomic $\mathrm{H}$ reaches $\sim
10^{-2}$ due to the photodissociation by the ISRF. As the self-shielding of
$\mathrm{H_{2}}$ is considered with the effective shielding length $l_{\rm
cap}\sim 1~{}\mathrm{pc}$ on a uniform medium (see also §2.2) with
$\rho=250~{}m_{p}~{}\mathrm{cm}^{-3}$ (or $n(\mathrm{H_{2}})\simeq
200~{}\mathrm{cm}^{-3}$) and $T=30~{}\mathrm{K}$, one gets the self-shielding
factor $f_{\mathrm{H_{2}}}\simeq 2.5\times 10^{-5}$ using the recipes in
Draine & Bertoldi (1996). This leads to a $\zeta_{\rm
diss,\mathrm{H_{2}}}\simeq 3\times 10^{-16}~{}\mathrm{s}^{-1}$
photodissociation rate, yielding $n(\mathrm{H})\simeq
1.9~{}\mathrm{cm}^{-3}\simeq 10^{-2}~{}n_{\mathrm{H}}$. The free electrons are
predominantly provided by the ISRF ionization of $\mathrm{C}$, whose abundance
almost equals the elemental abundance of carbon. The $\mathrm{H^{+}}$
concentration is eliminated by the radiative recombination (whose rate
coefficient is $\alpha_{\mathrm{H}}\simeq 2\times
10^{-11}~{}\mathrm{cm}^{3}~{}\mathrm{s}^{-1}$ at $30~{}\mathrm{K}$) at the
timescale $\tau_{\mathrm{H^{+}}}=(n_{e}\alpha_{\mathrm{H}})^{-1}\simeq
4.6\times 10^{4}~{}\mathrm{yr}$, in the absence of cosmic rays. With
$\zeta_{\rm CR}=10^{-17}~{}\mathrm{s}^{-1}$, radiative recombination maintains
$n(\mathrm{H^{+}})\simeq n_{\mathrm{H}}\zeta_{\rm
CR}\tau_{\mathrm{H^{+}}}\simeq 1.2\times 10^{-3}~{}\mathrm{cm}^{-3}$. These
results agree quantitatively with the test models.

The formation of $\mathrm{CO}$ follows a chain of reactions from
$\mathrm{C^{+}}$ to $\mathrm{CH}$ (see e.g., Draine 2011), starting from
$\mathrm{C^{+}}+\mathrm{H_{2}}\rightarrow\mathrm{CH_{2}^{+}}$. At
$30~{}\mathrm{K}$ and the 1 abundances, this reaction has a $\sim 4.6\times
10^{-15}~{}\mathrm{cm}^{-3}~{}\mathrm{s}^{-1}$ production rate of
$\mathrm{CH_{2}^{+}}$, a quarter of which is instantly converted to
$\mathrm{CH}$ by dissociative recombination with $e^{-}$ (equivalent to a
$\mathrm{CH}$ production rate of $\sim 1.2\times
10^{-15}~{}\mathrm{cm}^{-3}~{}\mathrm{s}^{-1}$). The final step
$\mathrm{CH}+\mathrm{O}\rightarrow\mathrm{CO}+\mathrm{H}$ competes with
$\mathrm{CH}+\mathrm{H}\rightarrow\mathrm{C}+\mathrm{H_{2}}$. The sum of these
two reactions consume $\mathrm{CH}$ at the timescales $\tau_{\mathrm{CH}}\sim
1.3\times 10^{11}~{}\mathrm{s}$, constraining the $\mathrm{CH}$ abundance to
be $n(\mathrm{CH})\simeq 1.6\times 10^{-4}~{}\mathrm{cm}^{-3}$, and leading to
a production rate of $\mathrm{CO}$ of $\sim 6\times
10^{-16}~{}\mathrm{cm}^{-3}~{}\mathrm{s}^{-1}$. Consider the ISRF
photodissociation rate $\zeta_{\mathrm{CO}}\simeq 4\times
10^{-11}~{}\mathrm{s}^{-1}$ (with $0.3~{}G_{0}$ and the total shielding factor
$\sim 0.57$ estimated using the recipes in Heays et al. 2017), such
$\mathrm{CO}$ production rate leads to $n(\mathrm{CO})\simeq 1.5\times
10^{-5}~{}\mathrm{cm}^{-3}$. The slight drop of $\mathrm{CO}$ abundance since
$t\sim 10^{5}~{}\mathrm{yr}$ is attributed to the increase of atomic
$\mathrm{H}$. In addition to confirming species abundances with Figure 13, we
also verify that the analytic estimates of reaction rates are quantitatively
consistent with the outputs given by the test simulations (not shown in the
figures).

Figure 13: Similar to the Panel (b) of Figure 1 but for the two test
simulation models under hydrostatics (viz. no turbulences or any kinematics).
The two panels present the results with and without cosmic rays, respectively.