The structure and statistics of interstellar turbulence

The following system of compressible MHD conservation laws for the mass, momentum, and energy complemented by the ideal induction equation is solved in a cubic domain of linear size L with triply periodic boundary conditions. The equations include external random forcing terms in the rhs of (1b) and (1d), as well as the ISM-specific generalized volumetric cooling function in the rhs of (1d):

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{u}})=0,\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&{\partial }_{t}(\rho {\boldsymbol{u}})+{\boldsymbol{\nabla }}\cdot \left[\rho {\boldsymbol{u}}{\boldsymbol{u}}-{\boldsymbol{B}}{\boldsymbol{B}}+\left(p+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{2}\right){\boldsymbol{I}}\right]={\boldsymbol{f}},\end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{B}}+{\boldsymbol{\nabla }}\cdot ({\boldsymbol{u}}{\boldsymbol{B}}-{\boldsymbol{B}}{\boldsymbol{u}})=0,\end{eqnarray} \tag{ 1c }$$

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal E }+{\boldsymbol{\nabla }}\cdot \left[\left({ \mathcal E }+p+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{2}\right){\boldsymbol{u}}-({\boldsymbol{B}}\cdot {\boldsymbol{u}}){\boldsymbol{B}}\right]={\boldsymbol{u}}\cdot {\boldsymbol{f}}-\rho { \mathcal L }(\rho ,T).\end{eqnarray} \tag{ 1d }$$

Here ρ is the gas density, ${\boldsymbol{u}}$—velocity, pressure is given by the ideal gas equation of state $p\equiv (\gamma -1)e\rho$, where $\gamma =5/3$ is the specific heats ratio and e—the specific internal energy density. The total energy density includes kinetic, internal and magnetic components ${ \mathcal E }=\rho {u}^{2}/2+\rho e+{B}^{2}/2$ and $\langle { \mathcal E }\rangle =K+U+M$ is one of the invariants of an ideal system (with zero rhs). The generalized heat-loss function $\rho { \mathcal L }={n}^{2}{\rm{\Lambda }}(T)-n{\rm{\Gamma }}$, where $n=\rho /{m}_{{\rm{H}}}$ is the Hi number density, is introduced to mimic uniform photoelectric and cosmic ray heating rates and local radiative cooling of the gas under typical interstellar conditions [7]. The magnetic field strength, ${\boldsymbol{B}}$, is subject to the usual solenoidal constraint ${\boldsymbol{\nabla }}\,{\boldsymbol{\cdot }}\,{\boldsymbol{B}}=0$. The following initial conditions: ${\rho }_{0}+\delta \rho$ (where $\delta \rho$ represent 20% random density perturbations), ${p}_{0}$, ${{\boldsymbol{u}}}_{0}=0$, and ${{\boldsymbol{B}}}_{0}=({B}_{0},0,0)$ are applied. Random forcing implementation assumes zero momentum input ${\boldsymbol{f}}\equiv \rho {\boldsymbol{a}}-\langle \rho {\boldsymbol{a}}\rangle$ and applies a fixed-in-time large-scale solenoidal non-helical acceleration ${\boldsymbol{a}}(x)$ with proper normalization, ensuing a constant kinetic energy input rate. The forcing is turned on in the middle of the run at $t=18.9\,\mathrm{Myr}$, when a finite initial velocity boost $\delta {\boldsymbol{u}}=\tau {\boldsymbol{a}}$ is applied. Here τ is a coefficient setting the finite amplitude of the one-time velocity perturbation.

System (1a)–(1d) is integrated numerically using the piecewise-parabolic method on a local stencil (PPML [20, 21]), following the traditional implicit large eddy (ILES) approach [22]. This implies that only an effective Reynolds number ${Re}\,=\,{u}_{\mathrm{rms}}L/{\nu }_{\mathrm{eff}}$, Prandtl number ${\Pr }={c}_{{\rm{p}}}{\rho }_{0}{\nu }_{\mathrm{eff}}/{\kappa }_{\mathrm{eff}}$ and magnetic Prandtl number ${Pm}={\nu }_{\mathrm{eff}}/{\eta }_{\mathrm{eff}}$ can be defined (${\nu }_{\mathrm{eff}}$, ${\kappa }_{\mathrm{eff}}$, and ${\eta }_{\mathrm{eff}}$ are the effective kinematic viscosity, thermal conductivity, and magnetic diffusivity controlled by numerical dissipation). In the molecular ISM, ${Re}\,\gg \,1$ and ${Pm}\gg 1$ [21], while in the ILES there is no explicit control of dissipative processes and the magnetic Prandtl number is usually set by the numerical method at ${Pm}\sim 1$. The effective Reynolds number is ultimately controlled by the nature of the numerical method and by the grid resolution N, with higher order accurate methods and larger grids delivering higher Re. Thermal conduction can stabilize thermal instability (TI) on scales below the so-called Field length ${\lambda }_{{\rm{F}}}=\sqrt{\kappa T/{\rho }^{2}{\rm{\Lambda }}(T)}$ [23] and should be ultimately included in the models. Our simulations, however, do not sufficiently resolve ${\lambda }_{{\rm{F}}}$ [24] that varies somewhere between 0.1 and 0.001 pc, depending on local physical conditions. With such limited resolution, the heat transport will be dominated by turbulent diffusion.

Besides the effective Re, Pr, and Pm, two other nondimensional numbers are important for the considered system: the sonic Mach number ${{ \mathcal M }}_{{\rm{s}}}=\sqrt{\langle {(u/{c}_{{\rm{s}}})}^{2}\rangle }$, where ${c}_{{\rm{s}}}$ is the sound speed, controls the degree of compressibility; and the Alfvén Mach number ${{ \mathcal M }}_{{\rm{a}}}=\sqrt{\langle {(u/{v}_{{\rm{a}}})}^{2}\rangle }=\sqrt{\langle \rho {u}^{2}/{B}^{2}\rangle }$, where ${v}_{{\rm{a}}}=\sqrt{{B}^{2}/\rho }$ is the Alfvén speed, controls the mean kinetic-to-magnetic energy density ratio.

Each model is defined by five control parameters: the domain size L; the mean gas number density in the box n₀; the rms velocity ${u}_{0,\mathrm{rms}}$ associated with the driving force amplitude; the uniform magnetic field strength B₀; and the grid resolution N³.

Table 1 provides a summary of parameters for the five cases discussed below. All models assume the box size L = 200 pc, representing the scale height of the Galactic neutral Hi. Our choices for n₀ in the range from 2 to 5 cm⁻³ reflect local average conditions in a spiral arm-like setting with Hi number densities somewhat higher than the typical mean value of 1 cm⁻³. Likewise, the values of ${u}_{0,\mathrm{rms}}$ are chosen to bracket the realistic local ISM kinematics at the given scale L, matching the linewidth-size relation for local molecular clouds [25].

Table 1.  Input parameters.

Case	N	n0	${u}_{\mathrm{rms},0}$	B0	${\beta }_{\mathrm{th},0}$	${\beta }_{\mathrm{turb},0}$	${{ \mathcal M }}_{{\rm{a}},0}$
		cm−3	km s−1	μG
A	512	5	16	9.54	0.2	3.3	1.3
B	512	5	16	3.02	2	33	4.0
C	512	5	16	0.95	20	330	13
D	256	2	16	3.02	2	13.2	2.6
E	256	5	7	3.02	2	8.3	2.0

The table also provides the magnetic field parametrization in terms of the plasma beta ${\beta }_{\mathrm{th},0}={p}_{0}/({B}_{0}^{2}/2)$, representing the ratio of the initial thermal pressure of the gas p₀ to the magnetic pressure due to the applied uniform field B₀. The values of turbulent plasma beta (defined with the dynamic pressure $\rho {u}^{2}$ replacing thermal p in the classical definition) ${\beta }_{\mathrm{turb},0}={\rho }_{0}{u}_{\mathrm{rms},0}^{2}/({B}_{0}^{2}/2)$ and the initial Alfvénic Mach number ${{ \mathcal M }}_{{\rm{a}},0}={u}_{\mathrm{rms},0}/({B}_{0}/\sqrt{{\rho }_{0}})$ are also given for convenience. Overall, cases A–C trace a transition from trans-Alfvénic (${{ \mathcal M }}_{{\rm{a}},0}\sim 1$) to super-Alfvénic (${{ \mathcal M }}_{{\rm{a}},0}=4,13$) magnetization levels, while cases D and E are only moderately super-Alfvénic (${{ \mathcal M }}_{{\rm{a}},0}\sim 2\mbox{--}3$). Finally, the total gas mass in the domain for cases A–C and E is ${\mathfrak{M}}=9.8\times {10}^{5}\,{{\mathfrak{M}}}_{\odot }$, while in case D the mass is 2.5 times smaller, ${\mathfrak{M}}=3.9\times {10}^{5}\,{{\mathfrak{M}}}_{\odot }$.

