Gravity versus Magnetic Fields in Forming Molecular Clouds

Magnetic fields are dynamically important in the diffuse interstellar medium. Understanding how gravitationally bound, star-forming clouds form requires modeling of the fields in a self-consistent, supernova-driven, turbulent, magnetized, stratified disk. We employ the FLASH magnetohydrodynamics code to follow the formation and early evolution of clouds with final masses of 3–8 × 103 M ⊙ within such a simulation. We use the code’s adaptive mesh refinement capabilities to concentrate numerical resolution in zoom-in regions covering single clouds, allowing us to investigate the detailed dynamics and field structure of individual self-gravitating clouds in a consistent background medium. Our goal is to test the hypothesis that dense clouds are dynamically evolving objects far from magnetohydrostatic equilibrium. We find that the cloud envelopes are magnetically supported with field lines parallel to density gradients and flow velocity, as indicated by the histogram of relative orientations and other statistical measures. In contrast, the dense cores of the clouds are gravitationally dominated, with gravitational energy exceeding internal, kinetic, or magnetic energy and accelerations due to gravity exceeding those due to magnetic or thermal pressure gradients. In these regions, field directions vary strongly, with a slight preference toward being perpendicular to density gradients, as shown by three-dimensional histograms of relative orientation.


INTRODUCTION
Understanding the relative importance of magnetic fields and self-gravity in the for-man splitting (e.g.Heiles 2004), which gives the line-of-sight field strength, and starlight polarization (e.g.Mathewson & Ford 1970;Heiles 1976), which traces the projected orientation of the field in the plane of the sky (Lazarian 2007;Andersson et al. 2015).The ionization of the diffuse gas is high enough even in neutral regions for the magnetic field to be effectively coupled to the gas (Tielens 2010;Draine 2011;Klessen & Glover 2016;Girichidis et al. 2020).As a result, uniform gravitational collapse will result in rapid growth of the field, preventing further collapse long before stars form (Mestel & Spitzer 1956;Mouschovias & Spitzer 1976).This led to the proposal that collapse and star formation could only occur through the mediation of ambipolar diffusion (ion-neutral drift; Mestel & Spitzer 1956;Mouschovias 1977;Shu 1977).
Molecular Zeeman splitting observations reveal magnetic field strengths for density ranges depending on the observed tracer molecule (e.g.Crutcher et al. 1975;Falgarone et al. 2008).The inferred values are typically insufficient to support the clouds against collapse, contradicting the expectation from ambipolar diffusion (Crutcher 1999(Crutcher , 2012)).Padoan & Nordlund (1999) showed that the observational characteristics of dense clouds were consistent with super-Alfvénic turbulence within them, leading to the gravoturbulent model of magnetized turbulence balancing gravitational collapse (Mac Low & Klessen 2004).More recent observations have emphasized the filamentary nature of the clouds (e.g.André et al. 2014) and have investigated the role of magnetic fields in filament formation (e.g.Planck Collaboration et al. 2016;Monsch et al. 2018;Fissel et al. 2019;Soler 2019), raising the question of whether the fields align with the filaments, as might be expected for weak fields, or are perpendicular, as expected for strong fields.
The observed turbulent motions cannot be driven from within the clouds (Brunt et al. 2009) but are instead likely associated with the very process of cloud assembly from the tenuous ISM (Klessen & Hennebelle 2010).Ibáñez-Mejía et al. (2016, hereafter Paper I) demonstrated that the supernova-driven turbulence in the diffuse ISM provides insufficient momentum to the dense clouds to drive the observed level of turbulence within them.Instead, they found that the observed motions can be driven by gravitational collapse, as suggested by Ballesteros-Paredes et al. (2011), leading to the conclusion originally proposed by Zuckerman & Palmer (1974).This leads to the scenario of molecular cloud formation during global hierarchical collapse (e.g.Elmegreen 2000;Vázquez-Semadeni et al. 2019), promptly followed by destruction of the clouds by stellar feedback (e.g.Dale et al. 2015;Haid et al. 2019;Wall et al. 2020;Grudić et al. 2021).Mac Low et al. (2017) demonstrated that Toomre instability of the gas disk (Goldreich & Lynden-Bell 1965) can form dense gas sufficiently quickly to replace clouds promptly destroyed by feedback in order to maintain the observed dense gas fraction in a steady state.
What magnetic field structure does global hierarchical collapse produce?Is this consistent with observed field properties?In this paper we study these questions using high-resolution, zoom-in simulations of individual clouds drawn from the kiloparsec-scale simulations of Paper I. The dynamics of accretion and collapse in these objects was described in Ibáñez-Mejía et al. (2017, hereafter Paper II), while their velocity structure functions were analyzed by Chira et al. (2019).Similar zoom-in simulations were described by Seifried et al. (2017) and Girichidis (2021), while their field structure was investigated by Seifried et al. (2020) and Girichidis (2021).
We describe the simulations in Sect. 2. We determine the physical state of the clouds in Sect. 3 using measurements of the sonic and Alfvénic Mach number, the different components of the energy, and the forces acting on each parcel of gas, all as a function of density.We then compare our models to observations in Sect.4, focusing on the relationship between line-of-sight magnetic field and density, and the alignment between magnetic field directions and density gradients as a function of density.
The analysis scripts and intermediate data products leading to the results presented here have been deposited in the digital repository of the American Museum of Natural History (Ibáñez Mejía et al. 2021).The underlying simulation results required to run the scripts are available from the same repository (Chira et al. 2018), and are also available through the Catalogue of Astrophysical Turbulence Simulations (Burkhart et al. 2020).

