# Calibrating the Davis-Chandrasekhar-Fermi Method with Numerical Simulations: Uncertainties in Estimating the Magnetic Field Strength from Statistics of Field Orientations

Junhao Liu (刘峻豪)¹,²,³,⁴✉, Qizhou Zhang³✉, Benoît Commerçon⁵, Valeska Valdivia⁶, Anaëlle Maury³,⁶, and Keping Qiu¹,²✉

¹School of Astronomy and Space Science, Nanjing University, 163 Xianlin Avenue, Nanjing, Jiangsu 210023, People's Republic of China; kpqiu@nju.edu.cn, liujh@smail.nju.edu.cn

²Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing, Jiangsu 210023, People's Republic of China
³Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

⁴East Asian Observatory, 660 N. A'ohōku Place, University Park, Hilo, HI 96720, USA; liujunhao42@outlook.com

⁵Univ Lyon, Ens de Lyon, Univ Lyon1, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69007, Lyon, France

⁶Laboratoire AIM, Paris-Saclay, CEA/IRFU/SAP—CNRS—Université Paris Diderot, F-91191 Gif-sur-Yvette Cedex, France

Received 2021 January 20; revised 2021 May 26; accepted 2021 June 17; published 2021 September 28

## Abstract

The Davis-Chandrasekhar-Fermi (DCF) method is widely used to indirectly estimate the magnetic field strength from the plane-of-sky field orientation. In this work, we present a set of 3D MHD simulations and synthetic polarization images using radiative transfer of clustered massive star-forming regions. We apply the DCF method to the synthetic polarization maps to investigate its reliability in high-density molecular clumps and dense cores where self-gravity is significant. We investigate the validity of the assumptions of the DCF method step by step and compare the model and estimated field strength to derive the correction factors for the estimated uniform and total (rms) magnetic field strength at clump and core scales. The correction factors in different situations are cataloged. We find the DCF method works well in strong-field cases. However, the magnetic field strength in weak-field cases could be significantly overestimated by the DCF method when the turbulent magnetic energy is smaller than the turbulent kinetic energy. We investigate the accuracy of the angular dispersion function (ADF; a modified DCF method) method on the effects that may affect the measured angular dispersion and find that the ADF method correctly accounts for the ordered field structure, the beam smoothing, and the interferometric filtering, but may not be applicable to account for the signal integration along the line of sight in most cases. Our results suggest that the DCF methods should be avoided to be applied below ~0.1 pc scales if the effect of line-of-sight signal integration is not properly addressed.

Unified Astronomy Thesaurus concepts: Magnetic fields (994); Radiative transfer simulations (1967); Magnetohydrodynamical simulations (1966); Star formation (1569)

## 1. Introduction

Magnetic fields play an important role in the star formation process (Crutcher 2012). It may provide support against gravitational collapse (Shu et al. 1987), prevent or regulate the fragmentation of dense clumps<fnref n="7" /> and cores (Commerçon et al. 2011; Zhang et al. 2014), and affect the isotropy of turbulence (Goldreich & Sridhar 1995). However, the relative importance of the magnetic field compared to the turbulence and the gravity in star formation is still uncertain (Crutcher 2012). Correctly estimating the strength of the magnetic field is crucial to compare the magnetic energy with the turbulent energy and the gravitational energy and to distinguish between star formation theories with a strong field (Mouschovias et al. 2006) or with a weak field (Mac Low & Klessen 2004).

There are several ways to observationally determine the magnetic field strength in star-forming regions. Zeeman splitting observations of spectral lines can provide information on the direction and strength of the line-of-sight magnetic field component and is the only technique to directly measure the magnetic field strength (Crutcher & Kemball 2019). However, Zeeman splitting observations require particular line tracers and are very time consuming if the line is weak. The recently developed velocity gradient technique proposes to trace the

plane-of-sky magnetic field orientations from spectral line observations with the assumption that the velocity gradient is perpendicular to the magnetic field (González-Casanova & Lazarian 2017). On the other hand, polarized thermal dust emissions at millimeter/submillimeter wavelengths have been proven to be a promising way to trace the orientation of the plane-of-sky magnetic field if the shortest axis of irregular dust grains is aligned with magnetic field lines (Davis & Greenstein 1949; Hildebrand 1988; Lazarian 2007; Lazarian & Hoang 2007; Andersson et al. 2015). Once the plane-of-sky magnetic field orientation is obtained for the observed region either from dust polarization or velocity gradient of molecular lines, the Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951; Chandrasekhar & Fermi 1953) or the polarization-intensity gradient method (Koch et al. 2012) can be applied to indirectly estimate the magnetic field strength.

Specifically, the DCF method and its modified forms have been the most widely used method to determine the magnetic field strength in star-forming regions. The method is based on the assumptions that there is a prominent mean-field component, the perturbation on the mean-field component is due to turbulence, and the magnetic field perturbation can be traced by angular dispersions on polarization maps. It also assumes that the turbulence is isotropic and there is equipartition between the turbulent components of the kinetic energy and the magnetic energy. With these assumptions, the magnetic field strength can

<footnote n="7">⁷Following the nomenclature in Zhang et al. (2009), we refer to a molecular clump as an entity of ~1 pc and a dense core as an entity of ~0.1 pc.</footnote>
be derived from the DCF method once the angular dispersion on polarization maps, the turbulent velocity dispersion, and the gas density are estimated.

The DCF methods are subject to uncertainties from various effects that may affect the measured angular dispersion and lead to underestimation or overestimation of the magnetic field strength. These effects include the signal integration and averaging along the line of sight (Zweibel 1990; Myers & Goodman 1991), the contribution from the turbulent magnetic field to the angular dispersion at scales smaller than the turbulent correlation scale (Hildebrand et al. 2009; Houde et al. 2009), and the contribution from the ordered magnetic field to the angular dispersion (Girart et al. 2006; Falceta-Gonçalves et al. 2008; Hildebrand et al. 2009; Pillai et al. 2015; Pattle et al. 2017). There are also additional observational uncertainties from the beam-smoothing effect (Heitsch et al. 2001; Falceta-Gonçalves et al. 2008; Houde et al. 2009) and the interferometric filtering effect (Houde et al. 2016). Moreover, nonsatisfaction of the assumptions of the DCF method, i.e., the turbulence is not isotropic, the angular dispersion does not correctly trace the perturbation of the magnetic field, or the energy is not in equipartition, could also lead to inaccurate results.

There were several numerical studies that incorporated one or several of the above-mentioned effects (Heitsch et al. 2001; Ostriker et al. 2001; Padoan et al. 2001; Falceta-Gonçalves et al. 2008; Skalidis & Tassis 2021). However, none of these simulations have explored physical conditions comparable to high-density and small-scale star-forming regions where self-gravity is dominant. On the other hand, the angular dispersion function method (ADF; Falceta-Gonçalves et al. 2008; Hildebrand et al. 2009; Houde et al. 2009, 2016) analytically takes into account all of the above-mentioned effects that affect the measured angular dispersion but has not been tested in comparison with numerical simulations. Thus, it is of great importance to test the raw DCF method and the ADF method with synthetic observations of small-scale self-gravitating MHD simulations.

In this work we perform high-resolution MHD simulations and radiative transfer simulations of massive star-forming clusters. We study the accuracy of the assumptions of the DCF method step by step and apply the raw DCF method and the ADF method to synthetic observations of the simulations. We describe our simulations in Section 2. We present the results of the simulation and the application of the DCF method in Section 3. In Section 4 we discuss the uncertainties of the DCF method and compare our work with previous simulation works. Section 5 is a summary of this work.

## 2. Simulations

We adapt several simulations from Fontani et al. (2018) and perform a set of new 3D ideal radiative magnetohydrodynamic (MHD) numerical simulations to simulate the gravitational collapse of massive star-forming regions with the RAMSES code (Teyssier 2002; Fromang et al. 2006; Commerçon et al. 2011, 2014). RAMSES uses a second-order finite-volume scheme to integrate the MHD equations, combined with a constrained transport scheme for the induction equation (Fromang et al. 2006), and an implicit solver for radiative transfer (Commerçon et al. 2011, 2014). The initial conditions of our simulations are very similar to those in Hennebelle et al. (2011), Commerçon et al. (2011), and Fontani et al. (2016). We consider dense spherical structures with an initial mass of

300 $M_{\odot}$ and an initial radius of $r_0$. The initial Plummer-like density profile is given by $\rho(r) = \rho_c/(1 + (r/r_c)^2)$, where $\rho_c$ is the central density and $r_c$ is the extent of the central plateau. The density contrast between the center and the border of the structure is 10. The initial magnetic field is along the x-axis, and its strength is proportional to the column density. The initial magnetic critical parameter $\mu$, which is defined as the mass-to-flux ratio over its critical value (Mouschovias & Spitzer 1976), ranges from 1.2 to 200 in our models. We include initial (super)sonic turbulence by seeding an initial turbulent velocity field with a random-phase Kolmogorov-like power spectrum (see Hennebelle et al. 2011). We vary the initial turbulence level by changing the amplitude of the velocity dispersion, but we use the same velocity perturbation power spectrum at initialization. The magnetic fields and density are not initialized with perturbations, but the onset of turbulence and gravitational collapse modifies their distribution. We allow the turbulence to decay (the initial freefall time is always shorter than the turbulent crossing time). To follow the gravitational collapse, we employ adaptive mesh refinement based on the local Jeans length (see Truelove et al. 1997). We use a Jeans criterion that ensures mass greater than the Jeans mass is always resolved by at least eight cells. We refine the mesh until we reach our maximum resolution, which corresponds to level 14 in RAMSES. Depending on the initial thermal support (either set with initial temperature $T=10$ K or 20 K), the maximum resolution achieved is either 12.5 au or 25 au (as indicated in Table 1). When the maximum level of refinement is reached, sink particles (i.e., protostars) are introduced to avoid further refinement of the grid to follow the collapse. We use the implementation presented in Bleuler & Teyssier (2014), with slight modifications on the creation and accretion schemes. The sink particles accrete the gas that sits in their accretion volumes (sphere of radius four cells at the maximum level of refinement) and that is Jeans unstable. We consider that half of the mass accreted into the sink particles actually goes into stellar material. The luminosity, radius, and effective temperature of the protostars are then computed using the mass–radius and mass–luminosity empirical relations of main-sequence stars (e.g., Weiss et al. 2004). This protostellar luminosity is then injected within the accretion volume in the computational domain as a source term (e.g., Krumholz et al. 2009). We do not account for accretion luminosity. The initial physical parameters of different simulation models are listed in Table 1. The physical parameters would evolve temporally and locally and deviate from the initial condition during the simulation. We analyze the simulations at two different time snapshots: when the first sink forms and when the star formation efficiency (SFE) reaches 15%.

Based on the results of the MHD simulations, we perform radiative transfer modelings and recompute the temperature using the POLARIS code (Reissl et al. 2016). The radiative torque (RAT) alignment theory is assumed to be the dominant mechanism for the dust grain alignment. The sinks as in the RAMSES output are set as the central heating sources of the POLARIS input. Other setups of the radiative transfer modelings are identical to those used in Valdivia et al. (2019). An interstellar radiation field of strength $G_0 = 1$ (Mathis et al. 1983) is included as the external radiation field. The gas-to-dust ratio is assumed to be 100, and the dust is composed of 62.5% astronomical silicates and 37.2% graphite grains (Mathis et al. 1977). The dust size distribution is described by
<table><caption>**Table 1**
Initial Parameters of the MHD Models</caption><thead><tr><th>Model</th><th>T<sup>a</sup><br>(K)</th><th>Resolution<br>(au)</th><th>M<sup>b</sup></th><th>M<sub>A,peak</sub><sup>c</sup></th><th>μ</th><th>(10<sup>-19</sup> g cm<sup>-3</sup>)</th><th>r<sub>c</sub><br>(pc)</th><th>r<sub>0</sub><br>(pc)</th></tr></thead><tbody><tr><td>T10M1MU1</td><td>10</td><td>25</td><td>1</td><td>0.13</td><td>1.2</td><td>1.9</td><td>0.17</td><td>0.5</td></tr><tr><td>T10M3MU1</td><td>10</td><td>25</td><td>3</td><td>0.38</td><td>1.2</td><td>1.9</td><td>0.17</td><td>0.5</td></tr><tr><td>T10M3MU20</td><td>10</td><td>25</td><td>3</td><td>3.2</td><td>20</td><td>1.9</td><td>0.17</td><td>0.5</td></tr><tr><td>T10M6MU2<sup>d</sup></td><td>10</td><td>25</td><td>6.4</td><td>1.44</td><td>2</td><td>1.9</td><td>0.17</td><td>0.5</td></tr><tr><td>T10M6MU200<sup>d</sup></td><td>10</td><td>25</td><td>6.4</td><td>163</td><td>200</td><td>1.9</td><td>0.17</td><td>0.5</td></tr><tr><td>T20M3MU2<sup>d</sup></td><td>20</td><td>12.5</td><td>3</td><td>0.14</td><td>2</td><td>15</td><td>0.085</td><td>0.25</td></tr><tr><td>T20M3MU5</td><td>20</td><td>12.5</td><td>3</td><td>1.68</td><td>5</td><td>15</td><td>0.085</td><td>0.25</td></tr><tr><td>T20M3MU200<sup>d</sup></td><td>20</td><td>12.5</td><td>3</td><td>63</td><td>200</td><td>15</td><td>0.085</td><td>0.25</td></tr><tr><td>T20M6MU2<sup>d</sup></td><td>20</td><td>12.5</td><td>6.4</td><td>1.4</td><td>2</td><td>15</td><td>0.085</td><td>0.25</td></tr><tr><td>T20M6MU5</td><td>20</td><td>12.5</td><td>6.4</td><td>3.6</td><td>5</td><td>15</td><td>0.085</td><td>0.25</td></tr><tr><td>T20M6MU200<sup>d</sup></td><td>20</td><td>12.5</td><td>6.4</td><td>144</td><td>200</td><td>15</td><td>0.085</td><td>0.25</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Initial gas and dust temperature.

<sup>b</sup> Initial Mach number, which is the ratio between the turbulent and thermal velocity dispersions.

<sup>c</sup> Initial peak Alfvén Mach number, which is the ratio between the nonthermal velocity dispersion of the whole simulation and the Alfvén velocity of the center magnetic field value (maximum amplitude). Because the initial amplitude of the magnetic field scales with the column density, the initial Alfvén velocity is not uniform throughout the cloud. So, we only present the value of the initial peak Alfvén Mach number, which is the minimum value of the initial Alfvén Mach number.

<sup>d</sup> Adapted from Fontani et al. (2018).

$dn(a) \propto a^{-3.5}da$ (Mathis et al. 1977), where $a$ is the effective radius and $n(a)$ is the number of dust grains. The minimum and maximum of the grain size are set to be 5 nm and 1 µm, respectively.

Based on the results of the radiative transfer modelings, we produce synthetic 1.3 mm Stokes I, Q, and U maps of the simulated dust emission observed at three orthogonal directions (corresponding to the xy, zy, and xz planes) at a distance of 2 kpc with POLARIS. We postprocess these synthetic I, Q, and U images through the CASA tasks simobserve and simanalyze to produce synthetic interferometer observations. We adopt the ALMA Cycle 7 configurations<fnref n="8" /> (Configuration 1–5, hereafter AC1–AC5) for these CASA tasks. Table 2 lists the synthesized beam size and the maximum recoverable scale of the ALMA configurations. Thermal noises are not considered in these synthetic ALMA observations. The polarized intensity $PI$, polarization percentage $P$, and polarization position angle $\theta$ are estimated with $PI = \sqrt{Q^2 + U^2}$, $P = PI/I$, and $\theta = \frac{1}{2} \arctan(U/Q)$, respectively.

### 3. Results

#### 3.1. Results of Simulations

Our simulation models are state of the art and accurately reproduce the key outcomes of massive star formation. In Appendix A we show examples of the synthetic images for a model initially with a strong magnetic field and weak turbulence (T10M1MU1) and a model initially with a weak magnetic field and strong turbulence (T10M6MU200) at the time when SFE = 15%. In this work we focus on the ability to measure the magnetic field strength using the dust polarization signal as observed from realistic 3D simulation models. We refrain from analyzing the difference in magnetic field morphology among different models in this work and will present them in detail in a separate study. Here we briefly describe the outcome of the simulations.

<table><caption>**Table 2**
Beam Size and Maximum Recoverable Scale of Different ALMA Configurations at 1.3 mm</caption><thead><tr><th>Configuration</th><th>Synthesized<br>Beam (")</th><th>Maximum Recoverable<br>Scale (")</th></tr></thead><tbody><tr><td>AC1</td><td>1.5</td><td>12.6</td></tr><tr><td>AC2</td><td>1.0</td><td>9.8</td></tr><tr><td>AC3</td><td>0.61</td><td>7.0</td></tr><tr><td>AC4</td><td>0.40</td><td>4.9</td></tr><tr><td>AC5</td><td>0.23</td><td>2.9</td></tr></tbody></table>