For the sake of simplicity, we adopt an analytical approximation for ${\rm{\Lambda }}(T)$ [26, 27] based on a calculation [28] of the thermal equilibrium gas temperature of the diffuse interstellar medium that includes photoelectric heating rate from small grains and polycyclic aromatic hydrocarbons and a detailed treatment of the ionization rates and heating due to the soft x-ray background and due to cosmic rays. The cooling function approximation assumes solar metallicity. Radiative cooling is turned off for $T\lt 26$ K, however there is no temperature floor and the gas temperature drops to ∼4 K in strong rarefactions adiabatically. The density-independent heating rate is set at ${\rm{\Gamma }}=2\times {10}^{-26}$ erg s⁻¹. The high temperatures extend beyond 32 000 K in low-density gas, but the volume fraction of this 'overheated' material is very low and we ignore changes in the molecular weight in the equation of state due to ionization processes. With these assumptions, the S-shaped thermal equilibrium ${\rm{\Gamma }}=(\rho /{m}_{{\rm{H}}}){\rm{\Lambda }}(T)$ in the pressure–density plane has two extrema separating the two stable branches: a maximum at ${T}_{\max }=5250$ K and a minimum at ${T}_{\min }=184$ K (figure 1). These two temperatures set the upper and lower limits to the TI range of the linear isobaric mode, assuming equilibrium conditions [23]. Even though TI away from the equilibrium (i.e. where ${ \mathcal L }(\rho ,T)\ne 0$) is determined by a generalized criterion [29], we take a simplified approach and define the warm and cold stable phases solely based on the temperature. The warm phase thus includes all gas with $T\gt {T}_{\max }$ and the cold phase assumes $T\lt {T}_{\min }$. The range of intermediate temperatures ${T}_{\min }\lt T\lt {T}_{\max }$ roughly captures the thermally unstable regimes. With these definitions, our warm phase includes both warm neutral medium (WNM) and warm ionized medium (WIM).

Zoom In Zoom Out Reset image size

All cases were evolved in two stages. First, we allowed TI to develop in response to the small initial density perturbations, so that the system undergoes a phase transition and forms cold dense clouds embedded in warm medium, generating weak subsonic turbulence [30]. Then, as the thermal phases equilibrate and the system transitions into a stationary state, we turn on the random forcing at t = 19 Myr. The force is normalized to reach the specific target velocity dispersion values ${u}_{\mathrm{rms},0}$ assigned to each case, see table 1.

Figure 2 illustrates the evolution of mean energy densities in the code units for all cases. The kinetic energy, $K=\langle \rho {{\boldsymbol{u}}}^{2}/2\rangle$, shown in green indicates modest growth followed by saturation after the phase transition occurs. After that, a jump due to the activation of the forcing leads to a final turbulence saturation. The magnetic energy, $M=\langle {{\boldsymbol{B}}}^{2}/2\rangle$ (blue), starts at a level determined by the mean field B₀, which remains essentially unchanged during the phase transition stage. The magnetic energy then grows and reaches equipartition $M\sim K$ in cases A, D, and E, while remaining subdominant in cases B and C. Case D is similar to B, but has a lower kinetic energy saturation level due to $2.5\times$ lower mean density, which also helps to establish equipartition $M\sim K$. Case E is also similar to B, but has weaker forcing, causing lower energy saturation levels and equipartition $M\approx K$. The internal energy, $U=\langle \rho e\rangle$ (brown), first drops on a time scale of a few Myr due to radiative losses, then remains stationary during the post-phase-transition stage, then gets a boost from turbulence dissipation after the forcing is activated, and then finally saturates. The total energy, $E=\langle { \mathcal E }\rangle =K+U+M$, is shown in red. It is initially ($t\lt 19$ Myr) dominated by M in case A, by $M\sim U$ in cases B and E, and by U in case C. In fully developed driven turbulence, however, the internal energy is subdominant in all cases, so that $U\lt M\lesssim K$.

Zoom In Zoom Out Reset image size

The particular way turbulence is initialized in the system does not really matter, as a test case (not shown) with forcing activated at t = 0 resulted in the same stationary turbulence after a different initial transient. The model thus keeps no memory of the initialization process and the resulting turbulence and phase structure do not depend on initial conditions, besides the generalized cooling function and global parameters listed in table 1.

The details of magnetic field growth are illustrated in figure 3, where thick and thin lines show the total field ${B}_{\mathrm{rms}}$ and fluctuations ${b}_{\mathrm{rms}}$ for each case. The magnetic field is not directly forced in these models, but it can receive energy via random compressions and stretching of the field lines through interaction with the velocity field. The nonhelical random force we use does not generate a mean field, but still leads to amplification of the small-scale magnetic energy. The mean field strength, B₀, together with the mean density, n₀, and the rms velocity, ${u}_{\mathrm{rms},0}$, control the level of magnetic fluctuations in the developed turbulence. When B₀ is sufficiently strong, the saturation would imply equipartition of the kinetic and magnetic energy—the most interesting and realistic situation to explore. In this case, one can identify the uniform field, B₀, with the regular field generated by the Galactic dynamo and spatially ordered on scales $\gt 100$ pc [31, 32]. Since for the Milky Way-like galaxy the observed ratio of the small-scale and large-scale field strengths, $1.7\leqslant {b}_{\mathrm{rms}}/{B}_{0}\leqslant 1.9$ [33], this case would fall in the range covered by cases A, D, and E (table 2). It is worth noting that in the Galactic dynamo models, the origin of the large-scale field is associated with the interplay of small-scale turbulence generated by the star formation processes [34], differential rotation of the disk, and the inverse cascade of small-scale helicity [35–37]. In our simple local ISM model, the small-scale field is instead reverse-engineered from the external mean field using external forcing. As soon as the right ratio ${b}_{\mathrm{rms}}/{B}_{0}$ is achieved, local magnetic field properties and cloud magnetization levels can be considered as realistic.

Zoom In Zoom Out Reset image size

Table 2.  Magnetic field properties and compressibility in developed turbulence.

Case	${B}_{\mathrm{rms}}$	${b}_{\mathrm{rms}}$	${b}_{\mathrm{rms}}/{B}_{0}$	$\langle {\beta }_{\mathrm{th}}\rangle$	$\langle {\beta }_{\mathrm{turb}}\rangle$	${n}_{\mathrm{rms}}/{n}_{0}$	${u}_{\mathrm{rms}}$	${v}_{{\rm{a}},\mathrm{rms}}$	${{ \mathcal M }}_{{\rm{s}},{\rm{v}}}$	${{ \mathcal M }}_{{\rm{s}},{\rm{m}}}$	${{ \mathcal M }}_{{\rm{a}}}$
	μG	μG					km s−1	km s−1
A	16.6	13.6	1.4	0.10	2.0	3.7	15.5	15.4	4.9	11.6	1.0
B	11.7	11.3	3.7	0.36	4.0	3.5	15.4	12.6	5.4	13.1	1.4
C	7.1	7.0	7.3	1.76	8.0	3.5	16.0	8.8	5.9	13.2	2.1
D	8.2	7.6	2.5	0.77	2.7	2.8	15.8	12.2	3.1	7.5	1.6
E	7.2	6.5	2.2	0.37	2.2	3.4	7.4	10.2	2.8	6.4	1.3