NUMERICAL SIMULATIONS
We analyze models of gravitationally collapsing clouds embedded in a kiloparsec-scale model of three-dimensional, stratified, magnetized, supernova-driven turbulence simulated in the FLASH v4.2.2 adaptive mesh refinement code (Fryxell et al. 2000;Dubey et al. 2012).The computational domain of the embedding simulation is 1 × 1 × 40 kpc 3 , with a background gravitational potential representing a stellar disk and dark matter halo centered vertically.Supernovae (SNe) are set off at the Galactic rate (Tammann et al. 1994) with distributions characteristic of Type Ia SNe, and both isolated and clustered core-collapse SNe.Diffuse far-ultraviolet photoelectric heating and radiative cooling are both included, as is an ideal treatment of the magnetic fields.The details of this simulation, including the initial conditions, are described in Paper I.
After establishing an equilibrium background flow at 0.95 pc resolution in the disk midplane (|z| < 300 pc), self-gravity was implemented using the FLASH multigrid solver (Daley et al. 2012).At this point we have a well-established, multiphase ISM, vertical gas stratification, and a galactic fountain up to ±20 kpc (Hill et al. 2012;Walch et al. 2015;Girichidis et al. 2016).So far, 7,515 SNe have exploded in the simulation, and they continue to be injected in the subsequent evolution presented in this work.Our analysis is performed from the moment t SG = 0 when self-gravity is turned on.
Three clouds were chosen for simulation at much higher resolution.The simulation was restarted prior to the formation of those clouds, with a refinement rule to increase resolution in any region where the Jeans length was resolved with less than four cells, down to a resolution of 0.06 or 0.1 pc.For the latter, this corresponds to resolving a density of 8 × 10 3 cm −3 for gas at 10 K (Eq.( 15) of Paper II).Our chosen clouds had masses in the simulation without self-gravity of 3 × 10 3 , 4 × 10 3 , and 8 × 10 3 M ; we identify them as M3, M4, and M8 hereafter.The first two were resolved to 0.06 pc, whereas M8 was only resolved to 0.1 pc.The clouds continue to accrete mass throughout their lifetimes.Further details of the high-resolution, zoom-in, re-simulation, and evolution of the hierarchically collapsing clouds are found in Paper II.
As illustration, Figure 1 shows a threedimensional rendering of all the gas above 100 cm −3 near the midplane, with the location of the high-resolution cloud M3, a close-up of the cloud region, and a second close-up of a dense, collapsing core within the cloud.The zoom-in view includes 125 magnetic field lines that are highly tangled and twisted, penetrating the density structure down to the massive, collapsing cores.
We stop the simulations at an evolutionary time of t SG ≈ 10 Myr, as it is expected that by this time massive stars must have already formed and be feeding back energy in the form of radiation, winds, and SN explosions, which should influence not only the cloud properties but its environment.As we do not include selfconsistent star formation and feedback in our simulations, running for longer would lead to unphysical results.

PHYSICAL ANALYSIS
In order to determine which physical processes dominate the dynamics of the gas in clouds and their envelopes during their formation, evolution, and collapse, we divide our analysis into three parts.First, in Sect.3.1 we consider the mean velocities at different densities to allow us to determine where transitions in the sonic and Alfvénic Mach numbers of the flow occur.Second, in Sect.3.2 we compare the distribution of the energy density of the gas in its different forms-magnetic, kinetic, thermal, and gravitational-as a function of the gas density, in and around the collapsing clouds.Third, in Sect.3.3 we compute the actual accelerations acting on cloud gas to directly determine which forces dominate the flow in each density regime, confirming the results suggested by the earlier analyses.
The cloud is defined by a number density contour of n c = 100 cm −3 and a gravitational binding criterion for nearby fragments (Paper II).A 100 pc cube containing the cloud is followed through the larger simulation domain and projected onto a uniform Cartesian grid for analysis using the yt toolkit Turk et al. (2011), as described in Paper II.

