Anchoring Magnetic Fields in Turbulent Molecular Clouds II - from 0.1 to 0.01 parsec

We (Li et al. 2009; Paper-I) compared the magnetic field directions inferred from polarimetry data obtained from 100-pc scale inter-cloud media (ICM) and from sub-pc scale molecular cloud cores. The highly correlated result led us to conclude that cloud turbulence must be sub-Alfvenic. Here we extend the study with 0.01-pc cores observed by interferometers. The inferred field directions at this scale significantly deviate from that of the surrounding ICM. An obvious question to ask is whether this high-resolution result contradicts the sub-Alfvenic picture concluded earlier. We performed MHD simulations of a slightly super-critical (magnetic criticality = 2) clouds with Alfvenic Mach number $M_A = 0.63$, which can reproduce the Paper-I results, and observed the development towards smaller scales. Interestingly, all subregions hosting cores with $n_H$$_2$>$10^{5}$/cc (the typical density observed by interferometers) possess $M_A = 2-3$. Not too surprisingly, these slightly super-Alfvenic cores result in B-field orientation offsets comparable to the interferometer observations. The result suggests that gravity can concentrate (and maybe also contribute to, which takes more study to confirm) turbulent energy and create slightly super-Alfvenic cores out from sub-Alfvenic clouds. The results of our simulations also agree with the observed velocity-scale (Kauffmann et al. 2013), mass-scale (Lombardi et al. 2010) and field strength-density (Li et al. 2015; Crutcher et al. 2010) relations.


Introduction
The critical role played by turbulence in star formation is solidly established in the past two decades; see, for example, the review by McKee & Ostriker (2007). The review noted that cloud magnetic fields (B-fields) are "not too weak, however" based on the ordered B-fields mapped from three clouds . We (Paper-I) further confirmed this picture of ordered cloud fields by showing the field-direction correlation between 25 cloud cores (1-0.1 pc) and their surrounding inter cloud media (ICM ~ 100 pc). The offsets between the core and ICM fields are shown in Fig. 1 (black). Afterwards, more evidence of dynamically important cloud B-fields have been detected; see our PPVI review .
Recent interferometer (SMA and CARMA) surveys showed that the core-ICM correlation in Bfield direction seems not to present within small cores at ~ 0.01 pc (Zhang et al. 2014;Hull et al. 2014). This raises several questions: (1) Is the situation the same for the cores collected in Paper-I? (2) If so, what does the situation imply? Does it mean that the turbulence changes from sub-Alfvenic, the conclusion from Paper-I, to super-Alfvenic at smaller scales? To answer question (1), we have collected data from Zhang et al. (2014) and Hull et al. (2014) (Z/H14 hereafter) for those cores in Paper-I and the result is presented in section 2, Table I and Fig. 1. To answer question (2), we have performed sub-Alfvenic molecular cloud simulations to first reproduce the result from Paper-I (above 0.1 pc) and then observe the circumstance at scales 0.1-0.01 pc. The results of the simulations are summarized in section 3, which are compared with observations in section 4.

Fig. 1
Core/ICM B-field offsets against the spatial scales of cores. The black symbols are from the data in Paper-I, and the red symbols are the interferometer data (Z/H14). Data from Paper-I with no interferometer counterparts are shown as hollow circles. The error bars indicate the interquartile ranges (IQRs) of the B-field orientation distributions. Above ~ 0.1 pc, most of the offsets are smaller than 45 degrees and the upper envelope grows with decreasing scale. Below 0.1 pc, there seems no direction preference at all.