To investigate how the MHD simulations are affected by the initial model parameter, we estimate the physical parameters within spheres (indicated by superscript “sph”) of different radii with respect to the most massive sink for all simulation models for comparison. For the volume density $n_{H_2}$, we adopt the volume-weighted value in the space concerned. Because the kinetic energy is proportional to the gas mass, we adopt the mass-weighted standard deviation of the velocity within the space concerned as the total kinetic velocity dispersion of the system: $\delta v = (v_{rms}^2 - \bar{v}^2)^{1/2}$, where $v_{rms} = (\sum(m^{cell}(v^{cell})^2)/\sum n^{cell})^{1/2}$ is the mass-weighted rms velocity, $m^{cell}$ is the mass of a cell, $v^{cell}$ is the velocity of a cell, and $\bar{v}$ is the mass-weighted mean velocity. The equivalent total (i.e., rms) magnetic field strength $B^{tot}$ within the space concerned is estimated as $B^{tot} = (\sum(V^{cell}(B^{cell})^2)/\sum V^{cell})^{1/2}$, where $B^{cell}$ is the total magnetic field strength of a cell and $V^{cell}$ is the volume of a cell. We also estimate the modified magnetic critical parameter as $\lambda_{mod} = M/M_B^{tot}$, where M is the total mass (including gas mass and sink mass) in the space concerned, and $M_B^{tot}$ is the total magnetic virial (critical) mass (Liu et al. 2020). $M_B^{tot}$ is given by

$$M_B^{tot} = \frac{\pi R^2 B^{tot}}{\sqrt{\frac{3(3-i)}{2(5-2i)} \mu_0 \pi G}}, \quad (1)$$

where R is the radius, $\mu_0$ is the permeability of vacuum, G is the gravitational constant, and i is the power-law index of the

<footnote n="8">⁸ https://almascience.eso.org/observing/prior-cycle-observing-and-configuration-schedule</footnote>
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Left: the relation between the 3D total magnetic field strength and the volume density. Right: the relation between the 3D total magnetic field strength and the magnetic critical parameter. The estimates are within different radii in all simulation models at two different time snapshots. The magnetic field strength is the volume-weighted value. Different colors correspond to different radii. Different symbols represent different initial $\mu$ values. Larger $\mu$ values correspond to weaker magnetic fields.</figcaption></figure>

<figure><img src="image_2.png" /></figure>

radial density profile. *i* is derived by fitting the radial density profile with a power law $\rho \propto r^{-i}$. We note that principally the magnetic critical parameter could be directly derived from the ratio between the gravitational energy and the total magnetic energy in the simulation space. However, directly estimating the gravitational energy is very time consuming for 3D space with a large number of grids and is beyond our computation ability. So we refrain from directly deriving the magnetic critical parameter from energy ratios.

### 3.1.1. Properties of Magnetic Fields

Figure 1(a) shows the relation between the 3D total magnetic field strength within spheres ($B_{3D}^{\text{tot,sph}}$) and the volume density. It seems that models with smaller initial $\mu$ value (i.e., stronger initial magnetic field) tend to have stronger total magnetic field strength. The total magnetic field strength is positively correlated with the volume density and negatively correlated with the radius.

Figure 1(b) shows the relation between the $B_{3D}^{\text{tot,sph}}$ and the magnetic critical parameter within spheres ($\lambda_{\text{mod}}^{\text{sph}}$). Models with smaller initial $\mu$ values have smaller $\lambda_{\text{mod}}^{\text{sph}}$, suggesting that the magnetic critical parameter of a simulation is determined by the initial magnetic critical parameter. In our simulations, the magnetic critical parameter does not show a strong relation with radius. We also find that if the sink mass is not included in the calculation of the magnetic critical parameter, there will be a fake relation that the magnetic critical parameter is positively correlated with radius for each simulation.

### 3.1.2. Properties of Velocity Fields

The measured velocity dispersion contains the contributions from outflows, gravitational collapse (i.e., infall), rotation, and turbulence. To properly estimate the turbulent velocity dispersion, the contribution from other mechanisms must be removed from the total velocity dispersion.

We first consider the contribution from outflows. Because the resolution of our simulations (12.5 au or 25 au) does not allow the detailed process of disk accretion to be modeled, our simulations do not produce the commonly observed

parsec-scale jet-like or wide-angle outflows (Arce et al. 2007) that are expected to be driven by accretion-related protostellar winds (e.g., X-wind or disk winds: Pudritz et al. 2006; Shang et al. 2006). However, there are still some materials ejected from the center that may be related to some MHD processes in our simulations. These “outflows” would still contribute to the total velocity dispersion. Because we do not track the motion of particles, we cannot discern low-velocity outflows from other low-velocity gas motions in our simulations. Observationally, low-velocity outflowing gas usually cannot be distinguished from the environmental gas either. The distinguishable outflows usually have 1D outflowing velocities of >5 km s⁻¹ in massive star-forming regions (e.g., Beuther et al. 2002; Zhang et al. 2005; Liu et al. 2018). So we only investigate the contribution from high-velocity outflows in the simulations. To remove the contribution from the high-velocity outflows, we exclude grids with 3D velocity >10 km s⁻¹ (approximately with 1D velocity >5.8 km s⁻¹) and measure the mass-weighted velocity dispersion of the remaining grids as the nonoutflow velocity dispersion.

Figure 2(a) shows the ratio between the 3D nonoutflow velocity dispersion ($\delta v_{3D,\text{noout}}$) and the 3D total velocity dispersion ($\delta v_{3D}$) within different radii in all simulation models. In most of the cases, the contribution from outflows to the total velocity dispersion is negligible. Only in some cases, mostly at small radii, could the ratio between the nonoutflow and total velocity dispersions be considerably lower than 1.

We then consider the contribution from infall. We adopt a simple spherical symmetric collapsing model for infall. The mass-weighted inward radial velocity with respect to the most massive sink is measured as the average radial infall velocity ($v_{\text{r,inf}}$) within spheres of different radii. We ignore the infall toward sinks of smaller mass. During the calculation, we also exclude grids with 3D velocity >10 km s⁻¹. The 3D noninfall velocity dispersion is estimated as $\delta v_{3D,\text{noinf}} = (\delta v_{3D}^2 - v_{\text{r,inf}}^2)^{0.5}$.

Figure 2(b) shows the ratio between the 3D noninfall velocity dispersion and the 3D total velocity dispersion within different radii in all simulation models. In most cases, the ratio between the noninfall and total velocity dispersions are greater than 0.71 (corresponding to an energy ratio of 0.5), which suggests that the infall energy is mostly smaller than half of the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Upper: the ratio between the nonoutflow velocity dispersion and the total velocity dispersion. Middle: the ratio between the noninfall velocity dispersion and the total velocity dispersion. Bottom: the ratio between the nonrotational velocity dispersion and the total velocity dispersion. The velocity dispersion is the mass-weighted value. Different symbols correspond to different initial $\mu$ values, and different colors represent different radii.</figcaption></figure>

total kinetic energy. The average ratio between the two velocity dispersions is ~0.87, which indicates that the infall does not affect the velocity dispersion too much.

Finally, we consider the contribution from the rotation. We adopt a simple rotation model where the rotation is along the

axis of the total angular momentum. The total angular momentum is estimated by $L = \sum \mathbf{r}^{\text{cell}} \times m^{\text{cell}} \mathbf{v}^{\text{cell}}$, where $\mathbf{r}^{\text{cell}}$ is the position vector from the most massive sink to the position of an individual cell. The tangential velocity at the plane perpendicular to the axis of the total angular momentum is measured as the rotational velocity. We estimated the total angular momentum and mass-weighted rotational velocity within different radii for all simulation models. Similarly, we exclude grids with 3D velocity >10 km s⁻¹ in the calculations. Once the average rotational velocity ($v_{\text{t,rot}}$) is determined, the 3D nonrotational velocity dispersion can be estimated as $\delta v_{3D,norot} = (\delta v_{3D}^2 - v_{\text{t,rot}}^2)^{0.5}$.

Figure 2(c) shows the ratio between the 3D nonrotational velocity dispersion and the 3D total velocity dispersion within different radii in all simulation models. The average ratio between the two velocity dispersions is ~0.71 (corresponding to an energy ratio of ~0.5), which indicates the rotation could significantly affect the measured velocity dispersion and contribute to the total kinetic energy. However, this may be because our simulations do not produce outflows that are comparable to the observed ones. So the excess angular momentum cannot be transferred out via outflows, and the rotation in our simulations would be stronger than in the actual star-forming regions.

As mentioned above, the infall and rotation in our simulations can both contribute to the total velocity dispersion. However, it is hard to include both mechanisms in a simple model to exclude their combined contributions. Because the contribution from the rotation is larger than that of the infall in the total velocity dispersion, we adopt the nonrotational velocity dispersion as the turbulent velocity dispersion $\delta v_{\text{turb}}$ in the following analysis.

Figure 3(a) shows the relation between the 3D turbulent velocity dispersion within spheres ($\delta v_{3D,turb}^{\text{sph}}$) and the volume density. Models with different initial Mach numbers do not show a distinct difference in ranges of turbulent velocity dispersion, suggesting that the initial turbulence environment might not significantly affect the magnitude of turbulence in star-forming regions if there is no continuous turbulence injection and/or when the self-gravity is dominant. In this situation, the turbulence is from the star formation activities. Generally, there is no strong relation between the turbulent velocity dispersion and the density/radius.

Figure 3(b) shows the relation between the $\delta v_{3D,turb}^{\text{sph}}$ and the magnetic critical parameter. With the same initial $\mu$ value, models with larger initial Mach numbers do not show larger velocity dispersion. This confirms the conclusion that the initial turbulence could decay quickly if there is no continuous turbulence injection.

## 3.2. The Principle of the DCF Method

The DCF method is widely used to estimate the plane-of-sky magnetic field strength from polarized dust continuum observations. The basic assumption of the DCF method is that there is a prominent mean (uniform/ordered) magnetic field component $B^u = \langle B^{\text{loc}} \rangle$, where $B^{\text{loc}}$ is the local plane-of-sky magnetic field vector. Then, it assumes Alfvénic magnetic field perturbations, i.e., the perturbation on the mean magnetic field is due to turbulence, or in other words, there is equipartition between the turbulent magnetic energy<fnref n="9" /> $E_B^\text{t} = (B^u)^2 V / (2\mu_0)$

<footnote n="9">⁹ All of the magnetic-related equations and terms in this work are in SI units or cgs units unless otherwise noted.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left: the relation between the 3D turbulent velocity dispersion and the volume density. Different symbols correspond to different initial Mach numbers and different colors represent different radii. Right: the relation between the 3D turbulent velocity dispersion and the magnetic critical parameter. Different symbols represent different initial $\mu$ values and different colors represent different initial Mach numbers. The estimations are within different radii in all simulation models at two different time snapshots. The velocity dispersion is the mass-weighted value. Larger initial $M$ values correspond to stronger initial turbulence.</figcaption></figure>

<figure><img src="image_5.png" /></figure>

and the turbulent kinetic energy $E_K^t = \rho V \delta v^2 / 2$, where $B^t$ is the equivalent turbulent magnetic field strength, $V$ is the volume, and $\rho$ is the gas density. The equivalent turbulent magnetic field strength $B^t$ is estimated as $B^t = ((B^{\text{tot}})^2 - (B^u)^2)^{1/2}$, where the equivalent total (i.e., rms) magnetic field strength $B^{\text{tot}}$ is estimated as $B^{\text{tot}} = \langle(B^{\text{loc}})^2\rangle^{1/2}$. Thus, at the orientation perpendicular to ($\perp$) the mean plane-of-sky field and in the plane of the sky, we have $B_{\text{pos}\perp}^t = \sqrt{\mu_0 \rho} \delta v_{\text{pos}\perp}$, where “pos” stands for the plane of the sky. If the turbulence is isotropic (see Appendix D for a discussion), i.e., $\delta v_{\text{pos}\perp}^2 = \delta v_{\text{pos}\parallel}^2 = \delta v_{\text{los}}^2 = \delta v_{3D}^2 / 3$, where “‖” stands for the orientation parallel to the mean plane-of-sky magnetic field, “los” stands for the line of sight, and “3D” stands for three dimensions (3D), we have

$$B_{\text{pos}\perp}^{t} = \sqrt{\mu_{0}\rho} \delta v_{\text{los}}. \qquad (2)$$

To estimate $B_{\text{pos}}^u$, Equation (2) can be transformed to

$$B_{\text{pos}}^{u} = \sqrt{\mu_{0}\rho} \frac{\delta v_{\text{los}}}{B_{\text{pos}\perp}^{t} / B_{\text{pos}}^{u}}. \qquad (3)$$

If we define angle $\phi$ as the difference between the local polarization angle $\theta$ and the mean polarization angle $\bar{\theta}$, the relation between the local $\phi$ and magnetic field components can be written as

$$\tan \phi = \frac{|\mathbf{B}_{\text{pos}\perp}^{\text{t,loc}}|}{|\mathbf{B}_{\text{pos}\parallel}^{\text{t,loc}} + \mathbf{B}_{\text{pos}}^{u}|}, \qquad (4)$$

where $B_{\text{pos}\perp}^{\text{t,loc}}$ and $B_{\text{pos}\parallel}^{\text{t,loc}}$ are the local turbulent field vector components with orientations perpendicular to and parallel to $B_{\text{pos}}^u$, respectively. When the turbulent magnetic field is small compared to the ordered magnetic field ($B_{\text{pos}\parallel}^t \sim B_{\text{pos}\perp}^t \ll B_{\text{pos}}^u$), there is $B_{\text{pos}\perp}^t / B_{\text{pos}}^u \sim \delta(\tan \phi)$, where $\delta(\tan \phi)$ is the dispersion of $\tan \phi$. Under these assumptions, the uniform plane-of-sky magnetic field strength is given by

$$B_{\text{pos}}^{\text{u,dcf}} \sim \sqrt{\mu_{0}\rho} \frac{\delta v_{\text{los}}}{\delta(\tan \phi)}. \qquad (5)$$

When $\phi$ is small, there is $B_{\text{pos}\perp}^t / B_{\text{pos}}^u \sim \delta\phi$. The uniform plane-of-sky magnetic field strength in small-angle approximation is

given by

$$B_{\text{pos}}^{\text{u,dcf,sa}} \sim \sqrt{\mu_{0}\rho} \frac{\delta v_{\text{los}}}{\delta\phi}. \qquad (6)$$

On the other hand, Equation (2) can be rewritten as

$$B_{\text{pos}}^{\text{tot}} = \sqrt{\mu_{0}\rho} \frac{\delta v_{\text{los}}}{B_{\text{pos}\perp}^{t} / B_{\text{pos}}^{\text{tot}}}. \qquad (7)$$

to estimate the equivalent plane-of-sky total magnetic field strength $B_{\text{pos}}^{\text{tot}}$. Similarly, we have

$$\sin \phi = \frac{|\mathbf{B}_{\text{pos}\perp}^{\text{t,loc}}|}{|\mathbf{B}_{\text{pos}}^{\text{loc}}|}. \qquad (8)$$

If $B_{\text{pos}}^{\text{loc}}$ does not vary too much, there is $B_{\text{pos}\perp}^t / B_{\text{pos}}^{\text{tot}} \sim \delta(\sin \phi)$. Thus, the $B_{\text{pos}}^{\text{tot,dcf}}$ could be estimated as

$$B_{\text{pos}}^{\text{tot,dcf}} \sim \sqrt{\mu_{0}\rho} \frac{\delta v_{\text{los}}}{\delta(\sin \phi)}. \qquad (9)$$

Under small-angle approximation, $B_{\text{pos}\perp}^t / B_{\text{pos}}^{\text{tot}} \sim \delta\phi$ and the plane-of-sky total magnetic field strength is given by

$$B_{\text{pos}}^{\text{tot,dcf,sa}} \sim \sqrt{\mu_{0}\rho} \frac{\delta v_{\text{los}}}{\delta\phi}. \qquad (10)$$

### 3.3. Angular Dispersions and Ratios between Magnetic Field Components

#### 3.3.1. Comparing Angular Dispersions with Ratios between Magnetic Field Components in the Model

The DCF method assumes the turbulent-to-ordered magnetic field strength ratio $B_{\text{pos}\perp}^t / B_{\text{pos}}^u$ and the turbulent-to-total magnetic field strength ratio $B_{\text{pos}\perp}^t / B_{\text{pos}}^{\text{tot}}$ are traced by the dispersion of magnetic field position angles. Here we test the reliability of this assumption with our simulation results.