Velocities
We compute three velocities as a function of number density n in the frame of the cloud on the uniform grid constructed from the original adaptive grid.These are the turbulent rms velocity v rms , the sound speed c s , and the Alfvén velocity v A .For each cell i we compute Here, v i , B i , ρ, and T i are the local velocity, magnetic field, density, and temperature in the cell, respectively, and v c is the center of mass velocity of the cloud under consideration.Furthermore, k B is the Boltzmann constant, and we take γ = 5/3 (consistent with either monatomic gas or molecular gas too cold to excite rotational and vibrational levels).We also adopt a mean mass per particle µ = 1.3 m p , with m p being the proton mass.Finally, we subdivide the number density range covered by the simulation into 100 logarithmic bins of size δ(log n) between n = 10 −2 cm −3 and 10 4.3 cm −3 .We obtain the average values as a function of number density n by summing over all N cells j in the range from log n − δ(log n)/2 to log n + δ(log n)/2, We also compute the 25th and 75th percentile values to indicate the amount of variation at each number density.
In this analysis we have neglected to remove any uniform rotation of the cloud (as was done, for example, by Federrath et al. 2016), although we have subtracted the bulk motion.Because the analyzed cube extends well beyond the cloud, subtracting off an average rotation across the cube is unlikely to capture the dynamics of the cloud itself.In order to assess how much of the rms velocity can be attributed to uniform rotation, we further compute the mass-weighted average rotation of all gas with n > 100 cm−3 in the analysis cube.
To do this, we start by computing, for each zone q with mass m q = ρ q δV q lying within the high-density region, the components (L x,q , L y,q , L z,q ) of the angular momentum vector where ijk is the Levi-Civita symbol.The net angular momentum of the cloud is then L = q L q .We similarly compute the mass M cl = q m q and the spherical radius The average rotation velocity is then The variation of these velocities as a function of density and time shows how the system reacts as the dense gas contracts.The two dimensionless values that determine its behavior are the sonic Mach number M = v rms /c s and the Alfvénic Mach number M A = v rms /v A .Supersonic turbulence with M > 1, characteristic of the ISM (Mac Low & Klessen 2004), leads to large density fluctuations, as isothermal shocks produce density increases ∆n ∝ M 2 .Supersonic but sub-Alfvénic flows can compress gas along field lines, but cannot strongly influence the field (see e.g., Beattie et al. 2021).Flows that are both supersonic and super-Alfvénic can also compress gas perpendicular to magnetic fields, leading to strong fluctuations of the magnetic field strength, as the field gets compressed by the shock along with the gas (see, e.g., Federrath et al. 2008Federrath et al. , 2010)).
Figure 2 shows the characteristic velocities of the system as a function of the density, at three evolutionary times of the cloud.At each time, the rms velocity exceeds the sound speed of the system at all but the very lowest densities in the hot, diffuse medium.The flow ranges from mildly supersonic (M ≈ 1) in the diffuse ISM at number densities 10 −2 < n < 1 cm −3 , to hypersonic, M ≈ 5-10 for n > 10 cm −3 .This reproduces the expected general behavior of the turbulent ISM with supersonic shocks permeating the gas at all scales and densities (Mac Low & Klessen 2004;Padoan & Nordlund 2011;Vázquez-Semadeni 2015;Klessen & Glover 2016;Girichidis et al. 2020).As density increases, the transition to hypersonic Mach numbers has been shown to occur at roughly the same density as the transition from atomic to molecular gas (Mandal et al. 2020).
The development of the Alfvén number, on the other hand, shows the changing nature of the cloud as self-gravity takes hold.At all times the rms velocity remains trans-Alfvénic with M A ≈ 1 at n < 10 2 cm −3 .At t = 0, the rms velocity at the highest densities remains flat, while the Alfvén velocity drops, suggesting constant field strength as density increases, consistent with compression primarily along field lines.However, as the density goes up and gravitational acceleration becomes stronger, the Alfvén velocity increases again.
The flow velocity grows more quickly, suggesting gravitational collapse has become the dominant process, producing flows strong enough to compress the field lines kinematically and increase the magnetic flux in the core region.Comparing the flow velocity to the rotational velocity shows that although rotation is present, it is not the dominant component of the rms velocity at the highest densities.The collapse results in strongly super-Alfvénic flows with M A 1 at n > 10 2 cm −3 .Such super-Alfvénic flows in the dense regions of clouds were already proposed by Padoan & Nordlund (1999) based on the morphology and dynamics of observed and simulated clouds.

