© ESO 2021

1 Introduction

Characterizing the magnetic field properties of the interstellar medium (ISM) is an important task for various fields in astrophysics. Dust polarization is the most widely used magnetic field tracer, since others such as the Zeeman effect (e.g Heiles 1997; Crutcher et al. 2010) or the Goldreich-Kylafis effect (Goldreich & Kylafis 1981, 1982) are usually dominated by systematic uncertainties. A major limitation of dust polarization is that it only directly probes the magnetic field morphology and not the strength of the field (e.g., Andersson et al. 2015). For this reason, various methods for the indirect estimation of the magnetic field strength from dust polarization data have been developed.

Davis (1951) and Chandrasekhar & Fermi (1953) (DCF) were the first to propose that polarization data can be used to estimate the ISM magnetic field strength. They suggested that magnetic field fluctuations should be imprinted on the polarization map. If there is an ordered magnetic field component which is stronger than the fluctuating one, then the magnetic field lines should be highly ordered, that is the dispersion of polarization angles (δθ) should be low. On the other hand, when fluctuations are stronger or comparable to the mean field component, then magnetic field lines should be highly dispersed and δθ large. They assumed that δθ are due to the propagation of incompressible magnetohydrodynamic (MHD) waves, known as Alfvén waves. The total magnetic energy density of a cloud is (1)

where is the mean field energy density which exists even if there are no fluctuations in the cloud and the other two terms correspond to the fluctuating magnetic field component. Alfvén-wave fluctuations are always perpendicular to B₀, and thus B₀ ⋅ δB = 0. DCF further assumed that gas kinetic energy is completely transferred to magnetic energy fluctuations and obtained that, (2)

Observationally, δθ is a proxy of δB∕B₀ (Zweibel 1996), and hence the above equation can be rearranged in the following form, (3)

where f is a constant factor inserted to account for various biases in the DCF method: only one of three Cartesian velocity components perturbing the field lines (Chandrasekhar & Fermi 1953); line-of-sight (LOS) averaging of the polarization signal (Zweibel 1990; Myers & Goodman 1991); averaging within the telescope beam (Falceta-Gonçalves et al. 2008; Houde et al. 2009), etc. To estimate an appropriate empirical correction for these effects, DCF has been calibrated against numerical magnetohydrodynamic (MHD) simulations (Ostriker et al. 2001; Heitsch et al. 2001; Padoan et al. 2001; Liu et al. 2021) and a factor f = 0.5 is usually applied to Eq. (3).

However, the ISM is highly compressible (Heiles & Troland 2003) which implies that the DCF approximation of turbulence incompressibility may not be generally applicable. Skalidis & Tassis (2021) (ST) relaxed this assumption and proposed an alternative method for estimating the magnetic field strength from polarization data in the presence of strongly magnetized and compressible turbulence. They suggested that in compressible but sub-, trans- Alfvénic turbulence the coupling term¹, B₀ ⋅ δB, is not zero and dominates in Eq. (1). As a result, it is the coupling term which should be comparable to the gas kinetic energy rather than δB²∕(8π), and hence (4)

Thus they derived that (5)

Equations (3) and (5) have a different scaling dependence on δθ. The strength of B₀ in the DCF equation is inversely proportional to δθ, while in ST it is inversely proportional to the square root of δθ. These scalings can be equivalently expressed as scalings of δθ with the Alfvén Mach number M_A, which is the ratio of δu over the Alfvén speed (V_A). For DCF, (6)

while for ST, (7)

This difference is due to the inclusion (ST) or omission (DCF) of the B₀ ⋅ δB term. Skalidis & Tassis (2021) compared the performance of Eqs. (3) and (5) against MHD simulations with M_A = 0.7, and they showedthat Eq. (5) is significantly more accurate than Eq. (3), even when the calibration factor f is included in the latter. However, they did not test the scaling of δθ with M_A, as this would require data from a larger number of simulations with varying values of M_A. The scope of this work is to test the δθ – M_A scaling, and thus the validity of DCF and ST assumptions regarding the contribution of B₀ ⋅ δB in the cloud energetics.

In addition, we address the following theoretical subtlety. Both DCF and ST are statistical methods and represent averaged quantities. From this point of view, one should compare the average kinetic energy density with the average magnetic pressure, (8)

For incompressible turbulence, one can directly derive the classical DCF equation, Eq. (3), with the substitutions and , with δθ probing . For compressible turbulence B₀ ⋅ δB≠0, since it is connected with the density fluctuations (Bhattacharjee & Hameiri 1988; Bhattacharjee et al. 1998). However, for a periodic signal the number of rarefactions (δB ⋅B₀ < 0) is equal to the number of compressions (δB ⋅B₀ > 0), and hence ⟨δB ⋅B₀⟩ = 0. This on a first reading implies that DCF (Eq. (8)) also applies to compressible turbulence. And yet, the raw DCF method (without an appropriately calibrated f factor) is found to be highly inaccurate when tested in compressible MHD simulations: the estimated B₀ values using Eq. (3) and f = 1 systematically deviate from the actual value by large factors. The method produces reasonable estimates only when f ≤ 1∕2 (Ostriker et al. 2001; Heitsch et al. 2001; Padoan et al. 2001; Skalidis & Tassis 2021; Liu et al. 2021), which means that DCF estimates without a fine-tuned f are more than 100% larger than the actual values. Even when combined with more sophisticated techniques, such as the dispersion function analysis (Hildebrand et al. 2009; Houde et al. 2009, 2013) which takes into account the LOS effects, DCF produces highly biased estimates (Skalidis & Tassis 2021). On the other hand, Skalidis & Tassis (2021) found that Eq. (5) is significantly more accurate in simulations with compressible turbulence than DCF, whether the latter includes or not the correction factor f. Therefore the argument that the B₀ ⋅ δB term cancels out in the averaged energetics of strongly magnetized and compressible turbulence does not hold.

