Turbulence in Zeeman Measurements from Molecular Clouds

The "textbook" argument for the B–n relation indices solely depends on gravity to concentrate gas, and the B-field can only play a role against gravity (Figure 1). In reality, however, thermal pressure should help support the cloud, and turbulence can either disperse or enhance the densities of gas and/or field lines. As a consequence, even the temporal B–n relation should not be interpreted by the textbook model before a careful examination of the model assumptions. Next, we check whether the gas flows, which shape the B–n relation, are converged only by self-gravity and hindered by B-field forces.

Only forces perpendicular to the B field, f⊥, can regulate the field strength. Positive/negative divergence of f⊥ can accelerate the gas to disperse/compress field lines, thereby decreasing/increasing the field strength. To determine what force is mainly compressing or dispersing core fields, we can integrate the divergence of f⊥ within the core; the results are shown in Figure 4, where, for simplicity, f⊥ is defined by the mean-field direction of a core.

Zoom In Zoom Out Reset image size

As expected, the gravitational force always compressed the field (with negative f_G ⊥ divergence in Figure 4), whereas the thermal pressure constantly opposed gravity with a positive divergence (note that plotted in Figure 4 is negative f_Th⊥ divergences to save space). What should be noticed is that, most of the time, the magnitude of the f_Th⊥ divergence is significantly greater than that of f_G ⊥. In other words, it is gravity that needs assistance against thermal pressure in order to hold a core together. This differs greatly from the textbook model (Figure 1), in which the thermal force is assumed to be negligible. Indeed, the divergence of the advection and B-field forces (f_Adv⊥ and f_B ⊥ in Figure 4, respectively) are negative for significant periods of time, when they confined the cores rather than supported the cores against gravity. It is worth noting that the simulation is isothermal, and the turbulence remains supersonic in the cores (Table 1). It is the density gradients set up by turbulent shocks that are responsible for the f_Th⊥. This is in agreement with the scenario that the 2/3 index is also due to the random turbulence compression. Only the densest Core 3 survived until violating the Truelove criterion (Truelove et al. 1997), while the other two cores are dynamic transient structures. B-fields play more of a stabilizing role in all cases. In Figure 4, the divergences of -f_Th⊥ and f_B ⊥ (the two green lines) share the same trend, i.e., when f_Th⊥ became more dispersive, f_B ⊥ turned more compressive, and vice versa. In addition, B-fields also reacted against the advection force, which can be observed from the fact that the higher-frequency variations of the divergences of f_Adv⊥ and f_B ⊥ are usually complementary. This is especially apparent for Core 3.

In assessing the gravitational boundness of molecular clouds, observers estimate the virial parameter (Bertoldi & McKee 1992), α = 2E_Tu /E_G (2 times the turbulence to gravitational energy ratio; see Appendix A.2 for detailed definitions of E_Tu and E_G ), which usually has a negative power-law relationship with the core mass (e.g., Kirk et al. 2017; Keown et al. 2017; Kerr et al. 2019). Cores 1–3 have α close to one, which are consistent with observations (Figure 5). Note that α ignores thermal and magnetic energy, so α <1 does not guarantee a contraction. Our simulation also replicated other observations, such as magnetic criticality = 1–2 (Table 1; Li et al. 2013, 2015; Myers & Basu 2021) and core density profiles n ∝ r^{−1.46±0.12} (Pirogov 2009; Kurono et al. 2013; Zhang et al. 2019) leading to strong pressure gradients (Figure 4). With the marginal virial parameter, magnetic criticality, and strong density gradient, cores are on the verge of contraction and expansion, as demonstrated by Cores 1 (expanding), 2 (bouncing), and 3 (contracting). This is another perspective to avoid using the ever-collapsing textbook model (Figure 1, left) to interpret observations.

Zoom In Zoom Out Reset image size

Finally, we note that turbulence may induce both ambipolar (Li & Houde 2008; Tang et al. 2018) and/or reconnection (Lazarian et al. 2012) diffusion, which are not simulated here but can potentially modify B–n relations. This study, however, is not intended to predict the exact index of B–n relations but rather aims to explore the rationale behind accessing temporal B–n relations via Zeeman measurements (spatial relations). Field diffusion should only flatten both the temporal and spatial B–n relations to some extent (see, e.g., Tsukamoto et al. 2015) instead of resolving their significant discrepancy (Figure 2).