Energetics
To more directly examine the relationship between thermal and magnetic pressure and gravitational potential, we measure the different forms of energy as a function of number density over time.The kinetic, thermal, and magnetic energy density in each cell i is To calculate the gravitational energy density of the cloud, we first must determine the background potential Φ 0 within which it is em-bedded.We calculate this as the mean cellweighted value of the potential in the cells with n i < 10 cm −3 in a 100 pc box centered on the cloud.The gravitational energy in cell i is then We find the distribution of the kinetic, thermal, magnetic, and gravitational energy densities in each number density bin log n k , similar to the velocities in Sect.3.1, by calculating the values in cells j in the range from log n k −δ(log n k )/2 to log n k + δ(log n k )/2. Figure 3 shows the evolution in time of the median value of each energy density along with the 25th and 75th percentage values as a function of density for each of the three clouds.This complements the picture offered by the mean velocities of the physical processes affecting the cloud dynamics.
The thermal energy peaks at a density of n ≈ 0.1 cm −3 and decreases at higher densities as temperatures drop.The magnetic energy is lower than the kinetic energy at all times, and is maintained at a fraction of 0.1 − 0.3 of the total kinetic energy, probably due to the saturation of the small-scale, turbulent dynamo in our simulation box (Balsara et al. 2004;Meinecke et al. 2014;Gent et al. 2021).The small-scale dynamo drives flux growth at wavelengths close to the dissipation scale, which is determined in these models by the numerical resistivity.
As the simulation evolves, the gas with n 100 cm −3 is dominated by gravitational potential energy, implying that the clouds are gravitationally bound and in a state of hierarchical contraction, not supported by turbulent, magnetic, or thermal energy.As the clouds contract, gas flows towards local centers of collapse.These increase in density, resulting in a rise of the gravitational potential energy density at higher densities, followed by a rise of the kinetic energy density as the gas gains velocity falling down the gravitational potential wells.Still, the kinetic energy density within dense regions remains at a fraction 0.2-0.5 of the gravitational potential energy density, closely following its mass density distribution.This is expected for gravitational collapse and supports our hypothesis (Paper I) that hierarchical contraction drives the fast, non-thermal, motions of the dense gas in the cloud.The magnetic field energy also increases as the cloud contracts, as field lines are compressed in local centers of collapse.However, magnetic energy density does not grow as fast as kinetic energy density, suggesting that a significant amount of gas moves along field lines.
At densities corresponding to the outer envelope of the cloud, i.e. n = 10-100 cm −3 , the kinetic energy contributes the majority of the total energy followed by the magnetic energy.This suggests that the accretion onto the cloud proceeds generally along the magnetic field lines without compressing them strongly.
At every time in the evolution of the simulation, the kinetic energy dominates the energy of the gas below n ≈ 10 cm −3 , as expected in the mildly supersonic and super-Alfvénic diffuse ISM.Large variations in the kinetic, thermal, and magnetic energies are observed.In most cases, the warm medium with 0.1 < n < 10 cm −3 has a higher magnetic energy than thermal energy as expected (e. g.Heiles & Troland 2005), while in the lower density hot medium thermal energy tends to dominate.Equipartition is sometimes reached, but is by no means universal, with order of magnitude deviations also occurring in this diffuse gas.This occurs because we perform our analysis in a (100 pc) 3  box, with open boundaries to a turbulent, multiphase ISM.The inflow and outflow of blast waves and rarefied gas from nearby SN explosions is the main cause for these large fluctuations.

Accelerations
Finally, in order to confirm our results from the previous subsections, we determine the acceleration a of the gas in each cell i from the pressure and gravitational gradients and from the magnetic force, with the mass density ρ i = µ n i .Similar to the velocity discussed in Sect.3.1 we determine the mean acceleration as function of number density n by averaging over all cells in the range log n − δ(log n)/2 to log n + δ(log n)/2.The result is presented in Figure 4, which shows that all three clouds evolve in qualitatively similar manner.At the lowest densities, magnetic and pressure forces roughly balance each other, while local gravitational forces are unimportant.At intermediate densities, where the gas has radiatively cooled, removing pres-sure support, magnetic forces dominate.At early times, the highest density regions have been assembled by turbulent flows.Gravitational forces already dominate the accelerations there, once self-gravity is turned on (by our choice of the region) driving the accretion and collapse of the cloud.The lack of hydrostatic equilibrium in cloud cores was also shown by Girichidis (2021) using lower-resolution models of more clouds.As time passes, the highest density regions become increasingly gravitationally dominated as the cloud continues to accrete and collapse.Thus, these results remain consistent with the picture of magnetically dominated envelopes and gravitationally dominated cores.

DISCUSSION
In this discussion we turn to quantities that we derive from the three zoom-in clouds that can more directly be compared to observations, even if they are not as easily interpretable physically.In Sect.4.1 we focus on the magnetic field strength as a function of density for the diffuse and dense gas in the ISM, with and without gas self-gravity, and compare our measurements with Zeeman observations of magnetic field strengths of MCs in the Galaxy.In Sect.4.2 we compute the relative orientation of the magnetic field to two other quantities: the velocity field (Sect.4.2.2), which establishes how the field influences the flows of gas in the diffuse and dense ISM, and the density gradient (Sect.4.2.1), which quantifies the participation of the field in the structure and collapse of dense MCs.We also address some of the caveats and shortcomings of the current approach.