In this work, we investigate where the argument fails, and establish that Eq. (5) of ST indeed provides a much superior estimate of the magnetic field strength whenever compressible turbulence is non-negligible, and a good estimate of the magnetic field strength (better than 50%) across all cases. In summary, we address the following questions: i) What is the actual scaling between δθ and M_A in MHD simulations?, ii) Does the B₀ ⋅ δB term contribute in the cloud energetics when averaging over the total volume of a cloud?, and iii) How accurately can we probe the contribution of B₀ ⋅ δB with polarization data?

To address these questions, we significantly expanded the set of numerical tests of both methods (DCF, ST) to a wide variety of forced ideal-MHD, isothermal simulations without self-gravity, with different values of Alfvénic Mach numbers M_A (ranging between 0.1 and 2.0) and sonic mach number M_s (ranging between 0.5 and 20), and with different driving of the turbulence. We produced synthetic polarization data from all simulations, and we investigated the scaling of δθ with M_A, the energetics of the simulations, and the behavior of averaged quantities. Finally, we applied the DCF and ST methods to these data, and we investigated the accuracy of the methods in estimating the strength of B₀.

This article is organized as follows. In Sect. 2 we present the numerical simulations used in this work. In Sect. 3, we test the DCF and ST scaling relations between δθ and M_A in synthetic data. In Sect. 4, we show why the omission of B₀ ⋅ δB in the energetics of compressible turbulence is physically wrong and verify our arguments with numerical data. We also test how accurately we can trace the B₀ ⋅ δB fluctuations from polarization data. In Sect. 5, we apply both methods (DCF and ST) in synthetic data and test their accuracy in simulated clouds covering a wide range of M_A and M_s. In Sect. 6, we present the main conclusions of this work.

2 Numerical simulations

We used data from the following simulations in our tests. “Cho-ENO” (Cho & Lazarian 2003; Burkhart et al. 2009; Portillo et al. 2018; Bialy & Burkhart 2020) simulations from the publicly available CATS database (Burkhart et al. 2020): These are ideal-MHD, isothermal simulations without self-gravity. Turbulence is driven in velocity Fourier space by injecting solenoidal modes only at scales equal to half the size of the simulated cube. Models are characterized by M_A = 0.7 and 2.0, while M_s ranges between 0.7 and 7.0. Skalidis & Tassis (2021) have tested the two methods in the M_A = 0.7 simulations of this dataset, but we also included them in our results for completeness. Simulation data are dimensionless and scale-free. A dimensionless sound speed, which is defined as , regulates the units. We assume that the sound speed is 0.91 km s⁻¹ for every model and follow Hill et al. (2008) in order to convert to cgs units. The resolution is 256³.

“AREPO” simulations from the CATS database (Burkhart et al. 2020): The simulations setup is presented in detail in Mocz et al. (2017); Burkhart & Mocz (2019). These are isothermal, ideal-MHD simulations run with the AREPOcode (Springel 2010). Turbulence is driven solenoidally until a quasi-static state was reached with M_s = 10 and then self-gravity is switched on. We used the model with M_A = 0.35 at a time step without self-gravity. The resolution of this model is 256³.

Simulations from Beattie et al. (2020, 2021). They solve the ideal-MHD equations without self-gravity and isothermal conditions using a modified version of the FLASH code (Fryxell et al. 2000; Dubey et al. 2008; Federrath et al. 2021). Turbulence is driven in Fourier space by injecting the same amount of power between compressible and solenoidal modes at large scales. Sound speed is c_s = 1 in every model and gas velocities are expressed in M_s units. The Alfvénic and sonic Mach numbers cover a wide range of the parameter space, M_A = 0.1–2.0 and M_s = 0.5–20. Models with M_s = 0.5 have resolution equal to 576³, while every other model 512³.

Simulations of Körtgen & Soler (2020). These are ideal-MHD, isothermal simulations without self-gravity, which were run with the FLASH code (Fryxell et al. 2000). These simulations are in cgs units with T = 11 K, that is sound speed is equal to 0.2 km s⁻¹, and n = 536 cm⁻³. We used the model with M_A = 0.5 and M_s = 7.5 driven solenoidally. The resolution is 512³.

Boundary conditions are periodic in every simulation presented here. In total we used 26 MHD numerical simulations with properties summarized in Table 1.

Fig. 1 Polarization angle dispersion as a function of the Alfvénic Mach number. Blue line: ST scaling; magenta line: DCF scaling. The two lines are normalized so that they pass through the data for MA = 1.0.

3 Testing DCF and ST scalings with numerical simulations

Skalidis & Tassis (2021) tested the DCF and ST methods in numerical simulations characterized by a unique M_A value (0.7) and five different M_s values. However, ISM turbulence spans a wide range of M_A and M_s values. There is overwhelming observational evidence that ISM turbulence is sub-, trans- Alfvénic (e.g., Mouschovias et al. 2006; Franco et al. 2010; Pillai et al. 2015; Panopoulou et al. 2016; Planck Collaboration Int. XXXV 2016; Cox et al. 2016; Tritsis & Tassis 2018) and highly compressible (e.g., Heiles & Troland 2003; Miville-Deschênes & Martin 2007; Brunt 2010; Burkhart et al. 2015; Orkisz et al. 2017; Nguyen et al. 2019; Beattie et al. 2019). This means that ISM clouds are statistically characterized by M_A ≤ 1 and M_s > 1. The numerical simulations employed in this work (Sect. 2) enable us to test the two methods to a large number of MHD simulations with a wide range of parameters consistent with observations.

The major difference between DCF and ST is the dependence of the magnetic field strength on δθ. In DCF, B₀ scales as δθ⁻¹ (Eq. (3)), while in ST B₀ scales as ~ δθ^−1∕2 (Eq. (5)). One can divide both equations with and obtain the magnetic field strength in velocity units (the Alfvénic speed, V_A). The scaling relations of DCF and ST expressed in terms of V_A are V_A ~ δθ⁻¹, and V_A ~ δθ^−1∕2, respectively.The Alfvénic speed is (9)

