Anisotropic Turbulence in Position–Position–Velocity Space: Probing Three-dimensional Magnetic Fields

The development of MHD turbulence theory is a long story. The earliest model of MHD turbulence is isotropic (Iroshnikov 1963; Kraichnan 1965). Later, this model was revised by a number of theoretical and numerical studies revealing that the MHD turbulence is anisotropic in sub-Alfvénic conditions and is isotropic in large-scale super-Alfvénic conditions (Montgomery & Turner 1981; Shebalin et al. 1983; Higdon 1984; Montgomery & Matthaeus 1995). Thereafter, a key concept of the "critical balance" condition, i.e., equating the cascading time (k_⊥ v_l )⁻¹ and the wave periods (k_∥ v_A )⁻¹, was given by Goldreich & Sridhar (1995), denoted as GS95. Here k_∥ and k_⊥ are the components of the wavevector parallel and perpendicular to the magnetic field, respectively. v_l is turbulent velocity at scale l, and v_A is Alfvén speed. By considering the Kolmogorov-type turbulence, i.e., v_l ∝ l^1/3, the GS95 anisotropy scaling be can easily obtained:

$$\begin{eqnarray}&&{k}_{\parallel }\propto {k}_{\perp }^{2/3},\end{eqnarray} \tag{ 1 }$$

which reveals the anisotropic nature of turbulence eddies, i.e., the eddies are elongating along the magnetic fields. Nevertheless, GS95's consideration is drawn in the global reference frame, i.e., the direction of wavevectors is defined with respect to the mean magnetic field. Due to the averaging effect, only the largest eddies are dominant in the global frame so that the observed anisotropy appears to be scale-independent (Cho & Vishniac 2000).

The scale-dependent anisotropy is later obtained from the study of fast turbulent reconnection (Lazarian & Vishniac 1999, denoted as LV99), which considers a local reference frame. The local reference frame is defined with respect to the magnetic field passing through the eddy at scale l. LV99 explained that only the motion of eddies perpendicular to the local magnetic field direction obeys the Kolmogorov law (i.e., ${v}_{l,\perp }\propto {l}_{\perp }^{1/3}$), because the magnetic field gives minimal resistance along this direction. Considering the "critical balance" in the local reference frame: ${v}_{l,\perp }{l}_{\perp }^{-1}\approx {v}_{A}{l}_{\parallel }^{-1}$, the scale-dependent anisotropy scaling is then:

$$\begin{eqnarray}&&{l}_{\parallel }={L}_{\mathrm{inj}}{\left(\displaystyle \frac{{l}_{\perp }}{{L}_{\mathrm{inj}}}\right)}^{\tfrac{2}{3}}{{\rm{M}}}_{{\rm{A}}}^{-4/3},{M}_{{\rm{A}}}\leqslant 1\end{eqnarray} \tag{ 2 }$$

where l_⊥ and l_∥ are the perpendicular and parallel scales of eddies with respect to the local magnetic field, respectively. L_inj is the turbulence injection scale. The corresponding anisotropy scaling for velocity fluctuation is:

$$\begin{eqnarray}\begin{array}{rcl}{v}_{l,\perp } & = & {v}_{\mathrm{inj}}{\left(\displaystyle \frac{{l}_{\perp }}{{L}_{\mathrm{inj}}}\right)}^{\tfrac{1}{3}}{{\rm{M}}}_{{\rm{A}}}^{1/3}\\ & = & {v}_{\mathrm{inj}}{\left(\displaystyle \frac{{l}_{\parallel }}{{L}_{\mathrm{inj}}}\right)}^{\tfrac{1}{2}}{{\rm{M}}}_{{\rm{A}}},{M}_{{\rm{A}}}\leqslant 1\end{array}\end{eqnarray} \tag{ 3 }$$

where v_inj is the injection velocity. This scale-dependent anisotropy in the local reference frame was numerically demonstrated (Cho & Vishniac 2000; Maron & Goldreich 2001) and in situ observations were made in the solar wind (Wang et al. 2016).

The observed intensity distribution in a given spectral line in PPV space is defined by both the emitters' density and their velocity distribution along the LOS. The LOS component of velocity v_los at the position (x, y) in the POS is a sum of the turbulent velocity v( r ) and the residual component due to thermal motions. This residual thermal velocity v_los − v( r ) has a Maxwellian distribution so that the gas distribution ρ(x, y, v) in PPV cubes and ρ( r ) in real space has the relation (Lazarian & Pogosyan 2004):

$$\begin{eqnarray}&&\rho (x,y,{v}_{\mathrm{los}})={\int }_{-S}^{S}\displaystyle \frac{\rho ({\boldsymbol{r}})}{\sqrt{2\pi {\beta }_{T}}}\exp \left[-\displaystyle \frac{{\left({v}_{\mathrm{los}}-v({\boldsymbol{r}})\right)}^{2}}{2{\beta }_{T}^{2}}\right]{dz}\end{eqnarray} \tag{ 4 }$$

where β_T = k_B T/m, m is the mass of atoms, k_B is the Boltzmann constant, and T is the temperature, which can vary from point to point if the emitter is not isothermal.

2.2.1. Dependence on M_A

In Figure 1, we consider a simple case that the magnetic fields are perpendicular to the LOS. The three isometric eddies (eddy₁, eddy₂, and eddy₃ in position–position–position, PPP, space) have different magnetic fields ( B ₁ > B ₂ > B ₃). For a given channel width Δv, for example Δv = 1 km s⁻¹, the amplitude of maximum LOS velocity v_los is then determined. However, here v_los is only contributed to by the perpendicular component of turbulent velocity. Hu et al. (2021) and Xu & Hu (2021) determined that the amplitude of velocity fluctuations for an isometric eddy is anisotropic, i.e., the maximum amplitude v_⊥ appears in the direction perpendicular to the local magnetic fields while the minimum amplitude v_∥ is in the parallel direction. In particular, the anisotropy of turbulent velocity (i.e., the ratio of ${v}_{\perp }^{2}/{v}_{\parallel }^{2}$) is a power-law relation with M_A:

$$\begin{eqnarray}{v}_{\perp }^{2}/{v}_{\parallel }^{2}=\left\{\begin{array}{lc}{\left({l}_{\parallel }/{L}_{\mathrm{inj}}\right)}^{-1/3}{{{\rm{M}}}_{{\rm{A}}}}^{-4/3}, & (\mathrm{local},{M}_{{\rm{A}}}\leqslant 1)\\ {{{\rm{M}}}_{{\rm{A}}}}^{-4/3}, & (\mathrm{global},{{\rm{M}}}_{{\rm{A}}}\leqslant 1)\end{array}\right..\end{eqnarray} \tag{ 5 }$$

Here "local" means the measurement is performed in the local reference frame, i.e., the parallel and perpendicular directions are defined by local magnetic fields. "global" means global mean magnetic fields define the global reference frame, i.e., the parallel and perpendicular directions.

Zoom In Zoom Out Reset image size

In our case, the LOS velocity v_los only consists of the turbulent velocity v_⊥, which is perpendicular to the magnetic field. For a given channel width, we have ${v}_{\mathrm{los}}={v}_{\perp }^{1}={v}_{\perp }^{2}={v}_{\perp }^{3}$. Consequently, a strong magnetic field ( B ₁ > B ₂ > B ₃) corresponds to a small value of v_∥ (i.e., ${v}_{\parallel }^{1}\lt {v}_{\parallel }^{2}\lt {v}_{\parallel }^{3}$, see Equation (5)). It indicates a more anisotropic turbulent velocity field v( r ). These three isometric eddies (in real PPP space) are then being mapped to the PPV space with identical channel width Δv. As shown in Equation (4), the observed intensity in a PPV channel is related to the distribution of turbulent velocity v( r ). ⁵ An anisotropic v( r ) would result in an anisotropic intensity distribution (Lazarian & Pogosyan 2000). Therefore, the observed intensity is more anisotropic for a strongly magnetized medium.

In a real scenario, the isometric eddies may incline to the LOS with angle γ. This inclination changes the observed anisotropy in the PPV channel. In Figure 2, we consider three magnetized isometric eddies (eddy₁, eddy₂, and eddy₃). The magnetic field strengths are identical (∣ B ₁∣ = ∣ B ₂∣ = ∣ B ₃∣), but B ₁ is perpendicular to the LOS, B ₂ inclines to the LOS with angle γ, and B ₃ is parallel to the LOS. We consider mapping the eddies to a given channel width Δv. For eddy₁, v_los is purely contributed by v_⊥. However, v_los of eddy₂ consists of both v_⊥ and v_∥. Considering the fact that v_∥ < v_⊥ and the ${v}_{\perp }\propto {l}_{\perp }^{1/3}$, the only way to achieve the same amplitude of v_los is to increase the eddy's size. Eddy₃ is facing a similar situation. As v_los of eddy₃ only comes from v_∥, the eddy must further increase its size to get the same v_los value. However, this increment of size results in a smaller anisotropy of turbulent velocity, i.e., a smaller ratio of ${v}_{\perp }^{2}/{v}_{\parallel }^{2}$ (see Equation (5)), as well as the observed intensity. This change of anisotropy induced by inclination angle is distinguishable from the one induced by magnetic field strength. It is clear that eddy₃'s anisotropy (i.e., γ = 0 case) is not observable in PPV space; however, this is not the case when only magnetic field strength gets changed (i.e., γ ≠ 0).

Zoom In Zoom Out Reset image size

The anisotropy of observed intensity in a PPV channel is thus constrained by both M_A and γ. By measuring two PPV channels' anisotropies, at least M_A and γ could be determined simultaneously. In the following, we will analytically show the anisotropy's dependence on M_A and γ through the second-order structure function.

The statistical description for a turbulent field within PPV space was first performed by Lazarian & Pogosyan (2000). Later, Lazarian & Pogosyan (2004) derived that, for optically thin lines, the intensity correlation function ξ_I (R, ϕ, Δv) is: ⁶

$$\begin{eqnarray}&&\begin{array}{l}{\xi }_{I}(R,\phi ,{\rm{\Delta }}v)\propto \displaystyle \frac{{\epsilon }^{2}\bar{\rho }}{2\pi }{\displaystyle \int }_{-S}^{S}{dz}[1+{\widetilde{\xi }}_{\rho }({\boldsymbol{r}})]{\left[{D}_{z}({\boldsymbol{r}})+2{\beta }_{T}\right]}^{-1/2}\\ \times {\displaystyle \int }_{-{\rm{\Delta }}v/2}^{+{\rm{\Delta }}v/2}{{dv}}_{\mathrm{los}}W({v}_{\mathrm{los}})\exp \left[-\displaystyle \frac{{v}_{\mathrm{los}}^{2}}{2({D}_{z}({\boldsymbol{r}})+2{\beta }_{T})}\right]\end{array}\end{eqnarray} \tag{ 6 }$$

where R is the two-dimensional separation of two points in the POS, ϕ is the position angle that R makes with the POS magnetic field, and Δv gives the channel width. is emissivity, $\bar{\rho }$ is the mean density, and S is the LOS distance. The integration along the LOS involves the overdensity correlation ${\widetilde{\xi }}_{\rho }{({\boldsymbol{r}})={\bar{\rho }}^{2}({r}_{0}/r)}^{\nu }$, the projected (along the LOS) velocity structure function D_z ( r ), and the thermal broadening term β_T = k_B T/m, where T is temperature, k_B is Boltzmann constant, and m is molecular/atomic mass . r = r ₁ − r ₂ means the three-dimensional separation of two points. The integration over velocities is described by the window function W(v_los). In terms of Alfvén waves, as they are incompressible and do not create any density fluctuations, the overdensity correlation ${\widetilde{\xi }}_{\rho }({\boldsymbol{r}})$ must be zero. The variables and parameters are summarized in Table 1.

Compared with the intensity correlation function, the intensity structure function D(R, ϕ, Δv) is better at describing intensity statistics at small scales (Lazarian & Pogosyan 2004). D(R, ϕ, Δv) can be expressed as:

$$\begin{eqnarray}&&D(R,\phi ,{\rm{\Delta }}v)=2[{\xi }_{I}(0,0,0)-{\xi }_{I}(R,\phi ,{\rm{\Delta }}v)].\end{eqnarray} \tag{ 7 }$$

The isotropic degree is defined as:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{iso}(R,\gamma ,{M}_{{\rm{A}}},{\rm{\Delta }}v) & = & \displaystyle \frac{D(R,0,{\rm{\Delta }}v)}{D(R,\pi /2,{\rm{\Delta }}v)}\\ & = & \displaystyle \frac{{\xi }_{I}(0,0,0)-{\xi }_{I}(R,0,{\rm{\Delta }}v)}{{\xi }_{I}(0,0,0)-{\xi }_{I}(R,\pi /2,{\rm{\Delta }}v)},\end{array}\end{eqnarray} \tag{ 8 }$$

where γ is the inclination angle of three-dimensional magnetic field with respect to the LOS, and M_A is the Alfvén Mach number. γ and M_A are implicitly included in D(R, ϕ, Δv) (see Appendix A). The dependence of iso(R, γ, M_A, Δv) will be shown below.

The analytical calculation of iso(R, γ, M_A, Δv) was carried by Kandel et al. (2016). Here we just briefly summarize it. Computing Equation (6) requires the knowledge of D_z ( r ), which can be obtained from the projection of the structure function tensor for the velocity field:

$$\begin{eqnarray}\begin{array}{rcl}{D}_{z}({\boldsymbol{r}}) & = & \langle ({v}_{i}({{\boldsymbol{r}}}_{1})-{v}_{i}({{\boldsymbol{r}}}_{2}))({v}_{j}({{\boldsymbol{r}}}_{1})-{v}_{j}({{\boldsymbol{r}}}_{2}))\rangle {\hat{z}}_{i}{\hat{z}}_{j}\\ & = & 2[(B(0)-B(r,\mu ))+(C(0)-C(r,\mu )){\cos }^{2}\gamma \\ & & -\ A(r,\mu ){\cos }^{2}\theta -2D(r,\mu )\cos \theta \cos \gamma ],\end{array}\end{eqnarray} \tag{ 9 }$$

which is determined by the coefficients A(r, μ), B(r, μ), C(r, μ), D(r, μ), and the angle θ between the LOS and r . We list the coefficients in Appendix A, and the derivation can be found in Kandel et al. (2016). By performing Legendre expansion for the coefficients up to the second order, i.e., A(r, μ) = ∑_n A_n (r)P_n (μ) ≈ A₀(r) + A₂(r)P₂(μ), and using the relation $\mu (\gamma ,\theta ,\phi )=\sin \gamma \sin \theta \cos \phi +\cos \gamma \cos \theta$, D_z ( r ) can be further simplified to:

$$\begin{eqnarray}&&{D}_{z}({\boldsymbol{r}})\approx {c}_{1}-{c}_{2}\cos \phi -{c}_{3}{\cos }^{2}\phi .\end{eqnarray} \tag{ 10 }$$

The parameters c₁, c₂, and c₃, which absorb the dependence on γ and M_A, are listed in Appendix A for the Alfvénic mode. As we can write $| {\boldsymbol{r}}| =\sqrt{{R}^{2}+{z}^{2}}$ and $\tan \theta =\tfrac{R}{z}$, the parameters only depend on R, z, γ, and M_A. By performing integration along the z-direction, the resulting isotropic degree iso(R, γ, M_A, Δv) is only the function of R, γ, M_A, and Δv.

Here we normalize Δv so that its maximum value is 1. In the limit cases of thin channel Δv = 0 and thick channel Δv = 1, defining:

$$\begin{eqnarray}\begin{array}{rcl}\alpha (R,\phi ) & = & {\left[{D}_{z}({\boldsymbol{r}})+2{\beta }_{T}\right]}^{-1/2},\\ \chi (R,\phi ) & = & \alpha (0,0)-\alpha (R,\phi ),\end{array}\end{eqnarray} \tag{ 11 }$$

the corresponding isotropic degree can be expressed as:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{iso}(R,\gamma ,{M}_{{\rm{A}}},{\rm{\Delta }}v=0)\\ =\displaystyle \frac{{\displaystyle \int }_{-S}^{S}\chi (R,0){dz}}{{\displaystyle \int }_{-S}^{S}\chi (R,\pi /2){dz}};\\ \mathrm{iso}(R,\gamma ,{M}_{{\rm{A}}},{\rm{\Delta }}v=1)\\ =\displaystyle \frac{{\displaystyle \int }_{-S}^{S}\chi (R,0){dz}{\displaystyle \int }_{-1/2}^{+1/2}W({v}_{\mathrm{los}})\exp \left[-\tfrac{1}{2}{v}_{\mathrm{los}}^{2}{\alpha }^{2}(R,0)\right]{{dv}}_{\mathrm{los}}}{{\displaystyle \int }_{-S}^{S}\chi (R,\pi /2){dz}{\displaystyle \int }_{-1/2}^{+1/2}W({v}_{\mathrm{los}})\exp \left[-\tfrac{1}{2}{v}_{\mathrm{los}}^{2}{\alpha }^{2}(R,\pi /2)\right]{{dv}}_{\mathrm{los}}}.\end{array}\end{eqnarray} \tag{ 12 }$$

Note that in the global reference frame, the anisotropy is scale-independent. This means the denominator and numerator of iso(R, γ, M_A, Δv) are both proportional to ∝ R^a , where a is a constant determined by the turbulence's properties (as D_z ( r ) ∝ r^ν ; see Appendix A and Kandel et al. 2016). Consequently, the measured iso(γ, M_A, Δv = 0) and iso(γ, M_A, Δv = 1) only depend on γ and M_A. The two values of iso(γ, M_A, Δv), therefore, are sufficient to determine γ and M_A for a given system.

Considering that the intensity structure function D(R, ϕ, Δv) is proportional to ${\cos }^{-1}\gamma$, ${\cos }^{-2}\gamma$, and ${M}_{{\rm{A}}}^{-2/3}$ (as D_z ( r ) depends on ${\cos }^{-2}\gamma$, ${\cos }^{-4}\gamma$, and ${M}_{{\rm{A}}}^{-4/3};$ see Appendix A), we separate the variables into: ⁷

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{iso}(\gamma ,{M}_{{\rm{A}}},{\rm{\Delta }}v) & = & ({a}_{1}+{a}_{2}\cos \gamma +{a}_{3}{\cos }^{2}\gamma )\\ & & \times \ ({b}_{1}+{b}_{2}{{\rm{M}}}_{{\rm{A}}}^{2/3})f({\rm{\Delta }}v)\end{array}\end{eqnarray} \tag{ 13 }$$

where a₁, a₂, a₃, b₁, and b₂ are parameters to be determined, and f(Δv) is a function of Δv. Note that the dependence on γ, which only involves in the projection of the structure function (see Equation (9)), is the same for all MHD modes, i.e., Alfvénic, fast, and slow modes. The M_A term comes from the amplitude of the power spectrum (see Appendix A). Nevertheless, the power spectrum of the slow mode is the same as that of the Alfvénic mode, and the fast mode is isotropic in the zeroth approximation (see Kandel et al. 2016 and Appendix C). Therefore, we expect the fitting function to work for the mixture of all MHD modes. For example, in Figure 3, at γ = π/2 and constant density, iso(γ, M_A, Δv) is negatively related to Δv, but positively related to M_A. In the following sections, we will perform numerical study to determine the parameters.

Zoom In Zoom Out Reset image size

The calculation of isotropic degree is performed by the second-order structure function, which is the SFA (Hu et al. 2021b; Xu & Hu 2021).

The first step is to determine the POS magnetic field. As illustrated in Figure 4, we choose a velocity channel that has an arbitrary width Δv. For this velocity channel, we calculate the intensity structure function D(R, ϕ, Δv):

$$\begin{eqnarray}&&D(R,\phi ,{\rm{\Delta }}v)=\langle | I({{\boldsymbol{X}}}_{1})-I\left({{\boldsymbol{X}}}_{2}\right){| }^{2}\rangle \end{eqnarray} \tag{ 14 }$$

where X ₁ and X ₂ denote the two-dimensional coordinates of two intensity data points located at position angle ϕ. Note that this calculation is preformed in the global reference frame, which means the anisotropy is scale-independent. Consequently at arbitrary R, D(R, ϕ, Δv) always exhibits a maximum value when ϕ is perpendicular to the POS magnetic field direction and a minimum value when ϕ is parallel to the POS magnetic field direction. Therefore, the POS magnetic field direction can be determined by varying the position angle from 0 to π. Note that in numerical simulations, there exists a fake dependence of the anisotropy on length scales, due to the isotropic driving and insufficient inertial range (Cho & Vishniac 2000; Yuen et al. 2018). To avoid this fake anisotropy, one should select R at sufficiently small scales away from the driving scale. In our case, we have R = 10 pixels (see also Figure 11).

Zoom In Zoom Out Reset image size

The second step is to figure out the LOS magnetic field direction and three-dimensional magnetization, i.e., the inclination angle γ and three-dimensional M_A. Here we define the isotropic degree $\mathrm{iso}(\gamma ,{M}_{{\rm{A}}},{\rm{\Delta }}v)=\max [D(R,\phi ,{\rm{\Delta }}v)]/\min [D(R,\phi ,{\rm{\Delta }}v)]$. To extract γ and M_A, one needs at least two measurements of iso(γ, M_A, Δv). For instance, one can choose to measure iso at normalized channel width Δv ≈ 1 and Δv ≈ 0. By solving Equation (13), one can get the values of γ and M_A.

In Section 2, we discussed the pure Alfvénic mode case neglecting the compressible components, i.e., fast and slow modes, as well density fluctuations. This simplification holds for subsonic turbulence, in which the density field, passively regulated by the velocity field, follows the statics of velocity fluctuation (Beresnyak et al. 2005; Xu et al. 2019). However, for supersonic turbulence, the turbulent compression is more significant, and the presence of shocks modifies the density field's statics (Kowal et al. 2007; Hu et al. 2020b; Xu & Hu 2021). The contribution from density can further change the anisotropy of the observed intensity structure. It is, therefore, essential to remove the density's contribution from channel maps.

In addition to the density effect, in a real scenario, one has to consider the thermal broadening effect in the velocity channel map, particularly for warm media. In this work, the temperature in simulations is set to 10 K, so the thermal broadening is marginal. Nevertheless, Yuen et al. (2021) recently developed a novel technique, i.e., the velocity decomposition algorithm (VDA), to reduce the density and thermal broadening effect in a given channel map. We briefly discuss the recipe of VDA here, and more details can be found in Yuen et al. (2021).

For a given intensity field ρ(x, y, v_los) in PPV space (see Equation (4)), we can define the integrated intensity map I(x, y) (i.e., the moment-0 map) and velocity channel map Ch(x, y) as:

$$\begin{eqnarray}\begin{array}{rcl}I(x,y) & = & {\displaystyle \int }_{-\infty }^{+\infty }\rho (x,y,{v}_{\mathrm{los}}){{dv}}_{\mathrm{los}}\\ {Ch}(x,y) & = & {\displaystyle \int }_{{v}_{0}-{\rm{\Delta }}v/2}^{{v}_{0}+{\rm{\Delta }}v/2}\rho (x,y,{v}_{\mathrm{los}}){{dv}}_{\mathrm{los}}\end{array}\end{eqnarray} \tag{ 15 }$$

where v₀ is the velocity of the averaged emission line maximum. Velocity fluctuation dominates the observed intensity fluctuations in the channel map when channel width Δv satisfies ${\rm{\Delta }}v\lt \sqrt{\delta ({v}_{\mathrm{los}}^{2})}$, where $\sqrt{\delta ({v}_{\mathrm{los}}^{2})}$ is the velocity dispersion (Lazarian & Pogosyan 2000). Accordingly, the velocity contribution Ch_v (x, y) and density contribution Ch_d (x, y) in Ch(x, y) can be extracted from (Yuen et al. 2021):

$$\begin{eqnarray}\begin{array}{rcl}{{Ch}}_{v}(x,y) & = & {Ch}-(\langle {Ch}\cdot I\rangle -\langle {Ch}\rangle \langle I\rangle )\displaystyle \frac{I-\langle I\rangle }{{\sigma }_{I}^{2}}\\ {{Ch}}_{d}(x,y) & = & (\langle {Ch}\cdot I\rangle -\langle {Ch}\rangle \langle I\rangle )\displaystyle \frac{I-\langle I\rangle }{{\sigma }_{I}^{2}}\end{array}\end{eqnarray} \tag{ 16 }$$

where ${\sigma }_{I}^{2}=\langle {(I-\langle I\rangle )}^{2}\rangle$ and 〈...〉 denotes the ensemble average.

Table 1. List of Notations Used in This Paper

Parameter	Meaning
r	three-dimensional separation
R	two-dimensional sky separation
ϕ	two-dimensional angle between R and the POS magnetic field
θ	angle between LOS $\hat{z}$ and r
γ	angle between LOS and mean magnetic field direction
μ	cosine of the angle between r and mean magnetic field direction
Δv	channel width
ξI (R, ϕ, Δv)	intensity correlation function
D(R, ϕ, Δv)	intensity structure function
Dz ( r )	z-projection of velocity structure function
W(v)	window function
βT	thermal velocity
MA	Alfvén Mach number
MS	sonic Mach number
vlos	LOS velocity
v( r )	turbulent velocity
v⊥	velocity fluctuation perpendicular to local magnetic field
v∥	velocity fluctuation parallel to local magnetic field
	emissivity
$\bar{\rho }$	mean density

Download table as:  ASCIITypeset image

The three-dimensional incompressible MHD simulations are produced from a pseudo-spectral code developed by Cho & Vishniac (2000). The code is a third-order-accurate hybrid employing an essentially non-oscillatory scheme. It uses hyperviscosity and hyperdiffusivity to solve the periodic incompressible MHD equations in a periodic box of size 2π:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t} & = & ({\rm{\nabla }}\times {\boldsymbol{v}})\times {\boldsymbol{v}}-{v}_{A}^{2}({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}+\nu {{\rm{\nabla }}}^{2}{\boldsymbol{v}}+{\boldsymbol{f}}+{\rm{\nabla }}P^{\prime} \\ \displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t} & = & {\boldsymbol{B}}\cdot {\rm{\nabla }}{\boldsymbol{v}}-{\boldsymbol{v}}\cdot {\rm{\nabla }}{\boldsymbol{B}}+\eta {{\rm{\nabla }}}^{2}{\boldsymbol{B}}\end{array}\end{eqnarray} \tag{ 17 }$$