We measure the angular dispersions and the strength ratios between magnetic field components within spheres of different radii with respect to the most massive sink and at three orthogonal projections (xy, xz, and zy). The orientation of
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Dispersions of (a) $\phi_{\text{sim}}^{\text{sph}}$ and (b) $\tan \phi_{\text{sim}}^{\text{sph}}$ vs. turbulent-to-ordered magnetic field strength ratio for all simulation models within different radii with respect to the most massive sink and at three orthogonal projections. Different symbols represent different projections. Different colors correspond to different initial $\mu$ values. Solid lines correspond to 1:1 relation. Dashed lines correspond to 1:2 and 2:1 relations.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

plane-of-sky magnetic fields projected in the xy, xz, and zy planes in the simulation grids is gridded with a cell size of 0.0013 pc (see Appendix B for a discussion about grid size). In the calculation of the angular dispersion within a sphere, the direction of $B_{\text{pos}}^{\text{u}}$ within the space concerned is adopted as the orientation of the mean plane-of-sky magnetic field position angle. We define the angle $\phi_{\text{sim}}$ to represent the difference between the local and mean plane-of-sky magnetic field orientation in the simulation space. Then, we measure the dispersion of $\phi_{\text{sim}}$, $\sin \phi_{\text{sim}}$, and $\tan \phi_{\text{sim}}$. Because the magnetic energy is proportional to the volume, we adopt the volume-weighted value of the magnetic field strength. $B_{\text{pos}}^{\text{u}}$ is estimated as $B_{\text{pos}}^{\text{u}} = |\langle B_{\text{pos}}^{\text{cell}} \rangle|$, where $B_{\text{pos}}^{\text{cell}}$ is the plane-of-sky magnetic field vector of a cell. At orientation in the plane of the sky and perpendicular to the mean-field orientation, we have $B_{\text{pos}}^{\text{u}} = 0$ and $B_{\text{pos}}^{\text{tot}} = B_{\text{pos}}^{\text{t}} = (\langle B_{\text{pos}}^{\text{cell}} \rangle^2)^{1/2}$. It should be noted that because $B_{\text{pos}}^{\text{cell}}$ is the rms plane-of-sky magnetic field strength of a cell, $\langle B_{\text{pos}}^{\text{cell}} \rangle$ should not be used to represent the $B_{\text{pos}}^{\text{u}}$ as adopted by some previous simulation works.

Figure 4 shows the angular dispersions estimated within spheres of different radii in the simulation space ($\delta\phi_{\text{sim}}^{\text{sph}}$ and $\delta(\tan \phi_{\text{sim}}^{\text{sph}})$) versus the turbulent-to-ordered magnetic field strength ratio ($B_{\text{pos}}^{\text{t}}/B_{\text{pos}}^{\text{u}}$) measured in the same space. The angular dispersions in the zy plane, where the initial magnetic field orientation is along the line of sight (i.e., the x-axis), are generally greater than those in the xy and xz plane, where the initial magnetic field orientation is in the plane of the sky. This trend is consistent with what was found in Falceta-Gonçalves et al. (2008). Figure 4 also indicates that models with stronger initial magnetic fields (i.e., smaller $\mu$) tend to have smaller turbulent-to-ordered magnetic field strength ratios and $\delta\phi_{\text{sim}}^{\text{sph}}$ are better correlated with ($B_{\text{pos}}^{\text{t}}/B_{\text{pos}}^{\text{u}}$) for models with stronger initial magnetic fields. On the other hand, models with different initial $\mu$ values do not show a significant difference in the value range of the angular dispersion. Figure 4(a) shows that the ratio between $\delta\phi_{\text{sim}}^{\text{sph}}$ and ($B_{\text{pos}}^{\text{t}}/B_{\text{pos}}^{\text{u}}$) is from 1 to a factor of a few when $\delta\phi_{\text{sim}}^{\text{sph}} \lesssim 25^\circ$ (~0.44), while $\delta\phi_{\text{sim}}^{\text{sph}}$ can underestimate ($B_{\text{pos}}^{\text{t}}/B_{\text{pos}}^{\text{u}}$) by

up to more than two orders of magnitude for large $\delta\phi_{\text{sim}}^{\text{sph}}$ values. Figure 4(b) shows that $\delta(\tan \phi_{\text{sim}}^{\text{sph}})$ are correlated with ($B_{\text{pos}}^{\text{t}}/B_{\text{pos}}^{\text{u}}$) for $\delta(\tan \phi_{\text{sim}}^{\text{sph}}) \lesssim 1$. For large $\delta(\tan \phi_{\text{sim}}^{\text{sph}})$ values, there is significant scatter and $\delta(\tan \phi_{\text{sim}}^{\text{sph}})$ can greatly overestimate ($B_{\text{pos}}^{\text{t}}/B_{\text{pos}}^{\text{u}}$) for $\delta(\tan \phi_{\text{sim}}^{\text{sph}})$ by several orders of magnitude. In summary, $\delta\phi_{\text{sim}}^{\text{sph}}$ and $\delta(\tan \phi_{\text{sim}}^{\text{sph}})$ are correlated with the turbulent-to-ordered magnetic field strength ratio only when the angular dispersions are small.

Figure 5 shows the angular dispersions estimated within spheres in the simulation space ($\delta\phi_{\text{sim}}^{\text{sph}}$ and $\delta(\sin \phi_{\text{sim}}^{\text{sph}})$) versus the turbulent-to-total magnetic field strength ratio ($B_{\text{pos}}^{\text{t}}/B_{\text{pos}}^{\text{tot}}$) measured in the same space. Both $\delta\phi_{\text{sim}}^{\text{sph}}$ and $\delta(\sin \phi_{\text{sim}}^{\text{sph}})$ correlate well with ($B_{\text{pos}}^{\text{t}}/B_{\text{pos}}^{\text{tot}}$) for $\delta(\sin \phi_{\text{sim}}^{\text{sph}})$ by a factor of 2. Models with different initial $\mu$ values do not show a significant difference in the value range of the turbulent-to-total magnetic field strength ratio or the angular dispersion. Models with different initial $\mu$ values also do not show a significant difference in the correlation between the angular dispersion and the turbulent-to-total magnetic field strength ratio.

### 3.3.2. Comparing Angular Dispersions in the Model with Angular Dispersions in Synthetic Polarization Maps

The DCF method assumes the plane-of-sky angular dispersions of magnetic fields in 3D spaces are represented by the angular dispersions obtained from 2D polarization maps. As first proposed by Myers & Goodman (1991) and Zweibel (1990), the effect of signal integration along the line of sight may underestimate the angular dispersion because of a complex ordered magnetic field structure or perturbation from multiple independent turbulent cells. Here we investigate the signal integration effect on the measured angular dispersions.

We measure the observed angular dispersions, $\delta\phi_{\text{obs}}$, $\delta(\tan \phi_{\text{obs}})$, and $\delta(\sin \phi_{\text{obs}})$, within circles with different radii with respect to the projected position of the most massive sink on the synthetic polarization maps at three orthogonal projections. There are more than 200 uniformly sampled pixels within the smallest analyzing radius (~0.01 pc) in the polarization map, so the estimated angular dispersion does
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Dispersions of (a) $\phi_{\text{sim}}^{\text{sph}}$ and (b) $\sin\phi_{\text{sim}}^{\text{sph}}$ vs. turbulent-to-total magnetic field strength ratio for all simulation models within different radii with respect to the most massive sink and at three orthogonal projections. Different symbols represent different projections. Different colors correspond to different initial $\mu$ values. Solid lines correspond to 1:1 relation. Dashed lines correspond to 1:2 and 2:1 relations.</figcaption></figure>

<figure><img src="image_9.png" /></figure>

not suffer from the uncertainty from sparse sampling. The way to estimate the observed angular dispersion from 2D synthetic polarization maps is similar to that in Padoan et al. (2001), which is to shift the origin of the observed position angle $\phi_{\text{obs}}$ to the approximate minimum of the angular distribution and measure the standard deviations of $\phi_{\text{obs}}$, $\tan\phi_{\text{obs}}$, or $\sin\phi_{\text{obs}}$. In principle, the observed polarization emissions arise from all the polarized dust grains in the volume along the line of sight if the emission is optically thin. In practice, the polarized gas in the outer low-density space only has negligible contribution to the synthetic polarized dust emission due to mass concentration, so the observed polarized emission is usually correlated with condensed structures (i.e., clumps and cores). It is unclear whether the synthetic polarized emission is dominated by the contribution from the central condensed structures or the material at outer low-density space also has nonnegligible contribution, so we compare the observed angular dispersions within circles in polarization maps with the angular dispersions within both cylinders (indicated by the superscript “cyl”) and spheres in the simulation space. All of the space and area concerned are centered at the position or projected position of the most massive sink.

Figure 6 shows the relation between the angular dispersions estimated within spheres and cylinders in the simulation space and the observed angular dispersion within the projected circles in polarization maps. As shown in Figures 6(c) and (d), the angular dispersion in the form of $\delta(\tan\phi)$ has large scatter. Thus, we avoid using $\delta(\tan\phi)$ in the following analysis. As shown in Figures 6(a), (b), (e), and (f), the observed angular dispersions and angular dispersions in the simulation space are positively correlated with each other, and the observed angular dispersions are overall smaller than the angular dispersions in the simulation space. This underestimation of the angular dispersion could be due to line-of-sight signal integrations through more than one turbulent cell or due to field tangling of complex field geometry along the line of sight. For some data points within small radii (radii $\lesssim 0.1$ pc), we observe a significant underestimation of the angular dispersion in the model 3D space in some of our synthetic observations. Generally, the ratio between the angular dispersion in

polarization maps and the angular dispersion in simulation space at radii $>0.1$ pc is $\gtrsim 0.5$. There is no significant difference in the correlation between the angular dispersions in the simulation space and in the polarization map for models with different initial $\mu$ values.

Table 3 summarizes the average ratios between the angular dispersions in polarization maps and the angular dispersion in the simulation space for data points with radii $>0.1$ pc. These average ratios statistically take into account the contribution from the ordered field structure and should be regarded as the correction factor for the effect of signal integration along the line of sight.

### 3.3.3. Comparing Angular Dispersions in Polarization Maps with Ratios between Magnetic Field Components in the Model