where c_s is the sound speed. Using the above equation with the corresponding δθ scaling dependence of the two methods we obtain (10)

This is the key difference between the two methods, and it is based on the different scaling relation of M_A with the magnetic fluctuations in the incompressible (Goldreich & Sridhar 1995) and compressible turbulence (Federrath 2016; Beattie et al. 2020). We tested the two scalings in synthetic polarization data. We computed synthetic Q and U Stokes parameters (Eqs. (20) and (21) in Sect. 4.6) and the polarization angles as θ = 0.5arctanQ∕U. Then, we computed the dispersion of the θ angles, δθ. All δθ values are shown in Table 1.

In Fig. 1, we show δθ as a function of M_A (we do not show the solenoidally driven models with M_A = 2.0 since their δθ is not representative of the actual fluctuating-to-ordered magnetic field ratio, Sect. 5.2). The blue line corresponds to the ST scaling, and the magenta line to the DCF scaling. Both lines are normalized so that they pass through the simulated data for M_A = 1.0. The ST scaling shows a remarkable consistency with the data. In contrast, the DCF scaling fails to representthe synthetic data. Even if we fine-tune the DCF relation with the use of a factor f à la Ostriker et al. (2001), agreement with the data is achieved only around the M_A values used for the tuning. The scaling slope of DCF is clearly inconsistent with simulation data, regardless of the presence of f.

4 Theoretical considerations

In the previous section, we have shown that ST outperforms DCF over a large range of Alfvénic Mach numbers. The simple theoretical arguments to the contrary we discussed in the introduction should, then, be incorrect. In this section we discuss why this is so, from a theoretical perspective.

4.1 Does the B₀ ⋅ δB term contribute to the average energetics of compressible turbulence?

The coupling term is by definition zero in the incompressible regime and the averaged total energy (kinetic and magnetic) in the perturbations/waves is, (11)

This is similar to the energy equation of a harmonic oscillator, where energy fluctuates between kinetic and potential (here magnetic) forms. The physical analogy in incompressible Alfvénic turbulence works well: the magnetic field oscillates harmonically about B₀. The lowest value of the potential energy is achieved at the equilibrium state where δB = 0. We can thus consider each fluid element as a harmonic oscillator perturbed around . According to the ergodic theorem, time averaging is equivalent to spatial averaging (equivalently, we can say that within the cloud there exist all possible oscillation phases); Eq. (11) therefore holds – but, in addition, as in the harmonic oscillator, the average values of kinetic and potential energy are equal, ρ₀ ⟨δu²⟩∕2 = ⟨δB²⟩∕8π (Eq. (2)), whence DCF Eq. (3) is obtained.

In compressible and strongly magnetized turbulence, the fluctuations are also periodic. The B₀ ⋅ δB term, although much higher in absolute value than δB², can be either negative or positive. Therefore, ⟨B₀ ⋅ δB⟩ = 0, Eq. (11) holds, and it would appear that once again the problem can be reduced to that of a harmonic oscillation. However, in this case the physical analogy is incorrect. There are two reasons for this.

First, unlike a harmonic oscillator, the equilibrium state (δB = 0) is not the state of lowest potential energy; the maximum rarefaction state (δB = −|δB|_max) is. The cross term, B₀ ⋅ δB, can be either positive or negative and this means that locally it can either add or remove magnetic flux from the fluid elements. When B₀ ⋅ δB > 0, then the mean energy density, , locally increases by 2|B₀ ⋅ δB|∕(8π) due to the compression of the magnetic field lines. On the other hand, in regions where B₀ ⋅ δB < 0 the mean energy density locally decreases by −2|B₀ ⋅ δB|∕(8π) due to the decompression of the magnetic field lines. In contrast, in a harmonic oscillator (as well as incompressible Alfvénic turbulence),any deviation from the equilibrium position will only increase the potential energy (the magnetic energy, in the case of our fluid elements).

Second, the dependence of the “potential energy” on the perturbation (here of the magnetic field) is linear: δϵ_{p,compressible} ∝ δB). In contrast, in a harmonic oscillator the dependence is quadratic: δϵ_{p, harmonic} ∝ δB². Therefore, neither δu nor δB of a fluid element will behave harmonically with time. It is thus clear that the harmonic oscillator is not the appropriate physical analog to our problem. Can we substitute it with a more appropriate mechanical analog to guide our intuition?

Indeed we can. Let us consider a shaft with depth |h| below the Earth’s surface. We let a basketball fall from height + h above the Earth’s surface into the shaft. In classical Newtonian mechanics, the ball will accelerate from a height + h down to the bottom of the hole, − h, where it willbounce upward. In the absence of energy losses, the ball executes the reverse motion as if the clock now runs backwards, and the ball will again reach + h, before moving once more downwards toward the bottom of the shaft, continuing these oscillations forever. Consider now that we make the choice of taking the zero point of gravitational potential energy to be at the Earth’s surface, at its midpoint between its highest value (achieved at + h) and lowest value (achieved at − h). Now the potential energy mgz is positive above the surface of the Earth (z > 0), negative below the surface of the earth (z < 0), and its average value over an entire cycle is ⟨mgz⟩ = 0. At the same time, the kinetic energy is always non-negative, ⟨mv²∕2⟩ > 0. It is obvious that in this case, unlike the harmonic oscillator, ⟨mv²∕2⟩≠⟨mgz⟩. Equating them would lead to an absurdity. The absurdity is resolved when we compare the absolute maximum potential energy (that is the difference in potential energy between highest and lowest points) with the maximum kinetic energy²: yields the correct relation between the maximum positive height h, and maximum velocity, v_max.