Table 2 further details magnetic field properties in quasi-stationary developed turbulence, providing average rms field values and plasma beta regimes for the whole domain. The mean ${\beta }_{\mathrm{th}}$ is mostly below unity, except in case C, indicating subdominant contribution of the internal energy compared to the magnetic one. The mean ${\beta }_{\mathrm{turb}}$ is above 2 in all cases, meaning that the dynamic pressure always dominates over magnetic pressure on average. The table also presents diagnostics further characterizing the overall statistical properties of turbulence in numerical models. The relative rms density fluctuations ${n}_{\mathrm{rms}}/{n}_{0}$ in all cases fall within a narrow interval from 3 to 4, indicating very strong compressibility supported by TI and high volume- and mass-weighted sonic Mach numbers ${{ \mathcal M }}_{{\rm{s}},{\rm{v}}}$ and ${{ \mathcal M }}_{{\rm{s}},{\rm{m}}}$. The rms Alfvén speeds are comparable to the actually measured rms velocities in cases with strong magnetization, resulting in $K\sim M$ energy equipartition and ${{ \mathcal M }}_{{\rm{a}}}\sim 1$. Weakly magnetized cases demonstrate slightly super-Alfvénic rms velocities. Volume-weighted sonic Mach numbers (mostly sensitive to the warmer temperature regimes) fall roughly between 4 and 6, while mass-weighted ones (more biased to the cold phase) are largely supersonic, $6\lt {{ \mathcal M }}_{{\rm{s}}}\lt 13$.

Scatter plots of the gas density versus thermal pressure (the so-called phase diagrams or 2D (p, n)-PDFs) illustrate imprints of the cooling and heating processes on the overall dynamics of multiphase turbulence. Figure 4 shows volume-weighted (left) and mass-weighted (right) phase diagrams for data snapshots taken at $t=42.3\,\mathrm{Myr}$ in cases B, D, and E. As in figure 1, thick black lines show thermal equilibrium and dotted lines are isotherms at given T-values. The contours indicate constant levels of the volume (mass) fraction for different regimes of the thermal pressure, p, and density, n, separated by factors of 2. An orange circle corresponds to initial conditions $({p}_{0},{n}_{0})$ for each case.

Zoom In Zoom Out Reset image size

The volume-weighted diagrams show that the warm phase together with thermally unstable gas have large filling factors, while the cold phase occupies a tiny fraction of the domain volume. The mass-weighted diagrams instead highlight the distribution of the mass between the warm/unstable gas and the cold phase (table 3).

Table 3.  Mass and volume fractions, Mach number regimes for different thermal phases.

Case	Fv			Fm			${{ \mathcal M }}_{{\rm{s}}}$			${{ \mathcal M }}_{{\rm{a}}}$
	W	U	C	W	U	C	W	U	C	W	U	C
A	25	68	7	5	44	51	1.8	4.0	13.5	0.6	0.9	2.5
B	23	70	7	5	44	51	1.7	4.2	15.2	1.2	1.6	4.3
C	16	76	8	4	44	52	1.6	4.6	15.3	2.9	3.9	8.3
D	35	63	2	14	63	23	1.6	3.1	12.5	1.2	1.6	4.7
E	13	77	10	1	35	64	0.8	2.1	6.6	0.6	1.0	3.0

F_v—volume fraction (%); F_m—mass fraction (%); ${{ \mathcal M }}_{{\rm{s}}}$ and ${{ \mathcal M }}_{{\rm{a}}}$—volume-weighted rms sonic and Alfvén Mach numbers, respectively; W, U, and C label warm, thermally unstable, and cold regimes, respectively.

The phase diagrams for case B show that turbulence supports a wide range of thermal pressures and also that p in the molecular gas ($n\gt 100$ cm⁻³) is higher than that in the diffuse ISM, even though self-gravity is ignored in the model. Most of the volume is filled with the warm gas, while the mass is roughly equally distributed between warm (including considerable abundance of thermally unstable component [30, 39–41]) and cold components. Volume (${F}_{{\rm{v}}}$) and mass (${F}_{{\rm{m}}}$) fractions of thermal phases in simulations given in table 3 are generally consistent with locally observed values [42], e.g. $\gtrsim 48$% by mass of the (WNM, $\gt 500$ K) lies at temperatures that belong to the thermally unstable regime; about 60% of all Hi is WNM, which fills very roughly 50% of the volume.

Phase diagrams for cases D and E illustrate the dependence of the resulting multi-phase structure on the mean density of the gas n₀ and on the rms velocity ${u}_{\mathrm{rms}}$ (table 2). With a $2.5\times$ smaller density, case D retains less than a half of the cold phase mass fraction, compared to case B. The maximum pressure and the maximum density of the cold phase also drop. Note that the rate of kinetic energy supply required to support ${u}_{\mathrm{rms}}$ at the same level as in case B, is lower because the mean density is lower. Hence case D illustrates the combined effects of a reduced density and an effectively $\sim 2\times$ weaker forcing in a low-density environment. Case E has the same n₀ as B, but $2\times$ smaller ${u}_{\mathrm{rms}}$, which corresponds to $\sim 4\times$ smaller energy injection rate. Since the mean density remains moderately high, but the turbulence is weaker, the cold phase mass fraction is higher compared to B, but the range of thermal pressures visibly shrinks.

The effects of magnetization can be studied using the A-B-C sequence of cases, which have the same n₀ and ${u}_{\mathrm{rms}}$, but different B₀. However, as far as the ${p}_{\mathrm{th}}-n$ phase diagrams are concerned, the differences are quite subtle. Hence we defer the discussion to section 3.5 with a primary focus on control of star formation.

Table 3 also lists average Mach numbers for the cold, warm, and unstable phases. The sonic Mach number ${{ \mathcal M }}_{{\rm{s}}}$ (a measure of compressibility) is transonic in the warm phase (weakly compressible), is supersonic in the unstable regime (mild-to-strong compressibility) and is highly supersonic in the cold phase (very strong compressibility). Hence, the velocity scaling is expected to be Kolmogorov-like in the warm gas and Burgers-like in the cold phase (section 3.11). The ${{ \mathcal M }}_{{\rm{s}}}-n$ correlation, shown in figure 5 (left), follows the prediction ${{ \mathcal M }}_{{\rm{s}}}\propto {n}^{1/2}$ for isobaric conditions at intermediate densities $n\in [1,10]$ cm⁻³ [43] reasonably well and flattens into ${{ \mathcal M }}_{{\rm{s}}}\approx {\rm{const}}$ in quasi-isothermal conditions at low and high density extremes (section 3.6).

Zoom In Zoom Out Reset image size

The Alfvén Mach number ${{ \mathcal M }}_{{\rm{a}}}$ (a measure of relative strength of magnetic and kinetic effects) is mostly trans-Alfvénic in the warm phase, except in case C, indicating approximate energy equipartition $K\sim M$. The unstable phase is transient from trans-Alfvénic to super-Alfvénic ($K\gtrsim M$), and the cold phase is mildly (cases A and D) or strongly (case C) super-Alfvénic ($K\gg M$). This creates a curious physical situation, where nonlinear relaxation of the system results in a statistically stationary state with moderately magnetized turbulence across most of the volume, enforcing overall kinetic-to-magnetic energy equipartition. At the same time, the cold, high density material potentially available for star formation remains in supersonic, super-Alfvénic regime with kinetic energy dominating magnetic by a large margin. Figure 5 illustrates a clear positive correlation between the Alfvén Mach number and density for case B with the bulk of the gas mass at $n\gt {10}^{2}$ cm⁻³ having local ${{ \mathcal M }}_{{\rm{a}}}\gt 1$. A similar trend is seen in isothermal MHD simulations [44] in supersonic and weakly super-Alfvénic regimes, which can be explained by the dependence ${{ \mathcal M }}_{{\rm{a}}}(B,\rho )$ and the relationship between B and ρ. Our multiphase cases combine a transonic, sub-Alfvénic behavior at low densities with a supersonic and super-Alfvénic one at high densities.

Our simulations do not include gravitational effects and all cases represent globally Jeans-stable settings. We can, however, extrapolate our results to cases with self-gravity by picking a sample of evolved snapshots with stationary statistics and asking a simple question: what happens if we turn on self-gravity? To address this, we have to return to the discussion of phase diagrams in figure 4, where each plot contains a dashed line in the upper-right corner, separating locally Jeans-unstable material.