We combine the results in Sections 3.3.1 and 3.3.2 and compare the angular dispersion in polarization maps with the turbulent-to-ordered or -total magnetic field strength ratio in simulation space (see Figure 7). In general, the measured angular dispersions in polarization maps are positively correlated with the turbulent-to-ordered or -total magnetic field strength ratio. Figures 7(a) and (b) indicate that $\delta\phi_{\text{obs}}$ would overall underestimate the turbulent-to-ordered magnetic field strength ratio. Data points with large $\delta\phi_{\text{obs}}$ (e.g., $>0.44, \sim 25^\circ$ or small radii (radii $\lesssim 0.1$ pc) show large scatter. Most of $\delta\phi_{\text{obs}}$ in the $zy$ plane are greater than $\sim 25^\circ$ and are poorly correlated with the turbulent-to-ordered magnetic field strength ratio. We also find that $\delta\phi_{\text{obs}}$ is better correlated with $B_{\text{pos}\perp}^{\text{t}}/B_{\text{pos}}^{\text{u}}$ for models with stronger initial magnetic fields. Figures 7(c)–(f) indicate that the measured angular dispersion would be significantly lower than the turbulent-to-total magnetic field strength ratio for some data points with a small radius (radii $\lesssim 0.1$ pc) and small angular dispersion. For data points with radii $>0.1$ pc, the measured angular dispersion is comparable to or only slightly lower than the turbulent-to-total magnetic field strength ratio. Overall, $\delta\phi_{\text{obs}}$ in the $zy$ plane are slightly greater than $B_{\text{pos}\perp}^{\text{t}}/B_{\text{pos}}^{\text{tot}}$, while $\delta\phi_{\text{obs}}$ in the $xz$ and $xy$ planes are slightly smaller than $B_{\text{pos}\perp}^{\text{t}}/B_{\text{pos}}^{\text{tot}}$. On the other hand, $\delta(\sin\phi_{\text{obs}})$ are generally smaller than the $B_{\text{pos}\perp}^{\text{t}}/B_{\text{pos}}^{\text{tot}}$ and
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Observed angular dispersion within circles vs. angular dispersions in the simulation space obtained in spheres (left panels) and cylinders (right panels) for all simulation models within different radii with respect to the most massive sink and at three orthogonal projections. Different symbols represent different initial $\mu$ values. Different colors correspond to different radii. Solid lines correspond to 1:1 relation. Dashed lines correspond to 1:2 and 2:1 relations.</figcaption></figure>
<caption>**Table 3**</caption>

Average Ratios between the Raw Angular Dispersions in Polarization Maps and Angular Dispersions in the Simulation Space for Data Points at Radii > 0.1 pc

<table><thead><tr><th>Average Ratios</th><th>Spheres<sup>a</sup></th><th>Cylinders<sup>a</sup></th></tr></thead><tbody><tr><td>δφ<sub>obs</sub> / δφ<sub>sim</sub></td><td>0.70(28%)</td><td>0.77(31%)</td></tr><tr><td>δ(sin φ<sub>obs</sub>) / δ(sin φ<sub>sim</sub>)</td><td>0.68(26%)</td><td>0.73(31%)</td></tr></tbody></table>

**Notes.** Values in parentheses are the relative uncertainty.

<sup>a</sup> Indicates whether the angular dispersion in the simulation space are estimated in spheres or cylinders.

δ(sin φ<sub>obs</sub>) in the zy plane are better correlated with $B_{\text{pos}\perp}^t / B_{\text{pos}}^{\text{tot}}$ than in the xz and xy planes. There is no significant difference in the correlation between the angular dispersion and the turbulent-to-total magnetic field strength ratio for models with different initial μ values.

Table 4 summarizes the average ratios between the angular dispersions in polarization maps and turbulent-to-ordered or -total magnetic field strength ratios for data points with radii > 0.1 pc. Due to large scatter, we exclude data points with δφ<sub>obs</sub> > 25° (~0.44) in the estimation of δφ<sub>obs</sub>/(B<sub>pos⊥</sub><sup>t</sup>/B<sub>pos</sub><sup>u</sup>). If other assumptions (i.e., isotropic turbulence and energy equipartition) of the DCF method are satisfied, the ratios reported in Table 4 should be regarded as the correction factors for the plane-of-sky uniform or total magnetic field strength derived from the DCF method. These correction factors statistically take into account the contribution from the ordered field structure to the measured angular dispersion at clump and core scales.

## 3.4. Uncertainties in Estimating the Angular Dispersion

### 3.4.1. The ADF Method

There are several effects that may affect the measured angular dispersion. The effect of measuring angular dispersion at scales smaller than the turbulent correlation scale (hereafter the turbulent correlation effect; Hildebrand et al. 2009), the beam-smoothing effect of telescopes, and the effect of line-of-sight signal integration can underestimate the angular dispersion and thus overestimate the magnetic field strength. On the other hand, the contribution from the underlying ordered (large-scale) field structure could overestimate the angular dispersion and thus underestimate the magnetic field strength. The spatial filtering effect of interferometers could also affect the measured angular dispersion and magnetic field strength.

There were several approaches to quantitatively account for these effects that may affect the measured angular dispersion. The ADF method (Falceta-Gonçalves et al. 2008; Hildebrand et al. 2009; Houde et al. 2009, 2016) is the only method to analytically and simultaneously account for the ordered field contribution, the beam-smoothing effect, the turbulent correlation effect, and the interferometric filtering effect. Here we investigate the accuracy of the ADF method with our simulation results.

In Appendix C, we tested the ADF only accounting for the ordered field structure with simple Monte Carlo simulations, and the ADF works well on taking into account this effect. Further, the ADF accounting for the ordered field structure and the plane-of-sky turbulent correlation effect is given by

Hildebrand et al. (2009) and Houde et al. (2009, 2016):

$$1 - \langle \cos[\Delta\Phi(l)] \rangle \approx a_2 l^2 + \left( \frac{\langle B_t^2 \rangle}{\langle B^2 \rangle} \right)_{\text{or,tc}}^{\text{adf}} \times [1 - e^{-l^2/2l_\delta^2}], \quad (11)$$

where ΔΦ(l) is the angular difference of two magnetic field line segments separated by a distance l, $a_2 l^2$ is the first term of the Taylor expansion of the ordered component of ADF, $l_\delta$ is the turbulent correlation length, $\langle(B_t^2)\rangle/\langle B^2\rangle_{\text{or,tc}}^{\text{adf}} = \langle(B_t^2)\rangle/\langle B_0^2\rangle_{\text{or,tc}}^{\text{adf}} / ((\langle B_t^2 \rangle/\langle B_0^2 \rangle)_{\text{or,tc}}^{\text{adf}} + 1)$ is the turbulent-to-total magnetic energy ratio, and $\langle(B_t^2)\rangle/\langle B_0^2\rangle_{\text{or,tc}}^{\text{adf}}$ is the turbulent-to-ordered energy ratio.

We derive the ADFs within different radii in the polarization maps for all simulation models. Calculation of the ADFs can be very time consuming if a large number of angle samples are concerned. Thus, we bin the synthetic Q and U maps with an appropriate bin size and regenerate the binned position angles for some of the calculations to make sure the number of angle samples is smaller than 3000. Using the first term of the Taylor expansion to represent the ordered magnetic field is only valid when the scale l is smaller than an appropriate displacement d, while the displacement d varies case by case. Thus, we fit each specific ADF over different maximum spatial scales and obtain the best fit for each scale via reduced χ² minimization. We compare the reduced χ² of the best fit for different spatial scales and adopt the one with the smallest reduced χ² as the best fit for each specific ADF.

Figure 8 shows the relation between the directly estimated angular dispersions (see Section 3.3.2) in synthetic polarization maps and the turbulent-to-ordered or -total magnetic field strength ratio derived from the ADF method. Generally, the directly estimated angular dispersions are greater than the turbulent-to-ordered or -total magnetic field strength ratio derived from the ADF method, which can be explained by the contribution from large-scale field structure. If the turbulent-to-ordered or -total magnetic field strength ratio correctly traces the true angular dispersion in polarization maps, Figure 8 indicates that the contribution from ordered field structure can overestimate the angular dispersion by a few factors when the intrinsic angular dispersion in polarization maps is small. The average δφ<sub>obs</sub>/(⟨B<sub>t</sub><sup>2</sup>⟩/⟨B<sub>0</sub><sup>2</sup>⟩)<sub>or,tc</sub><sup>adf</sup>, δφ<sub>obs</sub>/(⟨B<sub>t</sub><sup>2</sup>⟩/⟨B<sup>2</sup>⟩)<sub>or,tc</sub><sup>adf</sup>, and δ(sin φ<sub>obs</sub>)/(⟨B<sub>t</sub><sup>2</sup>⟩/⟨B<sup>2</sup>⟩)<sub>or,tc</sub><sup>adf</sup> are ~2.7, ~2.8, and ~2.1, respectively, suggesting that the contribution from the ordered field may overestimate the angular dispersion by a factor of 2–3 on average.

If the signal integration along the line of sight is considered, the turbulent-to-ordered energy ratio is given by $\langle(B_t^2)\rangle/\langle B_0^2\rangle_{\text{or,tc,si}}^{\text{adf}} = (\langle B_t^2 \rangle/\langle B_0^2 \rangle)_{\text{or,tc}}^{\text{adf}} \times N'_{\text{adf}}$, where $N'_{\text{adf}}$ is the number of turbulent cells along the line of sight. $N'_{\text{adf}}$ is given by

$$N'_{\text{adf}} = \frac{l_{\Delta}}{\sqrt{2\pi l_{\delta}}}, \quad (12)$$

where $l_\Delta$ is the effective thickness along the line of sight. Here an assumption is adopted that the line-of-sight turbulent correlation scale is identical to the plane-of-sky turbulent correlation scale. Similarly, the turbulent-to-total energy ratio $\langle(B_t^2)\rangle/\langle B^2\rangle_{\text{or,tc,si}}^{\text{adf}}$ can be derived from the turbulent-to-ordered energy ratio taking into account this signal integration effect.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 7.** Angular dispersion within circles in polarization maps vs. turbulent-to-ordered or -total magnetic field strength ratios in the simulation space obtained in spheres (left panels) and cylinders (right panels) for all simulation models within different radii with respect to the most massive sink and at three orthogonal projections. Different symbols represent different projections. Different colors correspond to different radii. Solid lines correspond to 1:1 relation. Dashed lines correspond to 1:2 and 2:1 relations.</figcaption></figure>
<caption>Table 4</caption>

Average Ratios between the Angular Dispersions in Polarization Maps and Turbulent-to-ordered or -total Magnetic Field Strength Ratios for Data Points at Radii > 0.1 pc

<table><thead><tr><th>Average Ratios</th><th>Spheres<sup>a</sup></th><th>Cylinders<sup>a</sup></th></tr></thead><tbody><tr><td>δφ<sub>obs</sub>/(B<sub>pos⊥</sub><sup>t</sup>/B<sub>pos</sub><sup>u</sup>)<sup>b</sup></td><td>0.25(46%)</td><td>0.25(49%)</td></tr><tr><td>δφ<sub>obs</sub>/(B<sub>pos⊥</sub><sup>t</sup>/B<sub>pos</sub><sup>tot</sup>)</td><td>0.81(43%)</td><td>0.82(41%)</td></tr><tr><td>δ(sin φ<sub>obs</sub>)/(B<sub>pos⊥</sub><sup>t</sup>/B<sub>pos</sub><sup>tot</sup>)</td><td>0.57(29%)</td><td>0.58(28%)</td></tr></tbody></table>

**Notes.** The ratios are averaged over three different planes (xy, xz, and zy). These ratios should be regarded as the correction factor for the estimated magnetic field strength if assumptions of energy equipartition and isotropic turbulence are satisfied. Values in parentheses are the relative uncertainty.

<sup>a</sup> Indicates whether the turbulent-to-ordered or -total magnetic field strength ratios are estimated in spheres or cylinders in the simulation space.

<sup>b</sup> Additional selection criterion of δφ<sub>obs</sub> ≲ 25° (~0.44) is adopted.

Following Houde et al. (2009), we adopt the width at half of the maximum of the normalized autocorrelation function of the integrated normalized polarized flux as the effective thickness $l_Δ$. However, for ~35% of our estimations, the minimum of the normalized autocorrelation function of the integrated normalized polarized flux is greater than half of the maximum, so the effective thickness cannot be derived in this way. Thus, these data are not included in the following ADF analysis.

The turbulent-to-ordered or -total magnetic field strength ratio derived from the ADF method should be compared to the turbulent-to-ordered or -total magnetic field strength ratio without the contribution from ordered field structure in the 3D simulation space. However, due to the lack of an appropriate method to account for the 3D field geometry, this cannot be done in our analysis. In principle, the ADF method can be applied to the 3D field to remove the 3D field geometry and derive the 3D turbulent-to-ordered or -total field strength ratio. In 3D space, there is no need to take into account the line-of-sight signal integration. However, the upper limit of the derivable turbulent-to-ordered or -total field strength ratio without accounting for the signal integration effect are ~0.76 and ~0.6, respectively, for the ADF method (see Appendix C). So the 3D turbulent-to-total field strength ratio could be underestimated, and the 3D turbulent-to-ordered field strength ratio could be greatly underestimated if there are significantly turbulent fields. Application of the ADF method on 3D space with large grids could also be very time consuming. On the other hand, the unsharp masking method (Pattle et al. 2017) may also be applied to remove the 3D field structure, but its correctness is still uncertain in 3D space because this method does not consider the turbulent correlation length. Thus, we refrain from deriving the 3D turbulent-to-ordered or -total field strength ratio without the contribution from the ordered field structure.

Alternatively, to investigate whether the ADF method can correctly account for the effect of signal integration along the line of sight, we correct the directly measured angular dispersion in polarization maps by a factor of $\sqrt{N'_{adf}}$ and compare the corrected angular dispersion with the angular dispersion in simulation grids (see Figure 9). For some data points, the angular dispersion after correction exceeds the value expected for a random field, suggesting that the number of turbulent cells along the line of sight is overestimated by the ADF method and the correction of $\sqrt{N'_{adf}}$ is inappropriate for

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 8.** Ratio between the directly estimated angular dispersions in synthetic polarization maps and the turbulent-to-ordered or -total magnetic field strength ratio derived by the ADF method. Different symbols represent different initial Mach numbers. Different colors correspond to different initial $\mu$ values.</figcaption></figure>

these data points. For other data points, the angular dispersion after correction generally agrees with the angular dispersions in simulation grids within a factor of ~2. The average ratios between the corrected angular dispersion and the angular dispersion in simulation space for these “good” data points are listed in Table 5. The average ratios are ~1, suggesting that the
<figure><img src="image_20.png" /><figcaption>Figure 9. Observed angular dispersion within circles in polarization maps corrected for $\sqrt{N'_{\text{adf}}}$ vs. angular dispersions in the simulation space obtained in spheres (left panels) and cylinders (right panels) for all simulation models within different radii with respect to the most massive sink and at three orthogonal projections. Different symbols represent different projections. Different colors correspond to different radii. Solid lines correspond to 1:1 relation. Dashed lines correspond to 1:2 and 2:1 relations. The horizontal dotted lines indicate the expected angular dispersion for random fields (~52° for $\delta\phi$ and ~0.58 for $\delta(\sin\phi)$).</figcaption></figure>

ADF method can correctly account for the line-of-sight signal integration effect. However, these good data points only make up 20%–30% of our total estimations, which suggest the ADF method might not be applicable to account for the line-of-sight signal integration in the majority of cases.

Once $(\langle B_t^2 \rangle / \langle B_0^2 \rangle)_{\text{or,tc,si}}^{\text{adf}}$ and $(\langle B_t^2 \rangle / \langle B^2 \rangle)_{\text{or,tc,si}}^{\text{adf}}$ are derived, the plane-of-sky uniform magnetic field strength can be derived by<fnref n="10" />

$$B_{\text{pos,or,tc,si}}^{\text{u,adf}} \sim \sqrt{\mu_0 \rho} \delta v_{\text{los}} \left( \left( \frac{\langle B_t^2 \rangle}{\langle B_0^2 \rangle} \right)_{\text{or,tc,si}}^{\text{adf}} \right)^{-\frac{1}{2}}, \quad (13)$$

and the plane-of-sky total magnetic field strength can be derived by

$$B_{\text{pos,or,tc,si}}^{\text{tot,adf}} \sim \sqrt{\mu_0 \rho} \delta v_{\text{los}} \left( \left( \frac{\langle B_t^2 \rangle}{\langle B^2 \rangle} \right)_{\text{or,tc,si}}^{\text{adf}} \right)^{-\frac{1}{2}}. \quad (14)$$

The plane-of-sky uniform magnetic field strength $B_{\text{pos,or,tc}}^{\text{u,adf}}$ and the plane-of-sky total magnetic field strength $B_{\text{pos,or,tc}}^{\text{tot,adf}}$ without accounting for the line-of-sight signal integration effect can be derived in similar ways.

### 3.4.2. The Cho & Yoo (2016) Method

As mentioned before, the amount of angular dispersion measured in polarization maps would be reduced by a factor of $\sqrt{N}$ compared to the intrinsic angular dispersion, where N is the number of independent turbulent cells along the line of sight (Hildebrand et al. 2009; Cho & Yoo 2016). The ADF method analytically solves for the plane-of-sky turbulent correlation scale from the distribution of plane-of-sky polarization angles

<footnote n="10">¹⁰ It should be noted that there is an approximation in Equation (26) in Houde et al. (2016), $B_0 \simeq \sqrt{4\pi\rho\sigma(v)(\langle B_t^2 \rangle / \langle B^2 \rangle)^{-1/2}}$ in Gaussian units, where the turbulent-to-total magnetic field strength ratio is used to derive the uniform magnetic field strength. The approximation is only valid when the turbulent field is relatively small compared to the uniform field. The form without approximation is $B_0 \simeq \sqrt{4\pi\rho\sigma(v)(\langle B_t^2 \rangle / \langle B_0^2 \rangle)^{-1/2}}$ or $B \simeq \sqrt{4\pi\rho\sigma(v)(\langle B_t^2 \rangle / \langle B^2 \rangle)^{-1/2}}$ in Gaussian units.</footnote>
and derives the number of line-of-sight turbulent cells with the assumption that the line-of-sight turbulent correlation scale is identical to the plane-of-sky turbulent correlation scale. We test the ADF method in Section 3.4.1 and find that the ADF method only correctly accounts for this effect for a small portion of the total estimations. Alternatively, Cho & Yoo (2016) proposed an approach (hereafter the CY16 method) to estimate *N* from the ratio between the average line-of-sight velocity dispersion (δv<sub>los</sub>) and the standard deviation of centroid velocities (δv<sub>c</sub>) in the considered region:

$$N_{\text{cy16}} = \left( \frac{\delta v_{\text{los}}}{\delta v_c} \right)^2. \qquad (15)$$

We apply the CY16 method on our simulation data at different radii at three different projections for all simulation models. We adopt the mass-weighted standard deviation of centroid velocities within cylinders as δv<sub>c</sub>. The δv<sub>los</sub> is also estimated within cylinders. We derive N<sub>cy16</sub> and use it to correct for the angular dispersion measured in polarization maps.

Figure 10 shows the relation between the corrected angular dispersions and the angular dispersions in the simulation space. The data points show large scatter at radii <0.1 pc. At radii >0.1 pc, the corrected angular dispersion and the angular dispersion in simulation grids are mostly consistent with each other within a factor of 2, and most corrected angular dispersions are below the value expected for random fields.

The average ratios between the angular dispersion corrected by the CY16 method and the angular dispersion in simulation space at radii >0.1 are listed in Table 5. We notice that the relative uncertainty of the average ratios corrected by the CY16 method is slightly lower than the relative uncertainty of the average ratios without correction (see Table 3), suggesting that the CY16 method provides a slightly better estimation than that without correction.

In Figure 11 we compare the number of line-of-sight turbulent cells derived from the ADF method and the CY16 method. For most data points at radii <0.1 pc, N<sub>cy16</sub> could be significantly larger than N'<sub>adf</sub>. On the other hand, the value of N<sub>cy16</sub> shows smaller variations than N'<sub>adf</sub> at radii >0.1 pc.

### 3.4.3. Smoothing and Filtering Effects of Synthetic Interferometric Observations

As revealed by some previous numerical studies (Heitsch et al. 2001; Falceta-Gonçalves et al. 2008), the effect of finite beam resolution of telescopes tends to underestimate the angular dispersion and thus overestimate the magnetic field strength. When it comes to interferometric observations, there arise additional issues about the spatial filtering of large-scale structures. Here we investigate how the combination of the beam-smoothing effect and the filtering effect of interferometers affects the measured angular dispersion. It should be noted that the analyses carried out in previous subsections (Sections 3.3, 3.4.1, and 3.4.2) are on perfect synthetic skies, while the analyses in this subsection are carried out in synthetic ALMA polarization observations. The synthetic ALMA polarization maps are gridded with pixel sizes of one-fifth of the beam sizes.

First, we directly measure the dispersion of magnetic field position angles in synthetic ALMA maps (δφ<sub>alma</sub>) in a similar way to the estimation of δφ<sub>obs</sub> in Section 3.3.2. Only position angles at locations where *I* > 0 are considered in the estimation.

<table><caption>**Table 5**</caption><caption>Average Ratios between the Corrected Angular Dispersions in Polarization Maps and Angular Dispersions in the Simulation Space for Data Points at Radii > 0.1 pc</caption><thead><tr><th>Average Ratios</th><th>Spheres<sup>a</sup></th><th>Cylinders<sup>a</sup></th></tr></thead><tbody><tr><td>√N'<sub>adf</sub> δφ<sub>obs</sub> / δφ<sub>sim</sub><sup>b</sup></td><td>1.00(35%)</td><td>1.06(37%)</td></tr><tr><td>√N'<sub>adf</sub> δ(sin φ<sub>obs</sub>) / δ(sin φ<sub>sim</sub>)<sup>c</sup></td><td>0.94(20%)</td><td>1.00(18%)</td></tr><tr><td>√N<sub>cy16</sub> δφ<sub>obs</sub> / δφ<sub>sim</sub><sup>b</sup></td><td>0.86(23%)</td><td>0.94(27%)</td></tr><tr><td>√N<sub>cy16</sub> δ(sin φ<sub>obs</sub>) / δ(sin φ<sub>sim</sub>)<sup>c</sup></td><td>0.87(20%)</td><td>0.94(23%)</td></tr></tbody></table>

**Notes.** Values in parentheses are the relative uncertainty.

<sup>a</sup> Indicates whether the angular dispersion in the simulation space is estimated in spheres or cylinders.

<sup>b</sup> Data points with corrected angular dispersion >52° are not included.

<sup>c</sup> Data points with corrected angular dispersion >0.58 are not included.

Figure 12 shows the relation between δφ<sub>alma</sub> and δφ<sub>obs</sub>. When the radius of the area concerned, *R*, is comparable to the synthesized beam size *l*<sub>beam</sub> (i.e., *R* ~ *l*<sub>beam</sub>), δφ<sub>alma</sub> is overall smaller than δφ<sub>obs</sub>, which means the synthetic ALMA observations would underestimate the angular dispersion when the area concerned is barely resolved. This could be easily explained by the loss of perturbation information because of beam smoothing. On the other hand, for a significant amount of data points, the δφ<sub>alma</sub> are larger than δφ<sub>obs</sub>. When θ<sub>obs</sub> is small, the synthetic ALMA observations could significantly overestimate the angular dispersion. As *R* approaches the interferometer filtering scale *l*<sub>filter</sub> (l<sub>filter</sub> ~ 8.5, 9.8, 11.6, 12.2, and 12.4 × *l*<sub>beam</sub> for AC1, AC2, AC3, AC4, and AC5, respectively), the overestimation of angular dispersion is more obvious. This overestimation might be because the filtering effect of interferometric observations has altered the large-scale magnetic field structure (see Figures 17 and 18 for examples), or because there is not enough samples within the area to achieve statistical significance due to the exclusion of some angle samples at locations where *I* < 0.

The ADF method is the only method that analytically takes into account the two observational effects. In the most complicated form of the ADF (Houde et al. 2016), it takes into account the contribution from ordered field structure, the turbulent correlation effect, the signal integration effect, the beam-smoothing effect, and the interferometer filtering effect. The ADF is given by Houde et al. (2016):

$$1 - \langle \cos[\Delta\Phi(l)] \rangle \approx a_2 l^2 + \left[ \frac{N_{\text{adf}}}{1 + \frac{N_{\text{adf}}}{\left(\frac{B_1^2}{B_0^2}\right)^{\text{adf}_{\text{or,tc,si,bs,if}}}}} \right] \\ \times \left\{ \frac{1}{N_1} \left[ 1 - e^{-l^2/2(l_1^2+2W_1^2)} \right] \right. \\ + \frac{1}{N_2} \left[ 1 - e^{-l^2/2(l_2^2+2W_2^2)} \right] \\ \left. - \frac{2}{N_{12}} \left[ 1 - e^{-l^2/2(l_2^2+W_1^2+W_2^2)} \right] \right\}, \quad (16)$$

where W<sub>1</sub> and W<sub>2</sub> are the standard deviations of the Gaussian profiles of the synthesized beam and the large-scale filtering (i.e., the FWHM divided by √8 ln 2), respectively. N<sub>adf</sub> is the number of turbulent cells along the line of sight probed by the
<figure><img src="image_21.png" /><figcaption>**Figure 10.** Observed angular dispersion within circles in polarization maps corrected for $\sqrt{N_{cy16}}$ vs. angular dispersions in the simulation space obtained in spheres (left panels) and cylinders (right panels) for all simulation models within different radii with respect to the most massive sink and at three orthogonal projections. Different symbols represent different projections. Different colors correspond to different radii. Solid lines correspond to 1:1 relation. Dashed lines correspond to 1:2 and 2:1 relations. The horizontal dotted lines indicate the expected angular dispersion for random fields (~52° for $\delta\phi$ and ~0.58 for $\delta(\sin\phi)$).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 11.** Comparison between $N_{cy16}$ and $N'_{adf}$. Different symbols represent different projections. Different colors correspond to different radii. Solid lines correspond to 1:1 relation. Dashed lines correspond to 1:2 and 2:1 relations.</figcaption></figure>

telescope beam given by

$$N_1 = \frac{(l_{\delta}^2 + 2W_1^2)l_{\Delta}}{\sqrt{2\pi}l_{\delta}^3}, \qquad (17)$$

$$N_2 = \frac{(l_{\delta}^2 + 2W_2^2)l_{\Delta}}{\sqrt{2\pi}l_{\delta}^3}, \qquad (18)$$

$$N_{12} = \frac{(l_{\delta}^2 + W_1^2 + W_2^2)l_{\Delta}}{\sqrt{2\pi}l_{\delta}^3}, \qquad (19)$$

$$N_{\text{adf}} = \left( \frac{1}{N_1} + \frac{1}{N_2} - \frac{2}{N_{12}} \right)^{-1}. \qquad (20)$$

Similar to the approach in Section 3.4.1, we derive the ADFs within different radii in the synthetic ALMA maps, fit the ADFs, and obtain the best fit via reduced $\chi^2$ minimization. When $R \sim l_{\text{beam}}$ (i.e., the area concerned is barely resolved), the fitting of the ADFs does not converge. Once the turbulent-to-ordered magnetic field energy ratio taking into consideration all
<figure><img src="image_23.png" /><figcaption>**Figure 12.** Comparison between the directly measured angular dispersion in raw synthetic polarization maps, $\delta\phi_{\text{obs}}$, and the directly measured angular dispersion in synthetic ALMA polarization maps, $\delta\phi_{\text{alma}}$. Different colors correspond to different radii in units of multiples of the synthesized beam size. Solid lines correspond to the 1:1 relation. Dashed lines correspond to the 1:2 and 2:1 relations. The beam size and the filtering scale of each configuration of ALMA are indicated in the title.</figcaption></figure>

the effects, $\langle(B_t^2)\rangle/\langle B_0^2\rangle_{\text{or,tc,si,bs,if}}^{\text{adf}}$, is derived, the turbulent-to-ordered magnetic field energy ratio without accounting for the signal integration along the line of sight can be derived from $\langle(B_t^2)\rangle/\langle B_0^2\rangle_{\text{or,tc,bs,if}}^{\text{adf}} = \langle(B_t^2)\rangle/\langle B_0^2\rangle_{\text{or,tc,si,bs,if}}^{\text{adf}}/N_{\text{adf}}$.

Figures 13(a)–(c) compare the $(\langle B_t^2 \rangle / \langle B_0^2 \rangle)_{\text{or,tc}}^{\text{adf}}$ derived by fitting ADFs of raw (unsmoothed) synthetic polarization maps (see Section 3.4.1) with the $(\langle B_t^2 \rangle / \langle B_0^2 \rangle)_{\text{or,tc,bs,if}}^{\text{adf}}$ derived by fitting ADFs of synthetic ALMA polarization maps. Generally the $(\langle B_t^2 \rangle / \langle B_0^2 \rangle)_{\text{or,tc}}^{\text{adf}}$ and $(\langle B_t^2 \rangle / \langle B_0^2 \rangle)_{\text{or,tc,bs,if}}^{\text{adf}}$ agree with each other within a factor of 2, which indicates that the ADF method does correctly take into account the beam-smoothing effect and the interferometric filtering effect.

Figures 13(d)–(f) compare the $(\langle B_t^2 \rangle / \langle B_0^2 \rangle)_{\text{or,tc,si}}^{\text{adf}}$ derived by fitting ADFs of raw synthetic polarization maps with the $(\langle B_t^2 \rangle / \langle B_0^2 \rangle)_{\text{or,tc,si,bs,if}}^{\text{adf}}$ derived by fitting ADFs of synthetic ALMA polarization maps. Generally the relation between $(\langle B_t^2 \rangle / \langle B_0^2 \rangle)_{\text{or,tc,si}}^{\text{adf}}$ and $(\langle B_t^2 \rangle / \langle B_0^2 \rangle)_{\text{or,tc,si,bs,if}}^{\text{adf}}$ shows large scatter, which indicates that the ADF method may not simultaneously correctly account for the line-of-sight signal integration, beam smoothing, and interferometric filtering effects. This may be because the beam resolution and the pixel

size are insufficient to trace the turbulent correlation scale (see Appendix B), the turbulent correlation scale is anisotropic (Goldreich & Sridhar 1995), the turbulent correlation scales at different depth are not uniform, the information of turbulent correlation scale is lost due to observational uncertainty, or because of some other reasons.

### 3.5. Equipartition between the Turbulent Kinetic Energy and the Turbulent Magnetic Energy

The DCF method assumes the equipartition between the turbulent kinetic energy $E_{K,3D}^t$ and the turbulent magnetic energy $E_{B,3D}^t$. Falceta-Gonçalves et al. (2008) proposed that the energy equipartition assumption is fulfilled in sub-Alfvénic (i.e., the turbulent kinetic energy $E_{K,3D}^t$ is smaller than the uniform magnetic field energy $E_{B,3D}^u$) cases and found that a super-Alfvénic (i.e., $E_{K,3D}^t > E_{B,3D}^u$) model with initial $E_{B,3D}^t < E_{K,3D}^t$ can increase its $E_{B,3D}^t$ as time evolves and reach a state close to the energy equipartition ($E_{B,3D}^t \sim E_{K,3D}^t$) after $0.3 \times$ the crossing timescale.
<figure><img src="image_24.png" /></figure>

a

<figure><img src="image_25.png" /></figure>

c

<figure><img src="image_26.png" /></figure>

e

<figure><img src="image_27.png" /></figure>

b

<figure><img src="image_28.png" /></figure>

d

<figure><img src="image_29.png" /></figure>

f

**Figure 13.** (a)–(c) Comparison between $((B_t^2)/\langle B_0^2 \rangle)_{\text{or,tc}}^{\text{adf}}$ and $((B_t^2)/\langle B_0^2 \rangle)_{\text{or,tc,bs,if}}^{\text{adf}}$. (d)–(f) Comparison between $((B_t^2)/\langle B_0^2 \rangle)_{\text{or,tc,si}}^{\text{adf}}$ and $((B_t^2)/\langle B_0^2 \rangle)_{\text{or,tc,si,bs,if}}^{\text{adf}}$. Different colors correspond to different radii in units of multiples of the synthesized beam size. Solid lines correspond to the 1:1 relation. Dashed lines correspond to the 1:2 and 2:1 relations.

We measure the energies in the 3D volume, $E_{\text{B},3\text{D}}^{\text{t}}$, $E_{\text{B},3\text{D}}^{\text{u}}$, $E_{\text{K},3\text{D}}^{\text{t}}$; the 3D turbulent velocity dispersion $\delta v_{3\text{D,turb}}$; and the 3D turbulent-to-total magnetic field strength ratio $B_{3\text{D}}^{\text{t}}/B_{3\text{D}}^{\text{tot}}$ within

spheres of different radii for all models. It should be noted that in our estimation of the turbulent magnetic field component, the contribution from large-scale field structure is not excluded, so
$E_{B,3D}^t$ and $B_{3D}^t$ could be overestimated and $E_{B,3D}^u$ could be underestimated.

Figure 14(a) shows the relation between the ratios $E_{B,3D}^t/E_{K,3D}^t$ and $B_{3D}^t/B_{3D}^{tot}$. The turbulent magnetic energy is overall less than the turbulent kinetic energy. For large $B_{3D}^t/B_{3D}^{tot}$ values, the ratio $E_{B,3D}^t/E_{K,3D}^t$ shows large scatter over several orders of magnitude. When $B_{3D}^t/B_{3D}^{tot} \lesssim 0.8$ ($B_{3D}^t/B_{3D}^u \lesssim 1.3$), the ratio $E_{B,3D}^t/E_{K,3D}^t$ is between ~0.3 and ~1.6, and the average $E_{B,3D}^t/E_{K,3D}^t$ is ~0.9, suggesting that the two energies tend to be in equipartition for small turbulent-total field strength ratios. Figure 14(a) also shows that models with stronger initial magnetic field (i.e., smaller initial $\mu$ value) tend to have smaller $B_{3D}^t/B_{3D}^{tot}$ and larger $E_{B,3D}^t/E_{K,3D}^t$ values. However, such relation is not found between the plane-of-sky ratio $B_{pos}^t/B_{pos}^{tot}$ and the initial $\mu$ value (see Figure 5). This might be a consequence of the uncertainty from the projection effect when measuring quantities on the plane of sky (see Appendix D).

Figure 14(b) shows the relation between the ratio $E_{B,3D}^t/E_{K,3D}^t$ and the $\delta v_{3D,turb}$. The same as mentioned above, models with a stronger magnetic field tend to be closer to an energy equipartition. There is no strong relation between the energy ratio and the turbulent velocity dispersion in general. However, for strong-field models ($\mu=1.2$ or $\mu=2$), data points with smaller turbulent velocity dispersion, mostly at large radii ($>0.1$ pc), seem to be closer to the energy equipartition. For strong-field models at radii $>0.1$ pc, the $E_{B,3D}^t/E_{K,3D}^t$ ranges from ~0.4 to ~1.3 and the average $E_{B,3D}^t/E_{K,3D}^t$ is ~0.9, which is quite close to an energy equipartition.

Figure 14(c) shows that the energy ratio $E_{B,3D}^t/E_{K,3D}^t$ is positively correlated with $E_{B,3D}^u/E_{K,3D}^t$. Models with stronger initial magnetic field tend to have larger $E_{B,3D}^t/E_{K,3D}^t$ and $E_{B,3D}^u/E_{K,3D}^t$ values. As $E_{B,3D}^u/E_{K,3D}^t$ increases, $E_{B,3D}^t/E_{K,3D}^t$ reaches unity when $E_{B,3D}^u/E_{K,3D}^t \sim 1$, which agrees with the prediction in Falceta-Gonçalves et al. (2008). Figure 14(c) also shows an increase of the energy ratio $E_{B,3D}^t/E_{K,3D}^t$ as time evolves. However, there is only a slight increase from the time the first sink forms to the time when SFE = 15%, which cannot bring the state to an energy equipartition for models with small $E_{B,3D}^t/E_{K,3D}^t$ values.

### 3.6. The Correction Factor $Q_c$

The assumptions of the DCF method may not be valid in some cases; therefore, the magnetic field strength estimated from the DCF method can deviate from the model magnetic field strength. The correction factor $Q_c$ is defined as the ratio between the model and estimated magnetic field strength and is used to correct for the derived magnetic field strength. Here we derive the correction factors for the magnetic field strength derived from the raw DCF method and the ADF method based on our simulation results. The observational effects (i.e., the beam-smoothing effect of telescopes and the large-scale filtering effect of interferometers) are not considered in the derivation of the correction factors here.

We derive the plane-of-sky mean/total magnetic field strength in the models within spheres or cylinders of different

<figure><img src="image_30.png" /><figcaption>**Figure 14.** The relation between the energy ratio $E_{B,3D}^t/E_{K,3D}^t$ and different physical parameters within spheres of different radii for all simulations models. Different symbols represent models with different initial $\mu$ values. In (a) and (b), different colors correspond to models with different initial Mach numbers. In (c), different colors correspond to different time snapshots.</figcaption></figure>
radii with respect to the most massive sink and divide them by the plane-of-sky mean/total magnetic field strength derived from the raw DCF method or the ADF method in the projected circles in the polarization maps to obtain the correction factors. Figure 15 relates the correction factors for $B_{\text{pos}}^{\text{u,dcf,sa}}$ (see Equation (6)), $B_{\text{pos}}^{\text{tot,dcf}}$ (see Equation (9)), $B_{\text{pos,or,tc,si}}^{\text{u,adf}}$ (see Equation (13)), $B_{\text{pos,or,tc,si}}^{\text{tot,adf}}$ (see Equation (14)), $B_{\text{pos,or,tc}}^{\text{u,adf}}$ and $B_{\text{pos,or,tc}}^{\text{tot,adf}}$, with the measured angular dispersion or the turbulent-to-ordered/-total magnetic field strength ratio derived from the ADF method. The correction factor for $B_{\text{pos}}^{\text{tot,dcf,sa}}$ (see Equation (10)) is very similar to the correction factor for $B_{\text{pos}}^{\text{tot,dcf}}$ and is not shown here.

As shown in Figure 15, overall the magnetic field strength is overestimated by the DCF method. As stated in Section 3.3.3, the angular dispersion in polarization maps do not or only slightly overestimate the turbulent-to-ordered or -total magnetic field strength ratio for radii greater than 0.1 pc, so the overestimation is mostly due to the energy nonequipartition. We find that models with a strong initial magnetic field and models with a weak initial magnetic field show similar ranges of angular dispersion or derived turbulent-to-ordered/-total magnetic field strength ratio. Thus, for a given angular dispersion or derived turbulent-to-ordered/-total magnetic field strength ratio, the correction factor shows large scatter over different models. So, the angular dispersion or derived turbulent-to-ordered/-total magnetic field strength ratio should not be used as indicators for the correction factor.

The correction factor for models with a stronger initial magnetic field is closer to 1, which is because models with a stronger magnetic field are closer to the energy equipartition. Because the DCF method can significantly overestimate the magnetic field strength for weak-field models, we focus on interpreting the correction factors for strong-field models ($\mu = 1.2$ or $\mu = 2$). For strong-field models, it seems that correction factors at radii <0.1 pc show a larger scatter than those at radii >0.1 pc. Overall, the scatter in the correction factor for the uniform magnetic field strength is larger than that for the total magnetic field strength. The correction factor for strong-field models seems to be closer to 1 for larger angular dispersions or turbulent-to-ordered/total magnetic field strength ratios. For strong-field models, the correction factor for $B_{\text{pos}}^{\text{u,dcf,sa}}$ show small scatter at small angular dispersions (e.g., $\delta\phi_{\text{obs}} \lesssim 25^{\circ}$) and radii >0.1 pc (see Figure 15(a)), while the correction factor for $B_{\text{pos}}^{\text{tot,dcf}}$ shows small scatter at radii >0.1 pc (see Figure 15(b)). On the other hand, even for strong-field models, the correction factors for the uniform magnetic field strength derived from the ADF method taking into account the signal integration effect show large scatter regardless of whether the derived turbulent-to-ordered magnetic field strength ratio is large or small, or the radius is large or small (see Figure 15(c)). For the total magnetic field strength derived from the ADF method taking into account the signal integration effect (see Figure 15(d)), the scatter in the correction factor is smaller than that in Figure 15(c) but still spans over one magnitude. The large scatter in Figures 15(c) and (d) may be because the ADF method does not work well on the effect of the line-of-sight signal integration (See Section 3.4.1). So, we also derive the correction factors for the uniform/total magnetic field strength estimated from the ADF method without taking into account the signal integration effect (see Figures 15(e) and (f)). The relations in Figures 15(e) and (f) are

very similar to those in Figures 15(a) and (b). The critical turbulent-to-ordered magnetic field strength ratio is ~0.1 in Figure 15(e), over which the correction factor for $B_{\text{pos,or,tc}}^{\text{u,adf}}$ for strong-field models at radii >0.1 pc can show large scatter. We also compare the correction factors within cylinders of different radii in the models with the corresponding angular dispersions or turbulent-to-ordered/-total magnetic field strength ratios and find that the relations are similar to what are shown in Figure 15.

As mentioned above, the correction factors for strong-field models (with initial $\mu = 1.2$ or $\mu = 2$) are close to 1 and can be confined within a small value range if appropriate radii or angular dispersion selection criteria are applied. The average correction factor for strong-field models at radii >0.1 pc are derived and listed in Table 6. We also directly divide the 3D magnetic field strength within spheres or cylinders in the model by the plane-of-sky magnetic field strength derived from the raw DCF or ADF method and list the results in Table 6. It should be noted that the correction factors in Table 6 are only appropriate when there are other evidence supporting a strong-field case. We notice that the uncertainty for the correction factor of the uniform magnetic field strength derived from the ADF method taking into account the signal integration effect is close to 100%, so we do not recommend using this method to derive the uniform magnetic field strength.

# 4. Discussion

## 4.1. Implication of Angular Dispersions

The DCF method assumes the turbulent-to-ordered or -total magnetic field strength ratios are traced by the angular dispersion of magnetic field orientations. With similar strengths of the turbulent field, the stronger the ordered or total field is, the smaller the turbulent-to-ordered or -total strength ratio would be. Thus, in principle, the angular dispersions would also be expected to be smaller for a stronger field. However, we find that strong-field models and weak-field models do not differ significantly in the value range of angular dispersions, i.e., models with a strong field can have large angular dispersions and models with a weak field can have small angular dispersions. This could probably be partly explained by the uncertainty from the projection effect. Moreover, weak-field and strong-field models do not produce the same probability distribution functions in the density of the structures formed. Because the magnetic field measured is the combination of different structures with different field strengths at different scales in the different models, there would be difficulty in simply associating large dispersions to weak fields and small dispersions to strong fields. Thus, one cannot distinguish between strong- and weak-field cases only with the measured angular dispersions.

The turbulent-to-ordered field strength ratio is expected to be traced by $\delta(\tan \phi)$ or $\delta\phi$ (in the small-angle approximation). When $\phi \sim 90^{\circ}$, $\tan \phi$ would be very sensitive to $\phi$. $\delta(\tan \phi)$ would be dominated by data points with $\phi \sim 90^{\circ}$ and show large scatter (see Section 3.3). Thus, $\delta(\tan \phi)$ generally does not correctly estimate the turbulent-to-ordered field strength ratio. The turbulent-to-ordered field strength ratio could be much greater than 1 for weak-field cases, while the $\delta\phi$ for a random field is ~52° (~0.91). In that situation, the measured $\delta\phi$ could greatly underestimate the turbulent-to-ordered field strength ratio and overestimate the uniform magnetic field strength.
<figure><img src="image_31.png" /><figcaption>a</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>f</figcaption></figure>

**Figure 15.** (a)–(b) Relation between the correction factor within spheres of different radii and the corresponding angular dispersion in polarization maps for all simulations models. (c)–(f) Relation between the correction factor within spheres of different radii and the corresponding turbulent-to-ordered or -total magnetic field strength ratio derived from the ADF method for all simulations models. Different symbols represent models with different initial *μ* values. Different colors correspond to different radii.
<table><caption>**Table 6**
Correction Factors for Strong Magnetic Field Models (Initial $\mu = 1.2$ or $\mu = 2$) at $R > 0.1$ pc</caption><thead><tr><th>Correction Factor</th><th>Spheres</th><th>Cylinders</th></tr></thead><tbody><tr><td>$B_{\text{pos}}^{\text{tot}}/B_{\text{pos}}^{\text{tot,def}}$</td><td>0.64(37%)</td><td>0.41(39%)</td></tr><tr><td>$B_{\text{pos}}^{\text{u}}/B_{\text{pos}}^{\text{u,def,saa}}$</td><td>0.33(25%)</td><td>0.22(40%)</td></tr><tr><td>$B_{\text{pos}}^{\text{tot}}/B_{\text{pos}}^{\text{tot,def,sa}}$</td><td>0.75(44%)</td><td>0.48(45%)</td></tr><tr><td>$B_{\text{pos}}^{\text{u}}/B_{\text{pos,or,tc,si}}^{\text{u,adf}}$</td><td>0.25(83%)</td><td>0.12(113%)</td></tr><tr><td>$B_{\text{pos}}^{\text{tot}}/B_{\text{pos,or,tc,si}}^{\text{tot,adf}}$</td><td>0.49(56%)</td><td>0.29(57%)</td></tr><tr><td>$B_{\text{pos}}^{\text{u}}/B_{\text{pos,or,tc}}^{\text{u,adf b}}$</td><td>0.11(32%)</td><td>0.08(41%)</td></tr><tr><td>$B_{\text{pos}}^{\text{tot}}/B_{\text{pos,or,tc}}^{\text{tot,adf}}$</td><td>0.25(45%)</td><td>0.16(45%)</td></tr><tr><td>$B_{\text{3D}}^{\text{tot}}/B_{\text{pos}}^{\text{tot,def}}$</td><td>0.84(47%)</td><td>0.56(54%)</td></tr><tr><td>$B_{\text{3D}}^{\text{u}}/B_{\text{pos}}^{\text{u,def,saa}}$</td><td>0.34(24%)</td><td>0.23(39%)</td></tr><tr><td>$B_{\text{3D}}^{\text{tot}}/B_{\text{pos}}^{\text{tot,def,sa}}$</td><td>0.99(54%)</td><td>0.66(60%)</td></tr><tr><td>$B_{\text{3D}}^{\text{u}}/B_{\text{pos,or,tc,si}}^{\text{u,adf}}$</td><td>0.41(79%)</td><td>0.24(82%)</td></tr><tr><td>$B_{\text{3D}}^{\text{tot}}/B_{\text{pos,or,tc,si}}^{\text{tot,adf}}$</td><td>0.63(57%)</td><td>0.38(58%)</td></tr><tr><td>$B_{\text{3D}}^{\text{u}}/B_{\text{pos,or,tc}}^{\text{u,adf b}}$</td><td>0.12(28%)</td><td>0.08(37%)</td></tr><tr><td>$B_{\text{3D}}^{\text{tot}}/B_{\text{pos,or,tc}}^{\text{tot,adf}}$</td><td>0.32(55%)</td><td>0.21(60%)</td></tr></tbody></table>

**Notes.** Values in parentheses are the relative uncertainty.

a An additional selection criterion of $\delta\phi_{\text{obs}} \lesssim 25^{\circ}$ ($\sim 0.44$) is adopted.

b An additional selection criterion of $(\langle B_r^2 \rangle / \langle B_0^2 \rangle)_{\text{or,tc}}^{\text{ADF}} \lesssim 0.1$ is adopted.

strength. On the other hand, the turbulent-to-total field strength ratio is expected to be traced by $\delta(\sin\phi)$ or $\delta\phi$ (in the small-angle approximation). Our estimations suggest that the $\delta(\sin\phi)$ or $\delta\phi$ can be well correlated with the turbulent-to-total field strength ratio within a factor of 2 in the simulation models.

Observationally, the angular dispersion is measured from the 2D polarization maps and could be affected by the effect of signal integration along the line of sight. Our results indicate that the $\delta\phi$ or $\delta(\sin\phi)$ in the polarization maps agrees with the $\delta\phi$ or $\delta(\sin\phi)$ in the simulation models within a factor of $\sim 2$, when the radius of the structure is greater than $\sim 0.1$ pc (approximately at the core scale), but the angular dispersion in the polarization maps could be significantly lower than the angular dispersion in the simulation models when radii $< 0.1$ pc. For the effect of signal integration along the line of sight, the ADF method analytically accounts for its effect on the measured angular dispersion, while Cho & Yoo (2016) proposed a measure of the correction factor for this effect from velocity variations. However, our results suggest that the ADF method is not applicable in accounting for this effect in the majority of cases. The CY16 method correctly takes into account this effect at radii $> 0.1$ pc, while it fails at radii $< 0.1$ pc. Because the measured angular dispersion can be significantly underestimated at radii $< 0.1$ pc, and neither the ADF method nor the CY16 method works well on the influence of the signal integration along the line of sight below the core scale, we suggest that the DCF methods should be avoided to be applied on scales smaller than the core scale.

The contribution from the ordered field structure can bias the measured angular dispersion in polarization maps toward larger values. If the field structure is highly ordered, the background field structures may be subtracted by fitting with specific field models (Girart et al. 2006; Myers et al. 2018). The attempts to universally remove the contribution from the large-scale ordered field structure include the spatial filtering method (Pillai et al. 2015), the unsharp masking method (Pattle et al. 2017), and the ADF method. The results of our simple Monte

Carlo simulations suggest that the ADF method can correctly account for the ordered field contribution (see Appendix C). Some comparative studies indicate that the magnetic field strengths derived from the unsharp masking method are larger than those derived from the ADF method in the same region (Houde et al. 2009; Pattle et al. 2017; Liu et al. 2019). There are two factors that may explain the systematic difference between the results of the two methods: the turbulent-to-ordered magnetic field strength ratio might be overestimated by the ADF method because of the overestimated number of turbulent cells along the line of sight, while the angular dispersion might be underestimated by the unsharp masking method because this method does not consider the beam-smoothing effect and the turbulent correlation effect.

The ADF method also is the only method to analytically take into account the turbulent correlation effect, the beam-smoothing effect, and the interferometric filtering effect. Our results suggest that the ADF method does correctly account for the beam-smoothing effect and the interferometric filtering effect, but do not provide evidence on whether the ADF method correctly accounts for the turbulent correlation effect.

## 4.2. Equipartition between the Turbulent Kinetic Energy and the Turbulent Magnetic Energy

The equipartition between the turbulent kinetic energy and the turbulent magnetic energy is a basic assumption of the DCF method. This assumption is expected to be satisfied in sub-Alfvénic ($E_{K,3D}^{\text{u}} < E_{B,3D}^{\text{u}}$) cases, where weak turbulence generates Alfvénic magnetic perturbations on a strong magnetic field (Falceta-Gonçalves et al. 2008). For super-Alfvénic ($E_{K,3D}^{\text{t}} > E_{B,3D}^{\text{u}}$) models, the perturbation on the magnetic field cannot follow the strong turbulent motions. Thus, the turbulent magnetic energy would be lower than the turbulent kinetic energy, and the uniform/total magnetic field strength would be overestimated (Heitsch et al. 2001; Falceta-Gonçalves et al. 2008).

In our calculations, we find that models with stronger magnetic fields tend to be closer to the energy equipartition. However, it could be hard to discern between the strong-field cases and the weak-field cases from observations. Due to the previously discussed uncertainties, the measured angular dispersion in polarization maps does not give us solid information to qualitatively distinguish between the strong- and weak-field cases. Thus, the measured angular dispersion in polarization maps cannot determine how close the state is to the energy equipartition. The polarization-intensity gradient method proposed by Koch et al. (2012) may help to determine whether the field is strong or weak from polarization maps. We will present the application of this method on our simulation results in a separate work. Comparing the large-to-small magnetic field orientations may also be useful to qualitatively discern the strong-field case from the weak-field case (Hull et al. 2014; Zhang et al. 2014). Strong-field cases and weak-field cases may also show observable differences in the magnetic field morphology and density structure.

Massive star formation regions are usually associated with strong turbulence and large velocity dispersions. Low-mass star-forming regions usually have weaker turbulence and smaller velocity dispersions, so they are expected to be more likely in energy equipartition if they have a similar magnetic field strength to that of massive star-forming regions. To investigate whether low-mass star formation region is more
likely in energy equipartition, we perform an additional simulation of a 10 $M_\odot$ spherical dense core with an initial radius of 0.067 pc. The initial temperature, $\mu$, and Mach numbers are 10 K and 1.2 and 1, respectively. We find that the ratio $E_{B,3D}^t/E_{K,3D}^t$ within 0.067 pc is ~1.2 when the first sink forms and ~0.8 when SFE = 15%. For comparison, the ratio $E_{B,3D}^t/E_{K,3D}^t$ within 0.067 pc is ~0.2 for the massive star formation model T10M1MU1 at both time snapshots. Thus, we suggest that the low-mass star-forming region would be closer to energy equipartition for similar magnetic critical parameters. More simulations of low-mass star-forming regions with different initial parameters would be useful to gain more solid conclusions.

One of the purposes for estimating the magnetic field strength is to compare the relative importance of the magnetic field energy and the turbulent energy. The energy equipartition assumption of the DCF method implicitly requires $E_{K,3D}^t = E_{B,3D}^t < E_{B,3D}^t + E_{B,3D}^u = E_{B,3D}^{tot}$, i.e., the turbulent kinetic energy is smaller than the total magnetic energy. If the question of energy equipartition is not properly addressed in the application of the DCF method, the estimated total magnetic field energy would always be larger than the turbulent kinetic energy. Then the relative importance between the magnetic field and the turbulence in the energy balance cannot be properly compared. It also should be noted that the total magnetic energy could be larger than the turbulent kinetic energy in a super-Alfvénic case if there is significant magnetic energy in turbulent forms. So, a super-Alfvénic state does not necessarily mean the turbulence is more important than the magnetic field in the energy balance.

Another main purpose for estimating the magnetic field strength is to compare the relative importance of the magnetic field energy and the gravitational energy. Similarly, if there is a significant amount of magnetic field energy in turbulent forms, the total magnetic energy, rather than the uniform magnetic energy, should be used to compare with the gravitational energy in the energy balance (i.e., replace the uniform magnetic field strength by the total magnetic field strength in the derivation of the mass-to-flux-ratio-to-critical-value or the magnetic virial parameter).

Most previous observational DCF studies applied the DCF method on polarization maps to estimate the uniform magnetic field strength, but none of them considered the problem of energy nonequipartition. In practice, if the estimated Alfvénic Mach number is greater than 1 (super-Alfvénic), the magnetic field strength must have been overestimated due to the effect of energy nonequipartition. On the other hand, the estimated magnetic field strength for strong-field cases and weak-field cases could be degenerate. Even if the estimated Alfvénic Mach number is smaller than 1 (sub-Alfvénic), it does not necessarily ensure that there is equipartition between the turbulent kinetic energy and the turbulent magnetic energy. It is still possible that the magnetic field strength is overestimated by the DCF method and the true situation is super-Alfvénic.

In summary, our simulations suggest that energy equipartition is only achieved in strong-field models (i.e., $\mu = 1.2$ or $\mu = 2$). In weak-field models, energy equipartition generally could not be reached, and the magnetic field strength would be overestimated by the DCF method. On the other hand, low-mass star-forming regions are closer to or more likely in the energy equipartition because they are associated with weaker turbulence than massive star-forming regions. In Section 3.6,

the correction factors accounting for the energy nonequipartition are derived for the strong-field models of massive star formation. If the strong-field scenario can be confirmed independently from observations, these correction factors can be applied to correct for the magnetic field strength estimated with the DCF method. We note that there is a lack of theoretical effort to quantify the energy equipartition problem from observations so far.

### 4.3. Comparison with Previous Simulation Works

There have been several studies on testing the reliability of the DCF method with numerical simulations in the literature (Heitsch et al. 2001; Ostriker et al. 2001; Padoan et al. 2001; Falceta-Gonçalves et al. 2008; Skalidis & Tassis 2021). Here we compare our results with those of the previous numerical works.

Ostriker et al. (2001) performed 3D compressible ideal MHD simulations to present the conditions of self-gravitating giant molecular clouds. Their simulations have a resolution of 256³, a fiducial average temperature of 10 K, and a fiducial average density of 100 cm⁻². Their models considered three different magnitudes of magnetic field strength, but they only derived the correction factor for the plane-of-sky uniform magnetic field strength for a strong-field model (with initial $B_{3D}^u = 14 \ \mu\text{G}$) that is sub-Alfvénic, magnetically subcritical, and with a box length of ~8 pc. They found that the correction factor $B_{\text{pos}}^u/B_{\text{pos}}^{\text{u,dcf,sa}}$ is ~0.5 if the angular dispersion in the polarization map (without removing the ordered field structure) is smaller than 25° (~0.44). Because the model is sub-Alfvénic, it is very likely that the energy equipartition assumption is fulfilled. So, the correction factor reflects the ratio between the angular dispersion and the turbulent-to-ordered magnetic field strength ratio. However, their correction factors are appropriate for parsec-scale molecular clouds at low densities. At smaller scales of molecular clumps and cores, gravity becomes prominent and there could be a significant contribution from the ordered field structure to the directly measured angular dispersion in polarization maps. Our simulation results indicate that statistically, the average ratio between the directly measured angular dispersion in polarization maps and the turbulent-to-ordered magnetic field strength ratio is ~0.25 (see Table 4) in clump and core scales if $\delta\phi_{\text{obs}} < 25°$, which is smaller than the value reported by Ostriker et al. (2001).

Padoan et al. (2001) performed 3D compressible MHD simulations of a super-Alfvénic and magnetically supercritical molecular cloud with a size of 6.25 pc and a numerical resolution of 128³. They selected three self-gravitating ~1 pc clumps (they used “cores” in their terminology) and estimated an average correction factor to be ~0.4 for the plane-of-sky uniform magnetic field strength ($B_{\text{pos}}^u/B_{\text{pos}}^{\text{u,dcf,sa}}$). Because their model is super-Alfvénic, it is likely that the energy is not in equipartition and the correction factor derived in their work includes the contribution from energy nonequipartition. Although their estimated $\delta\phi_{\text{obs}}$ are in the range of 0.37–0.73 (21°–42°), the correction factors derived by Padoan et al. (2001) show small scatter, which may be because they only have several estimations. As indicated by our simulation results, $\delta\phi_{\text{obs}}$ could be much smaller than $(B_{\text{pos}\perp}^t/B_{\text{pos}}^u)$ and the ratio between the two quantities could have large scatter when $\delta\phi_{\text{obs}} > 25°$. The DCF method should be avoided to derive the
uniform magnetic field strength when the angular dispersion in the polarization map is large.

Heitsch et al. (2001) presented a set of 3D ideal MHD scale-free simulations with constant turbulent energy input to simulate the interior of molecular clouds. They considered both super-Alfvénic and sub-Alfvénic models with resolutions ranging from 128³ to 512³. They used the same equation as Equation (5) in this paper to estimate the plane-of-sky uniform magnetic field strength. They found that the plane-of-sky uniform magnetic field strength derived from this equation is overestimated by a factor of 2–3 for strong-field models and shows significant scatter for weak fields. They also proposed that the 3D total (rms) magnetic field strength is estimated as $B_{3D}^{\text{tot,Hei01}} \sim \sqrt{\mu_0 \rho} (\delta v_{\text{los}} / \delta(\tan \phi))(1 + 3\delta(\tan \phi)^2)^{1/2}$ if the turbulent magnetic energy is isotropic and the mean magnetic field orientation is along the plane of the sky. They found that the 3D total magnetic field strength derived from their equation is overestimated by a factor of no more than 3 for strong-field models, but the field strength could be overestimated by up to 3 orders of magnitude for weak fields. After correction for energy nonequipartition, they found both the estimated uniform and total magnetic field strength match the model field strength by a factor between 1 and 1.5 for strong-field models but still yield large deviations from the weak-field estimates. For comparison, our simulation results show large scatter for estimations including $\delta(\tan \phi)$ regardless of whether the magnetic field is strong or weak, so we do not recommend using $\delta(\tan \phi)$ as an estimate of the turbulent-to-ordered magnetic field strength ratio. On the other hand, similar to the results in Heitsch et al. (2001), our Equations (9) and (10) also only slightly overestimate the total magnetic field strength for strong-field models and overestimate the total magnetic field strength by up to several orders of magnitude for weak-field models (see Figure 15). Unlike the Heitsch et al. (2001) equation where the total magnetic field strength is significantly overestimated for weak-field cases, our Equations (9) and (10) only slightly overestimate the total magnetic field strength regardless of whether the magnetic field is strong or weak (see Figure 7) if the energy equipartition assumption is fulfilled.

Falceta-Gonçalves et al. (2008) presented 3D ideal high-resolution (512³) MHD simulations with continuous turbulence injection to study the polarization emission in molecular clouds. They considered both super-Alfvénic and sub-Alfvénic simulations that are non-self-gravitating and scale-independent. They proposed that the plane-of-sky total magnetic field strength is estimated as $B_{\text{pos}}^{\text{tot,Fal08}} \sim \sqrt{\mu_0 \rho} \delta v_{\text{los}} / \tan(\delta\phi)$, where they assume $(B_{\text{pos}\perp}^{\text{tot}}/B_{\text{pos}}^{\text{tot,Fal08}} \sim \delta(\tan\phi) \sim \tan(\delta\phi)$. However, we note that they might have approximated the turbulent-to-ordered magnetic field strength ratio with the turbulent-to-total magnetic field strength ratio because $B_{\text{pos}\perp}^{\text{tot}}/B_{\text{pos}}^{\text{tot}} \sim \delta(\sin\phi)$ and $B_{\text{pos}\perp}^{\text{tot}}/B_{\text{pos}}^{\text{tot}} \sim \delta(\tan\phi)$ there. Due to this approximation, we refrain from discussing their results for the correction factors. We note that $(B_{\text{pos}\perp}^{\text{tot}}/B_{\text{pos}}^{\text{tot,Fal08}} \sim \tan(\delta\phi)$ is only valid when there is a small-angle approximation $\tan(\delta\phi) \sim \delta\phi$.

Recently, Skalidis & Tassis (2021) used 3D ideal scale-free sub-Alfvénic MHD simulations with no gravity to investigate the non-Alfvénic modes in the interstellar medium. In contrast to the assumption of the raw DCF method that there is energy equipartition between the turbulent kinetic energy and turbulent magnetic energy, they assumed that the turbulent kinetic energy equals the cross-term of the magnetic energy. They proposed that the uniform plane-of-sky magnetic field strength

could be estimated with $B_{\text{pos}}^{\text{u,Ska20}} \sim \sqrt{\mu_0 \rho / 2} \delta v_{\text{los}} / \sqrt{\delta\phi}$ in small-angle approximation. They claimed that their new equation can yield an accurate estimate of the true uniform field strength within a mean relative deviation of 17%. However, for different ratios between the angular dispersion and the turbulent-to-ordered field strength ratio and the different ratios between the turbulent kinetic energy and the cross-term of magnetic energy, the estimated correction factor may degenerate with each other. As they did not compare the two relations separately and they only have five estimates of the correction factors for models with an Alfvénic Mach number of 0.7, it is unknown whether their correction factor of $\sim 1$ is a coincidence because there might not be equipartition between the turbulent kinetic energy and cross-term magnetic energy in sub-Alfvénic models.

## 5. Summary

In this paper we present a set of 3D high-resolution ideal MHD simulations and radiative transfer simulations of self-gravitating clustered massive star-forming regions to study the accuracy of the DCF method, which is the most widely used method to estimate the plane-of-sky magnetic field strength with statistics of polarization position angles from polarization maps. There are several assumptions in the DCF method: there is a mean-field component that is prominent (i.e., small angle approximation), there is energy equipartition between the turbulent kinetic energy and the turbulent magnetic energy, the turbulence is isotropic, and the perturbation on the magnetic field is traced by angular dispersions in polarization maps. The estimated angular dispersion in polarization maps is subject to uncertainties from the line-of-sight signal integration, the turbulent correlation effect, the contribution from ordered field structure, the beam-smoothing effect of telescopes, and the additional large-scale filtering effect of interferometric observations. We tested the assumptions of the DCF method and the factors affecting the measured angular dispersion step by step with our simulation results. The main conclusions are as follows:

1. The magnetic critical parameter (mass-to-flux ratio over its critical value) of a simulation is determined by its initial value. There is no strong relation between the magnetic critical parameter and the spatial scale in our simulations. On the other hand, simulations with different initial Mach numbers show similar levels of turbulent velocity dispersion, which indicates the initial environmental turbulence decays very quickly and the turbulence is most likely generated from star formation activities when the self-gravity is dominant.

2. If the energy equipartition between the turbulent magnetic energy and the turbulent kinetic energy is fulfilled and the directly measured angular dispersion in the polarization map is smaller than 25°, the average correction factor between the model and estimated plane-of-sky uniform magnetic field strength would be $\sim 0.25$ at clump and core scales, which is smaller than the correction factor of $\sim 0.5$ found in previous numerical studies at larger scales. When the energy equipartition is met, the average correction factor between the model and estimated plane-of-sky total magnetic field strength would be $\sim 0.81$ or $\sim 0.57$, respectively, depending on whether the angular dispersion is in the forms of $\delta\phi$ or $\delta(\sin\phi)$.
3. Our simulations indicate that models of stronger magnetic fields are closer to energy equipartition. Because low-mass star-forming regions are associated with weaker turbulence, they are closer to or more likely in energy equipartition compared to massive star-forming regions with similar magnetic critical parameters. When the turbulent magnetic energy is smaller than the turbulent kinetic energy, the magnetic field strength would be overestimated by the DCF method. In our simulations, only strong-field models are close to energy equipartition. We estimate and catalog the correction factors for the uniform and total magnetic field strength accounting for the energy nonequipartition for strong-field massive star formation models with initial magnetic critical parameters of 1.2 or 2. These correction factors are applicable to observational data to determine the correct magnetic field strength at clump and core scales in strong-field cases.

4. The ADF method (Falceta-Gonçalves et al. 2008; Hildebrand et al. 2009; Houde et al. 2009, 2016) analytically takes into account the factors that affect the estimated angular dispersion. We test the accuracy of the ADF method with our simulation results and find that the ADF method does correctly account for the ordered field structure, the beam-smoothing effect, and the interferometric filtering effect. However, the ADF method might not be applicable to account for the line-of-sight signal integration effect from observational data in most cases.

5. We test the Cho & Yoo (2016) method and find that this method correctly accounts for the effect of line-of-sight signal integration at >0.1 pc, but fails at <0.1 pc.

6. The line-of-sight signal integration only underestimates the angular dispersion by a factor of no more than 2 at radii >0.1 pc, while the angular dispersion may be significantly underestimated at scales below ~0.1 pc (the core scale). Because the signal integration effect may not be correctly accounted for at <0.1 pc, the DCF method should be avoided from being applied to scales below the core scale.

7. The DCF method is valid for deriving the uniform magnetic field strength only when the directly measured angular dispersion in the polarization map is small (i.e., $\delta\phi_{obs} < 25^\circ$). There is no limitation on the angular dispersion for the DCF method to derive the total magnetic field strength. Overall, the ratio between the model and estimated uniform magnetic field strength shows larger scatter than the ratio for the total magnetic field strength. To properly compare the relative importance between the magnetic field and other forces (e.g., turbulence and gravity), the total magnetic field strength

should be used instead of the uniform magnetic field strength if there is a significantly turbulent magnetic field.

We are indebted to the anonymous referee whose constructive comments improved the presentation and clarity of the paper. J.L. thanks Prof. Martin Houde for helpful discussions on the principle of the angular dispersion function method and Dr. Robert Brauer for the help in the installation of POLARIS. J.L. also thanks Dr. Xing Lu, Mr. Yu Cheng, and Dr. Shanghuo Li for helpful discussions on the components of velocity dispersions. We thank Dr. Stefan Reissl for helpful discussions on the POLARIS code. K.Q. and J.L. are supported by the National Key R&D Program of China No. 2017YFA0402600. K.Q. and J.L. acknowledge the support from the National Natural Science Foundation of China (NSFC) through grants U1731237, 11590781, and 11629302. J.L. acknowledges the support from the program of the China Scholarship Council (No. 201806190134), the Smithsonian Astrophysical Observatory pre-doctoral fellowship, and the support from the EAO Fellowship Program under the umbrella of the East Asia Core Observatories Association. This study was enabled by the MagneticYSOs project, a European Research Council (ERC) funded collaboration, under the European Union Horizon 2020 research and innovation program (grant agreement N. 679937, PI: Maury). B.C. acknowledges support by the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU co-funded by CEA and CNES. The 3D RAMSES simulations were performed using HPC resources from GENCI-CINES (Grant A0080407247). This research made use of APLpy, an open-source plotting package for Python (Robitaille & Bressert 2012), Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013), and Matplotlib, a Python 2D plotting library for Python (Hunter 2007).

*Software:* APLpy (Robitaille & Bressert 2012), Astropy (Astropy Collaboration et al. 2013), Matplotlib (Hunter 2007).

## Appendix A
### Examples of Synthetic Images

Figure 16 shows the synthetic images of model T10M6MU200 at SFE = 15% as an example of the synthetic observations of a model with a weak magnetic field and strong turbulence. Figure 17 shows the synthetic images of model T10M1MU1 at SFE = 15% as an example of the synthetic observations of a model with a strong magnetic field and weak turbulence. Figure 18 shows the synthetic ALMA images of model T10M1MU1 at SFE = 15% and at the zy plane as an example of the simulated ALMA observations. All images are centered at the projected position of the most massive sink.
<figure><img src="image_37.png" /><figcaption>**Figure 16.** Synthetic observations of model T10M6MU200 at SFE = 15% and at three orthogonal projections. Gray scales indicate Stokes *I* emissions. Line segments represent the orientation of the magnetic field (rotating the orientation of the observed linear polarization by 90°). Line segments length are arbitrary.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>**Figure 17.** Synthetic observations of model T10M1MU1 at SFE = 15% and at three orthogonal projections. Gray scales indicate Stokes *I* emissions. Line segments represent the orientation of the magnetic field (rotating the orientation of the observed linear polarization by 90°). Line segments length are arbitrary.</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>**Figure 18.** Synthetic ALMA observations of model T10M1MU1 at SFE = 15% and at the zy plane. Gray scales indicate Stokes *I* emissions. Regions with negative *I* values are not shown. Line segments represent the orientation of the magnetic field (rotating the orientation of the observed linear polarization by 90°). Line segment lengths are arbitrary. (a) AC1. (b) AC3. (c) AC5. The size of the synthesized beam is shown in the lower-left corner of each panel. The size of the maximum recoverable scale is indicated in the lower-right corner of each panel.</figcaption></figure>
## Appendix B
### Effect of Grid Size

The DCF method assumes there are magnetic field perturbations around the mean-field component. If the perturbation only dominates at scales smaller than a turbulent correlation scale $l_\delta$, increasing the pixel size at scales smaller than $l_\delta$ could increase the angular dispersion from field perturbation and increasing the pixel size at scales larger than $l_\delta$ does not change the turbulent contribution to the angular dispersion (Heitsch et al. 2001; Hildebrandt et al. 2009). In the presence of field curvature, the pixel size may also affect the directly measured angular dispersion if the complex small-scale field structure is smoothed in large pixel size. Here we investigate how the grid (2D pixel or 3D cell) size can affect the directly measured angular dispersions from 2D polarization maps and 3D simulation models.

The 2D polarization maps and 3D simulation parameter value grids used in this paper are gridded with sizes of 0.13 (\~0.0013 pc or \~260 au). Due to the limitation of computer disk space and time, it is hard for us to perform the analysis with grid sizes smaller than 0.13. In a similar way to that in Section 3.3, we measure the angular dispersions from binned polarization maps and magnetic field parameter grids with bin sizes of 5, 13, and 21 and compare them with the angular dispersion from unbinned polarization maps and field parameter grids.

Figure 19(a) compares the angular dispersions within spheres of different radii measured from magnetic field parameter grids of different cell sizes for all simulation models. We do not estimate the angular dispersion if the number of cells is smaller than 50 within a sphere. For most of the estimations, the ratio between the binned angular dispersions and unbinned angular dispersions is 0.5–1, which could be an effect of smoothing below $l_\delta$ (i.e., $l_\delta > 0.0013$ pc) and/or because the ordered field structure is smoothed. For a small number of estimations, the binned angular dispersions are larger than the unbinned angular dispersions, which could be because of combined effects of $l_\delta \lesssim 0.0013$ pc and random measurement scatter.

Figure 19(b) compares the angular dispersions within circles of different radii measured from polarization maps of different pixel sizes for all simulation models. We do not estimate the angular dispersion if the number of pixels is smaller than 50 within a circle. The binned angular dispersions generally agree with the unbinned angular dispersions with small scatter. There is only a slight tendency for the binned angular dispersions to underestimate the unbinned angular dispersions. The scatter in Figure 19(b) is much smaller than that in Figure 19(a), which might be because the line-of-sight length is much larger than the turbulent correlation length and the polarization position angles in polarization maps are smoothed along the line of sight.

Calculation of the ADFs are very time consuming if there are large numbers of pixels included in the analysis, so we only test a few cases. We find that in the cases we test, the derived $\langle(B_1^2)\rangle/\langle(B_0^2)\rangle_{\text{or,te}}^{\text{adf}}$ for different pixel sizes agree with each other with small scatter. However, the derived $\langle(B_1^2)\rangle/\langle(B_0^2)\rangle_{\text{or,te,si}}^{\text{adf}}$ for different pixel sizes can disagree with each other by more than one order of magnitude, which might be because larger pixel sizes are insufficient to trace the turbulent correlation scale. Thus, we conclude that the ADF method may not correctly

<figure><img src="image_40.png" /><figcaption>Figure 19. Comparison between the binned angular dispersions and the unbinned angular dispersions in (a) simulations grids and (b) polarization maps. Different colors indicate different bin sizes. The solid lines indicate the 1:1 relation and the dashed lines indicate the 1:2 and 2:1 relations.</figcaption></figure>

account for the signal integration along the line of sight if an inappropriate pixel size is used.

## Appendix C
### Accuracy of the ADF Method on Accounting for the Ordered Field Structure

Here we independently investigate whether the ADF method can accurately account for the contribution from the large-scale field structure with simple Monte Carlo simulations of randomly generated Gaussian angular dispersions on modeled ordered fields. Similar to the approach in the Appendix of Liu et al. (2019), we first generate a set of underlying parabola models with the form

$$y = g + gCx^2. \qquad (C1)$$

We only investigate three different curvatures: $C=1$, $C=0.1$, and $C=0.01$. Then we apply different magnitudes of Gaussian angular dispersions in the forms of $\delta\phi$, $\delta(\sin\phi)$, or $\delta(\tan\phi)$ to the underlying field models. The angular dispersions for
<figure><img src="image_41.png" /><figcaption>Figure 20. The modeled B segments with an angular dispersion of $\delta\phi = 0.1$ are shown in blue. The interval of B segments is 1 pixel. Magenta curves denote the underlying parabola field models with $C = 0.01$. Line segments are of unit length. The circles mark the regions within which we derive the ADFs.</figcaption></figure>

randomly distributed $\phi$ and $\sin\phi$ are ~0.91 (~52°) and ~0.58, respectively. An example of the field model with $C = 0.01$ and $\delta\phi = 0.1$ (~6°) is shown in Figure 20.

The ADF only accounting for the contribution from the large-scale ordered field is given by Hildebrand et al. (2009) and Houde et al. (2009, 2016):

$$1 - \langle \cos[\Delta\Phi(l)] \rangle \approx a_2 l^2 + \left( \frac{\langle B_t^2 \rangle}{\langle B^2 \rangle} \right)_{\text{or}}^{\text{adf}}, \quad (\text{C2})$$

where $\langle(B_t^2)/\langle B^2 \rangle\rangle_{\text{or}}^{\text{adf}} = (\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}}/((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}} + 1)$ is the turbulent-to-total energy ratio and $\langle(B_t^2)/\langle B^2 \rangle\rangle_{\text{or}}^{\text{adf}}$ is the turbulent-to-ordered energy ratio. Because the average of $\cos[\Delta\Phi(l)]$ for a randomly distributed $\Delta\Phi(l)$ is ~0.64 and $a_2'$ is defined as positive values, the maximum values of the derived turbulent-to-ordered strength ratio $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ and turbulent-to-total strength ratio $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ are ~0.76 and ~0.6, respectively. We derive the ADFs within radii of 10, 15, and 20 pixels with respect to the center of each model and fit each ADF with Equation (C2) in a similar way to that in Section 3.4.1.

Figure 21(a) shows the derived turbulent-to-ordered strength ratio $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ compared to the injected angular dispersion in forms of $\delta\phi$ and $\delta(\tan\phi)$. As shown by Figures 21(a), the derived $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ may overestimate the injected angular dispersion by a factor of more than 2 when the $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ is smaller than ~0.2. Other than the overestimation at small angles, the derived $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ generally agrees with the injected $\delta\phi$ within a factor of 2. However, $\delta(\tan\phi)$ could be very large when $\phi$ approaches -90° or 90°. Thus, the derived $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ would become insensitive to $\delta(\tan\phi)$ and greatly underestimate $\delta(\tan\phi)$ as $\delta(\tan\phi)$ increases.

<figure><img src="image_42.png" /><figcaption>Figure 21. (a) Estimated turbulent-to-ordered field strength ratio vs. injected angular dispersion in forms of $\delta\phi$ (black) and $\delta(\tan\phi)$ (blue). The horizontal dashed line corresponds to $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2} = 0.86$. The vertical dashed line corresponds to $\delta\phi = 0.91$. (b) Estimated turbulent-to-total field strength ratio vs. injected angular dispersion in forms of $\delta\phi$ (black) and $\delta(\sin\phi)$ (blue). The horizontal dotted line corresponds to $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2} = 0.65$. The vertical dotted lines correspond to $\delta\phi = 0.91$ and $\delta(\sin\phi) = 0.58$. In both (a) and (b), the solid lines indicate the 1:1 relation, and the dashed lines indicate the 1:2 and 2:1 relations.</figcaption></figure>