We simulated the interior of a molecular cloud using the ideal MHD code ZEUS-MP (Hayes et al. 2006; Otto et al. 2017). Assuming isothermal, the code solves the following set of equations:

$$\begin{eqnarray*}&&\begin{array}{l}\qquad \displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot \left(\rho {\boldsymbol{v}}\right)=0\\ \rho \left(\displaystyle \frac{\partial }{\partial t}+{\rm{v}}\cdot {\boldsymbol{\nabla }}\right){\boldsymbol{v}}={\boldsymbol{J}}\times {\boldsymbol{B}}-{\boldsymbol{\nabla }}p-{\boldsymbol{\nabla }}\phi \\ \qquad \displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\boldsymbol{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}})\\ \qquad \qquad p=\rho {c}_{s}^{2}\\ \qquad \qquad {\boldsymbol{J}}={\boldsymbol{\nabla }}\times {\boldsymbol{B}}\\ \qquad \qquad {\boldsymbol{\nabla }}\cdot {\boldsymbol{B}}=0\\ \qquad \quad {{\boldsymbol{\nabla }}}^{2}\phi =4\pi G\rho ,\end{array}\end{eqnarray*}$$

where ρ and p are mass density and thermal pressure, respectively; v and B are velocity and magnetic field vector, respectively; the constant c_s ∼ 0.2 km s⁻¹ is the sound speed assuming a temperature of 10 K; and $G=4.3\times {10}^{-3}\ \mathrm{pc}\cdot {M}_{\odot }^{-1}\cdot {\left(\tfrac{\mathrm{km}}{s}\right)}^{2}$ is the gravitational constant.

The initial B field is 14.4 μG in the z -direction uniformly. The ratio between thermal pressure and magnetic pressure is β = 0.05. The uniform initial density is ρ = 1.2 × 10²¹ g cm⁻³ or n = 300 H₂/cc, assuming a mean molecular weight of 2.36. The size of the simulation domain is 4.8³ pc³, resolved by 960³ cells, and a periodic boundary condition is applied. Correspondingly, the mass of the whole domain is twice the magnetic critical mass (M_Φ = Φ/2π G^1/2, where Φ is magneticflux).