All gas to the right of the dashed line ($J\equiv {\rm{\Delta }}x/{\lambda }_{{\rm{J}}}=0.25$, where ${\lambda }_{{\rm{J}}}=\sqrt{\pi {c}_{{\rm{s}}}^{2}/G\rho }$ is the Jeans length) would fail the Truelove resolution constraint [38], which requires the Jeans length to be resolved with at least four grid zones to avoid artificial fragmentation. If we were to follow local gravitational collapses with adaptive mesh refinement in these simulations, the grid zones failing Truelove's condition would be flagged. The Truelove condition is also used to determine whether sink particles should be introduced in models which rely on this technology. In practice, the condition implies that (with some efficiency factor $\lesssim 1$) this material would be converted into stars within one free-fall time ${t}_{\mathrm{ff}}=\sqrt{3\pi /32G\rho }\approx \sqrt{({10}^{3}\,{\mathrm{cm}}^{-3})/n}$ Myr, where G is the gravitational constant [45].

This Jeans-unstable gas occupies a tiny fraction (0.01%) of the computational volume and its mass ${\rm{\Delta }}{\mathfrak{M}}$ comprises 1.7% of the mass in case B at $t=42.3\,\mathrm{Myr}$. The cold phase in this case constitutes on average 7% of the volume and 51% of the mass (table 3). In this case, the star formation rate per free-fall time [13, 46] ${\epsilon }_{\mathrm{ff}}\equiv {t}_{\mathrm{ff}}/{t}_{\mathrm{dep}}={\rm{\Delta }}{\mathfrak{M}}/{\mathfrak{M}}\lesssim 0.017$, where the gas depletion time ${t}_{\mathrm{dep}}\equiv {\mathfrak{M}}/\dot{{\mathfrak{M}}}$, is consistent with the recent Spitzer legacy survey of low-mass star-forming clouds near the Sun, which gives ${\epsilon }_{\mathrm{ff}}\sim 0.01-0.1$ for clouds with mean densities ${n}_{{{\rm{H}}}_{2}}\sim {10}^{3}$ cm⁻³ [47].

In case D with $2.5\times$ lower mean density compared to B (and proportionally lower kinetic energy injection at large scales), the cold phase constitutes only 2% by volume and 23% by mass. In this case, star formation is very inefficient. It is active in just 0.002% of the domain volume and involves 0.22% of the mass per ${t}_{\mathrm{ff}}$. In case E with a $2\times$ lower rms velocity, the cold phase constitutes 10% of the volume, 64% of the mass, and the upper limit for the star formation rate is set at 7.1% of the mass per ${t}_{\mathrm{ff}}$ in 0.14% of the volume. In this case, relatively low level of turbulence results in a high star formation rate. Note that cases D and E have a $2\times$ lower grid resolution compared to B, which could to some degree affect the phase fractions. The estimates of ${\epsilon }_{\mathrm{ff}}$ include the resolution dependence, which enters the Truelove constraint [38].

We have traced the dependence of ${\epsilon }_{\mathrm{ff}}$ on the mean density and rms velocity in the box, but more subtle magnetic effects remained unexplored. In order to see how magnetization changes the expected star formation efficiency, we calculated the time-averaged star formation rates per free-fall time for cases A–C, which have the same n₀ and ${u}_{\mathrm{rms}}$, but different B₀. The results: $\langle {\epsilon }_{\mathrm{ff}}\rangle =0.021\pm 0.005$ (A), 0.013 ± 0.002 (B), and 0.011 ± 0.002 (C) indicate an overall $\sim 2\times$ drop in $\langle {\epsilon }_{\mathrm{ff}}\rangle$ from A to C. These estimates, however, are based on the purely thermal Truelove condition [38], which does not fully take into account the magnetic effects even though the pressure and density fields that enter the Truelove condition are affected. Since strong magnetization is able to suppress fragmentation, higher densities would be required to trigger the collapse, relaxing the critical J value ${J}_{{\rm{M}}}=J/\sqrt{1+0.74/{\beta }_{\mathrm{th}}}$ [48]. If we apply this ${\beta }_{\mathrm{th}}$-dependent correction, the average rates drop by a factor of ∼10 in all three cases, but the statistics become poor because the maximum density is limited by the grid resolution. Thus, higher resolution modeling is required to measure how magnetization changes the expected star formation efficiency.

We note that the star formation rates discussed here differ from those in the context of idealized isothermal simulations of MC turbulence, where ${\epsilon }_{\mathrm{ff}}$ is determined primarily by the virial parameter of an individual cloud [49]. These rates may vary in a wide range of values, depending on cloud's dynamical state [50]. Since our multiphase models contain large statistical ensembles of realistic clouds, ${\epsilon }_{\mathrm{ff}}$ measures ensemble-average rates. Naturally, the scatter between different snapshots from a given run under statistically stationary conditions is reasonably small, with standard deviations of order $(15 \% \mbox{--}20 \% )$. A better understanding of the functional form of ${\epsilon }_{\mathrm{ff}}(L,{n}_{0},{u}_{\mathrm{rms}},{B}_{0})$ will help to design robust subgrid-scale prescriptions for star formation and feedback in galaxy formation simulations [51, 52].

Thus the two-phase turbulence has a potential to regulate the supply chain of star formation, permitting only so much cold gas, dripping through the upper-right 'nose' of the phase diagram per ${t}_{\mathrm{ff}}$, to undergo gravitational phase transition. Star formation rates per free-fall time ${\epsilon }_{\mathrm{ff}}\sim (0.001\mbox{--}0.1)$ we measured in cases B, E, and D are controlled by the problem parameters $(L,{n}_{0},{u}_{\mathrm{rms},0},{B}_{0},{\rm{\Gamma }})$, but this simple model does not explicitly include the feedback processes, nor does it include the mass source associated with the cosmological gas infall onto the Galaxy. Natural further steps in development of the model would include closing the energy feedback loop and allowing the mass flux through the system (with external gas accretion on the disk as a source and star formation as a sink). This would be an advanced variant of the galactic bathtub model [53]. Remarkably, ${\varepsilon }_{\mathrm{ff}}\sim 0.01$ that we get in case B is a reasonable estimate for the star formation rate on scales from entire galaxies to single clouds [13]. This suggests that closing the feedback loop might help to nail down a stable, statistically steady regime for the recycling of gas into stars by disk-like galaxies between major merger events.

The gas density PDF is one of the basic statistics of compressible turbulence that bears unique signatures of gas thermodynamics [54], gravitational instability [45], and the presence of shocks. For instance, in isothermal homogeneous supersonic turbulence, the PDF is lognormal [55–57]. In two-phase ISM turbulence with moderate velocity dispersions, however, the PDF is expected to have two peaks corresponding to thermally stable phases and a power-law excess at high densities because the cold phase thermal equilibrium can be approximated by an effective soft polytropic relation ${p}_{\mathrm{eq}}\propto {\rho }_{\mathrm{eq}}^{{\gamma }_{\mathrm{eff}}}$ with ${\gamma }_{\mathrm{eff}}\approx 0.7\lt 1$ [30, 58, 59]. Cases A through D assume a moderately strong turbulence driving at an injection scale ${\lambda }_{f}\sim 100$ pc, which effectively merges the usual double-peak distribution obtained at ${u}_{\mathrm{rms}}\sim 5\mbox{--}10$ km s⁻¹ into an asymmetric single-peak PDF. Density PDFs for cases A–C shown in figure 6 (left) display a minor tendency toward weaker rarefactions at higher magnetization levels moderated by stronger magnetic tension. Otherwise, the three PDFs are very similar and have characteristic lognormal shapes at the high end with no signature of any power-law extension. The lack of a thermal power-law tail is due to: (i) high levels of turbulence, supporting a wide range of thermal pressures and populating the unstable phase, which smear the soft effective polytropic law (figure 4) and (ii) flattening of the Mach number–density correlation ${{ \mathcal M }}_{{\rm{s}}}\propto {\rho }^{1/2}$, consistent with quasi-isobaric conditions [43], above $n\sim 30$ cm⁻³ (figure 5). It is ultimately the lack of ${{ \mathcal M }}_{{\rm{s}}}-n$ correlation that supports the lognormal shape of the PDF high end populated primarily with cold material under supersonic conditions. Note that the moderately strong forcing responsible for the lognormal shape here is not merely a random choice of a free model parameter, but is required to match the local observed distribution of thermal pressure (section 3.7) and local linewidth–size relationship (section 3.11).

Zoom In Zoom Out Reset image size

