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

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 $320 M_{\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.


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).
★ Contact e-mail: pillswor@mcmaster.caIt 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.

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.

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 (∼ 10 21  −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: This gives us a virial line mass, as discussed in Fiege & Pudritz (2000), defined as: 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 , and critical velocity dispersion  2  = 4  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 molec- 10 .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   ,   = 128 ⊙  −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    ∝  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 M ⊙ 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 × 10 21  −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 filamen- tary 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.

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.
In Fig. 1, Chevance et al. show the dust emission and 12  (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.

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. (2020Brucy et al. ( , 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.

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  −1 in simulation domains of size 0.5×0.5×2kpc.Dispersions of this magnitude are often associated with turbulence, as the typical turbulent rms velocity is    ≈ 5   , which becomes the dominant speed compared to typical sound speeds of   ≈ 0.2 − 1.0   .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.

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 × 0.5kpc × 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. 2016Girichidis et al. , 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. 2011Pineda et al. , 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  () of the turbulence, where k is the wavenumber: 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 C1 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.
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 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  −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 320 ⊙  −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(Wilson et al. , 2019)).

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  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 × 10 −3   −1  −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  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 × 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 10 4 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 ∼ 10 21 cm −2 for the diffuse ISM) due to self-shielding.In many cases the dissociation rate of H 2 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 H 2 , in order to have enough present that its destruction would be necessary, is on the order of 10 9 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 H 2 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 H 2 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 10 21 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.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 , in accordance with average values of ISM magnetic field strengths.  , and   , represent the velocity dispersion in the simulations on initialization and once sink particles have formed, respectively. , ℎ  and  ,  give the initial Jeans lengths of each simulation calculated using the average temperature and the velocity dispersion, respectively.

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 −3 and a minimum mass of 100 ⊙ , 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 −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   = 3 2 /.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 (  , ) and at the onset of sink formation (  , ) 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.

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 ∼ 12  −1 , whereas molecular gas may have velocity dispersions ∼ 5  −1 .Taking the average, assuming completely equal portions of HI and molecular hydrogen gives a 3D velocity dispersion of 8.5  −1 .What is the origin of these values?This question motivated our first set of numerical experiments.

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 where   is the vertical component of the velocity dispersion, Σ is the average column density of the gas, and  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 ∼ 5  −1 , and a maximum of ∼ 10km 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   −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 ∼ 4 − 6  −1 , with some regions reaching in 13 CO as high as 10  −1 .Averaging between HI and 13 CO gives dispersion ranges of 6-12  −1 .
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 ∼ 100, 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.
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 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 ∼ 10.We see higher temperatures (≤ 600) 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 ∼ 10 − 30 (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 ∼ 1, 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 10.The presence of a large range of temperatures (spanning 1 − 1 × 10 3 ), 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 Σ and  in Equation 5, we find   ≃ 30km 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  ∝  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. (2013Kim et al. ( , 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 ≈ 1−10 M ⊙ /pc 2 and an effective scale height of ∼ 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 Σ 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.

Molecular gas and cloud morphology
In our models we define our molecular gas at densities of 100  −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  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 ∼ 10 2  −3 and an average temperature of ∼ 10, lending an average sound speed of 0.3 −1 throughout the simulation.
at densities of 100, 1000 and 10000  −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 ∼ 50% of our gas mass lies at densities of 100  −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.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.
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 ∼ 7.57 × 10 4  ⊙ (Lombardi et al. 2011), whereas the mass of the molecular gas contained in this cut out is ∼ 5.78×10 4  ⊙ .Our entire zoom region is larger (150 pc) compared to the size of Orion A, which is ∼ 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 ∼ 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.
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 highdensity 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 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.
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. A1, 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.

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  parameter, which describes the concavity of the histogram of ().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 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 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 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 8 presents the structure of the magnetic field via its magni- tude 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.
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.
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 ) 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.
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.
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  −3 or higher, when we increase to 1000  −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  −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 100 −3 and 1000 −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.

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 ∼ 10 3  −3 , up to a maximum of ∼ 10 5  −3 with a minimum of ∼ 10 1  −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 ∼ 10 4  −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  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 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)).

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.
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 ≥ 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 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.
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 disper-sion 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.

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  ⊙ /, 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  ⊙ / and up to 10 2  ⊙ / (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 × 10 5  ⊙ , comparable to the 7.5 × 10 4  ⊙ 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  ⊙ 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 × 10 3  ⊙ , whereas its lowest mass cluster is roughly 10 2  ⊙ .Lowest and highest masses of each of our MHD models can be seen in Table 4.

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).
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;    ∝  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 nearcritical 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.

Magnetic Field-Density Relation
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,  ∝  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.  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  2/3 .At densities below 10 2  −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  −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.

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 ∼50% of the total gas mass in molecular gas of densities 100  −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 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 ≃ 10 3 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 ≥ 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 correspond to star-forming densities of 10 4  −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.  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 ∼ 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  ⊙ 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  ⊙  −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).line masses in order to show a map of the criticality throughout the filaments.
In Sec.4.6, we show the projected line masses of both of our 8 km/s models.In that case we have taken all the cells stacked on top of each other in the direction of the projection, and integrated the line mass field over the path length.Due to this integration, we then correct for the extra term of total path length by dividing the value in each cell by the length of the box.We note here that this is relatively easy to do given that our box is symmetrical, but non-symmetric boxes would of course require the definition of as many fields as projection directions one wishes to take.From this correction, we thus have both 3D and 2D planar information about the critical line masses of filaments, and we note the similarities in criticality between the slices above and the projections in the body of this paper.

APPENDIX C: TURBULENT ENERGY SPECTRA
We show here our turbulent energy spectra for our disp8MHD model at times of 1.90 Myr and 3.79 Myr.These show that our turbulent energy is not fully saturated to the Burger's relation until sink formation begins, therefore no appreciable decay happens in the timescales we study.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 2 .
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 3 .
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 ∼ 10 2  −3 and an average temperature of ∼ 10, lending an average sound speed of 0.3 −1 throughout the simulation.

Figure 4 .
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  •  −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  −3 to 10 23  −3 to mimic resolution range of Lombardi et al. (2014).Sink particles are denoted by filled black circles.

Figure 5 .
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.

Figure 7 .
Figure 7. Magnetic field orientations relative to filaments shown through  -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  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 .
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.

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

Figure 11 .
Figure 11.Mass growth of sink particles over time shown for MHD models disp8MHD and disp10MHD.

Figure 12 .
Figure12.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.

Figure 13 .
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 .
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   = 2 2  /, thereby indicating only thermal stability.Line mass is taken as cell mass/cell width and sampled for each cell, therefore providing a lower bound.

Figure 15 .
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  ∝  2/3 relation, matching observations from Crutcher et al. (2010).

Figure A2 .
Figure A2.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.

Figure A3 .
Figure A3.Density and temperature PDFs for all of our models.

Figure B1 .
Figure B1.Examples of slices through the cube of our disp8 and disp8MHD models showing the line mass in evaluated in the slice cell by cell.

Figure C1 .
Figure C1.Turbulent kinetic energy spectra for our disp8MHD model at a time of 1.90 Myr (top) and 3.79 Myr (bottom).The dotted line represents the Kolmogorov spectrum  ( ) ∝  −5/3 .The dash dotted line shows the Burgers spectrum  ( ) ∝  −2 .The maximum wavenumber (k) corresponds to a cell of 0.48 pc width (our highest resolution), but we note that some resolution effects still arise in large k values.These spikes in energy are likely due to small numbers of cells of the specific corresponding resolution. Table

Table 3 .
Comparison between time of first sink formation (  ) and 10% of the global crossing time ( ,10%