Density-Magnetic Field Relation
We investigate the correlation between the magnetic field strength and the density in our high-resolution zoom-in models of collapsing clouds using a B-n scatter plot of our model over time.Figures 5 and 6 show mass weighted plots of gas density for all cells within the zoom region versus the magnetic field strength, either x-component (parallel to the line of sight) or total, measured along an arbitrary line of sight lying in the galactic midplane.These figures also include Zeeman observations towards diffuse and dense clouds compiled by Crutcher et al. (2010a) as well as the upper envelope of the magnetic field strength inferred from the observations from that paper.
At t = 0, we show the distribution of magnetic field strength as a function of density in the diffuse ISM, as described in Paper I. The bulk of the gas has field strength |B x | systematically below 10 µG, the upper limit for the magnetic field strength derived from H i observations in the diffuse ISM (Crutcher 2012).For number densities n in the range from 1-300 cm −3 , the magnetic field does not correlate with the gas density, in agreement with the observations.The data show a change in behavior at n ≈ 300 cm −3 , where the derived maximum magnetic field strength begins to climb as |B| ∝ n 0.65±0.05(Crutcher et al. 2010b).
Indeed, we find a similar increase of the field strength with density.To further quantify this behavior, we fit both the average and maximum field strength for densities from 10 2 -2×10 4 cm −3 using a power law of the form As tabulated in Figures 5 and 6, we find exponents α < 0.4, which are shallower than expected theoretically for either diffusiondominated or freely-collapsing clouds, and which are also shallower than what is inferred from the observations by Crutcher et al. (2010b).This behavior could result from two effects.First, as we approach the Jeans resolution limit, numerical diffusion becomes important in limiting field growth.Higher resolution models refined with more cells per Jeans length (Sur et al. 2010) may yield higher field strengths at high  Mean accelerations from gravitational potential gradients (blue dash-dot), magnetic forces (green dash), and pressure gradients (orange solid) as a function of density in the gas within a (100 pc) 3 box, centered on the cloud's center of mass for clouds M3 (left), M4 (center), and M8 (right) for three snapshots at t = 0, 2.5, and 4.8 Myr since the time self-gravity has been turned on.Lines show median values in each density bin, while shading shows the central 50% of values in the bins.
densities.Second, our models do not include any treatment of radiative transfer or chemical abundance variations that could selectively emphasize particular regimes.Indeed the shallow slopes of the maximum field (fits to the green dashed lines in Figs. 5 and 6) result in part from maximum field values that are higher than predicted by the Crutcher et al. (2010b) analysis at intermediate densities around n ∼ 10 3 cm −3 .However, the observations are quite sparse in this regime due to a lack of chemical tracers, so the behavior in this regime is actually quite weakly constrained.
Models of magnetized cloud collapse suggest that the turbulence during collapse does drive a small-scale dynamo that produces field growth (Brandenburg & Subramanian 2005;Schleicher et al. 2010;Federrath et al. 2011a; Schober et al.  (Crutcher et al. 2010b).The black dashed line with error bars gives the mean and dispersion of the simulated values; a linear fit gives the power-law index listed to the upper left.The green dashed line fits the maximum magnetic field at each density between the transition point to power-law behavior (n = 10 2 cm −3 ) and the Jeans resolution limit (n = 10 4 cm −3 ), which also yields a power-law index.

2012
).However, Xu & Lazarian (2020) argue that at saturation of the dynamo, reconnection diffusion reduces further growth rates below the ideal MHD limit, possibly consistent with our simulation results.Further observational investigation of the intermediate density region with n 10 3 cm −3 might demonstrate this behavior.

Histogram of Relative Orientations
We examine the morphology of the magnetic field with respect to the cloud's density structure, velocity field and gravitational forces using the histogram of relative orientations (HRO) as introduced by Soler et al. (2013) to studies of the ISM (see also Soler et al. 2019

Mass [M ]
Figure 6.Same as Figure 5, but for total rather than line-of-sight magnetic field.

