Magnetic field strength from turbulence theory (I): Using differential measure approach (DMA)

The mean plane-of-sky magnetic field strength is traditionally obtained from the combination of polarization and spectroscopic data using the Davis-Chandrasekhar-Fermi (DCF) technique. However, we identify the major problem of the DCF to be its disregard of the anisotropic character of MHD turbulence. On the basis of the modern MHD turbulence theory we introduce a new way of obtaining magnetic field strength from observations. Unlike the DCF, the new technique uses not the dispersion of the polarization angle and line of sight velocities, but increments of these quantities given by the structure functions. To address the variety of the astrophysical conditions for which our technique can be applied, we consider the turbulence in both media with magnetic pressure larger than the gas pressure corresponding e.g. to molecular and the gas pressure larger than the magnetic pressure corresponding to the warm neutral medium. We provide general expressions for arbitrary admixture of Alfv\'en, slow and fast modes in these media and consider in detail the particular cases relevant to diffuse media and molecular clouds. We successfully test our results using synthetic observations obtained from MHD turbulence simulations. We demonstrate that our Differential Measure Approach (DMA), unlike the DCF, can be used to measure the distribution of magnetic field strengths, can provide magnetic field measurements with limited data and is much more stable in the presence of large scale variations induces of non-turbulent nature. In parallel, our study uncover the deficiencies of the earlier DCF research.