The analogy with B₀ ⋅ δB oscillations works very well if we replace the basketball with a fluid element and the gravitational potential energy with magnetic energy, where δB_∥ now plays a role similar to the height of the bouncing ball. The magnetic field of fluid elements oscillates around B₀, just as the ball height fluctuates around z = 0. For this reason B₀ ⋅ δB can be positive or negative, that is fluctuations in magnetic energy can be positive or negative if its zero point is defined at δB_∥ = 0 (B = B₀), just as fluctuations in the gravitational potential energy of the ball can be positive or negative if its zero point is defined at z = 0. If oscillations are periodic, we obtain ⟨B₀ ⋅ δB⟩ = 0. Still, we should not be concluding that the coupling term does not contribute to the average energy budget, any more than we should conclude that the potential energy of the bouncing ball does not contribute to its average energy budget.

4.2 The ST method revisited

For each fluid element energy oscillates between kinetic and magnetic forms. The total energy density of a fluid element is, (12)

where E_tot = E₀ + δϵ. In the unperturbed case the total energy of each fluid element is , which is the mean magnetic energy density. When |B₀|≫|δB|, the perturbed energy of each fluid element is (13)

in the fluid rest frame where u₀ = 0. The δB²∕8π term is second order, and hence it was neglected. The δϵ are energy fluctuations around E₀ and can be negative when δu → 0 and B₀ ⋅ δB < 0. If we assume undamped oscillations, then the kinetic energy is completely transferred to magnetic and backwards periodically. It is reasonable to assume that kinetic energy fluctuations are dominated by δu² fluctuations when the temperature is constant within a cloud implying ρ ≈ ρ₀.

To make further progress, we should identify the physically correct way to relate the kinetic energy term with the first-order magnetic energy term. As in the case of the bouncing ball, the maximum kinetic energy of the fluid element will be comparable to the absolute maximum magnetic energy (that is the difference between maximum compression and maximum rarefaction): (14)

Of course, neither δu_max nor δB_∥max can be probed observationally. The quantities that we do have access to from observations are the spatially-averaged kinetic energy fluctuations, ⟨δu² ⟩ and magnetic field fluctuations (Sect. 4.6), which, by virtue of the ergodic theorem, correspond to the time-averaged fluctuations over an entire period of the evolution. What we need then is a way to relate δu_max to ⟨δu² ⟩, and to δB_∥max.

In the case of the bouncing ball, the time evolution of its velocity v and height z are straightforward to obtain, so we can in fact calculate these relations between and ⟨δv² ⟩, and between z_max and . Before we do so, however, and use them in the problem at hand, we should investigate how far we can take the analogy between bouncing ball and fluid element in compressible, strongly magnetized turbulence. Would the time behavior of the fluid element have in fact the same functional form as in the bouncing ball?

The answer is “yes”, provided that we can write a formally equivalent Lagrangian for the two systems, and show that the boundary conditions of the problem are similar. If we take the magnetic field perturbation parallel to B₀, δB_∥, to be a generalized coordinate for the problem, then the generalized velocity would be , which, by virtue of flux freezing, is proportionalto δu_⊥ (the velocity of the fluid element perpendicular to the magnetic field)³. In the bouncing ball case, we have a potential energy term that is proportional to z, and a kinetic energy term that is proportional to ż². Similarly, for the fluid element responsible for a magnetic field compression or rarefaction, we have a potential energy term that is proportional to δB_∥, and a kinetic energy term that is proportional to . The Lagrangians of the two problems are thus formally equivalent. The boundary conditions of the problem are also similar: the presence of the bulk medium and its large-scale magnetic field acts as the “ground”, forcing the fluid element with increasing velocity undergoing a rarefaction to reverse course back toward increasing magnetic field with velocity decreasing in magnitude. As a result, the time evolution profiles of δB_∥ and will be similar to those of z and v, respectively, in the bouncing ball problem.

For the bouncing ball, over one period of the motion and ⟨z² ⟩ = 7h²∕15. Therefore, by eliminating in favor of ⟨v²⟩, and h in favor of , the relation between maximum kinetic and potential energies, , can be rewritten as (15)

Since the dynamics of the two systems are analogous, we can also express Eq. (14) for compressible fluctuations as⁴, (16)

from which we derive, (17)

4.3 Comments on the assumptions and approximations of the ST method

In this section, we discuss some assumptions entering the ST method which can limit its accuracy. Firstly, we assumed that oscillations are undamped in compressible turbulence, but in reality shocks lead to significant energy loss and nonideal effects (e.g., ambipolar diffusion) can induce significant loss of magnetic flux from the cloud. Even in ideal-MHD simulations there is significant energy dissipation due to the presence of shocks and numerical diffusion. In these cases, equipartition between the root-mean-square (rms) kinetic and magnetic energy is not guaranteed. In nature, equipartition can hold when there is a constant energy source at large scales injecting energy to the cloud. In numerical simulations this is achieved with the so called forcing, applied in most cases. Secondly, we assumed that magnetic pressure dominates over gas pressure. This is a reasonable approximation for clouds with M_A < 1 and M_S > 1, but when B₀ ⋅ δB < 0 the magnetic pressure may locally become comparable to the gas pressure (especially in cases when B₀ ⋅ δB is close to its negative minimum). However, an oscillator spends only a limited amount of time at its negative minimum, which through ergodicity means that gas pressure will be important only for a small volume fraction of the cloud. Thirdly, to arrive to Eq. (16) we assumed that gas density fluctuations are much smaller than u², and hence that density is approximately constant in the integral, ρ ≈ ρ₀. Thus, Eq. (17) omits any contribution from density fluctuations in the kinetic energy, ⟨δρu²⟩∕2.

Fig. 2 Relative ratio of kinetic over magnetic energy density using Eq. (8) (vertical axis) and Eq. (16) (horizontal axis). Colors correspond to different MA as shown in the legend. Dots correspond to simulations with Ms < 1, “x” to 1 < Ms ≤ 4 and triangles to Ms > 4.

4.4 How good are the energy equipartition assumptions of DCF and ST?

In order to assess the validity of Eqs. (8) and (16), we tested them in numerical simulations. We computed the kinetic energy term () and compared it against the two different magnetic energy terms.