The simulation was separated into two phases. During the first phase, turbulence was driven without self-gravity. When the turbulent energy was saturated (i.e., when the turbulent energy spectrum is stable), the Mach number of turbulence (M was 5.6 (Table 1) and the rms Alfvén Mach number (M_A was 0.6. Self-gravity was turned on in the second phase, while the turbulence driving was halted. The simulation was terminated after gravity had been turned on for 2 Myr, before any core can violate the Truelove criterion (Truelove et al. 1997).

A detailed description of turbulence driving is given in Otto et al. (2017). Briefly, we set up a vector field a ( k ) in the Fourier space with a zero mean and a variance $\propto {k}^{6}{\exp }(8k/{k}_{c})$. Therefore, the power spectrum peaks sharply at k_c , which we set as 2, corresponding to a driving scale of 2.4 pc. We then apply inverse Fourier transform on vector field a ( k ) with a random phase to obtain the driving velocity field v ( x ). Details of the saturated velocity spectrum are given in (Zhang et al. 2019; see their Figure 3). In this simulation, the turbulence driving is purely solenoidal.

As stated in the main text, our goal is to study the temporal B–n relation of individual molecular cloud cores. This requires the data of a core at different evolution stages. Defining a core in a single simulation snapshot is simple. The difficulty lies in tracing this core in other snapshots due to transposition and deformation. Here, we describe how to trace a core across snapshots.

Prior to core extraction, we first define "clumps" as dense clumps are the birth beds of cores. We define clumps using potential field contour surfaces, such that, within a clump C, $\tfrac{| {E}_{G}| }{{E}_{{Tu}}+{E}_{B}+{E}_{{Th}}}\gt 1/2$, which is a relaxed gravitational bounded condition. In the equation, ${E}_{{Tu}}=\tfrac{1}{2}M{\sigma }^{2}$ indicates the turbulent energy. The energies ${E}_{B}=\tfrac{1}{8\pi }{\int }_{V}{B}^{2}{dV}$ and ${E}_{{Th}}=\tfrac{3}{2}{{Mc}}_{s}^{2}$ are magnetic field and thermal energy, respectively. Finally, the ${E}_{G}=\tfrac{1}{2}{\int }_{V}\rho \phi {dV}$ is the potential energy, where the potential field ϕ is due to all of the mass under the periodic boundary condition. We assume the global maximum of potential as the zero potential point.

In accordance with the above criterion, we scanned through potential contour surfaces in each snapshot to look for clumps. From high to low, we scanned through 40 potential values evenly distributed between zero and the minimum potential. For each potential value, the connected component labeling (CCL) method (Silversmith 2021) was applied to index individual volumes confined by the contour surfaces. The clump criterion was applied to each of these volumes, and clumps were excluded from further examination with lower potential values. Two clumps in adjacent snapshots are considered to be temporally connected if they spatially overlap. The algorithm is summarized in Figure 6 below.

Zoom In Zoom Out Reset image size

Several sets of temporally connected clumps are found. Only clumps persisting for more than 1 Myr are retained for further analysis because only long-lived clumps can host long-lived cores for temporal study. Three clumps meet the requirements.

We define a core with fixed mass as the following. The mass of a core is typically 0.1 to 10 solar masses (Mac Low & Klessen 2004); we use two solar masses in this study. For each clump, the snapshot with the highest peak density is identified, and the core is defined by the density contour surface containing this peak and two solar masses. The algorithm to identify the same core in the earlier and later snapshots is detailed in Figure 7. Although a core can transpose the position and deform the shape, the difference in adjacent snapshots is expected to be small as long as the time difference is relatively short. For this purpose, we took frequent snapshots every 0.02 Myr. The algorithm in Figure 7 requires the center of mass (COM) of a core to always fall within the effective core scale from the projected COM position defined by the position and velocity of the COM in the previous snapshot; the effective core scale is defined by 2 times of $r={\left(\tfrac{3}{4\pi }V\right)}^{1/3}$, where V is the core volume in the previous snapshot.

Zoom In Zoom Out Reset image size

A major focus of our analysis is the evolution of the B-field strength, which is driven by force components perpendicular to the local mean field $\overline{{\boldsymbol{B}}}$:

$$\begin{eqnarray*}&&{{\boldsymbol{f}}}_{i,\perp }={{\boldsymbol{f}}}_{i}-{{\boldsymbol{f}}}_{i,\parallel }={{\boldsymbol{f}}}_{i}-{{\boldsymbol{f}}}_{i}\cdot \displaystyle \frac{\overline{{\boldsymbol{B}}}}{\left|\overline{{\boldsymbol{B}}}\right|},\end{eqnarray*}$$

where i = B, Adv, Th, or G, standing for the Lorentz, advection, thermal, and gravitational forces, respectively. More precisely, they are formulated as follows:

$$\begin{eqnarray*}&&\begin{array}{l}\qquad {{\boldsymbol{f}}}_{\mathrm{Adv}}=-\rho ({\boldsymbol{v}}\cdot {\rm{\nabla }}){\boldsymbol{v}}\\ {{\boldsymbol{f}}}_{B}={\boldsymbol{J}}\times {\boldsymbol{B}}=\left({\rm{\nabla }}\times {\boldsymbol{B}}\right)\times {\boldsymbol{B}}\\ \qquad \qquad {{\boldsymbol{f}}}_{\mathrm{Th}}=-{\rm{\nabla }}P\\ \qquad \qquad {{\boldsymbol{f}}}_{G}=-\rho {\rm{\nabla }}\phi .\end{array}\end{eqnarray*}$$

To determine whether, in general, the field lines were compressed or dispersed by f _i,⊥, we integrated the divergence within D_m , the core volume above the mean density:

$$\begin{eqnarray*}&&{\int }_{{D}_{m}}{\rm{\nabla }}\cdot {{\boldsymbol{f}}}_{i,\perp }{dV}.\end{eqnarray*}$$