# Necessary conditions for the formation of filaments and star clusters in the
cold neutral medium

Rachel Pillsworth1 , Ralph E. Pudritz1,2,3,4
1Department of Physics & Astronomy, McMaster University, 1280 Main St West,
Hamilton, ON L8S 4L8, Canada
2 Origins Institute, McMaster University, Hamilton, ON L8S 4M1
3 Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische
Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany
4 Max-Planck Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg,
Germany Contact e-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Star formation takes place in filamentary molecular clouds which arise by
physical processes that take place in the cold, neutral medium (CNM). We
address the necessary conditions for this diffuse ($n\approx 30$ cm-3), cold
(T $\approx$ 60 K), magnetized gas undergoing shock waves and supersonic
turbulence, to produce filamentary structures capable of fragmenting into
cluster forming regions. Using RAMSES and a magnetized CNM environment as our
initial conditions, we simulate a 0.5 kpc turbulent box to model a uniform gas
with magnetic field strength of 7 $\mu G$, varying the 3D velocity dispersion
via decaying turbulence. We use a surface density of $320M_{\odot}pc^{-2}$,
representative of the inner 4.0 kpc CMZ of the Milky Way and typical luminous
galaxies. Filamentary molecular clouds are formed dynamically via shocks
within a narrow range of velocity dispersions in the CNM of 5 - 10 km/s with a
preferred value at 8 km/s. Cluster sink particles appear in filaments which
exceed their critical line mass, occurring optimally for velocity dispersions
of 8 km/s. Tracking the evolution of magnetic fields, we find that they lead
to double the dense star forming gas than in purely hydro runs. Perpendicular
orientations between magnetic field and filaments can increase the accretion
rates onto filaments and hence their line masses. Because magnetic fields help
support gas, MHD runs result in average temperatures an order of magnitude
higher than unmagnetized counterparts. Finally, we find magnetic fields delay
the onset of cluster formation by $\propto 0.4$ Myr.

###### keywords:

MHD – turbulence – methods:numerical – stars:formation – ISM:clouds –
ISM:structure

††pubyear: 2022††pagerange: Necessary conditions for the formation of
filaments and star clusters in the cold neutral medium–C

##  1 Introduction

Molecular clouds represent the intermediate galactic scales (100s of pc) at
which large scale galactic structure and evolution gives way to more localized
star formation in forming star clusters. The multiple scale nature of this
process has been explored within molecular clouds by the Herschel space
observatory, and larger galactic scale studies of molecular cloud populations
by PHANGS-ALMA. Over a decade ago, Herschel made a major breakthrough by
revealing the ubiquitous filamentary networks that comprise molecular clouds
in which star-forming cores form on subpc scales (Alves de Oliveira et al.,
2014; André et al., 2010; Henning et al., 2010; Men’shchikov et al., 2010).
PHANGS is a recent addition to the research, and produced the largest
catalogue (100,000) of individually resolved molecular clouds across multiple
galaxies (Lee et al., 2022; Sun et al., 2022; Turner et al., 2022; Brunetti &
Wilson, 2022). It reveals the connection between molecular cloud properties
and their galactic environments across a large array of galaxies (Thilker et
al., 2021).

Molecular clouds are believed to condense out of a phase of the interstellar
gas known as the cold neutral medium (CNM) in which they are embedded (Klessen
& Glover, 2016; Chevance et al., 2022). It is therefore essential to
understand the dynamics and thermal evolution of the CNM - in particular under
what conditions molecular clouds and their star clusters may form in it. For
instance, Sun et al. (2022) use the PHANGS survey to show that molecular cloud
properties depend on local environmental properties such as gas surface
density and star formation density. On the galactic scale, they find it is
primarily galactic shear that affects the final properties of molecular
clouds. These recent surveys have given a much better picture of the
properties of entire populations of GMCs both within and outside the Milky
Way, although they are not yet sufficient to penetrate scales below the star
cluster.

Many of the galactic observations are well supported by theoretical and
computational work. As one example, the results in Sun et al. (2022) compare
well to the simulations of Jeffreson et al. (2020), in that both studies show
that cloud properties depend on some local and some galactic environmental
properties. Jeffreson et al. (2020) showed that rotation in the galactic disk
influenced the geometry and rotation of the molecular clouds which formed
within it. Conversely, their work finds that gravo-turbulent and star
formation properties of the molecular clouds are somewhat disconnected from
the overall galactic disk dynamics, supporting an agreement with Sun et al.
(2022) that the gas and star formation densities are tied to local
environment. Additionally, simulations from Grisdale (2021) showed that
different global star formation prescriptions produce different local effects,
especially along spiral arms in galactic disks, with turbulent star formation
prescriptions forming stars in high density pockets located in the arms. These
different environments also affect the final masses, radii and velocity
dispersions of the GMCs that form. Finally, observations from Soler et al.
(2022) studied giant HI filaments towards the Galactic plane finding that the
filamentary structure is closely linked to the properties of the Galactic
disk, being formed as direct consequences of galactic dynamics. The
interstellar medium (ISM) throughout a galaxy is also repeatedly shocked,
heated, and compressed by turbulence and expanding superbubbles (Watkins et
al., 2023a, b).

Molecular clouds host star clusters embedded in filamentary structures
(Krumholz et al., 2014; Grudić et al., 2021). These are the parent structures
of most of the stars (70-90 %) that form in GMCs (Lada & Lada, 2003) and makes
them a natural starting point for studies of star formation. On these parsec
and sub-parsec scales of clusters, the observational work has also been rich
and varied. Telescopes such as the Chandra X-Ray Observatory, Hubble and
ALMA/SMA have allowed us to peer into Milky Way clouds and their clusters, and
investigate star formation within them (Townsley et al., 2014; Massey &
Hunter, 1998; Galván-Madrid et al., 2013; He et al., 2022).

The ionizing radiation from massive stars will also contribute to the heating
in the star cluster, preventing low-mass stars from forming and effectively
setting the IMF (Klassen et al., 2012). Protostellar jets and outflows from
massive stars can relieve radiation pressure around the protostar and set
final masses, allowing them to have a higher final mass than without outflows
(Krumholz et al., 2005; Kuiper & Hosokawa, 2018). These active areas of star
formation can extend out many parsecs from the protostar, affecting the
environment they’re embedded within (Beuther et al., 2018) and, therefore, the
local star formation rate. Furthermore, Kuffmeier (2018) used zoom-in
simulations to investigate star formation and find that protostellar accretion
depends on the environment in which the protostar is located, firmly linking
cluster environment and protostar activity.

The feedback effects of these processes on GMC formation have also been
extensively studied. For instance, recent work has determined that the tides
tearing gravitationally bound clouds apart can come, in part, from massive
star-forming regions within the cloud (Ramírez-Galeano et al., 2022). Howard
et al. (2018) found that the most massive clusters tend to form in high
density filaments in the most massive GMCs that form, while others discuss
that radiation feedback plays an essential role in the lifetime of GMCs (see
Chevance et al., 2022, and references therein). Other works have discussed the
role of turbulence in GMC formation, finding it essential for filament
formation (Federrath et al., 2010a), and support against gravitational
collapse. Supernova explosions within clouds have been found to play a less
important role in GMC destruction, as the first supernovae explosions seem to
go off too late to affect their GMC, rather expanding into the local ISM
(Smith et al., 2021). Recently, STARFORGE (Grudić et al., 2022) has combined
many of these processes into one simulation starting from GMC and
investigating individual stars.

In considering regions in which most molecular clouds form in galaxies,it is
well-known that our own Milky Way does not contain many massive clouds (i.e.,
greater than $10^{6}$ M⊙). In fact, CO survey catalogues of the entire
Galactic plane have found only 1064 massive clouds (Rice et al., 2016), with
the outer disk’s upper mass limit reaching only $10^{6}$ M⊙. The same work
finds that those more massive clouds are found in the inner disk of the
Galaxy, closer to the dense, active Central Molecular Zone where higher line-
widths suggest a more direct role from the galactic disk in the formation of
clouds. Indeed, much higher surface densities, as can be found in luminous or
starburst galaxies, are needed in order to form the massive star clusters that
give rise to giant OB associations. These galaxies host 1-2 orders of
magnitude more molecular gas than our Milky Way and are capable of forming
molecular clouds 1000 times more massive than GMCs in the Milky Way (Solomon,
2001). Therefore the local neighbourhood in the Milky Way isn’t ideal for
simulating the processes that give rise to the most massive or numerous clouds
in the galaxy. For the massive end of cluster and cloud formation, one
requires a region with a much higher gas surface density than our own local
galactic environment can represent. Thus, many of the previous studies we have
cited likely fall short in investigating the typical massive cluster’s
formation.

In this paper, we investigate the sensitivity of the CNM to turbulence,
gravity, and magnetic fields towards building self-gravitating giant molecular
clouds for environments in which the bulk of massive molecular clouds are
known to form, such as luminous starburst galaxies or the very active CMZ. In
particular, we simulate a $500^{3}$ pc patch of the ISM in a high gas surface
density environment to investigate the scaling relations dominated by ISM
processes in these regions. The CNM on these scales is the starting point for
the transition from diffuse, cold atomic gas to star-forming molecular gas. It
is the initial state in our simulations that are designed to isolate the
effects of CMN properties on structure formation. We present the argument for
connecting the two scales in the following Sec.2. We present our models and
the previous works that guided our setup in Sec.3. A discussion of each
important result from our models follows in Sec.4. Finally, we present our
conclusions and the future of this problem in Sec.5.

##  2 Bridging scales

Computationally, a plethora of different initial conditions for any scale from
the 100 kpc-scale galactic disk to the pc-scale clusters have given us ample
opportunity to investigate many of the intricacies of structure formation from
different origins (see, for example, Grudić et al., 2022; Howard et al., 2018;
Brown & Gnedin, 2022; Rieder et al., 2022; Lane et al., 2022; Lue et al.,
2021; Kim & Ostriker, 2017, for simulations). Yet, even with all this work,
there are lingering questions. For instance, is there a point we can
disconnect the galactic disk and the GMC? How much influence does the disk
have on star cluster formation, or the formation of individual massive stars
within them? How many scales do protostellar jets, as well as feedback from
massive stars, affect? We argue here that many of these questions can not be
answered in isolation from the other, but rather that they require us to
connect scales.

###  2.1 Filamentary Structures on Many Scales

On molecular cloud scales, André et al. (2014) presented a comprehensive view
- involving observations, theory, and simulations - on the central role played
by filaments in star formation through their fragmentation into stars via
gravitational instabilities. This picture shares the point made by the
gravoturbulent theories of star formation (Klessen & Glover, 2016), that
supersonic turbulence creates the filamentary structure in the ISM (Mac Low &
Klessen, 2004; Larson, 1981). However, it differs in emphasizing that when
pushed to sufficiently high mass per unit length filaments become
gravitationally unstable and fragment into star forming cores (Fiege &
Pudritz, 2000; Molinari et al., 2010; André et al., 2014). The authors make
the argument for the strengthening link between star formation and the
structure of the ISM in galactic disks, motivated by this new paradigm of star
formation in filaments. More recent reviews, such as Chevance et al. (2022),
further support this view of star formation by exploring the filamentary
environments we see surrounding GMCs from both observational and theoretical
work. In either case, and among many other groups, there is general agreement
that star formation is connected to filamentary structures in the (cold) ISM.

Observational studies strongly support the existence of filamentary structures
across multiple scales. In both HI and CO, systems of filaments have been
visible in the ISM, with the cold phase in particular being highly filamentary
(see McClure-Griffiths et al., 2006; Falgarone et al., 2001; André et al.,
2010; Henning et al., 2010; Motte et al., 2010). The CNM provides a suitable
initial condition for this work in a few ways. The gas in the ISM follows a
cycle. Starting with the warm gas, thermal instabilities cause the gas to cool
into the CNM and create giant, dense filaments (Hacar et al., 2022). Out of
those cold filaments, molecular gas begins to form (Rathjen et al., 2021),
which we can trace observationally using dust extinction due to dust growth
being fastest and most abundant in cold gas (Klessen & Glover, 2016). The CNM
hosts most of the filamentary structures, especially those leading to and
within dense molecular clouds (Seifried et al., 2017). The molecular hydrogen
forms on dust grains when it has reached a sufficiently high column density
that it can shield itself from UV radiation. In our galaxy, the density
required ($\sim 10^{21}cm^{-3}$) is also the density at which self gravity
becomes important. Through gravitational collapse and instabilities discussed
previously, molecular clouds create the bound clusters which host stars.
Feedback from the stars heat the surrounding gas, completing the cycle and
taking the gas from cold and molecular to its hot or warm phase once again
(Girichidis et al., 2016). Furthermore, this filamentary structure appears
well before star formation starts (André et al., 2014). These filaments,
independent of the observational tracer used, appear to be coherent
structures, joining with each other to form dense pockets (referred to as
hubs). These hubs are associated with groups of young stellar objects, or
young clusters (Myers, 2009), showing a physical connection between the cold
ISM and the birth environments of massive clusters.

Magnetic fields may be playing an important role in the formation and
evolution of filaments. The dense filaments show more linear structure and
tend to be perpendicular to magnetic field orientation (Soler, 2019). On the
other hand, lower density filaments are more scattered in their directions and
tend to align parallel with magnetic field orientation (see also Kwon et al.
(2022)).

We can see small-scale, bound, star-forming structures form in the hubs at
filament junctions, as well as along the densest filaments from gravitational
instabilities. The gravitational instabilities in a fragmenting filament come
about from the transition to self-gravitation, and can best be determined
through filament line masses. When in the ISM though, we must also consider
the physical affects of mixing and turbulence. In this case, the critical line
mass includes both thermal and non-thermal contributions to the velocity
dispersion:

$\sigma_{tot}=\sqrt{c_{s}^{2}+\sigma_{turb}^{2}}$ (1)

This gives us a virial line mass, as discussed in Fiege & Pudritz (2000),
defined as:

$M_{vir}=\frac{2\sigma_{tot}^{2}}{G}.$ (2)

For magnetized filaments, (including the possibility of a helical magnetic
field wrapping around a constant filament, in addition to the poloidal
threading field) there is a magnetic correction to apply to the critical line
mass calculation. Following from Equations 27 and 28 from Fiege & Pudritz
(2000), the critical line mass becomes

$M_{mag,vir}=f_{B}*M_{vir}=\frac{1+(v_{A}/\sigma)^{2}}{1+(v_{A,\phi}/\sigma_{c})^{2}}M_{vir}$
(3)