Figure 21(b) shows the derived turbulent-to-total strength ratio $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ compared to the injected angular dispersion in the forms of $\delta\phi$ and $\delta(\sin\phi)$. Similarly, the derived $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ may overestimate the injected angular dispersion by a factor of more than 2 when the $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ is smaller than ~0.2. Other than the overestimation at small angles, the derived $((\langle B_t^2 \rangle/\langle B^2 \rangle)_{\text{or}}^{\text{adf}})^{1/2}$ generally agrees with the injected angular dispersion within a factor of 2.

We also measured the raw angular dispersions, $\delta\phi_{\text{obs}}$, $\delta(\sin\phi_{\text{obs}})$, and $\delta(\tan\phi_{\text{obs}})$, in a similar way to that in Section 3.3.2. Figure 22 compares the estimated angular dispersions with the injected angular dispersions. Generally, the estimated angular dispersions are greater than the injected angular dispersions because of the contribution from a large-scale field structure. As shown in Figures 22(a) and (b), the
<figure><img src="image_43.png" /><figcaption>Figure 22. Estimated angular dispersions compared to injected angular dispersions. (a) δφobs vs. δφ. (b) δ(sin φobs) vs. δ(sin φ). (c) δ(tan φobs) vs. δ(tan φ).</figcaption></figure>