where f is the random driving force, and $P^{\prime} =P+{\boldsymbol{v}}\cdot {\boldsymbol{v}}/2$ is the pressure. ν and η represent kinematic viscosity and magnetic diffusivity (η = ν = 6.42 × 10⁻³²), respectively. The magnetic field is considered to be B = B ₀ + δ B , where B ₀ is the uniform background field, and δ B is fluctuation. B ₀ initially is perpendicular to the LOS. v_A is the Alfvén speed of the uniform background. The code employs a pseudo-spectral method. We calculate the MHD equations in real space and transform them into Fourier space to obtain the Fourier components of nonlinear terms. The calculation of the temporal evolution is performed in Fourier space. Turbulence is solenoidally driven by 21 forcing components with $2\lt k\lt \sqrt{12}$ resulting in a peak of energy injection at k ≈ 2.5. Each forcing component has a correlation time of one. In our case, we vary the value of B ₀ to achieve M_A = 0.8 and M_A = 3.2 and stagger the simulation to 512³ cells/pixels. The respective parameters are listed in Table 2. We refer the reader to Cho & Vishniac (2000) and Cho (2010) for further details.

We generate three-dimensional compressible MHD simulations through ZEUS-MP/three-dimensional code (Hayes et al. 2006), which solves the ideal MHD equations in a periodic box.