In Fig. 2, we show in the vertical axis the ratio of kinetic over magnetic energy, with both terms calculated from Eq. (8) (DCF). In the horizontal axis, we show the ratio of kinetic over magnetic energy, with both terms calculated from Eq. (16) (ST). Different color points correspond to simulations with different M_A, while different shapes correspond to simulations with different M_s. Subsonic (M_s < 1) simulation results are shown with dots, supersonic models with 1 < M_s ≤ 4 are shown with an “×” and M_s > 4 models are shown with triangles. The vertical line indicates exact equipartition à la ST and the horizontal line exact equipartition à la DCF. Diagonal lines separate the regions where each method outperforms the other. Green-shaded regions are closer to the horizontal than the vertical line; there, it is the quadratic term in magnetic energy that dominates and is better comparable to the kinetic energy. White-shaded regions are closer to the vertical line; there, it is the coupling term in the magnetic energy that dominates, and it is that term that is closer to the kinetic energy.

There is a strong dependence of these results on M_A. The coupling term dominates the energetics of sub-Alfvénic simulations. This has been already shown in numerical simulations of strongly magnetized and compressible turbulence (Federrath 2016; Beattie et al. 2020). There is a slight difference between Eq. (16) and the one from Federrath (2016) and Beattie et al. (2020). We include only fluctuations which are parallel to B₀, while they included the total δB rms, and we include only the perpendicular velocity rms in the kinetic energy, while they included the total velocity rms; however, the overall behavior is the same.

Our results suggest that Eq. (8) (DCF) is highly inaccurate (by factors of several to 100) in sub-Alfvénic simulations, performs comparably to Eq. (16) (ST) in trans-Alfvénic simulations, and clearly outperforms Eq. (16) (ST) only in trans-, super- Alfvénic turbulence simulations driven solenoidally (see Fig. 3 and also discussion in Sect. 4.5). This shows the weakness of the incompressible approximation to accurately describe theenergetics of sub-Alfvénic turbulence and indicates that the DCF method is precarious to use when the M_A of a cloud is unknown. On the other hand, ST is fairly accurate (better than factor of 2) in the entire sub-, trans- Alfvénic regime, even for high M_s cases. Overall, the coupling term dominates in the energetics over the δB² term when M_A < 1 and cannot be ignored. This is consistent with our discussion in Sect. 4 and highlights the importance of this term when estimating the magnetic field strength.

Fig. 3 Same as in Fig. 2, but the dots correspond to models driven with a mixture of compressible and solenoidal modes and “×” denotes models driven solenoidally.

4.5 How does the forcing affect the energetics?

Turbulence in MHD simulations is driven in order to achieve the desired M_s. Driving is implemented by injecting compressible or incompressible or a mixture of modes in the cloud through a stochastic process in Fourier space. This process is supposed to mimic the driving mechanisms found in nature. We explored if there is any dependence of the results shown in Fig. 2 on the forcing mechanism.

In Fig. 3, we show the same data as in Fig. 2, but here dots correspond to models driven with a mixture of compressible and incompressible modes, and “×” represents solenoidally driven simulations. The driving mechanism of each simulation is shown in the second column of Table 1. The majority of simulations used in the current work were driven with an equal mixture of modes, but there are also a few driven solenoidally. Overall, our results in Fig. 3 are weakly affected by the forcing mechanism in the sub-Alfvénic simulations.

We compared the effect of forcing at models with M_A = 0.5 and M_A = 2.0, since we have an overlap of both solenoidally and mixed driven simulations. Blue points correspond to models with M_A = 0.5. The solenoidally driven model (shown with the blue “×”) has the largest kinetic over ratio. The reason is that forcing injects only incompressible modes in the cloud which are only traced by δB². On the other hand, when mixed forcing is used the kinetic energy is equally shared among and δB². For this reason, the blue “×” point is shifted toward larger and smaller values in the horizontal, and vertical axis respectively compared to the blue dots. However, the difference between the forcing mechanisms does not create significant deviations between these models. In sub-Alfvénic turbulence B₀ is much stronger than the forced perturbations and determines how energy is transferred among the modes.

The effect of forcing in super-Alfvénic simulations is more prominent. We compared the models with M_A = 2.0 (shown with cyan) and mixed forcing (denoted by dots) with the solenoidally driven models (denoted by ×). Models with mixed forcing are clustered in the bottom left corner of Fig. 3, while models with solenoidal forcing in the right bottom corner. The reason is that solenoidal modes are not represented by δB_x, and hence all the injected kinetic energy goes to δB²∕(8π). On the other hand, when mixed forcing is applied the injected energy is shared. In super-Alfvénic turbulence |B₀| < |δB|, and hence the dynamics of the cloud are determined by δB, instead of B₀. In these cases the mode of the forced fluctuations determines the cloud dynamics. Thus, there is a large separation between the solenoidal and mixed driven simulations with M_A = 2.0 in Fig. 3. Magnetic field lines are highly curved in the solenoidally driven simulation, and hence . On the other hand, in mixed driven simulations . This is in striking contrast with the sub-Alfvénic turbulence where fluctuations evolve independently of the forcing mechanism and fluctuations are dictated by B₀ and not by δB.

4.6 Can we trace B₀ ⋅ δB with polarization data?

The DCF and ST methods can be used to estimate the strength of B₀ using the incompressible and compressible approximation respectively, provided that we can estimate from observations ρ₀, , and (for DCF) or (for ST). Observationally, the gas volume density is one of the most uncertain parameters inserted in these equations. In the diffuse ISM, ρ₀ can be estimated using specific gas tracers, such as the 158 μm fine structure transition line of [C II] (e.g., Langer et al. 2010; Goldsmith et al. 2018). The rms velocity canbe estimated using the turbulent broadening of the dominant gas tracer at each cloud and the fluctuating-to-ordered magnetic field ratio from the dispersion of the polarization angle data, δθ.