using the poloidal Alfvèn speed $v_{A}=\frac{B_{\parallel}}{\sqrt{4\pi\rho}}$,
toroidal Alfvèn speed $v_{A,\phi}=\frac{B_{\perp}}{\sqrt{4\pi\rho}}$, and
critical velocity dispersion $\sigma_{c}^{2}=4\pi G\rho r^{2}$. Because of the
explicit presence of both poloidal and toroidal field orientations in the
above equation, one can see their respective effects. The poloidal magnetic
increases the critical line mass, because it lends more pressure support to
the gas. A toroidal field, on the other hand, squeezes the filament thereby
decreasing the critical line mass. As such, it is the geometry of the magnetic
field which affects line mass more so than its presence alone.

As derived, the thermal critical line mass scales in cold molecular clouds as
$M_{crit}\approx 16M_{\odot}pc^{-1}\times\bigg{(}\frac{T_{gas}}{10K}\bigg{)}$.
This gives a general minimum critical line mass, with the value changing
depending on turbulent motions within denser molecular gas. In higher
temperature environments such as the CNM environment, with an average
temperature of 80K, we expect the average critical line mass to be
$M_{crit,CNM}=128M_{\odot}pc^{-1}$, in atomic gas filaments - almost an order
of magnitude larger than a cold molecular cloud at 10K. This suggests that
larger scale structures more reminiscent of small cluster masses, can be
expected to grow unstable within a CNM filament that is becoming
gravitationally interesting, but not yet condensing into molecular gas.

Supersonic turbulence in the galaxy, whatever it’s source, obeys a size-line
width relation $\sigma_{turb}\propto L^{1/2}$ (Hacar et al., 2022). This is
identical to the scaling of Burgers turbulence that arises from shocks. This
scaling, when taken up to kpc scales, give line masses of the order
$10^{4}\rm{M_{\odot}pc^{-1}}$ indicating that massive atomic filaments can
form molecular clouds via gravitational instabilities (Zhao et al., 2023).

While in theoretical examples the critical line mass gives us conditions for
collapse, it is also supported observationally, making it a universally
helpful criterion for structure formation in filaments. On the smallest scales
characteristic of individual star formation, observations of Herschel
filaments suggest that supercritical filaments are in virial equipartition and
gravitationally bound. The subcritical filaments, on the other hand, are
unbound and have transonic or subsonic velocity dispersions (Arzoumanian et
al., 2013). In fact, subcritical filaments hold less than one third of the
bound prestellar cores found in a molecular cloud, and those cores tend to be
less massive and less dense than their counterparts forming in dense filaments
(Polychroni et al., 2013). Furthermore, the more massive bound cores found in
filaments tend to be closer to intersections of filaments, supporting the idea
of hub sections being the preferential site of massive cluster formation.

The hubs funnel the flow of gas, becoming dense and massive, and providing
conditions for clustered star formation. These dense clumps will be the
preferential sites for massive cluster formation, due to their connection to
high rates of gas inflow from multiple filaments and their already dense
environments (Myers, 2009). Moreover, the Herschel Gould Belt Survey found
that the deeply embedded protostars (usually those which form massive stars)
are found in filaments with column densities higher than $7\times
10^{21}cm^{-2}$ (André et al., 2014). Even without limiting our view to only
massive clusters and cores, we see that the vast majority ($>70\%$) of
Herschel-identified prestellar cores are found within filaments as opposed to
outside of filamentary structures, suggesting (as seen in Men’shchikov et al.,
2010) that cores form along filaments via cloud fragmentation.

In addition to gravitational instabilities giving rise to star-forming cores,
the filaments act as conduits for gas flow into smaller scales. Observational
works (see Smith et al., 2012) have outlined an accretion flow along filaments
onto cores, growing their mass over time. They have also suggested that a
primary role of filaments is to help focus the gas towards embedded cores
(Gómez & Vázquez-Semadeni, 2014).

However, recent simulations (see Zhao et al., 2023) have found higher
accretion rates onto the filament than along them indicating more significant
building up of material. This is also observed for accretion onto and flow
along the filament in which the Serpens South cluster is forming (Kirk et al.,
2015), where infall rates onto filaments exceed flow rates along them by
nearly a factor of 10. Regarding massive cluster formation, Peretto et al.
(2013) have found it possible that massive protostars accrete the majority of
their mass from larger (filamentary) scales as opposed to just from their
prestellar core. Dense cores also tend to share certain kinematic properties
with their filaments, furthering the argument of gas flows along the
structures. It is therefore clear that filamentary structures in the ISM
connect to multiple scales, from the cloud down to the single protostar, and
the expected scales of instability in the filaments are approximately four
times the scale length.

###  2.2 The transition of the CNM to GMCs

Given the importance of the filamentary network in the ISM, we now must ask
what exactly a GMC is. Is there a distinct difference between ISM environment
and GMC? Chevance et al. (2022) argue that GMCs should not be thought of as
discrete entities, but rather as an observationally defined feature in gas and
dust maps. Perhaps the greatest representation of this fuzzy distinction are
maps of the Orion A GMC.