$$\begin{eqnarray}&&\begin{array}{l}\partial \rho /\partial t+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0\\ \partial (\rho {\boldsymbol{v}})/\partial t+{\rm{\nabla }}\cdot \left[\rho {\boldsymbol{v}}{\boldsymbol{v}}+(P+\displaystyle \frac{{B}^{2}}{8\pi }){\boldsymbol{I}}-\displaystyle \frac{{\boldsymbol{B}}{\boldsymbol{B}}}{4\pi }\right]={\boldsymbol{f}}\\ \partial {\boldsymbol{B}}/\partial t-{\rm{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}})=0\end{array}\end{eqnarray} \tag{ 18 }$$

where f is a random large-scale driving force, ρ is the density, v is the velocity, and B is the magnetic field. We also consider a zero-divergence condition ∇ · B = 0, and an isothermal equation of state $P={c}_{s}^{2}{\rho }_{0}$, where P is the gas pressure, and c_s ≈ 0.192 (in numerical units) is the sound speed. The magnetic field and density field are considered to be B = B ₀ + δ b and ρ = ρ₀ + δ ρ, where B ₀ and ρ₀ are the uniform background fields. δ ρ and δ b represent fluctuations. Initially, B ₀ is assumed to be perpendicular to the LOS.

We consider single-fluid and operator-split MHD conditions in the Eulerian frame. The simulation is staggered to 792³ cells/pixels, and turbulence is solenoidally injected at wavenumber k ≈ 2 in Fourier space. The turbulence gets numerically dissipated at wavenumber k ≈ 100.