One concern raised about the ST method is that δθ is considered to be tracing only fluctuations perpendicular to the mean magnetic field (e.g., Zweibel 1996) and, as such, is dominated by the Alfvénic modes. However, in the ST method (Eq. 17) it is the term that is inserted in the energetics, rather than . In this section we explore if this indeed poses a problem for the ST method, or whether we can indeed probe parallel magnetic fluctuations from δθ.

Consider an ISM cloud permeated initially by an undisturbed and homogeneous magnetic field B₀ = (B₀, 0, 0). We perturb B₀ with δB = (δB_x, δB_y, δB_z) and we assume that δB is random, and hence ⟨δB⟩ = 0. The total magnetic field is B_tot =B₀ + δB, (18)

If the LOS is parallel to the z-axis, then the projected magnetic field morphology of B_tot as traced by dust polarization will be given by the Stokes parameters (Lee & Draine 1985), (19) (20) (21)

where |B| denotes the total strength of the field. The above equations hold when dust grain properties are uniform and temperature is constant throughout the cloud. For simplicity we assume that density is uniform along each LOS and that |B|² variations along the LOS are negligible (that is that the ρ₀ and terms strongly dominate ρ and |B|², respectively, and thus the latter can be pulled out of the integrals). From Eqs. (18) and (20), we derived for the Q Stokes parameter (22)

where . Using Eqs. (18) and (21), the U Stokes parameter can be expressed as (23)

If |B₀|≫|δB|, then we can drop the second-order terms, obtaining (24) (25)

The polarization angle is tan2θ = U∕Q. In the |B₀|≫|δB| regime θ is small and tan2θ ≈ 2θ. From Eqs. (24) and (25) we obtain, (26)

The δB_x term in the denominator of Eq. (26) is due to turbulence compressibility. In the limit where perturbations are limited only to Alfvén waves, this term is by definition zero and θ traces the perpendicular fluctuations of the magnetic field fluctuations (Zweibel 1996).

The dispersion of polarization angles is , since , where brackets here denote averaging in the x–y plane. For convenience we adopt the following notation:

Thus, we obtain, (30)

Since |B₀|≫|δB|, is of second order, and hence δθ is simplified to, (31)

The term is due to δB_x and represents the average coupling of B₀ with δB.

It is therefore clear that although ΔB_x does contribute to the dispersion of δθ, this contribution comes from a first-order term in the denominator of the average in the RHS of Eq. (31). When |B₀|≫|δB|, parallel fluctuations have but a limited contribution in δθ and can be neglected. Thus, to first order, δθ represents perpendicular fluctuations: .

However, in the ST method, it is parallel fluctuations that need to be estimated and used in Eq. (17). Only when would δθ be an adequate metric of parallel fluctuations. The question then becomes: do parallel and perpendicular fluctuations have similar dispersions in compressible turbulence?

The answer is “yes”, as we show in Fig. 4, where we have plotted the ratio as a function of M_A. We find that in all cases the dispersion of perpendicular fluctuations is comparable to that of parallel fluctuations, with deviations always smaller than a factor of 2. These results are consistent with the findings of Beattie et al. (2020). We note that inFig. 4 we are displaying the square root of the dispersions of the LOS averages of perpendicular and parallel fluctuations, respectively, since these are the relevant quantities in relating our observable (δθ) with the quantity of interest for the ST method (δB_∥).

Fig. 4 Ratio of perpendicular over parallel LOS averaged magnetic field fluctuations.

5 Applying the DCF and ST methods in synthetic data

In this section, we explore the accuracy of the DCF and ST methods in estimating the magnetic field strength. To this end, wecreated synthetic data from every simulation and applied the DCF (with an optimization factor f =0.5) and ST methods. Both methods can be significantly affected when the magnetic field is pointing toward the observer since the LOS angle fluctuations induces extra dispersion in δθ (Ostriker et al. 2001; Falceta-Gonçalves et al. 2008; Hensley et al. 2019). For this reason, we assumed that B₀ lies completely in the POS.

To obtain an estimate of , we created spectroscopic data in the form of a position-position-velocity (PPV) cube, as in Miville-Deschênes et al. (2003). We used the following equation, (32)

where v_los(x, y, z) is the LOS velocity component and v is the central velocity of each velocity channel. This equation assumes optically thin emission. In Eq. (32), σ(x, y, z) is due to thermal broadening and is equal to , where k_B is the Boltzmann constant. We then fitted Gaussian profiles to every I_v spectrum and derived an “observed” line spread (σ_obs) as a free parameter of the fitting. We computed the turbulent velocity by subtracting in quadrature the thermal broadening, (33)

To obtain an estimate of δθ, we applied the formalism of Sect. 4.6 (Eqs. (20) and (21)) to calculate the Stokes parameters for each LOS, we estimated polarization angles through tan 2θ = U∕2Q, and calculated the dispersion of θ over the entire cloud through (taking θ = 0 in the direction of the mean magnetic field).

We then applied the DCF method with f = 0.5, and the ST method, by inserting the synthetic δu and δθ in Eqs. (3) and (5), respectively. We divided both equations with in order to derive the estimated magnetic field strength in Alfvénic speed units. Finally, we compared the estimated Alfvénic speed () with the actual value () of each simulation and computed their relative error as, (34)

In Fig. 5, we show the relative error of DCF (ϵ_DCF) versus ST (ϵ_ST). Red color is used for points when M_A = 0.1, magenta for M_A = 0.35, blue for M_A = 0.5, green for M_A = 0.7, black for M_A = 1.0 and cyan for M_A = 2.0. We did not include models with M_A = 2.0 and solenoidal forcing, since the polarization angle distribution of these models is uniform, and thus uninformative. We discuss the effect of forcing of these simulations in more detail in Sect. 5.2. The green shaded region corresponds to a smaller error for DCF than for ST.