2020
).The HRO is a diagnostic derived from the machine vision technique of histograms of oriented gradients (Freeman & Roth 1994;Dalal & Triggs 2005).The original technique examined the distribution of orientations of gradients in two-dimensional images, while the HRO compares the orientation of a gradient to that of a vector field, possibly in three dimensions.We extend the standard analysis of the HRO method by working in three dimensions (position-position-position space), and by not only studying the relative orientation of the density gradient to the magnetic field, but also the orientation of the density gradient with respect to the velocity field.The density structure of the cloud can be characterized by the density gradient.Use of the gradient allows local comparison of the orientation of the magnetic field to iso-density contours.The HRO is computed by measuring the relative angle between the local density gradient ∇n and the magnetic field B as We can also measure the relative angle between the magnetic field and the velocity φ v B using the velocity vector in place of the density gradient in Eq. ( 17).Together with the total energy measurements and the characteristic velocities, this diagnostic can reveal whether gas flows are restricted to follow field lines, or whether mag-netic fields are just being advected by the turbulence.Values of φ v B consistently close to zero suggest that the velocities are constrained by the field, while a broad range of values suggests that field is being advected by the turbulence kinematically.The latter would also be likely if the kinetic energy exceeds the gravitational energy, and the rms velocity exceeds the Alfvén velocity.However then the scale of the turbulence determines the scale of magnetic field coherence.
Once the relative angles have been computed everywhere in the model, we collect this data into a histogram to determine if there are any preferred orientation angles.We proceed to bin the results with respect to the density, in order to explore the behavior of these relative orientations for a wide range of environments.Given that we perform our study of the HRO for physical scales over the range 0.1-100 pc in an AMR simulation, we do not force our density bins to have equal numbers of pixels, but rather force a minimum number of pixels per bin (> 10 3 ) to produce a well-sampled distribution of angles.
To characterize the shape of the HRO in each bin we employ the HRO shape parameter also used by Seifried et al. (2020) where A || is the area under the histogram for 1 > cos φ > 0.75 and A ⊥ is the area under the histogram for 0.25 > cos φ > 0. The relevant ranges are indicated in Figures 7-9 with blue and orange shading.A positive value of ξ corresponds to parallel dominance, while a negative value corresponds to a largely perpendicular orientation.We compute the variance in ξ within each density bin by taking the standard deviation σ || of the HRO values across the range of angles considered for A || , and a corresponding value for σ ⊥ .Then Finally, as we are interested in the dynamical evolution of the cloud as it collapses, we compute the time evolution of the HRO.

Orientation
We next use this characterization to describe the relative orientation of the magnetic field with respect to the density gradients in and around the cloud.Panels (b) of Figures 7-9 show the evolution over time of the HRO as a function of gas density.For n > 100 cm −3 the flow becomes super-Alfvénic as collapse proceeds.The relative orientation at these higher densities shows a transition from predominantly parallel with ζ nB > 0.5 to much more random or even perpendicular with ζ nB 0 and variance σ ζ > 1 as the field gets swept along by the flow.
This result agrees with previous applications of the HRO to simulated clouds by Soler et al. (2013), Soler & Hennebelle (2017), Seifried et al. (2020), and Girichidis (2021) using different approaches to cloud formation or structure.The general observation that fields typically shift from parallel to perpendicular during the collapse of dense structures has a long history (Heitsch et al. 2001;Ostriker et al. 2001;Li & Nakamura 2004;Nakamura & Li 2008;Collins et al. 2011;Hennebelle 2013;Soler et al. 2013;Chen & Ostriker 2015;Li et al. 2015;Chen et al. 2016;Zamora-Avilés et al. 2017;Mocz & Burkhart 2018;Chen et al. 2020;Seifried et al. 2020;Girichidis 2021).The critical density for the transition furthermore agrees with the density at which the field starts to increase with density (Figures 5 and 6), as would be expected from a collapse-dominated flow.

Orientation
In order to explain the flat behavior of the magnetic field-density relation for densities below n < 100 cm −3 , it has been suggested that

P(N)
A A HRO iso n- However configurations with the field perpendicular to the flow still constitute a substantial fraction of the volume, as even the lowest probability densities exceed 0.6, a quarter of the maximum value.The probability density of relative orientations is almost flat for number densities 0.5 < n < 50.This density range corresponds to transition gas in a thermally unstable phase of the ISM.Gas remains in this density range for only a short time, as it consists mostly of shocked gas, rapidly cooling towards the cold, dense phase of the ISM.
Large fluctuations in the relative distribution of the orientations are observed at t ∼ 4 Myr.This is likely caused by the explosion of a SN within the analyzed volume at t = 2.62 Myr, whose blast wave hits the cloud after a transit time of ∼ 1 Myr.The SN remnant expands into an inhomogeneous density distribution, producing turbulent motions both parallel and perpendicular to the magnetic field, resulting in the predominantly flat behavior of the probability seen.
The alignment between the velocity and the magnetic field might not only be caused by strong magnetic fields.As discussed by Padoan & Nordlund (1999) there are two opposite processes that can cause alignment: dynamical alignment, occurring when the field is strong enough to restrict the gas flows along field lines, and kinematic alignment, occurring when the field is swept up by the gas flows, forcing alignment.It is likely that for densities n < 100 cm −3 , the velocity and magnetic field are dynamically aligned, as the flows are trans-Alfvénic (see Sect. 3.1), at least when the system has not been recently perturbed by a nearby SN explosion.In this regime, the magnetic tension can restrict gas flows perpendicular to the magnetic field.
For densities n > 100 cm −3 , the angle φ v B shows large fluctuations as a function of time and density.This regime is less likely to be affected by random turbulence in the environment, but on the other hand is most likely to be affected by hierarchical gravitational contraction, as gravity is the dominant form of energy here (see Sect. 3.2).As the cloud contracts, it can compress the magnetic field.At t = 2 Myr, φ v B transitions at these densities from a preferentially aligned flow into a random distribution of the alignment.At t = 4 Myr, φ v B shows some alignment between the velocity and the magnetic field, which is most certainly caused by kinematic alignment of the collapsing gas, as the flow here is super-Alfvénic.Later, at t = 6 Myr, the relative distribution of the angles is very random, with more of the flow aligned with the field, but with large fractions of the gas having oblique relative orientations.This erratic behavior of the relative angles shows that when gravity controls the dynamics of the gas, the magnetic field is carried along with the contraction.