Figure 6 (right) shows PDFs of column densities (or N-PDFs) obtained by averaging projected density distributions in three coordinate directions. Interestingly, the convolutions involved in the projection operation completely remove the signature of two-phase density structure and return nearly perfect lognormal N-PDFs, particularly in case B, however with noticeable statistical variations at low column densities in cases A and C. Such lognormal N-PDFs should be observable with temperature-blind tracers in clouds with no signs of active star formation, where the effects of self-gravity are too weak to build a power-law tail at high column densities (e.g. the Draco cloud [60]). Note that the sonic Mach numbers for cases A–C are generally higher than 5 (table 2). At ${{ \mathcal M }}_{{\rm{s}},{\rm{v}}}\lesssim 2$, the lognormal property of the low tail of N-PDFs is lost [59].

Since temperature plays the role of an important third dimension in column-density mapping of star forming clouds [61], we proceed by computing N-PDFs conditioned on the temperature $\wp (\mathrm{log}N| {\rm{\Delta }}T)$, such that $\int \wp (\mathrm{log}N| {\rm{\Delta }}T){d}\,\mathrm{log}\,{\rm{N}}=1$ and ${\rm{\Delta }}T$ defines the allowed temperature interval. As an example, figure 7 shows a projected density image (left) and a synthetic map of the column density of the cold phase gas (right), showing a plethora of multi-scale filamentary structures created by turbulence with the aid of TI. Even though our models do not include any particular tracer's chemistry, simple conditioning on the temperature can yield new predictions for temperature-dependent N-PDF shapes. Figure 8 (left) shows case-B N-PDFs for the cold (purple, $4\lt T[{\rm{K}}]\lt 184$) and warm (red, $T\gt 5250$ K) phases, as well as the distribution including all temperature regimes (black solid line, $\forall \,T$). Three best-fit lognormal distributions are shown with black dashed and dotted lines together with their corresponding standard deviations: ${\sigma }_{{\rm{c}}}=0.29$, ${\sigma }_{{\rm{w}}}=0.17$, and ${\sigma }_{{\rm{a}}}=0.25$. Remarkably, both warm and cold phases bear extended low-end tails (that would be difficult measuring observationally), while the full N-PDF is purely lognormal. None of the lognormal distribution widths can be accurately predicted using the variance–Mach number relation ${\sigma }_{N}^{2}=0.11\mathrm{ln}(1+{b}^{2}{{ \mathcal M }}_{{\rm{s}}}^{2}$) based on simulations of isothermal turbulence [62], if we assume ${{ \mathcal M }}_{{\rm{s}}}$ values from table 3 and $b\approx 1/3$ due to a divergence-free nature of the forcing [63]. Similar discrepancy of the widths of fitted lognormal distributions with isothermal predictions was seen earlier in multiphase ISM simulations without a magnetic field [59], in GALFA-HI data for the Perseus molecular cloud [64] and in data for seven molecular clouds from the Leiden/Argentine/Bonn Galactic Hi survey [65]. In all cases, the widths of observational H_I N-PDFs are similar to our model predictions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Self-shielding of H₂ becomes effective around ${N}_{{\rm{H}}}\sim 2\times {10}^{20}$ cm⁻², while column densities in excess of ${N}_{{\rm{H}}}\sim 4\times {10}^{21}$ cm⁻² are required to form CO [67, 68]. The range of column densities in between is populated by the 'CO-dark H₂' gas that may account on average for $(30 \% \mbox{--}40 \% )$ of the Milky Way molecular mass [69, 70]. The shaded areas in both panels of figure 8 show the column density interval, where CO-dark H₂ is expected to reside. Figure 8 (right) details the temperature dependence of N-PDFs for the cold gas in case B. In all cases, except perhaps for the $T\lt 25$ K case that hits our resolution limits, the high tails of these distributions can be reasonably well represented as lognormal. Similar plots for cases A and C (not shown) display only subtle differences in the N-PDFs with case B, mostly limited to the low column density tails. This is consistent with figure 6 and indicates that corresponding observational diagnostics are not very sensitive to the ISM magnetization levels.

At high levels of turbulence, the distribution of thermal pressure for cases A–C spans about six decades leaving no room for the classical isobaric pressure-supported cloud picture in the violent ISM (figure 9, left). The PDFs are asymmetric in a very characteristic way, showing a positive skewness, reported in simulations earlier [71, 72–74]. The origin of this asymmetry can be traced back to the turbulence-regulated and radiative cooling-controlled distribution of thermal phases in the phase diagram (figure 4). The peak tends to shift slightly to higher pressures as the magnetization becomes stronger, approaching the mean ISM pressure of 3700 K cm⁻³ reported in [75] in the range of B₀ bracketed by cases A and B. The peak location is sensitive to pressure in the warm phase and in the thermally unstable regime, which together make up $\approx 93$% of the volume in cases A–C (table 3). Models with higher magnetization, however, tend to have a progressively more abundant warm phase and less populated thermally unstable regime. Since the stable warm phase has higher ${p}_{\mathrm{th}}$ on average, compared to the unstable gas (figure 4), the PDF's peak shifts to the right along the C-B-A sequence. The physics behind the shift apparently has to do with the stabilizing effect of the magnetic field on the TI [23]. Note that the overall shape of the PDF is not lognormal because at high turbulence levels each pressure bin includes contributions from a wide range of temperatures, while some parts of the high tail can still be approximated by a lognormal with reasonable accuracy.

Zoom In Zoom Out Reset image size

In order to compare the pressure PDF with observations that include line-of-sight convolutions, we use density-weighed PDFs compiled using sightlines parallel to the orthogonal coordinate directions of ∼70 data snapshots from each model. In addition, conditioning is used to match observations that can be blind to certain temperature regimes and are limited to sightlines below some extinction threshold. The synthetic PDFs shown are conditioned to mimic the HST archive data for the Solar Neighborhood obtained from fine structure Ci lines in UV stellar spectra sensitive to the cold neutral gas [75] ⁴ . We thus exclude lines of sight with $N({\rm{H}}{\rm{}}{\rm{I}})\gt 2.5\times {10}^{21}$ cm⁻² and mask out the gas with temperatures outside 40 K $\lt \,T\lt 183$ K. Besides these limits, there are no other free parameters involved, although we use arbitrary normalization vertically to fit the data. The effects of the column density and temperature thresholds can be seen in the synthetic map shown in figure 7.

Distributions for cases A through C with ${u}_{\mathrm{rms},0}=16$ km s⁻¹ match the characteristic mean pressure typical for the Milky Way disk at the solar radius, reproduce its overall asymmetry, and show only weak dependence on B₀. The width of the distribution, however, is sensitive to ${u}_{\mathrm{rms}}$ and n₀. Case D marginally reproduces the shape and the width of the observed distribution with some deficiency in the high tail. A lower turbulence level in case E yields too narrow a distribution, and hence should be rejected. We also formally applied the Kolmogorov–Smirnov (KS) test to compare the HST sample with our modeled cumulative distribution functions. Based on the KS test results, case B provides the best match to the HST data with the goodness of fit D = 0.066 and the significance level of this KS statistic $\wp =0.95$. Case B is only slightly better than A, but other cases returned poor results, as the rough visual comparison suggested. Overall, due to the conditioning by temperature and opacity, the mass-weighted PDF is substantially narrower than the actual PDF shown in the left panel. While both low and high ends are suppressed, the effects of conditioning are stronger at high thermal pressures.

These numerical experiments probe the levels of magnetic field strength in molecular clouds that form self-consistently in the magnetized turbulent diffuse ISM. Figure 10 shows scatter plots of the magnetic pressure ${p}_{\mathrm{mag}}\equiv {B}^{2}/2$ versus the thermal pressure ${p}_{\mathrm{th}}=(\gamma -1)\rho e$ (left) and versus the dynamic pressure ${p}_{\mathrm{dyn}}\equiv \rho {u}^{2}$ (right) for case B at 42.3 Myr. Since ${\beta }_{\mathrm{th}}\equiv {p}_{\mathrm{th}}/{p}_{\mathrm{mag}}=(\gamma -1)U/M$, ${\beta }_{\mathrm{turb}}\equiv {p}_{\mathrm{dyn}}/{p}_{\mathrm{mag}}=2K/M=2{{ \mathcal M }}_{{\rm{a}}}^{2}$, and ${\beta }_{\mathrm{turb}}/{\beta }_{\mathrm{th}}=\gamma {{ \mathcal M }}_{{\rm{s}}}^{2}=2/(\gamma -1)K/U$, one can use the two plasma beta parameters to access pairwise energy equipartition and Mach number regimes at different locations on these diagrams. For instance, the dashed line ${\beta }_{\mathrm{th}}=0.1$ in the left panel, crossing right through the 'summit' of the black contours, indicates that magnetic pressure is an order of magnitude higher than thermal in the bulk of the volume (same applies to magnetic and internal energies, $M\gg U$). Interestingly, the same remains true for the cold molecular gas ($T\lt 100$ K and $n\gt 100$ cm⁻³), as indicated by the color contour map in figure 10.

