The structure and statistics of interstellar turbulence

TL;DR

This paper uses advanced MHD simulations to analyze the structure, statistics, and observable properties of interstellar turbulence in the Milky Way's ISM, linking turbulence characteristics to star formation rates.

Contribution

It provides a comprehensive simulation-based analysis of multiphase, magnetized ISM turbulence, matching many observed features and revealing how turbulence influences star formation.

Findings

Reproduces observed molecular cloud properties

Predicts non-Gaussian magnetic field PDFs

Shows turbulence controls star formation rates

Abstract

We explore the structure and statistics of multiphase, magnetized ISM turbulence in the local Milky Way by means of driven periodic box numerical MHD simulations. Using the higher order-accurate piecewise-parabolic method on a local stencil (PPML), we carry out a small parameter survey varying the mean magnetic field strength and density while fixing the rms velocity to observed values. We quantify numerous characteristics of the transient and steady-state turbulence, including its thermodynamics and phase structure, kinetic and magnetic energy power spectra, structure functions, and distribution functions of density, column density, pressure, and magnetic field strength. The simulations reproduce many observables of the local ISM, including molecular clouds, such as the ratio of turbulent to mean magnetic field at 100 pc scale, the mass and volume fractions of thermally stable HI, the lognormal distribution of column densities, the mass-weighted distribution of thermal pressure, and the linewidth-size relationship for molecular clouds. Our models predict the shape of magnetic field probability density functions (PDFs), which are strongly non-Gaussian, and the relative alignment of magnetic field and density structures. Finally, our models show how the observed low rates of star formation per free-fall time are controlled by the multiphase thermodynamics and large-scale turbulence.

Tables

Table 1. Table 1: Input parameters.

Case	$N$	$n_0$	$u_{\mathrm{rms},0}$	$B_0$	${\beta }_{\mathrm{th},0}$	${\beta }_{\mathrm{turb},0}$	${ℳ}_{a,0}$
		cm-3	km/s	$\mu$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

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

Case	$B_{\mathrm{rms}}$	$b_{\mathrm{rms}}$	$b_{\mathrm{rms}}/B_0$			$n_{\mathrm{rms}}/n_0$	$u_{\mathrm{rms}}$	$v_{a,\mathrm{rms}}$	${ℳ}_{s,v}$	${ℳ}_{s,m}$	${ℳ}_a$
	$\mu$G	$\mu$G					km/s	km/s
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 3. Table 3: Mass and volume fractions, Mach number regimes for different thermal phases.

Case		$F_v$			$F_m$			${ℳ}_s$			${ℳ}_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

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