A scale-free MHD turbulence simulation is only characterized by the sonic Mach number M_S = v_inj/c_s and Alfvén Mach number M_A = v_inj/v_A , where v_inj is the isotropic injection velocity, and v_A is the Alfvén speed. We vary the value of B ₀ and ρ₀ to achieve different values of M_S and M_A. The parameters are listed in Table 2. In the text, we refer to the simulations by either the model name or the physical parameters.

We first examine the application of SFA to incompressible turbulence. Figure 5 presents the velocity channel maps (normalized Δv = 0.1) of incompressible simulations M_A = 0.8 and M_A = 3.2 at various inclination angles of mean magnetic fields. When the total mean magnetic field is perpendicular to the LOS (i.e., γ = π/2), the channel map's intensity structures (M_A = 0.8) are dominated by striations aligned with the POS magnetic field. The decreasing inclination angle, however, diminishes the anisotropic striation. When the total mean magnetic field is parallel to the LOS, the intensity structures become isotropic. As for super-Alfvénic turbulence, the intensity structures are always isotropic. In Figure 6, we calculate the intensity structure function D(R, ϕ, Δv) using the incompressible simulation M_A = 0.8 and choosing Δv = 0.1, R = 10 pixels, and γ = π/2. We vary the position angle ϕ from 0° to 180°. The maximum value of D(R, ϕ, Δv) appears at 0°, and 180°, while the minimum value appears at 90°. Note that there is a 180° ambiguity. Also, from the histogram of the POS magnetic field direction, we find the magnetic field direction concentrates at 90° with a mean value ≈8996. ϕ = 90°, which corresponds to the minimum D(R, ϕ, Δv), and therefore gives the direction of the mean POS magnetic field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Furthermore, we rotate the simulation cube so that the relative angle between the mean magnetic field and the LOS is γ. By varying the normalized channel width Δv, we plot the isotropic degree corresponding to different inclination angles in Figure 7. Its uncertainty is given by the standard error of the mean, which is negligible here due to the large sample size of the entire cube. We find the isotropic degree generally decreases when the channel becomes thick. The maximum and minimum values appear at normalized Δv ≈ 0.01 and Δv ≈ 1, respectively. This decrease can be understood as all thin channel emitters having similar LOS velocities, and anisotropy is suppressed. Also, when γ is smaller, i.e., the mean magnetic field is more parallel to the LOS, the observed anisotropy gets smaller as well. This decrease with respect to γ comes from the fact that the anisotropy is less projected onto the POS.

Zoom In Zoom Out Reset image size

In Figure 8, we fix the normalized channel width Δv to be 0.01, 0.10, and 1.00 but varying the value of γ, which is uniformly spaced in [0, 90°]. We can see the isotropic degree decreases when γ becomes larger, since more anisotropy is projected onto the POS. In the case of γ < 10°, the isotropic degree is 1, which indicates that the velocity channel is isotropic. In addition, we find the data points well fit the model $\mathrm{iso}(\gamma ,{M}_{{\rm{A}}},{\rm{\Delta }}v)=a{{\prime} }_{1}+a{{\prime} }_{2}\cos \gamma +a{{\prime} }_{3}{\cos }^{2}\gamma$. Note here that we already fixed Δv and M_A = 0.8. The fitting parameters are:

• 1.  
  Δv = 0.01: $a{{\prime} }_{1}=0.58\pm 0.04$, $a{{\prime} }_{2}=0.07\pm 0.18$, and $a{{\prime} }_{3}=0.33\pm 0.17;$
• 2.  
  Δv = 0.10: $a{{\prime} }_{1}=0.36\pm 0.03$, $a{{\prime} }_{2}=-0.06\pm 0.13$, and $a{{\prime} }_{3}=0.75\pm 0.12;$
• 3.  
  Δv = 1.00: $a{{\prime} }_{1}=0.19\pm 0.08$, $a{{\prime} }_{2}=-0.39\pm 0.38$, and $a{{\prime} }_{3}=1.24\pm 0.35.$

We find that $a{{\prime} }_{1}$ is small and $a{{\prime} }_{3}$ becomes large for a thick channel. This implies that the thick channel is more anisotropic and the γ has a more important role in regulating the thick channel's intensity structure.

Zoom In Zoom Out Reset image size

In this section, we test the SFA using compressible MHD simulations. Unlike the incompressible case, compressible slow and fast modes as well as the density field start to affect the anisotropy.