Zoom In Zoom Out Reset image size

The dynamic pressure, however, is of the same order as the magnetic one in the bulk of the volume, as the dash-dot line ${\beta }_{\mathrm{turb}}=1$ in the right panel shows. This also means that $K\sim M$ and ${{ \mathcal M }}_{{\rm{a}}}\sim 1$ in the bulk of the volume, consistent with data in table 3. It is worth noting, however, that these conditions do not extend to the cold molecular gas, as the color contours centering around ${\beta }_{\mathrm{turb}}=30$ indicate. The cold molecular gas in our simulations is instead super-Alfvénic with ${{ \mathcal M }}_{{\rm{a}}}\sim 8$ (see also figure 5, right) and $K\sim 10M$. Qualitatively, this holds even in the strongly magnetized case A.

A key to understanding the origin of this super-Alfvènic regime of strong MHD turbulence in the cold and dense molecular gas lies in the process of self-organization caused by stochastic flux freezing, which is in turn mediated by the reconnection diffusion process [76, 77]. Magnetic flux conservation in turbulent fluids at high magnetic Reynolds numbers is not valid in the classical sense, nor is it entirely broken. Instead it holds in a statistical sense associated with the 'spontaneous stochasticity' caused by turbulent Richardson-like diffusion of magnetic field lines, which leads to a breakdown of the classical flux freezing concept. This subject is still under development for compressible MHD turbulence and our simulations may not completely capture the small-scale dynamo effects in the cold dense gas due to limited resolution and because of the implicit numerical nature of dissipation resulting in ${Pm}\sim 1$ and uniform ${\nu }_{\mathrm{eff}}$ that do not properly approximate conditions in the multiphase ISM.

The statistical steady state explored by our simulations is characterized by a certain degree of alignment between the velocity and magnetic field lines. If molecular clouds form in the turbulent ISM via large-scale compression of the warm diffuse trans-Alfvénic Hi primarily along the field lines, then turbulence in the clouds will naturally be super-Alfvènic [78]. A closely related subject is the relative orientation of the magnetic field direction with respect to filamentary density structures in simulated turbulent molecular clouds [79, 80]. Observations of striations and filaments in local molecular clouds and their close environments [81–86] suggest that the filaments are preferentially oriented perpendicular to the magnetic field lines, while the low-density subcritical striations are parallel. Moreover, it is believed that the background cloud material is funneled along the magnetic field lines onto star-forming filaments through striations, consistent with the kinematic constraints from CO observations [87]. Here we check the likelihood of such a scenario by looking at the alignment of local magnetic field and velocity vectors in different density regimes.

Figure 11 shows the PDFs of the cosine of the magnetic field–velocity angle $\theta =\arccos ({\boldsymbol{B}}\cdot {\boldsymbol{u}}/\sqrt{{B}^{2}{u}^{2}})$ for cases A–C (left) and details the alignment properties as a function of density for case A (right). The ${\boldsymbol{B}}-{\boldsymbol{u}}$ alignment is most pronounced in the strongly magnetized case A and in the bulk of the volume at densities slightly below the mean n₀ (shown by the dash-dot line), where the turbulence is sub-Alfvénic. However, the alignment is very weak at both extreme ends of the density distribution: in rarefactions and in shocked material at very high densities (figure 11, right). Indeed, in the magnetically dominated (sub-Alfvńic) warm phase, the gas structure is expected to be preferentially aligned with the local magnetic field.

Zoom In Zoom Out Reset image size

For the dense cold phase with super-Alfvénic conditions, however, this dynamic alignment mechanism does not work, and the filament elongation becomes preferentially perpendicular to the local magnetic field. Case C, with supersonic and super-Alfvénic turbulence in all phases (table 3) shows no significant alignment between the velocity and magnetic field vectors. In this case, however, the kinematic alignment (between ${\boldsymbol{B}}$ and the eigenvectors ${\boldsymbol{\lambda }}$ of a symmetric part of the rate-of-strain tensor ${s}_{{ij}}=({\partial }_{j}{u}_{i}+{\partial }_{i}{u}_{j})/2$) is actually expected instead [88]. Filament orientation with respect to the magnetic field depends on the shear strain alignment with respect to the field lines. If the strain is strong and (anti)parallel to the magnetic field, the filament will be aligned with the field [89]. Gravitational effects can further complicate the alignment pattern in dense filaments through the amplification of the potential (curl-free) component of the gas velocity.

Most of what we know about the interstellar magnetic fields comes from observations that usually involve complex convolutions along the observed sightlines [90]. In order to decipher that convoluted information, additional assumptions must be made concerning possible symmetries involved (e.g., whether the turbulent fluctuations of the field are isotropic) and the shape of probability distributions of various magnetic field components (e.g., if they are Gaussian). Figure 12 shows PDFs of magnetic field fluctuations ${\boldsymbol{b}}={\boldsymbol{B}}-{{\boldsymbol{B}}}_{0}$ perpendicular (b_⊥) and parallel (b_∥) to the mean field ${{\boldsymbol{B}}}_{0}$ for cases A–C. All PDFs represent averages over 70–100 flow snapshots, uniformly covering the statistically stationary state. In the case of b_⊥, the results include both b_y and b_z components normal to ${{\boldsymbol{B}}}_{0}$ directed along the positive x-coordinate direction.

Zoom In Zoom Out Reset image size

The plots show that the PDFs are strongly non-Gaussian in all cases and that the B-field fluctuations in the direction of ${{\boldsymbol{B}}}_{0}$ are strongly asymmetric, i.e. they know fairly well about the direction of the mean field. In the strongly magnetized case A, the PDFs of perpendicular components show a plateau around zero with a width of order B₀. The case-A PDF is platykurtic with an excess kurtosis ${\gamma }_{2}\equiv {\mu }_{4}/{\sigma }^{4}-3=-0.4$. Moderately and weakly magnetized cases B and C instead show leptokurtic (i.e. more peaked than the normal distribution) PDFs with ${\gamma }_{2}=0.3$ and 3.6, respectively.

The PDF of the parallel component b_∥ is skewed to negative values with the skewness ${\gamma }_{1}\equiv {\mu }_{3}/{\sigma }^{3}=-1.0$ and −0.2 in cases A and B, respectively. Case C shows modest positive skewness ${\gamma }_{1}=0.2$. The models thus predict strongly non-Gaussian distributions for the magnetic field components both parallel and perpendicular to the mean field in cases with sufficiently strong magnetization.

Figure 13 shows PDFs of $| {\boldsymbol{B}}|$ for cases A–C; vertical black dashed lines indicate the values of B₀. All three PDFs have fat stretched-exponential tails extending to the field strengths exceeding dozens or even hundreds of the mean field value. Curiously, the likelihood to find extreme values in the range $B\in [100,120]$ μG in the weakly magnetized case C (green line) is approximately eight times higher than in the strongly magnetized case A (red line). The mean field B₀ thus controls the abundance of cold dense gas with extreme field values produced in strong intermittent (both in space and in time) 3D compressions which are naturally more likely under weak magnetization.

Zoom In Zoom Out Reset image size