Figure 1: Image from Chevance et al. (2022). Orion A GMC dust emission map
(grayscale) and ${}^{12}CO(1-0)$ emission (contour lines in position-position
space.

In Fig.1, Chevance et al. show the dust emission and ${}^{12}CO(1-0)$ emission
of the Orion A cloud together. Both tracers show what we would expect of a
GMC: higher aspect ratio (>1), filamentary (or clumpy, in the case of CO
contours) structure and overall connection to its environment. The largest
distinction is with the resolution between the two maps. While the contour
lines depict a specific contour demarcating the GMC, the dust emission shows
much better the connection between GMC and environment through filamentary
networks. This comparison is further complicated by the fact that, while CO
emission certainly guarantees that the region is in a molecular phase of gas,
it is not necessarily a tracer of molecular hydrogen in GMCs (Pineda et al.,
2008), making the CO emission view of a GMC an estimate of one based on
resolution and limitations in the tracer itself. Thus, when looking only at
gas emission, it is natural to assume a GMC as a structure which we can
isolate from its environments in numerical studies.

Dust emission, on the other hand, traces the neutral ISM, specifically the
column density. Dust maps, being done with higher resolution and at sub-
millimetre wavelengths where dust emission is optically thin, often allow us
to resolve small-scale structure in our GMCs. With the benefit of recent GAIA
surveys, dust emission is now also used to map out 3D structure of these
clouds in a way that gas tracers cannot (Rezaei Kh. et al., 2020; Zucker et
al., 2021). Thus, dust emission maps suggest that GMCs are inherently
connected to their surrounding environment and, therefore, should not be
isolated from their galactic disks to simulate their formation.

To summarize, the connection between filamentary networks and star clusters is
strong. Furthermore, the GMC as a structure in molecular hydrogen is certainly
a distinct step in the transition from neutral atomic gas to cold star-forming
gas. In galactic scale simulations of cluster formation, the cold phase of the
ISM gas is crucial in setting realistic cluster mass functions (CMF) on par
with the Milky Way (Reina-Campos et al., 2022). It is clear that dynamical
processes in the CNM are important for the formation and final properties of
massive star clusters and vice versa. This motivates our choice of conditions
in the cold ISM as the initial condition of our work.

##  3 Numerical Methods

Our simulations use Ramses, which is a magneto-hydrodynamics code (Teyssier,
2002) popularly used for cosmological simulations (for an example of star
formation research, see Calura et al., 2022; Decataldo et al., 2017; Zhao et
al., 2023; Brucy & Hennebelle, 2021). On the smaller ISM scale and below, the
code has become increasingly popular to use due to its capability for high
resolution and efficient run times, arising from its implementation of an
adaptive time step. For example Bellomi et al. (2020); Han et al. (2022);
Ntormousi & Hennebelle (2019); Brucy et al. (2020, 2023) all use Ramses at
parsec scales to investigate star formation with high resolution. Given the
increasing use of the code on smaller scales, and the link this paper has to
our larger galactic simulations (also done in Ramses), we choose to use it for
our ISM scale work, as it sits between the aforementioned scales already
tested. In order to achieve its high resolutions, Ramses uses adaptive mesh
refinement, an important mechanism that we implement in our work in order to
resolve our filamentary structure.

###  3.1 Physical Mechanisms

When investigating the CNM and molecular cloud formation, turbulence plays an
important role in creating the structure. As we have already discussed, it is
an essential addition to the consideration of critical line masses,
influencing the gravitational collapse of filaments and thus pushing the
subsequent structures that form to higher masses. Turbulent mixing from
supernovae explosions specifically is most crucial in the cold and warm phases
of the ISM (Kim & Ostriker, 2017), where we see velocity dispersions of 5-12
$kms^{-1}$ in simulation domains of size $0.5\times 0.5\times 2$ kpc.
Dispersions of this magnitude are often associated with turbulence, as the
typical turbulent rms velocity is $v_{rms}\approx 5\frac{km}{s}$, which
becomes the dominant speed compared to typical sound speeds of $c_{s}\approx
0.2-1.0\frac{km}{s}$. Large-scale turbulent motions have also been found to
have significant effects on the structure of the ISM at kpc scales, with the
turbulent power spectra of the filaments and molecular clouds containing
signatures from an imposed larger-scale power spectrum (Colman et al., 2022).

Generally, magnetic fields are difficult to measure, especially in GMCs, due
to the fact that most GMCs sit in or within only a few degrees of the Galactic
mid-plane which makes line-of-sight magnetic field observations next to
impossible with current technology (Pattle et al., 2022). Yet, despite this
fact, we have been able to uncover a substantial amount of information
regarding the effects of magnetic fields from theoretical simulations as well
as our relatively limited observations. Magnetic fields link gas across
multiple physical scales, being most dynamically important around the pc-scale
of a molecular cloud. On scales of 10 pc, molecular clouds have highly ordered
magnetic fields, with the orientation of their internal structure closely
aligning with the orientation of the fields (Pillai et al., 2020; Pattle et
al., 2022; Tahani et al., 2022). These highly ordered fields are indications
that the field is strong enough to resist distortions due to turbulence. They
also provide support against compressive shocks, as well as general pressure
support against collapse, thus delaying the onset of molecular cloud formation
by more than 20 Myr (Girichidis et al., 2018). On the other hand, though it
will take longer to form molecular clouds, fragmentation in the ISM will
happen sooner with strong magnetic fields due to the action of global Parker
instability modes, and support long filaments that extend in either radial or
azimuthal directions in the disk, as opposed to hydro cases which see
predominantly ring-like structures in the filaments (Körtgen et al., 2019).

Cluster formation will also be affected by radiation effects and galactic
shear. Radiative feedback and stellar winds from massive stars will control
the star formation rate through their influence on local gas properties, thus
forcing clusters that form after to be lower mass(Rathjen et al., 2021).
Galactic shear has a loose correlation with local gas properties around
molecular clouds (Sun et al., 2022), while also having a moderately strong
effect on the overall rotation of a cloud and, therefore, the angular momentum
available to transfer from cloud to cluster (Jeffreson et al., 2020).

While we recognize the importance of both of these effects, we neglect them in
this work for simplicity, as we focus on short times in the early formation
and evolution of molecular clouds such that neither shear nor radiation will
have significant effects on the dynamics of the gas. Our goal here is to focus
on the sensitivity of molecular cloud and cluster formation to dynamical
processes that shape the CNM - magnetic fields, turbulence, and gravity. The
level of turbulence in the galaxy may be regulated by star formation feedback
to a value that is able to offset global gravitational collapse (Kim &
Ostriker, 2007, 2015). The question we address here is what do different
levels of CNM turbulence imply for filament and molecular cloud formation?

In a separate study of molecular clouds in a galactic disk, presented in Zhao
et al. (2023), we include full galactic environmental effects, including
galactic shear, and investigate the role of galactic shear on the clouds. In
future work from that paper, we will further resolve the cluster scale and
compare to the present work, including the effects of radiation.

###  3.2 Simulations

####  3.2.1 Initial Conditions

We simulate the CNM via a turbulent box setup for two general cases; a
magnetic and non-magnetic case. We start with a 0.5kpc $\times$ 0.5kpc
$\times$ 0.5kpc box, to match the size of the PHANGS-ALMA hexagonal kpc-scale
observations (Sun et al., 2022) as well as being on par with similar
theoretical works (cf. Bellomi et al., 2020; Colman et al., 2022; Kim &
Ostriker, 2017). This size is chosen such that it can contain entire cloud
complexes, with GMCs on the scale of 70 pc, as well as the properties and
physical mechanisms acting in the local environment, such as those discussed
above. However, we do note that this will not account for the gravitational
potential in a galaxy due to its stars and dark matter, instead only allowing
a central collapse by assuming no initial distribution in the potential. The
origin is set in the center of the box, such that half of the z-dimension of
our box can be imagined as containing 250 pc both above and below the plane,
similar to pressure scale heights of CNM gas (Körtgen et al., 2019; Soler et
al., 2022).

We choose a domain size and surface density to align with the lower end of
Figure 2a from Wilson et al. (2019), such that a domain of half height 250 pc
places our CNM gas at the less luminous end, more similar to that which would
be present in a starburst galaxy. Scatter in these numbers show scale heights
of up to 316 pc, so our choice of 250 pc places us within a very reasonable
assumption of scale height. This idea is further supported in observational
works which find scale heights of 270 pc minimum for low SFE and full vertical
disk thicknesses greater than 0.5 kpc (Fisher et al., 2019; Bournaud et al.,
2009). Additionally, previous works discussed in Fisher et al. (2019) argue
that other methods of kinetic energy injection, such as supersonic turbulence,
would increase these scale heights even more. So our domain size is perfectly
placed within reasonable values for our gas surface density. We discuss our
choice of gas surface density further on.

While a larger box size would be possible with a more explicit disk setup in
the density profile, and simulations of collimated flows are popular for this
very reason (see, for example, the SILCC projects: Girichidis et al., 2016,
2018; Rathjen et al., 2021; Seifried et al., 2017), they require additional
manual setup in the initial conditions to the velocity and angular momentum.
This is something we wished to avoid in our project, as we wanted to isolate
some of the most basic physical processes that regulate the CNM and under
which conditions it gives rise to cluster forming filaments. Due to both the
possible vertical extents of gas and the simplicity we aimed for in our
initial conditions, we keep the full extent of our z-direction 0.5 kpc, making
it symmetric with the rest of the box.

Observations of pc scale star-forming filaments in molecular clouds have shown
typical widths of 0.1 pc, scaling with gas temperature (Arzoumanian et al.,
2011; Pineda et al., 2011, 2010). Given these results, we allow our box to
refine up to a highest resolution of 0.48 parsec per cell, which is achieved
consistently throughout our dense filaments, so that we highly resolve our
clouds but do not resolve individual star-forming filaments and our sink
particles can act to represent protoclusters. Our simulations employ periodic
boundary conditions. This allows for some effects of galactic gas flow through
the simulation volume, without additional computational expense. We do not
include shear across our box. Thus, local angular momentum arises purely from
the local angular momentum that comes about from turbulent mixing and shocks
moving through the medium.

Since we do not include shear or any effects of galactic rotation in this box,
we use a mixture of solenoidal and compressive turbulence to simulate general
motions and mixing in a patch of the ISM and allow for filament formation
within our box. The turbulence is initialized with a 1/3 compressive fraction.
We allow our turbulence only to decay, and keep driving forces turned off such
that the turbulence is never evolved past initialization, as motivated by
previous work such as Howard et al. (2018). Finally, the power spectrum
follows a Burger’s power-law shape as supported by previous studies of
supersonic turbulence in the ISM and GMCs (Federrath et al., 2010a; Kitsionas
et al., 2009), given by the energy spectrum $E(k)$ of the turbulence, where k
is the wavenumber:

$E(k)\propto k^{-2}$ (4)

wherein the velocity dispersion is $\sigma_{turb}^{2}(k)=\int E(k)dk\propto
k^{-1}$.

In setting up the turbulence in a simulation, RAMSES has effectively two
control parameters that need to be set for the strength of the initial
turbulence field: the source term of the kinetic energy (called the "RMS
forcing parameter"), and the decay time scale. By selecting various decay
timescales of our turbulence, the initial velocity dispersion measured in our
simulation can be fixed.

Given the use of decaying turbulence, one must be careful with the decay
timescale and the timescale of structure formation in the simulation. In order
to achieve the 3D velocity dispersions we find necessary in this work, the
decay timescales must be set sufficiently long enough to maintain the desired
dispersion for structure formation. We then run our simulations for a maximum
of 10% of their decay time in order to ensure turbulence has not significantly
decayed throughout formation, as well as to remain sufficiently early enough
in protocluster formation that our assumptions surrounding the lack of
radiative transfer can remain valid. As such, decay timescales are set to be
100, 65 and 50 Myr which create dispersions of 5 km/s, 8 km/s and 10 km/s in
our models, respectively.

Figure 20 in Appendix C shows the turbulent energy spectra for times of 1.90
and 3.79 Myr in our disp8MHD simulation. Early times show the turbulent
cascade has not reached Burger’s scaling for the majority of the scales and
oscillates around it at small scales. The spectrum at 3.79 Myr shows the model
has reached Burger’s scaling at large and small scales, while the middle of
the spectrum trends closer to Kolmogorov scaling. This transition to
Kolmogorov scaling is expected for later times as magnetic energy increases
(Schleicher et al., 2013). The differences between early and late time do not
correspond to a significant decay fraction in the energy spectra as there are
no appreciable differences in upper energy magnitudes nor any deviation away
from Burger’s scaling at the small scales with time. In reality, the
simulations run up to the point their oldest sink particle reaches 0.5 Myr in
age, in order to avoid feedback and radiation effects that are here neglected.
Finally, it is this decay time scaling which affects the velocity dispersions
in our runs. We stress that velocity dispersion is not a parameter in the
setup of our simulation explicitly, and the values we are studying are direct
consequences of the strength of the turbulence initialized in the simulations.

Run | Density($H/cm^{3}$) | Box Size (kpc) | Initial Temperature (K) |  $\Delta x_{start}$ (pc) |  $\Delta x_{final}$ (pc) |  $\sigma_{turb,initial}$ (km/s) |  $\sigma_{turb,sink}$ (km/s) | B field ($\mu G$) |  $\lambda_{J,thermal}$ (pc) |  $\lambda_{J,total}$ (pc)
---|---|---|---|---|---|---|---|---|---|---
disp5 | 30.0 | (0.5,0.5,0.5) | 58.8 | 1.95 | 0.48 | 2.72 | 5.51 | 0.0 | 17.44 | 48.23
disp8 | 30.0 | (0.5,0.5,0.5) | 58.8 | 1.95 | 0.48 | 6.96 | 6.99 | 0.0 | 17.44 | 123.3
disp10 | 30.0 | (0.5,0.5,0.5) | 58.8 | 1.95 | 0.48 | 10.5 | 11.5 | 0.0 | 17.44 | 186.2
disp5MHD | 30.0 | (0.5,0.5,0.5) | 58.8 | 1.95 | 0.48 | 2.72 | 3.2 | 7.0 | 12.16 | 48.23
disp8MHD | 30.0 | (0.5,0.5,0.5) | 58.8 | 1.95 | 0.48 | 6.96 | 7.0 | 7.0 | 12.16 | 123.3
disp10MHD | 30.0 | (0.5,0.5,0.5) | 58.8 | 1.95 | 0.48 | 10.5 | 10.26 | 7.0 | 12.16 | 186.2

Table 1: Parameters for all models computed. The first three simulations
represent our fiducial models with no magnetic fields and maximum resolution
of 0.48 pc. We set our magnetic field runs to have a field strength of 7 $\mu
G$, in accordance with average values of ISM magnetic field strengths.
$\sigma_{turb,initial}$ and $\sigma_{turb,sink}$ represent the velocity
dispersion in the simulations on initialization and once sink particles have
formed, respectively. $\lambda_{J,thermal}$ and $\lambda_{J,total}$ give the
initial Jeans lengths of each simulation calculated using the average
temperature and the velocity dispersion, respectively.

Table 1 lists the initial conditions we use in each of our models, and we
outline here the decisions behind these parameters. In order to mimic the
magnetized CNM we initialize a constant magnetic field of magnitude $7\mu G$
in only the y direction, allowing it to evolve with the gas over time. We
choose this setup because large scale ordered magnetic fields in galaxies are
toroidal in the plane, for which a small patch will look approximately
constant and oriented in one direction. The gas is initially isothermal and
set at a density of 30 $cm^{-3}$ and a temperature of 58 K, comparable to
median values for the CNM, as outlined in Table 1 of Klessen & Glover (2016).
We then evolve the chemistry as described in the following subsection. These
values give surface densities of $320M_{\odot}pc^{-2}$, which are not
dissimilar to surface densities within the first few hundred parsec of the
Milky Way disk, otherwise known as the central molecular zone (CMZ). Gas
surface densities in this portion of the disk can be on average
$10^{2}-10^{3}$ M⊙ pc-2 for the atomic gas, especially in its early evolution
(see Figures 5 and 9 of Tress et al., 2020). Additionally, this chosen surface
density reflects the more active starburst and luminous galaxies, in which one
expects higher gas densities, on the order of $10^{2}$-$10^{4}$ M⊙ pc-2, due
to the much more star-forming regions these galaxies represent (Wilson et al.,
2023, 2019).

####  3.2.2 Heating & Cooling

Through the chemistry code Grackle (Smith et al., 2017), we include non-
equilibrium chemistry with a 9-species network to form both atomic and
molecular hydrogen. Grackle includes chemical heating and cooling, as well as
metal line cooling, radiative cooling and $H_{2}$ formation on dust grains
(which introduces dust cooling and dust-gas heat transfer). Through the
addition of molecular hydrogen formation, Grackle also initializes a default
ISRF strength of $2.72\times 10^{-3}~{}erg~{}s^{-1}~{}cm^{-2}$, but we do not
include radiative transfer from ramses nor do we include a UV background. Due
to the lack of UV background and radiation in the simulation, explicit self-
shielding of $H_{2}$ is ignored.

A constant photoelectric heating rate is used such that at a hydrogen number
density of n=1 cm-3, the heating rate is $8.5\times 10^{-26}$ erg cm-3 s-1.
Because we do not include any radiative transfer due to computational
limitations, there is no photodissociation. Our molecular hydrogen abundance
therefore represents an upper limit, as destruction of the molecule is not
accounted for in the abundances.

In our simulations we argue this is reasonable, as the photodissociation rate
for gas at temperatures below 10000 K is sufficiently low as to be reasonably
ignored in the vast majority of our gas (Coppola et al., 2011). For gas below
densities of 104 cm-3, the probability of dissociation via UV is only 15% and,
even in those higher density areas, a sufficiently high column density can
contribute to an order of magnitude drop in the dissociation rate (such column
densities translate to $\sim$ 1021 cm-2 for the diffuse ISM) due to self-
shielding. In many cases the dissociation rate of H2 is not so high as to make
it completely unavoidable in cases where one is not studying in detail the
chemistry of the gas (Klessen & Glover, 2016).

Furthermore, given the inverse relationship between density and
photodissociation rate, combined with the timescale of destruction being on
the order of 1 Myr or higher for dense, optically thick gas (Abgrall et al.,
1992), we determine that the inclusion of photodissociation would not alter
the dense gas enough to significantly affect our results. The overall
timescale of formation of H2, in order to have enough present that its
destruction would be necessary, is on the order of 109 yr, with only an order
of magnitude decrease in supersonically turbulent environments (Klessen &
Glover, 2016). This is comparable to the crossing and dynamical timescales of
our simulations, both of which we do not model for any large fraction of their
scale.

Overall, self-shielding of H2 plays a large role in CNM clouds with no nearby
massive stars, such as we are modelling here, and Grackle does perform
estimations of self-shielding tied to H2 formation on dust grains, even
without radiative transfer or UV sources. We conclude that our molecular
hydrogen content acts sufficiently as an upper limit to the true abundance one
would have. To account for this being an upper limit, we henceforth discuss
molecular hydrogen through column density cuts of 1021 cm-2 and above. In the
case of volume density discussions, we have used the column density cut to
determine an average volume density of 100 cm-3 and higher for the ‘molecular
hydrogen’. Finally, we are tracking the formation of molecular hydrogen in
these simulations so we follow its rapid formation by starting with an
entirely atomic gas that converts to molecular quickly.

####  3.2.3 Clumpfinding & Sink Particles

At any time during the simulation run time, the gas can form both clumps and
sink particles. We define our clumps with a minimum density of 100$cm^{-3}$
and a minimum mass of 100$M_{\odot}$, such that they represent molecular gas
clumps and GMCs (with enough mass gained over time). These are formed through
the 3D clumpfinder algorithm in Ramses, outlined in Bleuler et al. (2015). As
a brief overview, this algorithm works by identifying density peaks in a data
set using a peak-based approach. Peaks with significant height compared to
their valleys, defined as topological relevance, are identified and labelled,
while nearby peaks with lower relevance are merged into the same clump as the
larger, given a certain saddle ratio between them. This algorithm allows us to
identify and label larger, extended structures that cannot be accurately
represented by a particle, and keep track of their properties. In our
simulations, this translates well to tracking the clouds in the complex we
form from our CNM, allowing for a matching between the clusters that form and
which clouds they form within.

The clusters form via sink particles, initialized at very low masses in order
to allow the sink’s mass to be a product of accretion using a density
threshold scheme, outlined in Bleuler & Teyssier (2014). Starting with very
low masses ensures that we follow the formation of the protoclusters as early
as possible, such that sink particle ages correspond very closely to the age
of the protocluster and we can accurately limit our runtimes to timescales
well before star formation and feedback would take place. Sink particles are a
useful catch-all particle that can easily represent any point-like density
peak, given the appropriate sub-grid physics is applied, so we have ample
freedom in our code to allow them to represent clusters. In general, the
agreed upon rules for sink formation are those outlined in Federrath et al.
(2010b). These rules outline conditions for density thresholds, refinement,
boundedness checks, gravitational potential minima, stability, and location.
Ramses adopts the same general conditions of sink particles, but simplifies
the rules. Sink particles are formed from the clumps found in clumpfinder,
which asserts they form within already dense structures, and we set our sink
formation density threshold to 5000$cm^{-3}$ to ensure we are only picking
cluster candidates out of the star-forming gas density peaks within a clump.

Once sink candidates are identified, they undergo 3 checks. First is the
virial check, in which the code checks if the gas that would be contained in
the sink is able to gravitationally collapse. More precisely, the virial check
analyzes the gravitational field at a possible site for sink formation,
ensuring it is universally compressive and strong enough to overcome internal
pressure support. Second, the collapse check verifies that the gradient of the
velocity is negative across all principle axes, ensuring that not only is the
gas collapsing, but it is contracting. The final check is the proximity check,
which ensures that gas which is being accreted by a sink particle cannot form
its own sink particle. Contrary to the outline of Federrath et al. (2010b),
Ramses does not carry out separate tests for boundedness, Jeans instability or
gravitational potential minima. We refer the reader to the original methods
paper for sink particles in Ramses for more detailed explanations (Bleuler &
Teyssier, 2014). We also note that discarding these tests allows for a more
general implementation of the existing tests. For example, while virial checks
give valuable information on gas collapse, the authors point out it is
possible for gas to exist in virial equilibrium but not be collapsing. They
then argue that their virial check approaches both collapse and boundedness by
being implemented more generally.

We split our simulation runs into two groups: those with and without magnetic
fields. Within these two groups, we further split it into three runs with
different turbulent velocity dispersions. Velocity dispersions are set via the
turbulent auto-correlation time (Federrath et al., 2010a) and the rms
acceleration (Schmidt et al., 2009), defined as $f_{RMS}=3v^{2}/L$. These
represent the parameters through which one can set an initial turbulent
velocity. The auto-correlation time sets the amount of time needed for one
wave or shock to cross half the length of the box, such that it is half of the
crossing time. The rms acceleration, on the other hand, sets the amplitude of
the turbulence by setting the acceleration of the material. Due to their
dependence on the rms velocity, the two must be set to determine the velocity
field within the simulation at startup. We set these to correspond to speeds
approximately double the desired dispersions, found through trial and error
for all other parameters of our box already set, such that we would reach the
desired dispersions at $10\%$ of the global crossing time, or 1 local crossing
time for a 100 pc patch. Actual dispersions achieved at initialization
($\sigma_{turb,initial}$) and at the onset of sink formation
($\sigma_{turb,sink}$) are given in Table 1. In the following section we
compare results from the different dispersion values as well as between hydro
and magnetic field runs.

##  4 Results

One of the most important physical attributes of any phase of the ISM is the
velocity dispersion that it supports. As is discussed in Klessen & Glover
(2016), there is a theoretical range of velocity dispersions (both in 1
dimension and 3) which can be present in the ISM, specifically depending on
the phase of the gas. From Table 4 of Klessen & Glover (2016), the velocity
dispersion of the gas decreases as we transition from atomic to molecular
hydrogen. Assuming a scale height of 0.5 kpc in the solar neighbourhood, the
HI gas can display velocity dispersions of $\sim 12km~{}s^{-1}$, whereas
molecular gas may have velocity dispersions $\sim 5km~{}s^{-1}$. Taking the
average, assuming completely equal portions of HI and molecular hydrogen gives
a 3D velocity dispersion of $8.5km~{}s^{-1}$. What is the origin of these
values? This question motivated our first set of numerical experiments.

###  4.1 Self consistent turbulent amplitudes in the CNM

One approach to understanding turbulent amplitudes is to assume that the ISM
achieves a hydrostatic balance between turbulence gas pressure resulting from
stellar feedback, and the weight of the gas (Kim & Ostriker, 2015). Too little
turbulence and gravitational collapse of the medium ensues - while too much
will heat and disperse it. The velocity dispersion predicted by this
hydrostatic balance condition is

$\sigma_{z}^{2}=\frac{\pi}{2}\frac{G\Sigma^{2}}{\rho}$ (5)

where $\sigma_{z}$ is the vertical component of the velocity dispersion,
$\Sigma$ is the average column density of the gas, and $\rho$ its average
volume density.

We first chose initial velocity dispersions near the predicted values of
Klessen & Glover (2016). We ran a large number of simulations for a box
without periodic boundary conditions. For this setup high velocity dispersion
strongly shock heats the gas and disperses gas out of the box, while too low a
velocity dispersion leads to general collapse. After some trial and error, we
find minimum 3D velocity dispersions of $\sim 5km~{}s^{-1}$, and a maximum of
$\sim 10\rm{km~{}s^{-1}}$, in order to create a filamentary ISM, similar to
Joung et al. (2009).

With these experiments completed, we then used these initial conditions in
periodic box simulations, which we show in all subsequent figures.
Additionally, to consider the formation of dense filaments and clusters, we
further limit that range between 8 and 10 $km~{}s^{-1}$, due to the fact that
our 5 km/s run experiences a very slow global collapse. This makes structure
transient, as the turbulent shocks are not strong enough to maintain sharp
filamentary structure. The formation of filaments and subsequent clusters is
therefore happening on much longer timescales. We find that the 5 km/s
simulation fails before reaching its burst of sink formation due to
overcooling of the completely empty edges as the entire simulation globally
collapses. For these reasons, we do not take the 5 km/s models into account
for any of our conclusions surrounding sink particles. We still include said
models in the discussion for completeness.

Our numerical results are also reflected in observational work such as HISA
observations of GMF38.1-32.4a (Wang et al., 2020). These cite velocity
dispersions in HI of $\sim 4-6km~{}s^{-1}$, with some regions reaching in 13
CO as high as 10$km~{}s^{-1}$. Averaging between HI and 13 CO gives dispersion
ranges of 6-12$km~{}s^{-1}$.

Figure 2: Column density projections of our hydro runs once the first sink
particle has formed. Top: Model for our 5 km/s dispersion case at 6.1 Myr.
Projections in both the x and z plane are shown, with the center of our box
centred at (0,0). Middle: Column density projection of the 8 km/s dispersion
case at a time of 3.7 Myr. Bottom: Projections for our 10 km/s dispersion case
at 2.8 Myr. Sink particles are shown as black filled circles.

Figure 2 shows the column density plots in both the x and z plane for each of
our hydrodynamic runs at the beginning of cluster formation. While the plots
of the 8 km/s and 10 km/s runs give very similar structure, the properties of
the structure are different. With a dispersion of 8 km/s, the structure is
sharper, due to the self-gravity of the filaments being strong enough to pull
them together, yet not so strong as to dominate completely and collapse all
the gas. In our 10 km/s run, the structure becomes more diffuse looking. In
our 8 km/s models, the ratio of kinetic to gravitational energy is 0.01,
indicating turbulent motions play a dynamic role in the formation of the
structure. On the other hand, as this global ratio reaches 0.03 for our 10
km/s models, we see turbulence overcome gravity on scales of $\sim 100pc$,
creating more transient structures. This higher proportion of turbulent energy
to gravitational potential energy contributes too much mixing throughout the
box, creating structure that cannot be pulled together by gravity as much as
in the 8 km/s case. A dispersion of 5 km/s, on the other hand, displays the
opposite. In this case, the gravity is far stronger than the turbulence can
be, not allowing it to create nearly as much structure. The gravitational
potential dominates the gas, causing a global collapse towards the centre of
the box, still creating structure, though with far shallower density
fluctuations than our more turbulent setups. Based on structure alone, it is
clear that a dispersion of 8 km/s can create a well balanced and realistic ISM
slice.

Figure 3: Slices through the center of the box of our disp8 model with density
(left) and temperature (right) plotted at 3.7 Myr, just before sink formation
begins. Sharp filamentary features are seen in the density, mirrored in the
temperature slice, where hot gas outlines the edges of filaments. This slice
has an average density of $\sim 10^{2}cm^{-3}$ and an average temperature of
$\sim 10K$, lending an average sound speed of $0.3kms^{-1}$ throughout the
simulation.

If one compares each of these cases to Figure 1, it can also be seen that the
8 km/s case creates structure much more akin to observations of Orion A,
specifically the structure visible in the 250$\mu m$ dust emission. Figure 3
shows densities and temperatures in a slice of the 8 km/s model, in which we
can see the densest filaments rest in temperatures $\sim 10K$. We see higher
temperatures ($\leq 600K$) at the boundaries of the filaments, from shock
waves creating higher temperature shock fronts. Our filament temperatures
compare well with observations of the Orion A GMC, which has temperature
ranges of $\sim 10-30K$ (Lombardi et al., 2014), while the hot, diffuse gas
surrounding the filaments would likely be too diffuse to be observed. We note
some cold gas of temperatures $\sim 1K$, tending to appear in low density
pockets of the gas. Given the placement of these cold pockets, we conclude
that their presence is due to shock waves emptying the area of gas, pushing it
into the hot boundaries of the filaments. The average temperature of the
entire simulation is $10K$. The presence of a large range of temperatures
(spanning $1-1\times 10^{3}K$), emphasizes the importance of cooling processes
that give rise to a multiphase ISM in simulations of filamentary structure,
most especially at the current scale. Additionally, in Fig. 2, we notice a
large structure just below the center on the 8 km/s density projections. The
size and components of gas present in this structure are similar to Orion A,
which is 90 pc in length and has approximately 45% of its gas sitting at high
column densities (Großschedl et al., 2018). We further discuss comparisons
between our models and observations of Orion A in the following Sec. 4.2.

Finally, we note that using our quoted values for $\Sigma$ and $\rho$ in
Equation 5, we find $\sigma_{z}\simeq 30\rm{km\>s^{-1}}$. This is at least a
factor of 3 higher than our own simulations suggest, yet matches well with
extragalactic observations such as Wilson et al. (2019) who find dispersions
of 30-80 km/s for unresolved galaxies. Given $T\propto\sigma_{tot}^{2}$, a
velocity dispersion of 30 km/s gives temperatures of more than $10^{4}$ K,
representing highly shocked gas that will quickly destroy the CNM. In fact,
our own simulations do see some of the shocked gas create high temperatures
over 1000 K in Fig. 3.

However, we note that all the simulations of Kim & Ostriker (2015); Kim et al.
(2013, 2011) also have velocity dispersions (see Figure 2 of their 2015 paper)
of 5-10 km/s, in excellent agreement with our own values. They use setups
typical of the local ISM with gas surface densities of $\approx
1-10~{}\rm{M}_{\odot}/\rm{pc}^{2}$ and an effective scale height of $\sim 60$
pc on average (Kim & Ostriker, 2015). However, they also report that these
velocity dispersions seem to be somewhat independent of the surface density
used; increasing $\Sigma$ by a factor of 10 (see Figure 11 in Kim et al.,
2013) made no difference to their resulting velocity dispersion. Our
simulations, sitting a factor of 10 higher in column densities than theirs
also agree with this finding. Thus, Equation 5 appears to be only a rough
guide to the link between velocity dispersion and gas column density with
large deviations possible. Taken together, all of these results point towards
a more general idea that velocity dispersions in the ISM are a consequence of
the gas phases themselves.

###  4.2 Molecular gas and cloud morphology

In our models we define our molecular gas at densities of 100 $cm^{-3}$ or
higher in order to best compare to dust emission maps of molecular clouds,
which contain no chemical information about the molecular gas. Using our 8
km/s models, we can consider the masses of gas at densities of 100, 1000 and
10000 $cm^{-3}$ to determine amount of molecular and star forming gas in our
model. Table 2 shows these masses. At the time of formation of the first sink
particle, we see $\sim 50\%$ of our gas mass lies at densities of 100
$cm^{-3}$ or higher, showing that a significant fraction of our gas is tied up
in dense structures that will eventually form stars. Our highest density cut,
representing gas which has surpassed star-forming densities and, therefore,
the threshold for our sink particle formation, contains only $0.2\%$ of our
total gas mass, showing a significant jump from molecular gas densities to
star-forming densities. From this, we can see a considerable amount of gas in
the CNM will become molecular gas, as we can also see in Fig. 4, but very
little of that will continue to climb in density to eventually form stars.

| Hydro | Magnetic Fields, B=7$\mu$G
---|---|---
$M_{T}(M_{\odot})$ | $9.3\times 10^{7}$ | $9.3\times 10^{7}$
$M_{100}(M_{\odot})$ | $4.9\times 10^{7}$ | $4.7\times 10^{7}$
$M_{1k}(M_{\odot})$ | $1.2\times 10^{7}$ | $1.4\times 10^{7}$
$M_{10k+sink}(M_{\odot})$ | $1.9\times 10^{5}$ | $4.7\times 10^{5}$
Table 2: Comparison of gas masses between the hydro and magnetic field cases
for our disp8 model. From top to bottom we give total gas mass, and total
masses for gas above densities of 100, 1000, and 10000 $cm^{-3}$,
respectively. For density cuts of 10000 $cm^{-3}$ we include total mass in
sink particles as well.

We also note a large complex of molecular clouds just below the center of the
box, extending roughly 100 pc in length and a maximum of 30 pc in breadth. We
find this size is comparable to that of the Orion A GMC in Fig. 1. The
similarities between the two also extend into the contours tracing molecular
gas. In both figures we see the majority of the molecular cloud complex
outlined by one contour, with the resolution of the contours unable to depict
any of the fine filamentary structure creating the complex. Additionally, this
is where we see many of the sink particles form, indicating it is the primary
area of star-forming gas. We emphasize that while there are similarities with
the Orion A cloud properties, the environment of our simulated cloud is very
different. Our point is just that the formation of GMCs with recognizable
properties is natural in a turbulent CNM setting.

Figure 4: Side-by-side comparison of the molecular gas structure for the disp8
model. Left: Z-plane projection of gas, cut to contain only that above
$10^{-21}g\cdot cm^{-3}$. Maximum resolution is 0.48 pc, chosen to approach
typical filament widths of 0.1 pc. Right: Same, with added contours to
simulate observations of molecular gas structures. Contours are drawn from
$10^{20}cm^{-3}$ to $10^{23}cm^{-3}$ to mimic resolution range of Lombardi et
al. (2014). Sink particles are denoted by filled black circles.

We provide a zoomed in look at the molecular cloud complex noted above in Fig.
5. While the contours break up into smaller clumps, they still closely trace
the filamentary structure and outline a cloud sized complex of molecular gas.
For the purposes of comparison, the Orion A GMC complex is $\sim 7.57\times
10^{4}M_{\odot}$ (Lombardi et al., 2011), whereas the mass of the molecular
gas contained in this cut out is $\sim 5.78\times 10^{4}M_{\odot}$. Our entire
zoom region is larger (150 pc) compared to the size of Orion A, which is $\sim
90$ pc (Großschedl et al., 2018). Additionally, we find we have a lower
percentage of high density gas (25%, compared to 45% in Orion A according to
Großschedl et al. (2018)), which can also be attributed to the young age of
our simulation. Our average velocity in the area is $13.9$ km/s, higher than
the assumed 10 km/s from Wilson et al. (1970), though not entirely outside an
acceptable range. The velocity dispersion in this region is $7.5$ km/s, more
than double the estimates of $\sim 2.5$ km/s for Orion A(Theissen et al.,
2022).

Overall, this molecular cloud candidate can be compared to observations of
molecular clouds, though it does not fully represent the Orion A cloud
complex. We attribute the differences to the short timescales we run our
simulations for, and intend to investigate the evolution of this area at later
times in future work. We find our simulations produce clouds similar to
observed molecular cloud structure and properties though we hasten to add that
a more detailed comparison requires a statistical sample of simulated clouds,
which is beyond our computational capability in this paper.

Figure 5: Zoom-in of contoured region in Fig. 4. Density projection and
contours are the same, though plotting has been recentered on the center of
the location and the domain plotted is 150 pc. Sink particles are represented
by filled black circles.

Fig. 6 shows the density and temperature PDFs of our 8 km/s dispersion model.
We use dotted vertical lines to show the average density of all the gas. We
also do a temperature cut of all the cold gas based on where our temperature
PDF flattens out. The dashed vertical lines in our density PDF mark the
average density of the colder gas. Temperature PDFs show an order of magnitude
difference in the peak temperatures of our models, corresponding with the
results discussed for Fig.10 in Section 4.3 below. Discussion of the
temperature differences follows there.

While the density PDFs between the two models are quite similar, we notice a
second peak beginning to form in our magnetic model. As this also corresponds
to a shift in our average density, we can relate this to the cold gas, such
that we are producing a double log-normal, with each peak corresponding to a
region in the temperature PDF (roughly split here between cold and warm gas).
We note that our peaks are quite subtle as we are limiting ourselves to very
early times in the gas, whereas the work of Joung et al. (2009) show much more
distinct peaks. However, we still see similar general trends between our work
and theirs, expecting that later stages of our simulation would evolve quite
similarly.

In the high density range of the PDF we notice a deviation from a normal
shape, indicating a power-law tail produced via self-gravity as explored in
works such as Kainulainen et al. (2014). The density PDF of our MHD model is
also wider than our hydro model, possibly due to the correlation between gas
density and temperature, as higher temperatures would be found in more diffuse
gas, and the high-density shocked gas.

Finally, the bottom panel of Fig. 6 shows that the temperature of the peak in
the density PDF in the hydro simulation is lower by nearly an order of
magnitude than is the case for the MHD run (a few degrees K in the former
versus 30 K in the latter case). This is an important consequence of the fact
that magnetic fields prevent the high gas compression that can occur in purely
hydrodynamic shocks. Lower densities imply lower cooling rates so more of the
gas remains warmer than in the pure hydro simulation. This is further
emphasized in our disp5 models shown in Fig. 16, which show that the exclusion
of magnetic fields allows the gas to overcool due to the global collapse of
the model. This brings a non-negligible amount of gas below 3 K and
effectively stops our run.

Figure 6: Density and temperature PDFs for the 8km/s models. Dotted vertical
lines on the density PDF mark the average density of all gas, whereas dashed
lines represent the average of cold gas in the simulation. The hydro model
shows a cooler temperature peak by an order of magnitude compared to the MHD
case, significant of the higher compression of the gas in these models.

###  4.3 Hydro vs. magnetic fields

To begin our analysis of the effects of magnetic fields on structure formation
in the CNM, we first look at its effects on the structure created in different
scenarios of turbulent strength. In figure 7, we show the orientation of the
magnetic field to the filament. We use Magnetar (Soler et al., 2013), which
analyzes the relative orientation of the magnetic field and the filament using
density gradients to determine the orientation of the filament. The cosine of
the angle between filament and magnetic field will be normally distributed for
a gas with parallel magnetic field alignment. A flat distribution signifies
random orientations with no net alignment. In Fig.7, we plot the $\xi$
parameter, which describes the concavity of the histogram of $cos(\phi)$. When
this value is positive, the magnetic field is aligned parallel to a filament
while for negative values, the magnetic field is aligned perpendicular. For
each of our MHD simulations, we plot this parameter against the density for
two times: one at which sink particle formation has just begun, and one
halfway between the start of the simulation and the onset of sink formation.
In this way, we are able to analyze not only the orientation of the magnetic
field, but the evolution of the orientation as well.

We begin by noting that the early time frame has a smaller density range as
the gas in our box is still evolving at this point. It is still highly
turbulent and has not yet formed any stable filaments, and therefore has not
yet reached high densities. Even with this smaller density range though, it is
still clear that there are significant differences in the orientation of the
field. In all three early time cases, magnetic fields are mostly parallel to
filamentary structure, or showing no real alignment, as signified by the line
hovering around zero in low densities. However, both the 8 and 10 km/s cases
show small peaks into perpendicular alignment at molecular gas densities, with
magnitudes of 0.2. These peaks of perpendicular alignment could be the
beginning of molecular filament formation in the simulations, indicating a
link between the slightly larger range of densities in these two compared to
the 5 km/s case and their perpendicular field alignments.

At the point of sink formation, we see a more evolved gas, as evidenced by the
larger range of densities. The magnetic field orientations have changed over
this time scale, showing higher magnitude concavity, indicating a stronger
alignment parallel to the filaments. For all three models, the magnetic field
begins to have a perpendicular alignment at densities of
$10^{2}~{}\rm{cm}^{-3}$, corresponding to our molecular gas filaments. Beyond
this point, the models’ similarities end.

Dispersions of 5 km/s show a weaker perpendicular alignment to filaments. At
densities of $10^{3}~{}\rm{cm}^{-3}$, the alignment of the magnetic field
switches to parallel. In the highest densities, the concavity approaches zero,
possibly showing a random alignment in the innermost regions of a filament and
around clusters. Dispersions of 8 and 10 km/s show stronger perpendicular
alignment. Once perpendicularly aligned, the 8 km/s model gradually climbs
back to random and even weakly parallel alignments, becoming parallel at
densities of approximately $3000~{}\rm{cm}^{-3}$. While the change from
perpendicular to parallel alignment is smoother in the 8 km/s model, the
alignment still definitively changes at the innermost high density portions of
filaments. Our 10 km/s model does not display any parallel alignment in the
gas though it is trending towards parallel alignment at high densities.
Instead, the gas remains at least weakly perpendicularly aligned throughout
the entire density range, indicating a higher density needed in filament
centres in order to sweep magnetic fields into parallel alignment.

As we reach the highest density filaments, we see that fields parallel to the
filament occur at the time that sink particles appear. This is the signal that
gravity now is stronger than magnetic field tension of the perpendicular
field. The gravity driven flow towards the clusters (i.e., sink particles)
bends the field into a parallel configuration (Klassen et al., 2017). This
alignment aids in funnelling material into dense clumps along the filaments.
We see evidence of this phenomenon in all three of our models, but notably the
density at which this occurs is not constant. Instead, the density necessary
to pull magnetic fields into a parallel alignment is connected to the amount
of turbulent support in the gas, where high velocity dispersions reach higher
densities in the centres of filaments before the magnetic field aligns
parallel with gas flows within the filament.

Figure 7: Magnetic field orientations relative to filaments shown through
$\xi$ \- the concavity parameter of the histogram of relative orientations
from Soler et al. (2013) for each of our MHD models. From top to bottom: 5
km/s dispersion, 8 km/s then 10 km/s. A $\xi$ value greater than 0 indicates
parallel alignment of magnetic field with the filament axes, whereas less than
0 indicates a perpendicular alignment. Evaluated for times of the onset of
sink formation (orange lines) and halfway between simulation start and sink
formation time (blue lines).

Figure 8 presents the structure of the magnetic field via its magnitude
through a slice of the simulation box. These maps reveal where our highest
magnetic field strengths lay. We can see strong magnetic field lines tracing
the high-density filaments, displaying a link between field strength and
density in our gas. However, though all three cases show the same link between
gas density and field strength, we see differences in the filaments traced in
each model.

In our least turbulent model, the 5 km/s dispersion case, the filaments of
strong magnetic field are thin, sharply defined and quite long. Many of the
filaments trail the length of the entire simulation domain, and shorter
filaments seem to quickly diffuse out into lower field strength (and lower
density) gas. We see some bending in the field around the area we see our
highest density structure forming, but the majority of the structure is
smooth.

In our 8 km/s dispersion model, we see more intense field variations. While
our filaments can still span the length of the domain, there are fewer short
filaments overall, and those that are present do not diffuse in the background
field as quickly. Additionally, we see that though our filaments are wider, we
have more structure variation throughout. Notably, we can distinguish the
location of our highest density structure from the complication of the field
lines around the area. As an example, in our z-plane view we know the
structure to exist just below the center point of the box, spanning across the
x direction. Likewise, in the slide of our magnetic field, we see a very dense
filament traced by highly magnetic gas breaking off into smaller filaments and
loops just below the center of the domain. We can as well see many filaments
turn around completely causing sharp u-bends and overlapping with other
filaments.

Finally, in our highest turbulence case, the 10 km/s dispersion model shows
many of the same features. Magnetic field strengths increase with gas density,
effectively showing strong magnetic fields in dense filament centres.
Filamentary areas where we see cluster formation in our gas also contain the
strongest field strengths along the direction of the filaments. The difference
between this and the 8 km/s model largely comes in through the filament
widths. We can see our more dynamic case lends wider and more diffuse looking
filaments than our 8 km/s case, not unlike our discussion of the densities,
where we find the 10 km/s models to be more diffusive due to the stronger
turbulence.

Figure 8: Magnetic field magnitude maps for each of our MHD models. From left
to right: x, z and y plane slices through the center. Sink particles are
represented by white filled circles.

Magnetic fields also delay structure formation (Banerjee et al., 2009). In our
simulations, we can see some ties between magnetic field strength, turbulent
velocity dispersion and the onset of structure formation. The time at which
the first sink particle forms is also linked to structure formation
timescales. Typical structure formation timescales of a turbulent medium
correspond to the crossing time of the gas, the time in which it takes for a
wave to cross the length of the domain. On the scales of structure necessary
for cluster formation, we consider the crossing time across the typical size
of a cloud (100 pc), corresponding to $10\%$ of the initial auto-correlation
time of the turbulence. However, we find that structure formation begins
rapidly in our simulations, and the formation of sink particles as clusters
can equally represent the fragmentation of structures. When comparing these
two, we find that each of our simulations form filaments well before they
reach $10\%$ of their crossing time, as can be seen Table 3.

|  $t_{cross,10\%}$ (Myr) |  $t_{sink}$ (Myr)
---|---|---
disp5 | 9.7 | 5.76
disp5MHD | 9.7 | 6.12
disp8 | 6.1 | 3.21
disp8MHD | 6.1 | 3.71
disp10 | 4.8 |  $\geq$ 2.7
disp10MHD | 4.8 | 2.77
Table 3: Comparison between time of first sink formation ($t_{sink}$) and
$10\%$ of the global crossing time ($t_{cross,10\%}$) for each of the six
simulations we run. Note the disp10 model, which does not form any sink
particles before resource limitations force us to end the run.

We see some distinct differences in timescales with the inclusion of magnetic
fields. While all of the models form sink particles well before reaching one
crossing time, the time of sink particle formation can vary. In the 5 km/s and
8 km/s models, we see sink particles form roughly 0.4 and 0.5 Myr sooner,
respectively, in the absence of magnetic fields. As such, there is an obvious
delay in cluster formation with the inclusion of magnetic fields demonstrating
the important dynamic role that fields can play in the ISM. As magnetic fields
support the gas and slow the formation of structure, so too will it delay the
onset of star formation.

However, as an outlier, we note that our 10 km/s model does not display this
same behaviour. In fact, while the hydro and magnetic cases are run for the
same amount of time, we see no sink formation start in our hydro model whereas
the inclusion of magnetic fields begins sink particle formation just before
the end of the run.

Figure 9: Mass-weighted pressure-density phase plot of the disp8 (top) and
disp8MHD (bottom) model, for densities above $1cm^{-3}$. The majority of the
mass is found in median pressures, molecular to star-forming gas.

Figure 10: Mass-weighted temperature-density relation for disp8 (top) and
disp8MHD (bottom) models, for a temperature range between $1-10^{3}K$ and
densities of $1-5\times 10^{5}K$.

The gas densities present in our simulation also vary between hydro and
magnetic cases. While we have discussed the breakdown of the gas in our hydro
models in Section 4.2 with Table 2, we also note the breakdown when including
magnetic fields. While we have the same total gas mass, we can see the
percentage of gas in each density cut compared to the total drops with each
higher cutoff. In our disp8MHD model, $50\%$ of the gas is molecular densities
of 100 $cm^{-3}$ or higher, when we increase to 1000 $cm^{-3}$ for our cutoff,
this drops to $15\%$ of the total gas mass, showing that most of the gas is
remaining at molecular gas densities or lower. When we consider the star-
forming gas, with density cutoffs of 10,000 $cm^{-3}$ or higher, only $0.49\%$
of the gas in the box reaches this level. While our magnetic field cases can
create an abundance of molecular gas, very little of that gas reaches star-
forming densities, leaving few clusters to form. Specifically comparing to the
hydro case, we see that the fractions vary slightly between the two. While
cuts for densities of $100cm^{-3}$ and $1000cm^{-3}$ contain similar gas
masses, we note a large difference in the star-forming gas cuts.

We have, therefore, discovered an important distinction between these two
cases. Even though the two models have similar amounts of molecular gas, our
hydro case has roughly half the amount of gas at star-forming densities,
indicating a lower number of possible sink particles that can form. This lower
fraction of high density gas would suggest that the inclusion of magnetic
fields does not significantly affect the cold gas, but can contribute to a
higher star formation rate as the magnetic fields allow more gas to funnel to
higher densities. Indeed, we also see this in the number of sink particles
each model forms, where our hydro cases formed only a couple cluster sinks but
our MHD models were able to form anywhere from 5-10 cluster sinks in less than
0.5 Myr.

###  4.4 Phase diagrams

Comparing the mass-weighted pressure-density relations in Figure 9, we note
the changes in slope of the relation for the highest masses. In particular, we
notice the hydro case contains a majority of its mass in molecular gas at
pressures of $\sim 10^{3}Kcm^{-3}$, up to a maximum of $\sim 10^{5}Kcm^{-3}$
with a minimum of $\sim 10^{1}Kcm^{-3}$. In the MHD models, on the other hand,
the gas sits in higher pressures by about an order of magnitude, with most
sitting at pressures of $\sim 10^{4}Kcm^{-3}$. Thus, even considering what
pressure at which the majority of the gas measures, the magnetic fields
contribute to overall higher pressure, a result which has also been found in
previous works such as Fiege & Pudritz (2000).

In Figure 10, we plot the temperature-density phase diagram of our two 8 km/s
dispersion models. Both have near constant slopes, indicating an isothermal
gas in both, though the average temperature is different by an order of
magnitude between the two. While the majority of the gas in our hydro model
sits around $10^{1}K$ cooler than our initial temperature, the MHD models are
supported at the initial temperature of our gas just by the inclusion of
magnetic fields. This net cooling present in the hydro models is likely due to
higher compression of the gas. Without magnetic fields, the gas experiences
less support against shock waves and compression, ultimately allowing for more
compressed gas than if magnetic fields are present. This compressed gas cools
faster, and the average temperature of our simulations will decrease as the
amount of highly compressed gas increases.

We see evidence of the increased compression in our hydro models in Figures 9
and 10. Note that high densities have their pressure and temperature spreads
significantly thinned, and we see fairly sharp relation with little to no
spread in both temperature and pressure for our hydro models. This indicates
the compression of gas, as high densities are pushed to a stronger pressure or
temperature relation. On the other hand, our MHD models do not display such a
sharp relation, especially so in our temperatures, where high densities follow
our isothermal line and are free to have larger spread in temperature. The
overall temperature trends, in tandem with our results from Table 2, indicates
that the presence of magnetic fields could be necessary to allow for the
creation of hot, gas as well as setting the star formation rate in the ISM.

These temperature spreads also indicate the presence of multiple phases. While
the majority of our gas exists in a CNM phase, we can also see a spread
towards warmer temperatures at similarly low densities, indicating a warm,
neutral phase of hydrogen present. This multiphase ISM structure -
specifically the heating to warm and even hot temperatures - comes about
purely from shock waves. However, we do note that an ionized medium phase is
missing because we lack supernova feedback in our simulations. In order to
push to these dense, hot regions of the temperature-density space it is likely
necessary to add those ISM heating mechanisms we know to be important at kpc
scales like these (i.e. photoelectric heating and UV radiation, as discussed
in Draine (2011); Klessen & Glover (2016)).

###  4.5 When do clusters form?

####  4.5.1 Timescales of formation

One can also consider the timescale of cluster formation, as well as how many
clusters form. In Figure 11, we show the mass function of the sink particles
that form across our disp8MHD and disp10MHD models. Disp5MHD is excluded as it
only forms one sink particle during the simulation. Mass growth of clusters is
rapid early on, quickly gaining more than 7 orders of magnitude of mass in
less than 0.2 Myr. The mass growth starts to exhibit a turnover, slowing down
significantly after the initial burst of growth, consistent with results from
Howard et al. (2018). Both the 8 and 10 km/s dispersions display similar mass
growth, indicating this pattern is not linked to the dynamics of the gas.

We additionally compare the formation time for sink particles across all our
simulations. From Table 3, we can see that sink formation starts well before
one crossing time has passed on a cloud scale and all of our runs begin
forming young cluster sinks at least 2 Myr before one local crossing time has
passed.

Figure 11: Mass growth of sink particles over time shown for MHD models disp8MHD and disp10MHD. | Lowest Mass ($M_{\odot}$) | Highest Mass ($M_{\odot}$)
---|---|---
disp5MHD | $4.92\times 10^{2}$ | —
disp8MHD | $1.32\times 10^{2}$ | $2.94\times 10^{3}$
disp10MHD | $2.28\times 10^{1}$ | $1.36\times 10^{3}$
Table 4: Lowest and highest mass sink particles for each of our MHD models. In
the case of our 5 km/s dispersion, we note only one sink particle forms and
therefore no upper bound is available for masses.

Both the 8 and 10 km/s simulations show rapid sink formation before reaching 4
Myr. In our non-magnetic case, the first sink to form for the 8 and 10 km/s
cases form at times of 3.33 and $\geq 2.7$ Myr, respectively. The two have
similar bursts of star formation, displaying a rapid onset of formation
immediately. The magnetic cases display remarkably similar behaviour. While
onset of sink formation occurs later, 3.71 and 2.77 Myr for the 8 and 10 km/s
cases respectively, the two cases still display early bursts in star
formation, forming many sinks within 0.1 Myr. Both these time frames are
consistent with findings from Li et al. (2017).

For our 5 km/s simulations, sink formation occurs significantly later and many
fewer sink particles form. In our non-magnetic case, sink formation begins at
5.76 Myr, and our magnetic case begins sink formation at 6.12 Myr. While this
is significantly later than the 8 and 10 km/s models, these times are still
well ahead of the local crossing time, unanimously indicating cluster
formation begins well before one crossing time. Additionally, both 5 km/s
cases produce very few sink particles, indicating we have not reached the
burst of formation for these models yet. Computational limitations prevent us
from running these models longer.

In the case of disp5MHD, we see only one sink form at 6.12 Myr, while the
simulation is cut off at 6.2 Myr. Similarly in the hydro case, our disp5 model
forms 1 cluster sink before the simulation is ended. With sink particle
formation beginning so late in the run, and the formation rate being
significantly slower than either the 8 or 10 km/s cases, the cloud is close to
the end of its lifetime before forming any massive clusters and may never form
more than a few. Further evolution of the cloud is needed to investigate the
later cluster formation for a burst as is seen in our 8 & 10 km/s simulations.
Additionally, though the gas is self-gravitating, the density fluctuations
formed via shock waves in the medium are not strong enough to create the extra
compression needed for supercritical filaments and their associated star
clusters. For these reasons, we conclude that a velocity dispersion of 5 km/s
is too low to create the extensive filamentary structure observed in the CNM.
We conclude that 8 km/s is ideal for both cluster formation onset and
structure formation.

####  4.5.2 Cluster accretion

Another important factor that governs the time scale for cluster formation is
the accretion rate onto the region. In Fig.12, the accretion rates of each
individual sink particle in our MHD simulations are plotted against the age of
the sink particle, for better sink-to-sink comparison. We can see that when
sinks are just forming, their accretion rates can range from
$10^{-6}-10^{-3}M_{\odot}/yr$, depending on the strength of turbulence and the
location of the sink particle. For example, our highest number of sink
particles formed occurs in our disp10 simulation, where turbulence is very
strong and the filamentary structure forms quickly. The sink particles in this
case (represented by the purple lines) show a spread in accretion rates, and
each sink’s rate stays steady, peaking upwards at later ages, when the sink is
more massive and more active. While the average accretion is low, we note that
these cluster sinks are still very young and we expect their accretion rates
to increase rapidly as time goes on.

Other published simulations of cluster formation show accretion rates of at
least $10^{-4}M_{\odot}/yr$ and up to $10^{2}M_{\odot}/yr$ (see, for example,
Reina-Campos et al., 2022; Bieri et al., 2022; Howard et al., 2018).
Furthermore, the average mass of a cloud in our simulations is $1.3\times
10^{5}M_{\odot}$, comparable to the $7.5\times 10^{4}M_{\odot}$ of Orion A,
but lower than average cloud masses from galactic simulations of GMC formation
(Grisdale, 2021). Given that our clouds are lower mass, the scaling relation
between cloud and cluster mass discovered in Howard et al. (2018) shows that
cluster masses would also be in the low mass end of cluster mass ranges,
particularly in the time frame for which we run our simulations.

We also note that over longer periods of time we would expect a varying
accretion rate, and the smoothness of the rates presented here is due to the
relatively short timescale on which we are looking at these clusters. An
analysis of these clusters at later times will follow in future work. No
mergers take place between sink particles to grow their mass, though the early
phases of a cluster’s growth will be dominated by mergers (Guszejnov et al.,
2022). Instead, we see mass growth only through initial formation and gas
accretion as the timescales discussed here are still shorter for the clusters
than discussed in the paper referenced. This result is in agreement with
simulations of massive isolated GMCs wherein only the most massive clusters
undergo significant mergers with smaller clusters (Howard et al., 2018). On
average, our sink particles reach masses of $10^{3}M_{\odot}$ before our
simulations end - i.e. an Orion Nebula Cluster type of mass scale. Our 8 km/s
dispersion model is responsible for the most massive cluster sink, reaching a
mass of $2.94\times 10^{3}M_{\odot}$, whereas its lowest mass cluster is
roughly $10^{2}M_{\odot}$. Lowest and highest masses of each of our MHD models
can be seen in Table 4.

Figure 12: Accretion rate vs. age of the sink particle for our three MHD
cases: disp5MHD, disp8MHD and disp10MHD. Age represents how long the cluster
sink has been around, not the physical time in the simulation.

###  4.6 Angular Momentum and Filament Stability

In this subsection, we analyze two other key factors governing structure
formation; angular momentum and gravitational stability. On the first issue,
we note that the local specific angular momentum distribution inside our
simulations is a consequence of both the initial conditions - the power
injected in solenoidal modes - as well as the subsequent interaction of
oblique shocks.

In Fig. 13, we plot each individual component of the specific angular
momentum, as well as the total magnitude, as a function of radius from the
center of the box. Most notably, one can notice the shallow, positive slope of
the total magnitude, indicating a lack of net rotation in the box. We further
support this conclusion with the flux of specific angular momentum, shown by
the dashed line. As the flux per radius is constant and low, we conclude that
global rotation is not present in the simulation. Instead, any local
rotational motion is essentially random across the box, brought on by the
interaction of shock waves moving through the domain (Jappsen & Klessen,
2004). The directional components of the specific angular momentum show
relatively even levels with each throughout most of the box, indicating they
are balanced and no dominant axis for rotation stands out. However, the outer
edge does show a dominance from the x component, which may indicate some
rotation along the outer edge of the box, but is not dynamically important as
we do not see any other evidence of the possibility of rotation.

We conclude that turbulence does not create any net rotation in our simulation
domain, which is a good check on the veracity of the RAMSES code since no net
initial rotation was put into our simulation box. Turbulence does produce
local fluctuations in angular momentum which will carry over into the
formation of protostellar disks on the smallest scales (not probed in our
simulations). We also note that the effects of galactic shear on our box could
be important, and this we investigate in Zhao et al. (2023).

Figure 13: Specific angular momentum profile of gas throughout the box, moving
radially outwards from the center, for our disp8MHD model. Both the individual
components and total magnitude of the gas’s specific angular momentum are
plotted.

Figure 14: Line mass ratio maps for projections of top: our disp8 model and
bottom: our disp8MHD model. Black X’s show the positions of sink particles.
Projections are taken for the x and z plane, as they showcase the most
structure. Critical line mass is calculated as $M_{crit}=2c_{s}^{2}/G$,
thereby indicating only thermal stability. Line mass is taken as cell
mass/cell width and sampled for each cell, therefore providing a lower bound.

We consider now the stability of the filamentary structure created in our 8
km/s dispersion models. Figure 14 depicts projections of the line mass ratio
of the gas in our models. We note that our simulations have considerable
turbulent amplitudes on the supra molecular scales, and these are much larger
than found within molecular clouds (recall the turbulent scaling relations;
$\sigma_{turb}\propto L^{1/2}$). André et al. (2010) considered nearly thermal
filaments in low mass molecular clouds. At the same time, the CNM has much
higher temperatures than the cold molecular clouds so that thermal speeds and
critical line masses in CNM filaments are higher. In order to make comparisons
then, we use only the thermal motions for our critical line mass as they do.
We do not include any turbulent motion or magnetic field corrections. As such,
our ratios in Figure 14 represent lower bounds on the true line mass ratios of
the gas.

In both the hydrodynamical and MHD models, we see largely subcritical thermal
line masses for gas throughout the domain, though the MHD model contains gas
that is more subcritical than our hydro model. Without magnetic fields, the
gas is allowed to approach much higher line mass ratios, the majority of the
gas sitting around a line mass ratio of 0.1, corresponding to transcritical
line masses. One can also notice that the structures in line mass projections
trace the structures we see in Fig. 2. Our densest filaments also contain the
gas closest to supercritical ratios. Specifically the region between z=-0.1
and z=0.0 in the top left plot of Fig. 14 contains gas approaching the
critical ratio, as indicated by the faint green colour in the filamentary
structure. Empty pockets of gas which do not make up the filamentary structure
are shown in dark pink, with very low ratios. In our hydrodynamical models,
there are very few of these areas, indicating most of this gas is still close
to collapse and that the structure is either overall more unstable or younger
and therefore has had less time to settle.

The bottom row of Fig. 14 shows the same projections for our disp8MHD model,
to compare the effects from magnetic fields. Here we see far less gas
approaching critical values, with most of the domain having significantly
subcritical line masses. The structure in these projections is therefore more
evolved than our hydrodynamical models, with only the densest areas of the
filaments being near fragmentation, and most of the gas having already
fragmented into clumps. Likewise, we see fewer filaments of transcritical line
mass in these projections. Instead, we note only the densest filaments of the
projections are close to approaching the critical value. Furthermore, we see
no elongated structures approaching supercritical values in green. Instead, we
see clumpier structure overall, and very small areas approaching supercritical
values and starting to fragment. We therefore conclude that the magnetic
fields have an overall stabilizing effect on the gas when only considering
thermal stability. In future work we intend to analyze the impact of
accounting for magnetic field affects in the critical line mass.

We can also see that, while the magnetic runs may be more efficient at
producing star-forming gas as indicated in Table 2, these projections do not
show near-critical or supercritical line masses in as large a fraction of the
total gas. Overall, the line masses of the filaments would suggest that our
non-magnetic models could have higher star formation rates due to the higher
concentration of near-critical filaments. Instead, we have shown throughout
our previous results that our magnetic models in fact have the higher star
formation and produce more cluster sinks than our non-magnetic models. In
fact, what is apparent from our line mass projections is that the magnetic
fields have acted to speed up the formation of filamentary structure. In this
way, filaments take less time to form, being able to start forming star
clusters more rapidly than the hydrodynamical models, despite unmagnetized
models starting to form clusters sooner than the magnetized.

We can conclude that cluster formation begins after filaments have reached
supercritical line masses, whereupon they begin to fragment into clumps that
ultimately host star clusters, lending support to beads on a string formation
scenarios of star clusters in molecular filaments.

###  4.7 Magnetic Field-Density Relation

Figure 15: Relation between magnetic field magnitude and volume density of the
disp8MHD model. Blue dashed line represents the simulation, and the solid
black line shows a $B\propto\rho^{2/3}$ relation, matching observations from
Crutcher et al. (2010).

It is important to track the evolution of the magnetic field in our simulation
as it is likely to evolve over time, changing it’s strength and the role it
plays on the gas. Magnetic fields may increase in strength simply due to the
compression of the gas which contains them. In this event,
$B\propto\rho^{2/3}$ discussed in Crutcher et al. (2010). On the other hand,
some cases may see magnetic field strength increase disproportionately to the
density increase. These instances are evidence of dynamo effects, where
turbulent motions in the dense gas may feed magnetic fields without needing to
increase density. In a dynamo case, one would see a deviation to steeper B vs.
$\rho$ relation from the well-accepted trend of Crutcher et al. (2010).

In Figure 15 we plot the relation between magnetic field magnitude and
density, and compare it to the scaling with $\rho^{2/3}$. At densities below
$10^{2}~{}cm^{-3}$, we have good agreement with the $2/3$ power-law relation,
but notice a significant deviation from the scaling relation at densities
above this value. As we deviate above the scaling from Crutcher et al. (2010),
we see a higher magnetic field strength. This may indicate the operation of a
turbulent dynamo effect that increases the field strength beyond values that
can be reached by compression alone. At the highest densities, the trend
flattens, with spiking present from the formation of sink particles
artificially decreasing the amount of gas at $10^{4}-10^{5}~{}cm^{-3}$ from
our formation criteria. Given the resolution limit and the formation of sink
particles, it is not possible to conclude definitively that this flattening is
a real physical result in our simulations, or a consequence of our resolution
limit.

##  5 Conclusions

In this paper we have presented simulations of the effects of the CNM
environment on molecular cloud and star cluster formation. We used initial
conditions for typical CNM densities and temperatures, although at higher
column densities than the local region - more in accord with typical massive
star-forming regions, such as the CMZ and luminous starburst galaxies. We used
the RAMSES code - that includes, MHD, cloud chemistry, turbulence, and self
gravity, but without stellar feedback - to set up MHD simulations for a range
of three different initial 3D velocity dispersions imposed on the gas but
allowed to decay. Since we chose the decay times to be much longer than our
simulation time scales (5 Myr), the turbulent amplitudes remain nearly
constant over our runs. We then analyzed the structure formed in the various
cases and found an optimal condition for structure and cluster formation. We
list our conclusions below.

Arguably our most important conclusion is that velocity dispersion in the CNM
plays a significant role in structure formation. A dispersion of 5 km/s in a
0.5 kpc cube is dominated by its self-gravity. Its weak turbulence, even
combined with thermal pressure, can barely support it. Dispersions of 8 and 10
km/s (initial energy ratios of 0.01 and 0.03 respectively) both allow for
turbulence to create filamentary structure, though 8 km/s is able to reproduce
some features that resemble observations of the Orion A GMC. Additionally, in
all three models we note that structure formation begins and is set well
before one crossing time has passed on the scale of individual clouds. The
formation of dense filaments therefore does not solely rely on the turbulent
crossing time in the gas, and we here stress the importance of gravity in
creating stable, dense structure as well.

Several important conclusions about the effects of magnetic fields are:

  * •

Although both 8 km/s models contain $\sim$50% of the total gas mass in
molecular gas of densities 100 $cm^{-3}$ or higher by the onset of sink
particle formation, the dense gas fraction between hydro and magnetic models
is very different. Just 0.2% of the total gas is converted to star-forming
densities in unmagnetized runs, whereas the inclusion of magnetic fields
doubles the fraction to 0.5%. The relative rarity of reaching star-forming
densities is mirrored in the sink particle formation, in which only a handful
of cluster sinks form before the end of our simulations, exhibiting a low star
formation rate in the early formation of clusters in the CNM.

  * •

Magnetic field orientations are primarily perpendicular to the gas above
densities of $100~{}\rm{cm}^{-3}$ once clusters have begun to form.
Dispersions of 5 and 8 km/s see magnetic field orientations return to be
parallel to the gas again above densities of $\simeq 10^{3}~{}\rm{cm}^{-3}$,
and dispersions of 10 km/s are trending towards parallel alignment at higher
densities. Our models see evidence for the effects of high density gas flows
in the centres of filaments dragging magnetic field lines into parallel
alignment, though higher dispersions may contain higher levels of turbulent
support, needing higher density flows to pull the field into parallel
alignment.

  * •

Magnetic fields also significantly contribute to the pressure of the ISM, with
average pressures an order of magnitude higher than in our non-magnetic
models. This results in an order of magnitude increase in temperature due to
the lower compressibility of the gas arising through magnetic support. Both
temperature and pressure plots showcase these compression differences in the
high-density ranges, where both are limited to a very small spread in
unmagnetized models. Additionally, Fig. 6 shows that the gas in the pure hydro
simulation overcools by an order of magnitude compared to the MHD case - a
clear indication that more of the shocked gas in our MHD models has lower
cooling rates and will be warmer. While neither hydro nor MHD models contain a
hot, ionized medium, the models do both recreate a multiphase ISM, containing
both cold and warm neutral gas, with significantly more cold gas in our hydro
models. In future work, we intend on introducing radiation effects in order to
add the hot, ionized medium and recreate the three-phase ISM.

  * •

Filament line masses are largely transcritical in unmagnetized simulations.
Including magnetic fields leaves the line masses at subcritical values in all
but the high-density gas. Additionally, magnetic fields can suppress the
amount of gas that fragments to form clusters, such that even though the
presence of fields correlates with more gas in high density regimes, the star
formation rate is suppressed from the stabilizing effects of the fields on the
line mass.

Other significant results are:

  * •

Temperatures of filaments are cold, sitting around 10-30 K, but they are
surrounded by hot gas at temperatures $\geq 600$ K, heated by the shock waves
that form the filaments. We note the intricate filamentary structure formed,
specifically the braided features that dense filaments display in density
projections, temperature slices and in magnetic field structure. These braids
are crucial for the interaction of filaments to create dense structures.

  * •

Density PDFs for models with dispersions of 8 km/s produce a standard log-
normal shape. The inclusion of magnetic fields supports the beginning of a
secondary peak, formed by the cold gas in the simulation. High densities show
evidence of power-law tails produced via self-gravity, regardless of presence
of magnetic fields.

  * •

Magnetic field strengths increase with density, where field strengths of up to
$100\mu G$ correspond to star-forming densities of $10^{4}cm^{-3}$, showing a
clear link between gas density and magnetic field density. The filamentary
structure of the ISM can be seen in the magnetic field strength, displaying
ridges of high magnetic field along the filaments, and the structure shows
perpendicular field lines funnelling these filaments to their high densities.
We also see, in our 8 km/s models, the majority of cluster formation happening
along a high density filament and, thus, in highly magnetized environments. We
link this to the delay we see in sink particle formation, as the magnetic
field strengths in the high density gas prove very dynamically important.

  * •

Turbulence produces no net, global rotation throughout the box, such that
rotation and shear in the ISM are likely created mostly via galactic dynamics.
The angular momentum magnitude per radius bin shows very low values,
indicating a low and constant momentum flux which supports our conclusions
towards low/no net rotation. Turbulence does produce local fluctuations in
angular momentum, and these are likely to be important for the formation of
disks on much smaller (and unresolved) scales.

  * •

We compare our models to Crutcher et al. (2010) through the B vs. $\rho$
relation. We note that low densities follow a standard $2/3$ power-law
relation, but high densities experience more compression, deviating above the
$2/3$ line. This may indicate some turbulent dynamo action in the denser gas,
though more work must be done to investigate thoroughly.

Finally, there are several interesting conclusions regarding the very earliest
phases of star cluster formation:

  * •

The differences in the sink particle formation onset times highlight the delay
of formation that occurs in our MHD models. Between all three dispersions, we
see a consistently later time to sink formation with the inclusion of magnetic
fields. The inclusion of magnetic fields in our simulations delays the onset
of sink particle formation by $\sim 0.4$ Myr in our 8 km/s dispersion models.
The delay does not seem to affect the amount of molecular gas formed, but
rather affects the high density star-forming gas, indicating a support role
from the magnetic fields in preventing very high density clumps from forming
as quickly. This result is further mirrored in the density PDFs, wherein our
magnetic runs exhibit a different shape at high densities than our hydro runs.

  * •

Cluster mass sink particles form quickly, and we surpass a total mass in sinks
of $10^{3}M_{\odot}$ within 0.4 Myr, indicating rapid growth. As in previous
works (see, for instance, Howard et al., 2018), mass growth is fast. After
initial rapid growth, our models see average accretion rates of
$10^{-5}~{}M_{\odot}~{}yr^{-1}$. While we only look at early evolution here,
there are already signs of a turnover in total cluster mass, signifying
cluster formation happens early on and finishes rapidly, as found in Li et al.
(2017).

## Acknowledgements

The authors thank an anonymous referee for a report that helped to clarify the
manuscript considerably. We would also like to acknowledge and thank the
following people for helpful and interesting discussions during the course of
this work: Noé Brucy, Juan Soler, Marta Reina-Campos, Hector Robinson, Claude
Cournoyer-Cloutier, Jeremy Karam, Davide Decataldo, Ralf Klessen, Alison
Sills, James Wadsley, Henrik Beuther. We would like to especially thank Romain
Teyssier for his help and support in implementing the RAMSES code. We thank
Juan Soler and Patrick Hennebelle for allowing us the use of Magnetar to
analyze magnetic field orientations. Computation was enabled through the
Graham and Cedar clusters of the Digital Research Alliance of Canada, as well
as support from their team. RP is grateful to the funding support of the
Ontario Graduate Scholarship and the Joanne and Joseph Lee Scholarship at
McMaster University. REP is supported by a Discovery Grant from NSERC Canada.

## Data Availability

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

## References

  * Abgrall et al. (1992) Abgrall H., Le Bourlot J., Pineau Des Forets G., Roueff E., Flower D. R., Heck L., 1992, Astronomy and Astrophysics, 253, 525
  * Alves de Oliveira et al. (2014) Alves de Oliveira C., et al., 2014, Astronomy & Astrophysics, 568, A98
  * André et al. (2010) André P., et al., 2010, Astronomy & Astrophysics, 518, L102
  * André et al. (2014) André P., Di Francesco J., Ward-Thompson D., Inutsuka S. I., Pudritz R. E., Pineda J. E., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 27 (arXiv:1312.6232), doi:10.2458/azu_uapress_9780816531240-ch002
  * Arzoumanian et al. (2011) Arzoumanian D., et al., 2011, Astronomy & Astrophysics, 529, L6
  * Arzoumanian et al. (2013) Arzoumanian D., André P., Peretto N., Könyves V., 2013, Astronomy & Astrophysics, 553, A119
  * Banerjee et al. (2009) Banerjee R., Vázquez-Semadeni E., Hennebelle P., Klessen R. S., 2009, Monthly Notices of the Royal Astronomical Society, 398, 1082
  * Bellomi et al. (2020) Bellomi E., Godard B., Hennebelle P., Valdivia V., Pineau des Forêts G., Lesaffre P., Pérault M., 2020, Astronomy and Astrophysics, 643, A36
  * Beuther et al. (2018) Beuther H., et al., 2018, Astronomy & Astrophysics, 617, A100
  * Bieri et al. (2022) Bieri R., Naab T., Geen S., Coles J. P., Pakmor R., Walch S., 2022, arXiv e-prints, p. arXiv:2209.06842
  * Bleuler & Teyssier (2014) Bleuler A., Teyssier R., 2014, Monthly Notices of the Royal Astronomical Society, 445, 4015–4036
  * Bleuler et al. (2015) Bleuler A., Teyssier R., Carassou S., Martizzi D., 2015, Computational Astrophysics and Cosmology, 2, 5
  * Bournaud et al. (2009) Bournaud F., Elmegreen B. G., Martig M., 2009, The Astrophysical Journal, 707, L1
  * Brown & Gnedin (2022) Brown G., Gnedin O. Y., 2022, Monthly Notices of the Royal Astronomical Society, 514, 280
  * Brucy & Hennebelle (2021) Brucy N., Hennebelle P., 2021, Monthly Notices of the Royal Astronomical Society, 503, 4192
  * Brucy et al. (2020) Brucy N., Hennebelle P., Bournaud F., Colling C., 2020, The Astrophysical Journal, 896, L34
  * Brucy et al. (2023) Brucy N., Hennebelle P., Colman T., Iteanu S., 2023, arXiv e-prints, p. arXiv:2305.18012
  * Brunetti & Wilson (2022) Brunetti N., Wilson C. D., 2022, Monthly Notices of the Royal Astronomical Society, 515, 2928
  * Calura et al. (2022) Calura F., et al., 2022, Monthly Notices of the Royal Astronomical Society, 516, 5914
  * Chevance et al. (2022) Chevance M., Krumholz M. R., McLeod A. F., Ostriker E. C., Rosolowsky E. W., Sternberg A., 2022, arXiv e-prints, p. arXiv:2203.09570
  * Colman et al. (2022) Colman T., et al., 2022, Monthly Notices of the Royal Astronomical Society, 514, 3670
  * Coppola et al. (2011) Coppola C. M., Diomede P., Longo S., Capitelli M., 2011, The Astrophysical Journal, 727, 37
  * Crutcher et al. (2010) Crutcher R. M., Wandelt B., Heiles C., Falgarone E., Troland T. H., 2010, The Astrophysical Journal, 725, 466
  * Decataldo et al. (2017) Decataldo D., Ferrara A., Pallottini A., Gallerani S., Vallini L., 2017, Monthly Notices of the Royal Astronomical Society, 471, 4476
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium
  * Falgarone et al. (2001) Falgarone E., Pety J., Phillips T. G., 2001, The Astrophysical Journal, 555, 178
  * Federrath et al. (2010a) Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2010a, Astronomy & Astrophysics, 512
  * Federrath et al. (2010b) Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010b, The Astrophysical Journal, 713, 269
  * Fiege & Pudritz (2000) Fiege J. D., Pudritz R. E., 2000, Monthly Notices of the Royal Astronomical Society, 311, 85
  * Fisher et al. (2019) Fisher D. B., Bolatto A. D., White H., Glazebrook K., Abraham R. G., Obreschkow D., 2019, The Astrophysical Journal, 870, 46
  * Galván-Madrid et al. (2013) Galván-Madrid R., et al., 2013, The Astrophysical Journal, 779, 121
  * Girichidis et al. (2016) Girichidis P., et al., 2016, Monthly Notices of the Royal Astronomical Society, 456, 3432
  * Girichidis et al. (2018) Girichidis P., Seifried D., Naab T., Peters T., Walch S., Wünsch R., Glover S. C. O., Klessen R. S., 2018, Monthly Notices of the Royal Astronomical Society, 480, 3511
  * Gómez & Vázquez-Semadeni (2014) Gómez G. C., Vázquez-Semadeni E., 2014, The Astrophysical Journal, 791, 124
  * Grisdale (2021) Grisdale K., 2021, Monthly Notices of the Royal Astronomical Society, 500, 3552
  * Großschedl et al. (2018) Großschedl J. E., et al., 2018, Astronomy and Astrophysics, 619, A106
  * Grudić et al. (2021) Grudić M. Y., Kruijssen J. M. D., Faucher-Giguère C.-A., Hopkins P. F., Ma X., Quataert E., Boylan-Kolchin M., 2021, Monthly Notices of the Royal Astronomical Society, 506, 3239
  * Grudić et al. (2022) Grudić M. Y., Guszejnov D., Offner S. S. R., Rosen A. L., Raju A. N., Faucher-Giguère C.-A., Hopkins P. F., 2022, Monthly Notices of the Royal Astronomical Society, 512, 216
  * Guszejnov et al. (2022) Guszejnov D., Markey C., Offner S. S. R., Grudić M. Y., Faucher-Giguère C.-A., Rosen A. L., Hopkins P. F., 2022, Monthly Notices of the Royal Astronomical Society, 515, 167
  * Hacar et al. (2022) Hacar A., Clark S., Heitsch F., Kainulainen J., Panopoulou G., Seifried D., Smith R., 2022, arXiv e-prints, p. arXiv:2203.09562
  * Han et al. (2022) Han D., Kimm T., Katz H., Devriendt J., Slyz A., 2022, arXiv e-prints, p. arXiv:2207.05745
  * He et al. (2022) He H., Wilson C., Brunetti N., Finn M., Bemis A., Johnson K., 2022, The Astrophysical Journal, 928, 57
  * Henning et al. (2010) Henning T., Linz H., Krause O., Ragan S., Beuther H., Launhardt R., Nielbock M., Vasyunina T., 2010, Astronomy & Astrophysics, 518, L95
  * Howard et al. (2018) Howard C. S., Pudritz R. E., Harris W. E., 2018, Nature Astronomy, 2, 725
  * Jappsen & Klessen (2004) Jappsen A. K., Klessen R. S., 2004, A&A, 423, 1
  * Jeffreson et al. (2020) Jeffreson S., Kruijssen J., Keller B., Chévance M., C.O. G. S., 2020, Monthly Notices of the Royal Astronomical Society, 498, 385*429
  * Joung et al. (2009) Joung M. R., Mac Low M.-M., Bryan G. L., 2009, The Astrophysical Journal, 704, 137
  * Kainulainen et al. (2014) Kainulainen J., Federrath C., Henning T., 2014, Science, 344, 183
  * Kim & Ostriker (2007) Kim W.-T., Ostriker E. C., 2007, The Astrophysical Journal, 660, 1232
  * Kim & Ostriker (2015) Kim C.-G., Ostriker E. C., 2015, ApJ, 815, 67
  * Kim & Ostriker (2017) Kim C.-G., Ostriker E. C., 2017, The Astrophysical Journal, 846, 133
  * Kim et al. (2011) Kim C.-G., Kim W.-T., Ostriker E. C., 2011, The Astrophysical Journal, 743, 25
  * Kim et al. (2013) Kim C.-G., Ostriker E. C., Kim W.-T., 2013, The Astrophysical Journal, 776, 1
  * Kirk et al. (2015) Kirk H., Klassen M., Pudritz R., Pillsworth S., 2015, The Astrophysical Joural, 802, 75
  * Kitsionas et al. (2009) Kitsionas S., et al., 2009, Astronomy & Astrophysics, 508, 541
  * Klassen et al. (2012) Klassen M., Pudritz R. E., Peters T., 2012, Monthly Notices of the Royal Astronomical Society, 421, 2861
  * Klassen et al. (2017) Klassen M., Pudritz R. E., Kirk H., 2017, Monthly Notices of the Royal Astronomical Society, 465, 2254
  * Klessen & Glover (2016) Klessen R. S., Glover S. C. O., 2016, Saas-Fee Advanced Course, 43, 85
  * Körtgen et al. (2019) Körtgen B., Banerjee R., Pudritz R. E., Schmidt W., 2019, Monthly Notices of the Royal Astronomical Society, 489, 5004
  * Krumholz et al. (2005) Krumholz M. R., McKee C. F., Klein R. I., 2005, The Astrophysical Journal Letters, 618, L33
  * Krumholz et al. (2014) Krumholz M. R., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. pp 243–266 (arXiv:1401.2473), doi:10.2458/azu_uapress_9780816531240-ch011
  * Kuffmeier (2018) Kuffmeier M., 2018, in Constraints from zoom-in simulations on the protostellar accretion process. Cambridge University Press, pp 56–60, doi:10.1017/S1743921319001510
  * Kuiper & Hosokawa (2018) Kuiper R., Hosokawa T., 2018, Astronomy and Astrophysics, 616, A101
  * Kwon et al. (2022) Kwon W., et al., 2022, The Astrophysical Journal, 926, 163
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, Annual Review of Astronomy & Astrophysics, 41, 51
  * Lane et al. (2022) Lane H. B., Grudić M. Y., Guszejnov D., Offner S. S. R., Faucher-Giguère C.-A., Rosen A. L., 2022, Monthly Notices of the Royal Astronomical Society, 510, 4767
  * Larson (1981) Larson R. B., 1981, Monthly Notices of the Royal Astronomical Society, 194, 809
  * Lee et al. (2022) Lee J. C., et al., 2022, The Astrophysical Journal Supplement Series, 258, 10
  * Li et al. (2017) Li H., Gnedin O. Y., Gnedin N. Y., Meng X., Semenov V. A., Kravtsov A. V., 2017, The Astrophysical Journal, 834, 69
  * Lombardi et al. (2011) Lombardi M., Alves J., Lada C. J., 2011, Astronomy and Astrophysics, 535, A16
  * Lombardi et al. (2014) Lombardi M., Bouy H., Alves J., Lada C. J., 2014, Astronomy and Astrophysics, 566, A45
  * Lue et al. (2021) Lue A., Guszejnov D., Offner S. S. R., Grudić M. Y., 2021, Research Notes of the American Astronomical Society, 5, 225
  * Mac Low & Klessen (2004) Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125
  * Massey & Hunter (1998) Massey P., Hunter D. A., 1998, The Astrophysical Journal, 493, 180
  * McClure-Griffiths et al. (2006) McClure-Griffiths N. M., Dickey J. M., Gaensler B. M., Green A. J., Haverkorn M., 2006, The Astrophysical Journal, 652, 1339
  * Men’shchikov et al. (2010) Men’shchikov A., et al., 2010, Astronomy & Astrophysics, 518, L103
  * Molinari et al. (2010) Molinari S., et al., 2010, Astronomy & Astrophysics, 518, L100
  * Motte et al. (2010) Motte F., et al., 2010, Astronomy & Astrophysics, 518, L77
  * Myers (2009) Myers P. C., 2009, The Astrophysical Journal, 700, 1609
  * Ntormousi & Hennebelle (2019) Ntormousi E., Hennebelle P., 2019, Astronomy and Astrophysics, 625, A82
  * Pattle et al. (2022) Pattle K., Fissel L., Tahani M., Liu T., Ntormousi E., 2022, arXiv e-prints, p. arXiv:2203.11179
  * Peretto et al. (2013) Peretto N., et al., 2013, Astronomy & Astrophysics, 555, A112
  * Pillai et al. (2020) Pillai T. G. S., et al., 2020, Nature Astronomy, 4, 1195
  * Pineda et al. (2008) Pineda J. E., Caselli P., Goodman A. A., 2008, The Astrophysical Journal, 679, 481
  * Pineda et al. (2010) Pineda J. E., Goodman A. A., Arce H. G., Caselli P., Foster J. B., Myers P. C., Rosolowsky E. W., 2010, The Astrophysical Journal, 712, L116
  * Pineda et al. (2011) Pineda J. E., Goodman A. A., Arce H. G., Caselli P., Longmore S., Corder S., 2011, The Astrophysical Journal, 739, L2
  * Polychroni et al. (2013) Polychroni D., et al., 2013, The Astrophysical Journal Letters, 777, L33
  * Ramírez-Galeano et al. (2022) Ramírez-Galeano L., Ballesteros-Paredes J., Smith R. J., Camacho V., Zamora-Avilés M., 2022, Monthly Notices of the Royal Astronomical Society, 515, 2822
  * Rathjen et al. (2021) Rathjen T.-E., et al., 2021, Monthly Notices of the Royal Astronomical Society, 504, 1039
  * Reina-Campos et al. (2022) Reina-Campos M., Keller B. W., Kruijssen J. M. D., Gensior J., Trujillo-Gomez S., Jeffreson S. M. R., Pfeffer J. L., Sills A., 2022, Monthly Notices of the Royal Astronomical Society,
  * Rezaei Kh. et al. (2020) Rezaei Kh. S., Bailer-Jones C. A. L., Soler J. D., Zari E., 2020, Astronomy & Astrophysics, 643, A151
  * Rice et al. (2016) Rice T. S., Goodman A. A., Bergin E. A., Beaumont C., Dame T. M., 2016, The Astrophysical Journal, 822, 52
  * Rieder et al. (2022) Rieder S., Dobbs C., Bending T., Liow K. Y., Wurster J., 2022, Monthly Notices of the Royal Astronomical Society, 509, 6155
  * Schleicher et al. (2013) Schleicher D., Schober J., Federrath C., Bovino S., Schmidt W., 2013, New Journal of Physics, 15, 023017
  * Schmidt et al. (2009) Schmidt W., Federrath C., Hupp M., Kern S., Niemeyer J. C., 2009, Astronomy and Astrophysics, 494, 127
  * Seifried et al. (2017) Seifried D., et al., 2017, Monthly Notices of the Royal Astronomical Society, 472, 4797
  * Smith et al. (2012) Smith R. J., Shetty R., Stutz A. M., Klessen R. S., 2012, The Astrophysical Journal, 750, 64
  * Smith et al. (2017) Smith B. D., et al., 2017, Monthly Notices of the Royal Astronomical Society, 466, 2217
  * Smith et al. (2021) Smith M. C., Bryan G. L., Somerville R. S., Hu C.-Y., Teyssier R., Burkhart B., Hernquist L., 2021, Monthly Notices of the Royal Astronomical Society, 506, 3882
  * Soler (2019) Soler J. D., 2019, Astronomy & Astrophysics, 629, A96
  * Soler et al. (2013) Soler J. D., Hennebelle P., Martin P. G., Miville-Deschênes M. A., Netterfield C. B., Fissel L. M., 2013, The Astrophysical Journal, 774, 128
  * Soler et al. (2022) Soler J. D., et al., 2022, Astronomy and Astrophysics, 662, A96
  * Solomon (2001) Solomon P. M., 2001, in Tacconi L., Lutz D., eds, Starburst Galaxies: Near and Far. p. 173 (arXiv:astro-ph/0101482), doi:10.48550/arXiv.astro-ph/0101482
  * Sun et al. (2022) Sun J., et al., 2022, arXiv e-prints, p. arXiv:2206.07055
  * Tahani et al. (2022) Tahani M., et al., 2022, Astronomy & Astrophysics, 660, L7
  * Teyssier (2002) Teyssier R., 2002, Astronomy and Astrophysics, 385, 337
  * Theissen et al. (2022) Theissen C. A., Konopacky Q. M., Lu J. R., Kim D., Zhang S. Y., Hsu C.-C., Chu L., Wei L., 2022, The Astrophysical Journal, 926, 141
  * Thilker et al. (2021) Thilker D. A., et al., 2021, arXiv e-prints, 1, arXiv:2106.13366
  * Townsley et al. (2014) Townsley L. K., Broos P. S., Garmire G. P., Bouwman J., Povich M. S., Feigelson E. D., Getman K. V., Kuhn M. A., 2014, The Astrophysical Journal Supplement, 213, 1
  * Tress et al. (2020) Tress R. G., Sormani M. C., Glover S. C. O., Klessen R. S., Battersby C. D., Clark P. C., Hatchfield H. P., Smith R. J., 2020, Monthly Notices of the Royal Astronomical Society, 499, 4455
  * Turner et al. (2022) Turner J. A., et al., 2022, Monthly Notices of the Royal Astronomical Society, 516, 4612
  * Wang et al. (2020) Wang Y., et al., 2020, Astronomy and Astrophysics, 634, A139
  * Watkins et al. (2023a) Watkins E. J., et al., 2023a, arXiv e-prints, p. arXiv:2302.03699
  * Watkins et al. (2023b) Watkins E. J., et al., 2023b, The Astrophysical Journal, 944, L24
  * Wilson et al. (1970) Wilson R. W., Jefferts K. B., Penzias A. A., 1970, The Astrophysical Journal Letters, 161, L43
  * Wilson et al. (2019) Wilson C. D., Elmegreen B. G., Bemis A., Brunetti N., 2019, The Astrophysical Journal, 882, 5
  * Wilson et al. (2023) Wilson C. D., Bemis A., Ledger B., Klimi O., 2023, Monthly Notices of the Royal Astronomical Society, 521, 717
  * Zhao et al. (2023) Zhao B., Pudritz R. E., Pillsworth R., Robinson H., Wadsley J., 2023, Filamentary Hierarchies and Supperbubbles: Galactic Multiscale MHD Simulations of GMC to Star Cluster Formation, submitted
  * Zucker et al. (2021) Zucker C., et al., 2021, The Astrophysical Journal, 919, 35

##  Appendix A Density and temperature PDFs

Figures 16 and 17 show the isolated density and temperature PDFs for the 5
km/s and 10 km/s models, respectively. Dotted vertical lines pinpoint the
average density for all gas in the simulation, while dashed vertical lines
pinpoint the average density when considering only the gas contained in or
below the peak temperature. We choose this temperature cutoff based on where
the temperature curve has begun to flatten out.

For both the 5 and 10 km/s models, these averages are approximately the same,
as evidenced by the overlapping vertical lines in the plots. Additionally,
while these models show different average densities between their magnetized
and unmagnetized runs, they have similar average densities to each other
despite their different density ranges. Our 10 km/s models show more gas at
low densities and reach maximum densities of $\simeq 10^{5}~{}\rm{cm}^{-3}$.
The 5 km/s models, on the other hand, do not contain densities much lower than
$\simeq 1~{}\rm{cm}^{-3}$, but reach high densities of
$10^{6}~{}\rm{cm}^{-3}$.

Though the high density range above $10^{5}~{}\rm{cm}^{-3}$ is affected by the
resolution limit of the simulation, we can still discern some general
information. First, the higher density limit in the 5 km/s case shows that the
lower total number of sink particles formed in the simulation is not due to
high density gas not being present at all, but rather due to that high density
gas never forming collapsing clumps. On the other hand, and to our second
point, the 10 km/s models are not necessarily missing sink particles because
of a lack of high density gas either, but more likely due to high turbulence
creating shocks which are actively tearing high density filaments apart,
preventing clumps from forming. Compared to our results for the 8 km/s model
in Section 4, these indicate that the balance between the gravity and
turbulence plays a key role in all structure formation and directly correlates
to cluster formation. We see no direct evidence of the magnetic field
affecting these balances in this case, though we also not that a more detailed
study of the magnetic field strength and its effects would be necessary to
make strong conclusions.

Figure 16: Density and temperature PDFs for our 5 km/s models. Dotted vertical
lines on the density PDF (top) mark the average density of all gas. Dashed
lines represent the average density of cold gas in the simulation. Drastically
low temperatures in the disp5 hydro model showcase the overcooling of the gas
brought on by significantly higher compression allowed due to the lack of
magnetic support, as discussed in Sec. 4.2.

Figure 17: Density and temperature PDFs for our 10 km/s models. Dotted
vertical lines on the density PDF (top) mark the average density of all gas.
Dashed lines represent the average density of cold gas in the simulation.

For completeness and direct comparison, we include in Figure 18 the density
and temperature PDFs of all six models together. For similar density PDFs, all
of our models have a unique temperature PDF. Th unmagnetized hydro models
follow an increasing peak temperature with increasing velocity dispersion, as