In Figure 9, we calculate the correlation of isotropic degree and normalized velocity channel width Δv at γ = π/2. We find for sub-Alfvénic turbulence, the isotropic degree is negatively related to Δv ≤ 0.3. When Δv ≥ 0.3, the isotropic degree starts increasing, which indicates smaller anisotropy. This can be understood as the density fluctuation dominating the thick channel so that the anisotropy is diluted (see also Figure 12 and Lazarian & Pogosyan 2000).

Zoom In Zoom Out Reset image size

In additional, in Figure 10, we further fix the normalized channel width Δv to be 0.01, 0.10, 0.30, and 1.00. We find it is clear that strong magnetic field cases exhibit more significant anisotropy. The super-Alfvénic case is closer to being isotropic, as the intrinsic turbulence is isotropic. We also fit the data points with the model $\mathrm{iso}(\gamma ,{M}_{{\rm{A}}},{\rm{\Delta }}v)=b{{\prime} }_{1}+b{{\prime} }_{2}{M}_{{\rm{A}}}^{2/3}$. Note here that we already fixed Δv and γ = π/2. The fitting parameters are:

• 1.  
  Δv = 0.01: $b{{\prime} }_{1}=0.52\pm 0.04$ and $a{{\prime} }_{2}=0.37\pm 0.06;$
• 2.  
  Δv = 0.10: $b{{\prime} }_{1}=0.16\pm 0.07$ and $b{{\prime} }_{2}=0.57\pm 0.10;$
• 3.  
  Δv = 0.30: $b{{\prime} }_{1}=0.08\pm 0.14$ and $b{{\prime} }_{2}=0.61\pm 0.19;$
• 4.  
  Δv = 1.00: $b{{\prime} }_{1}=0.19\pm 0.09$ and $b{{\prime} }_{2}=0.49\pm 0.13.$

The models' goodness of fits (see also Figure 8) with the data points confirm our theoretical expectation. We find that the fitted curve of Δv = 1.00 gets crossed with other curves. We expect that this comes from the density effect, as the Δv = 1.00 case keeps only density fluctuations. We will study the effect of the density contribution in the following section.

Zoom In Zoom Out Reset image size

5.2.1. Removing Density Contribution

From Equation (6), it is clear that when the channel is thick, i.e., Δv → ∞ , the observed intensity fluctuations only include the density's contribution (Lazarian & Pogosyan 2004):

$$\begin{eqnarray}&&\begin{array}{l}{\xi }_{I}(R,\phi )\propto {\displaystyle \int }_{-S}^{S}{dz}\displaystyle \frac{1+{\widetilde{\xi }}_{\rho }({\boldsymbol{r}})}{\sqrt{{D}_{z}({\boldsymbol{r}})+2{\beta }_{T}}}{\displaystyle \int }_{-\infty }^{+\infty }{{dv}}_{\mathrm{los}}\\ \exp \left[\displaystyle \frac{-{v}_{\mathrm{los}}^{2}}{2({D}_{z}({\boldsymbol{r}})+2{\beta }_{T})}\right]\propto {\displaystyle \int }_{-S}^{S}{dz}[1+{\widetilde{\xi }}_{\rho }({\boldsymbol{r}})]\end{array}\end{eqnarray} \tag{ 19 }$$

in which all of the velocity information is erased, and the density contribution plays a primary role in the observed intensity statistics. As velocity information is the most crucial in calculating the isotropic degree, we have to remove the density contribution in channel maps. To do so, here we use the VDA method (see Section 3).

In Figure 11, we use the compressible simulation M_S = 0.66 and M_A = 0.12, choosing Δv = 0.6 and Δv = 0.01 at γ = π/2. First, we plot the raw channel maps, and we find that the thin channel map (Δv = 0.01) displays more filamentary structures that are elongating along the mean magnetic field direction. We also calculate the intensity structure function with respect to the mean magnetic fields' parallel and perpendicular directions. The structure functions get shallower for the thin channel map as more small-scale structures appear. After that, we decompose the velocity and density contributions from the raw channel maps with the VDA method. We find for both the thick and thin channel maps that their corresponding density contribution maps are highly similar. Importantly, the (raw) thick channel map's intensity structures show more similarity with the density contribution map. In contrast, the (raw) thin channel map is similar to the velocity contribution map. This is exactly the theoretical prediction of Lazarian & Pogosyan (2000), i.e., the intensity fluctuations in the thick and thin channels are dominated by density fluctuations and velocity fluctuations, respectively.

Zoom In Zoom Out Reset image size

Also, the velocity contribution map exhibits more significant anisotropy in terms of the structure functions, i.e., the larger difference between the parallel and perpendicular components. The density contribution map's structure functions show higher amplitude for the thick channel, as the density fluctuation is more significant in the thick channel. In contrast, the velocity contribution map's structure functions for the thin channel, in which the velocity fluctuation is more important, have higher amplitude. We also analyzed the pure velocity caustic effect by setting a uniform density field when generating a PPV cube. ⁸ The VDA decomposed velocity contribution maps are highly similar to the pure velocity caustic maps, which confirm the validity of the VDA method. However, for the thick velocity caustic map, its structure function starts oscillating when the separation is larger than 20 pixels. We find that this comes from the fact that velocity information is marginal in a thick channel. For instance, for the pure velocity caustic case, because we set a uniform density field (i.e., ρ(x, y, z) = 1) for a 792³ MHD simulation, the intensity value at the density dominated position would be saturated to 792, which erases all velocity information. These saturated intensity values are excluded when calculating the structure function so that only velocity information is taken into account. The remaining velocity information does not guarantee an accurate structure function at a large scale from what we see. This insufficiency of velocity information is more significant if the channel width increases further.

In Figure 12, we make a comparison for the isotropic degrees calculated from the raw channel map, velocity caustic map, and VDA decomposed velocity contribution map of simulation M_S = 0.66 and M_A = 0.12. We find three crucial Δv ranges, at which the isotropic degree shows distinguishable behavior. For all three kinds of channel maps (i.e., raw channel map, velocity caustic map, and VDA decomposed velocity contribution map), the isotropic degree monotonically decreases when Δv gets larger until Δv ≈ 0.3. However, for the raw channel map, the isotropic degree starts increasing when 0.3 ≤ Δv ≤ 0.6. In contrast, the isotropic degrees of the velocity caustic map and VDA velocity map keep decreasing until Δv ≈ 0.6. As discussed above, the velocity caustic map and VDA velocity map contain only velocity information while the raw channel map includes both density and velocity contribution. The difference of the isotropic degrees in the range of 0.3 ≤ Δv ≤ 0.6 therefore comes from density contribution, which diminishes anisotropy. When 0.6 ≤ Δv, the isotropic degree of the raw channel map gets saturated since the projected density field's anisotropy is not related to channel width. The other two isotropic degrees, however, are dramatically divergent. This is because the velocity information in a very thick channel is not statistically sufficient for a structure function's calculation, i.e., the sample size is not large enough. For instance, in Figure 11, the structure function of the pure velocity caustic case fluctuates when the separation is larger than 20 pixels. When the channel becomes thicker, the remaining valid information of velocity fluctuations may drop down so that the isotropic degree at the 10 pixel scale (i.e., the numerical dissipation scale) gets diverged. Nevertheless, in observations, the inertial range and sample size are sufficiently large. One should perform the SFA at larger separations for thick channels to find sufficient velocity information. For supersonic turbulence, the density contribution is more significant. We discuss the case of supersonic turbulence in Appendix B.