Using the gas temperature and density as additional variables to condition the distributions, one can check which physical regimes contribute most to the fat tails or to the peak of the $| {\boldsymbol{B}}|$-PDFs. Figure 14 shows that in the strongly magnetized case A the extreme tail is due to the cold phase (left) and contains exclusively dense ($n\gt 30$ cm⁻³, most likely molecular) gas (right), whereas the peak of the distribution is mostly due to the gas at intermediate (thermally unstable) temperatures and at lower densities ($n\lt 2$ cm⁻³). The most likely rms B-field values for the cold phase (or for the dense gas) are only slightly higher than those for the full distribution, whereas the peak of the warm phase (or low density gas) is only slightly shifted to weaker field values. These distributions can predict general trends in expected magnetic field regimes accessible with observational tracers populating parts of molecular clouds with different temperatures and densities.

Zoom In Zoom Out Reset image size

Larson's first law [91, 92], also known as the linewidth–size relation, is an observable surrogate for first order (m = 1) structure functions of the velocity field that are defined by

$$\begin{eqnarray}&&{S}_{m}(r)=\langle {[u({\boldsymbol{x}}+{\boldsymbol{r}})-u({\boldsymbol{x}})]}^{m}\rangle \propto {r}^{{\zeta }_{m}},\end{eqnarray} \tag{ 2 }$$

where $\langle \cdot \rangle$ denotes the ensemble average over all possible point pairs separated by the lag $r=| {\boldsymbol{r}}|$. The velocity components here can either be parallel or perpendicular to vector ${\boldsymbol{r}}$, corresponding to longitudinal ${S}_{m,\parallel }(r)\propto {r}^{{\zeta }_{m}^{\parallel }}$ and transverse ${S}_{m,\perp }(r)\propto {r}^{{\zeta }_{m}^{\perp }}$ structure functions, respectively.

Figure 15 shows scaling results for ${S}_{1,\parallel }(r)$ and ${S}_{1,\perp }(r)$ based on 16 lag values in the range of $r\in [1.5,50]$ pc for cases A–C. Black triangles and black dashed lines show structure functions obtained from uniform sampling of the whole computational domain, where all temperature regimes are represented. In cases A and B with strong-medium magnetization, scaling exponents ${\zeta }_{1}^{\parallel }\approx {\zeta }_{1}^{\perp }=0.4$, while in weakly magnetized case C ${\zeta }_{1}^{\parallel }\approx {\zeta }_{1}^{\perp }=0.5$. The fact that case C shows an exponent close to that found in isothermal hydrodynamic turbulence simulations at high Mach numbers is not surprising since ${{ \mathcal M }}_{{\rm{s}}}\approx 13$ and ${{ \mathcal M }}_{{\rm{a}}}\approx 2$ (table 3), i.e. the turbulence is supersonic and super-Alfvénic on average in the whole volume. In all cases, the scaling range is about one dex, which implies the exponents are reasonably accurate.

Zoom In Zoom Out Reset image size

Observational measurements are usually done using a single molecular hydrogen tracer (e.g. CO) available only in certain density and temperature regimes. To mimic this situation, one can condition the ensemble averaging operator in (2) to only include point pairs in the density and temperature regimes appropriate for the given tracer. Technically, simulations with sufficient resolution will provide reasonably accurate statistics for the density and temperature regimes that are not overtly restrictive. This would help study trends in the tracer-specific bias in observed velocity statistics. For instance, figure 15 shows structure functions conditioned for the cold phase ($T\lt 184$ K, purple) and for the warm phase ($T\gt 5,250$ K, red). The warm phase shows Kolmogorov scaling ${\zeta }_{1}^{\parallel }={\zeta }_{1}^{\perp }\approx 1/3$ in all three cases. Indeed, the Mach number regimes for the warm phase are only mildly transonic (table 3) and thus imply ${S}_{1}(r)\propto {r}^{1/3}$ [93]. In contrast, the scaling of structure functions for the cold phase is steeper: ${\zeta }_{1}\approx 0.5$ in case A, 0.6 in case B, and 0.65 in case C. This is not surprising since ${{ \mathcal M }}_{{\rm{s}}}\approx 15$ for the cold phase (table 3) and ${\zeta }_{1}\sim 1/2$ are common for isothermal turbulence at ${{ \mathcal M }}_{{\rm{s}}}\sim 10$ [57]. Our simulations, thus, reproduce the observed linewidth–size relations for molecular clouds (including both the slope and the offset) within the uncertainty introduced by the projection effects.

We now turn to second order, two-point statistics and discuss the power spectra of primitive variables, such as density and velocities. The power spectrum is defined here in a usual way, including integration over spherical shells, e.g., for a scalar field $q({\boldsymbol{x}})$, the spectrum $P(q,k)\equiv \int | \widehat{q}({\boldsymbol{\kappa }}){| }^{2}\delta (k-| {\boldsymbol{\kappa }}| ){\rm{d}}{\boldsymbol{\kappa }}$, where k is the wavenumber, $\widehat{q}({\boldsymbol{\kappa }})$ denotes the Fourier transform of $q({\boldsymbol{x}});$ $\delta (k)$ is the Dirac delta function.

Figure 16 shows power spectra of the density ($P(\rho ,k)$, left) and projected density ($P({\rm{\Sigma }},k)$, right) for cases A–C compensated with ${k}^{0.2}$ and ${k}^{1.2}$, respectively so that the scaling range is approximately flat. Unlike in isothermal turbulence at sonic Mach numbers 5–10, where $P(\rho ,k)\propto {k}^{-1}$, the spectra in our multiphase simulations are substantially more shallow. This is in part due to the modest resolution of the simulations, which is responsible for very fragmented under-resolved cold phase at high densities. The spectra for moderately and weakly magnetized cases B and C are nearly identical, and the case-A spectrum has only a mild power offset and the same slope. Remarkably, the projected density spectra, which involve 2D Fourier transforms and integration over annuli in k-space, are very similar in shape and have a slope offset of −1, i.e. $P({\rm{\Sigma }},k)\approx {L}^{2}{k}^{-1}P(\rho ,k)/2$. Thus the Σ-spectrum is a very reliable proxy to the ρ-spectrum, assuming isotropy. This important property of the column density spectrum provides direct access to the scale-dependent gravitational potential energy distribution in self-gravitating turbulence (section 3.13 and [94]). We also note that the measured slope of the column density spectrum $P({\rm{\Sigma }},k)\propto {k}^{-1.2}$ implies a very weak dependence of Σ on scale, ${\rm{\Sigma }}\propto {r}^{0.1}$ for $r\in [5,60]$ pc.

Zoom In Zoom Out Reset image size

Figure 17 (left) shows compensated velocity power spectra $P({\boldsymbol{u}},k)$ for the same three cases. The spectra show similar scaling properties $P\propto {k}^{3/2}$ at $k/{k}_{\min }\in [10,30]$, which most likely does not reflect the inertial range scaling in the limit of ${Re}\,\gg \,1$ due to the limited resolution. The bump in the spectrum at $k\lt 3{k}_{\min }$ reflects the random forcing operating at these scales. The 'true' inertial range is most likely contaminated by overlapping effects of the forcing and numerical dissipation.

Zoom In Zoom Out Reset image size

Figure 17 (right) shows results of the Helmholtz decomposition of the velocity field into dilatational (curl-free) and solenoidal (divergence-free) components, ${\boldsymbol{u}}={{\boldsymbol{u}}}_{{\rm{d}}}+{{\boldsymbol{u}}}_{{\rm{s}}}$. The dilatational-to-solenoidal ratio $P({{\boldsymbol{u}}}_{{\rm{d}}},k)/P({{\boldsymbol{u}}}_{{\rm{s}}},k)$ starts below ∼0.3 in case A and increases to ∼0.4 in case C with decreasing magnetization, consistent with expectations for solenoidally driven supersonic isothermal MHD turbulence [95] and for supernovae-driven multiphase ISM turbulence [96]. Since our case includes the baroclinic effect otherwise missing in isothermal simulations, we quantified its contribution to the production of enstrophy in stationary multiphase turbulence, using the formalism based on the enstrophy equation [97]

$$\begin{eqnarray}&&{\partial }_{t}{\rm{\Omega }}={S}_{{\rm{s}}}+{S}_{{\rm{b}}}+{S}_{{\rm{c}}}+{S}_{{\rm{m}}}+{S}_{{\rm{d}}}.\end{eqnarray} \tag{ 3 }$$