Caveats
Our results apply in detail only to the earliest evolution of dense clouds, as we have no star formation or feedback models in our simulations.As a result, we do not follow the ex-pected prompt cloud destruction from feedback that must occur for gravitationally collapsing clouds not to exceed the observed star formation rate, as we discussed in Sect. 1.
Even during the period of collapse, we are not following the detailed chemistry and ionization structure of the clouds but rather relying on a tabulated temperature-dependent cooling prescription to follow the thermodynamics of cloud formation.For comparison, Seifried et al. (2017) and Girichidis (2021) included chemistry, but not early stellar feedback, while Grudić et al. (2021) included both, reaching overall similar results on collapse and field structure.Our approach does allow us to focus in detail on the force structure and energetics of the collapsing cloud.
Our most significant numerical limitation is that we only resolve Jeans lengths with a minimum resolution of four cells, as expressed by the Truelove et al. (1997) criterion (see also Heitsch et al. 2001).This may be inadequate to resolve the turbulence and field evolution within the collapsing regions of the clouds (e.g.Sur et al. 2010;Federrath et al. 2011b), which could limit the applicability of our conclusions about the dynamics of the densest regions of the clouds.It is worth noting that Girichidis (2021) did a resolution study from 2 pc to 0.25 pc grid resolutions confirming convergence of the major results on field alignment and gravitational collapse.
The study of field angles has all been done for three-dimensional structures, where the dynamics remains unobscured by projection effects.These have been shown by Soler et al. (2019), Seifried et al. (2020), and Girichidis (2021) using detailed radiative transfer models to significantly impact the ability to deduce information about field structure from observations.

SUMMARY AND CONCLUSION
We have presented here three-dimensional MHD simulations of dense clouds formed in a supernova-driven, stratified, turbulent, galactic environment, collapsing under the action of their own self-gravity.These calculations were originally presented in Paper I and Paper II.We find that: • Dense clouds with n > 100 cm −3 have super-Alfvénic rms velocity dispersions, caused by fast turbulent motions driven by hierarchical gravitational contraction, suggesting that the fields are not dynam-ically important in this regime.In contrast, diffuse gas with n < 100 cm −3 in the envelopes of MCs and the diffuse ISM exhibits trans-Alfvénic velocity dispersion (Boulares & Cox 1990).This suggests the medium is magnetically supported against gravitational collapse (Elmegreen 2007), and its flow is generally constrained to follow field lines (Padoan & Nordlund 1999;Heiles & Troland 2005), although substantial deviations to this simple picture are seen in our simulations.
• Examination of the relative energy in the clouds shows that kinetic, magnetic and thermal energy are in equipartition in the diffuse medium with n < 1 cm −3 , kinetic energy dominates in the range of roughly 1 < n < 100 cm −3 , whereas gravitational energy dominates at densities n > 100 cm −3 .This is consistent with collapse dominated clouds.
• Direct analysis of the acceleration terms contributing to the momentum equation emphasizes that thermal pressure gradients and Lorentz forces are equally important in the diffuse medium, magnetic fields dominate at intermediate densities, but self-gravity determines the flow properties at high densities.
• Random turbulent motions in the diffuse ISM and gravitational contraction of dense structures maintain the magnetic field strength at magnitudes consistent with those reported in observations (Crutcher 2012).The behavior of the magnetic field strength in relation to density also appears generally consistent, with no correlation between the field strength and gas density for densities n < 100 cm −3 , and a positive correlation above that density.However, the average B-n relation in our model has a significantly shallower slope than the observed relation between the maximum line-of-sight field strength and the derived density.
• Three-dimensional HROs of magnetic field vs. density gradient in our models indicate the transition from parallel fields in the diffuse medium and unbound cloud envelopes to random orientationshowever, with some tendendy towards perpendicular configurations-in bound high density cloud cores.
• While gas flows are dynamically aligned along the magnetic field at low densities, oblique supersonic shocks produce strong density fluctuations that lead to significant volumes with the field more perpendicular to the density gradient than expected for pure dynamical alignment.Software: yt (Turk et al. 2011), Flash (Fryxell et al. 2000), numpy (Harris et al. 2020), matplotlib (Hunter 2007)