We find that DCF is extremely inaccurate at low M_A ≤ 0.5: when the method fails, it fails by factors of several to tens. The accuracy of the method is improved for M_A ≥ 1.0, as expected, since it is in this regime for which the value of f we are using here (f = 0.5) has been optimized (Ostriker et al. 2001; Heitsch et al. 2001; Padoan et al. 2001). The overall trend of this figure is consistent with Fig. 2. The DCF method estimates are systematically biased toward larger values, because the kinetic energy of the cloud is much larger than the magnetic fluctuations, even for models with M_s = 0.5 and 0.7. The incomressible approximation employed by DCF is reached when M_s tends to zero, but even for weakly compressible flows (e.g., M_s ≈ 0.1) compressible terms can dominate the dynamics (Bhattacharjee et al. 1998). Only in trans-, super- Alfvénic cases DCF starts yielding reasonable estimates, because in this regime the contribution of the δB² ∕(8π) term in the energy increases.

In contrast, the error of the ST method remains low and uniform across the different models, because the kinetic energy remains comparable to the B₀ ⋅ δB fluctuations over a wide range of M_A (Fig. 2). The ST method overestimates the magnetic field strength at M_A = 0.1, while at M_A = 0.5 a transition happens. In the latter case half of the measurements overestimate the magnetic field strength, while the rest underestimate it. The underestimation is more prominent at large M_s, since the deviation between the estimated turbulent velocity and the true one becomes larger (Table 1). ST systematically underestimates the magnetic field strength at models with M_A ≥ 0.7. However, in all cases the error of ST is lower than 50%.

Fig. 5 Left panel: relative error of the DCF method versus the relative error of the ST method. Colors correspond to different MA as shown in the legend. Dots correspond to simulations with Ms < 1, “×” to 1 < Ms ≤ 4, and triangles to Ms > 4. The black dotted box marks the zoomed in region shown in the right panel. Right panel: zoomed region of the left panel.

5.1 Statistical properties of ϵ_DCF and ϵ_ST

In Fig. 6, we show the error distribution for the DCF method, ϵ_DCF, at different M_A, using kernel density estimation. The colorbar indicates the M_A of each distribution. Red, yellow, green, cyan and blue correspond to simulations with M_A of 0.1, 0.5, 0.7, 1.0, and 2.0, respectively. It is evident that as M_A increases ϵ_DCF decreases. Distributions become more symmetric and narrow at larger M_A since the properties of these models are more isotropic. The ϵ_DCF mean, median, and standard deviation of all the models with M_A ≥ 0.7 are 34%, 28%, and 49% respectively.This shows that the method estimates are systematically biased toward large values and the distribution is skewed to positive values. We also computed the same statistical quantities for the absolute values of ϵ_DCF for the same models. The mean, median, and standard deviation of |ϵ_DCF | is 47%, 37%, and 37%, respectively. DCF is completely inaccurate at lower M_A, and hence we do not report any statistics for these models.

The ϵ_ST distributions at different M_A are shown in Fig. 7. Colors are the same as in Fig. 6. The ST estimates are biased toward larger values at M_A = 0.1, while at M_A ≥ 0.7 the method estimates are systematically biased toward smaller values. The peak of the M_A = 0.5 distribution is close to zero and the probability of overestimating and underestimating the magnetic field strength is equal there. Distributions become more isotropic at larger M_A, which, similarly to ϵ_DCF, happens because the turbulent properties of these models are more isotropic. The ϵ_ST mean, median, and standard deviation of all the models are −2%, − 6%, and 24%, respectively and the distribution is close to a Gaussian. The same quantities for |ϵ_ST | are 20%, 18%, and 14%, respectively. The ϵ_ST distribution is more symmetric and peaks close to zero. Positive values are dominated by models with M_A ≤ 0.5, while negative by models with M_A ≥ 0.7. The right panel of Figs. 6 and 7 show the error distributions of the two methods (DCF and ST) on the same scale.

Fig. 6 Left panel: kernel density estimation of ϵDCF for different MA simulations. Right panel:zoomed in region of the left panel.

Fig. 7 Kernel density estimation of ϵST of simulations with different MA.

5.2 How does forcing affect the polarization data?

The dynamics of a sub-Alfvénic cloud are mainly determined by B₀ and forcing has a weak role. On the other hand, in super-Alfvénic turbulence forcing plays a dominant role in the cloud dynamics (Sect. 4.5). In Sect. 5, we could not apply the two methods (ST and DCF) in simulations with M_A = 2.0 that were driven solenoidally, but we did apply them to simulations with M_A = 2.0 and mixed forcing. The reason is that forcing has a strong impact on the distribution of polarization angles in super-Alfvénic simulations. In Fig. 8 we show the normalized polarization angle distributions of simulations with M_A = 2.0 and M_s = 2.0. The black histogram corresponds to a cloud driven solenoidally, while the blue histogram to a cloud driven with a mixture of modes. The black histogram is uniform and characterized by a spread equal to 48°, while the blue histogram has still a well defined mean with a spread equal to 28°. The dispersion in the black histogram (solenoidal driving) is so large it no longer encodes the magnetic fluctuations, due to the limited domain range, δθ ϵ [−90°, 90°]. In this case neither of the two methods can be applied.

6 Conclusions

The DCF and ST methods have been proposed to estimate the magnetic field strength from dust polarization data. DCF employs the incompressible approximation and infer that δθ ∝ M_A, while ST relaxes this assumption and infers that . In this work, we tested both scalings with synthetic data that we produced from ideal-MHD numerical simulations of isothermal clouds without self-gravity, spanning a wide range of M_s and M_A (M_s ϵ [0.5, 20] and M_A ϵ [0.1, 2.0]). In total we used 26 different models. We found that the synthetic polarization data can be fit very well with the ST scaling and not with the DCF scaling (Fig. 1). This means that the assumptions and approximations of the DCF method do not hold in compressible turbulence, while the ST assumptions are valid over the entire range of M_A studied here.

The major difference between the two methods is that in the DCF energy equation (Eq. (8)) by definition δB ⋅B₀ = 0, since only Alfvénic distortions are assumed to be present. On the other hand, ST relax this assumption and consider the more general case of compressible fluctuations where δB ⋅B₀≠0 (Eq. (17)).