INTRODUCTION
The role of magnetic fields in astrophysics is difficult to overestimate. Magnetic force is the second most important force in the present day Universe after gravity. The magnetic field plays an important role at different stages of star formation (e.g. Mestel & Spitzer 1956;Galli et al. 2006;Mouschovias et al. 2006;Johns-Krull 2007). In view of astrophysical flows with large Reynolds numbers the magnetic fields are turbulent (see Elmegreen & Scalo 2004;McKee & Ostriker 2007, McKee & Stone 2021. The evidence of turbulent magnetic field is coming from observations of density structure of the interstellar medium (e.g. Armstrong et al. 1995;Chepurnov & Lazarian 2009) and velocity fluc-alazarian@facstaff.wisc.edu kyuen@astro.wisc.edu pogosyan@ualberta.ca tuation studies (Larson 1981;Heyer & Brunt 2004;Chepurnov & Lazarian 2010, Chepurnov et al. 2016, Yuen et.al 2021Yuen et al. 2022c,a) and synchrotron polarization studies (Lazarian & Yuen 2018b;Zhang et al. 2020b). In particular the detailed properties of the turbulent velocity can be obtained using using the theory of mapping of velocity fluctuations in real space into the fluctuations of intensity in the Position-Position-Velocity (PPV) space in Lazarian & Pogosyan (2000) and subsequent publications (Lazarian & Pogosyan 2004, 2006. The theory is used in Yuen et.al (2021) to develop Velocity Decomposition Algorithm (VDA) approach that further improves the mapping of velocity statistics. The grand velocity cascade spanning from the large galactic scales (10 4 to 10 −2 pc reported in Yuen et al. (2022a) is the example of what sort of detailed information is becoming available combining the advances of theory with the modern spectroscopic measurements.
It is also important that improvement in resolution and sensitivity of instruments (see Cortes et al. 2021, Henley et al. 2021) allows to map polarization of hundreds of molecular clouds with unprecedented accuracy, obtaining from dozens to hundreds of polarization vectors per cloud. The distribution of these vectors reflects POS distribution of magnetic field. The latter is determined by the interplay of the randomizing effects of turbulence and aligning effect of magnetic tension. The quantitative description of the two effects should relate the magnetic field strength with the intensity of turbulent motions. The latter can be obtained from the measurements of turbulent velocity fluctuations. The availability of detailed information of the distribution of both magnetic field fluctuations and turbulent velocity fluctuations calls for the development of quantitative techniques that can provide the detailed distribution of magnetic field strength. The latter being essential for the quantitative understanding of the role of magnetic field in molecular clouds, including the processes of star formation (see McKee & Ostriker 2007).
The well-known attempt to get the information using the logic above, unfortunately, provides not a distribution of Bstrength on the plane of sky, but a single estimate for a map. In addition, this estimate is known to be of insufficient accuracy. The technique that employs the observed fluctuations of polarization together with the fluctuations of the Dopplershifted velocities of gas was proposed by Davis (1951) and, independently, Chandrasekhar & Fermi (1953). Their technique termed Davis-Chandrasekhar-Fermi technique (henceforth DCF), assumes that magnetic field fluctuations and those of the velocity are directly related through the averaged Alfvén relation (Alfvén 1942): where δB is assumed to be perpendicular to the mean magnetic field and ρ is the mean density, which can be obtained through independent observations. It is natural to assume that velocity fluctuations with δv < V A , where V A is the Alfvén velocity, induce the deviations of the underlying field direction by an angle δφ ≈ δB/B mean . If we associate δv and δφ with, respectively, the dispersion of velocities σ v and angles σ φ available from observations, Eq.(1) can be used to express the mean magnetic field strength through these observables: where f DCF is an adjustable factor that traditionally assumed to be a constant. Another derivation that does not directly refer to Alfvén modes, but assumes that magnetic turbulence and kinetic energies are at equipartition is provided in Li et al. (2021). In fact, this assumption also corresponds to Eq. (1). For nearly incompressible runs this assumption corresponds to numerical studies of sub-Alfvénic turbulence (Haugen & Brandenburg 2004) at the injection scale where the DCF measurements are done. However, even at this situation, one should ask a question what fraction of energy is being in the fluctuations of magnetic field parallel to mean magnetic field and which part corresponds to the fluctuations perpendicular to the mean magnetic field. In the DCF, only the dispersion of the latter is being measured. According to Cho & Lazarian (2003), the energies in the aforementioned components are nearly equal for incompressible turbulence, which gives the factor f DCF ≈ 1/2 in Eq.
(2). When the compressibility gets important, more energy is being transferred into gas compression.
As a result, f DCF is, in fact, a function that depends on the actual properties of turbulence. This oversimplification results in DCF having a reputation of an inaccurate technique that can provide only order of magnitude estimates (see Li et al. 2021). Naturally, this decreases the scientific output of polarization studies. Nevertheless, attempts to improve the DCF accuracy have been numerous (e.g. Heitsch et al. 2001;Crutcher 2004;Houde 2004;Girart et al. 2006;Falceta-Gonçalves et al. 2008), but the foundations of the technique stayed the same. More recent notable suggestions include Hildebrand et al. (2009) ;Cho & Yoo (2016) and Skalidis & Tassis (2021), and we shall discuss further in our paper why we do not believe that these attempts can solve the challenges that the DCF faces.
The inability of the DCF to capitalize on modern high resolution polarization and spectroscopic maps stems both from its use of dispersion of angle φ, which implies averaging over the large areas of sky. The inaccuracy of the DCF arises mostly from the inadequate model of turbulence that is adopted in the derivation of the technique. Indeed, the DCF assumes that the perturbations of magnetic field direction arise from the collection of Alfvén waves moving along the mean magnetic field. This assumption, however, disregards the actual properties of MHD turbulence (see Beresnyak & Lazarian 2019 and ref. therein). Anisotropy of MHD turbulence is its fundamental property and the disregard of this property is a serious deficiency of the DCF approach. 1 The issue of angle γ between the magnetic field and the line of sight is usually not discussed within the DCF technique and the numerical testing are done for γ = π/2. However, it is easy to see that the DCF expressions should change as γ changes. In particular, repeating the arguments that brought Eq. (2), but taking into account that the magnetic field is inclined in respect to the line of sight, we can obtain a modified expression for the DCF in the form: where B ⊥ is the POS component of magnetic field. This is still an oversimplification that does not account for the actual properties of MHD turbulence.
To deal with the aforementioned deficiencies of the DCF, in this paper we propose an alternative technique of measuring magnetic field strength. Our approach employs a re-alistic description of velocity and magnetic field turbulent fluctuations based on the modern theory of MHD turbulence (Lazarian & Pogosyan 2012;Kandel et al. 2017a) and suggests measurements of velocity and angle variations on small scales using the structure functions.
Apart from addressing the issue of magnetic field for the evolution of molecular clouds and star formation, the present study that opens avenues for providing the distribution of magnetic field strength is important for understanding the properties of galactic foregrounds. Indeed, the recent progress of cosmological studies raise the importance of the detailed information of polarized galactic foreground to an unprecedented level. This is, first of all, is motivated by the search of the polarization arising from the enigmatic B-modes related to the gravitational waves in the early Universe (Kaminkowsky & Kovetz 2016). An additional motivation comes from the attempts to study emission of polarized atomic hydrogen lines (e.g. Zhang et.al 2020Zhang et.al , 2021 potentially at higher redshifts. The latter emission is not polarized at the source, but the effects of polarized foreground, nevertheless, seriously complicates the procedure of the foreground removal. Although the latter effect is related to the instrumental polarization, it is extremely challenging to get the signal when the unknown polarized foreground is a million times stronger than the signal. In §2 we discuss the properties of turbulent interstellar medium to which we seek the application of our technique. In §3 we introduce the observables that we deal with, as well as a new way of obtaining of magnetic field strength employing structure functions measured at the scales less than the turbulence injection scale. We demonstrate the advantages of using our differential measures compared to the global measures employed in the DCF. In §4 we introduce our approach of using the theory of MHD turbulence to accurately calculate the structure functions of velocity and positional angle φ. There we also explain how the anisotoropy of turbulence induces "anomalous" scaling of structure functions of φ with media magnetization. The projection of basic MHD modes are discussed in §5. The exact expressions for obtaining B-strength for pure Alfvénic and weakly compressible MHD turbulence are provided in §6. Obtaining magnetic field strength in the media similar to molecular clouds is quantified in §7. In §8 we provide the input of our compressible MHD simulations to guide us in finding which of the theoretically discussed cases are applicable to actual interstellar settings. We provide a set of simplified equations and practical recommendations for the analysis of the observational data in §9. In §10 we compare our technique to other existing techniques. We discuss the advantages and prospects for the our new technique in §11 and in §12 we summarize our work. Technical details are collected in Appendices. In particular, our list of notations can be found in Table 4, while our extensive set of numerical simulations used to test our analytical derivations is presented in Appendix D. Table 1. The typical Ms, MA, β values for interstellar media and molecular clouds , nH , δv for molecular clouds from Draine (2011), WNM/CNM from Kalberla et.al (2018), Yuen et.al 2021;Yuen et al. 2022c, β is from .

MODELING INTERSTELLAR TURBULENT MEDIA
Similar to the original DCF, we adopt a model of isothermal turbulent media. The assumption of isothermality is applicable to the molecular clouds as well as to the idealized phases of the ISM (see Draine 2006) 2 . This does not apply to the Unstable HI phase that constitutes a significant portion of the interstellar hydrogen Yuen et al. 2022c).
The turbulence is assumed to be injected at a large scale L inj and cascades to the small dissipation scale l d determined by viscosity, ion-neutral collisions or resistivity, depending on the media we study. The scales between L inj and l d correspond to the turbulence inertial range and the properties of turbulence at this range do not depend on the physics of the dissipation at l d . However, the properties of turbulence are affected by the properties of the driving that is given by the sonic Mach number M s = V L /V s and Alfvén Mach number M A = V L /V A , where V L is the turbulent injection velocity at the scale L, while V s and V A the sound and Alfvén velocities respectively. The response of the media to turbulence depends on the ratio of the gaseous to magnetic pressure β ∼ V 2 s /V 2 A . The warm diffuse medium has β > 1, while molecular clouds typically have β < 1.
For magnetic field strength studies in this paper we are interested in M A ≤ 1, as for M A 1 the turbulence is super-Alfvénic and marginally affected by magnetic field at scales larger than L inj M −3 A (see Appendix B). In this paper we focus on studying magnetic field strength in sub-Alfvénic and trans-Alfvénic turbulence. This automatically entails that turbulence in high β media is subsonic with M s < 1, meaning that it can be approximated by nearly incompressible turbulence. The compressible turbulence in our studies takes place in low β, i.e. β < 1, media.
The compressible turbulence in isothermal MHD turbulence has been extensively studied both theoretically and numerically (see Lithwick & Goldreich 2001, Cho & Lazar-2 For warm and cold phases the equation of state (EoS) is monoatomic-adiabatic (See, e.g. Kritsuk et.al 2017, Yuen et.al 2021Ho et al. 2021. For turbulence with adiabatic EoS one could express the mean sonic Mach number c 2 s ∝ ∂P ∂ρ. Simulation community have studied adiabatic turbulence for decades, e.g. Nolan et.al 2015. ian , Kowal & Lazarian 2010, Makwana & Yan 2020, also a monograph by Beresnyak & Lazarian 2019). The backbone of the theory of MHD turbulence is the model of Alfvénic turbulence by Goldreich & Sridhar (1995, henceforth GS95), that was formulated for M A = 1 and then generalized for smaller M A in Lazarian & Vishniac 1999, henceforth LV99). In the latter study it was shown that for scales larger than l tr = L inj M 2 A the Alfvénic turbulence is in "weak" regime, while for scales smaller than l tr the turbulence is in the "strong" regime. The "weak" and "strong" characterize the rate of non-linear interactions of turbulent motions and are not related to the amplitude of the turbulent perturbations. The properties of Alfvénic turbulence in weak and strong regimes are very different as we discuss in Appendix B The Alfvénic turbulence is very anisotropic and it imprints its properties on slow modes both in high-β and low-β media (GS95, Cho & Vishniac 2000, Lithwick & Goldreich 2001, while the fast modes have their own cascade that is isotropic both in high-β medium (GS95) and in low-β medium .
The anisotropic scaling of Alfvénic turbulence that we provide in Appendix B is defined in the frame of reference that is termed "local frame". This is the frame related to the local direction of magnetic field surrounding the Alfvénic perturbations at hand. The understanding of local frame is most natural in the eddy model of turbulence proposed in LV99. 3 Indeed, due to fast turbulent reconnection, in strong Alfvénic turbulence, there exist magnetic field eddies that mix up magnetic field perpendicular to the direction of magnetic field in the vicinity of the eddies. The local direction of magnetic field varies through the 3D volume and therefore the statistics of magnetic field is different for the observer.
These properties MHD turbulence in the frame of the observer were described in Lazarian & Pogosyan(2012, henceforth LP12). The latter and the subsequent studies (Kandel et al. 2016(Kandel et al. , 2017a provide the theoretical framework for describing our efforts in obtaining magnetic field strength from observations. The energy injection for different modes depends on the properties of the turbulent driving at L inj . The driving is expected to significantly affect the dispersions of turbulent velocities δv and δφ, but general theoretical considerations (see Beresnayk & Lazarian 2019) suggest that for sufficiently extended cascade, as it is the case of the cascade of the interstellar medium, the equipartition of energy between different modes tends to establish. The numerical simulations, even the largest ones (see Federath et.al 2021), have more limited inertial range and therefore the properties of the properties of turbulence driving may have imprint over the entire range of scales. Note, that the most robust in the cascade are Alfvén modes, while slow and fast modes are subject to collision-less damping (see Yan & Lazarian 2002, Brunitti & Lazarian 2007. Therefore, in the presence of the collisionless damping the fraction of Alfvénic waves and their importance is expected to increase with the decrease of the scale. In our study we consider an arbitrary admixture of Alfvén, fast and slow modes. Numerical issues present a serious problem that have not been properly considered in the framework of DCF studies. As MHD turbulence presents different regimes, the interpretation of numerical simulations and their relation to astrophysical reality may not be trivial. Indeed, even for the simplest case of incompressible MHD turbulence, the scaling of velocity and magnetic field fluctuations are different for weak, strong and super-Alfvénic turbulence. The relation between numerical dissipation scale the scales of the transfer from one regime to another must be carefully accounted for. Indeed, while in typical astrophysical conditions the dissipation scale of turbulence is much smaller than the L inj M 2 A for the transfer from weak to strong MHD turbulence (see Appendix ), for numerical studies for sufficiently small M A < 1 this is not true. Thus relating the numerical results with astrophysical measurements may be misleading.
In addition, our study includes an extensive use of compressible MHD turbulence simulations that are listed in Appendix B. These simulations are employed to test our theoretical predictions and the expressions that we obtain for the magnetic field strength.
Similar to DCF, we assume that the observed polarization represent the variations of projected magnetic field. The polarization can be of different origin, not only from dust, as we discuss in Appendix A. Our technique employs two types of observables, the polarization to represent the magnetic field direction, and the velocity centroids for measuring turbulent velocities.
Polarization as a tracer of the magnetic field direction: The magnetic field acting on dust grains, molecules and atoms induces polarization of the emitted radiation (see Appendix A). In polarization observations the Stokes parameters for, e.g., thermal dust emission are given by where for the case of the dust polarization n is the number density of dust grains, while θ and γ are the POS positional angle of the magnetic field and the inclination angle of magnetic field with respect to the line of sight, respectively. Simi-lar expressions are available for the polarization arising from atoms and molecules (see Appendix A). The statistics of polarization fluctuations in turbulent media is described in LP12. This is the study results of which our paper heavily relies upon. Velocity centroids : The velocity information on astrophysical turbulent volume is available from spectroscopic observations in the form of Position-Position-Velocity (PPV) cubes where intensity Em(X, v) is given as a function of the plane of sky (POS) coordinate X and the Doppler-shifted line-ofsight (LOS) velocity v. The first velocity moment of the intensity distribution is where we depending on the choice of the integration limits a, b one can get different measures. For instance, integrating over the entire spectral line width one gets a measure known as a velocity centroid. 4 One of the advantages of using velocity centroid is that such measures are not affected by thermal broadening.
Changing the integration from the velocity to real space, it is possible to show (see Lazarian & Esquivel 2003) that centroids are can be written as where L is the optical depth of the turbulent object under study. In what follows, we disregard that the fluctuations of ρ, assuming that their contribution can be significantly reduced in observation using the recipe in Yuen et.al (2021). The velocity centroids have been used extensively for studies of turbulence statistics (see Scalo & Elmegreen 2004). A detailed numerical study of their properties can be found in Esquivel & Lazarian (2005). The quantitative analytical study of the statistics of velocity fluctuations was initiated in Lazarian & Pogosyan (2000) and continued in the subsequent papers. In our study we will rely on the results obtained for the analytical description of statistics of velocity centroids obtained in (Kandel et al. 2017a), henceforth KLP17. 5 Within the DCF study the use of centroids to measure the dispersion of velocities was suggested in Cho & Yoo (2016). Though using similar information as DCF, this paper presents a new technique that focuses on short scale turbulent information and does not require assumptions about behaviour of the magnetic field and velocities at scales approaching injection scale.

Sampling Small Scale Turbulence
The technique that we introduce in this paper relies heavily on the statistical properties of MHD turbulence in the observer's frame. Before getting to the detailed calculations, we discuss the advantage of using measurements of magnetic and velocity fluctuations at small scales, rather than employing the dispersion of these quantities as it is done in the DCF. For the sake of simplicity, we do this first in the DCF framework, i.e. without accounting for the actual anisotropic properties of MHD turbulence.
Let us initially consider a toy problem assuming that the mean magnetic field B mean,⊥ is perpendicular to the line of sight, i.e. γ = π/2. Consider the variations of the observed magnetic field direction within a volume with size L measured along the line of sight and the turbulence injection scale L inj . The fluctuations of the magnetic field angle are where the integration is done along the line of sight and where, without losing generality, we assumed that δB ⊥ is a perturbation in the magnetic field component in the plane of the sky. Later, in this paper we will abandon the assumptions both about γ and toy model turbulence. Instead of following the DCF path and considering the dispersions of angles and velocities at the turbulence injection scale, we focus on the turbulence statistics at small scales that can be retrieved by using differential measures. In particular, if we are interested in the spatial variations of the observed magnetic field directions at the scale l, e.g. angle dispersions of magnetic field, those can be obtained using the secondorder structure functions (D) of the polarization angle φ 6 : where X is a two dimensional vector on the Plane of Sky (POS), ... X denotes the averaging over the position X. For practical applications, this means averaging for different X over the area l 2 . Substituting Eq. (7) in Eq. (8) one gets where the D 2D {B}(R) is structure function of the POS projected magnetic field. The latter is related to the 3D structure of the magnetic field D 3D {b}(r) as where as it is usual in turbulence studies, is assumed to be magnetic field in the turbulent volume is homogeneous, i.e. dependent only on the 3D lag between the points r (Monin & Yaglom 1976): with x being the 3D position vector and z 1 and z 2 denoting the pair of lines of sight along which the integration of the structure function of 3D fluctuating magnetic field b is performed.
For the sake of simplicity, let us discuss structure functions averaged over the positional angle, which will make these functions only dependent on the line of sight distance l separating the points. Observing that fluctuations of turbulent field are accumulated along the line of sight L in a random walk fashion one gets the structure function of the polarization angle φ where l is the separation between the lines of sight 7 . In statistical sense, the turbulence has the axial symmetry around the direction of the mean magnetic field (see discussion in Lazarian & Pogosyan 2012). However, for the sake of simplicity within this simplified treatment we do not consider the complications arising from the anisotropy in Eq. (12). In fact, the problem at hand has 3 length scales -separation on the sky l, integration/cloud depth L and the injection scale L inj , which is also the line-of-sight correlation length. In our discussion, it is assumed that l < L inj , thus the random walk takes place on with the step l and L inj does not enter the problem. If L l, the random walk results in D 2D {B}(l) ∝ D 3D {b}(l)lL, which makes the slope of the 2D structure function steeper by unity compared to the 3D structure function of magnetic field. Expressing D 2D {B}(l) through the structure function for the polarization angle (see Eq. 9), one obtains where D 2D {φ}(l) is available from observations. Let us now look at velocity fluctuations. Here we only consider the simplest case of pure Alfvén turbulence. For Alfvénic turbulence the fluctuations of velocity and magnetic field are symmetric in 3D. Therefore δv turb = δB ⊥ / √ 4πρ , where ρ is the medium density. Notice that for pure Alfvén mode the medium density is constant since Alfvén wave does not contribute to any density fluctuations (see Biskamp 2003). This means that the structure function of velocity is related to the structure function of the magnetic field in Eq. (11) as With observational spectral line data, one can measure the (constant density) 8 structure function of velocity centroids: which presents the proxy of the structure function of the velocities, averaged along the line of sight. In this procedure, the addition of velocity fluctuations is similar to summing up of magnetic perturbations δb turb that we dealt with earlier.
As the result, the summation process of the velocity fluctuations is a random walk process, i.e.
Combining Eqs. (13), (14) and (16), one gets the expression for the mean magnetic field: where, similarly to DCF, in our toy problem the factor f is constant. The actual calculations of this factor, based on the properties of anisotropic MHD turbulence, is done further in the paper and demonstrate it dependence on the Alfvén Mach number and the structure of turbulent modes. Reader should recognize that, while our Eq. (17) might look similar to the DCF relation in the RMS of the form, the sampling requirement is very different for these two methods. The radical difference between the DCF approach and our present approach is that we measure the structure functions at scales l significantly smaller than L inj which allows to use data sets covering regions smaller than L inj . Whereas DCF method requires sampling over l ∼ L inj . More importantly, with the use of structure functions, we can connect the magnetic field estimation to the detailed statistical theory of MHD turbulence.
We would like to stress that Eq. (17) is applicable to situations when the turbulence injection scale L inj is either larger or smaller than the line-of-sight depth L, as long as l L. The only requirement is that L should be the same for the calculations of D 1/2 2D {C}(l) and D 1/2 2D {φ}(l). This requirement is automatically fulfilled if the same spectral line data is used for obtaining both the polarization and Doppler broadening. This is the case for velocity gradients (see  or Ground State Alignment (GSA) (see Yan & Lazarian 2007). The case of dust-induced polarization requires more care to be sure that the polarization is collected from the same column of gas that contributes to the line emission. For instance, if the used line is 13CO, it is necessary to make sure that the column density of gas associated with CO emission is much larger than the column density of the HI along the same line of sight .
Finally, we can consider a limiting case of l L, when we get D 2D (l) = D 3D (l)L 2 for both the structure functions of magnetic field direction and centroids. This is a special case of studies when only a shallow surface area of the turbulent volume is being observed. This case can be realized in the presence of strong dust absorption as discussed in Kandel et al. (2018). However, it is easy to see that Eq. (17) is also applicable to this limiting case.
3.3. Role of instrumental beam in DCF and DMA Both DCF and differential DMA approaches compare the statistical measures of two distinct sky maps -one for polarization and the other related to LOS velocity. It is important for the technique that the polarization and the spectroscopic data sample the same regions. The differences can emerge due to the difference in the volumes sampled along the line of sight as well as due to the differences arising from differences of telescope beams. Within a modification of the DCF suggested in Hildebrand et al. (2009) , the careful study of the effects of the telescope beam is provided in Houde et al. (2009). Kandel et al. (2016) analyzed the isotropizing effect of a beam on the measured anisotropy of the velocity structure functions. The current paper deals mostly with the monopole part of the structure functions and therefore the beam-related suppression of higher multipoles reported in Kandel et al. (2016) is not so relevant.
The issues related to the beam size affect the practical data handling. In general, the two data sets used in DMA or DCF are obtained in different experiments, with different instruments, and as such the sky maps that serve as an input are convolved with different instrumental beams. This difference must be accounted for when numerical comparison between two measures is used to deduce the value of the magnetic field. This can be achieved by a detailed modelling of the beams through the formalism, though a more practical remedy is to deconvolve the beams from the datasets, and then smooth the datasets again with an equal synthetic beam to suppress the noise and artifacts induced by the instrumental beam deconvolution. The benefit of this procedure is the ability to control the properties of the final map smoothing, in particular, to suppress possible anisotropy of the instrumental beams that interferes with tracking the anisotropic turbulence effects, and to have a sufficiently effective noise reduction at small scales. Utilizing simple isotropic Gaussian smoothing is one example of such a well-behaved synthetic beam.
The downside is, of course, the reduced resolution of the final maps, since the synthetic smoothing window must be larger than either of instrumental beams. For differential measures, this limits the scale at which DMA can be applied and the appropriate optimal balance needs to be found when analyzing the data. Overall, controlling the smoothing of sky maps is more critical to DMA, as it is sensitive to small scales, than to DCF, as it operates primarily on large scales. In what follows, we assume that the input data has been preprocessed accordingly and smoothed to the same resolution both for polarization and centroid maps. While the DCF employs global dispersion of velocities and magnetic field directions defined on the turbulence injection scale, the DMA uses the local measurements by structure functions that can be applied at much smaller scales. The first notable advantage of this approach is that the value of magnetic field can be obtained over smaller areas of the sky. For instance, in Fig. 1 we show the coefficient f in DCF (Eq. 17) as a function of the sample size normalized by the injection scale L inj , using "run-2 (M A = 0.6)" (see Table 6) numerical simulation. It is clear that the DCF significantly overestimates the strength of magnetic field in the resolutions dependent manner, which is not straightforward to calibrate away. Only when the sampling size of the data is larger than the injection scale of turbulence L inj does the correction factor become constant. In the interstellar medium this scale is around 100 pc (see Chepurnov & Lazarian 2010) and exceeds the sizes of molecular clouds. 9 Therefore, DCF is bound to overestimate magnetic field if limited area over a molecular cloud is used for observations. Naturally, using DCF it is impossible to provide the detailed distribution of the magnetic field strength over the image of a molecular cloud.
In contrast, Fig. 1 testifies that the factor f in Eq. (17) does not change much. This allows using relatively small patches of data to obtain the magnetic field strength, in particular get the distribution of the projected magnetic field strength over the molecular cloud map. Thus, one can productively use the improved resolution of the polarization and spectroscopic maps that are getting available.
Apart from the regular change of f , we observe a significant statistical dispersion of its values within the classical DCF approach. We explain this effect in Appendix C. The dispersion of values when the Eq. (17) is employed is significantly smaller, which reflects the fact that the structure functions are more focused on small scale statistics and less subject to large scale variations.
The overestimation of the magnetic field strength in patches less than L inj ×L inj is not the only reason why DCF is inaccurate. As we discuss in later sections, the DCF has no inclusion on the nature of MHD turbulence, in particular, the anisotropic character embedded in some of the fundamental modes in MHD turbulence. Accounting for this differences is very difficult within the DCF, as the technique deals with large scale dispersions. At the same time, the properties of MHD turbulence change along the cascade, e.g. due to the transition from weak to strong turbulence that takes Figure 1. The variations of f as a function of sampling size for the DCF as compared to the variations of f in Eq. (17) for the simulation "huge-0" (See Table 6). The block sizes are given as a fraction of the turbulence injection scale Linj. The DMA method that we employed is from Eq.71 at R = Linj/198. Here, both DCF and DMA are applied assuming f = 1 2 .
place at L inj M 2 A (Lazarian & Vishniac 1999). The DMA focuses only on the small scale turbulent fluctuations by choosing the lag of the structure functions less than this scale, i.e. R < L inj M 2 A . This makes the DMA results that we derive further in this paper robust. 10

Realistic Inhomogeneity of Data on Large Scales
On large scales the structure of observed magnetic and velocity field is strongly affected by regular shear, gravity, outflows and other galactic processes. Therefore the dispersion δv and δφ that enter Eq. (17) can be significantly distorted, decreasing the accuracy of the DCF technique. In comparison, the DMA employs structure functions, that are not sensitive to the poorly controlled shifts of the mean background values. In fact, within our approach it is easy to remove any regular magnetic field and velocity distortions by using multi-point structure functions as we demonstrate further.
If necessary, the removal of regular contributions can be practically obtained with higher order structure functions which for the astrophysical context are discussed in Lazarian & Pogosyan (2008), Chepurnov & Lazarian (2009) (see also Falcon et al. 2007;Lazarian & Pogosyan 2008;Cho 2019).
For our approach there is no explicit dependence on how many points should we use in computing the structure func-10 At the first glance, if turbulence is uniform and homogeneous, Eq. (17) is reminiscent of Eq. (1) as for R → ∞ the structure functions are proportional to the squared dispersion. However, due to the aforementioned transition from weak to strong turbulence the expected scaling of f will change with R and it will be different for R < L inj M 2 A and much larger R > L inj for which the transition of the structure functions to the dispersion is justifiable. tions for the estimation of magnetic field strength, as the magnetic field strength enters the expression via the Alfvénic relation between the perturbations of magnetic field and velocity. Therefore with the multi-point structure functions we can use the equations in the main text ( §5) to determine the magnetic field strength.
The three & four point second order structure functions are defined as Cho (2019) has discussed how the multi-point structure functions can possibly remove the shear velocity field and recover the turbulent part. In a similar manner to 2-point structure function that cancels a constant contribution, 3-point and higher order structure functions cancel any linear and respectively higher order regular contributions to the field. The cost of using them is an increased noise contribution when applied to noisy data.
Here we test whether this is the case. Fig. 2 shows how the 2-point and 3-point angular averaged structure functions of the projected velocities of incompressible cube is affected by a large-scale coherent shear. Here we only perform a simple test assuming that the velocity experienced a constant randomly oriented shear: where A is some random vector representing the shear strength. We compare both the structure functions before and after shear modifications for the case of 2-and 3-point. We can see from Fig. 2 that the modifications of velocity has no effect on the amplitudes on the 3-point structure functions. However, for the 2-point structure function the introduction of the shear changes the structure function dramatically, both the slope and the shape of it. Therefore, the multi-point structure functions are indeed not altered by the presence of constant velocity shear, confirming the suggestions from Cho (2019). The dispersions of the observables, δv and δφ, that are employed by the DCF are one point statistics. Such measures are less demanding in terms of the number of measurements compared to the two point statistics of structure functions. The multi-point structure functions require even richer statistics. Thus the DCF could be applicable to the poorly sampled data, but at the expense of having low reliability results.
3.5. Common Problem of Eq. (17) and DCF: Failure to Account for Anisotropy and mode composition of MHD Turbulence The obvious limitation of Eq. (17) is that, similar to the DCF approach, it ignores both the anisotropic properties of MHD turbulence. This anisotropy changes with the media magnetization M A , the angle γ between the mean magnetic field and LOS, as well as the ratio of energy in Alfvén, slow and fast modes of turbulent motions. One effect of the anisotropy is that the contribution of the basic MHD modes to the variations of polarization angles and the variations of velocities changes with γ.  17)) (orange dots) as a function of γ obtained with the supersonic trans-Alfvénic numerical simulation ("huge-3" in Table.6). The actual value of POS magnetic field is given by a grey line. Fig. 3 demonstrates that Eq. (17) provides a better approximation to the actual B ⊥ compared to the DCF. However, we see that the fit is not perfect. This is one of the consequence of the anisotropic character of MHD turbulence that is missing in the derivation of Eq. (17).
We are not aware of any attempts to quantitatively account for the anisotropy of the MHD turbulence. When accounting for the anisotropy effect, one should keep in mind that the Alfvénic turbulence changes from weak to strong regime at a particular scale L inj M 2 A , with a significant changes of its properties. Within the DCF approach one deals with the properties of turbulence at large scales, i.e. mostly with the weak turbulence. On the contrary, using the DMA we are focused on the motions at the small scales, i.e. dealing with Alfvénic turbulence at the strong regime.
In addition, as we discuss in Appendix B, the realistic MHD turbulence is compressible and contains Alfvén, slow and fast modes, each having its own anisotropy. In §5 we outline our approach to dealing with the contributions of the modes and further in the paper we obtain the functional dependence of f that enters Eq. (17).

DMA: INTRODUCTION TO THE FORMALISM
The high resolution of polarization and spectroscopic data sets ushers a new era where more sophisticated analysis that provides high quality output is possible. The goal of DMA is to increase the accuracy of B-strength predictions by accounting for the actual properties of MHD turbulence. Introducing the differential measures in Eq. (17) is the first step in our derivation of the DMA formalism. Our next step is to provide a realistic description of turbulence at sufficiently small scales that are sampled by structure functions that describe magnetic field fluctuations as well as the fluctuations of turbulent velocities.
In Lazarian & Pogosyan (2012, henceforth LP12) we described the statistics of magnetic fluctuations arising from MHD turbulence. In the subsequent study by Kandel et al. (2017a, henceforth KLP17) the fluctuations of velocities have been described following the approach in LP12. These papers provide the basis for our detailed calculations.

Structure Function of Polarization Angles
For angle of polarization signal that traces the magnetic field direction (synchrotron, dust polarization, synthetic polarization from gradient maps) we can generally write cos(2φ) = Q/(pI) , sin(2φ) = U/(pI) where I, Q, U are the Stokes parameters (See §3.1) and pI is the polarized intensity For these quantities we can construct the structure function for polarization angle φ: where indexes 1 and 2 refer to two LOS separated by the 2D vector R = X 1 − X 2 on the sky. Note, that our measure differs from the measure introduced for the polarization angle in Houde et al. (2009) by the multiplier 2 in the cosine argument, and a factor 1/2 in front. The difference stems from the nature of the polarization direction that has a period of π rather than 2π. Naturally, in the small angle approximation, i.e. for |φ 1 − φ 2 | 1, this difference does not play a role. In this small angle case, our expression, as well as that of Houde et al. (2009), transfers to the "typical version" of the structure function of angle that coincides with D 2D {φ} given by Eq. (8) and that was also previously explored in Falceta-Gonçalves et al. (2008). However, D φ in Eq. (22) is a more general expression that is better defined from observations. D φ is also applicable to the case when angle fluctuations are large, e.g., when the Alfvénic Mach number is large.
In the case of dust polarization, which arises due to dust grain alignment by the magnetic field, we can model sky signal of dust emission as where B x , B y are the magnetic field components in x,y directions and z is along the line of sight. Eq. (24) is a combination of the local polar angle cos 2θ = (B 2 x +B 2 y ), and the angle between the magnetic field and LOS, sin 2 γ = (B 2 x + B 2 y )/(B 2 x + B 2 y + B 2 z ). Eq. (24) relates the properties of statistics of polarization directions that is available from observations and the underlying statistics of magnetic field, in particular, allowing to study the statistics of magnetic turbulence (see . Let us consider the coordinate system where the mean field lies in in (x, z) plane, B x = B sin γ, B y = 0, B z = B cos γ. The average in Eq. (22) can be computed rigorously in a series expansion in δB/B x . The leading contributions to Stokes parameters from Eq. (24) are Q ∼ sin 2 γ dz n dust U ∼ 2 sin 2 γ dz n dust δB y B which tells us that and that the leading term in the structure function Eq. (22) is provided by the U part If we assume that fluctuations of the dust density are of the same order as the fluctuations of the magnetic field δn dust n dust ∼ O δB B , then we find and where is the regularized projection of the 3D structure function D yy (r) = (B y (r 1 ) − B y (r 2 )) 2 for the magnetic field ycomponent that is orthogonal to both LOS and the direction of the mean field. Remarkably, in this leading order perturbations of the dust density (assumed to be of the same order as perturbations in the magnetic field) do not affect the result. Our numerical tests show that when the field is perpendicular to the line of sight, Eq. (28) is accurate to under 1% for M A = 0.15 (b21 in Table 6). At higher M A = 0.66 (b15 in Table 6) accuracy varies from 1% when structure functions are measured small lag to ∼ 10% at R ≈ L inj , which points to better accuracy when measuring angle differences at short separation instead of the angle variance as in DCF. For M A = 1.1 (b42 in Table 6), which is near the limit of approximations in Eq. (28), the accuracy drops to ∼ 40%. Note, that for magnetic strength determination, the uncertainties are roughly half of the quoted ones, since D φ enters in a square root.
Our leading order approximation, that starts with Eq. (25), features a constant angle of the magnetic field sin γ to LOS, whereas its local value sin 2 γ = (B 2 x + B 2 y )/(B 2 x + B 2 y + B 2 z ) varies along the line-of-sight as reflected in Eq. (24). Further we will account for this next-order effect within the model of "magnetic field wandering" along the LOS.

Effect of Turbulence Anisotropy on the Angle Fluctuations
According to Eq. (28), the most important component that determines the polarization angle fluctuation is the ycomponent of magnetic field. The y-component magnetic field fluctuations contain both spatial power distribution and also anisotropic information of magnetic field fluctuations in the sky. Following description of magnetic field fluctuations in LP12, Appendix E contains general derivation of D yy (R). However, for our context here we focus on the case when the y-axis fluctuations are described by Alfvén mode power spectra only A = F = E, which, in particular, is the case for strong turbulence in high-β regime. We shall also average our structure functions over all directions of R, i.e we will measure the monopole of the structure functions D 0 (R) = dφ R D(R = (R, φ R )). Then, from Eq. (E19) : where γ is the line of sight angle cos γ =B ·ẑ. Let us compare Eq.(30) to the total power in POS components given, as seen from Eqs. (E19, E20), by which differs by the absence of cos 2 φ k factor. The presence of the cos 2 φ k is the most important reason why the f parameter should contain a natural dependence on M A . If the spectrum of turbulence was isotropic, we would have D yy (R) = 1 2 D B ⊥ (R). However, for a typical anisotropic MHD power distribution, the power is concentrated in the wave modes orthogonal to the mean magnetic field, with suppression of power for modes aligned with it. Thus, E(K, sin γ cos φ K ) is peaked at cos φ K = 0. In particular, for solenoidal magnetic field that its perturbations are orthogonal to the wave vector in 3D which results into the projection of cos 2 φ K factor in Eq. (30) for the ycomponent. This projection will be zero when the modes power E(K, sin γ cos φ K ) attains its maximum. As a result, the contribution of structure function fluctuations D yy (R) due to solendoial projection is suppressed, so as the polarization angle perturbations, in comparison to expectation from the general level of perturbations in δB ⊥ . Mathematically, we can express the polarization angle fluctuations as a function of D yy (R)/ D B ⊥ (R): From our argument above, we find that the factor Dyy(R) D B ⊥ (R) 1 2 , and is smaller, when the MHD power spectrum is more anisotropic. Moreover, its magnitude is turbulent model dependent, but as we shall see further A for many cases of sub-Alfvénic turbulence.

Multipole Expansion of Polarization Structure Function
Let us now return to a general exposition of the angle structure function. We follow the description of magnetic fluctuations in LP12 and Appendix E. In Appendix E.3, we obtain Eq. (E25) for D yy (R) for sub-Alfvénic turbulence in the strong regime at sufficient short scales R < L trans = L inj M 2 A (See Eq.(B1)) and where the LOS depth L R as well. Using it in Eq. (28) gives the following expression for the coefficients of the multipole expansion D φ (R, φ R ) = n D φ n (R)e inφ R of the polarization angle structure function where I n (R) are scaling functions defined in Eq. (E26). Coefficients E 2D p are POS angular harmonic decomposition of the projected power spectrum of the magnetic field which depend on the angle γ of 3D orientation of the mean field relative to LOS and the Alfvén Mach number M A . This additional dependence on M A reflects the change of the anisotropy of MHD turbulence with the change of media magnetization and, as we just argued for in the previous section, can significantly change the DMA approach compared to the naive treatment given by Eq. (17). The geometrical functions G (I,A,F ) (γ) depend on the mode structure of the turbulence, with several particular cases discussed in § 5. From KLP16, KLP17 the higher order multipoles encode the information of γ since γ and M A enters to the multipole contributions in different power. However for our purpose in this paper, we will only consider the case when n = 0 in the RMS of B-field estimations.

Multipole Expansion of Velocity Centroid Structure Functions
The next step in DMA is to evaluate the structure function of velocity centroids in the nominator of Eq. (17). The statistics of centroids was discussed in Kandel et al. (2017a). There, for the sake of theoretical convenience the definition of centroids was modified compared with the standard one given by Eq. (5). In particular, the density was set to be constant in the computation of velocity centroids. This approach is now possible to realize in observational data analysis, as the recently developed Velocity Decomposition Algorithm (VDA, Yuen et.al 2021) allows one to significantly mitigate the effects of the density fluctuations. Thus, the theoretical model in Kandel et al. (2017a) coupled with VDA becomes applicable to the observational data that is frequently strongly affected by density inhomogeneities. Therefore, we proceed with the same constant-density assumption in velocity centroids assuming that we can always remove the density effects from observation using the VDA.
Velocity centroids structure function has a very similar behaviour to the polarization angle one. The multiple moments of the structure function of centroids KLP17, normalized by the mean column intensity of the gas along the line of sight Em = Lρ, where is the emissivity coefficient, can be written in the form similar to Eq. (33) : where we have expressed the amplitude of velocity fluctua-tionsÂ p of KLP17 via the variance δv 2 . In Appendix E.4 we show that the scaling function I n (R) here is the same as in Eq. (33) and give the geometrical functions W s for the velocity centroid structure functions. Velocity centroids structure functions are less affected by the anisotropy of the power spectrum than the angle structure functions. None of geometrical weights W, as Eq. (E28) attests, have a suppression of the projected modes orthogonal to the magnetic field via cos 2 φ K factor that G (I,A,F ) functions contained. Thus, no additional small parameter ∝ M A arises in the evaluation, e.g., of D v 0 .

The Ratio of the Multipole Coefficients of Magnetic Field Angle and Velocity Centroid Structure Functions
Eq (33) and Eq. (34) are purposefully written in the form to highlight the common and distinct parts of the angle and centroid structure functions. From these two equations the ratio of the centroids and angle structure function multipole coefficients is clearly If we consider Alfvén and slow modes for β 1, or consider Alfvén and fast modes for β 1 , the velocities and magnetic field values are connected through Alfvénic relation given by Eq. (1) the variances are Then we can obtain for all even n: Within our assumption of the similarity of velocity and magnetic field scaling, the factor f n in Eq.(38) does not depend on the POS lag R. This is true for the velocity and magnetic field fluctuations for Alfvénic turbulence, turbulence of slow modes and the admixture of Alfvén and slow modes. We deferred the study of fast modes that have different scaling , Kowal & Lazarian 2009) to §7. In general, for an arbitrary mixture of 3 modes with different spectra, A, S, F, one can write: (39) where we have factorized the Alfvén mode power spectrum and have used the relation Eq. (36) for the ratio of the magnitudes of δB and v fluctuations in the Alfvén mode.
We note that the evaluation of the magnetic field strength depends via f n on M A as well as the orientation of the magnetic field with respect to the LOS. These dependencies arises due to difference of anisotropy of line-of-sight velocities and perpendicular magnetic field fluctuations in MHD turbulence. We explore the properties of f n further in § 5.
It is important to stress once more that in our studies we focus on the small scale asymptotic behavior of MHD turbulence. The description of MHD turbulence is much less certain in the vicinity of the injection scale and can be significantly affected by the turbulent energy injection mechanisms. We also assume that we sample turbulence at scales less than the transition from the weak to strong turbulence A . This condition is especially crucial to check with for simulations with very low M A . Therefore our calculations are not applicable to larger scales, e.g. beyond L inj , and the transition from the DMA expressions to the DCF ones is not straightforward.

PROJECTION OF MHD MODES
Before considering particular admixtures of Alfvén, slow and fast modes of MHD turbulence 11 , in this section we first discuss the observed properties common to all feasible magnetic fluctuations. More specifically, we will discuss first the basis of the frame vectors which decide the directions of the modes fluctuations. We will discuss how both the velocity and magnetic field angle fluctuations are accumulated along the line of sight. In particular, we discuss an important but easily overlooked suppression of the projected magnetic angle fluctuations when the depth of integration is significantly larger than the correlation scale of the magnetic field fluctuations.

Orthogonal Modes of Magnetic and Velocity Vector Fields in the Presence of a Globally Preferred Direction
In the system with the preferred directionλ 0 set, in our case, by the mean magnetic field direction, the fluctuations of any solenoidal vector field such as δb can be decomposed, for each Fourier mode k, into two orthogonal components; of A-type (Alfvén) along ζ A = ( k ×λ 0 )/| k ×λ 0 |, and of F -type along ζ F = k × k ×λ 0 /| k ×λ 0 |. In the context of MHD turbulence, Alfvénic magnetic field perturbations are of A-type; while both slow and fast modes induce the same F -type component in the magnetic field perturbations, though they differ in the power distribution among different wave vectors k.
A general vector field, such as the vector of turbulent velocity, also has the third component, orthogonal to both A and F ones, that is along ζ P = k, which we'll call P (potential) one (See Fig.4). This is where slow and fast turbulent modes structurally differ by generating velocities orthogonal to each other at each Fourier mode.
As we discussed in Section 2 for the description of MHD turbulence it is convenient to distinguish the case of mag-11 In a number of sources the combination of slow and fast modes is called compressible to distinguish them from an incompressible Alfvén mode. This is misleading as both fast and slow modes have compressible and incompressible components. The actual distinction can be done on the basis of the Helmholtz decomposition (see Kowal et al. 2010). Thus, we always put "compressible" in "..." while talking about fast + slow modes.
netically dominated media with magnetic pressure P mag larger than the gas pressure P gas , i.e. the case of β = P gas /P mag < 1, and the gas pressure dominated media, i.e. with β > 1. The properties of basic modes of MHD turbulence are different in low and high β cases (see Cho & Lazarian 2003).
In the limiting case of high β plasma, where β is the ratio of the gas and magnetic pressures, slow mode velocities are of F -type, while fast mode ones are purely potential P ones. In the opposite limit of low-β plasma, this pair is rotated, each being a mixture of F and P motions, with slow mode velocities aligned along the magnetic field, while fast modes perpendicular to the magnetic field. Alfvén mode velocities remain to be of A-type. . An illustration of the frame vector system in the global frame of reference. The vectorsζ A,F,P correspond to the A,F,Pmodes. Readers should be careful that A-mode do correspond to the Alfvén contributions, but F-mode is not necessarily the fast mode. In particular, for the displacement vector the fast mode actually lies in betweenζ F andζ P , while for magnetic field fluctuations both slow and fast modes are alongζ F . See Cho & Lazarian (2003); Lazarian & Pogosyan (2012) and Appendix E.
In the global reference frame given by the mean magnetic field direction λ 0 , one can always formally decompose a vector field into A, F , P components. In the physical sense the Alfvén, fast and slow modes in the turbulent cascade are defined with respect to the local direction of the magnetic field λ (see Lazarian & Vishniac 1999, Cho & Vishniac 2000, while large scale fluctuations in the magnetic field direction are viewed as fluctuation of the preferred directionλ. 12 Such locally defined modes are then found to be mixed between AF P types when viewed in the global reference frame. This effect is discussed in more details in Yuen et al. (2022b).

Projection Effects due to LOS Signal Accumulation
Applying the DMA to observations one deals with the signal integrated along the line of sight. Let us summarize how this projection affects the contribution of individual basic MHD modes on the quantities of interest -fluctuations of the direction of the polarization angle that in the leading order is given by δB y , and the LOS component of velocity in velocity centroids v z . As the DMA samples fluctuations at small scales R < L inj , the effects of the averaging along the line of sight L can be more important for the technique than for the DCF technique that deals with the dispersion of fluctuations at scales larger than L inj . Nevertheless, when the techniques are tested with the numerical simulations obtained with periodic boundary conditions, the integration is formally over infinite scale for both the DCF and the DMA.
By choosing lag R in the DMA analysis we identify the eddies of the scale L corr = R that we sample. At a minimum level, integrating over the LOS depth L beyond the correlation length L corr we expect a decay of the relative magnetic field fluctuations δB/B due to random walk as (L corr /L) 1/2 . However, our signals depend only on specific components δB y and v z . Moreover, the structure of global modes in Fourier space leads to an additional suppression that depends on the mode, angle of the mean magnetic field to LOS and the depth of integration L.
Namely, one can see that in the limiting case when the integration length L is larger than the POS scale R, the LOS projection results in averaging out of all Fourier modes with k z > L −1 . This asymptotically leaves only the Fourier components with wavevectors lying in POS, k = ( K, 0). Consider now a setting with the mean magnetic field perpendicular to the line of sightλ = (Λ, 0) and Alfvén mode fluctuations. Since A-mode excites only the field component simultaneously perpendicular to the wave vector and the mean field, ζ A is along LOS, i.e. after the projection of the A-mode in this configuration no δB y is remaining. Thus, in the leading order angle fluctuations of the projected magnetic field in POS are absent. As another example, F-mode lies in a plane spanned by the wave vector and the mean field, so since the projection restricts the wavevector to POS, the F-mode will not produce any δB y when the mean field is along LOS.
This illustrates the geometrical suppression of the contributions from A and F modes when the mean magnetic field is either perpendicular or parallel to the LOS. Table 2 summarizes the asymptotic limit of the projection effects. We Table 2. Contributions from A, F and P parts to fluctuations of velocity and magnetic field in the limiting case of infinite LOS integration depth and perpendicular or parallel orientation of the mean magnetic field. . The numerical factors are given by the fraction of the modes' total power that are seen after projection.
Mode δBy, γ = π/2 δBy, γ = 0 Vz, γ = π/2 Vz, γ = 0 also note that the P mode does not contribute to DMA in this regime, since magnetic field is solenoidal, and the velocity contribution is averaged out. In Table 2, the numerical factors are given by the fraction of the mode total projected power, that is carried by the component that affects the observable signal.
Let us investigate how the projection changes the variances of the angles and velocities as in DCF, when the signal accumulates over the depth L as it varies relative to the injection scale L inj . We perform our analysis numerically, on synthetic representation of the turbulent magnetic field and velocity motion.
In Fig. 5, we show the dependence of the standard deviation for the POS magnetic field angle and velocity centroid with the LOS integration depth in units of the injection scale. For L > L inj centroid rms fluctuations (divided by L) decrease as (L inj /L) 1/2 in accordance with the random walk expectation, that is shown in both upper and lower panels with dashed lines. However angle fluctuations decrease faster both for strong turbulence and even more so for purely Alfvén turbulence, due to projection effects. The last point on the graphs corresponds to full projection over the periodic simulation box which is 10L inj , which is formally equivalent to the infinite integration range. In this limit, the fluctuations of magnetic field angle due to Alfvén mode formally vanish.
The effect of the suppression of Alfvénic mode contribution has been missed, as far as we know, by the researchers studying the DCF. In fact, the suppression ∼ N −1 has a simple explanation. The fluctuations of Alfvén modes as viewed at γ = π/2 can be modelled by a simple sine wave sin(kL max ). This approximation is only valid for Alfvén modes in this geometry since if γ = π/2 the Alfvén mode has additional contributions that we discussed in the Appendix. In this geometry it is easy to see that the cumulative RMS fluctuation accumulated along the line of sight is k los =0 dk sin(kL max ) ∼ cos(kL max )/L max ∼ 1/N instead of N −1/2 .

PURE ALFVÉN CASE AND INCOMPRESSIBLE TURBULENCE
Alfvén and incompressible MHD turbulence are valuable idealization to be considered. In both cases the Alfvénic rela-tion between magnetic field and velocity fluctuations is valid. The case of Alfvén mode is useful to discuss for both low and high β cases as the Alfvén modes are known to be dominant in inducing variations in 3D magnetic field directions, as was shown in LV99. Our discussion will focus on geometric factors G for angle fluctuations and W for velocity centroids, that enter determination of the magnetic field strength given by Eq. (37). Here we will formally show that in 2D projection Alfvén mode contribution to angle fluctuations may be suppressed when the mean magnetic field is nearly perpendicular to the line of sight. However the Alfvén mode remains after projection if the magnetic field has significant LOS component.
In the global reference frame we will be speaking about A-mode as Alfvén mode. The contribution of A-mode to the structure function of the projected magnetic field angles with L R is represented by the geometrical weight in the structure function of δB y component, reproduced here from Eq. (E22) Contribution of the Alfvén modes to the (unnormalized) 13 velocity centroids is governed, in turn, by (KLP17) We note that G A n vanishes at γ = π/2, thus global A-type mode is not reflected in the fluctuations of the projected magnetic field angles when the mean magnetic field is perpendicular to LOS. At this orientation A-mode is projected out, or highly suppressed if we account for finite LOS depth L.
Similarly, W A n vanishes in the opposite limit when λ is along the line of sight, γ = 0. However, the case of γ = 0 is less of an issue since B ⊥ = 0 in this case and it is problematic to determine the magnetic field strength in this configuration anyway.
When using both Eq (40) and Eq. (41), n is to be taken even. For most angles both G A n and W A n decrease substantially for |n| > 2 (see Fig. 6 as well as LP12, KLP16, KLP17a). In particular, for π/8 < γ < 3π/8 the leading RMS fluctuations in both geometrical functions are the monopole (n = 0) and the quadrupole (n = 2):  . Geometrical factors G A n (blue) and W A n (orange) for n = 0 (solid), n = 2 (dashed) and n = 4 (dotted) as functions of the mean magnetic field orientation γ.
As γ → π/2 all orders of G A n become present but are equally suppressed, while W A n retains only the monopole

Power Spectrum in Strong Turbulent Regime
For sub-Alfvénic M A < 1 strong turbulent regime, the power spectrum can be concisely described by Eq. (E39) which gives the following ratio of the projected 2D multipoles in terms of modified Bessel functions I p of the variable with the latter expression being a reasonable fit even for M A⊥ exceeding unity.
At small M A⊥ we may utilize the expansion which shows highly anisotropic power distribution approaching δ-function behaviour that concentrates power in the modes, orthogonal to the projection of the mean field. This limit gives rise to all multipoles in the power spectrum being of the same magnitude with alternating signs. As M A⊥ increases, the ever lower multipoles become subdominant. For M A⊥ > 0.3 one can limit consideration to the monopole and quadrupole. The following expression is a reasonable fit for the quadrupole even for M A⊥ exceeding unity which we draw both Eq.44 and Eq.46 on the left of Fig.7. On the right of Fig. 7 shows this ratio of projected power multipoles as the function of M A for three values of γ as given by the exact form of Eq. (46).

Transforming Alfvén Modes from Local to Global Frame of Reference
Changing our focus to physical turbulent Alfvén mode defined locally, we follow LP12 model that accounting for the large scale magnetic field direction can be considered as wandering of λ direction along LOS, level of which is determined by the Alfvén Mach number M A . LP12 suggests that the result of wandering (See also Yuen et al. 2022b)can be described by replacing the anisotropic geometrical function by the weighted combination with the isotropic term as and the power spectrum E( k · λ) by the average E( k · λ) according to Eq. (E37).
We adopt a simple model for the weights 14 together with the angle distribution of the power in Alfvén spectrum given by Eq. (E38). This model reflects, on one hand, that in isotropic, M A → ∞, limit W I → 1 /2 since Alfvén mode has one degree of freedom whereas general isotropic solenoidal field would have two, and on the other hand, that in describes the level of isotropization arising due to change in the mean magnetic field direction along the LOS. The α = 1 case can be realized when there is an idealized Alfvén mode on the background of a fixed mean field for small M A . While α = 0 corresponds to completely isotropic tensor structure of the turbulence when M A > 1. The resulting geometrical structure of the local Alfvénic turbulent perturbations is that of a mix of A and F modes, for our purpose here represented as a mix of A and the isotropic I = A + F combination of the correlation tensors.
14 This model is more accurate for Alfvén modes than what was suggested in LP12. We use this opportunity to note an inconsistency in LP12 where W I as used there in Eq. (71) is twice the one introduced in Eq. (48). See also Yuen et al. (2022b).

Asymptotic Expressions for f 0 Factor
Let us start with explicit form of Eq. (38) for the monopole coefficient f 0 when the geometrical functions are given by Eq. (47) that takes into account the wandering of a local direction of the magnetic field along the LOS. Using α defined in Eq. (49) as the measure of wandering level, one obtains the following expression Remarkably, for small M A⊥ , when the power spectrum multipoles are given by approximation Eq. (45), one can perform the summations in Eq. (50) completely, by noting that Substitution of these formulae into Eq. (50) gives a compact expression for f 0 factor in Alfvénic turbulence This shows that magnetic field wandering eliminates degeneracies due to projections that arise due to idealized structure of formal global A-mode when γ = π/2. However, some suppression of Alfvénic perturbation in line-of-sight projection for perpendicular field remains to be a real effect. From Eq. (52) at γ = π/2 we find (1 + M 2 A ) 2 (53) which shows that at perpendicular orientation of the mean field to LOS in the strongly sub-Alfvénic regime one has f 2 0 ∼ M 4 A dependence that comes equally from a high anisotropy of the power spectrum and a small amount of mean field wandering, each effect contributing one M 2 A factor. Though this configuration is rather special, the case often appears in numerical simulations and one should be aware of additional suppression of angle fluctuations in this case.
Although the power spectrum treatment in Eq. (52) has been truncated to M 2 A order, for γ π/2 expanding Eq. (52) to M 4 A gives the correct asymptotic behaviour for f 0 since the omitted corrections are proportional to the product M 2 A cos γ 2 and small. Eq. (54) demonstrates that the de- In the limiting case of sufficiently small M A √ 2 cos γ, Eq. (52) gives which corresponds to neglected field wandering. Thus we can interpret the transition from f 0 ∼ M A to f 0 ∼ M 2 A as a transition to the regime when field wandering becomes important.
The dependence of f 0 as a function of M A calculated numerically in our model is presented in Fig. 8 for several orientations of the magnetic field γ. The figure demonstrates the linear f 0 ∝ M A behaviour at low M A for all orientations of the mean field not strictly perpendicular to LOS. The formula Eq. (52) works well at least for M A 0.3 for all γ's. At higher M A 0.3, the exact dependence becomes steeper than Eq. (52), which is expected since Eq. (52) uses power spectrum expansion truncated to the quadratic M 2 A⊥ order. At nearly perpendicular configurations (e.g., γ = 2π/5) we observe an almost quadratic f 0 ∝ M 2 A raise in 0.3 M A 0.6 range, before saturating at M A > 1. As γ → π/2, the transition point from linear to quadratic dependence slides to ever smaller M A following cos γ in accordance to Eq. (54), revealing for γ = π/2 the quadratic f 0 ∝ M 2 A scaling across the whole range M A 0.6, as accurately described by Eq. (53).
For configuration of the mean field more aligned with LOS, f 0 factor is close to unity for 0.1 M A 1, though still scaling linearly f 0 ∝ M A at small M A 's. Worth noting is that Mach number dependence may become non-monotonic as γ = π/8 case in Fig. 8 attests to. To test the predictions from Fig.8, we perform the same calculation in numerical simulations by extracting the Alfvén contributions using Cho & Lazarian (2003). We see from M A < 0.5 is due to the strong turbulence transition is more and more restrictive in our numerical simulations. (See Appendix B). Clearly, in general, the determination of magnetic field strength depends on our knowledge of the angle between the mean magnetic field and the line of sight γ. The statistical determination of this angle is possible (see Paper II), but in this paper it acts as a parameter. Note, that this presents an uncertainty both for DMA and the traditional DCF techniques. In our case, the functional dependence on this parameter is explicitly given.

High-β Nearly Incompressible MHD Turbulence
If the sonic speed c s is larger than Alfvén velocity V A , the media is high-β, as β ∼ c 2 s /V 2 A . Such media is widely presented in interstellar environments.
Gas in HII regions, warm interstellar and circumstellar gases can be approximated as weakly compressible high-β gas (See Tab.1). In this case fast MHD modes are very similar to the sound waves and their contribution vanishes due to the integration along the line of sight provided that R is significantly smaller than the integration scale L. 15 In such weakly compressible media, the slow and Alfvén modes have very similar power spectra, being two "polarizations" of a solenoidal wave and they have equal contributions to the observed fluctuations. The tensor of "Slow + Alfvén" fluctuations has simply the isotropic solenoidal form. In this case W I n−p = δ np while G I n−p = 1 2 δ np + 1 4 (δ n−2,p + δ n+2,p ). As a result we find Focusing on the monopole f 2 0 = 1 2 (1 + E 2D 2 / E 2D 0 ) as the simplest case to measure observationally, we notice that Eq. (56) is equivalent to Alfvén case for α = 0, since the monopole case have no effects to the anisotropic structure of the modes. In this case, however, it is an exact result that only the monopole and the quadrupole of the power distribution add contributions to f 2 0 . Since their ratio in high-β plasma remains to be given by Eq. (46), we obtain for sufficiently small M A,⊥ :f  (46) is not truncated to the leading order in M A⊥ . This result of high-β weakly compressible turbulence should be contrasted with that of the purely Alfvénic turbulence, especially in near perpendicular configuration at small M A . The origin of small value of f 0 in sub-Alfvénic regime in both cases is the suppression in projection of δB y that is responsible for angle variations of the observed magnetic field lines whereas velocity centroids fluctuations are not suppressed and dominated by the full power of Alfvén modes. However, while high-β result comes only from the anisotropy of power distribution, in purely Alfvénic turbulence we have an additional geometrical suppression due to specific anisotropic structure of Alfvén perturbations. Thus, ignoring projection effects in the patches of the sky where the mean field is perpendicular to LOS may lead to an overestimation of the magnetic field strength by a significant factor for low M A .
In the areas of the sky where the mean field direction is close to LOS, in weakly compressible turbulence Alfvén mode dominates POS magnetic field fluctuation, while slow mode is responsible for the LOS velocity centroids. Both contribute similar power and f 2 0 ∼ 1 2 . We remind however that in the patches of the sky where γ is below some critical angle γ crit , DMA (and DCF) measurements can not be interpreted as a measure of the mean magnetic field. 6.6. High-β Turbulence with Isotropic Driving One may notice that in low M A regime the high-β model with equipartition of energy between Alfvén and slow modes does not lead to isotropic distribution of velocities, due to highly anistropic power distribution. Indeed, as M A → 0 we get v 2 = v 2 ⊥ , where parallel component with respect to the magnetic field is dominated by the slow mode, while perpendicular components are given by the Alfvén mode. But this means that v 2 z = 2 v 2 x = 2 v 2 y , if z is along the magnetic field. As anisotropy of the spectrum decreases with M A increase, velocity distribution becomes more and more isotropic.
If, however, the turbulence is driven isotropically at all the energy in the slow mode will be less than that in the Alfén one, with the M A dependent fraction v 2 S / v 2 A , changing from one half as M A → 0 to unity as M A > 1. We can approximate this transition as v 2 , consistent with our simulations. In Figure 11 we observe that the decrease of power in slow modes relative to Alfvén ones, results, for small M A , in an increase of f 0 at small γ by a factor up to √ 2, and in a decrease of f 0 at γ ≈ π/2 by up to the same √ 2 (this limiting value corresponds to the case when slow mode contain half the energy of Alfvén ones). Dependence of the fraction of the slow modes on M A makes f 0 to depend separately on M A and γ rather than just on M A⊥ combination. Nevertheless, the fit which works very well for γ ≈ π/2, continues to work for all M A⊥ < 1. In this expression f 2 0 ∼ M 2 A⊥ is a rigorous asymptotic at small M A⊥ , including the prefactor, that stems from Eqs. (46,52). The last term is a fit for the transition to M A⊥ ∼ 1.
7. LOW-β MEDIUM CASE Low-β turbulence presents a complex case where all three turbulent modes play part in the estimation of the magnetic field. We shall focus on the limit β → 0, in which slow modes are the magnetic compression propagating along magnetic field lines and therefore only marginally contribute to 1.0 f 0 Figure 11. The factor f0 in high-β turbulence with isotropic driving (blue) as compared with the case of equal energy in Alfvén and slow modes (grey, same as in Figure 10), shown as a function of 3D angle of the mean magnetic field γ for, from top to bottom, MA = 1, 0.8, 0.6, 0.4, 0.2, 0.1.
the changes of the magnetic field direction. Also, numerical simulations show that the contribution of the fast mode is subdominant to the velocity of the medium. In this situation, Slow modes are expected to have similar spectrum as Alfvén modes with the same scaling I S 0 (R) = I A 0 (R) and angular dependence S 2D p = A 2D p and the relative power in two modes is constant. On the other hand fast modes may have a different scaling I F 0 (R) from Alfvén modes. We adopt as before the anisotropic model power spectrum for Alfvén mode given in Eq. (E39) and assume that the fast mode power is distributed isotropically, F p = F 0 δ 0p . Eq. (59) can then be significantly simplified for n = 0. Noting that W S n−p = δ np cos 2 γ (Eq. E29) and finding from Eq. (E22) that For practical applications the difference between the scaling I F 0 (R)/I A 0 (R) of fast and Alfvén modes (see Cho & Lazarian 2003, Kowal & Lazarian 2010) is weak enough not to be of secondary importance for the range of R involved in the DMA study. Thus, in the zeroth approximation, the difference in Alfvén and fast mode scaling may be disregarded.
For small M A⊥ < 1, the multipole summation for Alfvén modes that led to Eq. (52) can be used again. We also need to additionally take into account the field wandering effect on the fast modes, which is accomplished by replacing G F 0 → W I + W L G F 0 , while the wandering effect on slow modes is accounted for by G S 0 → W I + W L G S 0 and using the same averaged spectrum, as for the Alfvén mode. The result is where 2D monopole of the projected Alfvén angular spectrum A 2D 0 (M A , γ) is given by Eq. (E40) to be The last expression is the leading term of expansion at small and the dependence of f 0 on M A is lost due to the fact that the variations of angle φ arise in this case from isotropic fast modes.
The above requirement on the relation between M A and can be somewhat relaxed for γ → π/2. The factor f 0 in this limit exhibit the dependence only on the ratio of energies in fast and Alfvén modes: which corresponds to the notion that slow modes in low β regime are similar to sound waves that cannot contribute neither to the variations of φ nor to the variations of δv for γ = π/2. For γ → 0 and M A < 1, Assuming also that M A⊥ 1, Eq. (61) transforms to which demonstrate that in this regime the variations of φ arises from anisotropic Alfvénic modes and therefore f 0 gets ∼ M A⊥ . However, the requirement of M A⊥ 1 is very restrictive in this case.
In the opposite limit of negligible fast mode contribution, i.e. Let us first look at the balance of contributions to angle fluctuations. In Fig. 12 (left panel) we show the ratio of Alfvén to Fast modes to the projected angle structure function at equal power in both modes (i.e δB 2 F ≈ δB 2 A ). The curves are, thus, variance normalized, and reflect only the geometrical structure of the modes. They are to be multiplied by the ratio of powers when it is known at the scale of study. Fig. 12 focuses on the dependence of contributions as a function of orientation γ and for different Alfvén Mach numbers. Fast modes depend on M A due to field wandering only, while Alfvén modes depend on M A primary due to anisotropic power distribution, while the effect of field wandering dilutes this effect.
We observe that when the magnetic field is predominantly aligned with POS, the fast mode plays an important role in POS angle fluctuations of the field. The reason behind is based on the earlier discovered suppression of Alfvén contribution for such orientations due to the projection effects, as discussed in § 5.2. The effect is stronger when M A is smaller. At the same time, fast mode magnetic perturbations in POS projection, given by G F = G F (see Eq. (E22)), do not vanish at γ = π/2. Having isotropic power spectrum, low-β fast modes will perturb equally the magnetic field components parallel and perpendicular to the mean field, thus perpendicular POS perturbations will not be suppressed. Hence, for γ ≈ π/2 the observed variations of magnetic field and, correspondingly, the angle fluctuations will be dominated by the fast modes over the Alfvén as well as slow one. When γ is less than π/2, Alfvén modes begin contributing more to the projected angle and become dominant. Transition angle depends on M A and how much power in fast modes there is.
The balance of Alfvén and slow contributions to the velocity centroids sets the denominator in Eq. (59). With similar spectra between these modes, the balance is set by the scale independent power ratio v 2 S v 2 A . In the right panel of Fig. 12 we study M A and γ dependence for equal power between two modes. In this case we find that Alfvén mode dominates velocity fluctuations at γ > π/4 while slow mode takes over at γ < π/4, in full accordance with simple model Eq. (61).
The relative importance of slow and fast modes depends on driving. Our analysis of simulations in §8 indicate a relatively low fraction of energy in fast modes at the level of 20% relative to Alfvén ones and roughly equal contribution of slow and Alfvén modes. In Fig. 13 we show the modelled f 0 taking δB 2 F / δB 2 A = 0.3M

1/2
A and v 2 S / v 2 A = 0.6. This is a special case and we see that as M A increases the contributions of the modes to angles and velocities compensate each other to make f 0 nearly independent from γ. The most significant variations with γ we observe for low M A . Our analysis above provides crucial insight on how the statistics of velocity and angle φ behave according to the statistical theory of MHD turbulence. Our study in §7 demonstrates that the use of DMA expressions, and also the magnetic field estimation method via the DCF-like formalism, requires accurate knowledge of the media provided that magnetized turbulence is present and dominant. For instance, the required initial knowledge is different in the case of low and high β media. For instance, only slow and Alfvén modes are important for high β case. This makes the study rather robust for separations R that are smaller than LM 2 A . In particular, for high beta case we expect that the energy of Alfvén and slow modes are similar in order of magnitude or even equipartitioned, which is expected true for sufficiently extended inertial range of incompressible turbulence. In this case, the dependence of f 0 to γ is much simpler to analyze according to theory, as we can see in §6 in detail.
The low β case is more challenging for practical studies. This, however, a very important case as it includes studies of turbulence in molecular clouds. The relative proportion of the modes can change and their contribution to the fluctuations of angle and velocity changes with angle γ. This is reflected in §7. The main idea of our analysis in §7 is that, when β < 1, the velocity fluctuations mostly depend on the Alfvén and slow modes, while that of polarization angles depend on the Alfvén and fast modes. Therefore both the composition of turbulence and also the line of sight angle γ is crucial in estimating the magnetic field strength.
Both angle γ and composition of turbulence in terms of modes are parameters that can potentially be obtained from observations (see , Hu et al. 2022. However, in what follows we will numerically explore the output of numerical simulations to evaluate to what extend the DMA can be applied if these aforementioned parameters are not determined prior to the application of the technique. For our suite of low-β models, the distribution of power in slow, E S , and fast, E F , modes relative to the power E A in Alfvén modes is shown in Fig. 14. We find that the slow mode fraction is relatively constant at E S /E A ≈ 0.6 level, while the fraction of fast modes shows a weak dependence on  it was found that when the fast modes are generated through the conversion of energy from Alfvén modes starting from initial perturbations that are purely Alfvénic, the relative energy in the generated fast modes is proportional to M A . In our simulations, the setting is different, the fast modes are generated by the solenoidal driving and the dependence of the energy of fast modes on M A is found to be more shallow, i.e. ∼ M 1/2 A as seen in Fig. 14. 8.1.2. π/2 Case In Fig. 15 we show the numerically evaluated correction factor f 0 for a set of low-β sub and trans-Alfvénic simulations for γ = π/2. The latter choice is motivated by the fact, that all the numerical studies of the DCF and alternative models (see e.g. §10) are done for this choice of angle.
We observe that, while a significant scatter is present, the supersonic models which have sonic Mach number M s > 1 follow To explore further the nature of the M A dependence for f 0 , in Fig. 16 we plot separately the (square root of the monopole of) the angle structure function √ D φ and the velocity structure function √ D v , measured in the units of Alfvén velocity V A = B/ √ 4πρ. The ratio of these two functions give f 0 in Fig. 15.
The velocity structure function shows linear dependence on M A , √ D v ≈ 0.1M A which just reflects overall proportionality of the level of fluctuations in Alfvén modes to M A as we discussed in §6 and Appendix B. The angle structure function, however, has a steeper dependence, that is fit quite well, at least for supersonic case low-β turbulence, by To understand what the extra M

1/4
A dependence reflects and to compare these measurements with our theoretical picture in the previous section, we decompose the MHD turbulence into modes following the algorithm in Cho & Lazarian (2003) and look at the contributions of individual modes in the result.
In Fig. 17 we see that at π/2 orientation, the contribution of the Alfvén modes to the angle fluctuations have been suppressed while fast modes do not contribute to velocity fluctuations, as expected in theoretical model (See §3, Tab.2). We do not show a separate contribution of slow modes, as on its own it is subdominant in amplitude at π/2 to both angle and velocity measurements (See §7).
Velocity structure function is completely determined by the Alfvén mode, there is virtually no difference between right panels of Fig. 17 and Fig. 16. For angle fluctuations fast modes is the dominant source. Non-zero contribution of Alfvén modes that grows with a M A is the result of field wandering effect, and is absent if fluctuations of sin γ along LOS are not included in the projection. Thus, our theoretical model, resulting at γ = π/2 in Eq. (64), is on a firm basis.
Eq. (64) predicts that any observed dependence of f 0 on M A is due to the M A dependence of fast-to-Alfvén magnetic perturbation ratio δB 2 F / δB 2 A . In our simulations (see Fig. 14), we found that the relative energy in the fast modes in low β sub-Alfvénic plasma roughly scales as 0.3M A . Thus the total amplitude of fluctuations from fast modes, which in low β medium are compressions of magnetic field, is ∝ M A × M 1/4 A , which exactly explains the best fit for angle structure function from the dominant fast modes in Fig. 17. Note that Eq. (64) predicts the magnitude f 0 (M A = 1) ≈ 0.49, in agreement with measurements.
We notice, however, that taking the ratio At the same time, this ratio is rather noisy, and we have found that the better fit to f 0 is given by Eq. (67). We feel that determining the exact scaling of f 0 (M A ) for γ = π/2 may require more detailed numerical studies. Here we suggest that for low beta and γ = π/2 the dependence of f 0 is ∼ M α A , where α is between 1/2 and 1/4.

π/4 case
The contribution of Alfvén modes to projected magnetic angle fluctuations increases with γ and therefore we consider a case of intermediate angles, i.e. γ = π/4. The intermediate angles are usually not discussed within the DCF study where the composition of the turbulence in terms of modes is ignored together with the effects of Alfvén and slow mode anisotropy.
For γ = π/4 we know both theoretically (see §3) and numerically (Fig. 21) that the strong projection suppression of Alfvén mode is not present and as the result field wandering does not qualitatively change the results. At γ = π/4 Eq. (61) gives where κ(α) = 1+ For velocity fluctuations, in Fig. 18 we observe, as expected, that √ D v scales as M A and is determined almost in equal fraction by Alfvén and slow modes, as was predicted in Fig. 12, given that we have a similar energy content of two modes seen in Fig. 14.
On the other hand, for polarization angles, we see from Fig. 19 that the Alfvén mode contribution is larger than that of the fast modes. It is expected at large M A , as Eq. (68) Figure 14. The energy fraction of slow to Alfvén mode (left) and that of fast mode relative to the Alfvén mode (right) in our simulations with β < 1 and M>1 (Table 6)  predicts. However, according to Eq. (68), fast modes ere expected to take over at small M 2 A κ δB 2 F / δB 2 A which, if we take κ ≈ 0.2 and the fast mode fraction to be 0.3M

1/2
A , resolves to M A 0.15. This is not observed in D φ in Fig. 19 which shows Alfvén modes prevailing over the fast ones even at smaller M A . This originates in Alfvén contribution falling as, approximately, Correspondingly, instead of f 0 ∼ M A expected for the Alfvénic component, in Fig. 20 we observe f 0 ∼ M 1/2 A for supersonic simulations and f 0 ∼ const for subsonic ones.
We believe that here we may face a restriction related to the resolution of our numerical simulations. Let us remind the reader that the origin of f 0 ∝ M A result for Alfvénic mode comes from the anisotropy of the power spectrum, namely Eq. (45), that reflected the properties of the strong turbulence, see Eq. (E39). As we discuss in Appendix B, in sub-Alfvénic regime the Alfvén modes exhibit two different regimes of turbulence, the weak turbulence from the injection scale down to the scale LM 2 A and the strong one for scales less then the transition scale LM 2 A . All the calculations of the DMA are performed for the strong turbulence regime. However, for sufficiently small M A this regime cannot develop if LM 2 A < d diss , where d diss is the scale of numerical dissipation. Therefore for small M A we deal with the pure weak Alfvénic turbulence which has different anisotropies from our model of strong turbulence that we employ in DMA. For our current testing, we can only perform analysis of the strong Alfvénic regime up to M A ≥ 0.2. 16 Future studies with higher numerical resolution should test our prediction. For fast modes, there exist only one regime of turbulence and therefore where fast modes are dominant, we see the scaling of D φ that corresponds to the change of the fast mode energy with M A , unaffected by resolution. At large lags the structure functions saturate at the twice the dispersion level. In this limit we transfer to the case similar to that the DCF deals with. Our analytical study does not cover this limiting case, as the DMA is limited to sufficiently small scales for which the MHD turbulence is in the strong regime. To have a more complete picture of the evolution of φ fluctuations below we explore this case numerically. Fig. 21 illustrates the dispersion of angles σ φ arising from different MHD turbulence modes in low β turbulence. For simplicity we did not include the contributions of slow modes since in theory they have infinitesimal contribution to σ φ when β < 1 (See §7) As expected, for γ = 90 (see left panel) the fast mode dominates over the Alfvén modes . In this scenario, the total σ φ scales approximately as ≈ 0.1 M A . It is important to note that the magnitude of polarization angle fluctuations projected onto the sky is significantly less than M A , i.e. σ φ M A . This adds to the accuracy of approximations that were made in our theoretical approach. The contribution of Alfvén modes to σ φ fits the M 3 A scaling, be-coming more shallow as M A → 1, in accordance with our description of "field wandering" (see § 6.3-6.4). It would be interesting to compare the dependencies of σ φ and D φ (left of Fig. 17) to M A as the former is simply the "large scale" result of the latter. From the left of Fig. 17 it is indicative that the change of the fast modes energy with M A may be different at small and large scales, although it is difficult with the existing noisy data to definitively distinguish the change of power from M

5/4
A to M 1 A . The contribution of Alfvén modes becomes much more important when γ switches from π/2 to π/4 both for σ φ and D φ , which is expected as in §7. For the case of σ φ to M A , which we can see from the right panel of Fig.21, we see that Alfvén mode takes over fast modes in terms of the amplitude. Furthermore, we notice that the dispersion σ φ arising from Alfvénic modes scaling is best described as M

5/4
A , though the slope M

3/2
A as was D φ is also possible (see right panel of Fig. 19 ). This is again differs from M 2 A , which is the scaling for D φ expected theoretically for Alfvén modes (see §3). The difference between the two scaling laws is significant and we interpret the difference as resulting from the transfer of the Alfvénic turbulence to the weak regime at large scales. We do not have a theoretical prediction for this regime, but σ φ ∼ M A seems plausible. Furthermore, we note that the fast modes contribution to σ φ remains ∝ M A , but at π/4 its relative role decreases compared to the case of π/2.
In terms of practical separation into modes, a better accuracy can be achieved for the separation into Alfvén and the mixture of slow and fast modes. The latter is frequently called "compressible", although this does not really reflect the nature of the modes, which have both compressible and incompressible parts. Fig. 22 shows that the contribution of these "compressible" modes to σ φ ∝ M A . Notice that the scaling of σ φ to M A does not change as γ changes from 80 o to 40 o , however there is a slight amplitude change as γ changes. In the case of low β, this constancy of the slope is related to the isotropy of magnetic perturbations induced by fast modes. Therefore, at the range of fast mode domination, the σ φ part arising from "compressible" modes is independent of γ. For this case if the velocity field has typical scaling σ v , √ D v ∝ M A , then it is very probable to arrive to f 0 ∼ const, recovering the classical DCF expression given by Eq. (3). Notice also that in the case of Alfvénic modes dominance to σ φ it is troublesome to ignore the dependence of f DCF on γ. Our discussion above shows that the error can be significant for sufficiently small γ. We shall discuss in depth ( §9) on how to utilize the extra dependence of M A will alter the formula of magnetic field estimation that we used before (Eq. 71). However, before we proceed with the exact magnetic field estimation, let us discuss what we expect for the statistics of f 0 .
Due to the complex dependencies ofD and D φ , the statistics of f 0 obviously depends on the relative contribution of modes, the Alfvénic Mach number M A , the compressibility β and the line of sight angle γ. From what we discussed in §8.2 the contributions of Alfvén mode does not act the same in velocity and polarization angle as we discussed in Tab.2. In particular, for the low β case, the relative contribution of Alfvén mode increases as the angle γ decreases from π/2 to zero. Therefore, initially at γ = π/2 we expect the fluctuations of φ to be dominated by fast modes. The transition to Alfvénic-dominated variations of φ happens at some intermediate angle γ, the value of which depends on the ratio of the energy of modes that is affected by the properties of driving. Therefore, for strongly compressible driving the transitional angle γ is expected to be smaller than for the incompressible driving that we employ in this paper. In addition, for very small M A fast modes can dominate for a wider range of angles due to the dependence of the contribution the Alfvén modes to φ on M A (see Eq. 50).
In the situations that fast modes dominate both variations of velocity and φ, no dependence of f 0 on M A is expected. However, what we observed in Fig. 21 that while for γ = π/2 fast modes are dominant for the φ variations, the variations of velocity can be dominated by Alfvénic component. Thus we get the observed dependence of f 0 on M A .
The main takeaway of the current work is that the DCF scaling factor f 0 is a function of M A , γ and the relative composition of modes. The non-trivial dependence of f 0 and also its lack of systematic testing over the decades for all three items contribute to the reasons on why not until our work do we observe that the f 0 factor is not a constant. More importantly, since the contributions of different mode plays dramatically different role in terms of magnetic field estimations, whether a particular mode is important depends on the plasma compressibility β. For instance, in the case of molecular clouds with β 1, the scaling factor f 0 depends on the mode composition in the form of δv 2 A /δv 2 S and δB 2 F /δB 2 A . This discovery is non-trivial and has not been discussed before in literature.
Notice that the mode ratio is also a function of the line of sight angle γ. As we explained in §6 the structure functions of velocities D v (R) are expected to scale with M 2 A for all MHD modes, irrespectively of their anisotropy. At the same time, in §6. we showed that the turbulence anisotropy significantly modifies the behavior of the angle structure functions D φ (R) compared to the naive expectations employed in DCF. The effect is most prominent for Alfvén modes and it is also present for slow modes. For equal energy in Alfvén and slow modes, the fluctuations of angles induced by slow modes are important in high β regime and subdominant in low β regime.
To utilize our results from §6-8, we start with Eq. (37): where the structure functions of velocities D v (R) and angle D φ (R) should be calculated at sufficiently small scale R < L inj M 2 A in order to insure that the MHD turbulence is in the strong regime (see Eq. B1). Here the factor f 0 = f 0 (M s , M A , γ, mode fraction) should be computed according to the analysis that we provided in §6, 7. To proceed, we recognize that f 0 is generally in the form of f 0 ∝ CM n A for some C = C(γ) and some constant n. The idea of reformulation is to express M A back via the mean magnetic field substitute it into f 0 dependence and resolve now implicit expressions for B via the ratio of structure functions wrt B.
In this scenario, if f 0 ∝ CM n A⊥ , M A,⊥ = M A / sin γ 17 , one can write the square of the perpendicular magnetic field strength as where readers should be careful that V L is the injection velocity and D v = D v /V 2 L . It follows from Eq. (71) that our empirical study in §8 suggests that for low β medium f 0 ≈ 0.8M 1/2 A for supersonic 17 Here we are assuming that the γ value is not too small so that our expressions in §6, 7 do not fall into a regime that M A,⊥ > 1 despite M A,tot < 1.   Table 6. turbulence and γ = π/2 and therefore Naturally, the dependence for γ = π/2 does not persist for other γ as the contribution from different modes to φ and velocity fluctuations changes with γ. One may argue that the contribution of the modes may get to the universal equipartition between them for sufficiently extended cascade. However, this has not been proven and, even if true, it is not clear that for a given observation this is satisfied. At least near the injection scale, the energy in the modes depends on the properties of turbulence driving. For our choice of incompressible driving, we may see that Eq. (72) is applicable also to γ = π/4, provided that the corresponding prefactor changes from 0.4 to 0.8.
In Table. 3 we collect the limiting cases that we discussed in §6 and §7, and in Fig. 23 we provide the numerical comparison on our formulae in Table 3 as compared to both the "naive" formulation assuming f 0 = const in Eq. 17 and also the traditional DCF calculation. For reader's reference, if the curve is closer to 1, that means the estimation of magnetic field from the specific formula is closer to the actual magnetic field strength .Notice that if the DCF result is not plotted, that means it is significantly deviated from the ranges that we are plotting (i.e. only 0.1 < B DM A /B true < 10) is plotted. We can see that our complex formalism (blue curves on each panel in Fig. 23) performs significantly better than both the "naive" formalism and also the DCF method.

Effect of Density Fluctuations
In the present work so far we disregarded the effect of density fluctuations, assuming that we can use the mean density instead. The contribution of density fluctuations can be significantly reduced using the VDA approach (Yuen et.al 2021, Appendix H), which is particularly efficient for subsonic turbulence for both velocity and polarization angle fluctuations. Below, however, we do not apply any of these techniques and deal with the data as is given by synthetic observations where both the fluctuations of φ and velocity centroids calculated in the presence of turbulent density. Fig.24 shows a comparison of f 0 obtained with synthetic observations with constant ρ and with actual turbulent ρ. The overall behavior of f 0 is not changed. Withing the DMA the density over the averaging block should be employed for the calculation of the B-strength.

Obtaining V L
The velocity at the injection scale enters the our expressions for B-strength (see Eq. 71, Table.3) One way to obtain this velocity is via measuring of the total line width of the Doppler-shifted lines.
This approach is not applicable, however, in the presence of galactic shear. In this situation, V L can be obtained from the turbulence sonic Mach number M s = V L /c s , where c s is the sound velocity. There are a number of ways of obtaining M s from the statistical analysis of data (see Chepurnov & Lazarian 2009, Burkhart et al. 2014 ). Therefore, V L can be obtained as where c s is known if the temperature of the gas is known, or by kernel estimation method proposed in Yuen & Lazarian (2020a). More discussion of how to obtain V L will be provided in Lazarian, Yuen & Pogosyan (2022, Paper II).

Sub-Alfvénic and Super-Alfvénic Turbulence
Within this paper we consider the application of DMA to sub-Alfvénic turbulence. One might wonder how do we extend our analysis to super-Alfvénic turbulence? This is a totally valid question that we plan to address in an upcoming publication but first we should discuss how super-Alfvénic turbulence is different from that of sub-Alfvénic turbulence from the theoretical point of view.
The main concept that we would like to repeatedly emphasize in this paper about the turbulence system is that, the statistics of full 3D turbulence is dramatically different from that of the projected observables. In 3D turbulence, it is well known from theory (See, e.g. , or a recent review from Beresnyak & Lazarian 2019) that the seemingly super-Alfvénic turbulence will behave like a sub-Alfvénic turbulence as long as the scale of consideration is smaller than the cut-off scale L inj M −3 A . If we measure the statistics of a super-Alfvénic turbulence and a sub-Alfvénic turbulence in sufficiently small scale, they will both exhibit the Alfvénic behavior as predicted in theory and also verified in numerical simulations. This also suggests a main way of distinguishing the turbulence behavior as the sub-Alfvénic turbulence at large scale only reverts to the weak cascade, but that for super-Alfvénic turbulence reverts to hydrodynamic cascade. The differences of the statistics of these cascades can be measured via structure functions at large scale. The understanding of the correlation scale L inj M −3 A could also be put in this way: in 3D statistics, for scales smaller than A there will be Alfvénic cascade, and larger than that it should be hydrodynamic L inj M −3 A . Notice that M A is much more easier to obtain than the actual magnetic field strength, e.g. by measuring the anisotropy of turbulence (see Esquivel & Lazarian 2005, Lazarian & Pogosyan 2012, Kandel et al. 2016, 2017a. A powerful way of finding M A is based on exploring the distribution  or curvature of velocity gradients (Yuen & Lazarian 2020b).
However, could we interpolate the scaling argument about into 2D observables? If it could be, we can definitely apply our theory to observations as long as sufficient filtering is applied to both centroid and polarization map. Fig.25 shows how the structure functions of velocity and polarization angle behave in two selective simulations with sub-Alfvénic and super-Alfvénic turbulence, respectively. We can see from the LHS of Fig.25 that, in the case of sub-sonic turbulence, the structure functions for both velocity and polarization angle behave very similarly, as predicted in theory carrying the same R 1+m factor, and saturates at the exact same transitional scale of l tr = L inj M 2 A which defines the transition from weak to strong turbulence. However, for the case of super-Alfvénic turbulence, we can see from Fig. 25 that the slope of the polarization angle structure function is significantly flatter than that of the velocity structure function until R < L inj M −3 A .The difference between the two scales allows to estimate M A at least to the level to distinguish the case of sub-Alfvénic turbulence to which the DMA in its present form is applicable. Furthermore, in Fig. 26 it shows that the effective f is changing at the L inj M −3 A scale. An additional indication of the super-Alfvénic turbulence is that the slope of the spectral functions for φ gets more shallow compared to the slope velocity. Note, that the two slopes are similar for the sub-Alfvénic turbulence. 9.5. Obtaining the mode fraction from independent measurement Figure 23. The ratio between the estimated magnetic field over the true magnetic field B estimated /Btrue for different cases list in Table. 3 . The optimal method would have this ratio = 1 The figures right now are in ascending order as in  Some of our limiting case formulae (Table 3) depend on the mode fractions δB F /δB A or δv A /δv S . Observers might question how do we obtain the mode composition from observations? While we show in Fig.14 that these two ratios are simple functions of M A , it would be synergistic if we obtain these two parameters independently, instead of a derived variable of M A , so as to minimize the uncertainty of the magnetic field strength estimation.
While there is no widely accepted method in obtaining the actual ratios between three modes, recent progresses have been made via the newly developed Synchrotron Polarization Analysis method (SPA, Zhang et al. 2020b;Yuen et al. 2022b) allowing the retrieval of the "compressible"-to-Alfvén ratio from direct observations of polarization maps. Notice that the "compressible modes" that SPA method is discussing about is the F-type tensor we discuss in the current work, i.e. including the contributions of slow and fast mode. Therefore the SPA technique provides a upper bound for the slow or fast to Alfvén mode ratio. This value is exceptionally important when we are in the cases where the β is either very high (warm, diffuse ISM, see Ho et al. 2021 for the reason why colder ISM, including both the unstable and the cold phases, are not ideal candidates) or very low (molecular clouds). In either of these cases, we have the corresponding predictions outlined in the previous sections, e.g. high β case in §6 and low β case in §7 and 8. The additional knowledge of mode ratio from SPA will enable us to obtain a more accurate estimation of magnetic field through on our results in Table 3 Cho & Yoo (2016), henceforth CY16, considered the case of the injection scale of turbulence L inj being less than the extend of the line of sight L within the emitting turbulent volume. The authors noted that while the total line width employed as δv in Eq. (2) represents the full dispersion of velocity in the volume, the variation δφ in Eq. (2) is a result of a random walk of φ from one turbulent injection scale to another. As a result, the DCF technique overestimates the magnetic field unless L is less or equal to L inj . To remedy this problem CY16 proposed to measure δv using velocity centroids.
To explain the problem addressed by CY16, consider a setting with mean magnetic field being along x-direction in the plane of the sky and the Alfvénic fluctuations δB is along y-direction. If the 3D magnetic field is b reg , it is adds up linearly along the line of sight and therefore the observed B x is L b reg dx ≈ b reg L. On the contrary, the fluctuating magnetic field b turb with correlation scale L inj is added up in the random walk fashion with δB y providing L b turb dx ≈ b turb L inj L. 18 As a result an additional factor enters the δB y /B x ratio, namely, the observed fluctuation gets reduced by a factor ≈ L inj /L, which corresponds to a random walk suppression.

CY16 Expression
To account for this factor, CY16 considered the ratio of the line of sight velocity and the centroid velocity. The latter is given by Eq. (5), while the former is the usual δv los arising from the velocity dispersion at the scale L inj . The velocity measured by centroids is, on the contrary δC = L δv los dx/L ≈ δv los L inj /L. As a result, if δv los is substituted by the dispersion of Velocity Centroid δC that can be used both for the case of L ∼ L inj and L L inj . In other words, the expression with constant f CY 16 that can be obtained from numerical simulations. Eq.(74) has a wider range of applications than the original DCF expression as they show in the series of numerical works (Cho & Yoo 2016;Yoon & Cho 2019;Cho 2019). In particular, the magnetic field strength computed based on Eq.(74) would not depend on L inj /L ratio. However, as we see from the applications of the DCF technique, it is very rare that the objects that observers are considering are having the size comparable to the injection  scale. For instances, common molecular clouds like Taurus, Perseus have their size at the order of 1 − 10pc, but the turbulent injection scale is believed to be at the order of 100pc (Chepurnov et.al 2015, Yuen et al. 2022a). If we only consider the signals from the cloud along the line of sight, the ratio L inj /L is unlikely to be smaller than 1. In other words, the suppression of the dispersion δφ due to random walk arising from the addition of contributions from independent turbulent regions along the line of sight is rather unlikely. Naturally, as all the modifications of the techniques, Cho & Yoo (2016) assumes that the equipartition of the magnetic energy and the kinetic energy, which is known to be violated for the sub-Alfvénic turbulence (Haugen & Brandenburg 2004).

Comparison with DMA
The similarity of the DMA with CY16 is that both techniques employ centroids of velocities. The difference is that the DMA employs the structure functions rather than dispersion. This provides the differences between our approaches at the level of our Eq. (17), which is the starting equation for the derivation of the DMA formalism. As a result, our measurements can be performed locally providing the distribution of magnetic field strength in a molecular cloud rather than a single value. In addition, structure functions are less affected by the large scale non-turbulent distortions.
It is very important that the DMA treatment did not stop at Eq. (17) but, on the basis of MHD turbulence theory, provides analytical predictions for the functional dependence of f 0 on angle between the line of sight and mean magnetic field γ, as well as the Alfvén Mach number M A . In contrast, CY16 considers the constant f 0 from the perspective of number of turbulent eddies along the line of sight, irrespective to the properties and orientations of the eddies according to the theory of MHD turbulence theory.
An important question arises of whether the expressions for f 0 that we obtained within the DMA can be used with the CY16 approach. Indeed, the structure functions of the velocities and the angle saturate at the values equal to twice of dispersion of the corresponding observables. However, this is not true, as the properties of sub-Alfvénic MHD turbulence change along the cascade from strong regime at small scales to weak regime at larger scales. The DMA deals with the robust scaling of small scale strong turbulence, while at larger scales the turbulence gets into the weak regime that we repeated described in the main text and also Appendix B. The weak regime of turbulence is beyond the existing DMA description and therefore there is no way of direct utilizing the analytical expressions obtained in this study for improving the predictive abilities of CY16.

Hildebrand et al. (2009): Model of uncorrelated turbulence
The Angle Dispersion Function (ADF) method was introduced in Hildebrand et al. (2009) andHoude et al. (2009). Its similarity with the DMA is related to employing structure functions of φ. However, while the DMA employs the structure functions at scales R < L inj , Hildebrand et al. (2009), as we will argue, is only applicable when R > L inj .
10.2.1. Model of Turbulence Adopted Hildebrand et al. (2009) employs structure functions of the polarization angle, which is similar to the DMA, but uses the linewidth to measure the velocity dispersion. 19 In this work of the authors replace δθ → D 1/2 φ and investigate its properties with the traditional DCF assumption that MHD turbulence is isotropic (cf §B), which is one of the differences from the DMA.
Within these assumptions Hildebrand et al. (2009) propose that the variations of the angle φ can be modelled by the Taylor expansions of structure functions of polarization angles: where b, m are fitting factors. There b describes the contribution from turbulence, which is applicable in the situation when the structure of turbulence is not important and only large scale dispersion of φ matters. The factor m is related to the large-scale variations of mean magnetic field, which are assumed to be small making the Taylor expansion possible. Within the adopted model, it is then shown that the turbulent to regular magnetic field strength ratio is: The rest of the approach is based on the DCF method, with B 2 t / B 2 associated with δφ 2 which finally gives the estimate As Hildebrand et al. (2009) employs the velocity dispersion measured using the linewidth, the technique has to deal with the issue of estimating the number of independent fluctuations along the line of sight ∼ L/L inj . This is elaborated further in Houde et al. (2009) that estimate L/L inj in molecular clouds. In CY16 and the DMA this issue is solved by construction by using velocity centroids structure function.

Applicability of the Model
We shall point that the model of turbulence that amounts to a constant in the structure function (more accurately, a constant at all lags except R = 0) is the model of negligible correlation length in the turbulent fluctuations. Basically, turbulent fluctuations are treated as a white noise.
However, as we discussed in Appendix B, the turbulent cascade is correlated on scales up to the energy injection scale L inj (see Monin & Yaglom (1975)). For R < L inj where we sample turbulent motions over the inertial range, we expect 2D SF tur ∝ R m+1 with fractional power index, e.g., for Kolmogorov turbulence m = 2/3, that extends for multiple orders of magnitudes, as shown in the form of electron power spectra (Armstrong et.al 1995), Hα data (Chepurnov & Lazarian 2010) and also more recently with successive combinations of HI and CO data (Yuen et al. 2022a). In general, the expectations from turbulent fluctuations is that m + 1 < 2. Thus Hildebrand et al. (2009) assumption of the correlation scale of turbulence being negligible makes their model not applicable at R < L inj to turbulence that we expect to encounter in the ISM and molecular clouds (see Table 1). This is a significant difference with the DMA in which the turbulent scaling at R < L inj is taken into account.
It can be asked whether having observational or synthetic beam that exceeds the turbulence correlation length and, therefore, suppresses the correlated inertial part of the spectrum, will restore the validity of Eq. (75) if one samples the structure function at much smaller scales than the beam width. This is not so. Computing explicitly the structure function filtered with Gaussian window of width σ > L corr , one finds that D(R) ∼ 1 − exp(−R 2 /σ 2 ). Then at small lag D(R) ∼ R 2 , which is a general result, insensitive to the particular beam shape. 20 No constant contribution b appears, and Eq. (77) fails. Thus, the model adopted in Hildebrand et al. (2009) can be applicable only when the correlation scale of magnetic turbulence is smaller than the lag R between the points for which the structure function is calculated. For Kolmogorovlike turbulence this means R > L inj . Houde et al. (2009) have significantly extended the analysis of Hildebrand et al. (2009) with a detailed treatment of the telescope beam and the depth of observations as applied to molecular clouds, without using a restrictive model of Eq. (75). Nevertherless, their analytical progress relies on modelling the turbulence and polarization autocorrelation funcions in the specific form of a Gaussian, which width δ defines the correlation length. This model again does not reflect the actual scaling of the turbulence. In the language of structure functions, autocorrelation of a Gaussian form gives D(R) ∝ R 2 again at R < δ and D(R) ≈ const at R > δ, not the fractional power scaling that the real turbulence shows.
The analysis in Hildebrand et al. (2009) andHoude et al. (2009) is potentially applicable to turbulence that has most of energy at small scales. These "shallow" spectra with 3D k α , α > −3 have not been observed for magnetic or velocity turbulence, although we cannot exclude their presence in some special circumstances.

Our Explanation of the Observational Data
The approach described in Hildebrand et al. (2009 was applied to observational data to both molecular clouds (Houde et al. 2011;Chitsazzadeh et al. 2012;Houde et al. 2016) and also galactic disks . In view of other observational studies of galactic turbulence (see Armstrong et al. 1995, Chepurnov & Lazarian 2009, Chepurnov et al. 2010, Yuen et al. 2022a it is difficult to accept the setting given our arguments in the previous subsections. A possible explanation of some of the observational data in the aforementioned papers is that the small measured correlation scale of polarization angle fluctuations may correspond to the Alfvénic scale L inj M −3 A of super-Alfvénic turbulence. In fact, we observe this type of behavior in Fig.  25. For instance, in Chitsazzadeh et al. (2012) it was claimed that magnetic field fluctuations in OMC-1 have the correlation scales ∼ 9 " and ∼ 7 " for the Stokes Q and U parameters and ∼ 13 " for unpolarized intensity fluctuations. From the theory of super-Alfvénic turbulence (see Appendix A2), the first two numbers can be associated with the angular size associated with l A , and the third number with the injection size L inj . Using Eq. (B5) one can estimate M A for OMC-1 as (13/9) 1/3 to (13/7) 1/3 ≈ 1.13 − 1.22. This means that OMC-1 is a mildly super-Alfvénic object. The potential super-Alfvénic measurements suggest that the structure function analysis must be modified accordingly as the behavior of the structure function changes according to the theory of MHD turbulence. Note, in our paper we focus on the case of sub-Alfvénic turbulence. The super-Alfvénic case contains extra complication due to the non-trivial saturation on the slope of the polarization angle structure function. We shall discuss that case in later publications.

Issues with the Starting Equation
An attempt to take into account the effect of "compressible modes" within the DCF formalism was undertaken in Skalidis & Tassis (2021, henceforth ST21). Following the DCF the ST21 assumed the equi-partition of kinetic and magnetic energies at the injection scale, but attempted to take into account the parallel to mean magnetic field B 0 component of magnetic turbulent fluctuations δB ,turb present in the compressible media. Therefore the corresponding technique was termed in Li et al. (2021) "parallel-δB version of DCF". For the sake of brevity we shall refer to it as "parallel DCF" or Par-DCF".
In the Par-DCF, Skalidis et al. (2021, henceforth SX21) further suggested that when considering energy E ∝ B 2 , for δB 2 ,turb B 2 one can approximate and therefore the term 2BδB ,turb was associated with the fluctuation of magnetic energy, i.e.
However this is not right, since on average such δE mag will be zero as δB ,turb is not a positively defined, fluctuating in sign and, by the definition of the fluctuation, δB ,turb = 0. The multiplication of the fluctuation by B mean does not change the result and, naturally, B mean δB ,turb = 0.
To address the averaging to zero, Skalidis et al. (2021) proposed that the mean δB ,turb is to be substituted by root mean squared value δB 2 , which is mathematically not justified. It is obvious that the magnetic energy is ∼ δB 2 and not the cross product of the mean field and the fluctuation.
A recent publication of Beattie et al. (2022) comes up with a suggestion that we interpret as an attempt not to deal with magnetic energies as the DCF does, but consider a second moment of magnetic energy fluctuation. The physical meaning of this quantity is not clear, but it serves the purpose of preserving the desired cross term BδB . It is interesting that the presented numerical simulations are suggestive that the square root of the the second moment of the energy fluctuation is equal to the kinetic energy of the turbulence. This is a surprising result for which we do not see a physical justification. This result means that for subAlfvenic turbulence with B δB, the energy of magnetic fluctuations ∼ δB 2 are significantly lower than the energy of turbulent fluctuations ∼ δv 2 . This excludes the possibility that Alfvenic modes have significant contribution to turbulence and means that most of the energy must rest with the pure hydrodynamic motions that marginally involve magnetic field fluctuations. Such motions are compressions of fluid along magnetic field, which provides significant limitations on the type of media and driving that can result in such motions. This also at odds with earlier results e.g. in Haugen & Brandenburg (2004), , Kowal & Lazarian (2010), as well as the results in the present paper. Note, however, that the aforementioned numerical studies dealt with the incompressible driving of turbulence.
Future research should clarify numerous numerical issues that potentially can distort the results of MHD simulations. For instance, results in Haugen & Brandenburg (2004) indicate that the distribution between the energies of velocity and magnetic field fluctuations can be affected the ratio of the size of the numerical box to the injection scale. At low M A the effects of numerical box get more important and those can prevent establishing the regime of weak turbulence, as was demonstrated in Santos-Lima et al. (2021). If the relations used in Par-DCF can be justified from physical arguments, we do not believe that such arguments have been presented yet.

Par-DCF Expression for B-strength
This problem with the SX21's starting expression being problematic was first noticed in Li et al. (2021), where the authors numerically calculated the value B mean δB ,turb for their simulations and found that that the value of this term is very small, i.e. significantly smaller than δB 2 ,turb term. In particular in Appendix A3 of Li et al. (2021) it is written "In view of the questionable assumptions underlying their method, more work is needed to understand the physical basis for the method (of Skalidis & Tassis (2021))..." . This motivates us to further explore the SX21 expression for the magnetic field, even though, as we discussed earlier, we do not understand the physics behind this derivation.
Equating the kinetic energy with the expression of turbulent magnetic energy given by Eq. (79) ST21 obtained where, as in Eq.
(2) the fluctuations of velocities and angle are associated with the dispersions of these quantities that can be obtained from observations, i.e. σ v ≡ δv and σ φ ≡ δφ. Note that Eq. (80) is different from the DCF expression given by Eq.
(2). The most important is the difference between the two equations is that the square root of dispersion of σ φ enters in Eq. (80).
Our simulations employed for this paper are compressible MHD simulations with solendoial driving while SX21 employ both solendoial and compressive driving. However, SX21 claim that their finding do not depend on the difference in the compressibility of driving.
The study in SX21 is performed for low-β simulations and γ = π/2. The latter, as we have shown in this paper is a very special case for numerical testing. Indeed, at this γ due to the periodic boundary condition, there exist the strong suppression of the contribution of Alfvén modes to δφ. Unlike SX21, we perform the decomposition of the turbulence into Alfvénic, slow and fast modes. This allows us to quantify the effects of compressibility that the Eq (80) is supposed to reflect. Our results in Fig. 17,21 show that for low-β simulations the dispersion of σ φ is dominated by fast modes and scales as M A or M 5/4 A depending on the choice of the scales R. However, neither these dependencies correspond to σ φ ∼ M 2 A reported in SX21. Incidentally, SX21 claim that this difference of scaling of σ φ , i.e. being ∼ M A and ∼ M 2 A is "based on the different scaling relation with the magnetic fluctuations in the incompressible (Goldreich & Sridhar 1995) and compressible turbulence (Federrath 2016; Beattie et al. 2020)." Our numerical results contradict to this claim. It is Alfvénic and slow modes following Goldreich & Sridhar scaling that can show in some settings, e.g. for scales R L inj , ∼ M 2 A scaling. In contrast, compressible modes do not show σ φ ∼ M 2 A scaling in our simulations (See lower right corner of Fig.21). An additional discussion of our numerical results related to the applicability of Eq. (80) is provided in Appendix F.
We also can notice a resemblance of Eq. (80) and some of our expressions obtained for incompressible turbulence (see Table 3). Indeed, following the DMA approach, it is easy to see that one can rewrite Eq. (80) as However, this resemblance is coincidental. The SX21 study is focused on molecular clouds, which are low β medium. Therefore, the simulations that SX21 employ to test the correspondence with their formulae are compressible simulations in low β medium. For this setting our Table 3 provides f 0 ∼ M A for a very special case of M A 1 and γ → 0. However, the numerical study in SX21 was done for γ = π/2. In this case, the contribution of both Alfvén and slow modes are expected to be subdominant.
In Fig. 27 we provide the comparison of DCF results in the form of CY16 relations and Eq. (80) for a range of M A . Our results show that the CY16 show a more predictable behavior compared to Eq. (80). Indeed, the flat dependence in CY16 can be corrected by choosing a constant pre-factor on the basis of numerical simulations, as it done frequently in many DCF studies. Eq. (81), on the contrary, shows a dependence on M A which is indicative of a problem in the choice of the function in Eq. (81). Thus our testing favors the traditional DCF rather than Eq. (80).
Our prediction of f 0 ∼ M A corresponds to the case of nearly incompressible turbulence in high β medium (see Table 3). In this case the slow mode fluctuations induce changes of φ and, as a consequence, induce the f 0 ∼ M A dependence. This effect cannot be present in the SX21 setting.
We would like to state that our main point is that we object to the physical interpretation of Eq. (80) given in the literature. At the same time, we have no reason to question this equation as the consequence of empirically established fact. 11. DISCUSSION 11.1. DMA versus Davis-Chandrasekhar-Fermi Technique

Limitations of the DCF
The DCF is an empirical technique that is widely used in spite of the serious problems related to its accuracy in obtaining the value of magnetic field strength. In the past there were many attempts in fixing the DCF technique (see Ostriker et.al. 2002, Cho & Yoo 2016). However, these attempts were not founded based on the theory of MHD turbulence, in particular, the properties of the eigenmodes in MHD turbulence is usually disregarded or not properly considered. Figure 27. The ratio f = Btrue/B predicted that enters Eq. (81) as a function of MA for γ = π/2 for MHD simulations with Ms = 6. The simulations are from Table 6 for low-β.
The theoretical basis of the DCF is the equipartition of the kinetic energy and the energy of magnetic fluctuations, which is problematic in view of finding that the energy of magnetic fluctuations are significantly smaller than the kinetic energy in sub-Alfvénic turbulence (Haugen et al 2004). This equipartition is equal only at the regime of strong turbulence that we deal at small scales using the DMA.
Some of the problems of the DCF we addressed by introducing Eq. (17) that employs instead of global dispersion of δφ and δv the differential measures in the form of structure functions that can be measured locally. However, this does not solve a major problem of the DCF, namely, the DCF ignores the knowledge that we have about MHD turbulence, most importantly, the MHD turbulence anisotropy.
Reformulating the DCF on the basis of the modern MHD theory is difficult due to two different regimes of turbulence, weak and strong, that influence the dispersion of magnetic field and velocity. This problem is alleviated within the DMA that, by construction, deals only with strong turbulence sampling turbulent fluctuations at small scales.

DMA versus DCF
In this paper, we analytically derive an important relation that relates the squared mean magnetic field strength to the structure functions of the velocity centroid and polarization angle, with an extra f 0 factor that is related to the geometric factor of the MHD modes. Our general expression given by (Eq.(37)) describes obtaining magnetic field strength with any mixture of MHD modes. Later, we consider the limiting cases of turbulence corresponding to molecular clouds, i.e. low β medium, and diffuse warm gas, i.e. high β medium, analyze the effects of the change angle between the magnetic field and the line of sight. Based on the formalism in §5, We demonstrated that for the DMA can return accurate values of magnetic field strength (see Fig 23) for various astrophysical settings.
An additional advantage of the DMA compared to the DCF is that it can be successfully used for studying cases when line broadening is sub-thermal. The separation of the velocity components into the thermal and non-thermal part is rather complicated and causes additional errors (Eswaraiah et al. 2021). This limits the accuracy at which this process can be performed withing the DCF. At the same time the DMA does not require such separation, as the structure function of centroids is not sensitive to the thermal part of the line (see Lazarian & Esquivel 2003;Esquivel & Lazarian 2005;Kandel et al. 2017a). Moreover, The recently proposed Velocity Decomposition (VDA) approach (Yuen et.al 2021) allows the actual separation of non-thermal velocity fluctuations from channel maps. The fluctuations of densities induce changes to the statistics of the measured quantities, while the VDA allows to remove these interfering contributions. In Appendix H we provide an approach that allows decreasing the effects of density for the polarization studies.
In Paper II (Lazarian, Yuen, Pogosyan) we discuss another approach to finding magnetic field strength by employing Alfvén and sonic Mach numbers instead of the structure functions of magnetic field direction and the velocities that we use within the DMA. The two technique are new and the study of their synergy will be presented elsewhere.

Important Role of Various MHD Modes
The DCF technique disregards the fact that MHD turbulence consists of different MHD modes with slow and Alfvén modes being very anisotropic. To address this issue, in the DMA we study the turbulent velocity and magnetic field angle fluctuations at the small scales, i.e. at the scales at which we can apply our quantitative knowledge of compressible MHD turbulence (GS95, LV99, Lithwick & Goldreich 2001, Kowal & Lazarian 2010, see also Beresnyak & Lazarian 2019). We employ the statistics of MHD turbulence in the observer's frame mostly following the LP12 approach and somewhat extending it.
In the DCF the accuracy of the magnetic strength determination is increased by the empirical adjusting the expression prefactor f 0 by using numerical simulations. On the contrary, the DMA provide expression of f 0 as a function of M A , angle γ and properties of the media at hand. This significantly increases the precision of the new technique and explains the reported uncertainties of the DCF.

DMA and Recent Attempts to Modify DCF
The well-known limitations of the DCF induce numerous efforts to modify the technique. We have reviewed a couple of recent ones applying in sequence the following criteria: • Solid physical model behind the derivation.
• Model of turbulence applied for deriving the expression.
• Advantages of the derived expression compared to the original DCF.
For instance, we showed in §6.3 that the ST21 study does not satisfy the first criterion. Their approach if treated as based on the numerical finding should not be dismissed but subject to more numerical studies. The differences that we find interpolating our results can be due to the differences in the adopted driving.
The first criterion is satisfied for Angle Dispersion Function (ADF) method in Hildebrand et al. (2009);. Incidentally, the ADF, similar to DMA, employs structure functions for the polarization angle, but stays within the DCF approach with both using the line width for evaluation of velocity fluctuations as well as returning to the use of dispersions of magnetic field and velocities in the final expressions. On the contrary, the DMA uses the small scale structure functions of both polarization angle and the velocity measures.
At the same time the vital difference between the DMA and the ADF is that the latter implicitly assumed that the injection scale of the turbulent velocity L inj is much smaller than the separation of scales at which the structure functions are calculated. This assumption does not agree with what we know of turbulence in molecular clouds. For instance, according to McKee & Stone (2021) the injection of turbulence scale is at least two times the size of molecular clouds, while in Yuen et al. (2022a) it is shown that the turbulence spectrum of velocity fluctuations in Taurus smoothly transfers to much larger scales in diffuse media, making the turbulence injection scale significantly larger than the scale of the molecular clouds. As a result the applicability of the approach the applicability of the ADF is limited to special cases that should yet be identified by observations.
The improvement of the DCF technique that corresponds to all three criteria above is given in Cho & Yoo (2016). Similar to the DMA it employs centroids of velocities. Formally, our Eq. (17) transfers into Cho & Yoo (2016) expression in the limiting case of l = R L inj . However, the analytical expressions derived in the DMA approach cannot be transferred to be used in Cho & Yoo (2016), as the DMA deals with the well defined properties of MHD cascade at small scales, which are different from the properties of dispersion δv and δφ at the injection scale that Cho & Yoo (2016) employs. In addition, the distribution of magnetic field strength over the observed object is available with the DMA.

Obtaining Magnetic Field Strength in 3D
Using Dust polarization. GAIA data on the distances from the stars provide a valuable source of the 3D distribution of magnetic field in the galaxy. This information was used, for instance, in Gonsalvez-Casanova & Lazarian (2019) to confirm that the 3D distribution of the POS magnetic field obtained with the an innovative Velocity Gradient Technique (VGT) (see Yuen & Lazarian 2017ab, Lazarian & Yuen 2018) is consistent with the star light polarization data. The direct application of the DCF to the distribution of polarization directions from the stars with known distances is prohibited by the fact that the velocity dispersion of the nearby gas is dominated by the galactic rotation. This is not a problem for the DCF, however. For instance, the DMA could be applicable to studies of magnetic field using the 21 cm line of atomic hydrogen. This line is broadened by both thermal motions and also galactic rotation, but one can still use structure functions of velocities using Reduced Velocity Centroids (RVCs) introduced in Lazarian & Yuen (2018a).
In fact, the DMA can be used directly with the stars without using any diffuse media spectroscopy. Indeed, it was recently found that the statistics of the velocity of young stars reflects the statistics of turbulence (Ha et al. 2021). Therefore the structure functions of the velocities of stars can be obtained with GAIA and combined with the starlight polarization data. The advantage of the GAIA data, that the structure functions of velocity are available not only for the LOS component of velocity, but also for the POS velocity components. As a result, more detailed studies can be performed using the DMA formalism. 21 Using Spectral line polarization. Magnetic field direction can be obtained by measuring the polarization from spectral lines that arises from Goldreich-Kylafis effect and the Ground State Alignment (GSA) (see Appendix B). The polarization of spectral lines exhibit properties that makes them attractive from the point of view of 3D magnetic field strength measurements. First of all, one can ensure that the measurements of the magnetic field direction and the Doppler broadening arise from the same volumes of gas. Very importantly, however, that the emission of different spectral lines is localized to different regions. For instance, the GSA is being destroyed by collisions and thus it produces polarization in rarefied regions of interstellar medium with strong illumination by radiation sources. As a result, unlike dust polarization, the GSA is less affected by the line of sight averaging. In addition, spectral lines are affected by galactic rotation and this allows to get separate the magnetic field measurements at different distances along the line of sight. As we discussed above, the DMA can work successfully in the presence of Doppler broadening arising from galactic rotation.
Using magnetic field directions from VGT. Spectral lines can be used to obtain the magnetic field directions using the VGT. This does not require any polarization information and the same spectral line information can be used both for evaluating the structure functions of the magnetic field and the structure functions of velocity. The VGT, however, provides yet an alternative way of measuring magnetic field strength that is based on the statistical properties of gradients. We discuss this in Paper II (Lazarian, Yuen, Pogosyan). The DMA-VGT and the other technique are complementary in obtaining the 3D distribution of magnetic field strengths.

Synergy of theory and numerical simulations
In this paper we employed numerical simulations to both to test the theoretical expectations and provide some of the input parameters. For instance, our numerical studies showed how the energy is being distributed between different turbulent modes for different settings.
Being guided by theory, we expected the non-trivial change of the results with angle γ between the line of sight and mean magnetic field. Therefore our numerical testing was performed for a variety of M A and γ. In contrast, the DCF numerical simulations that we are aware of, are limited to γ = π/2, which our study shows to be a very serious limitation.
Theory provides important warnings for interpreting numerical simulations. We list of couple of them below.
For instance, numerical testing faces problems arising from the insufficient resolution of numerical simulations. This problem is very serious for Alfvén and slow modes. As the consequence of insufficient numerical range, these modes for sufficiently small M A can be simulated only in one regime, the weak MHD turbulence regime. In astrophysical reality, the inertial range is usually sufficiently extended and both regimes of weak and strong MHD turbulence are present. This entails a serious difference in the properties of actual and numerically simulated turbulence. In this situation the numerical testing of both DCF and DMA becomes poorly justified.
In addition, for γ = π/2 the periodic boundary conditions induce the effect of infinite integration that strongly suppresses the Alfvénic contribution. This is something that should be seriously considered during the DCF testing. At the same as γ changes from π/2, the effects of periodicity along the line of sight is altered and the suppression of Alfvén waves is removed. This is a subtle but important effect that complicates the testing. The periodicity is being disturbed for γ = π/2 and, as a result, the suppression of Alfvénic modes gets altered. Thus small changes in γ around π/2 can result in significant changes of the observed contribution of Alfvén modes. These effects have never been discussed within the numerical studies of DCF where setting of γ = π/2 is employed On the contrary, in the case when fast modes dominate the measured fluctuations of φ, there exist only one regime of fast mode turbulence and therefore we see the scaling of D φ that corresponds to the change of the fast mode energy with M A , i.e. while the fast mode energy increases with M A in proportion to M A . Indeed, we observed this dependence for γ = π/2. There exist uncertainties related to the known properties MHD turbulence that still require studies. For instance, we are not certain to what extend the turbulence at small scales gets independent on the conditions of turbulent driving in terms of its mode composition. Naturally, at the injection scale, the distribution of energy withing compressible and incompressible modes is determined by the turbulence driving. However, one may expect that due to the partial coupling of modes (Cho & Lazarian 2003), the energy can be redistributed between them as turbulence cascades to small scales. If so, for sufficiently extended turbulence range, one may get turbulence properties independent of the initial driving. This issue calls for further studies. 11.5. Prospects of the Technique and a Broader Impact of the present study The expressions of the DMA allow to derive the Bstrength with any required precision. However, they exhibit dependencies on on angle between the line of sight and magnetic field γ, medium magnetization β and mode composition.
Potentially, these parameters can be obtained from independent studies. For instance, in Yuen et al. (2022b) the procedures for obtaining the relative distribution of MHD modes are proposed and tested with synthetic observations. The plasma β does not need to be obtained precisely and its evaluations can be improved via the iterative application of the DMA. When these parameters are uncertain, the DMA expressions allow to evaluate the uncertainties of B-strength obtained with the technique.
The new technique is really timely these days where both polarimetry and velocity measurements can done with high spacial resolution. This allows to measure more detailed statistics compared to the earlier days. Thus, with the DMA one can get detailed information about the magnetic field and its distribution over the turbulent astrophysical volume.
The DMA samples motions at smaller scales at which the motions preserve its Alfvénic character. This it opens prospects for obtaining the magnetic field strength in super-Alfvénic turbulence that we will explore in our next publication (Yuen et al in prep).
The current study is based on the LP12 formalism describing MHD turbulence in the laboratory frame. The first application was the analytical description of synchrotron intensity fluctuations. The part of it related to velocities was further advanced in the series of papers that followed (Kandel et al. , 2018.
In the present paper, the LP12 approach was elaborated for describing the observed statistics of polarization directions. The importance of this advance goes beyond the particular problem in hand, i.e. the problem of obtaining magnetic field strength by combining polarization and spectroscopic data. For instance, the developed formalism allows to address the studies of magnetic turbulence from observations (see . In addition, we noticed the limitation of LP12 study in terms of dealing with high β media case. Indeed, the equipartition of energy in slow and Alfvén modes is assumed in our earlier studies. This is a reasonable assumption an idealized infinite inertial range, when the measurements of turbulence properties lose their dependence on turbulent driving. However, for the measurements sufficiently close to the in-jection scale, the properties of turbulence are influenced by the driving. This is definitely the case for present day numerical simulations. For instance, for the case of the conventional isotropic turbulence driving, at low M A , the energy of slow modes is expected to be ≈ 1/2 of the energy of Alfvén modes. The difference between the two is expected to decrease with the increase of M A , which we also confirmed by numerical simulations. Therefore, in the present study we extended the LP12 approach for the case of isotropic turbulence driving for high β case. This extension is important for the applications that are different from the DMA. 12. SUMMARY The paper introduces a new way to measure the strength of magnetic field by combining the spectroscopic data and the polarization measurements. The gist of our approach is to use the differential measures, i.e. structure functions at smallest available scales, to characterize the fluctuations of both the velocity and the polarization. A distinguishing feature of this study is the use the theory of MHD turbulence to describe the relevant small-scale fluctuations, at R significantly smaller than the injection scale L inj This is in contrast to the traditional Davis-Chandrasekhar-Fermi (DCF) where the dispersion of the above two quantities were employed and the anisotropic properties of MHD turbulence are disregarded. As a result, our Differential Measure Analysis (DMA) technique is different from the traditional DCF and its more recent modifications/improvements. We can briefly summarize our results in the following way: 1. The DMA can be applied to small patches of observational data and can successfully deal with data inhomogeneity and interfering processes not related to the turbulent cascade. As a result, the DMA can provide a detailed distribution of the POS component of magnetic field.
2. The anisotropic nature of MHD turbulence makes the DCF approach not accurate. Our study focuses on the small scale asymptotic behavior of basic modes of MHD turbulence and provides robust expressions for magnetic field. The generalization of our results to the dispersion that the DCF deals with is challenging, due to changes of the nature of MHD cascade at large scales. On the basis of our study we can state that the coefficient of proportionality between the magnetic field strength and the ratio of the velocity and polarization angle structure functions is not a constant, but, in general, a function of Alfvén Mach M A , angle between the mean magnetic field and the line of sight γ, as well as of relative fraction of basic MHD modes in the turbulent volume.
3. We obtained general expressions for the DMA both in interstellar medium with magnetic pressure larger than the gaseous pressure, i.e. low β medium, and for gaseous pressure larger than the magnetic pressure, i.e. high β medium. Starting with our general expressions, we derive a set of simplified expressions that are applicable to magnetic field studies in molecular clouds and diffuse media (see Table 3).
4. Our study of high β medium provides robust expressions that can be applied to observational data with minimal assumptions. For instance, in the case of equipartition of Alfvén and slow modes, no additional information related to the value of the angle γ is required.
5. Our study of low β medium demonstrates pronounced dependence on γ. Therefore the earlier numerical studies of the DCF case limited to γ = π/2 are not adequate. The composition of turbulence in terms of basic MHD turbulence modes can significantly alter the results.
Our study provides general expressions that can be used for obtaining magnetic field by combining polarization and spectroscopic observations. It testifies that a further increase of the accuracy of obtaining of magnetic field strength in molecular clouds can be achieved by employing additional information, e.g. the information on the composition of MHD turbulence in terms of Alfvén, slow and fast modes. This information can be obtained both from numerical simulations and observations.

A. WAYS OF MEASURING MAGNETIC FIELD USING POLARIZATION A.1. Polarization from Aligned Dust
Dust polarization arises from emission of non-spherical grains aligned with long axes perpendicular to the ambient magnetic field (see Andersson et al. 2015). Similarly, polarization of starlight arises from the differential extinction by aligned grains. The processes of dust alignment is generally believed to happen due to radiative torques (RATs) (see Dolginov & Mytrophanov 1976;Draine & Weingartner 1996 ). The theory of the RAT alignment have is based on the analytical model in Lazarian & Hoang (2007) and further studies e.g. in Hoang & Lazarian (2008, 2016. The RAT alignment theory at its present form (see Lazarian & Hoang 2019) can account for the major observational features of grain alignment. In particular, in typical conditions of diffuse ISM the silicate grains are nearly perfectly aligned, while in dense molecular clouds the degree of alignment depends on the grain illumination mostly by embedded stars. In other words, the existing grain alignment theory can evaluate in what conditions one should expect the polarization arising due to the aligned dust to trace magnetic fields. With more polarization measurements obtained using starlight and with more distances to stars measured there is a possibility to trace magnetic field in 3D.
A.2. Goldreich-Kylafis Effect Goldreich, & Kylafis (1981;1982, henceforth GK) effect provides a viable way of tracing magnetic fields in molecular clouds. The polarization arises due to the differences of the radiation transfer in the media with anisotropies or shear. The resulting polarization is either parallel or perpendicular to the magnetic field. In spite of this ambiguity, the effect has been successfully employed to trace magnetic field structure of molecular clouds (Girart et al 1999, Li et.al 2011 A more recent study reveals that there are extra complication in understanding the molecular line polarization based on GK effect. It is proposed that the use of circular polarization (Houde et.al 2013, Chamma et.al 2018 can remove the correlation between the polarization angle and the magnetic field (See also Hezareh et al. 2013) In any case, combining GK with velocity gradients (Yuen & Lazarian 2017ab, Lazarian & Yuen 2018a, Yuen et.al 2021 one can remove the 90 degree ambiguity in the magnetic field direction. A recent work in analysing the polarization in the emission of atomic or molecular (sub)millimeter lines has been developed by Lankhaar & Vlemmings (2020).

A.3. Ground State Alignment
A promising development in the RMS of magnetic field tracing is presented by the atomic/ionic ground state alignment (GSA) effect suggested and quantified for use in astrophysical conditions by (Yan & Lazarian 2006Zhang et al. 2015;. 22 The GSA employs atoms/ions with fine and hyperfine split levels. The atoms/ions get aligned in the ground or metastable state by external anisotropic radiation. The Larmor precession in the ambient magnetic field re-aligns the atoms/ions imprinting its direction on polarization. The atoms/ions stay in ground or metastable state long and thus they can trace very weak magnetic fields. The effect has been recently confirmed with observations (Zhang et al. 2020a), opening a wide avenue of applying it for tracing magnetic fields in various environments. The difference in distribution of atoms and conditions for atomic alignment in space provides a way to get the 3D distribution of magnetic field in diffuse medium. The technique is especially interesting for probing magnetic field direction near bright sources. Here we briefly summarize the scaling laws in the local frame of reference for compressible MHD turbulence as we did in . If the energy is injected with the injection velocity V L that is less than the Alfvén speed V A , the turbulence is sub-Alfvénic. In the opposite case it is super-Alfvénic. The illustration of turbulence scaling for different regimes can be found in Table 5. We briefly describe the regimes below. A more extensive discussion can be found in the review by Brandenburg & Lazarian (2013) is termed the weak Alfvénic turbulence. This type of turbulence keeps the l scale stays the same while the velocities change as v ⊥ ≈ V L (l ⊥ /L inj ) 1/2 (Lazarian & Vishniac 1999) The cascading results in the change of the perpendicular scale of eddies l ⊥ only. With the decrease of l ⊥ the turbulent velocities v ⊥ decreases. Nevertheless, the strength of non-linear interactions of Alfvénic wave packets increases (see Lazarian 2016). Eventually, at the scale l trans , the turbulence turns into the strong regime which obeys the GS95 critical balance. The situations when the l trans is less than the turbulence dissipation scale l diss require M A that is unrealistically small for the typical ISM conditions. Therefore, typically the ISM turbulence transits to the strong regime. If the telescope resolution is enough to resolve scales less than l trans then we should observe the signature of strong turbulence in observation.
The anisotropy of the eddies for sub-Alfvénic turbulence is larger than in the case of trans-Alfvénic turbulence described by GS95. The following expression was derived in LV99: where l and l ⊥ are given in the local system of reference. For M A = 1 one returns to the GS95 scaling. The turbulent motions at scales less than l trans obey: i.e. they demonstrate Kolmogorov-type cascade perpendicular to local magnetic field. For the magnetic field it is more natural to rewrite this expression as where the value of mean field B 0 enters explicitly.
In the range of [L inj , l trans ] the direction of magnetic field is weakly perturbed and the local and global system of reference are identical. Therefore the velocity gradients calculated at scales larger than l trans are perpendicular to the large scale magnetic field. While at scales smaller than l trans the velocity gradients follow the direction of the local magnetic fields, similar to the case of trans-Alfvénic turbulence that we discuss in the main text. Turbulent models are characterized by an inertial range over which power-spectrum is power-law. The inertial scale is bounded by the injection scale L inj above which the spectrum is assumed to sufficiently quickly flatten out, or even decrease, and short dissipation scale L dis below which there is a fast cutoff of the power. For sub-Alfvénic turbulence, characterized by Alfvénic Mach number M A < 1, magnetic field correlation length is similar to injection scale, however the cascade has two regimes, of weak turbulence between L inj and L trans ≈ L inj M 2 A and that of strong turbulence, for scales between L trans and L dis . For super-Alfvénic case, M A > 1, correlation scale is shorter L A ≈ L inj M −3 A , and the magnetic field spectrum is flattened before reaching L inj .

B.3. Super-Alfvénic Turbulence
If V L > V A , at large scales magnetic back-reaction is not important and up to the scale (see : the turbulent cascade is essentially hydrodynamic Kolmogorov cascade. At the scale l A , the turbulence transfers to the trans-Alfvénic turbulence described by GS95 scalings, i.e. anisotropy of turbulent eddies start to occur at scales smaller than l A . C. STATISTICAL VARIATIONS IN THE CASE OF SMALL DATA SETS The two point global structure function at R = l itself represents the dispersion of the observable with block size of l during computation averaged across the maps with length L, which we draw an illustrative figure in Fig.28. Suppose we have an observable X within the circular area L and we sample the dispersion of X by selecting a smaller area with length scale l ignoring all other interfering factors (e.g. shot noise, dissipation ranges etc, telescope beams smoothing). The individual dispersion values of X within the length scale l will not coincide with the structure function at R = l. However, if we take adequate averaging of the areas with R = l, the value of such averaging will return to the values as drawn in structure function (blue point).
What do we meany by "adequate" here? Beresnyak & Lazarian (2005) considered the structure function calculations by selecting a partial number of statistics in simulations, and illustrate that the "adequacy" in recovering the structure functions is just a few % of the data. Yuen et al. (2018) further showed that even punching ∼ 50% of the data the anisotropy of structure function would not be altered. This means that potentially we can divide our map into sub-regions and compute the structure functions in obtaining a distribution of values of δX for magnetic field strength estimations. D. NUMERICAL TESTING Most of the numerical data cubes are obtained by 3D MHD simulations that is from a single fluid, operator-split, staggered grid MHD Eulerian code ZEUS-MP/3D to set up a three dimensional, uniform, isothermal turbulent medium. To simulate the part of the interstellar cloud, periodic boundary conditions are applied. These simulations use the Fourier-space forced driving solenoidal driving. 23 For isothermal MHD simulation without gravity, the simulations are scale-free. If V inj is the injection velocity, while V A and V s are the Alfvén and sonic velocities respectively, then the two parameters, namely, the Alfvén Mach numbers M A = V inj /V A and sonic Mach numbers M s = V inj /V s , determine all properties of the numerical cubes and the resultant simulation is universal in the inertial range. That means one can easily transform to any arbitrary units as long as the dimensionless parameters M A , M s are not changed. The chosen M A and M s are listed in Table 6. For the case of M A < M s , it corresponds to the simulations of turbulent plasma with thermal pressure smaller than the magnetic pressure, i.e. plasma with 23 Our choice of force stirring over the other popular choice, i.e. of the decaying turbulence, is preferable because only the former exhibits the full characteristics of turbulence statistics, e.g power law, turbulence anisotropy, extended from k = 2 to a dissipation scale of 12 pixels in a simulation , and matches with what we see in observations (e.g. Armstrong et al. 1995;Chepurnov & Lazarian 2010). β/2 = V 2 s /V 2 A < 1. In contrast, the case that is M A > M s corresponds to the thermal pressure dominated plasma with β/2 > 1. To investigate the behavior of the incompressible case, we adopt the incompressible cube our previous work Lazarian et al. (2017). Further we refer to the simulations in Table 6 by their model name. The selected ranges of M s , M A , β are determined by possible scenarios of astrophysical turbulence from subsonic to supersonic cases.
As our derivations employed the properties of velocities and magnetic field, we first deal with the numerical testing of our predictions without accounting for the effects of density fluctuations. This is possible in observational applications due to the development of Yuen et.al (2021). Yuen et.al (2021)  Two point statistics of general random vector field, such as the turbulent velocity field, can described by the structure function that in Fourier space is given by three power-spectra, correspondent to orthogonal motions of the medium that can be considered statistically independent Which orientation of the orthogonal triad ξ A , ξ S , ξ F of unit vectors describes independent motions depends on the physical properties of the medium and velocity excitations. In the magnetized turbulence three independent modes are identified to be Alfvén, Slow, and Fast (which explains our notation with A, S and F calligraphic labels, see Table.4) which have particular orientations w.r.t. to the direction of the mean magnetic fieldλ. MHD motions are statistically anisotropic, with both the mode decomposition dependent structurally on the direction of the magnetic field, and power spectra, generally, dependent on the relative angle of the wave vector with the mean field axis.
In the high-β plasma, β 1, the Alfvén, Slow, Fast decomposition is and coincides with A, F, P decomposition into two solenoidal and one potential modes introduced in the main text. Namely, in this regime Fast motions are purely potential, while Slow mode together with Alfvén mode are two components of the solenoidal motion. For plasma β finite, and, in particularly in low-β plasma, Alfvén modes remain the same, but Slow and Fast directions are rotated relative to ξ F and ξ P ξ A =k ×λ where the rotation angle α β depends on the plasma constant in such a way that α β → 0 as β → ∞ and α β → arcsin(µ) as β → 0. The following formula gives a qualitatively correct fit In the low-β plasma, β → 0, the Alfvén, Slow, Fast decomposition is then

E.2. Magnetic Field Turbulent Fluctuations
Properties of the magnetic field fluctuations δB are linked to medium velocities v in MHD turbulence by the field frozen condition where the wave frequency ω(k) ∝ k ·λ for Alfvén and Slow modes. For the fast modes, ω(k) ∝ k and is angle independent.
Magnetic field is solenoidal in nature and contains only two degrees of freedom which are conveniently described via A-typê ξ A and F -type ξ F contributions, as in the following explicit expression 24 The frozen condition Eq. (E11) projects three, A, S, F modes of motion into two types of perturbations of the magnetic field, A and F . While Alfvén motions A give rise directly to A-type perturbations of B-field, both Slow and Fast motions cause F -type magnetic perturbations only. This is true for all values of plasma β, but what relative contributions of Slow and Fast modes are in F -type δB, depends on β.
Namely, in the high-β plasma, β 1, we obtain, using Eqs. (E11) and (E7) where v A is Alfvén velocity, which in the regime is also the speed of slow modes, while v F is the speed of fast modes.
In the low-β regime Slow modes do not perturb the magnetic field An important case of high-β turbulence is that of strong, nearly incompressible turbulence, when A v (k, µ) = S v (k, µ) and the magnetic perturbations arising from fast modes are negligible. 25 In this case the turbulence is described by a single power spectrum with isotropic tensor structure where Alfvén and Slow modes give two polarizations of a transverse solenoidal wave. The variance of the magnetic field fluctuations can be obtained as half of the trace of the structure function at large separations where A 0 (k) is the monopole of the power spectrum A, A 0 = 1 4π dΩ k A(k,k ·λ), and, correspondingly, F 0 (k) is the monopole of the spectrum F .

E.3. 2D Projected Structure Functions
The projected 2D structure function is obtained by integration over LOS z direction, where 2D vectors orthogonal to LOS are introduced in capitalized notation as k = (K, k z ) and r = (R, z). The integral over z translates into δ D (k z )-like behaviour if the LOS integration range exceeds the POS scale R under study. In this approximation, subsequent integration over k z amounts to setting k z = 0, and the result of LOS projection can be obtained from Eq. (E13) by substituting k → (K, 0), includingk → (K, 0), andλ = (sin γΛ, cos γ) whereΛ is the 2D direction of the mean field on the sky. Under this transformation µ = sin γ cos φ K where cos φ K =K ·Λ is the cosine of the angle between POS projections of the wavevector and the symmetry axis.
D ij (R) = 1 2π 2 d 2 K 1 − e iK·R A(K, sin γ cos φ K ) δ ij −K iKj + sin 2 γΛ iΛj − sin 2 γ cos φ K (K iΛj +K jΛi ) 1 − sin 2 γ cos 2 φ K + +F (K, sin γ cos φ K ) −K iKj +K iKj + sin 2 γΛ iΛj − sin 2 γ cos φ K (K iΛj +K jΛi ) 1 − sin 2 γ cos 2 φ K (E18) For this paper we are interested in the structure function of the perturbations of the magnetic field component perpendicular to the mean field. Considering projected mean field to be along the x-direction,Λ = (1, 0), this is the y-component for which D yy (R) = 1 2π 2 d 2 K 1 − e iK·R A(K, sin γ cos φ K ) cos 2 γ cos 2 φ K 1 − sin 2 γ cos 2 φ K + F (K, sin γ cos φ K ) sin 2 γ cos 2 φ K sin 2 φ K 1 − sin 2 γ cos 2 φ K (E19) For completeness we shall also quote D xx (R) = 1 2π 2 d 2 K 1 − e iK·R A(K, sin γ cos φ K ) cos 2 γ sin 2 φ K 1 − sin 2 γ cos 2 φ K + F (K, sin γ cos φ K ) sin 2 γ sin 4 φ K 1 − sin 2 γ cos 2 φ K (E20) D zz (R) = 1 2π 2 d 2 K 1 − e iK·R A(K, sin γ cos φ K ) + F (K, sin γ cos φ K ) cos 2 γ 1 − sin 2 γ cos 2 φ K In what follows we shall consider two specific models of the turbulence. One of them is model (I) of strong incompressible turbulence where Alfvén and slow modes are equally distributed. This corresponds to A(k) = F (k). The other, model (A), is a pure Alfvénic turbulence F (k) = 0. In addition, we also give the expressions for A(k) = 0 case (F) which represents magnetic field perturbations in both slow and fast modes in high and low β regimes. To simplify notation, let us use E(k) to designate the power spectrum generally, without specific reference to one of the modes. Then we can write D yy (R) = 1 2π 2 d 2 K 1 − e iK·R E(K, sin γ cos φ K )G (I,A,F ) (γ, φ K ) , where      G I = cos 2 φ K G A = cos 2 γ cos 2 φ K 1−sin 2 γ cos 2 φ K G F = sin 2 γ cos 2 φ K sin 2 φ K 1−sin 2 γ cos 2 φ K (E22) One can point out that the model of strong turbulence utilizes all two degrees of freedom that general representation allows for, while the purely Alfvénic model is degenerate, using only one degree of freedom. This leads to some artifacts in mathematical approximations, as for instance in our idealized axis-symmetric approximation the vanishing of the projected structure function if cos γ = 0, i.e when the mean magnetic field is exactly perpendicular to LOS. In reality this means that in such particular configuration Alfvénic mode will give subdominant contribution. We shall focus on the multipole coefficients D n yy (R) = 1 2π 2π 0 dφ R D yy (R, φ R )e −inφ R for which we obtain, after performing integration over φ K D n yy (R) = where E 2D n (K, sin γ) are coefficients of 2D multipole expansion on the sky of the projected k z = 0 power spectrum (see Lazarian & Pogosyan 2012;Kandel et al. 2017b) with respect to angle φ K , while G (I,A,F ) n (γ) are similar multipoles for geometrical functions in Eq. (E22), with both sets of coefficients potentially dependent on sin γ. Decomposition of the geometrical factor is particularly simple in the strong incompressible case, G I n = 1 2 δ n0 + 1 4 (δ 2n + δ −2n ). In the following Section E.5 we discuss that the level of anisotropy of the power spectrum can be often considered scaleindependent, i.e the dependence of the power on the direction of the wave-vector can be separated as With this factorization, Eq. (E22) can then be presented in a general form where the formal expression for the scaling functions I n (R) is I n (R) = π KdK (δ n0 − i n J n (KR)) E 0 (K) k 2 dkE 0 (k) (E26)

E.4. 2D Structure Function of LOS Turbulent Velocities
In velocity centroid studies (see Kandel et al. 2017a for details) one looks at the projection of the line-of-sight component of the velocity, dzv z (z). In this projection the potential mode vanishes, so the relevant 2D structure function has the form D v zz (R) = 1 2π 2 d 2 K 1 − e iK·R × × A v (K, sin γ cos φ K ) 1 − cos 2 γ 1 − sin 2 γ cos 2 φ K + F v (K, sin γ cos φ K ) cos 2 γ 1 − sin 2 γ cos 2 φ K (E27) In high-β case, with velocity spectra A v = F v as usual for strong nearly incompressible regime (I) or, individually, for A and F modes, we get that differs from Eq. (E22) only in geometrical functions W. We note that in LOS projection, the F mode in the high-β regime is provided completely by the slow mode motions, while the fast mode motions, that are potential, are projected out.
In low-β case, while the form Eq. (E28) for the D v zz (R) remains the same, the relevant geometrical functions are the Alfvén W A that is unchanged from Eq. (E28) and, now separately, slow and fast mode contributions to the F part W S = (1 − µ 2 )W F = cos 2 γ . W F = µ 2 W F = cos 2 γ sin 2 γ cos 2 φ K 1 − sin 2 γ cos 2 φ K (E29) that could be easily found from Eq. (E10). Velocities have the same scaling as the magnetic field fluctuations, E v (k) ∝ E(k). Following the steps of the previous section, the multipole representation of D v zz is then with the same I n (R) scaling function as in Eq, (E25). Here we have dropped the zz label as we are always dealing only with z component of the velocity. In the strong incompressible case W I n = δ n0 .

E.5. Angular Anisotropy of the Turbulent Power Spectrum
Let us focus on sub-Alfvénic case in strong turbulence regime. We can model E(k) power spectrum as E(k, µ = k · λ) = A trans (kL trans ) −3−m E(k, µ), kL trans > 1 with m being index of power scaling, and angular dependence factored out in E(k, µ) that is defined normalized to have 3D monopole equal to unity, 1 2 1 −1 dµ E(k, µ) = 1. The spectral amplitude A trans at the transition scale can be linked to the amplitude at injections scale via weak turbulence scaling with index m w ≈ 1.
A trans = A inj (L trans /L inj ) mw (E32) A inj , in turn, determines the variance of the magnetic field thus E(k, µ = k · λ) = δB 2 π 2 m w L 3 inj (L trans /L inj ) mw (kL trans ) −3−m E(k, µ), kL trans > 1 The level of anisotropy, in particular the quadrupole of E(k, µ) can be scale dependent. The strong turbulence relation Eq. (B2) for m = 2/3 gives ratio l /l ⊥ of individual eddies in configuration space to be scale dependent, increasing for smaller ones ∝ l −1/3 . However, CLV02 has argued that if one considers the ensemble averaged power spectrum in the volume L 3 c < L 3 trans , the averaging of orientation of individual small eddies results in the scale independent anisotropy of the spectrum, determined by the shape of the largest eddies at scale L c , which satisfies Eq. (B2). Switching to the angle coordinates of the wavevectors relative to the magnetic field, so that k = L −1 c , l = k −1 = L c /µ and l ⊥ = k −1 ⊥ = L c / 1 − µ 2 , we have If L c is raised to L trans = L inj M 2 A , Eq. (E35) suggests the following form for the angular dependence of the power spectrum This form is valid for each separate L 3 trans volume, however each different volume will have different direction of the volumeaveraged magnetic field, i.e., the axis from which µ angle is measured. Variations of this mean axis are determined by long fluctuating modes k < L trans which are in the weak turbulence regime. Global definition of the power spectrum will average over such long-wave "wandering" of the preferred axis and will describe the residual anisotropy of the spectrum with respect to the globally mean magnetic field.
The effect of long wave field wandering can be estimated by averaging the power spectrum over the distribution of λ directions. We will model the distribution of λ direction by assuming that the perturbations in the magnetic field δb are Gaussian and over large scales distributed in 3D isotropically around λ 0 global mean direction with the variance δb 2 /B 2 0 = M 2 A , which gives 26 where θ λ ∈ (0, π/2) and φ λ ∈ (0, 2π) are spherical angles of the headless vector λ with respect to λ 0 axis. Let us note that in the isotropic limit M A → ∞, cos 2 θ λ → 1/3 with respect to an arbitrary axis. In the opposite limit, M A → 0, P(θ λ ) → δ D (θ λ ) and λ → λ 0 with sin 2 (θ λ ) ∼ M 2 A . For M A < 1, the result of the averaging can be approximated as Let us also quote a multipole expansion of the 2D projection of this function E 2D p = (−1) p/2 2 exp − 1 2 M −2 A⊥ I p/2 F. RELATIVE ROLE OF AND ⊥ MAGNETIC FIELD FLUCTUATIONS FOR LOW β CASE While the main body of the paper deals with the scaling of the projected fluctuations of magnetic field, here we will discuss the properties of 3D fluctuations. These fluctuations are known to be very anisotropic in MHD turbulence and it is only perpendicular to the mean magnetic field component of magnetic fluctuations that contributes to the observed fluctuations in angle. 27 For our MHD simulations with solenoidal driving, we show our results for the relative amplitude of magnetic fluctuations in Fig. 14. We see that for low β supersonic M s = 6 and M A = 0.1 turbulence, the fraction of energy in fast modes is around 6% 26 This expression retains qualitative features but is an improvement on LP12 27 This point was initially missed in ST21 and corrected in SX21. There it was assumed that turbulent fluctuations are statistically isotropic, i.e. δB = δB ⊥ for all M A < 1. This is a very strong assumption that contradicts to what we know about MHD turbulence. For instance, Alfvénic fluctuations are perpendicular to the local direction of magnetic field. The contribution of slow modes is mostly to δB and it is marginal for turbulence in low-β plasmas, while it is only fast modes that are isotropic.
of Alfvén mode energy. At the same time, while the total fraction of energy in slow modes can be comparable or even exceed the energy in Alfvén modes, the contributions of the slow modes in the 3D fluctuations of magnetic field is also marginal. Thus, in our simulations for turbulence in low β medium the dispersion of δB ⊥ scales as ∝ M A in agreement with the DCF approach. Fig.29 shows how the 3D magnetic field fluctuations relative to the guide field vary as a function of the global M A in the numerical cubes. For low β case we can observe the fluctuations arising from slow modes being almost parallel to the ambient field. The contributions of these fluctuations to the variations of magnetic field direction is negligible and for simulations at hand is on the order of accuracy of the decomposition technique. At the same time, we observe that the fluctuations from Alfvén modes that are perpendicular to the mean field, i.e. responsible for the variations of φ, scale ∝ M A . Figure 29. A figure showing how the 3D magnetic field fluctuations for each fundamental MHD mode (in the global frame of reference) varies as a function of MA for low β simulations. Blue points: Alfvén mode, green points: slow mode, red points: fast modes. The trend lines are added for readers' references. The size of the marker is ∝ β , where we provide a reference marker of size β = 0.1 at the legend.

G. SIMILARITIES AND DIFFERENCES OF SYNCHROTRON AND DUST POLARIZATION STATISTICS
In the current paper we discuss about mainly the statistics of dust polarization, which has been covered in §4 in extended detail. However readers might recognize that the fundamental theory of Lazarian & Pogosyan (2012) relies on the "synchrotron" formalism. Readers might wonder how similar and different are these formalism in our development of DMA?
Let us recall how the dust and synchrotron polarization depends is theoretically written in the form of P = Q + iU : P dust ∝ dzn dust e 2iθ B sin 2 γ P synch ∝ dzn re B 2 tot e 2iθ B sin 2 γ where we have to remind our readers that tan θ B = B y /B x , γ is the line of sight angle, and n re is the density of relativistic electrons. Notice that the relativistic electron density is well known to be less correlated to the ISM density, while the dust density is believed to the highly correlated to the ISM density (See Draine 2011) and a justification is required to disregard the dust fluctuations (either via Yuen et.al (2021), or via the treatment in §4). The fundamental differences between the dust and synchrotron polarization formulae (Eq.G41) is the emergence of the sin 2 γ term in the dust polarization formula, while the strong ∝ B 2 tot dependence as discussed in Lazarian & Pogosyan (2012). In our current paper, we are not concerned with the actual statistics of Q and U in general, but the polarization angle φ = tan −1 2 (U/Q). The easiest way to compare their performance is to compute the dispersion of polarization angle using the two formula (Eq.G41). Fig.30 shows how the dispersion of polarization angles vary as a function of M A . We can see that when M A 1 the difference between the dust and synchrotron formalism is infinitesimal, which is expected by theory (See Lazarian & Pogosyan 2012). However, when M A > 1, they start have have slight divergence. When M A 1, the synchrotron polarization formula could be approximated as P synch ≈ B 2 tot dz... which the remaining factor within the integral carries the same factor as the dust formula, and the B 2 tot term cancels when computing φ. However, when M A > 1 the aforementioned relation does not hold, and therefore we expect deviations between the two expressions. However as one can see from Fig.30 that the difference between them is rather minor, we expect the scaling in the main text would change only slightly. Here we discuss how the corresponding procedure can it be implemented for polarization data and what is the constrains the procedure accuracy. The form of the decomposition equation that is employed for the polarization is the same its velocity version. In Yuen et.al (2021) it was shown that the employed linear-algebra construction works best for subsonic media. We illustrate this effect using our numerical simulations. From Yuen et.al (2021), we can write the velocity caustics p v as where p is the channel map, and I is the total column density map. Very similarly, we can also compute the Stokes variant version: where I here is the unpolarized Stokes intensity. Notice that the formulae (Eq.H43) has to be tested to see whether they work. Fig.31 shows the decomposition of Q and U in subsonic turbulence while that of Fig.32 shows the same decomposition in Stokes parameters synthesized from supersonic turbulence. We can see that the result that we obtained here is consistent with the spectroscopic equivalent in Yuen et.al (2021): The subsonic decomposition can recover the constant density Stokes parameter with the Normalized Correlation Coefficient (NCC, see Yuen et.al 2021) to be very close to 1, while that of the supersonic case we can somehow recover the structure but the NCC drops below 0.6. This simple illustration suggests that we can remove the density contribution easily and with high accuracy in the case of subsonic turbulence.
Notice that our technique also applies to any single-frequency synchrotron emissions without significant Faraday rotation, because the density dependence in both dust and synchrotron emissions are the same. However, without multi-frequency emissions it is impossible for us to obtain the density-free information in supersonic turbulence.