Here ${\rm{\Omega }}=\langle {\omega }^{2}/2\rangle$, is the enstrophy $\langle \cdot \rangle$ denotes averaging over the triply periodic domain, ${\boldsymbol{\omega }}={\boldsymbol{\nabla }}\times {\boldsymbol{u}}$ is the vorticity, and the dissipative term ${S}_{{\rm{d}}}$ cannot be evaluated explicitly in ILES simulations. Following [97], we also dropped the contribution from random forcing, which is small since the large-scale acceleration ${\boldsymbol{a}}$ is mostly uncorrelated with small-scale vorticity ${\boldsymbol{\omega }}$. The remaining four terms of interest here include the vortex stretching term ${S}_{{\rm{s}}}=\langle {\boldsymbol{\omega }}\,{\boldsymbol{\cdot }}\,({\boldsymbol{\omega }}\,{\boldsymbol{\cdot }}\,{\boldsymbol{\nabla }}){\boldsymbol{u}}\rangle$, the baroclinic term ${S}_{{\rm{b}}}=\langle {\boldsymbol{\omega }}\,{\boldsymbol{\cdot }}\,({\boldsymbol{\nabla }}\rho \times {\boldsymbol{\nabla }}p)/{\rho }^{2}\rangle$, the compression term ${S}_{{\rm{c}}}=-\langle {\omega }^{2}{\boldsymbol{\nabla }}\,{\boldsymbol{\cdot }}\,{\boldsymbol{u}}/2\rangle$, and the magnetic term ${S}_{{\rm{m}}}=\langle {\boldsymbol{\omega }}\,{\boldsymbol{\cdot }}\,({\boldsymbol{\nabla }}\rho \times {\boldsymbol{\nabla }}{p}_{\mathrm{mag}})/{\rho }^{2}+{\boldsymbol{\omega }}\,{\boldsymbol{\cdot }}\,{\boldsymbol{\nabla }}\times ({\boldsymbol{B}}\,{\boldsymbol{\cdot }}\,{\boldsymbol{\nabla }}{\boldsymbol{B}}/\rho )\rangle$.

Table 4 summarizes the relative contributions of these four sources in statistically stationary conditions for cases A–C. The magnetic term ${S}_{{\rm{m}}}$ strongly dominates enstrophy production in all cases, but contributes slightly less under weaker magnetization. The baroclinic effect ${S}_{{\rm{b}}}$ is generally weak and contributes even less in presence of a stronger mean field. Vortex stretching ${S}_{{\rm{s}}}$ is always the second largest contributor, but magnetic tension progressively reduces its effect in cases with stronger magnetization. Vorticity amplification by compression ${S}_{{\rm{c}}}$ (e.g. in shocks, where a net alignment of the enstrophy gradient and velocity vectors is common [97]) is not sensitive to the magnetization level. The overall levels of enstrophy given in the first column of table 4 correlate positively with the magnetization level, being 1.5 times larger in case A, compared to C. Meanwhile, the characteristic time of enstrophy production by all four sources, $\tau ={\rm{\Omega }}/({S}_{{\rm{s}}}+{S}_{{\rm{b}}}+{S}_{{\rm{c}}}+{S}_{{\rm{m}}})$, is roughly the same in all three cases. Since this time scale $\tau \approx 0.25\,\mathrm{Myr}$ is short compared to the system lifetime, it naturally has to do with the small-scale processes in the turbulent ISM.

Table 4.  Enstrophy, source terms, and associated time scales.

Case	Ω	${S}_{{\rm{s}}}$	${S}_{{\rm{b}}}$	${S}_{{\rm{c}}}$	${S}_{{\rm{m}}}$	τ
	(Myr−2)	(%)	(%)	(%)	(%)	(Myr)
A	31.7	12.1	2.4	11.0	74.5	0.23
B	28.1	14.7	2.7	11.1	71.5	0.26
C	21.8	18.3	5.4	11.9	64.4	0.28

The total energy conserved by the ideal system without dissipation, heating, cooling, and forcing includes the kinetic, internal, and magnetic energy contributions: $E=K+U+M$, respectively. Realistic ISM turbulence at length scales ∼100 pc, is dominated by the kinetic and magnetic energy constituents, hence $U\ll K\sim M$. The individual spectral densities are defined by $K(k)=P({\boldsymbol{j}},{\boldsymbol{u}};k)/2$ and $M(k)=P({\boldsymbol{b}},k)/2+\delta (k){B}_{0}^{2}/2$. Here the kinetic energy spectral density is represented by a cospectrum⁵ of the momentum density ${\boldsymbol{j}}\equiv \rho {\boldsymbol{u}}$ and velocity ${\boldsymbol{u}}$. Both spectral densities satisfy Parseval's theorem: $K={\int }_{0}^{\infty }K(k){\rm{d}}{k}$ and $M={\int }_{0}^{\infty }M(k){\rm{d}}{k}$ [94]. An alternative approach to calculating K(k) for compressible turbulence, using the power spectrum of $\sqrt{\rho }{\boldsymbol{u}}$ [98], returns spectra with a large artificial bottleneck at high wavenumbers, originating due to strong autocorrelation of $\sqrt{\rho }$, which has a broad distribution with a width of ∼3 dex in these models (figure 6-left). Our choice of the cospectrum $P({\boldsymbol{j}},{\boldsymbol{u}};k)$ is motivated by recently derived exact relations for compressible isothermal turbulence [99] and numerical experiments [100].

Figure 18 shows the kinetic and magnetic energy spectra as well as their sum, which scale $\propto {k}^{-1.5}$ at this resolution, for cases A–C. If the inertial range is sufficiently resolved, detailed kinetic-to-magnetic energy equipartition $K(k)\approx M(k)$ should hold and both K(k) and M(k) are expected to have scaling exponents in a range from −2 to $-3/2$. All three cases indeed show that energy equipartition $K(k)\sim M(k)$ holds within a factor of two or less with magnetic and kinetic energy spectra running roughly parallel to each other over $\gt 1$ dex in k in case C and $\gt 1.5$ dex in case A.

Zoom In Zoom Out Reset image size

While at length scales ∼100 pc the internal energy is clearly subdominant, it may still play a role on smaller scales, e.g., where turbulence in the bulk of the volume becomes subsonic. To probe the scale-dependent part of U, we define the spectral density of internal energy, following a template for isothermal supersonic turbulence developed in [94], by $U(k)\equiv P(\rho ,e;k)/2+\delta (k)U/2$ ⁶ . Internal energy spectra computed for cases A–C are very shallow, $U(k)\propto {k}^{-1}$ (figure 19), most likely because thermal processes are not properly resolved and the cold phase is highly fragmented. Unlike in the isothermal case, where ρ and $e={c}_{{\rm{s}}}^{2}\mathrm{ln}(\rho /{\rho }_{0})$ are always positively correlated, the cospectra here are negative at large scales due to anti-correlation of density and temperature perturbations in the quasi-isobaric TI regime. The spectra then switch to positive on small scales (not shown) under quasi-adiabatic conditions. The shallow $U(k)\propto {k}^{-1}$ scaling in combination with the quasi-isobaric conditions in the multiphase gas create a depression in the total energy spectrum E(k) centered at ${k}_{{\rm{p}}}\sim 50{k}_{\min }$. The negative thermal pressure 'support' at the associated 'crossover' length scale ${\lambda }_{{\rm{p}}}=2\pi /{k}_{{\rm{p}}}$ would make objects of this size more prone to collapse, if their mean density is sufficiently high. Our models, however, do not resolve this scale properly. Thus we can only resort to extrapolation to get a rough estimate for ${\lambda }_{{\rm{p}}}$, which returns a fraction of a parsec (figure 19). Higher resolution simulations would show whether or not this crossover scale defines the characteristic width of star forming filaments measured by the Herschel satellite [101].

Zoom In Zoom Out Reset image size

Finally, in order to make predictions for self-gravitating turbulence, we can readily estimate the gravitational potential energy spectral density, corresponding to the density spectrum in figure 16, as $W(k)=-2\pi {{Gk}}^{-2}P(\rho ,k)\propto {k}^{-2.2}$ [94]. If we take the stationary state of developed turbulence as initial conditions for a simulation with self-gravity, the gravitational potential energy would initially have a slightly steeper scaling compared to the sum of kinetic and magnetic energy. Also, the (negative) spectral density of the potential energy is smaller in absolute value than E(k) in our cases with L = 200 pc, ensuring global gravitational stability of the ISM and hence purely local character of star formation.