Zoom In Zoom Out Reset image size

The POS magnetic field direction can be traced by varying the position angle used for calculating the intensity structure function (see Section 5.1). Combining the model of iso(γ, M_A, Δv) determined in this section, one can further access the inclination angle γ and the total Alfvén Mach number M_A. The three-dimension magnetic field information, including both orientation and strength, is achievable in PPV space.

As the isotropic degree iso(γ, M_A, Δv) for a given Δv depends only on γ and M_A, two measurements of iso(γ, M_A, Δv) are sufficient to determine a unique pair of γ and M_A, although multiple measurements could reduce uncertainty. Figure 13 considers a fitting model iso(γ, M_A, Δv) = $a+b\cos \gamma +c{\cos }^{2}\gamma$ + ${{dM}}_{{\rm{A}}}^{2/3}\cos \gamma +{{eM}}_{{\rm{A}}}^{2/3}+{{fM}}_{{\rm{A}}}^{2/3}{\cos }^{2}\gamma$, which is the expansion of Equation (13). We selected four normalized channel widths Δv = 0.01, 0.10, 0.30, and 0.60. We vary the values of γ for each compressible simulation M_S ≈ 0.6. By performing a two-variable fitting, we determine the coefficients and list them in Table 3. We find that the terms $\cos \gamma$, ${\cos }^{2}\gamma$, and ${M}_{{\rm{A}}}^{2/3}$ have higher weights in a thick channel. This means the anisotropy in a thick channel is more sensitive to $\cos \gamma$. Also, the weights of the ${M}_{{\rm{A}}}^{2/3}\cos \gamma$ term are close to zero when Δv > 0.01, which means their contribution is negligible. Note that density has an important role in the thick channel Δv ≥ 0.3. When calculating the isotropic degree, it is advantageous to take the thin channel width at the range of Δv ≤ 0.3.

Zoom In Zoom Out Reset image size

We test our fitting model in the compressible simulation M_S ≈ 0.89, M_A ≈ 0.54. We rotate the simulation so that the mean magnetic field is inclined to the LOS with angle γ. Following the recipe illustrated in Figure 4, we first determine the mean POS magnetic field direction. We use the alignment measure (AM; González-Casanova & Lazarian 2017):

$$\begin{eqnarray}&&{AM}=2{\cos }^{2}{\theta }_{r}-1\end{eqnarray} \tag{ 20 }$$

to quantify the relative angle θ_r , which is the difference between the orientation of the magnetic field inferred from the SFA and the real magnetic field. The AM value is in the range of [−1, 1]. AM = 1 means the difference between the two vector's orientations is zero, and AM = −1 indicates a 90° difference. To measure the mean POS magnetic field, we vary the position angle ϕ from 0 to 180° with resolution 1°. The position angle corresponding to the minimum value of the intensity structure function reveals the POS magnetic field direction. Also, here we select the velocity channel with widths Δv = 0.01 and Δv = 0.10.

In Figure 14, we find that the measured POS magnetic fields (using Δv = 0.01) have excellent agreement (AM ≈ 1) with the real POS magnetic fields when γ is larger than 18°. When γ is close to zero, the AM becomes negative. This can be understood from the fact that the POS magnetic field is isotropic in γ ≈ 0. Consequently, the intensity structure function gets similar results along all directions, which cannot determine the POS magnetic field. With the knowledge of the POS magnetic field direction derived from the SFA, one can further calculate the isotropic degree, which is defined as the ratio of the intensity structure functions measured in the direction parallel and perpendicular to the POS magnetic field. We calculate the isotropic degrees at Δv = 0.01 and Δv = 0.10 and then solve the fitting model iso(γ, M_A, Δv) = $a+b\cos \gamma +c{\cos }^{2}\gamma \,+{{dM}}_{{\rm{A}}}^{2/3}\cos \gamma +{{eM}}_{{\rm{A}}}^{2/3}\,+{{fM}}_{{\rm{A}}}^{2/3}{\cos }^{2}\gamma$. The uncertainties are given by the upper and lower limits of the fitting model listed in Table 3. As shown in Figure 14, the SFA measured γ gives AM ≥ 0.6 when the real γ is larger than 18°. In the case of γ < 18°, the AM dramatically drops down to negative and M_A is underestimated. The misalignment and underestimation come from the fact that the turbulent components dominate over the mean-field components at small γ. This effect is more significant in estimating total magnetic field strength. The bound for this underestimation theoretically is $\gamma \lt 4{\tan }^{-1}({{\rm{M}}}_{{\rm{A}}}/\sqrt{3})$ (see Lazarian et al. 2020). Also, we find the AM gets smaller than 0.8 when γ is in the range of [36°, 63°]. This comes from the deviation of the measured POS magnetic field, which affects the resulting isotropic degree. As the measured γ and M_A highly depend on the isotropic degree, a small fluctuation in the isotropic degree can lead to a significant deviation.

Zoom In Zoom Out Reset image size

Here we investigate the noise effect (i.e., the signal-to-noise ratio) on the measured isotropic degree. We use the compressible simulation M_S ≈ 0.6, M_A ≈ 0.12 as an example. We add Gaussian noise to each channel map. The noise amplitude varies from 10% to 50% of the channel's mean intensity value. In Figure 15, we find that the isotropic degree has two distinguishable behaviors. Starting from the thinnest case, the isotropic degree decreases. Then when Δv surpasses a turning point, the isotropic degree becomes positively related to Δv. The turning point, however, gets small for a strong noisy case, since more velocity fluctuation is overwhelmed by noise. Also, in all cases, the isotropic degree converges to a value of 1.0 at Δv = 1.0, which means the intensity map is isotropic. Nevertheless, even in the presence of noise, one should always choose a channel as thin as possible to calculate the isotropic degree.

Zoom In Zoom Out Reset image size

Probing magnetic fields in ISM is notoriously challenging. One option for measuring the three-dimensional magnetic fields is using polarized dust thermal emission (Chen et al. 2019). This method measures the inclination angle γ from:

$$\begin{eqnarray}&&{\sin }^{2}\gamma =\displaystyle \frac{{p}_{\mathrm{obs}}\left(1+\tfrac{2}{3}{p}_{0}\right)}{{p}_{0}(1+{p}_{\mathrm{obs}})}\end{eqnarray} \tag{ 21 }$$