direct estimation can significantly overestimate the angular dispersion for small injected angular dispersions. Figure 22(c) indicates that the estimated δ(tan φobs) generally would greatly

overestimate the injected angular dispersion not only because of the contribution from a large-scale field structure but also because δ(tan φ) is very sensitive to φ when φ ~ ±90°.

## Appendix D
### Projection Effect for Magnetic Fields and Velocity Fields

Ultimately we want to obtain the 3D magnetic field strength, while the DCF method only gives information about the plane-of-sky magnetic field component. The line-of-sight magnetic field strength can be determined by Zeeman observations. However, the tracers of Zeeman splitting are more sensitive in low-density materials (Crutcher 2012), and we lack the methods to determine the line-of-sight field component in high-density regions.

Alternatively, Crutcher et al. (2004) proposed that there is a statistical relation between the 3D uniform magnetic field strength $B_{3D}^u$ and the plane-of-sky uniform magnetic field strength $B_{pos}^u$:

$$B_{3D}^{u} = \frac{4}{\pi} B_{pos}^{u} \qquad (D1)$$

Similarly, we propose that there is also a statistical relation between the 3D total magnetic field strength $B_{3D}^{tot}$ and the plane-of-sky total magnetic field strength $B_{pos}^{tot}$:

$$B_{3D}^{tot} = \sqrt{\frac{3}{2} B_{pos}^{tot}} \qquad (D2)$$