We have explored whether the averaging over the total volume of a cloud can make this term vanish in the energetics, even if locally it is nonzero. We showed that ⟨δB ⋅B₀⟩ = 0 is obtained only if the zero-point of the magnetic “potential energy” is taken to be at δB = 0, rather than at the point of minimum potential energy, δB = −δB_max. By making a simple analogy with the problem of a bouncing ball in Newtonian gravity, we showed that for compressible fluctuations the correct energy conservation equation is obtained when comparing the kinetic energy with the rms of the B₀ ⋅ δB in the magnetic energy equation, which then naturally leads to the ST equation for the magnetic field strength (Eq. (17)). We have tested our theoretical arguments with numerical data and found that indeed the rms of B₀ ⋅ δB compares very well with the kinetic energy in all cases studied, with the exception of super-Alfvénic simulations with solenoidal forcing (see Fig. 2). On the other hand, when we omit this term, the kinetic energy can be up to two orders of magnitude larger than the magnetic energy. This explains why ST predict the right scaling between δθ - M_A, while DCF do not. This result highly supports the magnetic field fluctuations model of Federrath (2016) and Beattie et al. (2020).

We explored if B₀ ⋅ δB is imprinted in δθ. We showed analytically that when |B₀|≫|δB| the zeroth order approximation of δθ refers to perpendicular magnetic field fluctuations and the first order corrections correspond to compressible modes. However, the rms amplitude of parallel and perpendicular fluctuations is in all cases comparable, with deviations smaller than a factor of 2. As a result, δθ provides, indirectly, information about the amplitude of parallel fluctuations even to zeroth order.

We tested the accuracy of the two methods in estimating the magnetic field strength. We found that the DCF method with a “fudge factor” f = 0.5 failed completely in clouds with M_A ≤ 0.5 and only started producing reasonable estimates when M_A ≳ 0.7. The lowest errors for DCF were achieved for trans-,super- Alfvénic turbulence, M_A ≥ 1.0, because this is the regime for which the “fudge factor” f = 0.5 we have used was fine-tuned by Ostriker et al. (2001), Heitsch et al. (2001) and Padoan et al. (2001). Even in this regime (M_A ≥ 0.7 clouds), the relative error of DCF lied in the range [−37, 137]%. Over the entire M_A range, the error of DCF lied in the range [−37, 3500]%. The ST method on the other hand gave good results for all M_A examined here, without any fine tuning: the relative deviation from the true value lied in the range [−43, 51]% over the entire M_A. We did not find any strong dependence of the accuracy of the methods on M_s. Even in the cases where DCF would outperform ST, the ST method would still provide an adequate estimate of the magnetic field strength, while the reverse is not true.

Fig. 8 Polarization angle distributions for two different simulations with MA = 2.0 and Ms = 2.0. Black histogram corresponds to solenoidally driven cloud, while the blue to mixed forcing. Both histograms are normalized by dividing each bin with the total number of measurements. The legend shows the type of forcing.

Acknowledgements

We are grateful to Dr. B. Körtgen, Dr. S. Walch, Dr. D. Seifried, and Dr. P. Hennebelle for sharing their simulation data. We would like to thank C.F. McKee, P. André, Dr. E. Ntormousi, and Dr. V. Pelgrims for stimulating discussions, N.D. Kylafis, P. Goldsmith, and the anonymous referee for valuable comments on the draft. RS would like to thank Dr. K. Christidis for his constant support. This work was supported by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme under grant agreement No. 771282. V.P. acknowledges support from the Foundation of Research and Technology – Hellas Synergy Grants Program through project MagMASim, jointly implemented by the Institute of Astrophysics and the Institute of Applied and Computational Mathematics and by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “First Call for H.F.R.I. Research Projects to support Faculty members and Researchers and the procurement of high-cost research equipment grant” (Project 1552 CIRCE). J.R.B. acknowledges financial support from the Australian National University, via the Deakin PhD and Dean’s Higher Degree Research (theoretical physics) Scholarships, the Research School of Astronomy and Astrophysics, via the Joan Duffield Research Scholarship and the Australian Government via the Australian Government Research Training Program Fee-Offset Scholarship.

References

All Tables

All Figures

	Fig. 1 Polarization angle dispersion as a function of the Alfvénic Mach number. Blue line: ST scaling; magenta line: DCF scaling. The two lines are normalized so that they pass through the data for MA = 1.0.
In the text

	Fig. 2 Relative ratio of kinetic over magnetic energy density using Eq. (8) (vertical axis) and Eq. (16) (horizontal axis). Colors correspond to different MA as shown in the legend. Dots correspond to simulations with Ms < 1, “x” to 1 < Ms ≤ 4 and triangles to Ms > 4.
In the text

	Fig. 3 Same as in Fig. 2, but the dots correspond to models driven with a mixture of compressible and solenoidal modes and “×” denotes models driven solenoidally.
In the text

	Fig. 4 Ratio of perpendicular over parallel LOS averaged magnetic field fluctuations.
In the text

	Fig. 5 Left panel: relative error of the DCF method versus the relative error of the ST method. Colors correspond to different MA as shown in the legend. Dots correspond to simulations with Ms < 1, “×” to 1 < Ms ≤ 4, and triangles to Ms > 4. The black dotted box marks the zoomed in region shown in the right panel. Right panel: zoomed region of the left panel.
In the text

	Fig. 6 Left panel: kernel density estimation of ϵDCF for different MA simulations. Right panel:zoomed in region of the left panel.
In the text

	Fig. 7 Kernel density estimation of ϵST of simulations with different MA.
In the text

	Fig. 8 Polarization angle distributions for two different simulations with MA = 2.0 and Ms = 2.0. Black histogram corresponds to solenoidally driven cloud, while the blue to mixed forcing. Both histograms are normalized by dividing each bin with the total number of measurements. The legend shows the type of forcing.
In the text