where p_obs is the observed polarization fraction and ${p}_{0}=\tfrac{3{p}_{\max }}{3+{p}_{\max }}$, in which p_max is the maximum polarization fraction for a given cloud. It gives an integrated measurement of the magnetic field along the LOS.

To access the local three-dimensional magnetic fields, the anisotropic properties of MHD turbulence provide one solution. For instance, the fluctuation of synchrotron polarization has the imprint of magnetic fields. Using multiple-wavelength measurements of synchrotron polarization and taking the gradients of its wavelength derivative, the POS and LOS magnetic fields' orientation can be retrieved accordingly (Lazarian & Yuen 2018b). This method is called the synchrotron polarization derivative gradients (SPDGs).

Another possibility of tracing three-dimension magnetic fields through the anisotropy is enabled by young stars' three-dimensional velocity field, which is available in the GAIA survey (Ha et al. 2021). As the turbulent velocities measured in the directions perpendicular and parallel to the magnetic field can be significantly different, the value of the velocity's structure function could reveal the three-dimensional magnetic field (Hu et al. 2021a). Notably, the ratio of turbulent velocities measured in perpendicular and parallel directions of the magnetic field is proportional to ${M}_{{\rm{A}}}^{-4/3}$. The three-dimensional magnetic field's orientation and strength can therefore be determined simultaneously. This method is called the SFA.

This work further extends the applicability of SFA to the spectroscopic PPV cubes. Here we measure the intensity structure function in velocity channel maps with various channel widths. We find the isotropic degree, which is the ratio of intensity structure functions measured in parallel and perpendicular directions of the POS magnetic fields, satisfies iso(γ, M_A, Δv) = $a+b\cos \gamma +c{\cos }^{2}\gamma$+ $d\cos \gamma {M}_{{\rm{A}}}^{2/3}\,+{{eM}}_{{\rm{A}}}^{2/3}+f{\cos }^{2}\gamma {M}_{{\rm{A}}}^{2/3}$. The coefficients are listed in Table 3. By measuring the intensity structure function and isotropic degrees at two channel widths, one can simultaneously obtain the POS magnetic field, the magnetic field's inclination angle, and the total magnetization. Comparing with the dust polarization method and the SPDGs, the SFA gives complete measurements of three-dimensional magnetic fields, including both orientation and strength.

The SFA method is a revolutionary method that changes the stage of magnetic field measurements in ISM. It provides the possibility of simultaneously revealing the POS magnetic field, the LOS magnetic field, and total magnetic field strength using only a single spectroscopic PPV cube. This three-dimensional measurement of the magnetic field has benefited several studies.

The SFA is beneficial in modeling the foreground magnetic field and predicting the foreground dust polarization, which is indispensable in detecting the B-mode polarization of the inflationary gravitational wave. Previously, several efforts, including the rolling Hough transform (see Clark et al. 2015 and Clark & Hensley 2019) and the velocity gradients technique (VGT; see Hu et al. 2020c and Lu et al. 2020), have been made to predict the foreground dust polarization through atomic H I gas. These predictions, however, consider only the POS magnetic fields neglecting that inclination of the magnetic field, which is one of the significant sources of depolarization. The simplification, thus, results in a higher value of the predicted polarization fraction. However, our method can provide a considerably accurate foreground dust polarization map by incorporating the LOS magnetic fields and improving on these existing methods.

On the other side, the relative role of turbulence, magnetic fields, and self-gravity in star formation is a subject of intensive debate. Earlier measurements of the magnetic field are limited to two dimensions. The POS magnetic fields in molecular clouds are usually inferred from either far-infrared polarimetry (Andersson et al. 2015; Planck Collaboration et al. 2016) or VGT (Hu et al. 2019a, 2019b, 2021; Alina et al. 2020). The Zeeman splitting gives the signed LOS magnetic field strength (Crutcher 2012). Comparing to these methods, the SFA directly gives three-dimensional magnetic fields, which is crucial for an advanced understanding of star formation. The SFA can also be employed to study the magnetic fields in the galaxy's disk, where traditional far-infrared polarimetry suffers from distinguishing multiple clouds along the same LOS.

In ISM, turbulence is ubiquitous, and its injection scale L_inj is quite large. For example, Elmegreen & Scalo (2004) and Chepurnov & Lazarian (2010) reported L_inj ≈ 100 pc. Also, for diffuse ISM, the self-gravity is negligible so that turbulence and magnetic field are dominated there. These conditions well match the scale-free simulations of pure MHD turbulence used in this work. Therefore, our analysis is appropriate for the diffuse ISM regime.

At small scales, additional effects might become important. For instance, in star-forming regions, strong self-gravity, interstellar feedback, or outflows can be dominated. The stellar winds, jets, or bubbles can also change the picture of MHD turbulence. Consequently, additional studies are required for investigating turbulence at the scales or regions where these effects are important.

The VGT (see González-Casanova & Lazarian 2017; Yuen & Lazarian 2017; Lazarian & Yuen 2018a; Hu et al. 2018) is a new technique that has been developed to trace magnetic fields. As its theoretical foundation is the anisotropy of MHD turbulence, several properties of SFA can also immigrate to VGT. For example, the velocity gradient's direction indicates the eddy's semiminor axis, and the velocity gradient's dispersion is related to the eddy's size, i.e., the anisotropy. Since the eddy's size is regulated by M_A (see Figure 1), the velocity gradient's dispersion can reveal the magnetization (Lazarian et al. 2018). As we showed in this work, the anisotropy in PPV space is a function of channel width Δv, M_A, and inclination angle γ. The velocity dispersion is therefore correlated with these three variables. By analogy to SFA, we expect that two measurements of the velocity gradient's dispersion at two given channel widths can uniquely determine a pair of M_A and γ.

On the other hand, γ is crucial for achieving three-dimensional magnetic field strength. Lazarian et al. (2020) recently proposed a new technique to evaluate the POS magnetic field from two Mach numbers, i.e.,the sonic one M_S and the Alfvén one M_A:

$$\begin{eqnarray}&&B={\rm{\Omega }}{c}_{s}\sqrt{4\pi {\rho }_{0}}{{\rm{M}}}_{S}\ {{\rm{M}}}_{{\rm{A}}}^{-1},\end{eqnarray} \tag{ 22 }$$

where Ω is a geometrical factor (Ω = 1 corresponds to magnetic fields perpendicular to the LOS), c_s is sound speed, and ρ₀ is mass density. In the frame of VGT, the measurements of M_S and M_A come from velocity gradients' amplitude (Yuen & Lazarian 2020) and dispersion (Lazarian et al. 2018), respectively. Consequently, with the knowledge of γ, VGT can also achieve a measurement of three-dimensional magnetic fields, including direction and strength.