Figure 1 .
Figure 1.Three-dimensional rendering at t = 5 Myr of (Top to bottom) the galactic midplane for altitudes up to ±50 pc above and below the midplane (Paper I), close-up view of cloud M3 (Paper II), including magnetic field lines for a volume of 50 pc 3 , and a closer view of a collapsing core with the magnetic field lines crossing a plane with area 2 pc 2 centered in the core center of mass.

Figure 2 .
Figure2.Mean rms velocity v rms (red solid), sound speed c s (orange dashed), and Alfvén velocity v A (green dash-dotted) for clouds M3 (left), M4 (center), and M8 (right) as a function of gas density within a (100 pc) 3 box, centered in the cloud's center of mass.The average rotational velocity of the region with n > 100 cm −3 is also shown (blue dotted).The panels correspond to evolutionary times of 0, 2.5, and 4.8 Myr (from top to bottom) since the moment self-gravity was included.The lines show the volume-weighted medians, while the shaded areas show the central 50% of values in each density bin.

Figure 3 .
Figure3.Energy density in gravitational potential (blue), magnetic (green), thermal (dashed orange), and kinetic (red) forms as a function of density in the gas within a (100 pc) 3 box, centered on the cloud's center of mass for clouds M3 (left), M4 (center), and M8 (right) for three snapshots at t = 0, 2.4, and 4.8 Myr since the time self-gravity has been turned on.The lines show the volume-weighted medians, while the shaded areas show the central 50% of values in each density bin.
Figure 4.Mean accelerations from gravitational potential gradients (blue dash-dot), magnetic forces (green dash), and pressure gradients (orange solid) as a function of density in the gas within a (100 pc) 3 box, centered on the cloud's center of mass for clouds M3 (left), M4 (center), and M8 (right) for three snapshots at t = 0, 2.5, and 4.8 Myr since the time self-gravity has been turned on.Lines show median values in each density bin, while shading shows the central 50% of values in the bins.

Figure 5 .
Figure5.Line-of-sight magnetic field strength along the x-axis parallel to an arbitrary line of sight in the galactic midplane versus the gas density in a (100 pc) 3 volume centered in the cloud's center of mass for clouds M3 (left), M4 (center), and M8 (right).The panels show evolutionary times of 0, 2.9, and 4.8 Myr since the time self-gravity was included from top to bottom.A vertical dashed line marks the maximum resolution at which we resolve the Jeans length with at least four cells.Blue points show Zeeman observations of line-ofsight magnetic field strengths with their errors (see references inCrutcher et al. 2010b).The solid black line shows the most probable model of the maximum field strength along an observed line of sight, as a function of density(Crutcher et al. 2010b).The black dashed line with error bars gives the mean and dispersion of the simulated values; a linear fit gives the power-law index listed to the upper left.The green dashed line fits the maximum magnetic field at each density between the transition point to power-law behavior (n = 10 2 cm −3 ) and the Jeans resolution limit (n = 10 4 cm −3 ), which also yields a power-law index.

Figure 7 .Figure 8 .Figure 9 .
Figure7.Animations over the first 5 Myr after gravitational collapse for cloud M3 of (a) density slices with contours at 10 −2 , 1 and 100 cm −3 , and probability density functions of the relative orientations between the magnetic field B and (b) the local density gradient, ∇n and (c) the velocity v, as a function of the given gas density ranges, within a (100 pc) 3 box centered on the cloud's center of mass.The histogram shape parameter ξ given by Equation (18) with error bars from Equation (19) is shown as a function of density within narrower bins for the distribution.The regions contributing to ξ are colored for clarity: parallel dominance in blue and perpendicular dominance in orange.The still figure published is at the final time for the each model.The animation shows the evolution of the field in the highest density gas from parallel to neutral or perpendicular orientation to the density gradient as collapse proceeds.

Figure 10 .
Figure10.Probability density function of the relative orientations between the local velocity, v, and the magnetic field, B, as a function of gas density, within a (100 pc) 3 box centered on the center of mass of cloud M3.The panels correspond to evolutionary times of 0, 2, 4, and 6 Myr since the moment self-gravity was included.(At the earliest time there is no high density gas, so the final bin covers a wide range of densities.)