Simulation setup
The simulations are isothermal (10 K) and start with a uniform density (1.20 × 10 -21 g/cm 3 or 300 H2/cc assuming a mean molecular weight of 2.36) and a uniform B-field (14.4 µG) over a cube with a linear size of 4.8 pc, which is resolved evenly by 960 pixels. The boundary condition is periodic. The ratio of the total mass to the magnetic critical mass (Φ/2πG 1/2 , where Φ is magnetic flux), a.k.a. magnetic criticality, is 2. The simulations proceed in two stages. First, self-gravity is turned off, and the pure solenoidal turbulence is driven at 2.4 pc until the turbulent energy power spectrum is stable. At this point, the sonic Mach number is 5.7 and the overall Alfvenic Mach number, MA≣< V/vA>, as conventionally defined in the literature (e.g. Falceta-Gonçalves et al. 2008, Burkhart et al. 2009), is 0.74. This MA, however, is not observable, as both the Chandrasechar-Fermi and Zeeman methods measure the mean field strength of a targeted volume. To connect simulations with observations, here we define MA_obs = V/VA , where V and VA are similar to V and vA but the density weighted 3D velocity dispersion and Alfven velocity based on the density weighted mean field strength of the whole volume; MA_obs = 0.92. The viral parameter, 5 V 2 L/6GM, is 0.51, were M and L are, respectively, the total mass and total volume scale. During the second stage, self-gravity is turned on while the turbulence driving continues for another 1 Myr. At the end, the velocity spectrum remains Kolmogorov-type ( Fig. 3) with an overall MA_obs = 0.84 and MA = 0.63. More details about the turbulence driving can be found in the Appendix of Otto, Ji & Li (2017). The Jeans length is always resolved by at least 8 simulation pixels (Truelove et al. 1997). We repeated the simulation three times (Cube 1, 2 and 3; Table II) with different random seeds for turbulence generation.
As our goal is to understand observations, we will mainly use MA_obs hereafter. Note that, conventionally, only the initial MA of a simulation is reported in the literature (e.g., Li et al. 2015b, Mocz et al. 2017, but an initial condition is not comparable with observations. The energy redistributes between gravity, turbulence and B-fields as time goes on. Here we report MA_obs one million years after gravity is turned on, which is in the same order as the typical cloud age.

Overall field-density relation
As the first look of the simulation result, we simply survey the offset between local field orientation and the initial field direction and study how is this offset related to density (Fig. 2,upper panel). For "local" field orientation, we divided the 960 3 -pixel cube evenly into 30 3 -pixel subcubes and calculated the density weighted mean field direction of each subcube. The reason to use the 30pixel scale for a subcube is that the turbulent energy is artificially dissipated below 20-pixel scale ( Fig. 3; which is typical for numerical simulation; see Kritsuk et al. 2011), where the field structures due to turbulence are thus underestimated. The density-offset relation of the (960/30) 3 subregions of Cube 1 is shown in Fig. 2 (upper panel), where the upper-envelope of the offset systematically increases with the density. This is the first indication that our simulations might be able to explain Fig. 1, as the upper envelope of Fig. 1 also increases with decreasing scales, which implies increasing density.
Also shown in Fig. 2 (lower panel) is the field strength against density at the pixel scale in Cube 1. The slope of the upper envelop gradually changes from ~ 0.03 at the lower density (nH2 < Fig. 2 B-field versus density in Cube 1. Upper panel: Cube 1 is evenly divided into subcubes of 30 3 pixels, and each data point shows the mean density and the density weighted mean field offset of one subcube. Lower panel: The density is binned into groups with an even width of 200 H2/cc. For all the pixels with densities within one bin, their field strengths are averaged. This is a plot of the averaged field strength versus the central density of each bin. The reference slopes of 0.03 and 0.67 are shown as the grey solid lines.

Field-density relation of high-density cores
To study how the field can be further deviated in even higher densities, we need to zoom-in into even smaller scales, keeping in mind that once getting smaller than the 20-pixel scale, the field offset from our simulations should be treated as the lower limit because the turbulent energy is artificially dissipated (Fig. 3).
To ensure a density range that covers both single-dish and interferometer data, we first identify cores with nH2 > 10 5 /cc in our simulations; five cores reach this density. For each of these cores and their surrounding regions, contour surfaces of a series of nH2, as shown in Table II, are identified. The scale of the volume enclosed by each contour surface is estimated by the cube root of the volume. The mean field direction of each volume is calculated by the density-weighted mean of the field directions at all the pixels within the contour surface. This way, within a core region, density increases with decreasing scale, which imitates the observations in Paper-I and Z/H14. The result of core field orientations is shown in Table II and Fig. 3, where the 3D offsets are measured from the direction of the uniform field in the initial condition. Fig. 3 already looks inspirational, because of the large offset angles at smaller scales. Above 0.1 pc, which are mainly the scales probed by Paper-I, the angles are all smaller than 45 degrees. Below 0.1 pc, the offset increases but still within 90 degrees. Again, due to the unrealistic turbulence energy dissipation below the 20-pixel scale, the offsets from below 0.1 pc are just lower limits. Mocz et al. (2017) performed similar simulations with a much higher resolution such that artificial energy dissipation only happens at scales much smaller than 0.01 pc and still they found offsets not excessing 90 degree at 0.01 pc for MA = 1.2 and even 3.5. However, note that their MA is defined by the "input" kinetic energy and the initial B-field strength; at the snapshot they measured offsets, the MA should be significantly lower (see section 3.1).
Note that the B-nH2 relation "within each core" in Table II is much shallower than a power-law with index 0.67 (Crutcher et al. 2010). The index here (~ 0.3) is closer to the one observed by Li et al. (2015). Indeed, the B-nH2 relation of each core describes the condition within "affiliated structures", which is comparable to the observation carried out by Li et al. (2015), but different from Crutcher et al. (2010), where the B-nH2 relation is mostly between independent structures. Moreover, the cores in Table II extend to density as low as nH2 ~ 10 3.5 ; if one fits the upper envelope of nH2 > 10 3.5 in Fig. 2 (lower panel) with only one power low, the index will be around 0.3. The B-nH2 relation in Fig. 2 agrees with Li et al. (2015b) and Mocz et al. (2017) on that the 0.67 index is only for nH2 well above 10 4 . In other words, if one power low is fitted to a density range involving nH2 < 10 4 , the index has to be significantly less than 0.67. Yet the 0.67 index from Crutcher et al. (2010) covers nH2 > 300! This discrepancy between Zeeman measurements and simulations motivated us to revisit the Bayesian analysis of the Zeeman data (Jiang, Li & Fan 2019) and found that the index and the threshold density cannot be reliably derived from the data with large uncertainties in nH2.

Comparison between simulations and observations
Our simulations (Fig. 3) can be compared with the polarimetry observations ( Fig. 1) in two different ways. In both cases, we need to get the 2D projections of the field angles. Moreover, if a projected angle is obtuse, the supplementary angle should be adopted to imitate the fact that polarization vectors have the 180-degree ambiguity (headless), and the acute angle is utilized to describe the offset between two polarization vectors (Paper-I; Z/H14).

The overall trend of the field projections
We bin the 3D angles into three groups -"above 1 pc", "below 0.1 pc" and "in between". Each 3D angle in Fig. 3 is projected to 135 directions evenly distributed on a sphere centered at the angle. For each bin, all the projected angles are collected (with obtuse angles replaced by their supplementary ones), and their distribution is also shown in Fig. 3. The distributions are comparable to Fig. 1: above 0.1 pc, 90% of the projections are below 45 degrees. Below 0.1 pc scale, even though the 3D offsets here are only the lower limits, the distribution is already close to uniform within 0-90 degrees. In other words, the large field deviations (relative to the local ICM field orientations) observed by interferometers (Fig. 1) do not contradict the idea of an overall sub-Alfvenic cloud (MA = 0.63).

Fig. 3
The results of MHD simulations. The diamonds show the 3D core field offsets form the initial condition (left axis) against the core sizes (bottom axis; see Field-density relation of highdensity cores in section 3.2 for the definition of a core size). Each of these 3D offsets is projected along 135 directions evenly distributed over a 4π sr solid angle. The scales are binned into "above 1 pc", "below 0.1 pc" and "in between". The histograms present the distribution of the projected 2D offsets in each bin; the upper axis is the fraction. Note that the 2D histograms are "folded" at 90 degrees to emulate the disability of polarimetry observations on distinguishing supplementary offset angles. Also plotted is the turbulent 3D velocity spectrum (the dark blue line; the velocity is shown in the right axis) of Cube 1, which peaks at 2.4 pc, the turbulence driving scale, and is artificially dissipated below 0.1 pc. The dotted line (y=0.90 x 0.31 ) is a fit to the velocity spectrum between 0.1 and 1.25 pc.

Reproducing the field directions of affiliated structures
Besides the overall trend, here we show that the simulation results can also reproduce the observed multi-scale field orientations of individual cloud systems. Orion molecular cloud and Cube 1 are used as an illustration, for that they have more cores (OMC-1, 2 and 3; simulated core 1, 2 and 3) and each OMC core has at least three submillimeter observations at different scales (Table  I).
In Fig. 4, OMC-1, 2 and 3 are isolated, respectively, from Fig. 1. The smallest structures are OMC-1 KL NW/SE, OMC-2-FIR3, and OMC-3-MMS5/6 ( Table I). For each simulated core, we surveyed all the possible projections. If a projection falls within the error bar of any of the aforementioned five smallest cores and the error bars of its parental structures, we plot this particular projection in Fig. 4. Again, even with the underestimated field structures below the 20pixel scale, there is no problem for the simulation to reproduce the observed field variation over different scales.

The emergence of super-Alfvenic cores
Our MHD simulations can closely reproduce the field offsets from not only Paper-I but also Z/ H14. To further explore the reason of the field offsets and answer question (2), we study the relation between nH2 and MA. The subcubes under study are evenly distributed within Cube 1 with the smallest separation between two subcube centers equals to 100 pixels (thus there are overlaps between larger subcubes). The result is given in Fig. 5., where the MA_obs of core 1, 2 and 3 at ~100and ~150-pixel scales are also displayed. While the simulation cube is overall sub-Alfvenic, it can produce dense and super-Alfvenic subcubes and cores! Which is most likely due to the mass (and thus kinetic energy) concentration along the B-fields, that results in local kinetic (turbulent) energy enhancement without magnetic field compression. This anisotropic accretion channeled by B-fields is also observed in other simulations and discussed by Otto, Ji & Li (2017). Also, part of the gravitational potential energy may be converted into turbulent energy after the accretion (e.g. Heitsch 2013).
Note that the cores are only slightly super-Alfvenic with MA_obs ~ 3 ( Fig. 5; Table II). This level of kinetic energy is not enough to randomize the field orientation. The subcubes with nH2 ~ 10 3.5 in Fig. 2 (upper penal) possess MA_obs ~ 3 (Fig. 5), but their field direction offsets seldom go beyond 90 degrees. It will take an MA an order of magnitude higher to deviate the field to close to 180 degrees, i.e. to randomize the field, at 0.01-pc scale; see, for example, Fig. 3 of Mocz et al. (2017). Indeed, the 3D offsets shown in Fig. 3 are all within 90 degrees, which is far from random (evenly distributed between 0 and 180 degrees). Note that, however, polarization vectors, which have the ambiguity of 180 degrees, are not able to distinguish between randomness and 3D angles distributed within 0-90 degrees. Thus one should not conclude a random field based merely on Z/H14 data (red in Fig. 1).

Polarization holes
The trend of increasing field offset with decreasing scale (and thus with increasing density) explains the so-called "polarization holes" -the tendency of the decreasing fraction of submm polarization with growing column density (N). In Paper-I, we proposed that polarization holes can occur naturally due to more B-field structures along a line-of-sight with higher N, not necessarily due to the lower grain alignment efficiency in high-density regions as suggested in the literature (e.g. Podoan et al. 2001). The fact that Z/H14 see higher polarization fraction than Paper-I (Tang 2016) supports our proposal. Higher N usually implies a line-of-sight going through a higher density and thus more B-field direction dispersion can be expected (Fig. 2). When these richer field structures are not resolved, they will appear as a lower polarization fraction.

Summary
We extend the multi-scale study of B-field directions in Paper-I down to 0.01-pc scale using the data from Z/H14. The multi-scale field correlation decreases for scales below 0.1 pc. The directional correlation between 100-pc and 0.1-pc fields (Paper-I) has been understood as that a molecular cloud, as a whole, should be trans-or sub-Alfvenic. Our MHD simulations show that dense cores developing from sub-Alfvenic clouds (MA_obs = 0.84 and MA = 0.63) can be slightly super-Alfvenic (MA_obs ~ 3) to significantly deviate the core fields and thus explain the field offsets observed by interferometers (Z/H14), but not enough to randomize the field directions. The simulation results are also consistent with other observations, e.g., B ∝ n 0.32±0.08 for affiliated structures (Li et al. 2015), B ∝ n 0.67 for independent high-density structures (Crutcher et al. 2010), V ∝ l 0.31 , the linewidth-size relation (Kauffmann et al. 2013) and M ∝ R 1.54±0.12 , the core mass-size relation for core size R ∈ [0.1 1] pc (Lombardi et al. 2010).
Understanding observations through numerical simulations had become a common practice in modern astronomy. However, we note that MA ≡ < V/vA> in the initial condition, the conventional description of turbulence/B-field relative importance in a simulation, is not comparable with an observed MA. First, after the turbulence spectrum is stabilized in a simulation, the overall MA should be significantly lower than the initial condition (section 3.1). This is because part of the kinetic energy is converted to random B-field energy and contributed to the denominator of MA. Second, the observable, MA_obs ≡ V/VA, the ratio between the density weighted velocity dispersion and the Alfven velocity of the density weighted mean field, will not be affected by random B-fields and thus higher than the MA ≡ < V/vA> (section 3.1). Third, the localized MA, e.g. MA of a cloud core, can be significantly different from the overall value of the cloud (Fig. 5).
Finally, polarization offsets distributed over 0-90 degrees (the red data in Fig. 1) can be interpreted in two possible ways: (1) the 3D B-field offsets are simply within 90 degrees; (2) the 3D B-fields offsets are over 0-180 degrees, i.e. random, but polarization offsets cannot go beyond 90 degrees due to the headless nature of polarization orientations. The interpretation (2), however, requires the turbulence to be highly super-Alfvenic (e.g., MA=35 in Fig. 3 of Mocz et al. 2017), but Li et al. (2015b) showed that simulations with MA=10 already conflict with the observations in Paper-I (the black data in Fig. 1).
% Fig. 4 The possibility for our simulated sub-Alfvenic cloud (Cube 1) to reproduce OMC field structures (black). Over a scale range of 0.5 pc, the OMC field orientation can change by almost 50 degrees after projected on the plane of sky. Here we survey all the possible projections of core 1-3 in Cube 1 and show that the reproduction of OMC fields is possible. For any affiliated structures (e.g. Orion KL-SE ⊂ Orion KL ⊂ OMC-1), we can always find a simulated core (red, green or blue) with a proper projection to lay the projected field orientations within the error bars from any scale of the affiliated structures.   Table II. Properties of the cloud cores resulted from our sub-Alfvenic MHD simulations * MA_obs drops when the scale goes below 0.1 pc for all cores. This is due to the limitation of our MHD simulations, which have to artificially damp the turbulence below 0.1 pc; see Fig. 3. # The density (n) -scale (l) relation derived from these five cores is n ∝ l -0.98±0.25 and the core mass ~ (165.1±25.3)×l 1.54±0.12 M ☉ ✝ The field strength -density relation derived from these five cores is B ∝ n 0.32±0.08 ; see the discussion in section 3.2