because $(B_{3D}^{tot})^2 = (B_{los}^{tot})^2 + (B_{pos\perp}^{tot})^2 + (B_{pos\parallel}^{tot})^2$.

Figure 23(a) shows the ratio between the 3D uniform magnetic field strength and the plane-of-sky uniform magnetic field strength within spheres of different radii in all simulation models. It turns out that the ratio between the 3D and plane-of-sky uniform field strength is below $4/\pi$ and only slightly greater than 1 if the line of sight (i.e., y or z) is perpendicular to the initial magnetic field orientation (i.e., x). On the other hand, if the line of sight is along the initial magnetic field orientation, the ratio between the 3D and plane-of-sky uniform field strength could be greater than 10 and the ratio seems to be higher for larger radii. Thus, we conclude that this statistical relation for the uniform magnetic field strength could be very unreliable if the mean magnetic field orientation is along the line of sight.

Figure 23(b) shows the ratio between the 3D total magnetic field strength and the plane-of-sky total magnetic field strength within spheres of different radii in all simulation models. For most estimations, the ratio between the 3D and plane-of-sky total field strength is around $\sqrt{3}/2$. The strength ratio for lines of sight perpendicular to the initial magnetic field orientation tends to be slightly lower than $\sqrt{3}/2$. The strength ratio for the line of sight along the initial magnetic field orientation tends to be higher than $\sqrt{3}/2$, with an upper value of ~2.3 (smaller than $2 \times \sqrt{3}/2 \approx 2.4$). Thus, we conclude that the statistical relation for the total magnetic field strength can give reliable estimates of the 3D total magnetic field strength that is accurate within a factor of 2.