Case	$\Omega$	$S_s$	$S_b$	$S_c$	$S_m$	$\tau$
	(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

Equations

$${\partial }_t\rho$$

$${\partial }_t(\rho \mathbf{u})$$

$${\partial }_t\mathbf{B}$$

$${\partial }_t\mathcal{E}$$

$$S_m(r)=⟨[u(\mathbf{x}+\mathbf{r})-u(\mathbf{x}){]}^m⟩\propto r^{{\zeta }_m},$$

$${\partial }_t\Omega =S_{\mathrm{s}}+S_{\mathrm{b}}+S_{\mathrm{c}}+S_{\mathrm{m}}+S_{\mathrm{d}}.$$

Full text

1 Introduction

Turbulence is ubiquitous in the interstellar medium (ISM) [1]. It organizes the ISM structure to optimize and direct energy fluxes across an enormous range of length scales within the Milky Way disk. Interstellar turbulence, however, differs from the familiar Kolmogorov homogeneous and isotropic incompressible case [2, 3] in several important respects: it is neither isotropic nor homogeneous [4]; the interstellar gas is highly compressible and magnetized; and its thermal energy is not strictly conserved. Instead, the ISM-specific radiative cooling and volumetric heating functions determine an equilibrium multiphase nature of the neutral ISM [5, 6, 7]. Nonlinear advection couples with nonlinear radiative cooling, further complicating the treatment of the ISM. We hence deal with a classical example of a hierarchical nonequilibrium dissipative system, receiving energy primarily at large scales ($\sim 100$ pc) from stellar feedback, extragalactic gas accretion, large-scale shear, and gravitational instabilities within the disk and loosing energy through small-scale dissipation ($10^{-4}$ pc) and local radiative cooling [8].

The concept of a hierarchy of ‘eddies’ loosely defined as coherent patterns in the velocity field of incompressible turbulence, is replaced in the case of the ISM by an even more indistinct hierarchy of cloud complexes, clouds, clumps and cores [9]. Observational cloud definitions are usually tied to specific tracers and are essentially two-dimensional –due to complex line-of-sight convolutions– making straightforward physical interpretation of observations extremely difficult [10]. Since interstellar turbulence controls the structure, dynamics and chemistry of the ISM, it is a key building block of any future successful theory of star formation in molecular clouds [8, 11, 12, 13].

To explore this complexity, we develop self-consistent models based on very idealized numerical experiments, linking together scales relevant to molecular cloud formation and fragmentation. We exploit the concept of self-organization in nonequilibrium nonlinear dissipative systems [14] in application to the ISM, which is long overdue [15]. We treat interstellar turbulence as an agent that imposes ‘order’ in the form of coherent structures and correlations between flow fields emerging in a simple periodic box simulation when a statistically stationary state fully develops. In this case, the details of initial conditions are no longer important. Instead the steady state provides realistic initial conditions for star formation. Unlike various flavors of the popular converging flow scenario [16, 17, 18, 19], this approach allows us to reproduce basic ISM observables, including properties of local molecular clouds, with only a few control parameters. The model can be readily extended to include more physics (e.g., radiative transfer and chemistry), to augment the range of resolved scales, and to close the star formation feedback loop.

The paper is organized as follows. Section 2 contains the details of our numerical models, including equations, initial and boundary conditions, and parameters. Section 3 presents the results of statistical analysis and comparison with observations for a number of ISM diagnostics. § 3.1 describes the transition to turbulence and global properties of the statistically stationary state. § 3.2 summarizes the evolution of magnetic field under the action of random large-scale external forcing. In § 3.3 we discuss thermodynamics and phase structure of the ISM emerging in the simulations. § 3.4 accesses the turbulent Mach number regimes of various thermal phases. § 3.5 presents predictions for star formation rate controlled by the turbulence and thermodynamic properties of the ISM. §§ 3.6 and 3.7 discuss density, column density and thermal pressure distributions. §§ 3.8, 3.9, and 3.10 present statistics of the magnetic field. §§ 3.11, 3.12, and 3.13 discuss scaling properties of multiphase turbulence and energy spectral densities. We conclude with a summary of results in Section 4.

2 Numerical models

2.1 Governing equations, initial and boundary conditions

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

[TABLE]

Here $\rho$ is the gas density, $\bm{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 ${\cal E}=\rho u^{2}/2+\rho e+B^{2}/2$ and $\langle{\cal E}\rangle=K+U+M$ is one of the invariants of an ideal system (with zero r.h.s.). The generalized heat-loss function $\rho{\cal L}=n^{2}\Lambda(T)-n\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, $\bm{B}$, is subject to the usual solenoidal constraint $\bm{\nabla}\bm{\cdot}\bm{B}=0$. The following initial conditions: $\rho_{0}+\delta\rho$ (where $\delta\rho$ represent 20% random density perturbations), $p_{\rm 0}$, ${\bm{u}}_{0}=0$, and ${\bm{B}}_{0}=(B_{0},0,0)$ are applied. Random forcing implementation assumes zero momentum input $\bm{f}\equiv\rho\bm{a}-\langle\rho{\bm{a}}\rangle$ and applies a fixed-in-time large-scale solenoidal non-helical acceleration $\bm{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$ Myr and, when a finite initial velocity boost $\delta\bm{u}=\tau\bm{a}$ is applied. Here $\tau$ 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_{\rm rms}L/\nu_{\rm eff}$, Prandtl number $Pr=c_{p}\rho_{0}\nu_{\rm eff}/\kappa_{\rm eff}$ and magnetic Prandtl number $Pm=\nu_{\rm eff}/\eta_{\rm eff}$ can be defined ($\nu_{\rm eff}$, $\kappa_{\rm eff}$, and $\eta_{\rm 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$. Heat conduction can stabilize thermal instability (TI) on scales below the so-called Field length $\lambda_{\rm F}=\sqrt{\kappa T/\rho^{2}\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$ pc 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 ${\cal 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 ${\cal 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.

2.2 Input parameters

Each model is defined by five control parameters: the domain size $L$; the mean gas number density in the box $n_{0}$; the rms velocity $u_{\rm 0,rms}$ associated with the driving force amplitude; the uniform magnetic field strength $B_{0}$; and the grid resolution $N^{3}$.

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_{0}$ in the range from $2-5$ cm*-3* reflect local average conditions in a spiral arm-like setting with Hi number densities somewhat higher than the typical mean value of 1 cm*-3*. Likewise, the values of $u_{\rm 0,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].

The table also provides the magnetic field parametrization in terms of the plasma beta $\beta_{\rm th,0}=p_{0}/(B_{0}^{2}/2)$, representing the ratio of the initial thermal pressure of the gas $p_{0}$ to the magnetic pressure due to the applied uniform field $B_{0}$. The values of turbulent plasma beta (defined with the dynamic pressure $\rho u^{2}$ replacing thermal $p$ in the classical definition) $\beta_{\rm turb,0}=\rho_{0}u^{2}_{\rm rms,0}/(B_{0}^{2}/2)$ and the initial Alfvénic Mach number ${\cal M}_{a,0}=u_{\rm rms,0}/(B_{0}/\sqrt{\rho_{0}})$ are also given for convenience. Overall, cases A, B, and C trace a transition from trans-Alfvénic (${\cal M}_{\rm a,0}\sim 1$) to super-Alfvénic (${\cal M}_{\rm a,0}=4,13$) magnetization levels, while cases D and E are only moderately super-Alfvénic (${\cal M}_{\rm a,0}\sim 2-3$). Finally, the total gas mass in the domain for cases A, B, 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 $\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 (PAHs) 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<26$ K, however there is no temperature floor and the gas temperature drops to $\sim 4$ K in strong rarefactions adiabatically. The density-independent heating rate is set at $\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 $\Gamma=(\rho/m_{\rm H})\Lambda(T)$ in the pressure-density plane has two extrema separating the two stable branches: a maximum at $T_{\rm max}=5250$ K and a minimum at $T_{\rm min}=184$ K (Fig. 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 ${\cal L}(\rho,T)\neq 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>T_{\rm max}$ and the cold phase assumes $T<T_{\rm min}$. The range of intermediate temperatures $T_{\rm min}<T<T_{\rm max}$ roughly captures the thermally unstable regimes. With these definitions, our warm phase includes both warm neutral medium (WNM) and warm ionized medium (WIM).

3 Results

3.1 Energy evolution

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_{\rm 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\bm{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\bm{B}^{2}/2\rangle$ (blue), starts at a level determined by the mean field $B_{0}$, 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{\cal E}\rangle=K+U+M$, is shown in red. It is initially ($t<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<M\lesssim K$.

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.

3.2 Magnetic field amplification and global turbulence characteristics

The details of magnetic field growth are illustrated in Fig. 3, where thick and thin lines show the total field $B_{\rm rms}$ and fluctuations $b_{\rm 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_{0}$, together with the mean density, $n_{0}$, and the rms velocity, $u_{\rm rms,0}$, control the level of magnetic fluctuations in the developed turbulence. When $B_{0}$ 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_{0}$, with the regular field generated by the Galactic dynamo and spatially ordered on scales $>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\leq b_{\rm rms}/B_{0}\leq 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, 36, 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_{\rm rms}/B_{0}$ is achieved, local magnetic field properties and cloud magnetization levels can be considered as realistic.

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_{\rm th}$ is mostly below unity, except in case C, indicating subdominant contribution of the internal energy compared to the magnetic one. The mean $\beta_{\rm 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_{\rm 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 ${\cal M}_{\rm s,v}$ and ${\cal M}_{\rm s,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 ${\cal 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<{\cal M}_{\rm s}<13$.

3.3 Thermodynamics of the turbulent two-phase ISM

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 phase diagrams for data snapshots taken at $t=42.3$ Myr in cases B, D, and E. As in Fig. 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.

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).

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>100$ cm*-3*) 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 [40, 41, 30, 42]) 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 [43], e.g. $\gtrsim 48$% by mass of the warm neutral medium (WNM, $>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_{0}$ and on the rms velocity $u_{\rm 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_{\rm 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_{0}$ as B, but $2\times$ smaller $u_{\rm 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_{0}$ and $u_{\rm rms}$, but different $B_{0}$. However, as far as the $p_{\rm th}-n$ phase diagrams are concerned, the differences are quite subtle. Hence we defer the discussion to § 3.5 with a primary focus on control of star formation.

3.4 Mach number regimes of thermal phases

Table 3 also lists average Mach numbers for the cold, warm, and unstable phases. The sonic Mach number ${\cal 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 (§ 3.11). The ${\cal M}_{\rm s}-n$ correlation, shown in Fig. 5 (left), follows the prediction ${\cal M}_{\rm s}\propto n^{1/2}$ for isobaric conditions at intermediate densities $n\in[1,10]$ cm*-3* [44] reasonably well and flattens into ${\cal M}_{\rm s}\approx const$ in quasi-isothermal conditions at low and high density extremes (§ 3.6).

The Alfvén Mach number ${\cal 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. Fig. 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>10^{2}$ cm*-3* having local ${\cal M}_{\rm a}>1$. A similar trend is seen in isothermal MHD simulations [45] in supersonic and weakly super-Alfvénic regimes, which can be explained by the dependence ${\cal M}_{\rm a}(B,\rho)$ and the relationship between $B$ and $\rho$. Our multiphase cases combine a transonic, sub-Alfvénic behavior at low densities with a supersonic and super-Alfvénic one at high densities.

3.5 Control of star formation

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 Fig. 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\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_{\rm ff}=\sqrt{3\pi/32G\rho}\approx\sqrt{(10^{3}\;{\rm cm}^{-3})/n}$ Myr, where $G$ is the gravitational constant [46].

This Jeans-unstable gas occupies a tiny fraction (0.01%) of the computational volume and its mass $\Delta{\mathfrak{M}}$ comprises 1.7% of the mass in case B at $t=42.3$ 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 [47, 13] $\epsilon_{\rm ff}\equiv t_{\rm ff}/t_{\rm dep}{=\Delta{\mathfrak{M}}/{\mathfrak{M}}}\lesssim 0.017$, where the gas depletion time $t_{\rm 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_{\rm ff}\sim 0.01-0.1$ for clouds with mean densities $n_{\rm H_{2}}\sim 10^{3}$ cm*-3* [48].

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_{\rm 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_{\rm 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_{\rm ff}$ include the resolution dependence, which enters the Truelove constraint [38].

We have traced the dependence of $\epsilon_{\rm 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, B, and C, which have the same $n_{0}$ and $u_{\rm rms}$, but different $B_{0}$. The results: $\langle\epsilon_{\rm ff}\rangle=0.021\pm 0.005$ (A), $0.013\pm 0.002$ (B), and $0.011\pm 0.002$ (C) indicate an overall $\sim 2\times$ drop in $\langle\epsilon_{\rm 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_{\rm th}}$ [39]. If we apply this $\beta_{\rm th}$-dependent correction, the average rates drop by a factor of $\sim 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_{\rm 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_{\rm 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-20)$%. A better understanding of the functional form of $\epsilon_{\rm ff}(L,n_{0},u_{\rm 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_{\rm ff}$, to undergo gravitational phase transition. Star formation rates per free-fall time $\epsilon_{\rm ff}\sim(0.001-0.1)$ we measured in cases B, E, and D are controlled by the problem parameters $(L,n_{0},u_{\rm rms,0},B_{0},\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_{\rm 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.

3.6 Density and column density PDFs

The gas density PDF is one of the basic statistics of compressible turbulence that bears unique signatures of gas thermodynamics [54], gravitational instability [46], and the presence of shocks. For instance, in isothermal homogeneous supersonic turbulence, the PDF is lognormal [56, 57, 58]. 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_{\rm eq}\propto\rho_{\rm eq}^{\gamma_{\rm eff}}$ with $\gamma_{\rm eff}\approx 0.7<1$ [30, 59, 60]. 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_{\rm rms}\sim 5-10$ km s*-1* into an asymmetric single-peak PDF. Density PDFs for cases A, B, and C shown in Fig. 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 (Fig. 4) and (ii) flattening of the Mach number–density correlation ${\cal M}_{\rm s}\propto\rho^{1/2}$, consistent with quasi-isobaric conditions [44], above $n\sim 30$ cm*-3* (Fig. 5). It is ultimately the lack of ${\cal 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 (§ 3.7) and local linewidth–size relationship (§ 3.11).

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. in the Draco cloud [61]). Note that the sonic Mach numbers for cases A, B, and C are generally higher than 5 (Table 2). At ${\cal M}_{\rm s,v}\lesssim 2$, the lognormal property of the low tail of N-PDFs is lost [60].

Since temperature plays the role of an important third dimension in column-density mapping of star forming clouds [62], we proceed by computing N-PDFs conditioned on the temperature $\wp(\log N|\Delta T)$, such that $\int\wp(\log N|\Delta T)d\log N=1$ and $\Delta T$ defines the allowed temperature interval. As an example, Fig. 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<T[{\rm K}]<184$) and warm (red, $T>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\ln(1+b^{2}{\cal M}_{\rm s}^{2}$) based on simulations of isothermal turbulence [63], if we assume ${\cal M}_{\rm s}$ values from Table 3 and $b\approx 1/3$ due to a divergence-free nature of the forcing [64]. Similar discrepancy of the widths of fitted lognormal distributions with isothermal predictions was seen earlier in multiphase ISM simulations without a magnetic field [60], in GALFA-Hi data for the Perseus molecular cloud [66], and in data for seven molecular clouds from the Leiden/Argentine/Bonn Galactic Hi survey [67]. In all cases, the widths of observational Hi N-PDFs are similar to our model predictions.

Self-shielding of H2 becomes effective around $N_{\rm H}\sim 2\times 10^{20}$ cm*-2*, while column densities in excess of $N_{\rm H}\sim 4\times 10^{21}$ cm*-2* are required to form CO [68, 69]. The range of column densities in between is populated by the ‘CO-dark H2’ gas that may account on average for $(30-40)$% of the Milky Way molecular mass [70, 71]. The shaded areas in both panels of Fig. 8 show the column density interval, where CO-dark H2 is expected to reside. Fig. 8 (right) details the temperature dependence of N-PDFs for the cold gas in case B. In all cases, except perhaps for the $T<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 Fig. 6 and indicates that corresponding observational diagnostics are not very sensitive to the ISM magnetization levels.

3.7 Thermal pressure PDFs

At high levels of turbulence, the distribution of thermal pressure for cases A, B, and C spans about 6 decades leaving no room for the classical isobaric pressure-supported cloud picture in the violent ISM (Fig. 9, left). The PDFs are asymmetric in a very characteristic way, showing a positive skewness, reported in simulations earlier [55, 73, 74, 75]. 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 (Fig. 4). The peak tends to shift slightly to higher pressures as the magnetization becomes stronger, approaching the mean ISM pressure of 3700 K cm*-3* reported in [72] in the range of $B_{0}$ 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, B, and 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_{\rm 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.

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 $\sim 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 [72].111The HST data points provided by Dr. Ed Jenkins in 2010 reflect an earlier version of the analysis that relied on Wolfire et el. (1995) description of thermal processes [28], which is consistent with our simulations (§ 2.2) and differs from [72], where the Wolfire et al. (2003) [7] version was used. We thus exclude lines of sight with $N({\rm H}{\rm\sc I})>2.5\times 10^{21}$ cm*-2* and mask out the gas with temperatures outside $40$ K$<T<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 Fig. 7.

Distributions for cases A through C with $u_{\rm 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_{0}$. The width of the distribution, however, is sensitive to $u_{\rm rms}$ and $n_{0}$. 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 modelled 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.

3.8 Magnetic, dynamic, and thermal pressure

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_{\rm mag}\equiv B^{2}/2$ versus the thermal pressure $p_{\rm th}=(\gamma-1)\rho e$ (left) and versus the dynamic pressure $p_{\rm dyn}\equiv\rho u^{2}$ (right) for case B at 42.3 Myr. Since $\beta_{\rm th}\equiv p_{\rm th}/p_{\rm mag}=(\gamma-1)U/M$, $\beta_{\rm turb}\equiv p_{\rm dyn}/p_{\rm mag}=2K/M=2{\cal M}_{\rm a}^{2}$, and $\beta_{\rm turb}/\beta_{\rm th}=\gamma{\cal 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_{\rm 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<100$ K and $n>100$ cm*-3*), as indicated by the color contour map in Fig. 10.

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_{\rm turb}=1$ in the right panel shows. This also means that $K\sim M$ and ${\cal 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_{\rm turb}=30$ indicate. The cold molecular gas in our simulations is instead super-Alfvénic with ${\cal M}_{\rm a}\sim 8$ (see also Fig. 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_{\rm eff}$ that do not properly approximate conditions in the multiphase ISM.

3.9 Filaments, striations, and the alignment of magnetic and velocity fields

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, 82, 83, 84, 85, 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(\bm{B}\cdot\bm{u}/\sqrt{B^{2}u^{2}})$ for cases A, B, and C (left) and details the alignment properties as a function of density for case A (right). The $\bm{B}-\bm{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_{0}$ (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 (Fig. 11, right). Indeed, in the magnetically dominated (sub-Alfvénic) warm phase, the gas structure is expected to be preferentially aligned with the local magnetic field.

For the dense cold phase with super-Alfveńic 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 $\bm{B}$ and the eigenvectors $\bm{\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 amplification of the potential (curl-free) component of the gas velocity.

3.10 Magnetic field PDFs

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 $\bm{b}=\bm{B}-\bm{B}_{0}$ perpendicular ($b_{\perp}$) and parallel ($b_{\parallel}$) to the mean field $\bm{B}_{0}$ for cases A, B, and C. All PDFs represent averages over $70-100$ flow snapshots, uniformly covering the statistically stationary state. In the case of $b_{\perp}$, the results include both $b_{y}$ and $b_{z}$ components normal to $\bm{B}_{0}$ directed along the positive $x$-coordinate direction.

The plots show that the PDFs are strongly non-Gaussian in all cases and that the $B$-field fluctuations in the direction of $\bm{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_{0}$. 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_{\parallel}$ 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 $|\bm{B}|$ for cases A, B, and C; vertical black dashed lines indicate the values of $B_{0}$. 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]$ $\mu$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_{0}$ 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.

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 $|\bm{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>30$ cm*-3*, 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<2$ cm*-3*). 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.

3.11 Velocity structure functions

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

[TABLE]

where $\langle\cdot\rangle$ denotes the ensemble average over all possible point pairs separated by the lag $r=|\bm{r}|$. The velocity components here can either be parallel or perpendicular to vector $\bm{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, B, and 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 ${\cal M}_{\rm s}\approx 13$ and ${\cal 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.

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 (1b) 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 excessively 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<184$ K, purple) and for the warm phase ($T>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 ${\cal M}_{\rm s}\approx 15$ for the cold phase (Table 3) and $\zeta_{1}\sim 1/2$ are common for isothermal turbulence at ${\cal M}_{\rm s}\sim 10$ [58]. 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.

3.12 Velocity and density power spectra

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(\bm{x})$, the spectrum $P(q,k)\equiv\int|\widehat{q}(\bm{\kappa})|^{2}\delta(k-|\bm{\kappa}|)d\bm{\kappa}$, where $k$ is the wavenumber, $\widehat{q}(\bm{\kappa})$ denotes the Fourier transform of $q(\bm{x})$; $\delta(k)$ is the Dirac delta function.

Figure 16 shows power spectra of the density ($P(\rho,k)$, left) and projected density ($P(\Sigma,k)$, right) for cases A, B, and 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(\Sigma,k)\approx L^{2}k^{-1}P(\rho,k)/2$. Thus the $\Sigma$-spectrum is a very reliable proxy to the $\rho$-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 (§ 3.13 and [97]). We also note that the measured slope of the column density spectrum $P(\Sigma,k)\propto k^{-1.2}$ implies a very weak dependence of $\Sigma$ on scale, $\Sigma\propto r^{0.1}$ for $r\in[5,60]$ pc.

Figure 17 (left) shows compensated velocity power spectra $P(\bm{u},k)$ for the same three cases. The spectra show similar scaling properties $P\propto k^{-3/2}$ at $k/k_{\rm 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<3k_{\rm 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.

Figure 17 (right) shows results of the Helmholtz decomposition of the velocity field into dilatational (curl-free) and solenoidal (divergence-free) components, $\bm{u}=\bm{u}_{\rm d}+\bm{u}_{\rm s}$. The dilatational-to-solenoidal ratio $P(\bm{u}_{\rm d},k)/P(\bm{u}_{\rm s},k)$ starts below $\sim 0.3$ in case A and increases to $\sim 0.4$ in case C with decreasing magnetization, consistent with expectations for solenoidally driven supersonic isothermal MHD turbulence [94] and for supernovae-driven multiphase ISM turbulence [95]. 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 [96]

[TABLE]

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

Table 4 summarizes the relative contributions of these four sources in statistically stationary conditions for cases A, B, and 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 [96]) 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=\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$ Myr is short compared to the system lifetime, it naturally has to do with the small-scale processes in the turbulent ISM.

3.13 Energy spectral densities

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 $\sim 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(\bm{j},\bm{u};k)/2$ and $M(k)=P(\bm{b},k)/2+\delta(k)B_{0}^{2}/2$. Here the kinetic energy spectral density is represented by a cospectrum222The cospectrum $P(\bm{j},\bm{u};k)$ is the Fourier transform of the symmetric part of the cross-covariance function, ${\cal C}_{ju}(\bm{r})=\langle\bm{j}\bm{\cdot}\bm{u}^{\prime}+\bm{j}^{\prime}\bm{\cdot}\bm{u}\rangle/2$, integrated over spherical shells, i.e. $P(\bm{j},\bm{u};k)\equiv\int\widehat{{\cal C}_{ju}}(\bm{\kappa})\delta(k-|\bm{\kappa}|)d\bm{\kappa}$. of the momentum density $\bm{j}\equiv\rho\bm{u}$ and velocity $\bm{u}$. Both spectral densities satisfy Parseval’s theorem: $K=\int_{0}^{\infty}K(k)dk$ and $M=\int_{0}^{\infty}M(k)dk$ [97]. An alternative approach to calculating $K(k)$ for compressible turbulence, using the power spectrum of $\sqrt{\rho}\bm{u}$ [99], 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 $\sim 3$ dex in these models (Fig. 6-left). Our choice of the cospectrum $P(\bm{j},\bm{u};k)$ is motivated by recently derived exact relations for compressible isothermal turbulence [100] and numerical experiments [101].

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, B, and 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 $>1$ dex in $k$ in case C and $>1.5$ dex in case A.

While at length scales $\sim 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 [97], by $U(k)\equiv P(\rho,e;k)/2+\delta(k)U/2$.333Here, as in the case of kinetic energy spectral density, the cospectrum $P(\rho,e;k)$ is the Fourier transform of the symmetric part of the cross-covariance function ${\cal C}_{\rho e}(\bm{r})=\langle\rho e^{\prime}+\rho^{\prime}e\rangle/2$ integrated over spherical shells in $k$-space. Also note the $k=0$ component, which is similar the mean-field energy in the expression for $M(k)$, and that Parseval’s theorem $U=\int_{0}^{\infty}U(k)dk$ holds in this case as well. Internal energy spectra computed for cases A, B, and C are very shallow, $U(k)\propto k^{-1}$ (Fig. 19), most likely because thermal processes are not properly resolved and the cold phase is highly fragmented. Unlike in the isothermal case, where $\rho$ and $e=c_{\rm s}^{2}\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_{p}\sim 50k_{\rm min}$. The negative thermal pressure ‘support’ at the associated ‘crossover’ length scale $\lambda_{p}=2\pi/k_{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_{p}$, which returns a fraction of a parsec (Fig. 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 [98].

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 Fig. 16, as $W(k)=-2\pi Gk^{-2}P(\rho,k)\propto k^{-2.2}$ [97]. 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.

4 Conclusions

We explored the statistics of multiphase, magnetized, ISM turbulence and the formation mechanisms of molecular clouds with idealized numerical simulations, using the PPML solver [20]. Our periodic box models employed a large-scale solenoidal force to stir the turbulence up to the rms velocities consistent with observations. After a few dynamical times of forcing, the fully developed turbulence reached a steady state subject to our statistical analysis. A fraction of the kinetic energy supplied by the external forcing is stored in the form of induced turbulent magnetic energy, which tends to saturate at equipartition with the kinetic energy on small scales. Due to the Alfvén effect, a strong dynamic field alignment develops in the quasi-stationary state in most of the computational domain, at characteristic ISM densities $\sim 1$ cm*-3*. The models indicate that this factor may be responsible for the super-Alfvénic nature of turbulence in the cold and dense parts of molecular clouds. Being a consequence of self-organization in highly compressible MHD turbulence, this result does not depend on the way the turbulence is initiated or fed in our models. The simulations capture basic physics and help to constrain the range in parameter space, where the model is overall successful in reproducing the following observables of the local ISM, including molecular clouds:

• •

the ratio of the turbulent magnetic field component versus the regular field measured at length scales $\sim 100$ pc;

• •

the mass and volume fractions of thermally stable neutral Hi and the abundance of the unstable regimes;

• •

the variety of sonic and Alfvénic Mach number regimes in different thermal phases;

• •

the low rates of star formation per free-fall time controlled by the multiphase thermodynamics and large-scale turbulence;

• •

the lognormal distribution of column densities and deviations from lognormality for some specific tracers;

• •

the mass-weighted distribution of thermal pressure, including the mean value and the characteristic asymmetric shape;

• •

the overall hierarchical filamentary morphology of the molecular gas and the alignment of filaments with respect to magnetic field lines;

• •

the linewidth-size relationship for molecular clouds, including the slope and the offset;

• •

the ratio of solenoidal-to-compressive velocity power.

Our models also provide predictions for the shape of magnetic field PDFs, which are strongly non-Gaussian. We apply new probabilistic tools developed for self-gravitating compressible turbulence elsewhere [97] to assess the spectral energy distributions. These help to disentangle the effects of turbulence, thermodynamics, and self-gravity and single out characteristic length scales. By design, our models are limited to the local ISM and molecular clouds. Exploring cases with more active star formation, starbursts, mergers, and other non-steady phenomena would require more sophisticated experiments with properly closed feedback loops. Higher dynamic range modeling is required to further harden the conclusions.

Acknowledgments

This research was supported in part by the National Science Foundation through grants AST-0607675, AST-0808184, AST-0908740, AST-1109570, and AST-1412271; used XSEDE allocation MCA07S014, and DOE Office of Science INCITE-2009 and DD-2010 awards allocated at NCCS (ast015/ast021). A.K. and M.N. acknowledge significant contribution to this work by the late Dr. Sergey D. Ustyugov, who developed PPML and ran the simulations. Christoph T. Lee carried out statistical data analysis to produce Figure 15. The authors would like to thank both anonymous referees for constructive reports, which helped to improve the quality of presentation.