One of the assumptions of the DCF method is isotropic turbulence, which requires

$$\delta v_{3D,turb} = \sqrt{3} \delta v_{los,turb} \qquad (D3)$$

The turbulence may be highly anisotropic if the magnetic field is strong (Goldreich & Sridhar 1995). As stated in Section 3.1,
<figure><img src="image_44.png" /><figcaption>**Figure 23.** (a) The ratio between the 3D uniform magnetic field strength and the plane-of-sky uniform magnetic field strength within spheres of different radii for all models. (b) The ratio between the 3D total magnetic field strength and the plane-of-sky total magnetic field strength within spheres of different radii for all models. Different symbols correspond to different projections. Different colors represent models with different radii.</figcaption></figure>

the rotation can have a significant contribution to the total velocity dispersion. So we remove the contribution from the rotation in the velocity dispersion. The 3D turbulent velocity dispersion is estimated as $\delta v_{3D,turb} = (\delta v_{3D}^2 - v_{\text{t,rot}}^2)^{0.5}$ and the line-of-sight turbulent velocity dispersion is estimated as $\delta v_{\text{los,turb}} = (\delta v_{\text{los}}^2 - v_{\text{t,rot}}^2/3)^{0.5}$ (see Section 3.1 for the derivation of $v_{\text{t,rot}}$). We show the ratio between the 3D turbulent velocity dispersion and the line-of-sight turbulent velocity dispersion within spheres of different radii in all simulation models in Figure 24(a). For most estimations, the 3D to the line-of-sight ratio of turbulent velocity dispersion is between ~1 and ~3.5. Only for a few estimations is the line-of-sight turbulent velocity dispersion much smaller than the 3D turbulent velocity dispersion. In our simulations, the models with initial $\mu = 1.2$ or $\mu = 2$ are trans-Alfvénic and other models are super-Alfvénic. However, there is no significant difference between the 3D to line-of-sight velocity dispersion

<figure><img src="image_45.png" /><figcaption>**Figure 24.** (a) The ratio between the 3D turbulent velocity dispersion and the line-of-sight turbulent velocity dispersion within spheres of different radii for all models. (b) The ratio between the 3D turbulent magnetic field strength and the line-of-sight turbulent magnetic field strength within spheres of different radii for all models. Different symbols correspond to different lines of sight. Different colors represent models with different initial $\mu$ values.</figcaption></figure>

ratios for models with different initial $\mu$ parameters. Nor do we see a difference in the ratios of different projections. Thus, we conclude that the assumption of isotropic turbulence generally holds in our simulations, and the relation $\delta v_{3D,turb} = \sqrt{3} \delta v_{\text{los,turb}}$ is accurate within a factor of ~2.

The turbulent magnetic field would also be expected to be isotropic if the turbulence is isotropic. We show the ratio between the 3D turbulent magnetic field strength and the line-of-sight turbulent magnetic field strength within spheres of different radii in all simulation models in Figure 24(b). The relation is quite similar to what is found for the velocity dispersion but shows slightly less scatter.

ORCID iDs

Junhao Liu (刘峻豪) https://orcid.org/0000-0002-4774-2998
Qizhou Zhang https://orcid.org/0000-0003-2384-6589
Keping Qiu https://orcid.org/0000-0002-5093-5088
References

Andersson, B.-G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501

Arce, H. G., Shepherd, D., Gueth, F., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. Arizona Press), 245

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Beuther, H., Schilke, P., Sridharan, T. K., et al. 2002, A&A, 383, 892

Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015

Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 113

Cho, J., & Yoo, H. 2016, ApJ, 821, 21

Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11

Commerçon, B., Hennebelle, P., & Henning, T. 2011, ApJL, 742, L9

Crutcher, R. M. 2012, ARA&A, 50, 29

Crutcher, R. M., & Kemball, A. J. 2019, FrASS, 6, 66

Crutcher, R. M., Nutter, D. J., Ward-Thompson, D., et al. 2004, ApJ, 600, 279

Davis, L. 1951, PhRv, 81, 890

Davis, L., & Greenstein, J. L. 1949, PhRv, 75, 1605

Falceta-Gonçalves, D., Lazarian, A., & Kowal, G. 2008, ApJ, 679, 537

Fontani, F., Commerçon, B., Giannetti, A., et al. 2016, A&A, 593, L14

Fontani, F., Commerçon, B., Giannetti, A., et al. 2018, A&A, 615, A94

Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371

Girart, J. M., Rao, R., & Marrone, D. P. 2006, Sci, 313, 812

Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763

González-Casanova, D. F., & Lazarian, A. 2017, ApJ, 835, 41

Heitsch, F., Zweibel, E. G., Mac Low, M.-M., et al. 2001, ApJ, 561, 800

Hennebelle, P., Commerçon, B., Joos, M., et al. 2011, A&A, 528, A72

Hildebrand, R. H. 1988, ApL&C, 26, 263

Hildebrand, R. H., Kirby, L., Dotson, J. L., Houde, M., & Vaillancourt, J. E. 2009, ApJ, 696, 567

Houde, M., Hull, C. L. H., Plambeck, R. L., et al. 2016, ApJ, 820, 38

Houde, M., Vaillancourt, J. E., Hildebrand, R. H., Chitsazzadeh, S., & Kirby, L. 2009, ApJ, 706, 1504

Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, ApJS, 213, 13

Hunter, J. D. 2007, CSE, 9, 90

Koch, P. M., Tang, Y.-W., & Ho, P. T. P. 2012, ApJ, 747, 79

Krumholz, M. R., Klein, R. I., McKee, C. F., et al. 2009, Sci, 323, 754

Lazarian, A. 2007, JQSRT, 106, 225

Lazarian, A., & Hoang, T. 2007, MNRAS, 378, 910

Liu, J., Qiu, K., Berry, D., et al. 2019, ApJ, 877, 43

Liu, J., Qiu, K., Wyrowski, F., et al. 2018, ApJ, 860, 106

Liu, J., Zhang, Q., Qiu, K., et al. 2020, ApJ, 895, 142

Mac Low, M.-M., & Klessen, R. S. 2004, RvMP, 76, 125

Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 500, 259

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

Mouschovias, T. C., & Spitzer, L. 1976, ApJ, 210, 326

Mouschovias, T. C., Tassis, K., & Kunz, M. W. 2006, ApJ, 646, 1043

Myers, P. C., Basu, S., & Auddy, S. 2018, ApJ, 868, 51

Myers, P. C., & Goodman, A. A. 1991, ApJ, 373, 509

Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980

Padoan, P., Goodman, A., Draine, B. T., et al. 2001, ApJ, 559, 1005

Pattle, K., Ward-Thompson, D., Berry, D., et al. 2017, ApJ, 846, 122

Pillai, T., Kauffmann, J., Tan, J. C., et al. 2015, ApJ, 799, 74

Pudritz, R. E., Rogers, C. S., & Ouyed, R. 2006, MNRAS, 365, 1131

Reissl, S., Wolf, S., & Brauer, R. 2016, A&A, 593, A87

Robitaille, T., & Bressert, E. 2012, APLpy: Astronomical Plotting Library in Python, Astrophysics Source Code Library, ascl:1208.017

Shang, H., Allen, A., Li, Z.-Y., et al. 2006, ApJ, 649, 845

Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23

Skalidis, R., & Tassis, K. 2021, A&A, 647, A186

Teyssier, R. 2002, A&A, 385, 337

Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJL, 489, L179

Valdivia, V., Maury, A., Brauer, R., et al. 2019, MNRAS, 488, 4897

Weiss, A., Hillebrandt, W., Thomas, H.-C., & Ritter, H. 2004, Cox and Giuli's Principles of Stellar Structure (Cambridge: Cambridge Scientific Publishers)

Zhang, Q., Hunter, T. R., Brand, J., et al. 2005, ApJ, 625, 864

Zhang, Q., Qiu, K., Girart, J. M., et al. 2014, ApJ, 792, 116

Zhang, Q., Wang, Y., Pillai, T., et al. 2009, ApJ, 696, 268

Zweibel, E. G. 1990, ApJ, 362, 545